|Home | About | Journals | Submit | Contact Us | Français|
Individuals following stroke exhibit altered muscle activation and movement patterns. Improving the efficiency of gait can be facilitated by knowing which muscles are affected and how they contribute to the pathological pattern. In this paper we present an electromyographically (EMG) driven musculoskeletal model to estimate muscle forces and joint moments. Subject specific EMG for the primary ankle plantar and dorsiflexor muscles, and joint kinematics during walking for four subjects following stroke were used as inputs to the model to predict ankle joint moments during stance. The model’s ability to predict the joint moment was evaluated by comparing the model output with the moment computed using inverse dynamics. The model did predict the ankle moment with acceptable accuracy, exhibiting an average R2 value ranging between 0.87 and 0.92, with RMS errors between 9.7% and 14.7%. The values are in line with previous results for healthy subjects, suggesting that EMG-driven modeling in this population of patients is feasible. It is our hope that such models can provide clinical insight into developing more effective rehabilitation therapies and to assess the effects of an intervention.
Abnormal movement patterns are associated with disability following stroke [1, 2]. Knowledge of internal forces and moments during movement is important for developing better rehabilitation programs for this population. Unfortunately, in vivo muscle force measurements are invasive and therefore not practical. For this reason, optimization techniques have been employed to predict muscle forces using a variety of cost functions [3, 4]. Other studies using forward dynamics simulations optimize muscle excitation (i.e., activation) patterns to reproduce recorded movement kinematics and to estimate muscle forces. In these studies, the cost function was the error between measured and simulated joint kinematics and ground reaction forces [5, 6]. These cost functions may have difficulty predicting co-contraction and cannot account for different muscle synergy patterns during various tasks . EMG-driven biomechanical models are not subject to the same limitations since the models use recorded EMG to determine muscle force [7–11]. That is, these models use EMG as inputs rather than attempting to predict how muscles are activated to produce a given movement.
Although EMG-driven models have been used to investigate movement dynamics for healthy individuals during a variety of tasks, they have not been used to model muscle forcesproduced by people following stroke. Post-stroke patients often have difficulty generating an effective push-off during late stance and they exhibit abnormal activation patterns . Abnormal and more variable muscle activity presents a challenging paradigm to test the robustness of our EMG-driven musculoskeletal model.
In this paper we describe our EMG-driven model and how it is implemented to estimate muscle forces and joint moments. Of specific interest was: (1) to provide a comprehensive description of the model structure, and (2) evaluate the model’s ability to predict the ankle joint moment for patients following stroke. The ability to predict joint dynamics and muscle forces in patients following stroke would benefit treatment planning and evaluation of therapies geared at improving patient function.
Our EMG-driven model uses subject recorded joint kinematics and EMGs as inputs to estimate individual muscle forces and joint moments [10, 11]. The EMG-driven model is comprised of three separate but related sub-models: (1) anatomical, (2) muscle activation dynamics, and (3) muscle contraction dynamics. Sub-models (2) and (3) require the specification of several muscle specific parameters. These parameters characterize the force generating potential of each muscle and define how the recorded EMGs are transformed to muscle activations. Choice of parameters greatly affects the model’s ability to predict joint moments. Stochastic optimizationis used to determine ideal values for the muscle parameters during a model-tuning process which is described in the section entitled “Model Tuning”.
The lower limb anatomical model was developed using SIMM (Motion Analysis Corp.) , and extended by Lloyd and Buchanan . This model includes the four major plantar and dorsiflexor muscles: tibialis anterior (TA), lateral gastrocnemius (LG), medial gastrocnemius (MG) and soleus (Sol). The anatomical model was used to determine individual muscle tendon lengths and moment arms for the primary ankle plantar and dorsiflexors.
The muscle activation dynamics model transforms raw EMG to muscle activation. Firstly, the raw EMG were high-pass filtered to remove movement artifact, then full wave rectified, normalized by peak rectified EMG values obtained during maximum voluntary contraction (MVC), and then low-pass filtered to account for the natural filtering property of muscle . A recursive filter was then used to determine neural activation, u(t), from the rectified, normalized and filtered EMG input, e(t). This process characterized the time varying features of EMG and was represented in the following discrete form [14, 15]:
where d is the electromechanical delay, and α, β1, and β2 are coefficients that define the second-order dynamics. To realize a positive stable solution of Equation (1), the following must hold true:
where |γ1| < 1, |γ2| < 1, as discussed by Buchanan et al. .
Woods and Bigland-Ritchie  reported a nonlinear force-EMG relationship at low levels of force for some muscles and a linear relationship at higher levels. To account for this, we used a logarithmic function for low levels of EMG and a linear function for high values, as shown below:
where u(t) is the neural activation and a(t) is the muscle activation. The coefficients c, d, m, and b can be solved and reduced to a single parameter, A, to characterize the curvature of the relationship . A varies from 0.0 to approximately 0.12.
The muscle activations and joint kinematics were used as inputs to a modified Hill-type muscle model to calculate individual muscle forces. The muscle-tendon unit was modeled as a muscle fiber in series with a tendon. The muscle fiber had a contractile element in parallel with an elastic element and a damping element (Figure 1). Pennation angle (3) is the angle between the lines of action of the tendons and the muscle fiber. lm is muscle fiber length, lt is total length of the tendons (each one is half of the length), and lmt is musculotendon length. Muscle-tendon force (F) was calculated as below:
Where, Ft is tendon force, and Fm is muscle fiber force that is comprised of the three parallel elements. Fmax is maximum isometric muscle force, m is normalized muscle fiber length, m is normalized muscle fiber velocity, a(t) is muscle activation. A (m) represents the active force-length relationship , V (m) represents the force-velocity relationship [19–21], P (m) represents the passive elastic force-length relationship , and bm is the damping factor . These parameters are normalized to maximum isometric muscle force, optimal fiber length (lom) and maximum muscle contraction velocity (vmax).
Pennation angle ((t)) was calculated as described by Scott and Winter , and is shown below:
where lm(t) is the muscle fiber length at time t, and o is the pennation angle at muscle optimal fiber length .
The coupling relationship between optimal fiber length and muscle activation (i.e., optimal fiber length increases with decreasing activation) reported by Guimaraes et al.  and Huijing  is implemented in our model using the mathematical description of Lloyd and Besier :
where λ is the percentage change in optimal fiber length, is the optimal fiber length at time t.
Tendon force varies with tendon strain ε only when tendon length lt is greater than tendon slack length ; otherwise the tendon force is zero as described by Zajac . Subsequently the normalized tendon force, t, is given by
where tendon strain was defined as
Muscle tendon length and activation data were used as input to the muscle model, and then muscle fiber lengths were calculated by forward integration of the fiber velocities using a Runge-Kutta–Fehlberg algorithm. After the muscle fiber lengths and velocities were determined, the muscle force was calculated using Equation (4). Once individual muscle forces were computed, they were multiplied by their respective muscle moment arms derived from the anatomical model to calculate the net joint moment.
In our muscle activation dynamics model, γ1, γ2, d, A are four parameters that account for the relationship between EMG and muscle activation. In our muscle contraction dynamics model, optimal fiber length (lom), tendon slack length ( ) are anatomical parameters for each muscle, and λ is the global percentage change in lom as a function of muscle activation. In addition, to account for differences in muscle strength due to variations in PCSA between people, we scale the maximum isometric force, Fmax, for the flexor and extensor muscles . Because these parameters are likely subject-specific, we begin by assigning each parameter an average value obtained from the literature. Initial values are allowed to vary within physiolgoical limits based on data reported in the literature Yamaguchi et al. .
Our calibration process uses a hybrid forward-inverse dynamics approach, which assumesthe joint moment estimated using the EMG-driven model should match the moment computed using inverse dynamics. A simulated parallel annealing algorithm  was used during the calibration process to tune the model parameters to minimize the difference between the EMG-driven model joint moment and the moment computed using inverse dynamics. The initial values of and lom for each muscle were constrained to change within ±15% and ±5% of their initial values respectively. Ranges for γ1, γ2 and A were described earlier; the range for electromechanical delay, d, was between 10 – 100 ms . Maximum isometric force was allowed to vary between 50% and 200% of the initial value for Fmax. It is important to note that Fmax for the LG, MG and Sol were constrained to scale by the same amount to prevent unrealistic scenarios in which plantarflexion is dominated by a single muscle. Initial values for γ1, γ2, d, A and Fmax were set at the midpoint of their upper and lower limits. Once the model was tuned by adjusting the parameters, it was used to predict the ankle moment using EMG and kinematics for other trials on which it was not tuned.
Four stroke patients who gave informed consent were included in this experiment (Table 1). The subjects walked at self-selected speeds. Motion analysis and force plate data were used to compute joint kinematics and moments for four walking trials. In addition, EMG was recorded from the tibialis anterior (TA), medial gastrocnemius (MG), lateral gastrocnemius (LG) and soleus (Sol) muscles. Maximum voluntary isometric contraction (MVIC) trials were collected fornormalization of the EMG. data. The raw EMG were first high-pass filtered using a forward and reverse pass 4th order Butterworth filter (cutoff frequency 50 Hz), then full wave rectified, normalized by peak rectified EMG values obtained from MVC trials, and then low-pass filtered using a forward and reverse pass 4th-order Butterworth filter (cutoff frequency 4 Hz).
The EMG-driven model was tuned to each of the four walking trial for each subject. The coefficient of determination (R2) between the model (i.e., EMG-driven) and calculated (i.e., Inverse dynamics) and the normalized root mean square (RMS) error (normalized to peak-to-peak joint moment) were calculated for each of the tuning trials. Once the subject’s muscle parameters were tuned to a particular trial, the model was given EMG and joint kinematics for the subject’s other three walking trials and the tuned-model was used to predict the ankle moment. For example, the model was tuned to trial 1 and then used to predict trials 2, 3 and 4. Similarly, the model was tuned to trial 2 and used to predict trials 1, 3 and 4.
Gait in people who have had strokes is markedly different than in healthy subjects. For impaired subjects in this study, circumduction was present, likely to compensate for dorsiflexion weakness and prevent toe-stubbing during swing. The EMG patterns were also notably different than normal . For example, the triceps surae (gastrocnemii and soleus) were active during initial contact, a time when the plantarflexors are normally silent. In addition, an increase in tibialis anterior was noted following mid-stance, with peak values occurring at roughly the same time as triceps surae during late stance (Figure 2). In addition, there was noticable variability in the activation profiles.
Despite these differences, the model was able to be tuned to the subjects to capture the ankle joint moments accurately (Figure 3). When average values obtained from the literature were used to characterize the muscle parameters, poor agreement was observed between the model output and the moment computed using inverse dynamics (i.e., measured). However, once the model parameters were tuned to fit the subjects, the model output was in close agreement to the measured moment. The R2 value for all the four subjects was high (Table 2), indicating the model did an excellent job of following the shape of the joint moment computed using inverse dynamics. Excellent agreement in shape and similarity in magnitude and times to peak were observed when comparing ankle moments estimated from the model with the experimentally derived values (Figure 3). The normalized RMS error for the tuning trials ranged between 1.9% – 6.9%. This improvement implies that initial values for the muscle parameters obtained from the literature were clearly not appropriate for the subject.
The model also did well in predicting joint moments (without re-tuning) given novel data. Once the muscle parameters were tuned to a particular subject for one trial, the model was then used to predict the ankle moment using EMG and kinematics for a different trial for that subject (Figure 4). Note that different activation patterns were used for the trial to which the model wastuned and the muscle activity for the trial that was predicted (Figure 2). In addition, notice the difference in magnitude between the trial used to tune the model (Figure 3) and the ankle moment that was predicted (Figure 4). Normalized RMS error and R2 for the predicted trial was 5.8% and 0.95 respectively. The averageR2 and normalized RMS errors ranged between 0.87 – 0.92, and 9.7% – 14.7% respectively.
The objectives of this paper were to describe in detail our EMG-driven model and discuss how it is implemented to predict joint moments and muscle forces. The second goal of this paper was to evaluate the model’s ability to predict joint moments in patients following stroke. Joint moments and muscle activation patterns are abnormal in this population and therefore present a challenging paradigm with which to test our EMG-driven model. Model results were verified by comparing the predicted and inverse dynamic joint moments during novel walking trials.
Overall, the model did predict the ankle joint moment for the four subjects in this study, despite the variability in muscle activation patterns and joint moments between trials. This may be associated with dysfunction caused by spasticity, reduced coordinative control, muscle weakness, and impaired sensory-motor control [12, 32, 33]. Results for the four patients were comparable to data reported for healthy individuals. Bogey et al.  reported calibration results on the ankle of healthy subjects during walking with a joint moment R2 of 0.97 and a normalized RMS error of 4.8%, but they did not use the calibrated model to predict novel trials. Buchanan et al.  had success predicting ankle joint moments with an R2 of 0.94 and a normalized RMS error of approximately 9% for a healthy subject during gait. In another study, Manal and colleagues reported R2 values ranging between 0.87 and 0.95 for three healthy subjects during stance . Taken together, these studies and our results suggest that EMG-driven modeling can be used to estimate muscle forces and joint kinetics in patients following stroke.
Lloyd and Besier used an EMG-driven model to estimate knee moments and muscle forces in healthy subjects during dynamic tasks . The model was tuned using a variety of dynamic activities, including jogging and cutting trials. This was done to expose the model to a wide range of muscle activiation patterns and differing joint moment profiles. The model was able to predict the knee joint moment for a range of tasks with good accuracy. This is a more robust test of the model than simply tuning it to a walking trial and then predicting another walking trial for the same subject. For example, Kadaba et al. have shown that joint kinematics, kinetics and muscle activation patterns are similar within subject from trial to trial . Thus, once a model is tuned to a walking trial for a healthy subject, the model will likely do well predicting another walking trial since the input data (i.e., EMG and kinematics) will be similar to the data used to tune the model.
Neurologically impaired subjects become fatigued easily and therefore testing them under a variety of experimental conditions is not feasible. In our study we tuned the model to a walking trial, and once tuned, predicted other walking trials for the subjects. This paradigm is a good test of the model because the EMG and moment data were variable from trial to trial, unlike data for healthy subjects.
Since our EMG-driven model can be used to predict novel trials, it may be possible to determine how muscle activation patterns for a stroke patient would have to be augmented to achieve a desired moment profile. These calculated corrective changes in muscle activation patterns may be used as reference to derive appropriate stimulation patterns, which can be applied in stroke patients’ gait training with the use of functional electrical stimulation.
There are several assumptions and limitations of this study. First, although model parameters are tuned for each subject, the model does not account for differences in musculoskeletal size. Second, hemiparetic muscles of post-stroke patients have been shown to have altered fiber type proportions  and fiber atrophy . Such changes of muscle architecture may influence the muscle force-activation relationship of these patients. Since such changes are difficult to model, we included adjustable parameters, such as maximal muscle strength, optimal fiber length and tendon slack length. These parameters were adjusted to account for the impaired musculature of each subject, and it was verified that our calibrated models could be used to predict novel trials.
In conclusion, our EMG-driven model did predict the ankle joint moment during stance for four patients following stroke. The ability to accurately predict the joint moment is a first step in estimating muscle forces and joint loads associated with pathological gait. We are continuing to refine this model for the study of stroke and other pathologies (e.g., osteoarthritis). The long term goal of this work is to apply such models and provide clinicians with insight regarding patient specific muscle function and joint loading.
This work supported, in part, by NIH R01-HD38582 and NIH R01-AR46386. We wish to thank Daniel Benoit, Wei-Li Hsu and Trisha Kesar for assistance with the data collection.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.