Validation of forward simulations to predict the effects of bilateral plantarflexor weakness on gait

Background: Bilateral plantarflexor muscle weakness is a common impairment in many neuromuscular diseases. However, the way in which severity of plantarflexor weakness affects gait in terms of walking energy cost and speed is not fully understood. Predictive simulations are an attractive alternative to human experiments as simulations allow systematic alterations in muscle weakness. However, simulations of pathological gait have not yet been validated against experimental data, limiting their applicability. Research question: Our first aim was to validate a predictive simulation framework for walking with bilateral plantarflexor weakness by comparing predicted gait against experimental gait data of patients with bilateral plantarflexor weakness. Secondly, we aimed to evaluate how incremental levels of bilateral plantarflexor weakness affect gait. Methods: We used a planar musculoskeletal model with 9 degrees of freedom and 9 Hill-type muscle-tendon units per leg. A state-dependent reflex-based controller optimized for a function combining energy cost, muscle activation squared and head acceleration was used to simulate gait. For validation, strength of the plantarflexors was reduced by 80 % and simulated gait compared with experimental data of 16 subjects with bilateral plan- tarflexor weakness. Subsequently, strength of the plantarflexors was reduced stepwise to evaluate its effect on gait kinematics and kinetics, walking energy cost and speed. Results: Simulations with 80 % weakness matched well with experimental hip and ankle kinematics and kinetics (R > 0.64), but less for knee kinetics (R < 0.55). With incremental strength reduction, especially beyond a reduction of 60 %, the maximal ankle moment and power decreased. Walking energy cost and speed showed a strong quadratic relation (R 2 > 0.82) with plantarflexor strength. Significance: Our simulation framework predicted most gait changes due to bilateral plantarflexor weakness, and indicates that pathological gait features emerge especially when bilateral plantarflexor weakness exceeds 60 %. Our framework may support future research into the effect of pathologies or assistive devices on gait.


Introduction
In many neuromuscular diseases, like Charcot-Marie-Tooth or polio, gait can be impaired due to bilateral plantarflexor weakness [1,2].
Typically, such weakness leads to a gait pattern characterized by excessive ankle dorsiflexion, reduced ankle moment, persistent knee flexion and reduced internal knee flexion moment during stance, along with reduced ankle push-off power [1][2][3]. These gait deviations reduce walking speed [4] and increase walking energy cost [4,5], thereby limiting individuals in performing daily activities [6,7].
The degree of gait deviations vary between individuals as severity of plantarflexor weakness ranges from mild to completely paralyzed [1,8]. Mild ankle plantarflexor weakness already increases walking energy cost [7] and imposes limitations in performing higher intensity activities in individuals with Charcot-Marie-Tooth [6]. However, this is not necessarily reflected by a change in gait kinematics or kinetics due to compensatory muscle activations [9,10]. In more severely affected individuals, necessary compensatory muscle activations exaggerate until walking is only possible with severe gait deviations requiring substantial effort [4,9].
The effects of incremental plantarflexor weakness on gait compensations and deterioration of walking energy cost and speed have not been demonstrated, while knowledge about the point where progressive weakness causes severe walking difficulties might be useful for understanding pathological gait. Establishing this effect is challenging as systematic alterations in muscle weakness cannot be imposed in human experiments. Computer simulations are not limited in this regard, and have been used to gain knowledge about the required muscle strength for normal gait [9] as well as for crouch gait [11] by tracking experimentally recorded gait patterns. However, due to the required input of gait data, these tracking-simulations are unable to predict the effect of incremental weakness on gait kinematics and kinetics.
Forward dynamic simulations that can predict changes in kinematics and kinetics have the potential to provide insights into how impairments affect gait, as previously used to explain gait of elderly [12]. Furthermore, Ong et al. showed that predictive simulations with bilateral plantarflexor weakness demonstrate gait adaptations, such as a heel-walking gait pattern [13]. However, they only simulated extreme cases of weakness (less than 25 % remaining force) while simultaneously reducing passive muscle forces, which in patients are not necessarily linked [14]. Additionally, the simulations were validated against unimpaired gait, but not against specific pathological gait data [13,15]. Consequently, whether forward simulations predict compensation strategies as used by humans is uncertain.
Therefore, our first aim was to simulate healthy gait and gait with 80 % bilateral plantarflexor weakness to validate the predictive simulation framework by comparing predicted gait against a large dataset of both healthy subjects and individuals with bilateral plantarflexor weakness (median weakness 78 %) consisting of joint kinematics, kinetics, ground reaction forces, walking energy cost and speed. Secondly, we simulated gait with a stepwise increase in bilateral plantarflexor weakness to evaluate how gait outcomes are affected by incremental weakness.
3D-gait data of barefoot walking at comfortable speed were collected using an 8-camera 100 Hz Vicon MX 1.3 system (VICON, Oxford, UK) and four force plates (1000 Hz, OR6-7, AMTI, Watertown, USA). Markers were placed according to the Plug-In-Gait model [17]. Marker data were processed in OpenSim, which is an open source neuro-musculoskeletal modelling program [18]. After scaling the generic OpenSim model using a static trial, joint angles and moments were calculated using the inverse kinematic and inverse dynamic toolboxes [19]. Joint powers were calculated from joint angles and moments using custom-made scripts.
Walking speed (m/s) and walking energy cost (J/kg/m) were assessed during a 6-minute walk test at comfortable speed on a 35-meter oval track. During the test, oxygen consumption (VO 2 ) and carbon dioxide production (VCO 2 ) were measured (Cosmed K4B 2 , Rome, Italy). A steady-state period of at least 60 s for which VO 2 and speed were constant was selected. The mean gross energy consumption during steadystate (in J/kg/s) was calculated according to Garby & Astrup [20] and divided by steady-state walking speed to obtain walking energy cost (in J/kg/m).
To assess the severity of plantarflexor weakness of individuals with bilateral plantarflexor weakness, maximal isometric plantarflexor strength was measured with a fixed dynamometer (Biodex type 3, Corp., Shirley, NY, USA). The highest value (in Nm) of three trials was divided by a gender-specific reference of our gait laboratory (age: 60 ± 6.3) to determine percentage weakness [21]. The median strength reduction  was 78 % (first quartile:59 %, third quartile: 87 %).

