Transtibial limb loss does not increase metabolic cost in three-dimensional computer simulations of human walking

Subjects

As a computer simulation study, the “subjects” in the present study were not specific human participants but were instances of the model described in ‘Pre-Limb Loss Model’–‘Post-Limb Loss Model’ with different parameter values. “Subject #1” used the default values of the pre-limb loss model’s parameters. To generate additional pre-limb loss subjects, a scaling factor between 0.75–1.25 was drawn randomly from a standard normal distribution and scaled all the model’s maximum isometric muscle forces to adjust the model’s overall strength, e.g., if the scaling factor was 1.12, every maximum isometric force was multiplied by 1.12. Smaller random scaling factors of 0.95–1.05 that varied for each individual parameter were then applied to the scaled maximum isometric muscle forces, tendon slack lengths, and body segment masses. The body segment masses were constrained to still sum to the model’s original mass and all parameter values were still required to be bilaterally symmetric, as in the original model. The ranges over which these parameters were set to produce a range of pre-limb loss metabolic costs of roughly 3.4–3.8 J/m/kg. These models were then converted to post-limb loss models (‘Post-Limb Loss Model’). For each subject’s post-limb loss model, a random fraction of 20–80% of the prosthetic limb’s shank mass was assumed to be biological mass to represent a range of residual limb sizes.

The muscle maximum isometric forces and tendon slack lengths were perturbed because Hill model output is very sensitive to these parameter values (e.g., Scovil & Ronsky, 2006). The segment mass distributions were perturbed because confidence in their values for a human subjects experiment is typically low. Other unmentioned parameters that could feasibly affect metabolic cost were not perturbed. The element of “subjects” in this study should essentially be viewed as a limited sensitivity analysis, hopefully providing some confidence that the results are due to the independent variable (limb loss) and not to the values assumed for parameters that are difficult to assign on a subject-specific basis with confidence. Similarly, the “subjects” here are perhaps better viewed as different parameter sets that may feasibly be assigned to an individual, or as different possible baseline fitness levels of an individual, rather than distinct individuals.

The minimum detectable effect of interest was a 2% change in metabolic cost, which is the minimum detectable change reported in testing of a modern portable pulmonary gas exchange unit due to measurement error alone, in the absence of any human biological variability between tests (Guidetti et al., 2018). The motivation for this choice of 2% is that smaller changes would be difficult to measure on human participants even if they were real changes. Simulation results from 20 subjects suggested at least 26 subjects were needed to detect the 2% effect as significant with error rates of 5% for both Type I and Type II errors. A total of 36 subjects were therefore used. In the present data, 2% equaled an effect size of 0.71, or ±0.068 J/m/kg in raw units of metabolic cost.

Overview

The models described above were used to simulate periodic strides of walking with stride duration T = 1.0 s and speed 1.45 m/s, the averages from healthy young adults measured by Miller et al. (2014). The simulations were posed as optimal control problems, finding the time-varying muscle excitations and associated model states that minimized a cost function J that was the weighted sum of gait deviations plus muscular effort: (3)${\displaystyle J=\frac{1}{T}{\int }_0^T\left[\frac{1}{N_1}\sum _{i=1}^{N_1}{\left(\frac{x_i\left(t\right)-{\mu }_i\left(t\right)}{{\sigma }_i}\right)}^2+\frac{w}{N_2\Delta x}\sum _{i=1}^{N_2}{m}_i{u}_i^2\left(t\right)\right]dt.}$

The first term on the right-hand side of Eq. (3) is the “gait deviations”: $x_i\left(t\right)$ is the value of model variable i at time t, ${\mu }_i\left(t\right)$ is the mean of the analogous variable from the experimental gait analysis data of Miller et al. (2014), and σ_i is the standard deviation between-subjects of the experimental data, averaged over the stride cycle. The N₁ = 37 variables included in this term were the six translations and rotations of the pelvis, the 25 joint angles, and the 3-D components of each ground reaction force (GRF). The square root of this term after dividing the integral by T gives the average tracking error, in multiples of SD. For example, a tracking error of 1.0 indicates the model’s gait mechanics were 1.0 SD away from the experimental means when averaged over all timesteps and all tracking targets. Hip internal rotation, subtalar motion, and non-sagittal shoulder motions did not track experimental data due to low confidence in accurately measuring these motions with skin marker-based motion capture. These motions tracked target angles of 0° to minimize their ranges of motion, which are typically small in walking. Assigning no tracking targets for these degrees of freedom tended to produce unusual movements to track other targets unnecessarily well, e.g., the GRF. The hands were constrained to not pass through the pelvis. Lumbar joint motions were also not tracked. Instead, deviations of the trunk angle from the global axes were minimized. The standard deviations of the experimental data were still used for σ_i as they were presumed to reflect the variability of the underlying motions.

The second term on the right-hand of Eq. (3) is the “muscular effort per unit distance”: $u_i\left(t\right)$ ranges from 0-1 and is the excitation of actuator i at time t, m_i is the mass of muscle i, and Δx is the horizontal displacement of the model during the stride. N₂ = 92 is the number of actuators in the model: 84 muscles and eight torque generators. This “muscular effort” term is a surrogate of metabolic cost and was used in the cost function instead of metabolic cost itself because (i) this version of the Moco software could not include metabolic cost in the cost function, and (ii) it is unclear if the mechanics of normal walking are consistent with a control policy that minimizes metabolic cost itself (Bertram, 2005; Ackermann & Van den Bogert, 2010). Smaller gait deviations will generally require greater muscular effort, assuming the model’s minimum-effort gait is not identical to the tracking targets. The weighting constant w determines the emphasis of the optimization on minimizing effort vs. deviations. Values of w = 50–200$∕\bar{m}$, where $\bar{m}$ is the average muscle mass, produced realistic metabolic costs in the range of 3.0–3.8 J/m/kg, typical of healthy young adults at comfortable walking speeds (Das Gupta, Bobbert & Kistemaker, 2019), with average tracking errors typically under 1.0 SD. A value of w = 100$∕\bar{m}$ was therefore used in all simulations. For the upper limb torque generators, m_i was set equal to $\bar{m}$. The same value of w = 100$∕\bar{m}$ was used in all simulations for all subjects, both pre- and post-limb loss; the value of w was not further adjusted to necessarily achieve a realistic metabolic cost in every simulation, nor to minimize the change in metabolic cost with limb loss, nor any particular pre-determined result. Keeping the same value of w in both the pre- and post-limb loss simulations aligns with the goal of this study in isolating the effect of limb loss.

The optimal control problems were converted to nonlinear programming problems and solved using a direct collocation method in Moco 0.4.0 software (Dembia et al., 2020). Briefly, the model’s state and control variables were discretized on a temporal grid of 101 nodal values spaced evenly over the stride duration. Simulation results in model development were invariant to finer grids tested up to 401 nodes per stride. These nodal values were optimized to minimize Eq. (3), subject to constraints of the skeletal equations of motion, muscle activation and contractile dynamics, and task constraints of kinematic periodicity and average speed. This method solves the same optimal control problem as the traditional forward dynamics approach in biomechanics (e.g., Neptune, Kautz & Zajac, 2001) but avoids integrating the model’s state equations forward in time, which allows for high-resolution muscle excitations and relatively fast computational speeds (Van den Bogert, Blana & Heinrich, 2011). The constraint tolerance for convergence of the optimizations was 10⁻⁴.

Validation simulations

Directly validating the model’s ability to predict changes in metabolic cost with limb loss is not possible since there are no longitudinal data on changes in metabolic cost pre- vs. post-limb loss. To gauge the model’s validity for predicting changes in metabolic cost, we compared its predictions to literature reports on how metabolic cost changes when altering lower leg mechanics. Schertzer & Riemer (2014) added 0.5, 1.0, and 2.0 kg mass to each foot (weights around the ankles) and reported average respective increases in metabolic cost of +2.9, 6.7, and 16.0%. We added the same masses to our pre-limb loss model and simulation approach above, using the subject whose metabolic cost was closest to Schertzer & Riemer’s (2014) mean cost with no added mass. The masses were added to the model by increasing the mass of the calcaneus segment. The cost function and simulation procedure was the same as described in ‘Overview’. The simulated metabolic costs at all three added mass levels (+3.5, 6.5, 14.7%) were within one standard deviation of the mean from the experimental data (Fig. 3A). We also evaluated the model’s prediction of metabolic cost with ankle bracing. Wutzke, Sawicki & Lewek (2012) used a unilateral ankle brace that reduced ankle range of motion by 13° and reported an average increase in metabolic cost of 5.8%, when assuming the same resting metabolic rate as the model. We simulated an ankle brace by adding a linear torsional spring to the model’s right ankle and generating a walking simulation, again using the cost function and simulation procedure described in ‘Overview’. With a spring stiffness of 1,000 Nm/rad, the ankle range of motion was reduced by 16°, and the metabolic cost increased by 6.2% (Fig. 3B).

Changing leg mass or ankle stiffness is, of course, not the same as walking with a prosthesis, but reasonably accurate responses of the model’s metabolic cost to these conditions gives some confidence in its ability to accurately predict small changes in metabolic cost for conditions that cannot easily be tested experimentally, such as pre- vs. post-limb loss.

Pre-limb loss simulations

Simulations with each subject’s pre-limb loss model described in ‘Pre-Limb Loss Model’ were performed using the approach described in ‘Overview’. The same tracking targets ${\mu }_i\left(t\right)$ were used in the cost function (Eq. (3)) for all subjects. Each simulation was performed with three different initial guesses generated in model development and the result with the lowest cost function score was retained. The three “initial guess” simulations were generated using the generic pre-limb loss model with different values of w = 10$∕\bar{m}$, 100$∕\bar{m}$, and 1000$∕\bar{m}$ specified in Eq. (3). The initial guess for these simulations was the means of the tracking targets for the kinematic states and zeroes for the muscle states and controls. These simulations were first performed on a coarse collocation grid (11 nodes/stride) which helped with convergence when using a bad initial guess. The coarse result was then interpolated to generate an initial guess on a finer grid, and this process was repeated until a converged solution was achieved on the final 101-node grid.

For readers unfamiliar with the terms and methods of optimal control: (i) an “initial guess” provides a starting point for the optimization routine and contains all the variables to be optimized when generating a simulation, in this case the discretized values of the model’s time-varying state and control variables; (ii) multiple initial guesses are used and the result with the lowest cost function score retained because the goal of the optimization was to minimize the cost function in an absolute sense. For example, a pre-limb loss simulation that converged on a local minimum with a high metabolic cost could bias the post-limb loss simulation towards having a comparatively low metabolic cost.

Post-limb loss simulations

Simulations with the post-limb loss models were performed identically to the pre-limb loss simulations, with the exceptions that (i) the prosthetic limb’s ankle, subtalar, and toe joints were excluded from the tracking term in the cost function (N₁ = 34 in Eq. (3)) because the tracking targets were able-bodied data, and (ii) the 13 muscles of the right ankle were removed (N₂ = 79 in Eq. (3)). To approximate the patient-specific prosthesis tuning process typically done in clinical practice, the stiffness and damping parameters of the three prosthetic joints (Eq. (1)) were optimized in Moco for each subject’s post-limb loss simulation, along with the muscle excitation controls and model states, to minimize the cost function. Stiffness was bounded on 100–1,000 Nm/rad and damping on 0.01–10 Nm/(rad/s).

For all post-limb loss simulations, the same subject’s pre-limb loss simulation result was used as the initial guess. The cost function in the post-limb loss simulations used the values of ${\mu }_i\left(t\right)$ from the able-bodied subjects of Miller et al. (2014) as tracking targets. The post-limb loss model was therefore not explicitly trying to track gait data from individuals with limb loss. The goal of these post-limb loss simulations was to walk with minimal deviations from a typical able-bodied gait with low muscular effort, a common goal of rehabilitation after limb loss.

Outcome variable and statistics

The outcome variable from each simulation was the gross metabolic cost, expressed as the metabolic energy consumed per unit distance traveled, divided by the biological body mass. To be clear, the pre-limb loss metabolic cost was scaled by the pre-limb loss biological mass, and the post-limb loss metabolic cost was scaled by the post-limb loss biological mass. The biological mass was 75.4 kg pre-limb loss and 72.8 ± 0.3 kg post-limb loss. Scaling by biological body mass (excluding the prosthesis mass in the post-limb loss case) was done to avoid favoring the hypothesis of similar metabolic costs pre- vs. post-limb loss, which would have been favored by scaling the post-limb loss metabolic costs by the larger total body mass that includes the prosthesis.

The effect of limb loss on metabolic cost was tested first with a standard Welch’s t-test to determine if the null hypothesis of no change could be rejected. Equivalence of the pre- vs. post-limb loss metabolic costs was then tested using two one-sided tests (Lakens, Schell & Isager, 2018): ${\displaystyle d\le -d_{min}d\ge d_{min}}$

where d is the effect size in the data for the post-limb loss change in metabolic cost, and d_min is the minimum effect size of interest, defined as a 2% change earlier. When both of these hypotheses can be rejected, the change in metabolic cost lies statistically within the equivalence bounds of −d_min ≤ d ≤ d_min, suggesting the pre- vs. post-limb loss metabolic costs are practically equivalent.