Human elbow motor learning skills of varying loads: Proof of internal model generation using joint stiffness estimation

This study presents a human elbow motor learning strategy responding to varying loads. Inspired by Kawato’s internal model theory, we suggest hypothesis that human minimize the internal model error by updating the joint stiffness to generate stable and robust motion during repetitive voluntary action with varying weight of load condition. We designed experimental robotics device to verify our hypothesis and the device is capable of precisely measuring human elbow joint stiffness very accurately. The subject was instructed to perform the prescribed elbow motion without notifying the weight of the load for neutral experimental condition and we recorded joint position, perturbation torque of actuator, reaction torque from torque sensor, and mean absolute value (MAV) of the surface EMG (sEMG) in forearm muscles and upper arm muscles as a reference criterion for elbow joint impedance modulation during motor learning. Modified ensemble-based system identification was applied to characterize the dynamic elbow mechanical impedance in transient state of moving loads. Experimental results show that subjects utilized high joint stiffness initially, but it decreases gradually and saturated to the level of 20%~60% of initial value after repetitive motion tests. The degree of saturation of motor learning varied with the weight of loads, this result supports the hypothesis that motor learning reduces joint stiffness by providing accurate internal model.


Introduction
Human motor learning skills play an important role in quality of performance in physical interactions with environments. Typically, for interaction with unknown objects, such as moving or carrying arbitrary object, agricultural labor, manufacturing tasks using tools, human adapt to these unstable conditions using dexterous motor control. Identifying the underlying mechanisms in motor learning provide fundamental insight into human-like robotic control, haptics, and rehabilitation fields. In previous studies, many researchers tried to investigate human motor learning hypothesis based on human experimental data with biomechanical and neurophysiological models using robotic technologies. The main focus of these studies was to develop motor control process how stable motion trajectory is generated in unstable conditions. Equilibrium point hypothesis (Burdet et al., 2001) were first introduced to explain stable motion is generated towards the equilibrium points using joint stiffness as a source of a feedforward torque without complicated computation of limb dynamics error using feedback control. However, it has been reported that in case of highly dynamic state such as fast motion, the target trajectory could be deviated from equilibrium point trajectory (Gomi et al 1997). From these results, others suggested internal model hypothesis for describing the source of error dynamics minimization with feedback learning process (Kawato et al 1999). They focused on mechanical joint stiffness and endpoint force during interacting with unstable environments to show the evidence of human motor learning process. Similarly, other research reports also supported this internal model hypothesis with and without force field perturbation in repeating planar arm reaching movements (Shadmehr and Mussa-Ivaldi et al., 1994). In order to identify the mechanisms in motor adaptation under unstable physical conditions, they investigated both human response after force field perturbation and after sudden removal of the force field by monitoring trajectory error corresponding each case. Characterization of human joint impedance provides intuitive information in describing the interaction between neuromechanical system and physical environment, connecting motor learning hypothesis with human experimental data. Especially for explaining general cases of dynamic motion, an in-depth understanding of joint characteristics during motion with transient state is essential. Multiple literatures (Lee et al., 2013(Lee et al., , 2014 investigated the joint impedance in several human joints such as elbow, shoulder, knee, ankle. Especially for elbow joint, much has been known about joint impedance in steady state (stationary) conditions, little is known for transient-state dynamic joint impedance as a general case of motion, due to lack of adequate tools to be used with high transparency. To best of our knowledge, only one study (Bennett et al., 1992) measured elbow joint stiffness in transient state. This study simply investigated joint impedance corresponding to repetitive motion without further experimental schemes for identifying motor learning skills in manipulating unknown objects with body dynamics. Some studies estimated joint stiffness in time-varying condition using surface EMG in human elbow joint (Li et al., 2019) and others studies the robotic skill learning using sEMG based stiffness estimation (Zeng et al., 2020).
In this context, this paper presents a human elbow motor learning skill from identifying dynamic elbow joint impedance in arm reaching motion with unknown varying loads. We first hypothesized that human minimize the trajectory tracking error while creating internal model in the brain. In course of repetitive motion, human optimize joint stiffness which represents viscoelasticity by updating the new parameters to generate stable and robust trajectory tracking. In order to test our hypothesis, we developed an experimental system capable of precisely measuring elbow joint stiffness when subject voluntarily rotate arm carrying unknown weight of loads. The detailed experimental setup, protocol and system identification methodology with data processing procedures are described. The schematic of this study is divided into the following two process. First, we checked motor learning process through repeating the motion multiple times ( = 1,2 ••• ) in three different loads, ( ) ( = 1,2,3), modulating joint stiffness, ( ) ( ) (for each different load, j in different trials, i), shown in (Fig. 1). Second process was to identify relationship between weight of load with joint time- Shin, Chang and Kim, Journal of Biomechanical Science and Engineering, Vol.16, No.3 (2021) [DOI: 10.1299/jbse.21-00088] varying stiffness among varying load experiments, ( ( ) ( )).The criterion for evaluating the level of motor learning was based on the extent of saturated stiffness and degree of target trajectory tracking error. Additionally, we checked the level of muscle co-contraction, evaluating the mean absolute value of the surface EMG data.

Internal model theory
Previous researches (Gomi et al., 1997) experimentally showed that in case of highly dynamic motion, human joint trajectory could be different from equilibrium point trajectory. From these experimental results, other research result asserted that inside the cerebellum of the brain, human build internal model of inverse dynamics through feedback learning process (Kawato et al 1999). Based on internal model, human create feed forward motor control input to generate stable motion. In order to minimize the trajectory tracking error, joint torque input was composed of torque from joint impedance modulation, −K(θ − θ d ) − Ḃ and feedforward torque, ̂̈ with updated limb dynamics from learned internal model as follows in 1-doF horizontal plane.

̈ + ( ,) =
(1) where g is unknown non-linear function approximation of net internal muscle and connective tissue torque, −̂+̈, K is joint stiffness, B is joint viscosity, I is an inertia of the human forearm. Here, the term, = −1 in (2) represents the joint torque from internal model, as feedforward-based term stored in cerebellum of the brain. The next term (2) shows the joint impedance modulated torque, as a result of viscoelasticity ( Fig. 2.). This viscoelasticity can be regarded as the peripheral feedback control gain, modulated by regulating muscle cocontraction levels and reflex gains (Franklin et al., 2008). Substituting (2) for in (1), we can arrange internal model error, ∆ with joint impedance related equation (5) and error dynamics (6) is like follows.

= −̂+̈
where ∆ = −̂ as an internal model error. From obtained equations (6), we could expect near the equilibrium points, joint stiffness is K inversely proportional to internal model error. Based on the internal model theory, internal

Stiffness-based trajectory error compensation
Based on current motor learning theory, we hypothesized human motor control strategy of joint stiffness optimization to generate feed forward torque. This torque minimized the trajectory tracking error from unknown high inertial loads from learned inverse dynamics of limb. According to internal model hypothesis, trajectory error (e) depends mostly on joint stiffness, damping with internal model error, simply described as following relation (7).
where ∆ is internal model error, is inertial term, and is a joint stiffness. This represents the tendency of trajectory error, without detailed process of time-varying joint stiffness with corresponding inverse dynamics model error. In the  Shin, Chang and Kim, Journal of Biomechanical Science and Engineering, Vol.16, No.3 (2021) [DOI: 10.1299/jbse.21-00088] following sessions, we experimentally measured the time-varying elbow joint stiffness during repetitive elbow motion with unknown varying loads. The tendency result of stiffness variation within each trail and among different trials could provide human joint stiffness as viscoelasticity adaptation strategy, compensating internal model error in trajectory following tasks.

Hardware Design
In this study, we developed the robotic measurement system which can (1) perturb the human limb with enough power to achieve measurable data, (2) excite the human elbow joint with proper level of frequency bandwidth for system identification (3) robust mechanical structure to support high weight load in axial direction (4) wearable type arm connecting design with multiple cases for unknown slots (5) axis alignment for efficient transmission of perturbation to human. The source of the actuator is frameless type BLDC motor (EC frameless 90 flat 260Watt, Maxon, Switzerland). This motor shows 0.968[Nm] nominal continuous torque output, which is adequate to transmit the high frequency region energy to both human elbow with varying loads, satisfying the power condition (1). For satisfying Nyquist sampling theorem, the frequency contents of the input signal should be at least more than twice of the system, especially human elbow in this research. According to previous literatures on human elbow joint natural frequency, they reported it operates in nearly 25Hz. Therefore, we investigated the power spectrum density (PSD) of the actuator in frequency bandwidth to check the performance of the system. The result of power spectrum of the mild random torque signal was flat to 53Hz, as represented in Figure 4. This satisfies the requirement of frequency bandwidths specification (2). The mechanical design of thrust bearing with load support structure could disperse the high load, on the inner radius, section of the bearing (Fig. 3. a). Also, the double bearing structures could prevent the bending of the forearm structures from the human arm with additional loads (Fig. 3. a, b). Three slot cases are connected to forearm supporter, from 100g to 2500g. Forearm supporter is wearable type structures to directly align the axis of elbow joint with actuator. Also, the shoulder joint is isolated to fixed position with support fixture for upper arm. These structural features of the system satisfy  1,2,3). For each load, subject trained without perturbation (Before instant A) input and trained after torque perturbation for impedance measurements, gradually learning the unknown dynamics B) Instance of the actual motor learning is started. The trajectory error level in this period is over 40%. C) Instance of the motor learning is almost done. The trajectory error level is near 10%. D) Instance of the trajectory error is saturated.

