Muscle Synergies Facilitate Computational Prediction of Subject-Specific Walking Motions

Introduction

Roughly one in six people worldwide will suffer a stroke at some point in their lifetime, with ~15 million people experiencing a stroke each year (World Stroke Organization, 2016). Due to improvements in acute stroke management, the majority of these individuals will survive their initial stroke, which has helped make stroke a leading cause of serious, long-term disability in adults worldwide (Go et al., 2013; World Stroke Organization, 2016). More than a third of stroke survivors experience significant physical disability (Lloyd-Jones et al., 2010), with walking dysfunction being among the greatest stroke-related limitations contributing to disability. While the majority of persons who suffer a stroke regain some level of ambulatory function, their gait is typically slow, asymmetrical, and metabolically inefficient (Olney et al., 1986; Roth et al., 1997). Diminished walking ability is tied to decreased quality of life, increased risk of depression, and increased risk of serious secondary health conditions (Blair et al., 1989; Mutikainen et al., 2011; Ostir et al., 2013). Restoration of walking function following a stroke is therefore both a high priority for rehabilitation and an important public health problem.

Despite recognition of the problem, current clinic-based neurorehabilitation methods produce only modest improvements in walking function for persons post-stroke (States et al., 2009; Bogey and Hornby, 2014; Winstein et al., 2016). For this reason, researchers and clinicians have explored a variety of neurorehabilitation approaches in search of an effective means to restore post-stroke walking function. These approaches include functional electrical stimulation (FES) (Popovic et al., 1999; Kesar et al., 2009, 2010; Sabut et al., 2013; Chung et al., 2014; O’Dell et al., 2014; Pilkar et al., 2014; Auchstaetter et al., 2015; Chantraine et al., 2016), ankle–foot orthoses (AFOs) (Ferreira et al., 2013; Tyson et al., 2013; Kobayashi et al., 2016), exoskeletons (Nilsson et al., 2014; Bortole et al., 2015; Buesing et al., 2015), partial body weight support (Ng et al., 2008; Lee et al., 2013) and split-belt treadmill training systems (Reisman et al., 2007; Malone and Bastian, 2014), and robotic gait trainers (Pennycott et al., 2012; Mehrholz et al., 2013; Bae et al., 2014; Hussain, 2014; Dundar et al., 2015). Each of these approaches has shown varying levels of promise for improving post-stroke walking function. However, a critical challenge is determining the treatment prescription – which approach to apply, how much of the approach to apply, and how the approach should be applied – that will maximize recovery of walking function for any particular individual. Furthermore, there is currently no way to identify whether a small amount of treatment provided by a combination of approaches might be dramatically more effective than a large amount of treatment provided by a single approach (Belda-Lois et al., 2011). Current treatment design methods based on trial-and-error and subjective clinical judgment cannot address these challenges, since they do not provide an objective means for predicting a patient’s walking function following a specified treatment or treatment combination.

One possible approach for overcoming this challenge is to use patient-specific neuromusculoskeletal models to predict post-treatment walking function for different neurorehabilitation technologies (alone or combined) under consideration. Such models should account for how the patient interacts with the treatment approach (Mooney and Herr, 2016) so that the optimal prescription can be recommended based on objective predictions of walking improvement. A number of studies have pursued such modeling efforts by simulating the effects of FES (Riener, 1999; Heilman and Kirsch, 2003; Zhang and Zhu, 2007; Shao and Buchanan, 2008; Nekoukar and Erfanian, 2013; Sharma et al., 2014; Alibeji et al., 2015), exoskeletons (Fleischer and Hommel, 2008; Afschrift et al., 2014; Farris et al., 2014; Sawicki and Khan, 2015), orthoses (Zmitrewicz et al., 2007; Crabtree and Higginson, 2009; Silverman and Neptune, 2012), and strength training (Goldberg and Neptune, 2007; Knarr et al., 2014) on lower extremity function and walking ability in the context of stroke, spinal cord injury, amputee, and general rehabilitation. Few of these studies focused on stroke (Goldberg and Neptune, 2007; Shao and Buchanan, 2008; Knarr et al., 2014), few used three-dimensional models (Fleischer and Hommel, 2008; Afschrift et al., 2014; Farris et al., 2014; Knarr et al., 2014; Sawicki and Khan, 2015), few used subject-specific models created by calibrating critical neuromusculoskeletal model parameters to movement data collected from an individual (Fleischer and Hommel, 2008; Shao and Buchanan, 2008; Knarr et al., 2014), and only one included modeling elements that accounted for subject-specific neural control capabilities and limitations (Alibeji et al., 2015). No study to date has predicted a stroke patient’s complete post-treatment walking motion and speed resulting from application of a specific neurorehabilitation intervention.

As a first step toward optimizing patient interaction with stroke neurorehabilitation technologies, this study describes the development and evaluation of a subject-specific synergy-controlled neuromusculoskeletal simulation framework that can predict three-dimensional walking motions for an individual post-stroke. The main question we address is whether actuating a subject-specific neuromusculoskeletal model with muscle synergy controls (5 per leg) facilitates generation of accurate walking predictions compared to actuating the model with muscle activation controls (35 per leg) or joint torque controls (5 per leg). We hypothesize that synergy controls will work the best since they combine a low number of control signals with a subject-specific representation of the coupling between muscle activations within each leg. We collect gait data from a stroke subject walking at 0.4, 0.5, 0.6, 0.7, and 0.8 m/s on an instrumented treadmill and use data from his self-selected speeds of 0.4–0.6 m/s to develop a subject-specific neuromusculoskeletal model. We incorporate the subject-specific full-body model into a direct collocation optimal control framework to predict new walking motions for the subject. To evaluate the framework and the potential benefits of using synergy controls, we predict how the individual will walk (including cadence and stride length) at 0.5 and 0.8 m/s (conditions for which experimental data are available for comparison) using joint torque, muscle activation, or muscle synergy controls and at 1.1 m/s (a condition for which no experimental data are available) using only synergy controls. With future simulation of different neurorehabilitation approaches, our subject-specific synergy-controlled neuromusculoskeletal simulation framework may help identify optimal neurorehabilitation prescriptions that maximize recovery of walking function on an individual patient basis.

Methods

Experimental Data Collection

To assist with development and evaluation of our subject-specific synergy-controlled neuromusculoskeletal simulation framework, we collected experimental walking data from one high-functioning hemiparetic male individual with chronic stroke-related walking dysfunction (age 79 years, LE Fugl-Meyer Motor Assessment 32/34 pts, right-sided hemiparesis, height 1.7 m, mass 80.5 kg). All study procedures were approved by the University of Florida Health Science Center Institutional Review Board (IRB-01) and the Malcom Randall VA Medical Center Research and Development Committee and included approval to study individuals with stroke-related disability. Study personnel obtained written informed consent prior to participant enrollment and involvement in study activities. Study procedures were conducted in accordance with the Declaration of Helsinki. Motion capture (Vicon Corp., Oxford, UK), ground reaction (Bertec Corp., Columbus, OH, USA), and electromyography (EMG) data (Motion Lab Systems, Baton Rouge, LA, USA) were collected simultaneously while the participant walked on a split-belt instrumented treadmill (Bertec Corp., Columbus, OH, USA) at speeds ranging from 0.4 to 0.8 m/s in increments of 0.1 m/s. 0.8 m/s was the fastest speed at which the subject could walk safely without assistance. This range of speeds included the participant’s self-selected walking speed of 0.5 m/s. More than 50 gait cycles were recorded at each walking speed. A static standing trial was collected for model scaling purposes. To facilitate subsequent creation of subject-specific foot–ground contact models, the participant wore Adidas Samba Classic sneakers, which have a flat sole and neutral midsole with no cushioning, and we collected additional static trials where we used a marker wand to trace the outline of each sneaker sole on the ground. Motion capture data were obtained using a modified Cleveland Clinic full-body marker set with additional markers added to the feet (Reinbolt et al., 2005). Marker motion and ground reaction data were filtered at a variable cut-off frequency of 7/tf Hz, where tf is the period of the gait cycle being processed, using a fourth-order zero phase lag Butterworth filter (Hug, 2011). This variable cut-off frequency would cause data collected at a normal walking speed to be filtered at ~6 Hz.

Electromyography data were collected from 16 muscles in each leg and processed using standard methods (Lloyd and Besier, 2003). A combination of surface and fine-wire EMG electrodes was used. Surface EMG data were collected for gluteus maximus and medius, semimembranosus, biceps femoris long head, rectus femoris, vastus medialis and lateralis, medial gastrocnemius, tibialis anterior, peroneus longus, and soleus. Fine-wire EMG data were collected for adductor longus, iliopsoas, tibialis posterior, extensor digitorum longus, and flexor digitorum longus. All EMG data were high-pass filtered at 40 Hz (Lloyd and Besier, 2003), demeaned, rectified, and then low-pass filtered at a variable cut-off frequency 3.5/tf Hz. Filtering was performed using a fourth-order zero phase lag Butterworth filter. EMG data from each muscle were normalized to the maximum value over all trials and resampled to 101 time points per gait cycle (heel strike to heel strike for the less involved left side) while keeping an additional 20 time points before the start of the cycle to permit modeling of electromechanical delay. In addition, each processed EMG signal was offset on a cycle-by-cycle basis so that the minimum value was zero.

Neuromusculoskeletal Model Development

The subject-specific neuromusculoskeletal model that served as the foundation for our simulation framework incorporated four modeling components to account for the unique neurophysiological and musculoskeletal characteristics of the subject: (1) a subject-specific lower-body kinematic model to simulate the subject’s skeletal motion, (2) subject-specific foot–ground contact models to simulate how the subject’s feet interact with the ground, (3) subject-specific EMG-driven muscle moment models to simulate the subject’s lower extremity joint moments, and (4) a subject-specific muscle synergy control model to simulate the subject’s neural control system. Below we describe each of these four components in further detail. Unless otherwise noted, we calibrated model parameters in each component using a single representative walking trial collected at the subject’s self-selected speed of 0.5 m/s.

Subject-Specific Lower-Body Kinematic Model

Our neuromusculoskeletal model creation process started with a generic full-body musculoskeletal model (Arnold et al., 2010; Hamner et al., 2010) developed in OpenSim (Delp et al., 2007). The original model included 37 degrees of freedom (DOFs) and 44 Hill-type muscle-tendon actuators in each leg. We locked the wrist and forearm pronation–supination angles to anatomically reasonable values for walking, leaving the following 31 DOFs: 6 DOF ground-to-pelvis joint, 3 DOF hip joints, 1 DOF knee joints, 1 DOF ankle joints, 1 DOF subtalar joints, 1 DOF toe joints connecting rear foot and toe segments, 3 DOF back joint, 3 DOF shoulder joints, and 1 DOF elbow joints. We also eliminated nine muscle-tendon actuators without related EMG data (extensor hallucis longus, flexor hallucis longus, gemelli, gracilis, pectineus, piriformis, quadratus femoris, sartorius, tensor fascia latae), leaving 35 muscles per leg that actuated hip flexion-extension, hip adduction-abduction, knee flexion-extension, ankle flexion-extension, and ankle inversion–eversion on each leg. We then scaled the modified model using the standing static trial marker data and the OpenSim “Scale Model” tool, where distances between markers placed over identifiable landmarks were averaged between the two sides to maintain bilateral symmetry following scaling.

