An Optimization Method Tracking EMG, Ground Reactions Forces, and Marker Trajectories for Musculo-Tendon Forces Estimation in Equinus Gait

Introduction

Walking is often considered to be the most important activity in daily living (Chiou et al., 1985). The ability to move without pain, fatigue, or major gait deviation is closely related to quality of life (Cuomo et al., 2007; van Schie, 2008). Many neuro-orthopedic pathologies (e.g., cerebral palsy, stroke) induce impairments (i.e., paresis, muscle overactivity, soft tissue contractures, and bone deformities) that compromise normal movement. Consequently, the goal of many therapeutic interventions is to minimize gait deviations in patients. In order to improve these interventions, understanding the possible mechanisms of these gait deviations, and being able to identify the causal motor impairments is of primary importance (Davids et al., 2004; Gough and Shortland, 2008; Wren et al., 2011). Currently, the relationship between motor impairments and gait deviations is unclear (Bonnefoy-Mazure et al., 2016), and there is a lack scientific evidence for these relationships due to the inherent complexity of the human neuromusculoskeletal system during dynamic tasks such as walking (Armand et al., 2017). Compared to existing approaches (e.g., pathologic models, experimental procedures with human subjects, robots with human-like gait), in silico neuromusculoskeletal simulations of normal and pathological gait could provide additional insight into gait deviations (Armand et al., 2017). The advantage of neuromuscular simulations as a method is that large numbers of simulations may be performed relatively quickly, and without the ethical issues involved with performing invasive and lengthy experiments in vulnerable patient groups such as those with neuromuscular deficits.

To date, simulations reported in the literature are often limited to the analysis of the consequences of isolated impairments on gait such as muscle weakness (van der Krogt et al., 2009; Thompson et al., 2013) or muscle spasticity (Jansen et al., 2014). Most of these studies only report possible muscular compensations (adaptations) that occur due to muscular redundancy–results are achieved by tracking the normal gait kinematics and then applying a perturbation to the model. Very few studies have been based on a numerical framework allowing kinematic adaptations in response to more varied perturbations. Within the TLEMsafe project (https://www.tlemsafe.eu), Fluit et al. (2014b) combined an optimized inverse model and a ground reaction force predictive model to simulate lower limb kinematics after the removal of the rectus femoris and the vastus lateralis from the model. Santos et al. (2017) also proposed a numerical framework based on direct collocation and an optimal control package to simulate lower limb kinematics after the introduction of a weakening of the triceps surae and the tibialis anterior, or after increasing the ankle joint stiffness. These two approaches represent a first step toward the simulation of pathological gait. However, neither were able to reproduce muscle co-contractions. While this capacity was not necessarily needed in these studies, this feature is essential to establish a pathological gait simulator that would be able to reproduce physiological gait adaptations biofidelically.

From a methodological point of view, inverse dynamics-based approaches (such as static optimization) are commonly used due to their computational efficiency (Erdemir et al., 2007), but are not appropriate for predictive simulations. Moreover, static optimization method underestimates or neglects antagonist co-contractions unless hybrid approaches are used (Brookham et al., 2011; Son et al., 2012). On the other hand, forward dynamics-based approaches are often criticized for being time-consuming–several studies report convergence times in the hundreds of hours range (Anderson and Pandy, 2001). Despite this disadvantage, these methods have the potential to predict new movements, such as an adaptation in response to a perturbation. For example, state-of-the-art algorithms used in conjunction with existing musculoskeletal models–like direct collocation (Santos et al., 2017) and direct multiple shooting (Bélaise et al., 2018a,b)—can be used to solve forward dynamics problems in a timely manner. Recently, Bélaise et al. (2018a,b) introduced an EMG-marker tracking optimization method to predict musculo-tendon forces in a co-contraction case. Based on simulated datasets of upper limb movements, the authors showed the importance of tracking both marker trajectories and EMG, in particular to reproduce muscle co-contractions. To the best of our knowledge, such an approach has never been applied on experimental gait records with muscle co-contractions.

The objective of our project is to establish a computational framework appropriate for predictive simulations of healthy and pathological gait, that is able to reproduce muscle co-contractions, and that allows for the introduction of various kind of perturbations on the model (e.g., therapy-related, surgery-related, pathology-related perturbations). This preliminary study represents a first step toward this project by adapting the computational framework proposed by Bélaise et al. (2018a,b) to a lower limb musculoskeletal model during gait. This framework has been tested for this purpose on a type of pathological gait known as equinus gait.

Methods

Lower Limb Musculoskeletal Model

A generic three-dimensional musculoskeletal model of the lower limb [Lower Extremity Model, OpenSim (Delp et al., 1990)] was adapted for our study (Figure 1). This model consists of five rigid segments: the pelvis, right thigh, patella, shank, and foot. Twenty-six markers were associated with these segments by virtual palpation to reproduce the experimental marker locations (Table 1; see section Dataset). To simplify the dynamic optimizations to a two-dimensional motion in this preliminary study, the original degrees of freedom (DoF) were reduced to three DoFs for the pelvis-ground joint (vertical translation, translation in the direction of walking, pelvis tilt) and one DoF (flexion-extension modeled as a hinge joint) at the hip, knee, and ankle joints. Joints were actuated by the muscle torques resulting from 17 muscle lines of action (Table 2), and the pelvis DoFs were actuated by three generalized forces applied on the pelvis. The path, optimal length, maximal isometric force, tendon slack length, and pennation angle of each muscle lines of action were derived from the original model (Delp et al., 1990).

FIGURE 1

Figure 1. Anterior (A), lateral (B), and posterior (C) views of the right lower-limb musculoskeletal model derived from OpenSim [Lower Extremity Model, OpenSim (Delp et al., 1990)] and adapted to the bioRBD musculoskeletal modeling package (https://github.com/pyomeca/biorbd). Red lines and pink dots represents the 17 Hill-type muscle lines of action and the 26 markers related to this model, respectively.

TABLE 1

Table 1. List of the 26 markers used in this study.

TABLE 2

Table 2. List of the 17 Hill-type muscle lines of action included in the model.

Segment lengths were scaled to the dataset used in this study (see section Dataset) using OpenSim 3.3 (Delp et al., 2007) by minimizing the distance between experimental and model markers placed on bony landmarks (Table 1). All components of the model that depend on bone lengths (e.g., muscle attachment points, optimal fiber length), segment masses, and inertial parameters were also scaled. The resulting scaled model was transferred to the bioRBD musculoskeletal modeling package (https://github.com/pyomeca/biorbd) based on the Rigid Body Dynamic Library (Felis, 2017). This model was defined by 29 states (six generalized joint positions and their six related velocities, 17 muscle activations) and 20 controls corresponding to the 17 muscle neural excitations plus three generalized forces driving the three pelvis DoFs.

Equations of Motion and Activation Dynamics

The generalized accelerations $\ddot{\text{q}}$ of the rigid multibody system were computed using a forward dynamics approach for given generalized joint positions q, joint velocities $\stackrel{\bullet }{\text{q}}$, and generalized forces τ:

$\begin{array}{l}\stackrel{\mathrm{..}}{q}=M{\left(q\right)}^{-1}\left(\tau \left(q,\stackrel{.}{q},a,e\right)+C{\left(q\right)}^{T}R-N\left(q,\stackrel{.}{q}\right)\stackrel{.}{q}-G\left(q\right)\right)\\ s.t.C\left(q\right)\stackrel{\mathrm{..}}{q}+\stackrel{.}{C}\left(q\right)\stackrel{.}{q}=0andC\left(q\right)\stackrel{.}{q}=0\end{array}$

where M is the inertia matrix, C is the external contact Jacobian matrix, R is the Lagrange multipliers vector corresponding to the ground reaction forces (GRF), N is the non-linear effects (Coriolis and centrifugal forces) vector, and G is the gravity effects vector. It was assumed that contact points have a null acceleration and velocity throughout the entire contact phase. In line with equinus gait, one fixed contact point was defined on the forefoot for the entire contact phase (see section Dynamic Optimizations). Generalized forces were divided into ${\tau }_{\mathbf{\text{1}}}={\left[{\tau }_{\mathbf{\text{11}}}{\tau }_{\mathbf{\text{12}}}{\tau }_{\mathbf{\text{13}}}\right]}^T$ driving the three pelvis DoFs, and ${\tau }_{\mathbf{\text{2}}}=\frac{\partial {\mathbf{\text{L}}}_{\mathbf{\text{mt}}}}{\partial \mathbf{\text{q}}}{\mathbf{\text{F}}}_{\mathbf{\text{mt}}}$ corresponding to the net joint torques due to the musculo-tendon forces F_mt, where $\frac{\partial {\mathbf{\text{L}}}_{\mathbf{\text{mt}}}}{\partial \mathbf{\text{q}}}$ is the generalized muscular lever arms matrix and L_mt the vector of muscle line of action lengths. F_mt were computed from muscle activations a using a Hill-type muscle model with a generic force-length-velocity relation f (Zajac, 1989):

$F_{mt}\left(q,\stackrel{.}{q},a\right)=af\left(F_{mt}^0,l_m,v_m\right),$

where ${\mathbf{\text{F}}}_{\mathbf{\text{mt}}}^{\mathbf{\text{0}}}$ is the maximal isometric forces vector, l_m is the muscle fiber lengths vector, and v_m is the muscle fiber velocities vector. Again, musculo-tendon forces were divided into F_mt1 (with related activations a₁ and excitations e₁), corresponding to the muscles for which electromyographic (EMG) records were available, and F_mt2 (with related activations a₂ and excitations e₂) where EMG measurements were unavailable. Muscle activation dynamics was implemented as a set of first-order differential equations (Buchanan et al., 2004):

$\stackrel{.}{a}\left(t,e\left(t\right),a(t)\right)=\{\begin{array}{c}\frac{\left(e(t)-a\left(t\right)\right)}{t_{act}\left(0.5+1.5a(t)\right)},e\left(t\right)>a(t)\\ \frac{e(t)-a(t)}{t_{deact}}\left(0.5+1.5a(t)\right),e\left(t\right)\le a(t)\end{array}$

where e(t) is the muscle neural excitations at time t. Time constants t_act and t_deact (for activation and deactivation) were set at 10 and 40 ms, respectively (Thelen et al., 2003).

Dynamic Optimizations

As proposed in Bélaise et al. (2018a), controls and state variables were simultaneously optimized using an EMG-marker tracking optimization process. Because the model has a reduced muscular redundancy and because a generic model was used as opposed to a subject-specific model, the optimal maximal isometric forces were also identified during this process.

The optimization consisted of the minimization of the differences between predicted M_p and measured M_m marker trajectories in the sagittal plane, and between predicted e_1p and measured e_1m (EMG envelop) muscle neural excitations (corresponding to the muscles for which EMG records are available). This tracking optimization was extended to the minimization of the differences between predicted R_p and measured R_m GRF in the sagittal plane. This tracking was necessary to impose physiologic generalized forces to the pelvis (τ₁), i.e., generalized forces that compensate for the missing upper part of the body and contralateral lower limb.

To predict the activity of the muscles for which tracking was not possible (i.e., muscles for which EMG records were not available), the objective function J was written to find the least squared muscle activations a₂ that produced the prescribed marker trajectories, muscle neural excitations, and GRF (during stance phase only):

$\begin{array}{l}J={\displaystyle \sum _1^{N_i}}\left(w_{\text{M}}\left|\right|{\text{M}}_p-{\text{M}}_m|{|}^2+w_e||{\text{e}}_{1p}-{\text{e}}_{1m}|{|}^2+w_{\text{R}}\left|\right|{\text{R}}_p-{\text{R}}_m|{|}^2\right)\\ +w_L{\displaystyle {\int }_0^{T_i}}{{\text{a}}_2(t)}^{2}dt\\ where{w}_{\text{R}}\left|\right|{\text{R}}_p-{\text{R}}_m|{|}^2=0wheni=2,i.e.,theswingphase.\end{array}$

where w_M, w_e, w_R, and w_L are weighting factors adjusted to the relative importance of each term, T_i is the duration of the current stage (see section Simulations) and N_i is the related number of time frames.

This objective function was minimized under three sets of constraints. First, boundary conditions were applied on the state and the control variables. In this study, the range of motion of each DoF and related velocities were set to physiologic values (Table 3), while activations and excitations were bounded between 0 and 1. Second, the velocity of the contact point was constrained to be null at the first frame and its acceleration to be null at each time frame (see section Equations of Motion and Activation Dynamics). Third, periodicity was ensured by constraining the first and last time point of the cycle to have similar values in terms of hip, knee, ankle joint angles, and velocities, pelvis velocities, muscle excitations, and GRF.

TABLE 3

Table 3. Boundaries constraints applied during the optimization process to each degree of freedom and related velocities.

Simulations

Each dynamic optimization was solved using a direct multiple shooting algorithm with MUSCOD-II (Leineweber et al., 2003). Three phases were defined in the gait step: (1) the stance phase (with an external contact between foot and ground), (2) the swing phase (no external contact), and (3) the first frame of the next stance phase following the impact between the foot and the ground. These stages were divided into 25, 25, and 1 multiple shooting intervals, respectively. For the sake of simplicity, the first stage started just after the collision impact between the foot and the ground. The duration of each stage was fixed to the measured value.

The initial guess was set to the measured values for the joint positions and velocities, 1% for all activations and excitations, and 0 for the controls corresponding to the generalized forces related to the pelvis DoFs. Weighting factors were set to w_M = 30 (except for the foot markers for which w_M = 50 to ensure the correct position of the contact point), w_e = 1, w_R = 0.05, and w_L = 1. These weighting factors were adjusted empirically to set values around 1, in order to ensure optimization convergence and produce simulation results close to the experimentally measured data.

Dataset

The previously defined method was evaluated on a dataset of emulated equinus gait. All data were recorded on a healthy volunteer (male, 35 years old, 165 cm, 66 kg) without any neuro-orthopedic conditions. This participant gave written informed consent prior to his inclusion and the protocol was conformed to the Declaration of Helsinki and approved by the National Research Ethics Committee of Luxembourg (201805/01).

The 3D trajectories of 26 reflective cutaneous markers (bilateral iliac anterior and posterior spines, right leg great trochanter, medial and lateral femoral epicondyles, peroneal head, tibial tuberosity, medial and lateral malleoli, 1st, 2nd, and 5th proximal and distal metatarsal heads, calcaneum, completed by a four-marker cluster on the thigh and on the shank) (Figure 1) were recorded using a 10-camera optoelectronic system (OQUS-4, Qualisys AB, Sweden) sampled at 200 Hz. Markers were placed by anatomical palpation (Table 1) following the recommendation of van Sint Jan (2007) by an experienced user. GRF and moments were recorded using two side-by-side force plates (OR6-5, AMTI, USA) sampled at 2,000 Hz. The EMG activity of nine right leg muscles (tibialis anterior, soleus, gastrocnemius medialis, vastus medialis, rectus femoris, semitendinosus, biceps femoris long head, gluteus medius, gluteus maximus) was collected with a wireless electromyographic system (DTS clinic, Noraxon, USA) sampled at 2,000 Hz. The EMG surface electrodes were placed following the recommended standard of the Surface EMG for a Non-Invasive Assessment of Muscles (SENIAM) project (Hermens et al., 2000).

All data were imported under Matlab (R2018a, The MathWorks, USA) using the ezc3d package (https://github.com/pyomeca/ezc3d). Marker trajectories were interpolated when necessary using a cubic spline and smoothed by a 4th order low-pass Butterworth filter with a cutoff frequency of 6 Hz. Generalized kinematics ($q,\stackrel{\bullet }{q},\ddot{q}$) were computed using an extended Kalman filter (Fohanno et al., 2014) following the segmental coordinate systems defined in the original generic three-dimensional lower limb musculoskeletal model [Lower Extremity Model, OpenSim (Delp et al., 1990)]. GRF were smoothed by a 4th order low-pass Butterworth filter with a cut-off frequency of 15 Hz. Raw EMG signals were band pass filtered (4th order) between 30 and 300 Hz, rectified, and EMG envelops were obtained by a 4th order low-pass Butterworth filter with a cut-off frequency of 25 Hz. EMG envelops were then normalized to their respective maximal voluntary activation (Gaudet et al., 2018).

The participant was asked to mimic an equinus gait by producing voluntarily controlled co-contractions of the muscles crossing the ankle joint to restrain ankle dorsiflexion. Eight trials were recorded and the related right steps were analyzed in this study.

Analysis

In order to evaluate the capacity of the model to reproduce the measured gait pattern and muscle excitations under the mechanical properties and constraints imposed to the model, a set of goodness-of-fit parameters were employed. Root mean square error (RMSE) and coefficient of determination (R²) were computed to assess the differences in intensity and shape, respectively, between measured and estimated excitations, joint angles and GRF.

Only estimated muscle excitations corresponding to the measured EMG envelops (Table 2) are presented in this analysis. The coefficient of determination (CC) (Giroux et al., 2013) was computed for the muscles for which EMG data was recorded. This method uses active/inactive state concordance between the estimated muscle excitations and normalized EMG envelopes to compute a coefficient of concordance defined as the percentage of concordance elements.

Results

The convergence time of the eight optimizations using MUSCOD-II was 3.88 ± 1.04 h on an Intel^® Core™ i5-3570 CPU @3.4 GHz. Estimated and measured muscle excitations, musculo-tendon forces, joint angles, and GRF are reported in Figures 2–5, respectively. Goodness-of-fit parameters (RMSE, R² and CC) are reported in Table 4.

FIGURE 2

Figure 2. Mean and standard deviation of normalized measured (EMG envelops) and estimated muscle excitations during gait cycle (EMG envelops have been adimensioned (adim) by maximal voluntary contraction). Abbreviations of muscle names are given in Table 2. For illustration purpose, EMG envelops of gluteus maximus, gluteus medius, semitendinosus, and vastus medialis are reported on plots R_GLUT_MAX1/R_GLUT_MAX2/R_GLUT_MAX3, R_GLUT_MED1/R_GLUT_MED2/R_GLUT_MED3, R_SEMIMEM/R_SEMITEN, and R_VAS_MED/R_VAS_INT/R_VAS_LAT, respectively.

FIGURE 3

Figure 3. Mean and standard deviation of estimated musculo-tendon forces during gait cycle. Abbreviations of muscle names are given in Table 2.

FIGURE 4

Figure 4. Mean and standard deviation of measured and estimated pelvic position/orientation and joint angles during gait cycle. Abbreviations of each degree of freedom are given in Table 3.

FIGURE 5

Figure 5. Mean and standard deviation of measured and estimated vertical (R_GRF_V) and anterior/posterior (R_GRF_AP) ground reaction forces during gait cycle.

TABLE 4

Table 4. Root mean square error (RMSE), coefficient of determination (R2) computed to assess the differences in intensity and shape, respectively, between measured and estimated excitations (adimensioned), pelvis position/orientation, joint angles, vertical ground reaction force (R_GRF_V), and anterior/posterior ground reaction force (R_GRF_AP).

Considering the tracked muscle excitations (i.e., muscles for which EMG records are available, see Table 2), the temporal muscle activity of the model was in good overall agreement with the experimental measurements, with an average CC of 76.8 ± 5.3%. RMSE values were generally low with an average value of 0.2 ± 0.1 (values were adimensioned between 0 and 1). However, RMSE was found to be higher for the gastrocnemius medialis (0.3 ± 0.1) and the tibialis anterior (0.4 ± 0.1). For all muscles, the correlation remained low with an average R² at 0.02 ± 0.52. Regarding all other model muscles, for which EMG was not tracked, the optimized muscle excitations were higher than those estimated for muscles for which EMG was tracked. This is the case for R_GLUT_MAX2 and R_GLUT_MAX3 compared to R_GLUT_MAX1, R_GLUT_MED1, and R_GLUT_MED3 compared to R_GLUT_MED2, R_SEMIMEM compared to R_SEMITEN, R_VAS_INT, and R_VAS_LAT compared to R_VAS_MED, and for R_GAS_LAT compared to R_GAS_MED. The same results are observed on the estimated musculo-tendon forces. These forces are ranged between 0 and 2,000 N, with the highest peak forces obtained for R_SEMITEN, R_VAST_INT, R_VAST_LAT, and R_SOLEUS.

With regard to the pelvis position/orientation and joint angles, the model estimations were generally in agreement with the experimental measurements. Average RMSE were 5.4 ± 3.3 mm for the pelvic translations, and 1.9 ± 1.3° for pelvic rotations and joint angles. However, RMSE was found to be higher for the ankle joint (4.0 ± 0.9°). Considering all degrees of freedom, the correlation remained very high with an average R² at 0.94 ± 0.09.

For the GRF, the model estimations were found to be in agreement with the experimental measurements (the average RMSE is 17.3 ± 5.5 N). For these forces, the correlation remained high with an average R² at 0.97 ± 0.03.

Discussion

The main objective of this study was to adapt an EMG marker tracking optimization process to solve a forward dynamics problem on a 3D musculoskeletal model of the lower limb during equinus gait. To the best of our knowledge, the use of a direct multiple shooting algorithm on a musculoskeletal model with the tracking of measured EMG, marker trajectories, and GRF has never been performed to date. As already shown by Bélaise et al. (2018a,b), this approach allows for an accurate reproduction of joint kinematics and ensures temporal fidelity in muscle activity with improved computational time compared to traditional forward dynamic approaches. We have also highlighted that the tracking of GRF could be performed accurately, even without the use of a complex foot/ground contact model.

Limitations

A primary limitation of this preliminary study is that it was based on a small number trials for a single task, performed by a single participant. As such, limited conclusions can be drawn from this paper.

A second limitation is that only one contact point was defined at the forefoot and it was only constrained to null velocity and acceleration during the whole contact phase. While this approach was in line with an equinus gait and was able to accurately reproduce the tracked GRF, this definition cannot be applied during normal gait trials, during which several contact points should be defined (Fluit et al., 2014a). Elastic contact elements (Peng et al., 2018), artificial muscle-like actuators (Fluit et al., 2014a) or distance and velocity-dependent force models (Jung et al., 2016) should be adapted to the present model to extend its use to normal gait.

The proposed musculoskeletal model also only consisted of the pelvis and the right lower limb. A forward dynamic approach was made possible by replacing the forces and moments produced by the opposite lower limb and the upper part of the body by a set of generalized forces acting on the pelvis. With such an approach, the individual contribution of the opposite lower limb and the upper part of the body to the muscle activity and joint contact forces of the right lower limb cannot be evaluated individually. It would thus be an important next step to complete the missing body segments of the present musculoskeletal model in order to obtain a full body musculoskeletal model. Several full body musculoskeletal models (Rajagopal et al., 2016; Bassani et al., 2017) have been proposed in the literature and could be transferred to the bioRBD musculoskeletal modeling package. These segments could be actuated by joint torques instead of muscles.

In addition to this, while we used a 3D musculoskeletal model, DoFs were reduced to only allow a two-dimensional motion in the sagittal plane. It is however, established that walking is a locomotion task that is performed in sagittal, coronal, and transversal planes (Perry and Burnfield, 2010). In particular, patients often develop compensatory movements in the coronal plane when pathological impairments result in reduced foot clearance capacity in the sagittal plane (Chantraine et al., 2016). Despite this simplification, the accuracy of kinematic tracking observed in the present results suggest that there is potential for the EMG marker tracking optimization process to perform 3D gait motion simulations.

Finally, most of the parameters of the Hill-type muscle model were kept with generic property definitions in this study. Only optimal maximal isometric forces were identified during dynamic optimizations, while muscle optimal lengths, tendon slack lengths, and maximal isometric forces are usually identified in similar studies (Sartori et al., 2014; Pizzolato et al., 2015). This may explain the high excitations observed in muscles contained in a group, for which tracking of excitation was applied to other muscles. The introduction of further muscle parameters could be implemented in future studies, and would only be expected to impact the convergence time of the simulations.

Muscle Activity

The present study support the results of Bélaise et al. (2018a,b), which demonstrated that EMG tracking could be an efficient way to reproduce measured muscle excitations in simulations. While some amplitude differences appeared for certain tracked muscles, the temporal activity was generally reproduced with good fidelity. This outcome is crucial to ensure the ability of the model to produce muscle co-contractions. Similar approaches have already been proposed in the literature–EMG-driven musculoskeletal models have also been used to accurately reproduce muscle excitation patterns observed on EMG records (Shao et al., 2009; Sartori et al., 2012). However, these models are constrained to have as many muscle lines of action as EMG available in the dataset. To overcome this limitation, some authors have proposed hybrid approaches that combine EMG-driven and static optimization methods (Lloyd and Besier, 2003; Moissenet et al., 2014; Sartori et al., 2014). The drawback with this strategy is that by minimizing the difference between the motor joint moments computed by EMG-driven and inverse dynamics, kinematics may not be accurately reproduced. In that sense, the EMG-marker tracking algorithm proposed by Bélaise et al. (2018a,b) is a novelty. As this method tracks joint kinematics based on marker trajectories rather than joint moments as in a forward dynamic approach, the error diffusion is minimized and the simulation outputs reproduce the experimental kinematics more faithfully.

In our trials of emulated equinus gait, co-contractions of the ankle dorsiflexors, and plantarflexors can be observed during early stance to stabilize the joint in this specific posture. It is interesting to observe that, while the gastrocnemius medialis and the soleus (muscles for which EMG records were tracked) were contracted during this phase, the gastrocnemius lateralis (a muscle for which EMG records were not measured) was not in our simulations. Although use of EMG is somewhat limited to available hardware, muscle locations and the signal quality, measurement should be prioritized toward muscles presumed to be active during the task being investigated. In our case, focusing on a greater number of muscles crossing the ankle joint would have brought more relevant information to the model. A similar recommendation has already been proposed by Sartori et al. (2014); these authors suggested to prioritize EMG use on muscles “that reflect the patient's non-physiological muscular behavior.”

Although the EMG-GRF-marker tracking algorithm was able to reproduce physiological muscle activity, two points must be considered. First, we observed that when EMG was not tracked, the optimized muscle excitations and musculo-tendon forces were higher than the ones estimated when muscle EMG was tracked. As pointed out in the limitations of the study, this over-estimation is perhaps related to the use of generic muscle model parameters. For example, if the parameters applied to a muscle group would tend to limit its capacity to produce a motor joint moment, a higher muscle excitation would be required to reproduce the experimental measurements. This effect is further exacerbated if the excitation of a muscle in this group is constrained to a low level in accordance with the experimental EMG tracking, as the excitation of the other muscles of the group will have to compensate for this reduced excitation in the other muscle. Second, due to equinus (results in increased plantarflexion), the capacity of the triceps surae to produce a plantarflexion moment is reduced (Delp et al., 1990), in particular during the push-off phase. Thus, hip flexor recruitment may be increased in this gait pattern to pull the leg forward (Romkes and Schweizer, 2015). However, in our model the primary hip flexors, i.e., the iliopsoas muscles, were not included. van der Krogt et al. (2012) showed that an increased activation of the rectus femoris may be developed to compensate for a weakness of the primary hip flexors. In our case, the absence of the iliopsoas muscles (equivalent to a complete reduction of strength in these muscles) may have induced the increased rectus femoris excitations observed during the simulations compared to the experimental measurements. Because the rectus femoris is a bi-articular muscle (i.e., hip flexor, knee extensor), an increased excitation of this muscle used to assist in hip flexion would simultaneously act to reduce knee flexion, which would have require compensation from knee flexors in order to maintain experimental kinematics (van der Krogt et al., 2012). This could explain the non-physiological activity of the triceps surae (ankle plantarflexors) during pre-swing and early swing observed in our simulations, and the increased tibialis anterior activity (ankle dorsiflexor) used to balance ankle flexion due to the increased activity of ankle plantarflexors. All these observations support the need for a more comprehensive, full body musculoskeletal model, as already discussed in section Limitations.

Kinematics and Ground Reaction Forces

Unlike inverse dynamics-based optimization approaches (i.e., static optimization), where measured kinematics and calculated joint torques are the input constraints to the optimization problem, in a forward dynamics approach it is essential to assess the accuracy of reproduced kinematics. This is generally assessed by tracking experimental kinematics (Erdemir et al., 2007; Chèze et al., 2015). In the present study, marker trajectories were tracked rather than joint kinematics, as proposed by Bélaise et al. (2018a,b). This tracking was able to produce accurate marker trajectories as highlighted in Bélaise et al. (2018a) with a tracking residual of 0.31 ± 0.32 cm during elbow flexion. This approach, used in the present study gave accurate kinematics with a maximum RMSE obtained at the ankle joint < 5°, a threshold recognized as critical for clinical interpretation (McGinley et al., 2009). The main errors appeared at the end of the stance phase in most of the DoFs. This issue may be associated with the non-physiological activity of the triceps surae and the rectus femoris during this phase, as discussed in section Muscle Activity. Our simulations tended to reduce the plantarflexion induced by the emulated equinus gait. By increasing ankle dorsiflexion, the triceps surae moment arm was increased (Delp et al., 1990) and a minimal ankle plantarflexion moment was kept.

Interestingly, GRF were estimated accurately without an advanced foot-ground contact model definition. The use of simple generalized forces applied on the pelvis, designed to compensate for the absence of the opposite lower limb and the upper part of the body in our model, acted as efficient reserve actuators (Modenese et al., 2013) to provide the forces required to track experimental GRF. This approach may thus present a valid means by which to manage external forces and moments in the dynamic equation when the interactions between the upper limb and/or the contralateral lower limb with the ipsilateral lower limb are not known. Otherwise, as already discussed in section Limitations, a full body musculoskeletal model would be recommended.

Clinical Perspectives

The simulation of equinus gait represents an important clinical issue in the context of toe walking, a common gait deviation observed in many pathologies such as cerebral palsy, myopathy, and neuropathy (Armand et al., 2007). Numerical simulations may present a useful tool in this context of identifying potential biomechanical causes of this deviation, such as: pre-tibial muscle weakness, inadequate ankle dorsiflexors activity, ankle plantarflexors contracture, and/or spasticity, excessive voluntary ankle plantarflexion in compensation for quadriceps weakness, knee flexor contracture caused by overactivity of the hamstring, combined spasticity of the hamstring and ankle plantarflexors, and leg length discrepancy (Armand et al., 2007; Perry and Burnfield, 2010). In more general terms, there is a need for a numerical framework allowing for the introduction of pathology (Santos et al., 2017), treatment, or surgical intervention-related (Fox et al., 2009) perturbations in a model, and the analysis of their impact on the structures of the musculoskeletal system during daily living activities. However, before clinical applications, models have to be evaluated and validated (Hicks et al., 2015). It will thus be necessary that we assess the capacity of our approach to produce physiological musculo-tendon forces and joint contact forces. Validation datasets, such as the ones made available by Bergmann et al. (2016) and Fregly et al. (2012), should thus be tested on our numerical framework in the future.

Conclusion

In conclusion, we have improved the recent EMG-marker tracking optimization method to a multiphase cyclic movement with GRF. This numerical framework was successfully tested on a dataset of equinus gait for which our approach was able to estimate lower-limb kinematics, GRF and muscle activity with reasonable accuracy.

Data Availability

The dataset used for this study is available from figshare (https://figshare.com/s/ddef4ef25acc2ea91c04, doi: 10.6084/m9.figshare.7869293). They are all stored in c3d file format (https://www.c3d.org).

Ethics Statement

The participant gave written informed consent prior to his inclusion and the protocol was conformed to the Declaration of Helsinki and approved by the National Research Ethics Committee of Luxembourg (201805/01).

Author Contributions

FM conceived of the presented idea and adapted the theory to the present project. MB, BM, and CB developed the theory. FM and EP performed the computations. MB and BM verified the methods. All authors discussed the results and contributed to the final manuscript.

Funding

This project was supported by the Luxembourg National Research Fund (INTER/MOBILITY/17/11756737) and by the NSERC Discovery grant (RGPIN-2014-03912).

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.

Acknowledgments

The authors gratefully acknowledge the Interdisciplinary Center for Scientific Computing of Heidelberg University (Germany) for the opportunity to work with MUSCOD-II.

References

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

Armand, S., Moissenet, F., de Coulon, G., and Bonnefoy-Mazure, A. (2017). Identifying and understanding gait deviations: critical review and perspectives. Mov. Sport Sci. Sci. Mot. 2017, 77–88. doi: 10.1051/sm/2017016

CrossRef Full Text | Google Scholar

Armand, S., Watelain, E., Roux, E., Mercier, M., and Lepoutre, F.-X. (2007). Linking clinical measurements and kinematic gait patterns of toe-walking using fuzzy decision trees. Gait Posture 25, 475–484. doi: 10.1016/j.gaitpost.2006.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassani, T., Stucovitz, E., Qian, Z., Briguglio, M., and Galbusera, F. (2017). Validation of the AnyBody full body musculoskeletal model in computing lumbar spine loads at L4L5 level. J. Biomech. 58, 89–96. doi: 10.1016/j.jbiomech.2017.04.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Bélaise, C., Dal Maso, F., Michaud, B., Mombaur, K., and Begon, M. (2018a). An EMG-marker tracking optimisation method for estimating muscle forces. Multibody Syst. Dyn. 42, 119–143. doi: 10.1007/s11044-017-9587-2

CrossRef Full Text | Google Scholar

Bélaise, C., Michaud, B., Dal Maso, F., Mombaur, K., and Begon, M. (2018b). Which data should be tracked in forward-dynamic optimisation to best predict muscle forces in a pathological co-contraction case? J. Biomech. 68, 99–106. doi: 10.1016/j.jbiomech.2017.12.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Bergmann, G., Bender, A., Dymke, J., Duda, G., and Damm, P. (2016). Standardized loads acting in hip implants. PLoS ONE 11:e0155612. doi: 10.1371/journal.pone.0155612

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonnefoy-Mazure, A., Sagawa, Y., Pomero, V., Lascombes, P., De Coulon, G., and Armand, S. (2016). Are clinical parameters sufficient to model gait patterns in patients with cerebral palsy using a multilinear approach? Comput. Methods Biomech. Biomed. Engin. 19, 800–806. doi: 10.1080/10255842.2015.1064112

PubMed Abstract | CrossRef Full Text | Google Scholar

Brookham, R. L., Middlebrook, E. E., Grewal, T., and Dickerson, C. R. (2011). The utility of an empirically derived co-activation ratio for muscle force prediction through optimization. J. Biomech. 44, 1582–1587. doi: 10.1016/j.jbiomech.2011.02.077

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

Chantraine, F., Filipetti, P., Schreiber, C., Remacle, A., Kolanowski, E., and Moissenet, F. (2016). Proposition of a classification of adult patients with hemiparesis in chronic phase. PLoS ONE 11:e0156726. doi: 10.1371/journal.pone.0156726

PubMed Abstract | CrossRef Full Text | Google Scholar

Chèze, L., Moissenet, F., and Dumas, R. (2015). State of the art and current limits of musculo-skeletal models for clinical applications. Mov. Sport Sci. Sci. Mot. 90, 7–17. doi: 10.1051/sm/2012026

CrossRef Full Text

Chiou, I. I., Burnett, C. N., Chiou, I. I., and Burnett, C. N. (1985). Values of activities of daily living. A survey of stroke patients and their home therapists. Phys. Ther. 65, 901–906. doi: 10.1093/ptj/65.6.901

PubMed Abstract | CrossRef Full Text | Google Scholar

Cuomo, A. V., Gamradt, S. C., Kim, C. O., Pirpiris, M., Gates, P. E., McCarthy, J. J., et al. (2007). Health-related quality of life outcomes improve after multilevel surgery in ambulatory children with cerebral palsy. J. Pediatr. Orthop. 27, 653–657. doi: 10.1097/BPO.0b013e3180dca147

PubMed Abstract | CrossRef Full Text | Google Scholar

Davids, J. R., Ounpuu, S., DeLuca, P., a, and Davis, R. B. (2004). Optimization of walking ability of children with cerebral palsy. Instr. Course Lect. 53, 511–522. doi: 10.2106/00004623-200311000-00028

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

Delp, S. L., Loan, J. P., Hoy, M. G., Zajac, F. E., Topp, E. L., and Rosen, J. M. (1990). An interactive graphics-based model of the lower extremity to study orthopedic surgical procedures. IEEE Trans. Biomed. Eng. 37, 757–767. doi: 10.1109/10.102791

CrossRef Full Text | Google Scholar

Erdemir, A., McLean, S., Herzog, W., and Bogert, A. J. (2007). Model-based estimation of muscle forces exerted during movements. Clin. Biomech. Bristol. Avon. 22, 131–154. doi: 10.1016/j.clinbiomech.2006.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Felis, M. L. (2017). RBDL: an efficient rigid-body dynamics library using recursive algorithms. Auton. Robots 41, 495–511. doi: 10.1007/s10514-016-9574-0

CrossRef Full Text | Google Scholar

Fluit, R., Andersen, M. S., Kolk, S., Verdonschot, N., and Koopman, H. F. J. M. (2014a). Prediction of ground reaction forces and moments during various activities of daily living. J. Biomech. 47, 2321–2329. doi: 10.1016/j.jbiomech.2014.04.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Fluit, R., Andersen, M. S., Verdonschot, N., and Koopman, H. F. J. M. (2014b). Optimal inverse dynamic simulation of human gait. Gait Posture 39S, S1–S141. doi: 10.1016/j.gaitpost.2014.04.060

CrossRef Full Text | Google Scholar

Fohanno, V., Begon, M., Lacouture, P., and Colloud, F. (2014). Estimating joint kinematics of a whole body chain model with closed-loop constraints. Multibody Syst. Dyn. 31, 433–449. doi: 10.1007/s11044-013-9366-7

CrossRef Full Text | Google Scholar

Fox, M. D., Reinbolt, J. A., Õunpuu, S., and Delp, S. L. (2009). Mechanisms of improved knee flexion after rectus femoris transfer surgery. J. Biomech. 42, 614–619. doi: 10.1016/j.jbiomech.2008.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Fregly, B. J., Besier, T. F., Lloyd, D. G., Delp, S. L., Banks, S. A., Pandy, M. G., et al. (2012). Grand challenge competition to predict in vivo knee loads. J. Orthop. Res. Off. Publ. Orthop. Res. Soc. 30, 503–513. doi: 10.1002/jor.22023

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaudet, S., Tremblay, J., and Dal Maso, F. (2018). Evolution of muscular fatigue in periscapular and rotator cuff muscles during isokinetic shoulder rotations. J. Sports Sci. 36, 2121–2128. doi: 10.1080/02640414.2018.1440513

PubMed Abstract | CrossRef Full Text | Google Scholar

Giroux, M., Moissenet, F., and Dumas, R. (2013). EMG-based validation of musculo-skeletal models for gait analysis. Comput. Methods Biomech. Biomed. Engin. 16, 152–154. doi: 10.1080/10255842.2013.815878

PubMed Abstract | CrossRef Full Text | Google Scholar

Gough, M., and Shortland, A. P. (2008). Can clinical gait analysis guide the management of ambulant children with bilateral spastic cerebral palsy? J. Pediatr. Orthop. 28, 879–883. doi: 10.1097/BPO.0b013e31818e197c

PubMed Abstract | CrossRef Full Text | Google Scholar

Hermens, H. J., Freriks, B., Disselhorst-Klug, C., and Rau, G. (2000). Development of recommendations for SEMG sensors and sensor placement procedures. J. Electromyogr. Kinesiol. 10, 361–374. doi: 10.1016/S1050-6411(00)00027-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Hicks, J. L., Uchida, T. K., Seth, A., Rajagopal, A., and Delp, S. L. (2015). Is my model good enough? Best practices for verification and validation of musculoskeletal models and simulations of movement. J. Biomech. Eng. 137:020905. doi: 10.1115/1.4029304

PubMed Abstract | CrossRef Full Text | Google Scholar

Jansen, K., De Groote, F., Aerts, W., De Schutter, J., Duysens, J., and Jonkers, I. (2014). Altering length and velocity feedback during a neuro-musculoskeletal simulation of normal gait contributes to hemiparetic gait characteristics. J. Neuroeng. Rehabil. 11:78. doi: 10.1186/1743-0003-11-78

PubMed Abstract | CrossRef Full Text | Google Scholar

Jung, Y., Jung, M., Ryu, J., Yoon, S., Park, S.-K., and Koo, S. (2016). Dynamically adjustable foot-ground contact model to estimate ground reaction force during walking and running. Gait Posture 45, 62–68. doi: 10.1016/j.gaitpost.2016.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Leineweber, D. B., Bauer, I., Bock, H. G., and Schlöder, J. P. (2003). An efficient multiple shooting based reduced SQP strategy for large-scale dynamic process optimization. Part 1: theoretical aspects. Comput. Chem. Eng. 27, 157–166. doi: 10.1016/S0098-1354(02)00158-8

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

McGinley, J. L., Baker, R., Wolfe, R., and Morris, M. E. (2009). The reliability of three-dimensional kinematic gait measurements: A systematic review. Gait Posture 29, 360–369. doi: 10.1016/j.gaitpost.2008.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Modenese, L., Gopalakrishnan, A., and Phillips, A. T. M. (2013). Application of a falsification strategy to a musculoskeletal model of the lower limb and accuracy of the predicted hip contact force vector. J. Biomech. 46, 1193–1200. doi: 10.1016/j.jbiomech.2012.11.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Moissenet, F., Chèze, L., and Dumas, R. (2014). Introduction of a set of EMG-based muscular activations in a multi-objective optimisation when solving the muscular redundancy problem during gait. Comput. Methods Biomech. Biomed. Engin. 17, 132–133. doi: 10.1080/10255842.2014.931542

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, Y., Zhang, Z., Gao, Y., Chen, Z., Xin, H., Zhang, Q., et al. (2018). Concurrent prediction of ground reaction forces and moments and tibiofemoral contact forces during walking using musculoskeletal modeling. Med. Eng. Phys. 52, 31–40. doi: 10.1016/j.medengphy.2017.11.008

CrossRef Full Text | Google Scholar

Perry, D. J., and Burnfield, D. J. (2010). Gait Analysis: Normal and Pathological Function. Thorofare, NJ: SLACK Inc.

Google Scholar

Pizzolato, C., Lloyd, D. G., Sartori, M., Ceseracciu, E., Besier, T. F., Fregly, B. J., et al. (2015). CEINMS: A toolbox to investigate the influence of different neural control solutions on the prediction of muscle excitation and joint moments during dynamic motor tasks. J. Biomech. 48, 3929–3936. doi: 10.1016/j.jbiomech.2015.09.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajagopal, A., Dembia, C. L., DeMers, M. S., Delp, D. D., Hicks, J. L., and Delp, S. L. (2016). Full-body musculoskeletal model for muscle-driven simulation of human gait. IEEE Trans. Biomed. Eng. 63, 2068–2079. doi: 10.1109/TBME.2016.2586891

PubMed Abstract | CrossRef Full Text | Google Scholar

Romkes, J., and Schweizer, K. (2015). Immediate effects of unilateral restricted ankle motion on gait kinematics in healthy subjects. Gait Posture 41, 835–840. doi: 10.1016/j.gaitpost.2015.02.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Santos, G. F., Gomes, A. A., Sacco, I. C. N., and Ackermann, M. (2017). Predictive simulation of diabetic gait: Individual contribution of ankle stiffness and muscle weakening. Gait Posture 58, 208–213. doi: 10.1016/j.gaitpost.2017.07.124

PubMed Abstract | CrossRef Full Text | Google Scholar

Sartori, M., Farina, D., and Lloyd, D. G. (2014). Hybrid neuromusculoskeletal modeling to best track joint moments using a balance between muscle excitations derived from electromyograms and optimization. J. Biomech. 47, 3613–3621. doi: 10.1016/j.jbiomech.2014.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Sartori, M., Reggiani, M., van den Bogert, A. J., and Lloyd, D. G. (2012). 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

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, Q., Bassett, D. N., Manal, K., and Buchanan, T. S. (2009). An EMG-driven model to estimate muscle forces and joint moments in stroke patients. Comput. Biol. Med. 39, 1083–1088. doi: 10.1016/j.compbiomed.2009.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Son, J., Hwang, S., and Kim, Y. (2012). A hybrid static optimisation method to estimate muscle forces during muscle co-activation. Comput. Methods Biomech. Biomed. Engin. 15, 249–254. doi: 10.1080/10255842.2010.522187

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

Thompson, J. A., Chaudhari, A. M. W., Schmitt, L. C., Best, T. M., and Siston, R. A. (2013). Gluteus maximus and soleus compensate for simulated quadriceps atrophy and activation failure during walking. J. Biomech. 46, 2165–2172. doi: 10.1016/j.jbiomech.2013.06.033

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Krogt, M. M., Delp, S. L., and Schwartz, M. H. (2012). How robust is human gait to muscle weakness? Gait Posture 36, 113–119. doi: 10.1016/j.gaitpost.2012.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Krogt, M. M., Doorenbosch, C. A. M., Becher, J. G., and Harlaar, J. (2009). Walking speed modifies spasticity effects in gastrocnemius and soleus in cerebral palsy gait. Clin. Biomech. 24, 422–428. doi: 10.1016/j.clinbiomech.2009.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

van Schie, C. H. M. (2008). Neuropathy: mobility and quality of life. Diabetes Metab. Res. Rev. 24(Suppl 1), S45–51. doi: 10.1002/dmrr.856

PubMed Abstract | CrossRef Full Text | Google Scholar

van Sint Jan, S. (2007). Color Atlas of Skeletal Landmark Definitions. Edinburgh; New York, NY: Churchill Livingstone/Elsevier. 1623.

Google Scholar

Wren, T. A. L., Gorton, G. E., Õunpuu, S., and Tucker, C. A. (2011). Efficacy of clinical gait analysis: a systematic review. Gait Posture 34, 149–153. doi: 10.1016/j.gaitpost.2011.03.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Zajac, F. E. (1989). Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Crit. Rev. Biomed. Eng. 17, 359–411.

PubMed Abstract | Google Scholar