Simulation framework
To reduce the number of optimization parameters and redundancy of the neuromuscular system, we simplified the OpenSim gait2392 model, originally consisting of 23 degrees of freedom and 92 muscles. We created a planar model with seven segments (trunk-pelvis, and bilateral thigh, shank and foot), and nine degrees of freedom; three around the trunk-pelvis segment and one for hip, knee and ankle flexion [18]. This planar model was justified as most prominent compensations of plantarflexor weakness occur in the sagittal plane [3]. Nine hill-type muscles (Millard-Equilibrium muscle model [22]) per leg were modelled, including the Tibialis Anterior, Soleus, Gastrocnemius medialis, Vastus intermedius, Rectus Femoris, Semitendinosus, Biceps Femoris short head, Gluteus Maximus and Iliopsoas. Muscle path, optimal fiber length, pennation angle and tendon slack length were set as in the Gait2392 model. Peak isometric forces were summed for all muscles performing the same movement similar to Ong et al. [13]. Knee ligaments were modelled as a rotational spring (2 Nm/degree) and damper (0.2 Nm/degree/s)), activated when the knee moved beyond 120 degrees of flexion or extended beyond 5 degrees of flexion.
Ground contact was modelled using two viscoelastic Hunt-Crossley contact spheres on each foot [23]. The spheres' location and radius were set to match the rounding of the heel and metatarsal joints. The spheres' static, dynamic and viscous friction were set to 1, while stiffness and dampening were optimized by imposing healthy lower leg kinematics and minimizing the error between sphere forces and experimental ground reaction forces. See Appendix A for an overview of all model parameters.
Muscles were activated by a gait state-dependent reflex-based controller, based on Geyer & Herr [24] ( Table 1). The controller consisted of constant signals and muscle force and length feedback to generate muscle excitations. For the hamstrings, iliopsoas and gluteus maximus a proportional-derivative feedback loop was implemented during stance to stabilize the trunk. The original controller was adjusted on a few aspects. Reflex gains could differ between five gait states, similar to previous research [25,26], as reflexes are known to be modulated during gait [27]. This allowed us to add a force reflex for the soleus and gastrocnemius during early stance, which led to a more gradual increase in activation and better match with experimental EMG. Additionally, similar to Ong et al. [13], we added the biceps femoris short head and rectus femoris to the Geyer & Herr controller to have both uni-articular muscles around the knee and bi-articular muscles surrounding the knee and hip, respectively. Feedback gains within each state, transition thresholds between the phases and the initial joint angles were optimized using SCONE, which is open-source software to run predictive neuromuscular simulations [28].
In total, 103 parameters were optimized by minimizing a cost function consisting of walking energy cost (JCost), muscle activation squared (JMuscle) and head acceleration (JHeadAcc) per meter, while walking without falling down (JFall) and avoiding extreme joint angles at the ankle and knee (JAng). JCost penalized metabolic cost of walking calculated according to Uchida et al. [29]. JMuscle and JHeadAcc were calculated by dividing the sum of the muscle activation squared and head acceleration by the horizontal distance travelled. JFall penalized a center-of-mass lower than 0.85 of its starting value, representing the beginning of a fall. JAng penalized ankle angles exceeding 60 degrees of dorsiflexion or plantarflexion and when knee ligaments provided more than 5 Nm during the simulation time. Weightings of the different components were manually tuned to best match unimpaired gait and set as follows for all simulations: 1E8 * JFall + 0.1 * JAng + 0.15 * JCost + 15 * JMuscle + 0.1JHeadAcc. In the final simulation outcomes, JFall and JAng did not contribute to the score of the cost function.   Gait was simulated for 10 s without constraining walking speed, using SCONE. Control parameters were optimized using the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [30]. Each generation consisted of 17 simulations. The optimization was terminated when the average reduction of the cost function score in the last 500 generations was smaller than 0.001 %. For each condition three optimizations with different random seeds and a standard deviation (SD) of 0.05 for all parameters were performed. For better comparison ground reaction forces were filtered using a low-pass 6 Hz Butterworth filter similar to experimental ground reaction forces.