Once the model was scaled, we calibrated joint and marker positions in the torso, pelvis, and lower-body portion of the OpenSim model using marker data from a representative walking trial. The calibration approach was similar to one described previously (Reinbolt et al., 2005, 2008) except that it was performed on the scaled OpenSim model using the OpenSim-MATLAB Application Programing Interface and included modifications to maintain correct bone geometry positions within the body segments (Charlton et al., 2004). To facilitate the calibration process, we created marker plates on the torso, pelvis, thighs, and shanks to which all markers on the respective OpenSim body segments were attached. To perform the actual calibration, we used non-linear least squares optimization (lsqnonlin) in MATLAB to adjust joint (knee, ankle, and subtalar in both legs) and marker plate (torso, pelvis, thighs, and shanks) positions and orientations in their respective body segments such that marker errors from repeated OpenSim inverse kinematic analyses were minimized. The optimization cost function included penalty terms that prevented large changes in joint and marker plate positions and orientations that would produce only small improvements in marker tracking. Modification of the two hip joint center locations was achieved by modifying the position and orientation of the rigid marker plate on the pelvis. For joint centers and orientations, symmetry between left and right sides of the body was maintained during the kinematic calibration process. Markers on the feet were not adjusted since their locations were well defined. The position and orientation of the toe axis in each foot and of the back, shoulder, and elbow joints was maintained from the scaled OpenSim model.

Subject-Specific Foot–Ground Contact Models

Following kinematic calibration, we created a subject-specific foot–ground contact model for each foot of the OpenSim model using recently developed methods (Jackson et al., 2016). The elastic foundation contact models were developed in MATLAB and used a grid of contact elements that spanned the rear foot and toes segments of each foot. To create the element grid, we started with the shoe outlines obtained from the static trial marker data and used principal component analysis to identify the principal axes of each foot (rear foot and toes segments together). Using these axes, we constructed an 11 (anterior-posterior) × 8 (medial-lateral) grid of rectangular contact elements for the left foot, where the edges of the grid extended 2.5 mm beyond the edge of the shoe outline in both directions. Forty-seven elements whose centers fell within the shoe outline were retained in the contact model, while 41 elements whose centers fell outside the shoe outline were removed. Given the locations of the MATLAB contact element centers relative to the foot markers from the static left shoe outline trial, we calculated the locations of the element centers on the OpenSim rear foot and toes segments. We then projected the left toes axis of the OpenSim model onto the contact element grid. Elements whose centers were posterior to the axis were assigned to the rear foot segment, while elements whose centers were anterior to the axis were assigned to the toes segment. The complete MATLAB/OpenSim contact element grid for the left foot was mirrored to the right foot by aligning the principal axes of the mirrored grid with those of the right foot.

Each contact element in the foot–ground contact models generated normal force using a linear spring with non-linear damping and shear force using a continuous stick-slip friction model. For any contact element i, the required time-varying inputs for contact force calculations performed in MATLAB were the penetration into the floor y_i, the normal penetration rate ${\dot{y}}_i$, and the shear slip rate $v_{s_i}$ of the element center in the laboratory coordinate system as calculated by OpenSim. The normal contact force F_i for element i was calculated as (Hunt and Crossley, 1975)

$F_i=k_i(y_i-y_0)(1+c{\dot{y}}_i)$ (1)

where k_i is the spring stiffness unique to each spring, y₀ is the spring resting length common to all springs on the same foot and essentially adjusts the height of the floor, and c (= 1e–2) is a non-linear damping coefficient common to all springs. The linear spring also generates a small amount of force when the foot is off the floor, and this negligible force transitions in a smooth and continuous way to the large force produced when the spring is in contact with the ground (Anderson and Pandy, 2001; Ackermann and van den Bogert, 2010). The non-linear damping ensures that the normal contact force does not exhibit damping-related discontinuities when a spring enters or leaves contact. The shear contact force f_i for element i was calculated using a simple continuous and differentiable friction model (Ackermann and van den Bogert, 2010)

$f_i=F_i\left[\mu \tanh \left(\frac{v_{s_i}}{v_l}\right)\right]$ (2)

where μ [= 1 (Ackermann and van den Bogert, 2010)] is a dynamic friction coefficient common to all springs and v_l (= 5 cm/s) is a latching speed common to all springs that defines the edge of a linear transition region between zero slip rate and the start of dynamic friction. Shear contact force f_i was applied to the element center in the direction opposite to the slip velocity vector. The direction calculation included a small constant value of 1e−4 in the denominator to avoid division by a small number when the slip speed was near zero. Once F_i and f_i were calculated for each contact element, the net contact force and torque due to all contact elements in the rear foot segment were calculated with respect to the rear foot origin, and similarly for the toes segment using the toes origin (Kane and Levinson, 1985). These net contact forces and torques were then applied to their corresponding segments in the OpenSim model. This approach allowed rear foot and toes contributions to total ground reaction force to be predicted by the model.

We calibrated the spring stiffness k_i of each contact element in both feet and the spring resting length y₀ for all contact elements in each foot using marker and ground reaction data from a representative walking trial. We made two assumptions about the spring stiffness distribution across the bottom of the shoe to simplify the calibration process. First, we assumed that the mirrored stiffness distribution was the same for both feet. Second, we assumed that the stiffness distribution across the entire shoe bottom could be approximated by a three-dimensional parabolic surface, which possesses only six unknown coefficients rather than 47 unknown independent spring stiffness values and prevents neighboring springs from having dramatically different stiffnesses. To calibrate these six coefficients and the two resting lengths, we formulated a direct collocation optimal control problem that tracked experimental marker, ground reaction, and inverse dynamic joint torque data for the entire body with higher weight placed on matching marker position and ground reaction data for the two feet. Tracked ground reaction quantities for each foot included the three ground reaction force components and three ground reaction torque components calculated about the midfoot marker projected onto the floor (see Optimal Control Walking Predictions below for further details).

Subject-Specific EMG-Driven Muscle Moment Models

To generate subject-specific joint moments from predicted muscle activations, we calibrated lower extremity EMG-driven muscle moment models for both legs to a large number of walking trials collected from the subject. Complete details of our EMG-driven model calibration process, and a full assessment of its ability to predict joint moments accurately for the same subject walking at multiple speeds, can be found in Meyer et al. (2016). In brief, the model calibration process used experimental walking data collected from the subject at 0.4, 0.5, and 0.6 m/s, bracketing his self-selected speed of 0.5 m/s. Ten trials from each speed were used for calibration. We adjusted three types of model parameter values in our calibration process: (1) EMG-to-activation parameter values, (2) Hill-type muscle-tendon model parameter values, and (3) surrogate musculoskeletal geometry parameter values. Below we describe each category of adjusted model parameter values in greater detail.

For our EMG-to-activation model, parameter values adjusted during calibration included electromechanical delays, EMG scale factors, activation time constants, and muscle non-linearity constants. A single electromechanical time delay between 0 and 100 ms was used for all muscles in the same leg, allowing the two legs to have different electromechanical delays. An EMG scale factor was found for each muscle in each leg, resulting in 70 different scale factors. An activation time constant for a first-order activation dynamics model (He et al., 1991) was found for each muscle, with time constants assumed to be identical for the same muscle in both legs, resulting in 35 different time constants. Deactivation time constants were assumed to be four times larger than corresponding activation time constants (Thelen, 2003; De Groote et al., 2012; Millard et al., 2013). Finally, a non-linear constant defining the conversion of neural activation to muscle activation was found for each muscle (Buchanan et al., 2004), with non-linear constants assumed to be identical for the same muscle in the two legs, resulting in 35 different non-linear constants.

For our Hill-type muscle-tendon model, parameter values adjusted during calibration included optimal muscle fiber lengths and tendon slack lengths. We used a custom Hill-type muscle-tendon model with a rigid tendon, as a recent study has shown that use of a compliant tendon model for simulations of walking has little effect on predicted muscle activations and forces (De Groote et al., 2016). This model was implemented in MATLAB to facilitate customization of model properties. Initial Hill-type model parameter values (optimal muscle fiber length, tendon slack length, pennation angle) for each muscle were taken from the literature (Arnold et al., 2010) and assumed to be the same for both legs. Optimal fiber length and tendon slack length values were pre-calibrated to reproduce passive hip, knee, and ankle flexion moment data reported in the literature (Silder et al., 2007). Peak isometric force for each muscle was defined using regression equations for muscle volume reported in the literature (Handsfield et al., 2014) along with each muscle’s optimal fiber length and a maximum muscle stress of 61 N/cm² (Arnold et al., 2010). Maximum shortening velocity for each muscle was defined to be 10 optimal fiber lengths per second.

For our surrogate musculoskeletal geometry, parameter values adjusted during calibration included coefficients of polynomial functions defining muscle-tendon lengths as a function of spanned joint angles. For each muscle-tendon actuator in our kinematically-calibrated OpenSim model, we first created a surrogate model of muscle-tendon length using a cubic polynomial function of all spanned joint angles. Some muscles required a cubic function of only one joint angle (e.g., vastus medialis), while other muscles required a cubic function of multiple joint angles (e.g., gluteus maximus and gastrocnemius medialis). We then created corresponding surrogate models of muscle-tendon velocity and moment arms by defining muscle-tendon velocity as the first derivative of the muscle-tendon length polynomial with respect to time and each muscle moment arm as the negative of the first derivative of the muscle-tendon length polynomial with respect to the corresponding spanned joint angle (An et al., 1984). In this way, the polynomial functions defining muscle-tendon lengths, velocities, and moment arms shared common coefficients (Menegaldo et al., 2004; Sartori et al., 2012a).

To generate an initial polynomial fit for each muscle, we sampled muscle-tendon lengths and moment arms from our kinematically calibrated OpenSim model using a wide range of lower extremity joint angle combinations. The maximum and minimum value of each joint angle were allowed to go well beyond the values achieved by the subject during walking. Sampling was performed using 1000 different model poses specified using a Latin hypercube design, with muscle-tendon lengths and moment arms being calculated by an OpenSim “Muscle Analysis.” After outliers were removed corresponding to situations with muscle wrapping problems, we used linear least squares regression to fit muscle-tendon length and moment arms for each muscle as a cubic polynomial function of spanned joint angles.

We used sequential quadratic programing (SQP) optimization (fmincon) in MATLAB to adjust the parameter values described above such that EMG-driven models for both legs matched lower extremity inverse dynamic and passive joint moment curves as closely as possible. Inputs to the optimization were the subject’s processed experimental EMG data, joint kinematics from OpenSim “Inverse Kinematics” analyses, and joint moments from OpenSim “Inverse Dynamics” analyses from 30 selected walking trials, along with the published passive joint moment data described earlier (Silder et al., 2007). For modeled muscles without experimental EMG data (e.g., vastus intermedius) or with multiple compartments (e.g., gluteus maximus), EMG data from anatomically related muscles were used but with a separate scale factor for each muscle/compartment (Sartori et al., 2012b). Outputs were the model parameter values and predicted inverse dynamic and passive joint moments. The subject’s hip internal–external rotation moment from walking was not included in the calibration process, since EMG data were not collected from primary hip external rotator muscles. Thus, the EMG-driven model for each leg matched five inverse dynamic moments (hip flexion–extension, hip adduction–abduction, knee flexion–extension, ankle flexion–extension, and ankle inversion–eversion). During calibration, optimal muscle fiber length and tendon slack length values were allowed to vary within 25% of the values produced by pre-calibration, while changes in surrogate musculoskeletal geometry were strongly penalized so that such changes would be made only if they resulted in significant improvements in joint moment matching. The calibrated EMG-driven model for each leg was verified by predicting lower extremity joint moments for walking trials withheld from calibration, including trials from faster walking speeds. These models were used in all subsequent activation- and synergy-driven optimal control simulations.

Subject-Specific Muscle Synergy Control Model

