Muscle weakness but also contractures contribute to the progressive gait pathology in children with Duchenne muscular dystrophy: a simulation study

Abstract

Background

Muscle weakness and contractures cause gait deficits in children with Duchenne muscular dystrophy (DMD) but their relative contributions are poorly understood and hence it is unclear whether contractures should be treated. Therefore, we aimed to differentiate the effect of muscle weakness in isolation from weakness and contractures combined on the gait patterns.

Methods

We used computer simulations that generate gait patterns based on a musculoskeletal model (without relying on experimental data) to establish the relationship between muscle impairments and gait deviations. We previously collected a longitudinal database of 137 repeated measurements in 30 boys with DMD and found that the data measured through 3D gait analysis could be clustered in three gait patterns. We estimated weakness based on data from fixed dynamometry, and contractures based on goniometry and clinical measures. Foot deformities were modeled by reducing the height of all foot segments and decreasing the strength of the intrinsic foot muscles. We created musculoskeletal models that either represented (1) the mean weakness; (2) the mean weakness and contractures; or (3) the mean weakness, contractures and foot deformities, in each gait pattern.

Results

Simulations based on models with both weakness and contractures captured most (but not all) experimentally observed gait deviations, demonstrating the validity of our approach. While muscle weakness was primarily responsible for gait deviations, muscle contractures and foot deformities further contributed to gait deviations. Interestingly, the simulations predict that the combination of increasing weakness and contractures rather than increasing weakness alone causes loss of ambulation for the most affected gait pattern.

Conclusions

Predictive simulations have the potential to elucidate causal relationships between muscle impairments and gait deviations in boys with DMD. In the future, they could be used to design targeted interventions (e.g. stretching, assistive devices) to prolong ambulation.

Background

Duchenne muscular dystrophy is a severe X-linked neuromuscular disorder caused by mutations in the dystrophin gene, leading to progressive muscle degeneration characterized by loss of contractile tissue and its replacement by fat and fibrotic tissue [1,2,3]. Subsequently, progressive muscle weakness, stiffness and contractures manifest, altering posture and gait, with loss of ambulation between 7.1 and 18.6 years of age (mean age: 12.7 years) [1, 4, 5]. Since there is no cure, rehabilitation, orthopedic, and orthotic treatments target muscle stiffness and contractures in order to prolong ambulation and slow disease progression [1]. However, conflicting results of orthopedic and orthotic treatments on gait decline have been reported [6,7,8,9,10]. Additionally, disease-modifying treatments show promise in further delaying the loss of ambulation and slowing disease progression. However, progress and clinical implementation have been hampered by the lack of sensitive outcome measures [11,12,13,14]. Improved understanding of how impairments contribute to gait pathology is necessary to enhance clinical decision-making, advance current rehabilitation, orthopedic and orthotic treatments, and support the development of novel treatments. Yet, insights into the causal relationships between underlying weakness and contractures, and gait pathology in DMD are still limited.

It is poorly understood how muscle weakness and stiffness interact in progressive DMD gait pathology. Although muscle weakness is considered the most important contributor to pathological gait in DMD [1, 5], the role of contractures in conjunction with weakness remains unclear. It has been postulated that boys with DMD could use the passive forces caused by contractures to compensate for muscle weakness [15], but also that increased contractures could lead to loss of ambulation if they undermine successful compensation mechanisms for muscle weakness [1]. We recently provided objective evidence that muscle weakness and contractures are associated with, and interact longitudinally with, gait pathology in DMD [16]. However, the co-occurrence and simultaneous decline of multiple impairments hamper disentangling the causal contribution of weakness from that of contractures to gait pathology in DMD based on experimental data alone.

Three gait patterns have been identified in growing children with DMD: a mildly affected gait pattern, a tiptoeing gait pattern, and a flexion gait pattern with distinct anterior pelvic tilt and posterior trunk leaning [17]. Both the tiptoeing and flexion gait patterns deviated more severely from the typically developing (TD) gait pattern than the mildly affected gait pattern, but they differed in the location (proximal versus distal) of the most severe deviations. The tiptoeing gait pattern exhibited the most severe deviations around the ankle, while the flexion gait pattern exhibited the most severe deviations around the proximal joints. Notably, abnormally large midfoot motion (i.e. midfoot break) was observed in boys with DMD exhibiting the flexion gait pattern.

Physics-based predictive simulations of gait are a powerful tool to unravel complex causal relationships between underlying impairments and pathological gait patterns. Such simulations generate novel movement patterns based on a musculoskeletal model (i.e., a mathematical description of the musculoskeletal system) and the assumption that gait is achieved by optimizing performance, without relying on experimental data [18, 19]. They therefore enable computation of the direct effect of a change in the musculoskeletal system on the gait pattern by adjusting model parameters, allowing to study both the isolated effects as well as interaction effects between underlying impairments on the gait pattern. Falisse et al. [20] improved the computational efficiency of forward simulations of 3D gait [20,21,22,23] allowing us to perform many simulations with complex 3D models [24]. D’Hondt et al.[23, 25] demonstrated that incorporating a dynamic three-segment foot model with flexible longitudinal arch into a complex 3D model (based on OpenSim’s gait2392 model) improved gait simulations and especially the agreement between simulated and experimentally observed ankle–foot kinematics. In addition, the more detailed foot model allows for modeling foot deformities. Since boys with DMD exhibit severe gait deviations affecting all planes and could present a midfoot break, complex 3D models are necessary to accurately simulate DMD gait.

Previous predictive simulations studies have examined the effects of weakness and/or stiffness on gait, but the interaction between these impairments has thus far only been investigated for the ankle. Ong et al. [26] showed that the isolated effect of severe plantar flexor weakness resulted in calcaneal gait without crouch, while the isolated effect of severe plantar flexor contracture resulted in equinus gait. Waterval et al. [27] found that a simulated gait pattern with 80% of plantar flexion weakness closely matched the experimental hip and ankle kinematics and kinetics of participants with bilateral plantar flexion weakness, thereby validating their predictive simulations. By investigating the isolated effect of incrementing plantar flexion weakness, they observed that pathological gait features particularly emerged when bilateral plantar flexion weakness exceeded 60%. However, these studies employed 2D musculoskeletal models, whereas gait deviations in DMD also occur in the frontal plane. In contrast, Falisse et al. [20] performed predictive simulations with complex 3D models and showed that increased hip weakness caused increased trunk sway, increased step width and decreased hip moments, and that increased plantar flexion weakness caused increased knee flexion, increased ankle dorsiflexion and decreased stride lengths that reduced ankle moments. Previous studies have thus been limited to examining the isolated effects of weakness or contracture in separate muscle groups and have not yet investigated the interaction between weakness and contractures on the gait pattern, except for the ankle. Additionally, there is still a need to evaluate whether predictive simulations of 3D gait could capture the key features of DMD gait and whether modeling increasing levels of impairments would result in the different DMD gait patterns, based on severity and location of the impairments [17].

The first aim of this study was to investigate the ability to simulate DMD gait using models that reflect the DMD-specific weakness and contractures. Specifically, we aimed to evaluate whether we can capture the mildly affected gait pattern, the tiptoeing gait pattern, and the flexion gait pattern [17] by modeling the experimentally measured muscle weakness and contractures of the children who exhibit these gait patterns. Given that the interaction between weakness and contractures is only poorly understood in DMD, the second aim of this study was to differentiate the effect of weakness alone from the effect of weakness in combination with contractures on the simulated gait patterns. Lastly, the third aim of this study was to investigate the effect of progressive impairments on the pathological gait pattern by modeling increasing levels of impairment severity.

Methods

We performed a simulation study based on our previously collected longitudinal database[17, 28, 29].

Experimental data

We used experimental data to create models that reflect the muscle impairments in children with DMD and to evaluate how accurately the simulations (see "Personalizing musculoskeletal model" section) capture the gait deviations.

Our collected longitudinal dataset consists of instrumented strength assessments, clinical examinations, and 3D gait analyses of 30 boys with DMD. These boys were repeatedly measured between 2015 and 2022 at multiple time points (median number of repeated sessions = 4; range number of repeated sessions = 1–10; mean time between repeated sessions = 0.74 ± 0.49 years) at the Clinical Motion Analysis Laboratory in the University Hospital Leuven campus Pellenberg, resulting in a dataset of 137 sessions. The data and collection procedures have been extensively described in our previous work [17, 28, 29].

Data collection was approved by the local ethics committee (Ethical Committee UZ Leuven/KU Leuven; S61324) under the Declaration of Helsinki. The parents or participants’ caregivers provided written informed consents. Participants aged 12 years or older provided informed assents. All methodology adhered to the relevant regulations and guidelines.

Body mass, height, and lower limb segment lengths were measured at each observation.

Hip extension, flexion and abduction, knee extension and flexion, and ankle plantar flexion and dorsiflexion muscle weakness were assessed with an instrumented strength assessment [30, 31]. To avoid fatigue and ensure participant cooperation, the instrumented strength assessment was performed unilaterally. If asymmetry was detected on the manual muscle testing during the clinical exam (not common), the weakest side was assessed. If muscle weakness was symmetrical, the assessed side was randomly chosen by flipping a coin. For each child, the same leg was assessed consistently across all longitudinal assessments. Participants performed maximal voluntary isometric contractions (MVIC) in standardized test positions on a fixed dynamometer (MicroFet, Hogan Health Industries, West Jordan, UT United States). The mean maximal joint torque per muscle group was calculated by multiplying the mean maximal force over one to three representative MVIC trials with its lever arm with respect to the joint. The strength outcomes of the children with DMD were expressed as percentages of the mean strength of TD children with a similar body mass and height (based on a dataset of 153 children [32]), reflecting strength deficits.

At each observation a standardized clinical examination was performed. Passive range of motion (ROM) of hip extension (modified Thomas test [33]), hip adduction (with extended hip and knee of the assessed leg and hip and knee flexed in 90° of the contralateral leg), the hamstrings (true popliteal angle [33]), and ankle dorsiflexion (with knee extended and knee flexed in 90° [33]) was assessed using goniometry. The age-related normative values of Mudge et al. [33] and Sankar et al. [34] were used as the reference. Differences between the ROM measurements of the children with DMD and the reference ROM measurements of TD children with a similar age were calculated, reflecting ROM deficits. Muscle stiffness of hip flexors, hip abductors, rectus femoris, hamstrings, plantar flexors (with knee extended and flexed) was also measured using a clinical stiffness scale during passive elongation (0–3 point ordinal scale: 0 = no increased resistance; 1 = minimal increased resistance at the end of ROM; 2 = increased resistance; 3 = highly pronounced resistance). This scale was inspired by the Ashworth scale to assess joint hyper-resistance in neurological diseases for which validity and reliability have been established [35, 36]. Given that there is no contribution from abnormal tone to increased resistance in DMD, we removed references to neurally induced increases in resistance such as ‘catch’. Both sides were assessed, but only the values from the side measured with the instrumented strength assessment were included for further analysis. Weakness of the abdominal and back muscles was evaluated with the manual muscle testing [37]. Strength of the abdominal and back muscles was converted to strength deficits, expressed as a percentage relative to the maximal score of 5.

Gait was measured by 3D gait analysis. The boys with DMD were instrumented with reflective markers (Plug-In Gait Full-body marker model with simple foot model (two foot markers); marker diameter: 14 mm), whose trajectories were recorded with a 10–15 Vicon camera system (Vicon-UK, Oxford, UK; sampling frequency of 100 Hz; built-in Woltring filter with mode MSE and smoothing of 15 mm²) during barefoot walking at self-selected speed on a 10-m walkway. Two embedded force plates (AMTI, Watertown, MA, USA; sampling frequency: 1500 Hz) captured ground reaction forces. Nexus software (Nexus 2.10. Vicon-UK, Oxford, UK) was used to define gait cycles and estimate trunk, pelvis and lower limb kinematic waveforms, lower limb kinetic waveforms, and spatiotemporal parameters. Ten gait cycles with kinematic data, of which three to five with kinetic data, were collected bilaterally. Quality of collected gait cycles was checked in a custom-made MATLAB software (The Mathworks Inc., Natick, M.A., 2016 and 2021b). We averaged the kinematic waveforms, kinetic waveforms and spatiotemporal parameters of the selected gait cycles with good quality per observation separately for the right and left sides. Each observation was classified into the mildly affected gait pattern, the tiptoeing gait pattern, or the flexion gait pattern according to our previously introduced gait classification [17]. Gait data of nineteen TD boys (Plug-In Gait Lower-body marker model; marker diameter: 14 mm) with a similar age range as the boys with DMD was selected from our reference database at the University Hospital Leuven. Data for the trunk segment was not collected in the TD boys, as only the lower-body marker model was applied.

The entire database of strength assessments, clinical examinations, and 3D gait analyses was divided into three DMD groups based on the label of the gait classification [17]: “DMD1”, i.e., the mildly affected gait pattern, “DMD2”, i.e., tiptoeing gait pattern, and “DMD3”, i.e., flexion gait pattern (Fig. 1). The selected TD gait dataset was used to represent the TD group. For each group, the mean and standard deviation (SD) were calculated for the gait data. For the DMD groups, the mean and SD were also calculated for the strength deficits, the ROM deficits and the clinical stiffness scale. These values—mean, one SD above the mean and one SD below the mean—were later used as input to the musculoskeletal models.

Fig. 1

Experimental data used to personalize musculoskeletal models. A Anthropometry (mean ± 1 standard deviation). B Muscle strength deficits. C Passive range of motion. D Muscle stiffness. AB, abdominal; BA, back; DF, dorsiflexors; DMD1, DMD group with the mildly affected gait pattern; DMD2, DMD group with the tiptoeing gait pattern; DMD3, DMD group with the flexion gait pattern; GAS, gastrocnemius; HA, hip abductors; HAM, hamstrings; HE, hip extensors; HF, hip flexors; KE, knee extensors; KF, knee flexors; PF, plantar flexors; RF, rectus femoris; SOL, soleus; TD, typically developing

Personalizing musculoskeletal model

Generic musculoskeletal model

We used the musculoskeletal model proposed by D’Hondt et al. [23]. This model has 33 degrees of freedom (dofs) (pelvis-to-ground: 6 dofs, hip: 3 dofs, knee: 1 dof, ankle: 1 dof, subtalar: 1 dof, midtarsal: 1 dof, metatarsophalangeal: 1 dof, lumbar: 3 dofs, shoulder: 3 dofs, and elbow: 1 dof). 94 Hill-type muscle–tendon units (92 muscles according to OpenSim’s gait2392 model and the right and left intrinsic foot muscles) actuate the lower limb and lumbar joints. A right and left plantar fascia is included. Eight ideal torque actuators actuate the shoulder and elbow joints. Five Hunt-Crossley contact spheres per foot model the foot–ground contact. Passive joint torques with exponential stiffness and damping [38] were added to the lower limb and lumbar joints to represent the effects of unmodeled passive structures [20]. An adjustment we made compared to D’Hondt et al. [23] was fitting the exponential stiffness of the lumbar joint to experimental data of Panjabi et al. [39]. Muscle excitation-activation coupling was described by Raasch’s model [40, 41] and muscle–tendon interaction and the dependence of muscle force on fiber length and velocity by a Hill-type muscle model [42, 43]. Skeletal motion was modeled with Newtonian rigid body dynamics [44, 45].

Scaling of generic musculoskeletal model to child’s dimensions

For each DMD group and the TD group, the generic musculoskeletal model was linearly scaled to the mean anthropometry of that group, using the OpenSim scaling tool. For the TD group, the mean length and body mass of the TD group were therefore used. Since the boys with DMD have an increased body mass due to extra fat tissue, we first scaled the model parameters to the mean length of the DMD group and the expected mean body mass of TD children of the same length (i.e., expected TD body mass). The additional mass of the children with DMD, i.e. difference between the mean body mass of the DMD group and expected TD body mass was then divided with 2/3 over the trunk and pelvis segments, and with 1/3 over the lower limb segments. This decision was based on the observation that the additional mass was predominantly distributed around the abdomen and buttocks, with secondary accumulation in the legs due to increased intramuscular and subcutaneous fat. Notably, no excess mass was observed in the feet. Subsequently, we scaled the model parameters that are also dependent on anthropometry but were not scaled by the OpenSim scaling tool (e.g. muscles’ maximal isometric force) according to geometric similarity (Additional file 1).

Modeling impairments

The experimental data of the selected side was used as input for both the left and right sides of the musculoskeletal models to create symmetrical models. We chose this approach as muscle impairments typically manifest symmetrically in DMD [46].

The muscles in the model were represented as Hill-type muscle–tendon units (for a detailed description see [42, 43]). The muscle–tendon unit consists of an active contractile element in parallel with a passive element, which is in series with a tendon. The muscle force arises from both the active contractile component and the passive elastic element as:

$${F}_{m}={F}_{m}^{max}*\left({f}_{m}^{act}({\widetilde{l}}_{m} , {\widetilde{v}}_{m} ,a)+ {f}_{m}^{pass}({\widetilde{l}}_{m})\right)$$

(1)

where ${F}_{m}$ is muscle force, ${F}_{m}^{max}$ is maximal isometric force, ${f}_{m}^{act}$ is the active muscle force-length-velocity characteristic, ${\widetilde{l}}_{m}$ is normalized fiber length, ${\widetilde{v}}_{m}$ is normalized fiber velocity, $a$ is muscle activation, and ${f}_{m}^{pass}$ is the (non-linear) passive muscle force-length characteristic (see [43] for mathematical expression and visual representation of the characteristics).

The most common parametrization of this model [42] assumes that maximal isometric force (${F}_{m}^{max}$), and passive muscle (see Eq. (1)) and tendon stiffness are coupled. Therefore, they all scale with ${F}_{m}^{max}$. However, in DMD, active and passive muscle forces do not decrease simultaneously. The loss of contractile tissue is accompanied by its replacement with fat and fibrotic tissue, resulting in a decline in active muscle force while passive muscle stiffness increases. Therefore, we modeled muscle weakness by scaling only the active force component (${f}_{m}^{act}$), rather than scaling ${F}_{m}^{max}$ that also scales the passive elements:

$${f}_{m, DMD}^{act}={f}_{m, TD}^{act}* \left(\frac{{MVIC}_{DMD}}{{MVIC}_{TD}}\right)$$

(2)

where MVIC_DMD/MVIC_TD is either the mean, one SD above the mean, or one SD below the mean of the experimentally measured strength deficit.

In DMD, contractile tissue is not only lost but also replaced by fat and fibrotic tissue, resulting in increased muscle stiffness and eventually leading to contractures. We modeled this by shifting the passive muscle force-length relationship to shorter fiber lengths through a reduction in the fiber length at which passive muscle force begins to develop. Neither slack length nor fascicle length was altered; only the position of the passive force-length curve was adjusted. We used the ROM measurements and clinical stiffness scale to estimate this shift. For the ROM measurements, we estimated the difference in fiber length at which the muscle starts to develop passive force between TD and DMD from the difference in joint angle at the end of ROM. The joint angle at the end of ROM of TD children was based on age-related reference data reported by Mudge et al. [33] and Sankar et al. [34]. To estimate the corresponding difference in fiber length, we multiplied the difference in measured joint angle at end ROM between TD and DMD (in radians) with the moment arm of the muscles in the anatomical position. This difference in fiber length was normalized to optimal fiber length to compute the shift of the passive force-length relationship. For the clinical stiffness scale, the normalized fiber length at which passive force starts to develop was assumed 1 when the clinical stiffness score was 0 (no increased resistance), 0.83 when the score was 1 (minimal increased resistance), 0.67 when the score was 2 (increased resistance), and 0.5 when the score was 3 (highly pronounced resistance) corresponding to a shift of respectively 0, 0.17, 0.33, and 0.5. We shifted the passive force-length relationship by the mean of the shifts estimated based on the ROM and clinical stiffness score.

Boys in DMD3 presented a midfoot break on the videos, a feature associated with reduced arch height and arch stiffness. To model this foot deformity, we reduced the height of all foot segments with 10% and decreased the strength of the intrinsic foot muscles with 50% [23], since activation of the intrinsic foot muscles stiffens the arch and DMD is primarily a muscle disorder.

Predictive simulation of gait

Our previously described optimal control framework was used to simulate gait [20, 21]. We formulated simulations of gait as optimal control problems. Gait can then be simulated by solving for muscle controls that minimize a cost function, while imposing task constraints and the dynamics of a musculoskeletal model, without relying on experimental data. The task constraints for walking were the average forward speed of the pelvis and right-left symmetry. For each symmetrical model, we imposed the average experimentally measured gait speed. As contact between segments was not explicitly modeled, distance constraints were used to prevent segments to penetrate each other [20]. We scaled the imposed distances based on body height (Additional file 2).

We used a previously determined cost function, i.e. the integral of the weighted sum of squared metabolic energy rate $(\dot{E})$, muscle activations (a), joint accelerations (u_a), and passive joint torques (T_p):

$$J=\frac{1}{d}{\int }_{0}^{{t}_{f}}\left({w}_{1}*{\dot{E}}^{2}+{w}_{2}*{a}^{2}+{w}_{3}*{u}_{a}^{2}+ {w}_{4}*{T}_{p}^{2}\right)*dt$$

(3)

where d is the distance traveled, t_f is half gait cycle duration, t is time, and w_1-4 are the weight factors. Weight factors w₁ and w₄ were scaled based on anthropometry (Additional file 2).

This cost function was selected based on the agreement between simulated and experimental gait patterns at self-selected speed in a healthy individual and evaluated in other conditions (range of gait speeds; weakness; prosthesis use) [20]. The model of Bhargava et al. [47] was used to calculate the metabolic energy of the muscles and was made continuously differentiable by approximating conditional statements with a hyperbolic tangent.

Our framework is implemented in MATLAB (The Mathworks Inc., USA). Skeletal dynamics is formulated using OpenSim (based on version 4.3), CasaDi [48] is used to formulate the optimal control problem and calculate derivatives, and IPOPT [49] to solve the optimal control problem. The code and documentation of the user friendly PredSim code that was used for this study can be found here [22]: https://github.com/KULeuvenNeuromechanics/PredSim. We used the gait data of a TD individual (available on github “IK_Bounds_Default.mot”) as the initial guess for all simulations.

Analyses

Forward simulations were conducted based on a model that represented the mean weakness and contractures (± SD) in each group (and thus not for individual patients). For each DMD group, simulations were performed based on six models to study how weakness and contractures with different severities impact the walking pattern: (1) mean weakness; (2) weakness one SD below the mean; (3) weakness one SD above the mean; (4) mean weakness and mean contractures; (5) weakness and contractures one SD below the mean; and (6) weakness and contractures one SD above the mean. For DMD3, one additional simulation was performed, incorporating the midfoot break in addition to the mean weakness and mean contractures. For the TD group, a single predictive simulation on a model without impairments was conducted, serving as the reference. Contractures and foot deformities were not simulated in isolation as muscle weakness is the primary clinical symptom of DMD and the aim was to improve understanding of its interaction with contractures and foot deformities.

To validate whether the predictive simulations can capture DMD gait, we compared the experimental differences in gait kinematics and kinetics between DMD and TD with the differences in gait kinematics and kinetics simulated based on the models with impairments and the TD model. We focused on comparing experimental and simulated differences between TD and DMD kinematics and kinetics, as some of the differences between experimental and simulated gait patterns might be due to the use of different kinematic models.

To study the interaction between weakness and contractures, we compared the simulation based on the mean weakness model with the simulation based on the mean weakness and contractures model, for each DMD group. For DMD3, we also compared the simulation of both models to the simulation based on the mean weakness, contractures and midfoot break model to evaluate the additional effect of modeling foot deformities.

To evaluate the effect of progressive impairments, we compared the simulations based on the -1SD (i.e., one SD below the mean) and + 1SD (i.e., one SD above the mean) models to the simulations based on the mean models.

Results

Simulations based on models that reflect DMD-specific impairments capture DMD-specific gait impairments

For DMD1, the predictive simulation based on a model reflecting the mean weakness and contractures in this group captured the minor gait deviations between DMD1 and TD (Figs. 1–4). The experimental gait pattern of DMD1 deviated only minimally from TD gait. The main differences in the sagittal plane kinematics were increased anterior pelvic tilt, reduced hip extension at the end of stance, increased hip flexion during swing, reduced knee flexion motion during loading response, delayed knee flexion and increased peak knee flexion during swing, and reduced dorsiflexion during swing (Fig. 2). The increased anterior pelvic tilt and increased hip flexion during swing were minimally captured in the simulations. The reduced knee flexion motion during stance was captured, but the simulations predicted more flexion.

Fig. 2

Comparison of simulated and experimental sagittal plane kinematics for DMD1. A Simulated kinematics when only modeling weakness. B Simulated kinematics when modeling weakness and contractures. C Experimental gait kinematics. D Differences in kinematics between the TD model and the mean weakness model (light green), the TD model and the mean weakness + contractures model (dark green), and the experimental TD and DMD1 data (blue). Ant, anterior; DF, dorsiflexion; Ext, extension; Flex, flexion; PF, plantar flexion; Post, posterior; TD, typically developing

The main differences in the frontal and transverse plane kinematics were increased pelvic upwards position at initial contact, increased hip adduction during stance, and increased hip abduction during swing (Fig. 3). These differences were not fully captured in the simulations.

Fig. 3

Comparison of simulated and experimental frontal and transverse plane kinematics for DMD1. A Simulated kinematics when only modeling weakness. B Simulated kinematics when modeling weakness and contractures. C Experimental gait kinematics. D Differences in kinematics between the TD model and the mean weakness model (light green), the TD model and the mean weakness + contractures model (dark green), and the experimental TD and DMD1 data (blue). Abd, abduction; Add, adduction; Ex, external; In, internal; TD, typically developing

The main differences in joint moments were decreased hip extension moment, premature hip flexion moment (i.e., an early negative hip extension moment during stance), decreased knee extension and flexion moments, and reduced ankle plantar flexion moment (Fig. 4). The premature hip flexion moment was not captured in the simulations. The reduced ankle plantar flexion moment was captured, but the experimental difference was larger.

Fig. 4

Comparison of simulated and experimental kinetics for DMD1. A Simulated kinetics when only modeling weakness. B Simulated kinetics when modeling weakness and contractures. C Experimental gait kinetics. D Differences in kinetics between the TD model and the mean weakness model (light green), the TD model and the mean weakness + contractures model (dark green), and the experimental TD and DMD1 data (blue). Abd, abduction; Add, adduction; DF, dorsiflexion; Ext, extension; Flex, flexion; PF, plantar flexion; TD, typically developing

For DMD2, the predictive simulation based on a model reflecting the mean weakness and contractures in this group captured many but not all deviations between DMD2 and TD (Fig. 1; Figs. 5–7). The experimental gait pattern of DMD2 was characterized by a tiptoeing pattern and deviated more from TD gait than DMD1. The main differences in the sagittal plane kinematics were continuous ankle plantar flexion, reversed first ankle rocker, horizontal second rocker, large increased plantar flexion during swing, absent knee flexion motion during loading response, delayed and increased peak knee flexion during swing, increased hip flexion at the end of stance and during swing, and increased anterior pelvic tilt (Fig. 5). The ankle pattern was captured in the simulations, yet the experimental difference in stance was larger. The absent knee flexion motion was not captured in the simulations. The increased hip flexion and anterior pelvic tilt were minimally captured in the simulations.

Fig. 5

Comparison of simulated and experimental sagittal plane kinematics for DMD2. A Simulated kinematics when only modeling weakness. B Simulated kinematics when modeling weakness and contractures. C Experimental gait kinematics. D Differences in kinematics between the TD model and the mean weakness model (light green), the TD model and the mean weakness + contractures model (dark green), and the experimental TD and DMD2 data (yellow). Ant, anterior; DF, dorsiflexion; Ext, extension; Flex, flexion; PF, plantar flexion; Post, posterior; TD, typically developing

The main differences in the frontal and transverse plane kinematics were increased range of pelvic upwards/downwards and external/internal motion, and increased hip abduction in stance and swing (Fig. 6). The increased range of pelvic upwards/downwards motion was not fully captured in the simulations. The increased range of external/internal motion was captured, but the experimental difference was larger.

Fig. 6

Comparison of simulated and experimental frontal and transverse plane kinematics for DMD2. A Simulated kinematics when only modeling weakness. B Simulated kinematics when modeling weakness and contractures. C Experimental gait kinematics. D Differences in kinematics between the TD model and the mean weakness model (light green), the TD model and the mean weakness + contractures model (dark green), and the experimental TD and DMD2 data (yellow). Abd, abduction; Add, adduction; Ex, external; In, internal; TD, typically developing

The main differences in joint moments were reduced hip extension moment, premature hip flexion moment, decreased knee extension and flexion moments, absent dorsiflexion moment, premature plantar flexion moment, and reduced peak plantar flexion moment (Fig. 7). Only the premature hip flexion moment and reduced knee extension moment were not captured in the simulations.

Fig. 7

Comparison of simulated and experimental kinetics for DMD2. A Simulated kinetics when only modeling weakness. B Simulated kinetics when modeling weakness and contractures. C Experimental gait kinetics. D Differences in kinetics between the TD model and the mean weakness model (light green), the TD model and the mean weakness + contractures model (dark green), and the experimental TD and DMD2 data (yellow). Abd, abduction; Add, adduction; DF, dorsiflexion; Ext, extension; Flex, flexion; PF, plantar flexion; TD, typically developing.

For DMD3, the predictive simulation based on a model reflecting the mean weakness, mean contractures and midfoot break in this group captured most but not all deviations between DMD3 and TD (Fig. 1, 8, 9, 10). The experimental gait pattern of DMD3 was characterized by a flexion pattern with distinct anterior pelvic tilt and posterior trunk leaning and deviated more from TD gait than DMD1. The main differences in the sagittal plane kinematics were increased anterior pelvic tilt, continuous increased hip flexion, delayed peak knee flexion during swing, absent knee flexion motion during loading response but increased knee flexion compared to DMD2, reduced dorsiflexion during swing and at initial contact (Fig. 8). The increased anterior pelvic tilt was not captured and the increased hip flexion was only minimally captured, mainly at the end of stance in the simulations. The delayed knee flexion during swing was captured, but simulations predicted also increased peak knee flexion. The reduced dorsiflexion during swing and at initial contact was captured, but simulations predicted more plantar flexion.

Fig. 8

Comparison of simulated and experimental sagittal plane kinematics for DMD3. A Simulated kinematics when only modeling weakness. B Simulated kinematics when modeling weakness and contractures. C Simulated kinematics when modeling weakness, contractures and midfoot break. D Experimental gait kinematics. E Differences in kinematics between the TD model and the mean weakness model (light green), the TD model and the mean weakness + contractures model (dark green), the TD model and the mean weakness + contractures + midfoot break model (medium green), and the experimental TD and DMD3 data (red). Ant, anterior; DF, dorsiflexion; Ext, extension; Flex, flexion; PF, plantar flexion; Post, posterior; TD, typically developing

The main differences in the frontal and transverse plane kinematics were increased range of pelvic upwards/downwards and external/internal motion, but less pronounced than in DMD2, and increased hip abduction in swing (Fig. 9). The simulations predicted more hip abduction during stance than experimentally observed.

Fig. 9

Comparison of simulated and experimental frontal and transverse plane kinematics for DMD3. A Simulated kinematics when only modeling weakness. B Simulated kinematics when modeling weakness and contractures. C Simulated kinematics when modeling weakness, contractures and midfoot break. D Experimental gait kinematics. E Differences in kinematics between the TD model and the mean weakness model (light green), the TD model and the mean weakness + contractures model (dark green), the TD model and the mean weakness + contractures + midfoot break model (medium green), and the experimental TD and DMD3 data (red). Abd, abduction; Add, adduction; Ex, external; In, internal; TD, typically developing

The main differences in joint moments were absent hip extension moment, premature hip flexion moment, reduced peak hip flexion moment, decreased knee extension and flexion moments, absent dorsiflexion moment, premature plantar flexion moment, and reduced peak plantar flexion moment (Fig. 10). The premature hip flexion moment was not captured in the simulations and the simulation predicted a decrease in hip abduction moment that was not observed experimentally.

Fig. 10

Comparison of simulated and experimental kinetics for DMD3. A Simulated kinetics when only modeling weakness. B Simulated kinetics when modeling weakness and contractures. C Simulated kinetics when modeling weakness, contractures and midfoot break. D Experimental gait kinetics. E Differences in kinetics between the TD model and the mean weakness model (light green), the TD model and the mean weakness + contractures model (dark green), the TD model and the mean weakness + contractures + midfoot break model (medium green), and the experimental TD and DMD3 data (red). Abd, abduction; Add, adduction; DF, dorsiflexion; Ext, extension; Flex, flexion; PF, plantar flexion; TD, typically developing

Contractures and foot deformities contribute to gait impairments in DMD

For DMD1, mainly weakness captured the differences between the mildly affected gait pattern and TD gait (Figs. 2–4). Yet, by adding the contractures to the weakness, additional gait deviations were explained such as reduced dorsiflexion in swing and at initial contact, and the delayed knee flexion and increased peak knee flexion during swing.

For DMD2, weakness already captured most differences between the tiptoeing gait pattern and TD gait (Figs. 5–7). The pattern at the ankle was effectively reflected by only modeling muscle weakness. However, additional gait deviations were captured by modeling contractures in addition to weakness, such as increased hip flexion, knee flexion, and plantar flexion during swing. The loss of knee flexion motion during loading response was not adequately captured by either muscle weakness or the combination of weakness and contractures.

For DMD3, weakness already explained the majority of the differences between the flexion gait pattern and TD gait (Figs. 8–10). In contrast to DMD2, the loss of knee flexion motion during loading response was effectively captured by muscle weakness in DMD3. Some additional gait differences were reflected by adding contractures to weakness, such as increased hip abduction during swing. However, a tiptoeing pattern (i.e., the main characteristic of DMD2) at the ankle with increased knee flexion was predicted by combining weakness and contractures, though this was not observed experimentally. By modeling a midfoot break in addition to the weakness and contractures, this pattern shifted towards increased dorsiflexion at the end of stance and reduced knee flexion during loading response, which better reflected experimental observations.

The effect of progressive impairments

For DMD1, simulations based on the less affected −1SD models, reflecting muscle impairments one SD below the mean, closely resembled simulations based on the mean models, both when modeling weakness in isolation and weakness and contractures (Figs. 2–4). In contrast, simulations based on the more affected +1SD models, reflecting muscle impairments one SD above the mean, clearly differed from simulations based on the mean models and exhibited some gait features that were observed in more severely affected children with DMD, such as tiptoeing at the ankle.

For DMD2, simulations based on −1SD and +1SD models differed from the simulations of the mean models (Figs. 5–7). Simulations based on the −1SD models had less severe deviations, while simulations based on +1SD models had more severe deviations. The tiptoeing gait pattern at the ankle was consistently present across all simulations based on these ±1SD models, both when modeling weakness in isolation, and weakness and contractures.

For DMD3, simulations based on the less affected −1SD models closely resembled simulations based on the mean models (Figs. 8–10). The simulations based on the more affected +1SD weakness model differed slightly more from simulations based on the mean models, with more severe gait deviations such as increased drop foot during swing, increased plantar flexion at initial contact with a reversed first ankle rocker, increased pelvic downwards/upwards range of motion, and increased hip abduction during swing. Interestingly, simulations based on the more affected +1SD weakness and contractures model could not find a feasible solution.

Discussion

Most but not all key features of the three DMD gait patterns were captured in the predictive simulations based on models reflecting DMD-specific muscle weakness and contractures. Overall, the simulated differences between TD and DMD gait patterns more closely matched the experimental differences when stiffness/contractures were modeled alongside muscle weakness. This suggests that not only muscle weakness but also muscle stiffness/contractures contribute to the pathological gait pattern. Interestingly, the tiptoeing gait pattern (DMD2) was observed when only modeling weakness, suggesting that weakness causes tiptoeing even in the absence of contractures. For DMD3, the experimental differences were even better captured when foot deformities (midfoot break) were included in addition to weakness and stiffness/contractures. Our current results highlight the potential of predictive simulations to study the effect of progressive weakness and contractures. Not only could the predictive simulations capture a mildly affected gait pattern and two more severely affected gait patterns based on the severity and location of the muscle impairments, but they also predicted less affected gait features for less severe impairments (−1SD) and more affected gait features (DMD1 and DMD2) or even the inability to walk (DMD3) for more severe impairments (+1SD).

Despite the assumptions and the limited accuracy of the clinical measurements, the models and predictive simulations were able to capture the main features of DMD gait. This suggests that models can be personalized on a population-based level using clinical measurements. This is especially promising given that clinical measurements are sparse (e.g. no assessment of individual muscles), and depend on the cooperation of the boys with DMD and/or the interpretation of the clinician. The instrumented strength assessment is proven reliable [31] but depends on the child’s cooperation and might therefore not reflect the maximal force [50]. The clinical examination, including the range of motion measurements and the clinical stiffness scale, is subjective and depends on the assessor. Additionally, these clinical measurements were expressed as deficits relative to TD children, which could introduce an additional error if these measurements in TD children were also not performed accurately. Personalizing the active muscle force based on the experimentally measured strength outcomes is relatively straightforward. However, several assumptions were made when modeling contractures. In particular, it is unclear how clinical stiffness scores relate to alterations in the muscle’s passive force-length behavior. In a Hill type model, passive muscle force is a non-linear function of muscle length. Here, we chose to model contractures as a shift of the passive force-length relationship, i.e., a decrease in muscle length at which passive force starts to develop. This captures the increased resistance at smaller muscle stretches but it is possible that the tissue stiffness also increases. Future studies should experimentally measure passive torque-angle relationships and use this information to validate and if needed adapt the proposed modeling workflow. For accurate predictive simulations, it is sufficient that we capture the passive torque across postures. This can be obtained with a Hill model but a better understanding of how passive tissue properties change is needed to inform personalization of the passive force-length relationship. On the other hand, more complex muscle models [54,55,56] could give valuable insights into the relationship between changes in muscle microstructure and alterations in the muscle’s overall mechanical behavior [55, 56], and when informed with experimental data on muscle microstructure could outperform Hill-type muscle models.

To further evaluate the robustness of our findings against modeling errors, we evaluated the sensitivity of the differences between simulated TD and DMD kinematics and kinetics to the amount of muscle weakness and contractures separately (Additional files 3–5). The simulations in which weakness and contractures were simultaneously altered by 1SD already give an indication of the sensitivity of our simulations to model parameters. Additionally, we performed three further analyses. In the first, muscle weakness was held constant at its mean value, while contractures were progressively increased by shifting the passive force-length curve to 0.9, 0.8, and 0.7 of the normalized fiber length (Additional file 3). The results indicate that shifting the passive curve by 0.1 had only a small effect on the differences between simulated TD and DMD kinematics and kinetics. Given the non-linear relationship between model parameters and gait kinematics and kinetics, we varied contractures around their mean value by ±1 SD while holding muscle weakness constant at its mean value (see Additional file 4). The results showed minimal effects of varying the shift of the passive force-length curve by 0.04 to 0.1 on the differences between simulated TD and DMD kinematics and kinetics. Together, these analyses suggest that small deviations (i.e., a shift about 0.1) in contracture estimation do not substantially affect the outcomes. Only increasing contractures by +1SD for DMD1 had a larger impact and resulted in larger differences between the simulated and experimental differences as the simulated gait pattern began to resemble the gait pattern of DMD2. This suggests that the proposed estimates of the shift are more realistic than larger shifts introducing larger contractures. Finally, in the third analysis, we varied muscle weakness around its mean value by ±1SD while holding contractures constant at their mean value (Additional file 5). Here, the impact of altering muscle strength by 7% to 35% was greater than that of varying the shift of the passive force-length curve by 0.04 to 0.1 (Additional file 4), indicating that the simulations are more sensitive to muscle weakness than contractures. However, we also believe that our muscle weakness estimates are more accurate due to the use of quantitative data.

The proposed approach can also be used to create personalized models for individual patients. However, further validation is required as measurement errors (assessments of muscle weakness and contractures) might be larger for individual patients than on a population level where we could use a mean value across many patients. Whereas asymmetry is typically limited, it has been observed in some patients. Patient-specific simulations could further elucidate the effect of asymmetric weakness or contractures on the gait pattern.

Simulations based on a cost function that was designed for simulations of unimpaired gait captured key features of DMD gait, suggesting that motor control is largely unaffected by DMD. The cost function used here was selected because it best captured unimpaired gait at self-selected speed [20]. Our simulation results confirm the underlying assumption that motor control of children with DMD is unimpaired. However, it is possible that we need to consider additional factors in the cost function, such as avoiding eccentric contractions or pain, to predict fine details of DMD gait. For example, in DMD2, we were unable to replicate the loss of knee flexion motion during loading response because this group still had sufficient strength in the knee extensors, unlike DMD3. Therefore, it is possible that children with DMD adopt this pattern to avoid eccentric contraction or pain [57] in the knee extensors, despite having sufficient strength. Further consideration of these factors in the cost function may enhance the accuracy and applicability of the simulations.

The difference in mass and mass distribution between the TD group and DMD groups could have influenced the simulated gait patterns [21]. We aimed to reflect the altered anthropometry of boys with DMD by distributing the additional mass compared to TD with 2/3 over the trunk and pelvis segments (reflecting that additional fat tissue is predominantly distributed around the abdomen and buttocks) and with 1/3 over the lower limb segments (reflecting secondary fat accumulation due to increased intramuscular and subcutaneous fat in the legs). The mass and mass distribution were therefore not constant across the different models. Falisse et al. [21] investigated the influence of mass distribution on simulations and found that a lighter torso and heavier legs produced slightly larger knee flexion angles, knee extension moments, and knee extensor activations during stance. Consequently, the different mass and mass distribution of the models could have influenced the differences between simulated gait patterns. We adjusted only the mass and inertia. However, the accumulation of fat tissue around the abdomen shifts the trunk’s center of mass forward. Not accounting for this shift in the center of mass may have contributed to the predictive simulations’ limited ability to capture the increase in posterior trunk leaning and premature hip flexion moment. Future studies should investigate the impact of mass, mass distribution, and trunk center of mass on the simulated gait patterns.

Muscle weakness but also contractures contribute to DMD gait. Previously, muscle weakness was considered the most important factor influencing pathological gait in DMD [1, 5], whereas the effect of contractures, in addition to muscle weakness, was poorly understood. Our predictive simulations indicate that muscle weakness is the main contributor to DMD gait deviations but some gait features were only predicted by modeling the interaction between muscle weakness and contractures. Therefore, the role of contractures cannot be underestimated, as they also contribute to the pathological gait patterns.

Altered foot properties might have an important contribution to the flexion gait pattern (DMD3). When the midfoot break was not modeled, a tiptoeing gait pattern was predicted, which did not match the experimental data. Videos of some boys with DMD showed a premature heel rise with abnormally increased motion through the midfoot (i.e., a midfoot break), leading to a dorsiflexion overestimation in the experimental data due to the applied simplistic marker foot model. The midfoot break might be caused by more flexible feet with larger arch drops and lower arch rigidity, secondary from obesity, combined with the plantar flexion contracture restricting the ankle dorsiflexion [17]. When a midfoot break was included in the model, there was still tiptoeing gait at initial contact, but dorsiflexion increased again at terminal stance, and knee flexion decreased during loading response, better corresponding to the experimental gait. This highlights the importance of using detailed, complex 3D foot models and modeling foot deformities to accurately capture pathological gait patterns.

Predictive simulations have the potential to study the effect of progressive impairments on the DMD gait pathology. Our exploration based on models that were less or more impaired yielded plausible results and led to testable hypotheses. For example, when increasing weakness but not stiffness and weakness in DMD3 a gait pattern could still be found thereby suggesting that increases in stiffness are more threatening for gait function than increases in weakness. Future experimental research should test if increases in stiffness indeed play a crucial role in the loss of ambulation. Furthermore, future work is needed to examine the individual effects of the progressive weakness and contractures in different muscle groups to gain insights into potential supportive interventions, such as the use of exoskeletons.

The differences between experimental and simulated gait patterns arise from two factors: the use of different models and the inaccuracies of the predictive simulations. As demonstrated in a previous study [58], the Plug-in Gait and OpenSim models are not fully comparable with kinematic differences up to 13°. The majority of these discrepancies are attributed to differences in the anatomical models, including different anatomical segment frames and joint constraints, rather than differences in computational methods (direct versus inverse kinematics). The largest difference is observed in pelvis anterior/posterior tilt angles. This can be attributed to different definitions of neutral pelvis tilt. In the OpenSim model neutral pelvis tilt is aligned with the anatomical position, whereas in the Plug-in Gait model neutral pelvis tilt is based on the pelvic markers, resulting in a 13° offset between the two reference systems. We did account for this offset between Nexus and OpenSim, but our marker set did not allow us to perform inverse kinematics in OpenSim for the experimental data based on the model we used for simulations. In particular, we did not have enough markers on the feet to estimate detailed foot and ankle kinematics. The use of two different models primarily affects kinematics. We, therefore, focused on comparing the differences between DMD and TD kinematics but this comparison may still be slightly influenced by the use of different kinematic models. However, not all discrepancies between experimental and simulated gait can be attributed to differences in modeling joint kinematics. Predictive simulations do not fully capture experimental kinematics and kinetics. For example, the simulations underestimate peak knee flexion during swing. This might be due to errors in modeling quadriceps properties in the generic model leading to unrealistically large passive forces when the knee is deeply flexed. Additionally, the simulations failed to capture the large increases in anterior pelvic tilt. We used an existing model of musculoskeletal geometry that is most frequently used (OpenSim’s gait2392 model [51]). The trunk model is simplified with only one joint, and the muscle anatomy at the hip/pelvis is not very detailed. For instance, the psoas muscle in the musculoskeletal model [51] inserts at the pelvis instead of the lumbar vertebrae. This might explain the underestimated increase in anterior pelvic tilt and posterior trunk leaning, and the absent premature hip flexion moment in the simulations. Furthermore, the segment lengths of the generic adult model were scaled to those of children but the differences in bone shapes were not yet considered. For example, the femoral anteversion angle decreases from 40° to 15°, and the neck-shaft angle from 140° to 125°, from birth to skeletal maturity [52, 53], which affects the moment arms of the muscles. Nevertheless, this model provides the most realistic predictive simulations of gait [25] and we deemed it sufficiently realistic to study the effect of weakness and contractures on the gait pattern. However, refining the anatomical details and incorporating more accurate passive force-length relationship measurements might enhance the predictive accuracy of these simulations.

Conclusion

Predictive simulations based on personalized models capture key features of DMD gait patterns. While muscle weakness was primarily responsible for gait deviations, muscle stiffness/contractures and foot deformities further contributed to gait deviations. These results highlight the potential of predictive simulations to improve our understanding of the causal relationships between progressive impairments and pathological gait patterns in boys with DMD. Increased understanding of the differential effect of weakness and stiffness might guide treatment selection.