Validation
For validation, we compared simulations of unimpaired gait and gait with bilateral ankle plantarflexor weakness with experimental data. To model plantarflexor weakness, we reduced maximal force of the soleus and gastrocnemius muscle by 80 %, so 20 % of the original force remained, in order to match median weakness (78 %) of our subjects. In diseases maximal force and passive force are not related as in healthy subjects [14] and as passive stiffness for the included subjects was not known, we chose to adapt passive fiber and tendon force-length curves such that they matched those of the unimpaired model.
We quantified the shape-agreement between the simulated and experimental ground reaction forces, joint angles, moments and powers for the unimpaired and bilateral plantarflexor weakness group using time-normalized cross-correlations (R). A correlation coefficient <0.3 was considered poor, between 0.3− 0.5 fair, between 0.5− 0.7 moderate and above 0.7 strong [31]. Differences between simulated and experimental gait trajectories were quantified by the root mean square error (RMSE) normalised to the standard deviation (SD) of the experimental data. The higher SD for patient-data must be considered when interpreting RMSE-values. We tested whether simulated walking speed and walking energy cost differed significantly from experimental data using a one-sample t-test.

Effect of incremental plantarflexor weakness on gait
To study the effect of incremental bilateral plantarflexor weakness on gait, we reduced maximal isometric force of the soleus and gastrocnemius muscle by 20 %, 40 %, 60 %, 70 %, 80 % and 90 %. Consistent changes in joint kinematics and kinetics were qualitatively described. To assess the relation between walking speed and walking energy cost with    incremental bilateral plantarflexor weakness, we tested both a linear and quadratic fit.

Model validation
The unimpaired simulations showed strong cross-correlations for all joint kinetics, kinematics and ground reaction forces (R > 0.82), except hip power (moderate; R = 0.68) ( Fig. 1 & Table 2). All joint angles had an RMSE within 2 SD, while for the kinetics the ankle moment (2.56 SD), knee power (2.17 SD) and hip power (2.59 SD) had an RMSE above 2.
High RMSE values for both vertical (2.29 SD) and horizontal (4.09 SD) ground reaction forces were found.
Walking energy cost of the unimpaired simulation (4.21 J/kg/m) was significantly higher compared to reference data of healthy individuals (3.66 ± 0.47 J/kg/m, p < 0.001). Walking energy cost of the simulation with 80 % plantarflexor weakness (5.21 J/kg/m) did not differ significantly from experimental data of individuals with bilateral plantarflexor weakness (5.04±1.06 J/kg/m, p = 0.541).