The standard method of performing muscle synergy analysis uses only processed experimental EMG data. Following processing, EMG data for one complete walking cycle are organized into a matrix, and non-negative matrix factorization is applied via an optimization approach to decompose the high dimensional set of processed EMG signals into a lower dimensional set of time-varying signals (which we will call “synergy controls”) with associated sets of muscle weights (which we will call “synergy vectors”) (Lee and Seung, 1999; Tresch et al., 1999). The synergy vector weights specify how each synergy control contributes to each processed EMG signal. Each time a synergy analysis is performed, the number of synergies to fit must be specified a priori. Typically only three to six synergy controls are needed to reconstruct a much larger number of processed EMG signals with a high “variability account for” (VAF), typically above 90% (Ivanenko et al., 2005; Clark et al., 2010).

There are at least two drawbacks to the standard synergy analysis approach. First, the absolute amplitude of each processed EMG signal remains unknown. Though EMG signals are commonly normalized using data from a maximum voluntary contraction trial or the movement trial with maximum signal amplitude, maximal M-wave measurements reveal that these methods do not yield the true maximum EMG value (Clark et al., 2006; Simonsen et al., 2012; Racinais et al., 2013; Cronin et al., 2015). Since EMG amplitudes affect the results of a muscle synergy analysis, this issue makes it difficult to use experimentally derived synergy information for musculoskeletal modeling purposes. Second, since standard synergy analysis only uses experimental EMG data, it does not provide any information on how the inter-muscle coupling quantified by the synergy vectors affects an individual’s ability to perform specific movement tasks. Thus, a gap exists between subject-specific neural control information provided by standard muscle synergy analysis and the functional consequences of that information.

To bridge this gap while also addressing the EMG normalization issue, we perform muscle synergy analysis for our subject within the larger context of producing a dynamically consistent full-body walking motion using a subject-specific neuromusculoskeletal model (Sharif Razavian et al., 2015). With standard synergy analysis, the goal is to find a specified number of synergy controls and vectors that best match a larger set of normalized muscle EMG signals. With our approach, the goal is to find a specified number of synergy controls and vectors that best match experimental joint motions, ground reactions, lower-body inverse dynamic joint moments, and scaled EMG signals produced by our EMG-driven models. Thus, the breadth of data matched by our approach is much larger than that of a standard muscle synergy analysis. We describe our expanded approach as “dynamically consistent synergy analysis.” Use of EMG-driven models within our neuromusculoskeletal model provides a unique way to address the EMG normalization issue, while finding dynamically consistent full-body walking motions with our neuromusculoskeletal model bridges the gap between neural control information and its functional consequences.

To find synergy controls and vectors that could reproduce our subject’s walking data at his self-selected speed of 0.5 m/s, we followed a two-step process. First, we performed standard synergy analysis on the 35 muscle activations for each leg produced by the EMG-driven model calibration process for a representative walking trial. We performed this step on muscle activations rather than scaled EMG signals since we omitted activation dynamics from our final neuromusculoskeletal model to reduce model complexity. We incremented the number of synergies found by standard synergy analysis until the total VAF was greater than 95% and the VAF for each muscle was greater than 85%. We chose these high values since our goal was not simply reconstruction of activations but also reproduction of a dynamically consistent walking motion. Five synergies were required to achieve the target VAF values. Second, we performed a tracking optimization using direct collocation optimal control where all muscles in each leg were driven by five synergy controls. The synergy controls and associated synergy vectors were unknowns to be found by the optimization. The optimization tracked ground reactions, muscle activations, lower-body joint torques, and upper body joint motions while producing a dynamically consistent walking motion (see Optimal Control Walking Predictions below for further details).

Optimal Control Walking Predictions

We used the subject-specific neuromusculoskeletal model described above and direct collocation optimal control to predict the subject’s walking motion at 0.5 and 0.8 m/s given walking data from the most periodic trial collected at his self-selected speed of 0.5 m/s. A key advantage of collocation methods over shooting methods is that they use implicit rather than explicit simulation. Repeatedly during the non-linear programing (NLP) solution process, shooting methods perform explicit simulation to solve the system dynamics sequentially for one time frame at a time via numerical integration. This process is often unstable because either errors accumulate with each integration step or the system being simulated is inherently unstable, as with the human body during walking. By contrast, as part of the NLP solution process, collocation methods perform implicit simulation to solve the system dynamics for all time frames simultaneously with no notion of time stepping. Consequently, instabilities arising from accumulated integration errors or inherent system instabilities are eliminated, facilitating the use of gradient-based optimization for predicting motion. Another advantage of collocation is that feedback control, which is artificially introduced in explicit simulation to maintain system stability, is unnecessary since time stepping is not performed by implicit simulation. An overview of various numerical methods for solving optimal control problems, along with a brief discussion of the advantages and disadvantages of each, can be found in Limebeer and Rao (2015).

We investigated how well the predictions worked using three different control situations: joint torque controls (5 per leg – termed “torque-driven”), muscle activation controls (35 per leg – termed “activation-driven”), and muscle synergy controls (5 per leg – termed “synergy-driven”). To control the motion of each leg, the torque-driven problems used 5 joint torques rather than 35 muscles, the activation-driven problems used 35 muscles controlled by 35 independent muscle activations, and the synergy-driven problems used 35 muscles controlled by 5 independent muscle synergy controls that were linearly combined to produce 35 muscle activations. For all three control situations, hip internal–external rotation and toes flexion–extension in both legs along with the three pelvis rotations were found by tracking the corresponding joint angles from the 0.5 m/s periodic trial.

To generate our predictions, we performed a sequence of three categories of optimizations: (1) calibration optimizations, (2) tracking optimizations, and (3) prediction optimizations (see Table 1 for overview). This sequence was needed since large-scale direct collocation optimal control problems are often sensitive to the initial guess, making it helpful to increase the complexity of the problems being solved in a gradual and systematic fashion. Furthermore, for each category of optimization, initial problems were solved where the skeletal dynamic constraints (as quantified by pelvis residual loads) were not enforced, and these constraints were gradually tightened in subsequent problems until the desired tolerance was met. Below we describe concepts common to all three optimization categories, integration of OpenSim functionality into the optimal control framework, and details for the three optimization categories.

TABLE 1

Table 1. Sequence of direct collocation optimal control problems solved using GPOPS-II to predict patient-specific walking motions at 0.5, 0.8, and 1.1 m/s.

Common Concepts

The walking predictions reported in this study were generated using GPOPS-II, a direct collocation optimal control toolbox for MATLAB (Patterson and Rao, 2014). GPOPS-II solves for the state x(t), control u(t), and static parameters p that minimize the cost functional:

$J=\varphi \left(x(t_0),t_0,x(t_f),t_f,p\right)+{\displaystyle {\int }_{t_0}^{t_f}g(x(t),u(t),p)}dt$ (3)

subject to the constraints

$\dot{x}(t)=f(x(t),u(t),t,p)\text{,\hspace{0.17em}}\left(\text{dynamic\hspace{0.17em}constraints}\right),$ (4)

$c_{\min }\le f\left(x(t),u(t),t,p\right)\le c_{\max },\left(\text{algebraic\hspace{0.17em}constraints}\right),$ (5)

$b_{\min }\le b(x(t_0),t_0,x(t_f),t_f,p)\le b_{\max },\left(\text{boundary\hspace{0.17em}conditions}\right).$ (6)

The state and control are parameterized using variable-order Gaussian quadrature orthogonal collocation methods and formulated into a NLP problem.

Within the toolbox, two NLP solvers can be utilized: SNOPT (Gill et al., 2005) or IPOPT (Biegler and Zavala, 2009). SNOPT employs a quasi-Newton SQP active set method where the inverse of the Hessian is approximated using a Broyden–Fletcher–Goldfarb–Shanno (BFGS) update. At each step in the optimization, all linear constraints are satisfied and the active set is estimated. SNOPT achieves convergence via a merit function that is the objective function plus the sum of constraint infeasibilities. IPOPT employs an interior-point method where the goal is to satisfy the constraints via a barrier function. IPOPT does not distinguish between linear and non-linear constraints; all constraints are treated the same. IPOPT operates in two different modes. In first derivative mode, IPOPT employs a quasi-Newton method where the inverse of the Hessian is estimated using a BFGS update. In second derivative mode, IPOPT employs a full-Newton method where GPOPS-II provides a sparse representation of the lower triangular part of the Hessian of the Lagrangian of the NLP problem. Based on convergence and computation time considerations, we chose IPOPT in first derivative mode as the initial NLP solver for all optimal control problems. Though GPOPS-II contains an adaptive mesh refinement algorithm, we used a fixed mesh of 50 collocation points, divided into 10 intervals, over the entire gait cycle to reduce computation time (Ackermann and van den Bogert, 2010).

In this study, we developed a jerk-controlled inverse dynamic problem formulation using the 31 DOF skeletal dynamic equations generated by OpenSim (Delp et al., 2007). This approach was used since it exhibited improved convergence properties and solution smoothness over a directly controlled forward dynamic problem formulation (van den Bogert et al., 2011), but at the cost of added joint jerk controls. It also allowed us to violate the skeletal dynamics when starting from a poor initial guess, facilitating finding a feasible solution. For this problem formulation, we defined the control u(t) as the third time derivative of the generalized coordinates q(t) (i.e., joint jerk). The problem state consisted of the coordinate positions q(t), velocities v(t), and accelerations a(t), simplifying the dynamics of the optimal control problem to:

$\left[\begin{array}{c}\dot{q}\\ \dot{v}\\ \dot{a}\end{array}\right]=\left[\begin{array}{c}v\\ a\\ u\end{array}\right]$ (7)

To make the joint jerk controls unique, we added minimization of joint jerk to every optimal control problem. Jerk minimization terms in the cost function were scaled by $t_f^6$, where t_f is the specified final time, since jerk magnitude is proportional to $t_f^3$ based on analysis of analytic functions, and joint jerk squared is minimized in the cost function. Since joint jerk magnitudes change with final time, we did not want jerk minimization to affect the final times predicted by our optimizations.

For our inverse dynamic problem formulation, enforcement of the skeletal dynamics was achieved using algebraic path constraints. Two types of path constraints were used for this purpose. The first type ensured that residual forces and torques acting on the pelvis were eliminated. Each time an inverse dynamic analysis was performed with the OpenSim model, six pelvis residual loads R_pelvis(x(t)) were calculated. To enforce dynamic consistency, we added the algebraic path constraint

$R_{\text{min}}\le R_{\text{pelvis}}\le R_{\text{max}}$ (8)

to constrain the six pelvis residual loads to be within a specified tolerance. The second type of constraint enforced consistency between net muscle moments calculated from additional controls and lower extremity joint moments calculated from inverse dynamics. The additional controls were 5 joint torques per leg for torque-driven problems, 35 activations per leg for activation-driven problems, 5 synergies per leg for synergy-driven problems. For all three control situations, we added algebraic path constraints to ensure that the additional controls balanced five lower extremity inverse dynamic joint moments in each leg (i.e., hip flexion–extension, hip adduction–abduction, knee flexion–extension, ankle flexion–extension, and ankle inversion–eversion). Convergence tolerances were set to 1 N and 0.1 Nm for residual forces and moments, respectively. When this tolerance was met, we considered the resulting motion to be dynamically consistent, requiring negligible fictitious external loads to balance the dynamic equations. All optimal control solutions presented in this study utilized this inverse dynamic problem formulation. However, the cost function, controls, and constraints varied depending on the category of optimization problem being solved.

OpenSim Integration

Development of dynamically consistent full-body walking predictions required integrating OpenSim functionality into the MATLAB environment in a computationally efficient manner (Figure 1). To perform the integration, we took advantage of two aspects of our optimal control solution process. First, use of joint jerk controls made complete OpenSim model state and state derivative information available for any time point being analyzed, making it possible to evaluate the OpenSim skeletal dynamic equations in an inverse sense and avoiding the need to calculate OpenSim model state derivatives for forward solutions. Second, use of a direct collocation solution approach allowed each time point to be analyzed independently from all others, making it possible to employ parallel processing methods on a time-point-by-time-point basis.

