Abstract
Wearable assistive technology for the lower extremities has shown great promise towards improving gait function in people with neuromuscular injuries. But common secondary impairments, such as hyperreflexia, have been often neglected. Adding hyperreflexia prediction to the control loop would require expensive or complex measurement of muscle fiber characteristics. In this study, we explore a clinically accessible biomechanical predictor set that can accurately predict rectus femoris (RF) reaction after knee flexion assistance in pre-swing by a powered orthosis. We examined a total of 14 gait parameters based on gait kinematic, kinetic, and simulated muscle-tendon states from 8 post-stroke individuals with Stiff-Knee gait (SKG) wearing a knee exoskeleton robot. We independently performed both parametric and non-parametric variable selection approaches using machine learning regression techniques. Both models revealed the same four kinematic variables relevant to knee and hip joint motions were sufficient to effectively predict RF hyperreflexia. These results suggest that control of knee and hip kinematics may be a more practical method of incorporating quadriceps hyperreflexia into the exoskeleton control loop than the more complex acquisition of muscle fiber properties.
Introduction
Gait disorders are common in stroke survivors. More than 80% of stroke survivors have varying degrees of gait abnormalities11,21, and about 25% have a residual impairment that requires full physical assistance, despite rehabilitation efforts17. Recent technological advances have introduced wearable robots as potential solutions to mitigate physical barriers for people with impaired mobility. There are various commercialized and research purpose lower-body exoskeletons showing promise in performance augmentation, mobility assistance, and gait therapy following neurological injuries12,41. Assistive strategies using exoskeletons incorporate force, torque, electromyography (EMG), and gait kinematics into the control loop22,41. More recently, human-in-the-loop approaches have been introduced to maximize the benefit of the assistive device with feedback of metabolic energy measures in a real-time optimization43. However, existing exoskeletons have only addressed weakness and ignored critical neuromuscular impairments, such as hyperreflexia, defined as a hypersensitive stretch reflex, governed by musculotendon fiber stretch velocity42.
Stiff-Knee gait (SKG) is a common abnormal gait pattern following stroke, characterized by diminished knee flexion during the swing phase of the gait cycle. We previously developed a lightweight, remotely actuated powered knee orthosis38 to assist post-stroke individuals with SKG. This robotic knee exoskeleton was expected to improve swing-phase knee flexion kinematics. While kinematics improved somewhat, hyperreflexia in the rectus femoris (RF), represented by excessive early swing-phase muscle activation, primarily followed exoskeletal assistance, but occasionally without assistance as well2. Musculoskeletal modeling and simulation revealed that increased RF fiber stretch velocity preceded increases in RF muscle activation2. In a new cohort of individuals with SKG, Akbas et al. found that RF reflex excitability was highly associated with reduced swing phase knee flexion angle1. Taken together, this evidence indicates that robotic knee perturbations on post-stroke individuals with SKG could elicit counterproductive quadriceps hyperreflexia, a reaction likely spurred on by suprathreshold RF fiber stretch velocity. Thus, if excessive RF fiber stretch velocity could be avoided, it may be possible to enable the full benefits of exoskeletal assistance in those with post-stroke SKG. However, real-time measurement of fiber stretch velocity requires difficult to access equipment such as ultrasound and a real-time analysis and control system29. Another approach could be to identify clinically accessible correlates of fiber stretch velocity to avoid the hyperreflexia elicited by exoskeletal assistance.
The objective of this study was to determine to what degree kinematic and kinetic features could substitute for RF fiber stretch velocity in predicting hyperreflexive RF muscle activation following exoskeletal assistance. To achieve this, we performed parametric and nonparametric multivariate regression analyses for variable selection on a knee exoskeleton gait data from 8 post-stroke participants with SKG38. We hypothesized that hip and knee flexion kinematics would predict RF hyperreflexia due to the biarticular characteristics of the RF. This work represents a novel statistical and machine learning regression approach to identify predictors of hyperreflexia. Obtaining clinically accessible predictors for hyperreflexia will lead to more effective human-in-the-loop assistive strategies that account for neural impairments.
Materials and Methods
Experimental Data
We obtained previously collected data from 8 chronic, hemiparetic participants with post-stroke SKG who gave written informed consent using procedures approved by the local Institutional Review Board37. Inclusion criteria for the hemiparetic participants included their peak paretic knee flexion during swing was at least 15° less than their peak nonparetic knee flexion and the ability to walk for 20 min without rest at 0.55 m/s on a treadmill. A lightweight, powered knee orthosis was used to provide knee flexion torque perturbations during the pre-swing phase without affecting the remainder of the gait cycle38. Analog data, including GRF from the instrumented split-belt treadmill (Tecmachine, Andrez Boutheon, France) and applied torque by the powered orthosis, were acquired at 1 kHz. Motion capture data (Motion Analysis, Santa Rosa, CA) on the lower limbs was collected at 120 Hz. Electromyography (EMG) (Delsys Inc., Boston, MA) from the RF muscle was measured at 1 kHz.
The experimental protocol involved steps with and without pre-swing knee flexion assistance. Steps with knee flexion assistance ranged from 10 Nm to 40 Nm. A mean of less than 1 Nm of resistance was measured on the knee for the remainder of the gait cycle. The weight of the device did not alter walking kinematics. In this study, a total of 406 gait cycles were extracted including 260 gait cycles with knee flexion assistance (Mean ± SD: 23.32 ± 5.85 Nm; Min: 11.40 Nm; Max 34.51 Nm). Previous work describes further details of the protocol37.
Musculoskeletal Modeling and Simulation
In order to determine muscle-tendon states (i.e., RF muscle fiber stretch velocity at each time instance), we employed musculoskeletal modeling and simulation through OpenSim 4.310 using a method validated in previous work2,4. To summarize, we condensed upper body segments in the gait 2392 model of OpenSim to the pelvis segment to account for a lower body marker set, so that it had 18 degrees of freedom and 90 muscle-tendon actuators. The modified model was scaled to match the anthropometry of each subject. Marker data were fed into the inverse kinematics tool, which generated joint kinematics with the least square fit of marker trajectories. A residual reduction algorithm (RRA) adjusted the model to minimize dynamic inconsistencies between the experimental GRFs and body segment kinematics. Based on the adjusted model after RRA, computed muscle control (CMC)39 estimated muscletendon states and excitation patterns that reproduced the motion while minimizing the sum of excitations squared. Unmeasured handrail forces were estimated as previously described3. The dynamic consistency of the simulations was evaluated by examining the resulting residual forces and moments and verifying them with the OpenSim guidelines18.
Dependent Variable
The dependent variable for the regression analysis was RF activity. The raw RF EMG signals were filtered with a fourth-order band-pass Butterworth filter with cutoff frequencies of 20–400 Hz to remove artifacts, demeaned, rectified, and low-pass filtered with a 4th order Butterworth filter at 10 Hz. To reduce inter-individual variability, we used the mean of the EMG envelope in the trial as the normalization reference8. The processed signals were divided into gait cycles based on a paretic limb heel-strike event using vertical GRFs and then normalized into 100-time frames. A reflex response would occur within 120ms after stimulus onset through mono- or polysynaptic mechanisms31. This time interval has been employed to detect quadriceps stretch reflexes following mechanical knee perturbations during gait27. In addition to the 120ms window, we accounted for timing errors from computation of RF fiber stretch velocity and evaluated windows of 90ms and 150ms. We defined the initiation of reflex as the timing of the peak pre-swing muscle fiber stretch velocity estimated from the musculoskeletal simulation. We then numerically integrated the RF EMG envelopes following stimulus onset to obtain integrated EMG (iEMG) measures representing an involuntary response. Figure 1 visualizes RF iEMG for involuntary responses.
Gray line indicates simulated RF fiber velocity. Red solid line represents the external knee torque assistance profile by the powered knee exoskeleton device. Vertical bars show iEMG of RF. The involuntary response is captured from the peak simulated RF fiber velocity until 120 ± 30ms following (yellow).
Independent Variables
A total of 14 predictors were extracted from each gait cycle based on muscle-tendon states as well as gait kinematic and kinetic data. One muscle-tendon relevant variable was selected from the musculoskeletal simulation: the peak pre-swing RF fiber stretch velocity. Three variables from hip and knee joint kinematics in the pre-swing phase (i.e., peak hip flexion and knee flexion velocity, and relative peak knee flexion velocity to peak hip flexion) were chosen due to the biarticular physiology of the RF and its associations with post-stroke SKG15. Correspondingly, we used body kinematics which are velocities and accelerations based on the configuration (center of mass position and orientation) of each body segment. Specifically, we used 3 body kinematic variables from the thigh and shank segments (i.e., peak thigh and shank angular velocity, and relative peak shank angular velocity to peak thigh) in the pre-swing phase. Additionally, we used the pelvic segment acceleration in the anteroposterior plane as a proxy for body acceleration7. The remaining six variables were acquired from propulsive/brake impulse and ground reaction forces. Specifically, we computed the paretic limb’s propulsive, braking, and net impulse during the pre-swing phase 6. GRF signals were processed to obtain the peak anterior/posterior GRF in the pre-swing phase, peak medial/lateral GRF in the stance phase, and peak vertical GRF in the stance phase (Table 1).
Variable Selection
Variable selection is a procedure to select appropriate variables from a complete list of variables by removing those that are irrelevant or redundant. This step could be achieved by either parametric or nonparametric regression. The parametric approaches provide generalizable and interpretable statistical inference but require assumptions of linearity and normal distribution. Conversely, nonparametric machine learning methods do not need stringent distribution assumptions, but their models are less interpretable and risk overfitting without optimal tuning of hyper-parameters. In this study, we aimed to achieve a comprehensive variable selection by utilizing both parametric and nonparametric regression methods to address unknown variable relationships.
As a parametric approach, we deployed the least absolute shrinkage and selection operator (LASSO) method which is a penalized regression technique that uses an L1-norm penalty on the regression coefficients40. By the addition of the regularization term on linear regression, the LASSO regression produces a sparse and interpretable linear model with important variables automatically selected. In this study, we used glmmLasso package in R statistical software16 for modeling generalized linear mixed model using LASSO. We set a random intercept model with subjects as random effects and all predictors as fixed effects. We then built repeated 5-fold cross-validation to avoid a bias in the evaluation regression model. In data splitting, stratified sampling was applied based on subjects to obtain a sample population that best represents the entire population. In every cross-validation loop, the parameter λ was tuned by grid search using Bayesian Information Criteria (BIC). To quantify the importance of each predictor variable in regression, we ranked variables based on the selection proportions defined by frequencies of each selected variable in trained LASSO models from repeated cross-validation. The selection proportions have a score ranging from 0 (always shrunk, unimportant feature) to 1 (never shrunk, essential feature) by averaging variable selection frequencies so that we could infer individual feature importance in the regression model. The same ranking method was used in a previous regression study for post-stroke clinical outcomes25. In this study, we repeated 50 cross-validation repetitions to obtain variable selection proportions from the LASSO regression analysis. Variables above the threshold were determined as critical parameters for prediction. The threshold was 0.95 meaning 95% of the chance for selection across repeated cross-validation.
As a nonparametric approach, we deployed Bayesian Additive Regression Trees (BART), which is a nonparametric, ensemble method of regression tree models9. BART combines the advantages of the Bayesian model, incorporating past information about a parameter and forming a prior distribution for future analysis, and ensemble methods. It recently has gained popularity due to its superior prediction performance over other machine learning techniques (i.e., random forests, gradient boosting model, neural networks) in various study settings9,19. To evaluate the relative importance of each predictor variable from BART, Chipman et al. proposed the variable inclusion proportion corresponding to the proportion of times each variable is selected for splitting nodes across Markov chain Monte Carlo iterations in the sum-of-trees model9. Intuitively, variables with higher inclusion proportions are the more important variables in prediction. In this study, we implemented BART and extracted the variable inclusion proportions as a feature importance metric by using the BART package publicly available in R statistical software36. We performed a repeated 5-fold cross-validation with a total of 50 repetitions for unbiased regression outcomes. The stratified split dataset across repetitions was the same as the one used for glmmLASSO for a lateral comparison. Variables above the threshold, denoted by 1 / (total number of predictors), were determined as a crucial subset. This threshold stands for the probability of equal inclusion for all predictors in the model.
Reduced Model Selection
With key variables resulting from both parametric and nonparametric variable selection, we fit reduced models using a linear mixed and BART regression. For linear mixed model regression, we formulated a random intercept model with nested random effects in terms of subjects and exoskeletal assistance. Parametric and nonparametric models were compared to each other based on goodness-of-fit measures, such as adjusted R2 and Root Mean Square Error (RMSE) acquired by 50 repetitions of 5-fold cross-validations. Based on the reduced fit model, we additionally performed a model selection through backward stepwise elimination. This stepwise regression eliminated insignificant variables resulting in the final best fit model with minimum numbers of predictors. Additionally, the relationships between variables were examined by Spearman’s rank correlation coefficient. The overall workflow was summarized in Figure 2.
AIC, Akaike information criterion; Adj. R2, Adjusted R2.
Results
Parametric and non-parametric models showed similar results (Figure 3). Four parameters from kinematic and muscle-tendon state variables in the pre-swing phase remained (>95% of selection proportion) in the LASSO model across all 50 iterations: peak shank velocity, relative peak knee velocity with respect to hip, peak RF fiber stretch velocity, and peak knee flexion velocity. Using BART, six predictor variables from kinematic and muscle-tendon state variables in the pre-swing phase were determined as crucial variables. These were four common important variables between the parametric and nonparametric variable selection: the peak shank velocity, relative peak knee velocity with respect to the peak hip flexion, the peak RF fiber stretch velocity, and the peak knee flexion velocity. All parameters above the significance threshold using LASSO were also above the significance threshold using BART. Both variable selections resulted in low importance scores for all variables related to GRF and impulse.
The square markers represent the average, and the solid lines are the standard deviation of total 50 repetitions. The red-colored variables are essential variables showing above threshold (vertical dashed line) denoted by 95% of selection proportions and 1/total number of variables for glmmLASSO and BART, respectively. The common essential variables between methods denoted by superscript symbol of ‡.
To ensure the model was both representative and not overly complex, we fitted a reduced model using all six significant predictors from Figure 3, and then evaluated the goodness-of-fit of models. We used the same random effects and the same hyperparameter setting as the variable selection. Table 2 summarizes accuracy measures for full and reduced regression models from repeated 5-fold cross-validation. The reduced model showed slightly degraded goodness-of-fit (i.e., increased RMSE and decreased adjusted R2) compared to the full model. Except for the parametric regression model’s out-of-bag test, all accuracy measurements differed significantly between full and reduced models (p < .001).
Backward stepwise elimination additionally cut down the number of variables from the reduced model of six important predictors. The stepwise elimination on the parametric model was based on the Akaike information criterion (AIC). The adjusted R2 was used for the nonparametric model’s stepwise elimination. Stepwise regression using the parametric linear mixed model determined the last four crucial variables including two joint kinematic and another two body kinematic variables: peak hip flexion velocity, relative peak knee velocity to peak hip flexion, peak shank velocity and the relative peak shank velocity to peak thigh velocity. Table 3 represents the summary of the final regression model by a random intercept mixed model with nested random effects with respect to subject and assistance. Small marginal R2 (0.074) but moderate conditional R2 (0.505) accounted for the most variance of this model. The adjusted R2 of this final model was 0.533. Backward elimination based on the nonparametric model, BART, revealed the same last four critical predictors. Figure 4 illustrates the goodness-of-fit changes through backward elimination using BART. Once the model consisted of less than four predictors, its adjusted R2 started decreasing significantly meaning that variables in the gray shaded area of Figure 4 were the minimum crucial predictor set. The adjusted R2 of BART model with four variables was 0.638. The results were robust to differences in timing (Supplementary material).
The variable names shown in the bottom axis are the predictors removed at each step. For example, “KneeFlex_Vel” was removed at the first step of elimination from a model with six variables. The gray shaded area indicates the steps that the adjusted R2 was significantly affected by the elimination.
The relationships between selected kinematic predictors and muscle-tendon state, pre-swing rectus femoris fiber stretch velocity, were examined by Spearman’s rank correlation coefficient. Table 4 summarizes the coefficients.
Discussion
Exoskeletal assistance has the potential to unload some of the effort from clinicians and caregivers but does not yet account for other neuromuscular impairments such as hyperreflexia1,37. The aim of this study was to determine accessible biomechanical predictors of quadriceps hyperreflexia based on walking data from those with post-stroke SKG with and without knee flexion exoskeletal assistance. We used both parametric and non-parametric regression techniques to obtain two different perspectives of the key kinematic, kinetic and muscle-tendon predictors of rectus femoris hyperreflexia. We found four kinematic predictors obtained during pre-swing phase were shared between the regression techniques: peak hip flexion velocity, the relative peak knee to peak hip flexion velocity, peak shank angular velocity, and the relative peak shank to peak thigh angular velocity. These parameters can be measured by two or three wearable sensors. The implication of these findings is that these parameters can be used to predict hyperreflexia and control during exoskeletal assistance.
Spasticity is often characterized as velocity-dependent hyperexcitability of the stretch reflex23. Since it was widely believed that the spindle proprioceptive receptors encode information for muscle length and velocity changes20,26, spasticity modeling has been based on muscle length and velocity feedback24. Our group’s previous musculoskeletal simulation analysis confirmed the strong association between simulated RF muscle fiber stretch velocity and RF excitability2. In this study, both parametric and nonparametric regression techniques for variable selection resulted in the same finding that kinematic variables are predictive of an RF reflex. Indeed, there were strong correlations (ρ ≃ 0.90) between some of the selected predictors (relative peak knee flexion velocity to peak hip flexion velocity, peak shank angular velocity) and RF fiber stretch velocity (Table 4), as we expected. We did not expect that these parameters would be more predictive of a reflex than RF fiber stretch velocity, the physiological mechanism of a stretch reflex. The fiber stretch velocity was computed using generic Hill-type muscle models and mathematically driven cost functions39, rather than direct measurement, a model validated in our previous work2. Yet, it remains possible that RF fiber stretch velocity could have played a stronger statistical role if the true value could have been obtained. Regardless of the level of significance between variables, it is notable that out of all 14 parameters, the ones found to be most predictive of a reflex were ones related to hip and knee kinematics as well as RF fiber stretch velocity. This finding is consistent with our hypothesis that more easily measurable kinematic variables can serve as a proxy for fiber stretch velocity in exoskeletal control.
Over the past decade, the exoskeleton community has moved towards human-in-the-loop design43, including the incorporation of real-time estimation of musculotendon mechanics32, and ultrasound imaging13,29,35. While these approaches provide an exciting potential for the future of personalized robotics, the extra equipment remains limited to a few research institutes and adds significant complexity. In this work we found that joint kinematics accurately predict reflex responses in the RF. We estimate that kinematic predictors from this study could be measured by only two or three portable motion sensors such as inertial measurement units. This is significant because it questions the basis for using complex imaging or computational technology in lieu of simpler proxies. Enabling the use of ubiquitous sensor systems would allow exoskeleton engineers and rehabilitation researchers to design practical and effective devices.
It is possible that this work could lack a wide range of variables. Due to the exponential increase in the sampling volume, adding extra dimensionality is not always beneficial in regression analysis, specifically when there are sparse relevant predictors compared to the total number of predictors or the fundamental relationships are nonlinear (so called “curse of dimensionality”)14. To avoid these adverse effects, we limited the scope of parameter sources among previously assessed associations with post-stroke SKG. Additionally, we did not include parameters that were unable to be directly measured, such as joint kinetics, except for simulated muscle fiber stretch velocity. This omission was due to the practical difficulty in measuring real time joint kinetics. Therefore, in this work, we aimed to balance comprehensiveness, practicality and overfitting.
This study has several limitations. Regression approaches based on observations from only 8 participants may affect generalization of the results. Despite the small sample size, the repeated 5-fold cross-validation34 took into account the randomness in sampling resulting in unbiased and trustworthy variable and model selection results. This study does not provide causal relationships between kinematic parameters and quadriceps hyperreflexia. However, it is notable that the predictors of RF hyperreflexia (i.e. related to knee-hip velocities) found in this work were all physiologically related to the biarticular span of the RF suggesting that this statistical approach may help reveal mechanisms of hyperreflexia. For instance, the relative peak angular velocity between the knee and hip was the most correlated variable with RF fiber stretch velocity (ρ = 0.90 in Table 4). Still, the relative peak angular velocity between the shank and thigh had a weak negative correlation (ρ = -0.18 in Table 4) to RF fiber stretch velocity. While thigh-shank segment and knee-hip joint velocities are not equivalent, the large difference in correlations is difficult to explain mechanistically and other approaches may be needed to investigate further. Lastly, our findings were based on two types of regression techniques among numerous statistical and machine learning approaches. Although there exists an advanced penalized regression approach improving the limitations of LASSO33 (i.e., Elastic Net44), glmmLASSO in this study could result in a better model due to subjects’ random-effects modeling. There were also popular tree-based machine learning algorithms, such as gradient boosting28 and random forest methods5. It has been reported that BART has outperformed these other methods9. Therefore, the use of these two different regression methods was sufficient.
Conclusions
Despite the rapid growth of wearable robot technology, there is still a lack of consideration regarding neuromuscular impairments such as hyperreflexia into device design. In this study, we introduced regression-based variable selection for predicting rectus femoris hyperreflexia following pre-swing knee flexion exoskeletal assistance. We found that four kinematic variables relevant to knee and hip joint motions were sufficient to effectively predict hyperreflexia in people with post-stroke SKG. These results suggest that effective monitoring and control of these accessible knee and hip kinematics may be more productive at regulating hyperreflexia than the more complex acquisition of muscle fiber stretch velocity. This work represents a novel statistical approach to identifying biomechanical associations with reflex function. This information could be used to further understand the mechanisms of hyperreflexia and its associated clinical interventions.
Data Availability
All data produced in the present study are available upon reasonable request to the authors
Funding
This work was supported by the National Center for Medical Rehabilitation Research at the National Institute of Child Health and Human Development under the award number R01HD100416.
Supplementary Material
We performed a sensitivity analysis based on two additional reflex response windows by changing time windows for computing involuntary RF muscle activity: 90 ms and 150 ms. Figures A1, A2 and Table A1 summarize variable selection for 90 ms time window for a reflex response and Figures A3, A4 and Table A2 are for 150 ms. The last four key selected variables were invariant regarding the timing error of RF fiber stretch velocity estimated by OpenSim simulation.
References