Effect of incremental plantarflexor weakness on gait
With increasing ankle plantarflexor weakness, several consistent changes in the gait pattern were found. In stance, maximal ankle dorsiflexion angle was reached later in the gait cycle (from 46 % for the unimpaired simulations to 58 % in case of 90 % weakness), while the knee demonstrated less flexion in early stance and the hip moved less towards extension (Fig. 2). Regarding gait kinetics, maximal ankle moment (from to 1.44 to 0.21 Nm/kg) and power (from 1.57 to 0.07 W/ kg) decreased, especially when weakness exceeded 60 % (at 60 % weakness: max. ankle moment = 1.33 Nm/kg, max. ankle power = 1.19 W/kg). At the knee level, normal knee moment curves were found until 80 % weakness, from where onwards the model walked with a neutral knee moment. The first peak of the vertical ground reaction force showed no clear trend (varied between 0.99 -1.05 body weight), while the second peak decreased (from 1.09 to 0.97 body weight) with incremental weakness.

Discussion
Our predictive simulation framework can generate human-like walking and predict most gait changes due to bilateral plantarflexor weakness, as the majority of simulated kinematics and kinetics showed a good agreement with experimental gait data. In addition, our simulation framework indicated that in case bilateral plantarflexor weakness is 60 % or more large changes in gait emerge, leading to rapid deterioration of walking speed and increment of walking energy cost.
Despite the simplifications of a 2D model, unimpaired walking simulations demonstrated good agreement with experimental data for all sagittal plane gait kinetics and kinematics. Our simulations showed knee flexion in early stance improving the cross-correlation compared to previous simulations [13]. This improvement is likely explained by incorporation of muscle activation squared in our cost function, as it allows for activation of the large vasti muscles in stance and was previously shown to result in a better knee flexion during stance [12,32]. Song et al. were able to simulate knee flexion without including a muscle activation term, although they incorporated an additional control layer defining foot placement and heel-strike potentially explaining the knee flexion [26]. However, introducing such an additional layer limits possible foot-placement compensations, such as forefoot landing. Additionally, the muscle activation term was essential to avoid fatigue from unrealistically high muscle activations [15], although at the cost of the lower walking speed (1.05 m/s compared to 1.25 in previous work) [13]. In simulations of elderly gait, it has been shown that especially including a fatigue term explains reductions in speed while an energy cost term does not [12]. This may also hold in people with plantarflexor weakness, where fatigue of the plantarflexors could be an additional factor contributing to a slower walking speed. The slower simulated speed might have reduced the match with experimental gait kinetics and kinematics collected at higher speeds. This could explain the reduced knee flexion in swing and relatively higher ground reaction force during mid-stance, causing the high RMSE for ground reaction force variables. Additionally, the lower speed might explain why simulated energy cost was higher compared to experimental data.
Together with the lower walking speed, a lower ankle moment in mid-stance and reduced plantar flexion during push-off were the most prominent discrepancies between unimpaired simulations and reference experimental data. Potentially, both are caused by a late onset of activating the gastrocnemius and soleus muscles (see appendix B). This late onset could be a limitation of our controller, which is primarily based on reflexes and only uses a force-based reflex for the gastrocnemius and soleus. The importance of such reflexes during healthy gait is well established [33], but how different reflex-types within the various muscles contribute to impaired gait remains unknown. Consequently, our controller may be too compromised to establish a timely onset of the plantarflexor muscles.
Using the same framework, most effects of bilateral ankle plantarflexor weakness on gait could be predicted. For the ankle joint, a larger maximal dorsiflexion angle at a later time-point in the gait cycle and lower maximal ankle moment and power were predicted [3,4]. Previous studies simulating bilateral plantarflexor weakness found gait patterns where the ankle remained in an excessive dorsiflexion throughout the gait cycle [13,15], which is not typical for people with plantarflexor weakness [3,4]. The main difference with our study is that those studies reduced passive muscle forces together with maximal force, as these are coupled in the default OpenSim model, thereby also limiting the energy recoil of the Achilles tendon in pre-swing. This indicates that to accurately predict pathological gait, not only the maximal force but also passive forces should be assessed and individualized as they are not necessarily matched in case of pathologies [14].
Effects of plantarflexor weakness on the knee were less accurately predicted compared to the ankle and hip. The simulation framework predicted a reduced knee flexion in early stance combined with a neutral moment throughout stance, while our experimental data show a persistent knee flexion and internal knee extension moment throughout stance. Although the gait pattern of our experimental data is most common in people with plantarflexor weakness [3], gait with extended knees has been reported in people with unilateral calf muscle weakness and slow walking speed [34]. Indeed, in our subject group the slow walkers tended to have less knee flexion in early stance. To test the hypothesis that walking speed affects knee joint kinetics and kinematics in this population, we performed an additional simulation with the 80 % weakness model at a fixed speed of 0.85 m/s (speed of experimental data). Indeed, simulations at a higher speed demonstrated more knee flexion (see Appendix C). Additionally, besides muscle weakness, in patients the neuromuscular control can be affected, such as a reduced maximal neural drive or neural conduction speed [21]. These changes were not modelled but potentially affect gait compensations and might explain the lower match regarding knee kinematics and kinetics.
When simulating incremental ankle plantarflexor weakness, we found that gait deteriorated markedly when muscle force was reduced beyond 60 %, i.e. less than 40 % remaining force. Previous simulations were only able to track a healthy gait pattern with more than 60 % strength reduction by having unrealistically high muscle activations [9]. Similarly, our framework predicted that beyond 60 % weakness walking energy cost and speed worsened dramatically, coinciding with marked reductions in ankle power and shift towards a neutral knee moment throughout stance. Although the 60 % reduction depends on what is considered "normal" strength, it can be concluded, based on these simulations that there is a threshold after which gait deteriorates noticeably. Potentially this may coincide with the timing people seek assistance and start using orthoses [35].

Future research
Despite the limitations of our framework, its predictive ability may support future research on how impairments and assistive devices affect gait [36,37]. Additionally, research should focus on predicting gait in case of unilateral weakness to extend generalizability of the framework. Lastly, research should further develop and personalise the models, controller and cost function in order to increase its predictive value and allow for valid, personalized simulations for individual subjects, which is an essential next step towards further clinical application to evaluate the effect of pathologies or assistive devices on individual patient gait.

Conclusion
Our simulation framework showed an acceptable agreement between simulated and experimental kinetics and kinematics in people with bilateral ankle plantarflexor weakness. Based on the simulations, including its assumptions, gait kinematics and kinetics change substantially after a threshold of around 60 % weakness or more, thereby greatly reducing walking speed, and increasing walking energy cost. In the future and with necessary improvements, our framework could be used to predict individual pathological gait and help understand how assistive devices may improve gait.

Data availability statement
Data are available upon request.

Declaration of Competing Interest
None to declare.

Acknowledgement
This study was supported by Amsterdam Movement Sciences, under the Innovation Call 2018, Tenure Development grant to MK.

Appendix A
The foot-ground contact model See Table A1 In order to set the parameters of the Hunt & Crossley foot contact spheres, we optimized them to best match with both experimental foot kinematics and GRFs. For this purpose, we used experimental data of a subset of six subjects. For these subjects we performed inverse kinematics with the same musculoskeletal model except that instead of using the pelvis we used the calcaneus as free body with the ground. Foot kinematics were averaged over the subjects. To set the contact-sphere parameters we imposed the average foot kinematics, which were the joint angles and velocities expressed as Bsplines, to the musculoskeletal model. First, the spheres location and radius were set to match the rounding of the heel and metatarsal joints. The spheres' static, dynamic and viscous friction were set to 1 to avoid slipping. Subsequently, the stiffness and dissipation parameters of both spheres were optimized by minimizing the root mean square error (RMSE) between the experimental GRF and the GRF calculated by the contact model.
See Table A2 Appendix B Muscle activations of the unimpaired model (black) and plantarflexor model (red). Experimental data are from Bovi et al. [38]. No experimental data of the iliopsoas are available. Exp. = Experimental; Sim. = simulation.