FIGURE 1

Figure 1. Flowchart demonstrating the interaction of MATLAB and OpenSim functions for synergy-driven prediction optimizations performed within the GPOPS-II MATLAB environment. For each iteration of the NLP solver, GPOPS-II provides a new guess for the model state and control, which are used by a series of MATLAB and OpenSim functions to calculate path constraints, terminal constraints (not shown – used to enforce near-periodicity in motion and ground reactions), and the cost function. The cost function requires knowledge of state and control quantities that are being tracked, which are taken from the results of the corresponding tracking optimization. Path constraints are used to satisfy skeletal dynamics and to make lower-body joint moments calculated by an OpenSim inverse dynamic analysis match corresponding joint moments calculated by a MATLAB EMG-driven model. Synergy-driven prediction optimizations use five synergy controls per leg as controls. Activation-driven prediction optimizations use 35 muscle activations per leg as controls and do not require the MATLAB muscle synergy model. Torque-driven prediction optimizations use five joint torques per leg as controls and do not require the MATLAB muscle synergy model or EMG-driven model. Calibration optimizations include additional path constraints to match additional experimental quantities (e.g., ground reaction forces and moments).

To take advantage of these two aspects, we created MATLAB MEX functions in C++ that parallelized two OpenSim tasks: (1) calculation of positions and velocities of all contact elements on the bottom of each foot, and (2) calculation of inverse dynamic joint torques. We parallelized OpenSim point kinematic calculations using OpenMP and OpenSim inverse dynamic calculations using MPI. For parallelized point kinematic calculations, inputs were the OpenSim model state; while for parallelized inverse dynamic calculations, inputs were the OpenSim model state and calculated external ground reactions.

Using these parallelized functions, we employed an efficient series of calculation steps each time GPOPS-II needed to evaluate the OpenSim skeletal dynamic equations across all time points. First, we called our parallel C++ point kinematics function to calculate the position and velocity of every contact element in the model. Second, we used this point kinematic information to calculate the contact force acting on each element and the net contact force and torque to be applied to the rear foot and toes segments of each foot. Third, we evaluated our surrogate musculoskeletal geometry using the current joint positions and velocities to find muscle-tendon lengths, velocities, and moment arms for all muscles in the model. Fourth, we used our Hill-type muscle model with the current muscle activations and muscle-tendon lengths and velocities to calculate all tendon forces. Fifth, we calculated muscle contributions to all joint moments by multiplying tendon forces by their associated moment arms. Sixth, we called our parallel C++ inverse dynamics function to apply the calculated ground reactions to the rear foot and toes segments of both feet as OpenSim external forces and calculate inverse dynamic joint torques. Seventh, we calculated joint torque errors for the five muscle-controlled lower-body joints in each leg by subtracting the joint torques calculated in OpenSim using inverse dynamics from the joint torques calculated in MATLAB using GPOPS-II additional controls and moment arms where needed. All of these calculations were performed in MATLAB apart from the two parallel C++ functions that called OpenSim functionality.

Calibration Optimizations

We performed two calibration optimizations to determine parameter values in the foot–ground contact models and to provide a starting point for subsequent tracking optimizations (see Table 1, top section). For both optimizations, the cost function tracked experimental marker, ground reaction, and inverse dynamic joint torque data from the periodic trial while also minimizing joint jerk. The constraints enforcing dynamic consistency and joint jerk were the only controls. IPOPT was used in first derivative mode to generate an initial guess, and then SNOPT was started at that initial guess to refine the solution. For the first optimization, static parameters were included for the six coefficients that defined the parabolic distribution of spring stiffness across the bottom of each foot and for the common spring resting length for each foot. For the second optimization, all static parameters were removed from the problem, the parameter values were fixed to those found by the first optimization, and the optimization repeated to verify that identified foot–ground contact model parameter values could closely reproduce the experimentally measured ground reactions for both feet. Root-mean-square (RMS) errors in ground reactions produced by the second optimization were within 10 N for forces and 5 Nm for moments calculated about the midfoot marker projected onto the floor.

Tracking Optimizations

Starting from the results of the second calibration optimization, we performed three tracking optimizations for the 0.5 m/s walking speed to provide a starting point for subsequent prediction optimizations (see Table 1, middle section). Each tracking optimization utilized a different method to control 5 lower-body DOFs (hip flexion–extension, hip adduction–abduction, knee flexion–extension, ankle flexion–extension, and ankle inversion–eversion) in each leg. The first tracking optimization was torque-driven (5 independent joint torque controls per leg), the second was activation-driven (35 independent muscle activation controls per leg), and the third was synergy-driven (5 independent synergy controls per leg used to construct the 35 muscle activations per leg). For all three optimizations, the cost function tracked lower-body joint torques (apart from hip internal-external rotation and toes flexion–extension) and upper body, hip internal-external rotation, and toes joint angles from the second calibration optimization while minimizing joint jerk controls. We did not track joint angle and joint torque curves for the same joint to avoid having related terms for the same joint in the cost function. We selected joint angle tracking over joint torque tracking for the toes and upper body joints since we found that large changes in toes or arm motion required only small changes in the corresponding joint torques. We also tracked pelvis angles so that the model would maintain the proper orientation in the laboratory. Path constraints satisfied skeletal dynamics and bounded joint angle errors relative to the second calibration optimization and ground reaction errors relative to experimental data. For the torque-driven optimization, the controls were joint jerk and 5 joint torques for each leg, and algebraic path constraints were used to match 5 inverse dynamic joint torques per leg with the 5 joint torque controls per leg. The activation-driven optimization replaced the 5 torque controls with 35 lower extremity muscle activations per leg, added tracking of 35 activation controls per leg from the EMG-driven models to the cost function, and used algebraic path constraints to match 5 inverse dynamic joint torques per leg with activation controls per leg. The synergy-driven optimization was identical except that the 35 activation controls were replaced with 5 synergy controls per leg, inverse dynamic joint torque matching used synergy-constructed muscle activations, and static parameters were added to allow identification of the corresponding 5 sets of synergy vector weights. In addition, a constraint was added to force each synergy vector to have unit magnitude, making the synergy solutions unique. IPOPT in first derivative mode was used to solve all tracking optimization problems, with gradients calculated using central differencing.

Prediction Optimizations

Starting from the results of the tracking optimizations, we performed three prediction optimizations for the 0.5 m/s walking speed as a “sanity check,” three for the 0.8 m/s walking speed as a predictive evaluation, and one using only synergy controls for a 1.1 m/s walking speed as a challenge to our simulation framework. For the 0.5 m/s and 0.8 m/s predictions, the first tracking optimization was torque-driven, the second was activation-driven, and the third was synergy-driven. For the 1.1 m/s prediction, only a synergy-driven model was used to make predictions. The goal of these optimizations was to see whether each type of control could predict not only a realistic walking motion with realistic ground reactions but also the correct period for one gait cycle, which decreases with increasing walking speed. Since the 0.5 m/s speed was the same as the tracking optimization speed, matching the experimental period closely provides a “sanity check” on the solution. Due to problems we encountered with solving these optimal control problems using free final time, we performed each prediction optimization using five values of final time: the average experimental value rounded to the nearest tenth of a second plus or minus 0.1 and 0.2 s. We then fit a parabola to the final cost function values plotted as a function of final time and took the minimum value as the predicted final time. All results reported were taken from the simulations whose final times were closest to the parabola’s minimum value.

For each control situation, the optimal control problem formulation for both speeds was similar to the corresponding tracking optimization except that several cost function and constraint terms were removed and no static parameters were utilized, allowing new walking motions and ground reactions to be predicted. For all three optimizations, the cost function tracked upper body, hip internal–external rotation, and toes flexion–extension joint angles from the corresponding tracking optimization solution while minimizing joint jerk controls, and the path constraints satisfied the skeletal dynamics. The torque-driven problem added tracking of 5 lower extremity joint torques per leg found by the corresponding tracking optimization to the cost function. The activation-driven problem added 35 lower extremity muscle activations per leg to the controls, tracking of 35 activations per leg found by the corresponding tracking optimization to the cost function, and matching of 5 inverse dynamic joint torques per leg with corresponding joint torques produced by muscle activations to the algebraic path constraints. The synergy-driven problem added 5 synergy controls per leg to the controls, tracking of 5 synergy controls per leg found by the corresponding tracking optimization to the cost function, and matching of 5 inverse dynamic joint torques per leg with corresponding joint torques produced by the synergy-constructed muscle activations to the algebraic path constraints. The synergy-driven problem used the synergy vectors found by the corresponding tracking optimization. The implicit assumption in this prediction approach is that when the subject walks at any speed, he will choose controls (joint torques, muscle activations, or synergy controls) that are “close” to those he uses at his self-selected speed of 0.5 m/s. IPOPT in first derivative mode was again used to solve all prediction optimization problems.

Results

Tracking optimizations using all three types of controls closely reproduced the subject’s experimentally measured walking motion at 0.5 m/s. Simulated lower and upper body joint angles were within experimental ranges from multiple walking cycles (Figure 2 – first column), as were simulated lower-body joint torques (Figure 3 – first column). However, joint torques from the torque-driven model were generally less smooth than those from the activation-driven and synergy-driven models, especially for the two hips. Simulated ground reaction forces were also within experimental ranges (Figure 4 – first column). For the activation-driven and synergy-driven models, simulated activations were within the experimental ranges determined by the EMG-driven models (Figures 5 and 6), with only small changes in activations needed to produce dynamically consistent walking motions.

FIGURE 2

Figure 2. Sagittal plane lower extremity joint angles from inverse kinematics (gray lines for individual trials), torque-driven optimizations (red lines), activation-driven optimizations (green lines), and synergy-driven optimizations (blue lines). First column: results from tracking optimizations at 0.5 m/s. Second column: results from prediction optimizations at 0.5 m/s. Third column: results from prediction optimizations at 0.8 m/s. Top three rows: hip, knee, and ankle angles from left non-paretic leg. Bottom three rows: hip, knee, and ankle angles from right paretic leg.

FIGURE 3

Figure 3. Sagittal plane lower extremity joint torques from inverse dynamics (gray lines for individual trials), torque-driven optimizations (red lines), activation-driven optimizations (green lines), and synergy-driven optimizations (blue lines). First column: results from tracking optimizations at 0.5 m/s. Second column: results from prediction optimizations at 0.5 m/s. Third column: results from prediction optimizations at 0.8 m/s. Top three rows: hip, knee, and ankle torques from left non-paretic leg. Bottom three rows: hip, knee, and ankle torques from right paretic leg.

FIGURE 4

Figure 4. Ground reaction forces from the experimental data (gray lines for individual trials), torque-driven optimizations (red lines), activation-driven optimizations (green lines), and synergy-driven optimizations (blue lines). First column: results from tracking optimizations at 0.5 m/s. Second column: results from prediction optimizations at 0.5 m/s. Third column: results from prediction optimizations at 0.8 m/s. Top three rows: anterior (Fx), superior (Fy), and lateral (Fz) forces from left non-paretic leg. Bottom three rows: anterior (Fx), superior (Fy), and lateral (Fz) forces from right paretic leg.

FIGURE 5

Figure 5. Left leg muscle activations from the EMG-driven model (gray lines for individual trials), activation-driven tracking optimizations (green lines), and synergy-driven tracking optimizations (blue lines) at 0.5 m/s. Muscle names are indicated in the figure.

FIGURE 6

