Modeling Neuromotor Adaptation to Pulsed Torque Assistance During Walking

Multiple mechanisms of motor learning contribute to the response of individuals to robot-aided gait training, including error-based learning and use-dependent learning. Previous models described either of these mechanisms, but not both, and their relevance to gait training is unknown. In this paper, we establish the validity of existing models to describe the response of healthy individuals to robot-aided training of propulsion via a robotic exoskeleton, and propose a new model that accounts for both use-dependent and error-based learning. We formulated five state-space models to describe the stride-by-stride evolution of metrics of propulsion mechanics during and after robot-assisted training, applied by a hip/knee robotic exoskeleton for 200 consecutive strides. The five models included a single-state, a two-state, a two-state fast and slow, a use-dependent learning (UDL), and a newly-developed modified UDL model, requiring 4, 9, 5, 3, and 4 parameters, respectively. The coefficient of determination (R2) and Akaike information criterion (AIC) values were calculated to quantify the goodness of fit of each model. Model fit was conducted both at the group and at the individual participant level. At the group level, the modified UDL model shows the best goodness-of-fit compared to other models in AIC values in 15/16 conditions. At the participant level, both the modified UDL model and the two-state model have significantly better goodness-of-fit compared to the other models. In summary, the modified UDL model is a simple 4-parameter model that achieves similar goodness-of-fit compared to a two-state model requiring 9 parameters. As such, the modified UDL model is a promising model to describe the effects of robot-aided gait training on propulsion mechanics.


I. INTRODUCTION
Robot-assisted gait training has been implemented in rehabilitation due to an assumed advantage of greater repeatability and consistency compared to standard approach [1].In robot-assisted gait training, the torque/force assistance to be provided during training may be modified based on the specific needs of a participant.By establishing a relationship between the participant response and assistive parameters, we can predict the optimal torque/force profile that yields a desired participant response.
To obtain the optimal torque/force assistance profile for an individual, it is important to consider variability among participants and to balance between number of observations and accuracy of modeling the participant's response [2].To accomplish both objectives at once, human-in-the-loop (HIL) optimization method was introduced [3].HIL and G. Kim  successfully implemented with various optimization methods including gradient descent [4], covariance matrix adaptive evolution strategy (CMA-ES) [5], and Bayesian optimization [6], [7].Previous research has also shown that incorporating additional information about the relationship between control parameters and participant responses reduces the number of iterations to achieve desired response during training [8].This result suggests that additional information may boost optimization, thus reducing the number of observations required to construct a precise model.However, previous optimization methods tested in HIL gait training did not account for the participant's adaptation during gait training.The optimizers assumed that participants would yield the same response when the same torque/force intervention was applied, which limits extensions of these method to rehabilitation applications.Therefore, there is a need to develop models of neuromotor adaptation in response to robot-assisted gait training.
Motor adaptation models describe how the internal model used by our central nervous system to control our movements changes in response to changes in task dynamics.To describe how the central nervous system refines forward models that are used to control movements, error-based learning models are widely used.In an error-based learning model, the internal model is updated by two processes: a feed-forward signal predicted from internal dynamics and an error-based component that updates the model based on the perceived error.Errors are often calculated based on the difference between the actual and planned movement or force.The twostate error-based learning model is one of motor adaptation models most widely used in previous research that focused on upper-extremity movements [9], [10].To address usedependent effects in the human response, where individuals often rely on past movements to decide the current movement, a use-dependent learning (UDL) model was introduced in previous research [11].However, the previously UDL model may describe adaptation during training but are not able to address remaining changes in responses after training, which is the ultimate goal of rehabilitation after neurological injury.Also, the previous models have not been validated within the context of walking mechanics.
In this work, we identified the relationship between hip and knee pulse torque applied to participant and output measurements during gait training using available data previously collected by our lab [12].A modified use-dependent learning (modified UDL) model was constructed by introducing a reference state update equation in the existing UDL model to model possible after-effects of gait training.State space models including single-state, UDL model, two-state fast and slow model, two-state model, and modified UDL model were compared in terms of goodness of fit to the participant-and group-level responses during and after gait training.In our analysis, we tested the hypothesis that the modified UDL model better describe participant responses, both during and after gait training, compared to existing models.

A. Data Collection and Processing
In previous research study, sixteen healthy participants were exposed to pulse of torque to the hip and knee joint while walking on the treadmill to modulate propulsion mechanics [12] (Fig. 1a).Eight different torque pulse pattern used in this research were generated based on the combination of pulse timing and amount of pulse torque to the hip and knee joint (Fig. 1b).In each trial of experiment, each participant was exposed to the same torque pulse intervention pattern for a continuous 200 strides.100 strides of baseline without torque pulse were conducted before torque intervention session, and same length of after-effect session without torque pulse was conducted after torque intervention session.Two outcome measurements were collected during previous research: Hip extension angle (HE) at peak anterior ground reaction force and normalized propulsive impulse (NPI).
Outcome measurements for both group averages and individual participant were processed separately for each participant, pulse condition, and outcome measurement.Each responses of HE or NPI was processed by filling missing data through linear interpolation.To reduce the impact of noise in HE and NPI, response at each stride was averaged with two preceding and two subsequent responses, separately for each session: baseline, torque intervention, and aftereffect.Given that participants were not expected to fully adapt to the experiment setting during the early baseline, responses from the 81st to the 100th strides in the baseline were averaged and subtracted from the overall responses to establish a reference response of zero value.The overall responses were normalized by dividing them by the standard deviation of responses from 81th to the end of strides of each experiment.

B. Modeling Use-Dependent Learning in the Response to Pulsed Torque Assistance
Five different motor adaption models were used to describe participant responses data from Sec. II-A.To describe gradual and non-exponential changes of participant response during and early after training, as well as the positive association between effects of training and after-effects of training previous observed in our lab [12], we selected the UDL model [11].However, because the original UDL model could not account for asymptotic after-effects, we modified the UDL model by introducing a reference state update equation.In addition to these UDL models, standard models of neuromotor adaptation were also included.Specifically, a two-state fast and slow model was selected since this model was often used in previous upper extremity research where a is the retention coefficient (i.e., it indicates how the previously planned movement affects the next movement), b is the error update coefficient (i.e., it indicates how the previously experienced error affects the next movement), and c and d indicate how the observed output y (i.e., HE or NPI in our data) is affected by the movement plan and by the applied perturbation u.Error e(n) is calculated as the difference between the planned and actual movement, i.e., e(n) = y(n) − x(n).The equation above can be transformed in the standard discrete-time state-space form as where and d).
2) Two-State Fast and Slow Model: In this model, a movement plan x(n) is defined from the linear combination of two separate states: a fast and a slow state [9].Each state has different dynamics in that the fast state reacts to a perturbation and returns to a reference value faster than the slow state.Since the slow state contains error feedback effects induced from the perturbation, the system may exhibit temporary after-effects when training is completed.The model equations are expressed as: where e(n) is described as error between movement plan state x(n) and output observed y(n), e(n) = y(n) − x(n) in this model.Since the slow state (x s ) needs to react slower than fast state (x f ), the retention and error update coefficients are constrained to verify the following equations (a s > a f and b f > b s ).This model can be converted in discrete-time state space form as: The two-state fast and slow model is defined via 5 parameters (a f , b f , a s , b s , and d).When the inertia coefficient d is equal to zero, the model becomes identical as described in previous implementation in upper extremity experiments [9], [10].
3) Two-State Model: The two-state fast and slow model is an instance of a two-state model of motor adaptation.Using the same error definition, procedure shown in Sec.II-B.1, a two-state model is expressed in the most general form as: The model is defined via 9 parameters, which are the component of state matrices A, B, C, and D. 4) Use-dependent Learning Model: The use-dependent learning (UDL) model is a single-state of motor adaptation.UDL model describes the use-dependent aspect of the human response observed during perturbation.This mechanism is observed simultaneously with error based learning and prediction of movement [11].The model equation is: Parameter a u accounts for retention, while b u accounts for use-dependent learning (i.e., how the previous output affects the next planned movement).In this model, x 0 is the reference movement in absence of any perturbation, and it is assumed to be a constant.Thus x 0 would account for the reference movement before, during, and after training.As the response data in our experiment results were processed to set the reference movements to zero, x 0 was set to be zero.Therefore, Eq. 13 can be transformed in standard state-space form as: The UDL model is defined via 3 parameters (a u , b u , and d u ).

5) Modified Use-dependent Learning Model:
The modified use-dependent learning (modified UDL) model is a two-state model describing error-based and use-dependent learning occurring during motor adaptation.The modified UDL model is derived by modifying the existing UDL model, introducing an update equation for the reference state, which will allow the model to describe residual effects induced by training.The model equation is c mu is the reference state retention parameter, bounded on [0 1]. a mu is the previous movement retention parameter, and b mu is the use-dependent learning term.The standard statespace form for the modified UDL model with state vector where (a mu , b mu , c mu , and d mu ) define the modified UDL model.All constraints and number of parameter to be estimated for each motor adaptation model are summarized in Table I.

