Physics-based predictive simulations to explore the differential effects of motor control and musculoskeletal deficits on gait dysfunction in cerebral palsy: a retrospective case study

Model-based simulations of walking have the theoretical potential to support clinical decision making by predicting the functional outcome of treatments in terms of walking performance. Yet before using such simulations in clinical practice, their ability to identify the main treatment targets in specific patients needs to be demonstrated. In this study, we generated predictive simulations of walking with a medical imaging based neuro-musculoskeletal model of a child with cerebral palsy presenting crouch gait. We explored the influence of altered muscle-tendon properties, reduced neuromuscular control complexity, and spasticity on gait function in terms of joint kinematics, kinetics, muscle activity, and metabolic cost of transport. We modeled altered muscle-tendon properties by personalizing Hill-type muscle-tendon parameters based on data collected during functional movements, simpler neuromuscular control by reducing the number of independent muscle synergies, and spasticity through delayed muscle activity feedback from muscle force and force rate. Our simulations revealed that, in the presence of aberrant musculoskeletal geometries, altered muscle-tendon properties rather than reduced neuromuscular control complexity and spasticity were the primary cause of the crouch gait pattern observed for this child, which is in agreement with the clinical examination. These results suggest that muscle-tendon properties should be the primary target of interventions aiming to restore a more upright gait pattern for this child. This suggestion is in line with the gait analysis following muscle-tendon property and bone deformity corrections. The ability of our simulations to distinguish the contribution of different impairments on walking performance opens the door for identifying targeted treatment strategies with the aim of designing optimized interventions for neuro-musculoskeletal disorders.

that the altered muscle-tendon properties rather than the control impairments alone caused a crouch gait 126 pattern. As an additional analysis, we investigated whether the child's impairments impede a walking 127 pattern similar to TD walking or rather make such a walking pattern less optimal. To this aim, we extended 128 the performance criterion of the predictive simulations with a tracking term that penalized deviations from 129 a TD walking pattern. We found that the musculoskeletal impairments did not prevent an upright walking 130 pattern resembling TD walking but that upright walking was less optimal than walking in crouch.

MATERIAL AND METHODS
The overall process to evaluate the effects of impairments on walking performance through predictive 132 simulations is outlined in Figure 1B. The following sections provide details of this process.

133
Experimental data 134 We collected data from one child with diplegic CP (male; age: 15 years; height: 143 cm; mass: 33.1 135 kg). The data collection was approved by the Ethics Committee at UZ Leuven (Belgium) and written 136 informed consent was obtained from the child's parents. The child was instrumented with retro-reflective 137 skin mounted markers whose 3D trajectories were recorded (100 Hz) using a motion capture system (Vicon,138 Oxford, UK) during overground walking at self-selected speed. Ground reaction forces were recorded 139 (1000 Hz) using force plates (AMTI, Watertown, USA). EMG was recorded (2000 Hz) using a telemetric 140 Zerowire system (Cometa, Milan, Italy) from eight muscles of each leg (rectus femoris, biceps femoris 141 short head, semitendinosus, tibialis anterior, gastrocnemius lateralis, vastus lateralis, soleus, and gluteus 142 medius). EMG from the rectus femoris and vastus lateralis was of poor quality and excluded from the 143 analysis. 144 On the same day as the gait analysis, spasticity of the right medial hamstrings and gastrocnemii was 145 assessed using an instrumented passive spasticity assessment (IPSA; described in detail by Bar-On et al. 146 (2013)). Hamstrings and gastrocnemii were passively stretched by moving knee and ankle, respectively, 147 This is a preprint This MRI-based model as well as experimental data collected during walking and instrumented passive spasticity assessments (IPSA) are inputs to optimization procedures providing personalized estimates of Hill-type muscle-tendon parameters characterizing altered muscle-tendon properties and personalized feedback gains characterizing spasticity. The framework for predictive simulations generates gait patterns by optimizing a cost function, describing a walking-related performance criterion, subject to the muscle and skeleton dynamics of the MRI-based musculoskeletal model. We investigated the effects of impairments on predicted gait patterns (dotted arrows): Q i we evaluated the effect of altered versus unaltered muscle-tendon properties by using personalized versus generic muscle-tendon parameters in the muscle dynamics; Q ii we assessed the influence of reducing the neuromuscular control complexity by imposing a reduced number of muscle synergies; Q iii we explored the impact of spasticity on walking performance. Details on how we modeled these impairments are described in the methods. As an additional analysis, Q iv , we evaluated how well the model was able to reproduce the gait pattern of a typically developing (TD) child by adding a term in the cost function penalizing deviations between predicted gait pattern and measured gait data of a TD child. All these analyses can be combined as well as performed in isolation. Details are provided in section "model-based analyses". one at a time from a predefined position throughout the full range of motion (ROM). The stretches 148 were performed at slow and fast velocities. EMG was collected from four muscles (semitendinosus, 149 gastrocnemius lateralis, rectus femoris, and tibialis anterior) using the same system and electrode placement 150 as used for gait analysis. The motion of the distal and proximal segments were tracked using two inertial 151 measurement units (Analog Devices, ADIS16354). The forces applied to the segment were measured using 152 a hand-held six degrees of freedom (DOFs) load-cell (ATI Industrial Motion, mini45). The position of the 153 load-cell relative to the joint axis was manually measured by the examiner.
Muscle strength, selectivity, and ROM were evaluated (Table 1)  We processed the experimental gait and IPSA data, used as input for the estimation of muscle-tendon   (six between the pelvis and the ground; three at each hip joint; one at each knee,   170 ankle, and subtalar joint; and three at the lumbar joint), 86 muscles actuating the lower limbs (43 per leg), 171 three ideal torque actuators at the lumbar joint, and four contact spheres per foot (Delp et al. (1990,2007)). 172 We added passive torques to the joints of the lower limbs and the trunk to model the role of the ligaments where gait2392 refers to the OpenSim gait2392 model (Delp et al. (1990,2007)). 197 The personalized parameters reflect the muscle force generating capacity of the subject. Only optimal 198 fiber lengths and tendon slack lengths were personalized as gait simulations have been shown to be the problem identifies muscle excitations that reproduce joint torques underlying a given movement while 203 minimizing a performance criterion (e.g., muscle effort). We augmented this formulation in different ways.