Figure 6. Right leg muscle activations from the EMG-driven model (gray lines for individual trials), activation-driven tracking optimizations (green lines), and synergy-driven tracking optimizations (blue lines) at 0.5 m/s. Muscle names are indicated in the figure.

Prediction optimizations using all three types of controls also closely reproduced the subject’s experimentally measured walking motion, gait period, and stride length at 0.5 m/s (Table 2 – top rows). Simulated lower and upper body joint angles (Figure 2 – second column), lower-body joint torques (Figure 3 – second column), and ground reaction forces (Figure 4 – second column) were extremely similar to those produced by the corresponding tracking optimizations. The most noticeable minor differences were for the torque-driven model, where the predicted hip torques were again less smooth than for the other two control types, as were the predicted medial–lateral ground reaction forces. For the activation-driven and synergy-driven models, simulated activations were also very similar to those found by the corresponding tracking optimizations (Figures 7 and 8). All three models predicted a realistic gait period and stride length (Figure 11, first column; Table 2, top rows).

TABLE 2

Table 2. Comparison of gait period and stride length between the median experimental walking motions (±experimental ranges) and corresponding predicted walking motions for 0.5 and 0.8 m/s gait speeds.

FIGURE 7

Figure 7. Left leg muscle activations from the EMG-driven model (gray lines for individual trials), activation-driven prediction optimizations (green lines), and synergy-driven prediction optimizations (blue lines) at 0.5 m/s. Muscle names are indicated in the figure.

FIGURE 8

Figure 8. Right leg muscle activations from the EMG-driven model (gray lines for individual trials), activation-driven prediction optimizations (green lines), and synergy-driven prediction optimizations (blue lines) at 0.5 m/s. Muscle names are indicated in the figure.

Only prediction optimizations that used activation controls or synergy controls closely reproduced the subject’s experimentally measured walking motion, gait period, and stride length at 0.8 m/s (Table 2 – bottom rows). For these two control types, simulated lower joint angles were within or just outside experimental ranges from multiple walking cycles, while for torque controls, every simulated lower-body joint angle went beyond the experimental ranges (Figure 2 – third column), and the optimization predicted hiking of the paretic hip to compensate for reduced knee flexion on that side, which is a biologically plausible strategy. Interestingly, for all three control types, simulated upper body joint angles were extremely different from those observed experimentally (see Figure 12 for a visual comparison between experimentally measured and computationally predicted full-body motions). Similar to the joint angle results, simulated lower-body joint torques (Figure 3 – third column) and ground reaction forces (Figure 4 – third column) produced by activation and synergy controls generally remained within or at the edge of experimental ranges, while those produced by torque controls generally went beyond them, especially for the ground reaction forces. All three control types predicted an increase in anterior–posterior ground reaction force that was comparable to the increased range measured experimentally to achieve an increase in walking speed. For the activation-driven and synergy-driven models, simulated activations remained within the experimental ranges determined by the EMG-driven models (Figures 9 and 10). The synergy-driven model predicted the most realistic gait period and stride length, which were 1.16 s compared to an experimental median of 1.14 s for gait period and 0.93 m compared to an experimental median of 0.91 m for stride length (Figure 11, second column; Table 2, bottom rows). The gait period and stride length predictions for the activation-driven model were right at the boundary of the measured experimental ranges, while the predictions for the torque-driven model were outside those ranges.

FIGURE 9

Figure 9. Left leg muscle activations from the EMG-driven model (gray lines for individual trials), activation-driven prediction optimizations (green lines), and synergy-driven prediction optimizations (blue lines) at 0.8 m/s. Muscle names are indicated in the figure.

FIGURE 10

Figure 10. Right leg muscle activations from the EMG-driven model (gray lines for individual trials), activation-driven prediction optimizations (green lines), and synergy-driven prediction optimizations (blue lines) at 0.8 m/s. Muscle names are indicated in the figure.

FIGURE 11

Figure 11. Normalized cost function values as a function of specified gait period (i.e., final time) as produced by torque-driven (top row – red circles), activation-driven (middle row – green circles), and synergy-driven (bottom row – blue circles) prediction optimizations performed at 0.5 m/s (first column) and 0.8 m/s (second column) walking speeds. The thick black line in each plot indicates a quadratic fit to the cost function values. The thick gray bars indicate the range of experimentally measured gait periods, and the black hash marks indicate the median gait period. For the activation-driven prediction at 0.8 m/s, the normalized cost function value for a gait period of 0.9 s (not shown) was much higher than expected, likely due to entrapment in a local minimum, and was omitted from the quadratic fit.

FIGURE 12

Figure 12. Animation strips showing walking motions at 0.8 m/s obtained from the most periodic experimentally measured walking trial (first row), the torque-driven prediction optimization (second row), the activation-driven prediction optimization (third row), and the synergy-driven prediction optimization (fourth row). Each column represents a different point in the gait cycle starting at 0% on the far left and going to 100% on the far right.

The final prediction optimization used synergy controls to simulate walking at 1.1 m/s, a condition for which no experimental data are available, and predicted a gait period of 1.12 s. Assuming a linear relationship between walking speed and stride length (Jordan et al., 2007), as defined by the subject’s median speed and stride length for 0.5 and 0.8 m/s, the estimated experimental gait period for this faster walking speed would be 1.06 s. Thus, the synergy-driven model predicted not only the best gait period for the 0.8 m/s walking speed but also a physically realistic gait period with physically realistic walking motion (not shown) for an even faster walking speed without available experimental data.

Discussion

This study developed and evaluated a subject-specific synergy-controlled neuromusculoskeletal simulation framework that predicted three-dimensional walking motions for an individual post-stroke. We investigated whether actuating a neuromusculoskeletal model of the subject with muscle synergy controls (5 per leg) facilitated generation of accurate walking predictions compared to actuating the model with muscle activation controls (35 per leg) or joint torque controls (5 per leg). We found that walking predictions generated for both 0.5 and 0.8 m/s were more accurate (in terms of joint motions, joint torques, ground reactions, and final time) and converged more easily for the activation-driven and synergy-driven models than for the torque-driven model. Furthermore, the accuracy of the walking predictions at both speeds was comparable for the activation-driven and synergy-driven models, even though the synergy-driven model used 30 fewer controls per leg. When the synergy-driven model was used to predict the subject’s walking motion at 1.1 m/s, a condition for which no experimental data were available, the predicted stride length was close to that calculated by the linear speed–stride length relationship fitted to experimental data from 0.5 and 0.8 m/s. Overall, these findings are encouraging and suggest that the current simulation framework could provide a useful foundation for predicting how a patient will interact with different neurorehabilitation approaches (e.g., FES, AFO, exoskeleton, robotic gait trainer, strength training) so that an optimal neurorehabilitation prescription can be identified.

We had to overcome a number of practical challenges to generate the walking predictions presented in this study. One significant challenge was obtaining convergence of our optimal control problems, which we addressed using several strategies. First, we used sequences of optimizations that gradually increased the complexity of the problem being solved and provided a good initial guess for the next level problem. Second, we compared IPOPT in first derivative mode, IPOPT in second derivative mode, and SNOPT to determine which option worked best for our problems. Third, we tracked joint angles rather than joint torques for joints not controlled by muscles, since tracking joint torques permitted large changes in toes and arm motions using only small changes in joint torques. Fourth, we replaced free final time problems with a sequence of five fixed final time problems and then fitted a parabola to the cost function values to determine the final time of the corresponding free final time problem. Fifth, we replaced an explicit dynamics formulation requiring no additional controls with an implicit dynamics formulation requiring additional jerk controls to obtain smooth motion and control predictions. The other significant challenge was obtaining reasonable computation times, which we resolved by parallelizing the computationally costly skeletal dynamics and point kinematics calculations performed by OpenSim. Once a reasonable initial guess was available, torque-driven and synergy-driven problems typically required about 30 min of CPU time to converge, though the variation in convergence time was wide, while activation-driven typically required about an hour of CPU time. On a practical basis, however, we normally perturbed the most recent solution and re-ran each optimization to help avoid entrapment in a local minimum, making estimation of total CPU time difficult. More research is needed to determine how to improve the convergence properties of these problems.

Given that numerous previous studies have generated muscle-actuated full-body forward dynamic simulations of walking (Gerritsen et al., 1998; Anderson and Pandy, 2001; Thelen et al., 2003; Ackermann and van den Bogert, 2010, 2012; Geyer and Herr, 2010; McGowan et al., 2010; Allen et al., 2013; Knarr et al., 2013, 2014; Kia et al., 2014; Dorn et al., 2015), it is worth considering the unique aspects of our approach. Only a small number have predicted new walking motions for which experimental data are not available (Anderson and Pandy, 2001; Ackermann and van den Bogert, 2010, 2012; Dorn et al., 2015). Few studies have used direct collocation optimal control methods to solve for all time points simultaneously (Ackermann and van den Bogert, 2010, 2012), few have used subject-specific (Knarr et al., 2014) rather than scaled generic or simplified musculoskeletal models, and few have modeled individuals with neurological impairment (Allen et al., 2013; Knarr et al., 2013). No previous muscle-actuated full-body walking simulation studies have combined subject-specific EMG-driven modeling with subject-specific synergy controls to define the neural control structure of the model. Furthermore, no previous studies have calibrated joint and ground reaction force parameters in the model to match walking data collected from a specific subject. Thus, the most unique feature of our study was combining all of the various modeling and optimization elements listed above into a single comprehensive simulation framework capable of predicting new walking motions.

Only the activation-driven and synergy-driven models were able to predict the correct walking motion and gait period closely at both 0.5 and 0.8 m/s. Since the only difference between the torque-driven model and the other two models was the control method, it is likely that numerical issues related to the use of pure torque controls were the source of the problem. This hypothesis is supported by results from two previous simulation studies. Risher et al. (1997) showed that even small inconsistencies in inverse dynamic solutions, such as those introduced by spline fitting, can produce large motion errors when the calculated joint torques are used to control a forward dynamic simulation intended to reproduce the original motion. Our approach requires spline fitting of joint torque and joint angle data so that the optimal control solver can obtain values at any desired collocation point. Gerritsen et al. (1998) showed that muscle force–length–velocity properties provide proportional-derivative-like feedback control properties that help stabilize forward dynamic simulations of walking and prevent drift away from a desired motion. Based on the finding of these two studies, it is less surprising that the activation-driven and synergy-driven models performed better than the torque-driven model.

Although both our synergy-driven and activation-driven models generated accurate walking predictions, our synergy-driven model still possesses several distinct advantages. The primary advantage is the significantly reduced number of controls compared to the activation-driven model. By using muscle synergy controls, we were able to predict highly realistic subject-specific walking motions using the same number of controls as in the torque-driven model but without any of the problems encountered by that model. Significantly fewer controls (30 fewer per leg) reduce computational cost and complexity and make the optimal control solution process less sensitive to poor initial guesses. Though not demonstrated in our study, we believe that the synergy-driven model also has the best potential for simulating individuals with neurological impairment, since subject-specific synergy information limits how a subject can coordinate his muscles. For example, it would be interesting to eliminate one synergy at a time from the paretic leg in our model and predict the functional impact on our subject’s walking pattern. Would the model’s ability to reproduce the subject’s walking motion break down if a lower number of synergies were used? Synergy controls within a neuromusculoskeletal model provide an excellent avenue for simulating the functional consequences of reduced neural complexity (Allen et al., 2013).