C. Motor Adaptation Model Fitting
The five motor adaptation models from Sec.II-B were used to describe the data collected as described in Sec.II-A.As participant requires time to adapt to the experimental setting while walking on the treadmill, responses from the 81st stride to the end of the experiment were used to fit each motor adaptation model for each torque pulse condition.Constrained state space model estimation was conducted using MATLAB (MathWorks, Natick, Massachusetts, USA) to find optimal values of coefficients parameters in each motor adaptation model equations described in Sec.II-B in the least squares sense.Model fits were constrained to achieve a stable response, and to satisfy the additional validity constraints listed in Table I.To avoid local minima, fitting was repeated 100 times with different initial parameter values for the two-state fast and slow, UDL, and modified UDL models.Initial values for the single-and two-state models were determined based on the results of UDL and modified UDL models, respectively.Group-averaged response data was provided for each torque pulse condition for both HE and NPI.The best coefficient parameter estimations across the repeated runs, minimizing the sum of squared errors between the motor adaptation model and response, were selected for each motor adaptation model.The coefficient of determination (R 2 ) and Akaike information criterion (AIC) [13] were calculated based on each motor adaptive model estimations to quantify the goodness of fit and the understand the trade-off between model complexity and goodness of fit for the set of models considered.The same procedure was conducted for fitting motor adaptation models in individual participant responses.

D. Statistical Analysis
R 2 and AIC values from all motor adaptation models were analyzed through statistical software JMP (SAS Institute, Cary, North Carolina, USA).The Shapiro-Wilk test was conducted to check normality of the outcome measures from each motor adaptation models.If the outcome measures were normally distributed, a paired t-test conducted, otherwise a Wilcoxon signed rank test was conducted.To correct inference for multiple comparisons, the paired tests were conducted at a false positive rate of α < 0.05/10, using a Bonferroni correction for 5 motor adaptation models (10 pairs).Comparison between results of two motor adaptation models was conducted for each measurement condition: HE, NPI, both HE and NPI.

A. Model Fitting to the Group-averaged Response
Group-averaged responses were generated by averaging the response of sixteen participants for each torque pulse and outcome measurement condition.These group-averaged responses were then processed in Sec.II-A.Five different motor adaptation models were used to describe group-averaged responses for each torque pulse and outcome measurement condition.Four examples of motor adaptation model fitting to the group-averaged response of HE and NPI are shown in Fig. 2. Fitness values are presented in Table II

B. Model Fitting to the Participant-level Response
Five motor adaptation model introduced in Sec.II-B were utilized to describe responses of sixteen participants in terms of HE and NPI across different pulse conditions and outcome measurements (HE and NPI).
1) Statistical Analysis Results: All fitness values, including R 2 and AIC, from five motor adaptation models showed normal distribution.Paired t-test using matched pair test   II.
was conducted to check statistical difference between each fitness value from two different motor adaptation model, with corrected false positive standard.R 2 values showed significant difference for all pairs made from the 5 different motor adaptation models, except for the pair of single-state with UDL model.Statistical significance were maintain in all outcome measurement cases: HE, NPI, and both HE and NPI.Overall, two-state and modified UDL models exhibited significantly the highest and second-highest R 2 values compared to other motor adaptation models.
AIC values for fitting HE responses of participants indicated that all results were significantly different, except modified UDL with two-state model.For NPI case, UDL with two-state fast and slow model showed no statistical difference.For all response of HE and NPI, all pairs showed statistical difference.Overall, two-state and modified UDL model showed significantly lower AIC values compared to other motor adaptation models.
In all cases of participant-level response fitting, modified UDL and two-state model demonstrated better performance compared to single-state, UDL, and two-state fast and slow model.R 2 values indicated that two-state model showed significantly better fitness performance compared to modified UDL model (Fig. 3).However, when considering the number of coefficient parameters required to describe the response through each motor adaptation model, modified UDL model showed no statistical different AIC values when compared to AIC values from two-state model to describe participant responses of HE.

IV. DISCUSSION AND CONCLUSIONS
In this work, we compared performance of five different motor adaptation models describing group-averaged and individual participant responses collected in previously conducted experiment in our lab [12].Based on the goodness-offit values (R 2 and AIC), modified UDL model is expected to be the best model for describing group-averaged responses.
For addressing participant-level responses, both modified UDL and two-state model are suitable for response of HE, and two-state model is applicable for response of NPI.
Results of fitting to the group-averaged responses did not show statistical differences between fitness values from different motor adaptation models.Given that there were only 16 cases (8 pulses and 2 outcome measurements) in the group-averaged fitting results, the number of data points was not sufficient to demonstrate significant difference between fitness results.In most of cases from participantlevel and group-averaged fitting results, motor adaptation models were found to be more suitable to describe response of HE than NPI (R 2 values for HE: 0.61±0.27,NPI: 0.40±0.25).Kinematic response, specifically in HE, often showed larger changes in responses compared to NPI right after the introduction and removal of external perturbation in our experiment results.Therefore, in describing kinematic outcome measurements, motor adaptation models may be expected to be more effective than implementing them for kinetic outcome measurements.Not only that response of NPI shows less changes right after introduction and removal of perturbation compared to response of HE, NPI may contains more noise which was not eliminated via data processing procedure.Since group-averaged response were expected to contain less noise compared to participant-level responses, which may effect the fitness for describing responses using motor adaptation models.Effects on describing responses were demonstrated by the differences in R 2 values (groupaveraged: 0.72±0.25,participant-level: 0.49±0.28).
Considering capability of motor adaptation models to describe remaining after-effects, single-state and UDL model will eventually exhibit no after-effects when a large number of strides for after training session are provided in experiment.Therefore, for further implementation in gait training targeting after-effect modulation, modified UDL and two-state modes should be considered.When considering fitness values to describe both HE and NPI from group-level responses, modified UDL model performed best and superior compared to two-state model.However, this relationship was reversed when addressing participant-level response of NPI.This result suggests the possibility that two-state model performs better with noisy responses, while modified UDL model may describe better in absence of noise.For further research targeting kinematic variables such as HE or trailing limb angle during and after gait training using HIL optimization method, modified UDL model is expected to perform best with the additional noise elimination procedure during gait training.The modified UDL model will be implemented to reduce number of strides required for obtaining desired motion or to predict remaining after-effect, targeting remaining response changes after gait training.
is with the Department of Mechanical Engineering, University of Delaware, Newark, DE 19716, USA; secretgh@udel.eduF. Sergi is with the Department of Biomedical Engineering, University of Delaware, Newark DE, 19713, USA; fabs@udel.edu*This work is supported in part by NSF-CMMI-1934650 and in part by NIH-R01HD111071.

Fig. 1 :
Fig. 1: (a) Participant with the Active Leg Exoskeleton II (ALEX II) standing on the instrumented split-belt treadmill.(b) The torque pulse combinations applied to participants.

Fig. 2 :
Fig. 2: Model fitting for group-averaged responses in hip extension (HE) and normalized propulsive impulse (NPI) during experiments with selected torque pulse conditions.The gray area corresponds to the torque intervention session, where participants experienced torque pulses while walking on the treadmill.Corresponding fitness values (R 2 and AIC) are denoted in TableII.

Fig. 3 :
Fig. 3: Paired comparisons of goodness-of-fit parameters R 2 and AIC, reported as scatter plots (top-right), and paired line plots (bottom-left).Each row/column corresponds to one of the five models, and paired comparisons between models are reported in the intersecting cell.For the scatter plots, the unity line is provided in black, with R 2 values provided as red dots, and AIC values in blue.Large departures from the unity line are an indication of large differences of goodness-of-fit (note that smaller values of AIC indicate better fit).Asterisks indicate significant differences at the corrected false-positive rate.For the paired line plots, randomly undersampled ten outcomes between AIC values of 650 to 750 and R 2 values of 0.4 to 0.6 from all outcome measurement cases (HE and NPI) are displayed for improved figure clarity.Also, the AIC scale is inverted to align with the changes in R 2 values between different models.

TABLE I :
Model parameters and constraints and III, with the best values of goodness-of-fit highlighted in bold text. 1) Statistical Analysis Results: With the exception of R 2 values obtained from the single-state, UDL, and two-state fast and slow model, all other values, including AIC and R 2 values, exhibit non-normal distributions.Wilcoxon signed rank test was preformed to assess statistical difference in the distribution of R 2 value distribution among each model for different output measurement cases (HE, NPI, and both HE and NPI).The results of Wilcoxon signed rank test indicate that there was no significant difference between R 2 values derived from motor adaptation models in output measurement of HE or NPI cases.While comparing AIC values, two pairs, single-state with UDL and modified UDL with two-state model, showed significant difference with corrected false positive rate of 0.005.The Modified UDL model exhibited best AIC values

TABLE II :
Fitness values for group-averaged HE responses motor adaptation models in fifteen out of all sixteen cases.Additionally, both modified UDL and two-state model demonstrated the best R 2 values in fifteen out of sixteen cases.Two-state fast and slow model showed inferior R 2 values in three cases and higher AIC values in five cases compared to single-state or UDL model.

TABLE III :
Fitness values for group-averaged NPI responses