204
First, we added optimal fiber lengths and tendon slack lengths as optimization variables. Second, we 205 introduced a term in the cost function minimizing the difference between muscle activations and scaled 206 EMG signals where scale factors were included as optimization variables. Third, we assumed that muscles 207 operate around their optimal fiber lengths, and that maximal and minimal fiber lengths across movements should hence be larger and smaller, respectively, than their optimal fiber lengths. Fourth, we assumed 209 that resistance encountered when evaluating the ROM during the clinical examination may be, at least 210 in part, attributed to passive muscle forces. Hence, we included a term in the cost function minimizing 211 the difference between fiber lengths at these extreme positions of the ROM and reference fiber lengths 212 generating large passive forces. Finally, we minimized optimal fiber lengths, assuming that children with 213 CP have short fibers (Barrett and Lichtwark (2010)). The problem thus consisted in identifying muscle 214 excitations and parameters that minimized a multi-objective cost function: Passive forces in extreme positions force rate) to feedback muscle activations a s through a first order differential equation: where T s is a feedback threshold, g s is a feedback gain, and τ s = 30 ms is a time delay. 236 We determined the threshold for force feedback as the value 20 ms before the EMG onset (Staude 237 and Wolf (1999)) and used a zero threshold for force rate feedback. We identified the personalized 238 feedback gains that minimized the difference between EMG and feedback muscle activations during fast 239 passive stretches (IPSA measurements). We performed such optimization for the right medial hamstrings 240 (i.e., biceps femoris long head, semitendinosus, and semimembranosus) and for the right gastrocnemii 241 (i.e., gastrocnemius lateralis and medialis). We used semitendinosus EMG to drive the three hamstrings and gastrocnemius lateralis EMG to drive both gastrocnemii. We normalized EMG using scale factors 243 identified when estimating the personalized muscle-tendon parameters. We described the optimization 244 process in detail in previous work (Falisse et al. (2018)). We incorporated the spasticity model with 245 personalized feedback gains in our framework for predictive simulations (Figure 1). Since we only had 246 IPSA measurement for the right leg, we used feedback gains and thresholds identified with right leg data 247 for left leg muscles. Gait EMG data and spasticity, as clinically assessed (Table 1), were comparable for 248 both legs.

249
Muscle synergies 250 We modeled the reduced neuromuscular control complexity through muscle synergies. These synergies where a has dimensions N m × N f . Importantly, we did not impose subject-specific synergies when  Problem formulation 262 We predicted gait patterns by optimizing a gait-related cost function, independent of experimental data, 263 based on the MRI-based musculoskeletal model described above. In addition to optimizing performance, 264 we imposed average gait speed and periodicity of the gait pattern. We optimized for a full gait cycle to 265 account for asymmetry of CP gait. We solved the resultant optimal control problem via direct collocation.

266
The problem formulation and computational choices are detailed in previous work (Falisse et al. (2019)).