To use our neuromusculoskeletal simulation framework for actual clinical treatment design, we would need to determine how to incorporate different treatment approaches into the framework. Simulating the effects of strength training could be achieved by increasing the peak isometric strength of individual muscles or groups of muscles. Simulating the effects of FES could be achieved by adding controls to the optimal control problem that augment the activation of one muscle or several muscles. Simulating the effects of an AFO, exoskeleton, robotic gait trainer, or exercise device could be achieved by adding a model of the device to the patient’s OpenSim model [e.g., Fregly et al. (2015)] and adding static parameters and controls to the optimal control problem that account for modifiable design features of the device. The biggest challenge with simulating any of these treatment approaches is defining an appropriate optimization cost function. When predicting how the patient’s neural control system will respond to a treatment, should the cost function minimize an absolute quantity, such as metabolic cost, should it minimize neural control changes away from some experimentally measured baseline situation (as we have done here), or should it hold the patient’s neural control strategy constant and change only static parameters and controls related to the treatment? The issue of predicting how a patient will interact with a treatment is a critical one for researchers to explore in the future.

Despite the high level of subject specificity in our neuromusculoskeletal model, our study possesses a number of limitations that can help inform future research efforts. First, we modeled only a single hemiparetic subject. Our goal for the present study was to develop and evaluate the initial implementation of our neuromusculoskeletal simulation framework, which required use of only a single subject. In the future, we plan to test the framework further using walking data collected from additional hemiparetic subjects. Second, we performed all model calibration steps (lower-body kinematic model, foot–ground contact model, EMG-driven model, muscle synergy model) using only a static trial and walking data, primarily from the subject’s self-selected speed of 0.5 m/s. Though this limitation was planned to simplify the model calibration process, use of a wider variety of calibration movements could improve the predictive capabilities of the model. Third, we did not model any neural feedback mechanisms (e.g., from muscle spindles and Golgi tendon organs). Though the extent to which feedback mechanisms contribute to the control of walking remains controversial, it is possible that inclusion of neural feedback models could have a significant impact on our predicted walking motions (Geyer and Herr, 2010). Fourth, we evaluated our predicted walking motions using only the calibration speed and a single faster non-calibration speed, where the faster speed was only 0.3 m/s faster. A more thorough evaluation would involve a wider range of speeds (our subject was unable to walk comfortably for an extended period of time above 0.8 m/s) and movement tasks. Fifth, while we controlled most lower-body joints with muscles, we tracked experimentally measured joint motions to determine the motion of the toes, hip internal–external rotation, and all upper body joints. Better prediction of toes and upper body motion in particular would likely be achieved if these joints could be controlled by muscles as well. Finally, even though our subject was hemiparetic, we assumed that most model parameter values were the same on both sides of the body, the exceptions being electromechanical delays, EMG scale factors, and muscle synergy control structure. When we allowed model parameter values to be different on the two sides, the accuracy with which our EMG-driven modeling process could predict joint moments improved little, suggesting that our current assumptions about bilateral symmetry were reasonable.

In conclusion, we have presented a novel combination of neuromusculoskeletal and computational modeling approaches that permits prediction of subject-specific walking motions for individuals with neurological impairment. The primary strengths of the approach are the high level of subject-specificity in the neuromusculoskeletal model calibration process, including a subject-specific representation of neural control limitations and capabilities via muscle synergy controls, the flexibility of the framework for changing the optimal control problem formulation and the characteristics of the OpenSim model, the ability to tie muscle synergies to their functional consequences, and the ability to predict new walking motions. With future developments, it may be possible to use our framework to simulate different neurorehabilitation interventions and ultimately to optimize treatment prescription so as to maximize recovery of walking function on an individual patient basis.

Author Contributions

AM assisted with experimental data collection, processed all experimental data, performed all neuromusculoskeletal modeling work, formulated and ran all optimal control problems, and participated in writing the manuscript draft. IE integrated OpenSim dynamics functionality into the GPOPS-II optimal control environment and helped with revising the manuscript draft. JJ assisted with development and calibration of the foot–ground contact models and helped with revising the manuscript draft. AR assisted with development of optimal control problem formulations, helped troubleshoot GPOPS-II optimal control issues, and helped with revising the manuscript draft. CP recruited the experimental subject, ran the experimental data collection session, provided clinical assessment of predicted walking motions, and helped with revising the manuscript draft. BF planned and supervised the entire project, organized the experimental data collection session, assisted with experimental data collection, directed collaborative efforts between AM, IE, and JJ, assisted with formulation of optimal control problems, evaluated optimal control results, and participated in writing the manuscript draft.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer JS declared a shared affiliation, though no other collaboration, with the authors AM, IE, JJ, AR, CP, and BF to the handling Editor, who ensured that the process nevertheless met the standards of a fair and objective review.

Funding

This study was funded by the National Science Foundation under grants CBET 1404767, CBET 1159735, and CBET 1052754 and by the National Institutes of Health under grants R01 EB009351 and U54 GM072970.

Supplementary Material

MATLAB code and sample experimental inputs for performing direct collocation optimal control predictions of walking will be made available at https://simtk.org/projects/synwalkingpred.

References

Ackermann, M., and van den Bogert, A. J. (2010). Optimality principles for model-based prediction of human gait. J. Biomech. 43, 1055–1060. doi: 10.1016/j.jbiomech.2009.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Ackermann, M., and van den Bogert, A. J. (2012). Predictive simulation of gait at low gravity reveals skipping as the preferred locomotion strategy. J. Biomech. 45, 1293–1298. doi:10.1016/j.jbiomech.2012.01.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Afschrift, M., De Groote, F., De Schutter, J., and Jonkers, I. (2014). The effect of muscle weakness on the capability gap during gross motor function: a simulation study supporting design criteria for exoskeletons of the lower limb. Biomed. Eng. Online 13, 111. doi:10.1186/1475-925X-13-111

PubMed Abstract | CrossRef Full Text | Google Scholar

Alibeji, N. A., Kirsch, N. A., and Sharma, N. (2015). A muscle synergy-inspired adaptive control scheme for a hybrid walking neuroprosthesis. Front. Bioeng. Biotechnol. 3:203. doi:10.3389/fbioe.2015.00203

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, J. L., Kautz, S. A., and Neptune, R. R. (2013). The influence of merged muscle excitation modules on post-stroke hemiparetic walking performance. Clin. Biomech. (Bristol, Avon) 28, 697–704. doi:10.1016/j.clinbiomech.2013.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

An, K. N., Takahashi, K., Harrigan, T. P., and Chao, E. Y. (1984). Determination of muscle orientations and moment arms. J. Biomech. Eng. 106, 280–282. doi:10.1115/1.3138494

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, F. C., and Pandy, M. G. (2001). Dynamic optimization of human walking. J. Biomech. Eng. 123, 381–390. doi:10.1115/1.1392310

PubMed Abstract | CrossRef Full Text | Google Scholar

Arnold, E. M., Ward, S. R., Lieber, R. L., and Delp, S. L. (2010). A model of the lower limb for analysis of human movement. Ann. Biomed. Eng. 38, 269–279. doi:10.1007/s10439-009-9852-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Auchstaetter, N., Luc, J., Lukye, S., Lynd, K., Schemenauer, S., Whittaker, M., et al. (2015). Physical therapists’ use of functional electrical stimulation for clients with stroke: frequency, barriers, and facilitators. Phys. Ther. 96, 995–1005. doi:10.2522/ptj.20150464

CrossRef Full Text | Google Scholar

Bae, Y.-H., Ko, Y. J., Chang, W. H., Lee, J. H., Lee, K. B., Park, Y. J., et al. (2014). Effects of robot-assisted gait training combined with functional electrical stimulation on recovery of locomotor mobility in chronic stroke patients: a randomized controlled trial. J. Phys. Ther. Sci. 26, 1949–1953. doi:10.1589/jpts.26.1949

PubMed Abstract | CrossRef Full Text | Google Scholar

Belda-Lois, J.-M., Mena-del Horno, S., Bermejo-Bosch, I., Moreno, J. C., Pons, J. L., Farina, D., et al. (2011). Rehabilitation of gait after stroke: a review towards a top-down approach. J. Neuroeng. Rehabil. 8, 66. doi:10.1186/1743-0003-8-66

PubMed Abstract | CrossRef Full Text | Google Scholar

Biegler, L. T., and Zavala, V. M. (2009). Large-scale nonlinear programming using IPOPT: an integrating framework for enterprise-wide dynamic optimization. Comput. Chem. Eng. 33, 575–582. doi:10.1016/j.compchemeng.2008.08.006

CrossRef Full Text | Google Scholar

Blair, S. N., Kohl, H. W., Paffenbarger, R. S., Clark, D. G., Cooper, K. H., and Gibbons, L. W. (1989). Physical fitness and all-cause mortality. A prospective study of healthy men and women. J. Am. Med. Assoc. 262, 2395–2401. doi:10.1001/jama.1989.03430170057028

PubMed Abstract | CrossRef Full Text | Google Scholar

Bogey, R., and Hornby, G. T. (2014). Gait training strategies utilized in poststroke rehabilitation: are we really making a difference? Top. Stroke Rehabil. 14, 1–8. doi:10.1310/tsr1406-1

CrossRef Full Text | Google Scholar

Bortole, M., Venkatakrishnan, A., Zhu, F., Moreno, J. C., Francisco, G. E., Pons, J. L., et al. (2015). The H2 robotic exoskeleton for gait rehabilitation after stroke: early findings from a clinical study. J. Neuroeng. Rehabil. 12, 54. doi:10.1186/s12984-015-0048-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Buchanan, T. S., Lloyd, D. G., Manal, K., and Besier, T. F. (2004). Neuromusculoskeletal modeling: estimation of muscle forces and joint moments and movements from measurements of neural command. J. Appl. Biomech. 20, 367–395. doi:10.1123/jab.20.4.367

PubMed Abstract | CrossRef Full Text | Google Scholar

Buesing, C., Fisch, G., O’Donnell, M., Shahidi, I., Thomas, L., Mummidisetty, C. K., et al. (2015). Effects of a wearable exoskeleton stride management assist system (SMA^®) on spatiotemporal gait characteristics in individuals after stroke: a randomized controlled trial. J. Neuroeng. Rehabil. 12, 69. doi:10.1186/s12984-015-0062-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Chantraine, F., Schreiber, C., Kolanowski, E., and Moissenet, F. (2016). Control of stroke-related genu recurvatum with prolonged timing of dorsiflexor functional electrical stimulation: a case study. J. Neurol. Phys. Ther. 40, 209–215. doi:10.1097/NPT.0000000000000137

PubMed Abstract | CrossRef Full Text | Google Scholar

Charlton, I. W., Tate, P., Smyth, P., and Roren, L. (2004). Repeatability of an optimised lower body model. Gait Posture 20, 213–221. doi:10.1016/j.gaitpost.2003.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Chung, Y., Kim, J.-H., Cha, Y., and Hwang, S. (2014). Therapeutic effect of functional electrical stimulation-triggered gait training corresponding gait cycle for stroke. Gait Posture 40, 471–475. doi:10.1016/j.gaitpost.2014.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, D. J., Condliffe, E. G., and Patten, C. (2006). Activation impairment alters muscle torque-velocity in the knee extensors of persons with post-stroke hemiparesis. Clin. Neurophysiol. 117, 2328–2337. doi:10.1016/j.clinph.2006.07.131

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, D. J., Ting, L. H., Zajac, F. E., Neptune, R. R., and Kautz, S. A. (2010). Merging of healthy motor modules predicts reduced locomotor performance and muscle coordination complexity post-stroke. J. Neurophysiol. 103, 844–857. doi:10.1152/jn.00825.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Crabtree, C. A., and Higginson, J. S. (2009). Modeling neuromuscular effects of ankle foot orthoses (AFOs) in computer simulations of gait. Gait Posture 29, 65–70. doi:10.1016/j.gaitpost.2008.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Cronin, N. J., Kumpulainen, S., Joutjärvi, T., Finni, T., and Piitulainen, H. (2015). Spatial variability of muscle activity during human walking: the effects of different EMG normalization approaches. Neuroscience 300, 19–28. doi:10.1016/j.neuroscience.2015.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

