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

In the context of neuro-orthopedic pathologies affecting walking and thus patients' quality of life, understanding the mechanisms of gait deviations and identifying the causal motor impairments is of primary importance. Beside other approaches, neuromusculoskeletal simulations may be used to provide insight into this matter. To the best of our knowledge, no computational framework exists in the literature that allows for predictive simulations featuring muscle co-contractions, and the introduction of various types of perturbations during both healthy and pathological gait types. The aim of this preliminary study was to adapt a recently proposed EMG-marker tracking optimization process to a lower limb musculoskeletal model during equinus gait, a multiphase problem with contact forces. The resulting optimization method tracking EMG, ground reactions forces, and marker trajectories allowed an accurate reproduction of joint kinematics (average error of 5.4 ± 3.3 mm for pelvis translations, and 1.9 ± 1.3° for pelvis rotation and joint angles) and ensured good temporal agreement in muscle activity (the concordance between estimated and measured excitations was 76.8 ± 5.3 %) in a relatively fast process (3.88 ± 1.04 h). We have also highlighted that the tracking of ground reaction forces was possible and accurate (average error of 17.3 ± 5.5 N), even without the use of a complex foot-ground contact model.


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 cocontractions. 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 dynamicsbased 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 cocontractions. 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 cocontractions, and that allows for the introduction of various kind of perturbations on the model (e.g., therapy-related, surgeryrelated, 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.

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).
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).  (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. 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 accelerationsq of the rigid multibody system were computed using a forward dynamics approach for given generalized joint positions q, joint velocitiesq, and generalized forces τ: 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 τ 1 = [τ 11 τ 12 τ 13 ] T driving the three pelvis DoFs, and τ 2 = ∂L mt ∂q F mt corresponding to the net joint torques due to the musculo-tendon Palpation details used to place experimental reflective cutaneous markers (as well as virtual markers by virtual palpation) and the related segment are also mentioned.
forces F mt , where ∂L mt ∂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 Hilltype muscle model with a generic force-length-velocity relation f (Zajac, 1989): where F 0 mt 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 1 and excitations e 1 ), corresponding to the muscles for which electromyographic (EMG) records were available, and F mt2 (with related activations a 2 and excitations e 2 ) where EMG measurements were unavailable. Muscle activation dynamics was implemented as a set of first-order differential equations (Buchanan et al., 2004): 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 (τ 1 ), 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 2 that produced the prescribed marker trajectories, muscle neural excitations, and GRF (during stance phase only): , the swing phase.
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.  −300 300 Frontiers in Neurorobotics | www.frontiersin.org

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 neuroorthopedic 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,q,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 2 ) 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 R Core TM 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 2 and CC) are reported in Table 4.
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 2 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 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.  Table 2.
Frontiers in Neurorobotics | www.frontiersin.org  Table 3. 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 2 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 2 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 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). 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 preswing 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.

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.