267
The cost function represents the goal of the motor task. We modeled this task-level goal as a weighted 268 sum of gait-related performance criteria including metabolic energy rate, muscle fatigue, joint accelerations, 269 passive joint torques, and trunk actuator excitations: Metabolic energy rates + w 2 a 10 10 Muscle fatigue Trunk actuator excitations where t f is unknown gait cycle duration, d is distance travelled by the pelvis in the forward direction,Ė are  2016)). We manually adjusted the weight factors until we found a set of weights that predicted human-like walking: w 1 = (25/86/body mass) × 10 −2 , 277 w 2 = 25/86 × 10 2 , w 3 = 50/21, w 4 = 10/15 × 10 2 , and w 5 = 1/3 × 10 −1 . We added several path 278 constraints enforcing a prescribed average gait speed corresponding to the child's average gait speed 279 (d/t f = 1 m s −1 ), imposing periodic states over the complete gait cycle (except for the pelvis forward 280 position), and preventing inter-penetration of body segments.

281
Model-based analyses 282 We investigated the differential effects of altered muscle-tendon properties, reduced neuromuscular 283 control complexity, and spasticity on gait patterns predicted with the MRI-based musculoskeletal model 284 ( Figure 1). In particular, we compared predicted joint kinematics and kinetics, muscle activity, and stride 285 lengths to their experimental counterparts. We also evaluated how impairments affected the metabolic cost 286 of transport (COT), defined as metabolic energy consumed per unit distance traveled.

287
First, we tested the influence of altered versus unaltered muscle-tendon properties by using personalized 288 versus generic muscle-tendon parameters in the muscle dynamics (Q i in Figure 1). In this initial analysis, 289 we did not include spasticity, nor imposed synergies. isolation as well as in combination with altered muscle-tendon properties.

296
Finally, we evaluated the effect of spasticity in the three medial hamstrings and two gastrocnemii of both 297 legs (Q iii in Figure 1). We modeled muscle activations as the sum of reflex muscle activations determined 298 based on the personalized spasticity model and feedforward muscle activations: where a f f are feedforward muscle activations, and a F t and a dF t are muscle activations from muscle force 300 and force rate feedback, respectively. We only tested the effect of spasticity based on the model with 301 personalized muscle-tendon parameters, since these parameters were used to estimate the feedback gains. 302 We tested the effect of spasticity in combination with selective control (i.e., no synergy constraints) as well 303 as with a reduced number of muscle synergies.

305
As an additional analysis, we investigated whether the child adopted an impaired crouch gait pattern 306 because of neuro-mechanical constraints or because it was more optimal (Q iv in Figure 1). To this aim, we 307 added a term in the cost function that penalized deviations from measured kinematics of a TD child: whereq are measured joint positions of a TD child and w 6 = 100/20 is a weight factor. We generated 309 these simulations with personalized parameters as well as with and without synergies. We did not include accounting for synergies, we added initial guesses derived from simulated kinematics with the lowest 320 optimal costs produced without synergies and with more synergies (e.g., with three synergies, initial 321 guesses were derived from the best kinematic solutions with four synergies and without synergies). For 322 simulations accounting for spasticity, we added initial guesses derived from simulated kinematics with 323 the lowest optimal costs produced without spasticity. In all cases, initial guesses for muscle, trunk, and 324 synergy variables were constant across time and not informed by experimental data. Initial guesses for 325 synergy weights were constant across muscles and independent of experimental data.

Gait analysis 327
The child walked with a pronounced crouch gait pattern characterized by bilateral knee extension deficits 328 with reduced knee ROM during swing, a lack of right ankle dorsiflexion at the end of swing, excessive left 329 ankle dorsiflexion, excessive and deficient right and left hip adduction, respectively, and excessive bilateral 330 hip internal rotation (Figures 2 and S1).