De Groote, F., Demeulenaere, B., Swevers, J., De Schutter, J., and Jonkers, I. (2012). A physiology-based inverse dynamic analysis of human gait using sequential convex programming: a comparative study. Comput. Methods Biomech. Biomed. Engin. 15, 1093–1102. doi:10.1080/10255842.2011.571679

PubMed Abstract | CrossRef Full Text | Google Scholar

De Groote, F., Kinney, A. L., Rao, A. V., and Fregly, B. J. (2016). Evaluation of direct collocation optimal control problem formulations for solving the muscle redundancy problem. Ann. Biomed. Eng. doi:10.1007/s10439-016-1591-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Delp, S. L., Anderson, F. C., Arnold, A. S., Loan, P., Habib, A., John, C. T., et al. (2007). OpenSim: open-source software to create and analyze dynamic simulations of movement. IEEE Trans. Biomed. Eng. 54, 1940–1950. doi:10.1109/TBME.2007.901024

PubMed Abstract | CrossRef Full Text | Google Scholar

Dorn, T. W., Wang, J. M., Hicks, J. L., and Delp, S. L. (2015). Predictive simulation generates human adaptations during loaded and inclined walking. PLoS ONE 10:e0121407. doi:10.1371/journal.pone.0121407

PubMed Abstract | CrossRef Full Text | Google Scholar

Dundar, U., Toktas, H., Solak, O., Ulasli, A. M., and Eroglu, S. (2015). A comparative study of conventional physiotherapy versus robotic training combined with physiotherapy in patients with stroke. Top. Stroke Rehabil. 21, 453–461. doi:10.1310/tsr2106-453

CrossRef Full Text | Google Scholar

Farris, D. J., Hicks, J. L., Delp, S. L., and Sawicki, G. S. (2014). Musculoskeletal modelling deconstructs the paradoxical effects of elastic ankle exoskeletons on plantar-flexor mechanics and energetics during hopping. J. Exp. Biol. 217, 4018–4028. doi:10.1242/jeb.107656

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferreira, L. A. B., Neto, H. P., Grecco, L. A. C., Christovão, T. C. L., Duarte, N. A., Lazzari, R. D., et al. (2013). Effect of ankle-foot orthosis on gait velocity and cadence of stroke patients: a systematic review. J. Phys. Ther. Sci. 25, 1503–1508. doi:10.1589/jpts.25.1503

PubMed Abstract | CrossRef Full Text | Google Scholar

Fleischer, C., and Hommel, G. (2008). A human – exoskeleton interface utilizing electromyography. IEEE Trans. Robot. 24, 872–882. doi:10.1109/TRO.2008.926860

CrossRef Full Text | Google Scholar

Fregly, B. J., Fregly, C. D., and Kim, B. T. (2015). Computational prediction of muscle moments during ared squat exercise on the international space station. J. Biomech. Eng. 137, 121005. doi:10.1115/1.4031795

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerritsen, K. G., van den Bogert, A. J., Hulliger, M., and Zernicke, R. F. (1998). Intrinsic muscle properties facilitate locomotor control – a computer simulation study. Motor Control 2, 206–220. doi:10.1123/mcj.2.3.206

PubMed Abstract | CrossRef Full Text | Google Scholar

Geyer, H., and Herr, H. (2010). A muscle-reflex model that encodes principles of legged mechanics produces human walking dynamics and muscle activities. IEEE Trans. Neural Syst. Rehabil. Eng. 18, 263–273. doi:10.1109/TNSRE.2010.2047592

PubMed Abstract | CrossRef Full Text | Google Scholar

Gill, P. E., Murray, W., and Saunders, M. A. (2005). SNOPT: an SQP algorithm for large-scale constrained optimization. SIAM Rev. 47, 99–131. doi:10.1137/S0036144504446096

CrossRef Full Text | Google Scholar

Go, A. S., Mozaffarian, D., Roger, V. L., Benjamin, E. J., Berry, J. D., Borden, W. B., et al. (2013). Executive summary: heart disease and stroke statistics – 2013 update: a report from the American Heart Association. Circulation 127, 143–152. doi:10.1161/CIR.0b013e318282ab8f

CrossRef Full Text | Google Scholar

Goldberg, E. J., and Neptune, R. R. (2007). Compensatory strategies during normal walking in response to muscle weakness and increased hip joint stiffness. Gait Posture 25, 360–367. doi:10.1016/j.gaitpost.2006.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamner, S. R., Seth, A., and Delp, S. L. (2010). Muscle contributions to propulsion and support during running. J. Biomech. 43, 2709–2716. doi:10.1016/j.jbiomech.2010.06.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Handsfield, G. G., Meyer, C. H., Hart, J. M., Abel, M. F., and Blemker, S. S. (2014). Relationships of 35 lower limb muscles to height and body mass quantified using MRI. J. Biomech. 47, 631–638. doi:10.1016/j.jbiomech.2013.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

He, J., Levine, W. S., and Loeb, G. E. (1991). Feedback gains for correcting small perturbations to standing posture. IEEE Trans. Automat. Contr. 36, 322–332. doi:10.1109/9.73565

CrossRef Full Text | Google Scholar

Heilman, B. P., and Kirsch, R. F. (2003). “Selection of an optimal muscle set for a standing neuroprosthesis using a human musculoskeletal model,” in Proceedings of the 25th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (IEEE Cat. No.03CH37439) (Cancun: IEEE), 1535–1538. doi:10.1109/IEMBS.2003.1279645.

PubMed Abstract | CrossRef Full Text | Google Scholar

Hug, F. (2011). Can muscle coordination be precisely studied by surface electromyography? J. Electromyogr. Kinesiol. 21, 1–12. doi:10.1016/j.jelekin.2010.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunt, K. H., and Crossley, F. R. E. (1975). Coefficient of restitution interpreted as damping in vibroimpact. J. Appl. Mech. 42, 440. doi:10.1115/1.3423596

CrossRef Full Text | Google Scholar

Hussain, S. (2014). State-of-the-art robotic gait rehabilitation orthoses: design and control aspects. NeuroRehabilitation 35, 701–709. doi:10.3233/NRE-141174

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanenko, Y. P., Cappellini, G., Dominici, N., Poppele, R. E., and Lacquaniti, F. (2005). Coordination of locomotion with voluntary movements in humans. J. Neurosci. 25, 7238–7253. doi:10.1523/JNEUROSCI.1327-05.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Jackson, J. N., Hass, C. J., and Fregly, B. J. (2016). Development of a subject-specific foot-ground contact model for walking. J. Biomech. Eng. 138:091002. doi:10.1115/1.4034060

PubMed Abstract | CrossRef Full Text | Google Scholar

Jordan, K., Challis, J. H., and Newell, K. M. (2007). Walking speed influences on gait cycle variability. Gait Posture 26, 128–134. doi:10.1016/j.gaitpost.2006.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Kane, T. R., and Levinson, D. A. (1985). Dynamics: Theory and Applications. New York: McGraw-Hill.

Google Scholar

Kesar, T. M., Perumal, R., Jancosko, A., Reisman, D. S., Rudolph, K. S., Higginson, J. S., et al. (2010). Novel patterns of functional electrical stimulation have an immediate effect on dorsiflexor muscle function during gait for people poststroke. Phys. Ther. 90, 55–66. doi:10.2522/ptj.20090140

PubMed Abstract | CrossRef Full Text | Google Scholar

Kesar, T. M., Perumal, R., Reisman, D. S., Jancosko, A., Rudolph, K. S., Higginson, J. S., et al. (2009). Functional electrical stimulation of ankle plantarflexor and dorsiflexor muscles: effects on poststroke gait. Stroke 40, 3821–3827. doi:10.1161/STROKEAHA.109.560375

PubMed Abstract | CrossRef Full Text | Google Scholar

Kia, M., Stylianou, A. P., and Guess, T. M. (2014). Evaluation of a musculoskeletal model with prosthetic knee through six experimental gait trials. Med. Eng. Phys. 36, 335–344. doi:10.1016/j.medengphy.2013.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Knarr, B. A., Kesar, T. M., Reisman, D. S., Binder-Macleod, S. A., and Higginson, J. S. (2013). Changes in the activation and function of the ankle plantar flexor muscles due to gait retraining in chronic stroke survivors. J. Neuroeng. Rehabil. 10, 12. doi:10.1186/1743-0003-10-12

PubMed Abstract | CrossRef Full Text | Google Scholar

Knarr, B. A., Reisman, D. S., Binder-Macleod, S. A., and Higginson, J. S. (2014). Changes in predicted muscle coordination with subject-specific muscle parameters for individuals after stroke. Stroke Res. Treat. 2014, 321747. doi:10.1155/2014/321747

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobayashi, T., Orendurff, M. M. S., Singer, M. L. M., Gao, F., Daly, W. K., and Foreman, K. B. (2016). Reduction of genu recurvatum through adjustment of plantarflexion resistance of an articulated ankle-foot orthosis in individuals post-stroke. Clin. Biomech. 35, 81–85. doi:10.1016/j.clinbiomech.2016.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, D. D., and Seung, H. S. (1999). Learning the parts of objects by non-negative matrix factorization. Nature 401, 788–791. doi:10.1038/44565

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, H.-J., Cho, K.-H., and Lee, W.-H. (2013). The effects of body weight support treadmill training with power-assisted functional electrical stimulation on functional movement and gait in stroke patients. Am. J. Phys. Med. Rehabil. 92, 1051–1059. doi:10.1097/PHM.0000000000000040

PubMed Abstract | CrossRef Full Text | Google Scholar

Limebeer, D. J. N., and Rao, A. V. (2015). Faster, higher, and greener: vehicular optimal control. IEEE Control Syst. 35, 36–56. doi:10.1109/MCS.2014.2384951

CrossRef Full Text | Google Scholar

Lloyd, D. G., and Besier, T. F. (2003). An EMG-driven musculoskeletal model to estimate muscle forces and knee joint moments in vivo. J. Biomech. 36, 765–776. doi:10.1016/S0021-9290(03)00010-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Lloyd-Jones, D., Adams, R. J., Brown, T. M., Carnethon, M., Dai, S., De Simone, G., et al. (2010). Heart disease and stroke statistics – 2010 update: a report from the American Heart Association. Circulation 121, e46–e215. doi:10.1161/CIRCULATIONAHA.109.192667

CrossRef Full Text | Google Scholar

Malone, L. A., and Bastian, A. J. (2014). Spatial and temporal asymmetries in gait predict split-belt adaptation behavior in stroke. Neurorehabil. Neural Repair 28, 230–240. doi:10.1177/1545968313505912

PubMed Abstract | CrossRef Full Text | Google Scholar

McGowan, C. P., Neptune, R. R., Clark, D. J., and Kautz, S. A. (2010). Modular control of human walking: adaptations to altered mechanical demands. J. Biomech. 43, 412–419. doi:10.1016/j.jbiomech.2009.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehrholz, J., Elsner, B., Werner, C., Kugler, J., and Pohl, M. (2013). Electromechanical-assisted training for walking after stroke: updated evidence. Stroke 44, e127–e128. doi:10.1161/STROKEAHA.113.003061

CrossRef Full Text | Google Scholar

Menegaldo, L. L., de Toledo Fleury, A., and Weber, H. I. (2004). Moment arms and musculotendon lengths estimation for a three-dimensional lower-limb model. J. Biomech. 37, 1447–1453. doi:10.1016/j.jbiomech.2003.12.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer, A. J., Patten, C., and Fregly, B. J. (2016). Lower extremity EMG-driven modeling of walking with automated adjustment of musculoskeletal geometry. PLoS ONE (in review).