Experimental Setup
We designed elbow joint impedance measurement system for motor learning experiment. A forearm of the subject is located in parallel to manipulator, fixed using strap (Fig 5. (a)). Also, the upper arm is fixed to the supporting structure to isolate the elbow motion, minimizing the intervention of shoulder. Surface electromyography (EMG) sensors (Trigono, Delsys Inc. USA) were used to monitor the muscle activation during the elbow motion. Each of four surface EMG sensors were attached to the muscle belly of both forearm and upper arm muscles: Brachioradialis (BR), Flexor carpi ulnaris (FC), Biceps (BC) and Triceps (TC) each of muscle corresponding to (1), (2), (3), (4) in Fig 5.(a). Angular data from and perturbation torque data were acquired in LabVIEW FPGA system (NI compact RIO), sampled at 1kHz. EMG data were also sampled at 1kHZ and band-pass filtered between 20Hz and 450Hz.

Experimental Protocol
Seven healthy male subjects (right-handed male subjects) participated in this study. Subjects were between age of 21-26[mean (SD): 23.83(1.47)] with no history of biomechanical or neurological disorders. Following protocol was approved by KAIST institutional review board, informed consent from all participants.
For the first step, we attached four surface EMG electrodes in each of four muscle belly and instructed subjects to follow the trajectory in the visual panel without torque perturbation using only elbow joint to rotate without intervening the shoulder joint. The magnitude of random mild torque perturbation of each subjects was decided after monitoring the trajectory following error. When the degree of trajectory error reached the tolerance (±10%) and confirmed that participant could move the elbow joint without aware of the perturbation signal, we selected the corresponding torque signal (~0.93Nm) for our input. Next, subject rotated elbow joint watching the reference signal in visual panel to follow the sinusoidal trajectory with 20° magnitude and 0.6Hz frequency (Fig. 4.). A band-limited Gaussian white noise was applied to produce random torque perturbations transmitted to human elbow joint. For one session, we fixed the carrying load and connected it to robotic manipulator not revealing the weight to subject and repeated the same motion nearly 100 times. The detailed experimental design is described in (Fig.5 (b)).

System Identification Algorithm
Modified ensemble-based linear time-varying system identification technique was used in this study to characterize the dynamic elbow joint impedance in transient state (Lortie et al., 2001).
We set input of the system as a perturbation torque of the system. The input satisfies the noise-free condition. The output of each realizations in ensemble set was the difference between measured elbow position( ) and nominal position ( 0 ) from multiple ensemble set of the trials. The input-output relationship is determined using discrete convolution equation as follows (9).
where ĥ( , ) is an impulse response function (IRF) estimate (Fig. 7), one type of time-varying weightings for each time instants i, with a finite lag length = + 1, long enough for the IRF estimate to reach down almost to zero. By multiplying both sides of (9) with ( − ), summing over all realizations and averaging both sides, we obtain the following equation (10).
Left side is an input-output cross-correlation function estimate (∅ ) at time and − . Right side is an input auto- Fig. 7 Top is the obtained impulse response function (IRF) ̂( ) in two variables, time instant i and lag j. Bottom graph represents the IRF after smoothing process.
Shin, Chang and Kim, Journal of Biomechanical Science and Engineering, Vol.16, No.3 (2021) [ The auto-correlation function estimate(Ф ) is M × M matrix and cross-correlation function estimate (Ф ) and impulse response function (IRF) estimate ̂( ) is a M × 1 vector, rewritten as follows.
Based on the previous study that human elbow stretch reflex response is in excess of 40ms (Koike et al., 2007), we set the length of the reflex loop delay to 40ms (Sampling rate of the system is 1kHZ). For each time , the optimal timevarying parameters were estimated by minimizing the mean square error, ( ) between smoothed IRF, ̂( ) based estimated torque with measured torque ( ) (Fig. 6).

Muscle Co-Contraction during Motor Learning
We evaluated the muscle co-contraction in forearm and upper arm muscle using mean absolute value (MAV) of recorded surface EMG. It provides the degree of muscle activation during elbow motion. The MAV of surface EMG during motor learning represents the reference index of viscoelasticity, for validating the obtained joint impedance modulation history in our research. Each of the MAV in four muscle near elbow joint is like following equation (16).
where i is N is the length of data, is the raw EMG signal. Shin, Chang and Kim, Journal of Biomechanical Science and Engineering, Vol.16, No.3 (2021) [DOI: 10.1299/jbse.21-00088]

Dataset Construction & Processing
We reconstructed the measured position and torque data with criterions applied for statistical ensemble-based system identification method. To satisfy the assumption that every realization undergoes the same behaviors, we normalized the temporal and spatial size of the data and eliminated the outlier outermost 10% of realizations from the nominal dataset (Fig. 6). Also, the obtained impulse response function ̂( ) which represents the system characteristics was smoothed using time window and showed gentle shape near the zero lag regions. This shows the good result of matrix computation of correlation functions and optimizations.

Reliability of the Model Approximation
The reliability of the parameter estimation was calculated and shown in Fig. 8 using NRMSE between estimated (̂) and measured position( ) from ensemble set. The approximation of statistically obtained impulse response function at particular time instants, including multiple neighborhood data of time lag, with mechanical second order system underwent the optimization process using nonlinear least squares and time-windowing process to smooth the dataset. The overall values of the NRMSE were less than 2% within trials (Fig. 8). This result shows the goodness-of-fit of the model approximations we applied.

Joint Impedance Transition after Motor Learning
Motor learning skill is described with the estimated time-varying mechanical impedance, especially joint stiffness (main contribution of output torque, more than 85% of total torque. Fig. 9 (a) represents the trend of stiffness transition after trials. The stiffness transition trend in five main positions commonly showed initially continued for some time, gradually decrease to lower level compared to start and finally trials for joint stiffness saturation among different weight loads is in Fig. 9. (c), (d) showing the stiffness ratio and speed of motor learning after repetitive motion. Also, with multiple trials, the trajectory error decreases and saturated to nearly 10% as shown in Fig. 9 (b). Shin, Chang and Kim, Journal of Biomechanical Science and Engineering, Vol.16, No.3 (2021) [DOI: 10.1299/jbse.21-00088]

Evaluation of Viscoelasticity using surface EMG
MAV of surface EMG in both forearm and upper arm initially increased, gradually decreased and saturated to relatively lower level compared to initial state. The reduced rate of MAV after motor learning was nearly 2% ~ 4% (Fig .10). It supports our hypothesis that human modulate joint stiffness (viscoelasticity) to minimize internal error.

Discussion and Conclusion
This paper suggests human elbow motor learning strategy of carrying varying loads. Based on previous internal model theory, we hypothesized that human minimize the internal model error during repetitive voluntary action with varying weight load condition by updating the joint stiffness to generate stable. We validated our hypothesis by characterizing the time-varying elbow joint stiffness using designed transparent actuating system.
We identified that human modulate joint stiffness, initially relatively high value, gradually reduced the level nearly 20 ~ 60% (different among weight of loads) and trajectory error also showed similar trends. Also, we checked the mean absolute value (MAV) of the surface EMG to compare the tendency of joint stiffness transition with actual human muscle activation, as a reference index for viscoelasticity of the joint. From these results, we concluded that accuracy of internal model is modified during repetitive motion after trials and human stabilize the motion in unknown dynamic condition modulating joint stiffness by updating the optimized value from internal model error. Our motor learning results provide insight about human joint control of unknown dynamics. Motivated by motor learning skill of human from this research, we can simply modulate joint stiffness of robot joint like human and realize the stable motion (Chang et al., 2017). simply updating the ratio of joint stiffness, without complicated and heavy load calculation of body dynamics and optimization in each instant during motion (Chang et al., 2019).