332
Using personalized versus generic muscle-tendon parameters resulted in a crouch (i.e., excessive knee 333 flexion) versus a more upright gait pattern (Figures 2 and S1; Movies S1-2). Personalized optimal fiber 334 lengths and tendon slack lengths were generally smaller and larger, respectively, than their generic 335 counterparts (Tables S1-2). The use of personalized parameters resulted in decreased deviations (smaller

342
Predicted stride lengths were larger than the average stride length of the child but were within two standard 343 deviations.

344
Influence of the synergies with generic muscle-tendon parameters 345 Reducing the number of synergies in combination with generic muscle-tendon parameters did not induce 346 the amount of crouch that was experimentally measured in the child, although it altered muscle coordination and increased COT (Figures 3 and S2; Movie S3). The right knee flexion angles increased during stance 348 with the reduction of the neuromuscular control complexity but were still smaller than experimentally Figure 2. Influence of the muscle-tendon parameters on the predicted walking gaits. Variables from the right leg are shown over a complete gait cycle; left leg variables are shown in Figure S1 (Supplementary Material). Vertical lines indicate the transition from stance to swing. Experimental data is shown as mean ± two standard deviations. Experimental EMG data was normalized to peak activations. GRF is for ground reaction forces; BW is for body weight; COT is for metabolic cost of transport; lh is for long head.      Figure S4 (Supplementary Material). When accounting for spasticity, total activations (green) combine spastic (solid black) and nonspastic (dotted black) activations. Vertical lines indicate the transition from stance to swing. Experimental data is shown as mean ± two standard deviations. Experimental EMG data was normalized to peak activations. Lh is for long head.

384
Tracking the TD kinematics while using personalized muscle-tendon parameters produced an upright gait 385 pattern when not incorporating synergies, but decreased the overall gait performance (Figures 6 and S5; 386 Movie S6). Specifically, the simulated gait had a similar COT (3.46 J kg −1 m −1 ) as the crouch gait pattern 387 predicted without such tracking term but the contribution of most terms in the cost function increased, 388 suggesting that walking upright is not prevented by mechanical constraints (i.e., aberrant musculoskeletal 389 geometries and altered muscle-tendon properties) but is less optimal, due to these mechanical constraints, 390 than walking in crouch for this child. The contribution of the muscle fatigue term increased by 29%, in 391 part driven by higher activations of the glutei. The contribution of the joint acceleration, metabolic energy 392 rate, and passive joint torque terms increased by 15, 15, and 36%, respectively, when walking upright.

393
Similarly, passive muscle forces increased when walking upright for the iliacus and psoas (hip flexors), 394 and biceps femoris short head (knee flexor). Knee flexion increased when adding synergies but did not 395 reach the angle that was experimentally measured in the child ( Figure S6). Nevertheless, this suggests that 396 reduced neuromuscular control complexity may contribute to crouch gait. The gastrocnemius lateralis and 397 soleus (ankle plantarflexors) were also activated earlier during stance with synergies. Imposing synergies 398 increased the COT (4.12 and 4.05 J kg −1 m −1 with four and three synergies, respectively).
399 Figure 6. Influence of tracking the TD kinematics on predicted walking gaits. Variables from the right leg are shown over a complete gait cycle; left leg variables are shown in Figure S5 (Supplementary Material). Vertical lines indicate the transition from stance to swing. Experimental data is shown as mean ± two standard deviations. Muscle fatigue is modeled by activations at the tenth power. Passive muscle forces are normalized by maximal isometric muscle forces. Sh is for short head. The influence of synergies on predicted walking gaits is depicted in Figure S6 (Supplementary Material).

DISCUSSION
We demonstrated the ability of predictive simulations to explore the differential effects of musculoskeletal to the crouch gait pattern. These simulations thus suggest that adopting an upright gait pattern for this 407 child might produce an early onset of fatigue, which might explain in part why the child walks in crouch.

408
Importantly, not only fatigue, but also joint accelerations, passive joint torques, and metabolic energy rates were also reported to be spastic (Table 1). Including these other muscles may have an influence on walking 441 performance. Second, experimental data from the spasticity assessment was only collected for the right leg, 442 whereas bilateral spasticity was reported (Table 1). We optimized the feedback parameters using that data 443 but used the resulting parameters for both legs, which might affect our predictions. Third, we used feedback 444 parameters optimized from passive stretches to predict spasticity (i.e., reflex activity) during gait, assuming Other factors might have contributed to the deviations between predicted and measured movements.

478
First, the musculoskeletal model had generic rather than subject-specific (i.e., MRI-based) geometries 479 for feet and tibias. Since the child later underwent a surgery that included bilateral tibia derotation, muscle-tendon parameters in order to obtain valid parameter estimates (i.e., representative of the subject).

496
In this study, the available experimental data was limited to walking trials and passive stretches from one 497 leg. Hence, it is likely that some parameters were calibrated to fit the experimental data but did not truly 498 reflect the force-generating capacities of the child. When used in conditions different from the experiments, 499 these parameters may hence result in non-representative force predictions. A challenge for upcoming 500 research will be the design of experimental protocols to collect experimental data that contains sufficient

DATA AVAILABILITY STATEMENT
All data, code, and materials used in this study will be made available at https://simtk.org/ projects/predictcpgait upon publication.
result in a small number of muscle synergies during gait. Frontiers in Computational Neuroscience 8,