Google Scholar

Millard, M., Uchida, T., Seth, A., and Delp, S. L. (2013). Flexing computational muscle: modeling and simulation of musculotendon dynamics. J. Biomech. Eng. 135, 021005. doi:10.1115/1.4023390

PubMed Abstract | CrossRef Full Text | Google Scholar

Mooney, L. M., and Herr, H. M. (2016). Biomechanical walking mechanisms underlying the metabolic reduction caused by an autonomous exoskeleton. J. Neuroeng. Rehabil. 13, 4. doi:10.1186/s12984-016-0111-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Mutikainen, S., Rantanen, T., Alén, M., Kauppinen, M., Karjalainen, J., Kaprio, J., et al. (2011). Walking ability and all-cause mortality in older women. Int. J. Sports Med. 32, 216–222. doi:10.1055/s-0030-1268506

PubMed Abstract | CrossRef Full Text | Google Scholar

Nekoukar, V., and Erfanian, A. (2013). Dynamic optimization of walker-assisted FES-activated paraplegic walking: simulation and experimental studies. Med. Eng. Phys. 35, 1659–1668. doi:10.1016/j.medengphy.2013.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ng, M. F., Tong, R. K., and Li, L. S. (2008). A pilot study of randomized clinical controlled trial of gait training in subacute stroke patients with partial body-weight support electromechanical gait trainer and functional electrical stimulation: six-month follow-up. Stroke 39, 154–160. doi:10.1161/STROKEAHA.107.495705

PubMed Abstract | CrossRef Full Text | Google Scholar

Nilsson, A., Vreede, K. S., Häglund, V., Kawamoto, H., Sankai, Y., and Borg, J. (2014). Gait training early after stroke with a new exoskeleton – the hybrid assistive limb: a study of safety and feasibility. J. Neuroeng. Rehabil. 11, 92. doi:10.1186/1743-0003-11-92

CrossRef Full Text | Google Scholar

O’Dell, M. W., Dunning, K., Kluding, P., Wu, S. S., Feld, J., Ginosian, J., et al. (2014). Response and prediction of improvement in gait speed from functional electrical stimulation in persons with poststroke drop foot. PM R 6, 587–601;quiz601. doi:10.1016/j.pmrj.2014.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Olney, S., Monga, T. N., and Costigan, P. A. (1986). Mechanical energy of walking of stroke patients. Arch. Phys. Med. Rehabil. 67, 92–98. doi:10.1016/0003-9993(86)90109-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Ostir, G. V., Berges, I. M., Kuo, Y.-F., Goodwin, J. S., Fisher, S. R., and Guralnik, J. M. (2013). Mobility activity and its value as a prognostic indicator of survival in hospitalized older adults. J. Am. Geriatr. Soc. 61, 551–557. doi:10.1111/jgs.12170

PubMed Abstract | CrossRef Full Text | Google Scholar

Patterson, M. A., and Rao, A. V. (2014). GPOPS-II: a MATLAB software for solving multiple-phase optimal control problems using hp-adaptive gaussian quadrature collocation methods and sparse nonlinear programming. ACM Trans. Math. Softw. 41, 1–37. doi:10.1145/2558904

CrossRef Full Text | Google Scholar

Pennycott, A., Wyss, D., Vallery, H., Klamroth-Marganska, V., and Riener, R. (2012). Towards more effective robotic gait training for stroke rehabilitation: a review. J. Neuroeng. Rehabil. 9, 65. doi:10.1186/1743-0003-9-65

PubMed Abstract | CrossRef Full Text | Google Scholar

Pilkar, R., Yarossi, M., and Nolan, K. J. (2014). EMG of the tibialis anterior demonstrates a training effect after utilization of a foot drop stimulator. NeuroRehabilitation 35, 299–305. doi:10.3233/NRE-141126

PubMed Abstract | CrossRef Full Text | Google Scholar

Popovic, D., Stein, R. B., Namik Oguztoreli, M., Lebiedowska, M., and Jonic, S. (1999). Optimal control of walking with functional electrical stimulation: a computer simulation study. IEEE Trans. Rehabil. Eng. 7, 69–79. doi:10.1109/86.750554

PubMed Abstract | CrossRef Full Text | Google Scholar

Racinais, S., Maf, N. A., and Girard, O. (2013). M-wave, H- and V-reflex recruitment curves during maximal voluntary contraction. J. Clin. Neurophysiol. 30, 415–421. doi:10.1097/WNP.0b013e31829ddcf1

PubMed Abstract | CrossRef Full Text | Google Scholar

Reinbolt, J. A., Haftka, R. T., Chmielewski, T. L., and Fregly, B. J. (2008). A computational framework to predict post-treatment outcome for gait-related disorders. Med. Eng. Phys. 30, 434–443. doi:10.1016/j.medengphy.2007.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Reinbolt, J. A., Schutte, J. F., Fregly, B. J., Koh, B. I., Haftka, R. T., George, A. D., et al. (2005). Determination of patient-specific multi-joint kinematic models through two-level optimization. J. Biomech. 38, 621–626. doi:10.1016/j.jbiomech.2004.03.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Reisman, D. S., Wityk, R., Silver, K., and Bastian, A. J. (2007). Locomotor adaptation on a split-belt treadmill can improve walking symmetry post-stroke. Brain 130, 1861–1872. doi:10.1093/brain/awm035

PubMed Abstract | CrossRef Full Text | Google Scholar

Riener, R. (1999). Model-based development of neuroprosthesis for paraplegic patients. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 354, 877–894. doi:10.1098/rstb.1999.0440

PubMed Abstract | CrossRef Full Text | Google Scholar

Risher, D. W., Schutte, L. M., and Runge, C. F. (1997). The use of inverse dynamics solutions in direct dynamics simulations. J. Biomech. Eng. 119, 417. doi:10.1115/1.2798288

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth, E. J., Merbitz, C., Mroczek, K., Dugan, S. A., and Suh, W. W. (1997). Hemiplegic gait: relationships between walking speed and other temporal parameters. Am. J. Phys. Med. Rehabil. 76, 128–133. doi:10.1097/00002060-199703000-00008

PubMed Abstract | CrossRef Full Text | Google Scholar

Sabut, S. K., Bhattacharya, S. D., and Manjunatha, M. (2013). Functional electrical stimulation on improving foot drop gait in poststroke rehabilitation: a review of its technology and clinical efficacy. Crit. Rev. Biomed. Eng. 41, 149–160. doi:10.1615/CritRevBiomedEng.2013007621

PubMed Abstract | CrossRef Full Text | Google Scholar

Sartori, M., Reggiani, M., van den Bogert, A. J., and Lloyd, D. G. (2012a). Estimation of musculotendon kinematics in large musculoskeletal models using multidimensional B-splines. J. Biomech. 45, 595–601. doi:10.1016/j.jbiomech.2011.10.040

CrossRef Full Text | Google Scholar

Sartori, M., Reggiani, M., Farina, D., and Lloyd, D. G. (2012b). EMG-driven forward-dynamic estimation of muscle force and joint moment about multiple degrees of freedom in the human lower extremity. PLoS ONE 7:e52618. doi:10.1371/journal.pone.0052618

CrossRef Full Text | Google Scholar

Sawicki, G., and Khan, N. (2015). A simple model to estimate plantarflexor muscle-tendon mechanics and energetics during walking with elastic ankle exoskeletons. IEEE Trans. Biomed. Eng. 63, 914–923. doi:10.1109/TBME.2015.2491224

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, Q., and Buchanan, T. S. (2008). A biomechanical model to estimate corrective changes in muscle activation patterns for stroke patients. J. Biomech. 41, 3097–3100. doi:10.1016/j.jbiomech.2008.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharif Razavian, R., Mehrabi, N., and McPhee, J. (2015). A model-based approach to predict muscle synergies using optimization: application to feedback control. Front. Comput. Neurosci. 9:121. doi:10.3389/fncom.2015.00121

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, N., Mushahwar, V., and Stein, R. (2014). Dynamic optimization of FES and orthosis-based walking using simple models. IEEE Trans. Neural Syst. Rehabil. Eng. 22, 114–126. doi:10.1109/TNSRE.2013.2280520

PubMed Abstract | CrossRef Full Text | Google Scholar

Silder, A., Whittington, B., Heiderscheit, B., and Thelen, D. G. (2007). Identification of passive elastic joint moment-angle relationships in the lower extremity. J. Biomech. 40, 2628–2635. doi:10.1016/j.jbiomech.2006.12.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Silverman, A. K., and Neptune, R. R. (2012). Muscle and prosthesis contributions to amputee walking mechanics: a modeling study. J. Biomech. 45, 2271–2278. doi:10.1016/j.jbiomech.2012.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Simonsen, E. B., Alkjaer, T., and Raffalt, P. C. (2012). Reflex response and control of the human soleus and gastrocnemius muscles during walking and running at increasing velocity. Exp. Brain Res. 219, 163–174. doi:10.1007/s00221-012-3075-y

PubMed Abstract | CrossRef Full Text | Google Scholar

States, R. A., Pappas, E., and Salem, Y. (2009). Overground physical therapy gait training for chronic stroke patients with mobility deficits. Stroke 40, e627–e628. doi:10.1161/STROKEAHA.109.558940

PubMed Abstract | CrossRef Full Text | Google Scholar

Thelen, D. G. (2003). Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults. J. Biomech. Eng. 125, 70. doi:10.1115/1.1531112

PubMed Abstract | CrossRef Full Text | Google Scholar

Thelen, D. G., Anderson, F. C., and Delp, S. L. (2003). Generating dynamic simulations of movement using computed muscle control. J. Biomech. 36, 321–328. doi:10.1016/S0021-9290(02)00432-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Tresch, M. C., Saltiel, P., and Bizzi, E. (1999). The construction of movement by the spinal cord. Nat. Neurosci. 2, 162–167. doi:10.1038/5721

PubMed Abstract | CrossRef Full Text | Google Scholar

Tyson, S. F., Sadeghi-Demneh, E., and Nester, C. J. (2013). A systematic review and meta-analysis of the effect of an ankle-foot orthosis on gait biomechanics after stroke. Clin. Rehabil. 27, 879–891. doi:10.1177/0269215513486497

PubMed Abstract | CrossRef Full Text | Google Scholar

van den Bogert, A. J., Blana, D., and Heinrich, D. (2011). Implicit methods for efficient musculoskeletal simulation and optimal control. Procedia IUTAM 2, 297–316. doi:10.1016/j.piutam.2011.04.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Winstein, C. J., Stein, J., Arena, R., Bates, B., Cherney, L. R., Cramer, S. C., et al. (2016). Guidelines for adult stroke rehabilitation and recovery: a guideline for healthcare professionals from the American Heart Association/American Stroke Association. Stroke 47, e98–e169. doi:10.1161/STR.0000000000000098

PubMed Abstract | CrossRef Full Text | Google Scholar

World Stroke Organization. (2016). Available at: http://www.world-stroke.org/advocacy/world-stroke-campaign [Accessed January 1, 2016].

Google Scholar

Zhang, D., and Zhu, K. (2007). Modeling biological motor control for human locomotion with functional electrical stimulation. Biol. Cybern. 96, 79–97. doi:10.1007/s00422-006-0107-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Zmitrewicz, R. J., Neptune, R. R., and Sasaki, K. (2007). Mechanical energetic contributions from individual muscles and elastic prosthetic feet during symmetric unilateral transtibial amputee walking: a theoretical study. J. Biomech. 40, 1824–1831. doi:10.1016/j.jbiomech.2006.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar