Transfemoral limb loss modestly increases the metabolic cost of optimal control simulations of walking

Virtual subjects

The computational speed of modern optimal control methods allows for simulations of multiple subjects, producing datasets that can undergo formal statistical analysis and hypothesis testing similar to studies on human participants. As a computer simulation study, the subjects in the present study were instances of the model with different model parameter values. The pre-limb loss model had a standing height of 1.70 m and body mass of 75.3 kg. Post-limb loss, the total body mass was 70.8 kg, with 67.3 kg of biological mass and 3.5 kg of prosthesis mass. This instance of the model was defined as “Subject 00”. Additional subjects were generated using the OpenSim “Scale Model” tool to scale the Subject 00 pre-limb loss model to the height and mass of each of the 14 able-bodied control participants from Russell Esposito, Ràbago & Wilken (2018), all of whom were United States military service members (mean age 26 years, range 22–44 years). The dimensions of all axes of all body segments in the model were scaled linearly by the ratio of the target height to the Subject 00 height. The mass of each body segment was scaled linearly by the ratio of the target total body mass to the Subject 00 total body mass. The segment dimension scaling also adjusted the muscle optimal fiber lengths and tendon slack lengths. Muscle maximum isometric forces were then scaled linearly by the ratio of the target whole-body mass to the Subject 00 whole-body mass.

Including Subject 00 and the 14 subjects scaled to the heights and masses from Russell Esposito, Ràbago & Wilken (2018), there was a total of 15 subjects in the study. We emphasize that these subjects were not intended to represent any specific individuals, and did not represent an exhaustive sensitivity analysis of the model’s parameter values. The element of subjects in the present study is better viewed as a limited sensitivity analysis, providing some confidence that any reported effects were due to the independent variable (limb loss) and generalize beyond the original model’s mass and height. The sample of 15 subjects had a minimum detectable effect size of 0.68 in the paired directional t-test of the change in metabolic cost pre- vs. post-limb loss, using standard long-term error rates of 5% for false-positives and 20% for false-negatives. In the present data, this effect size equated to a minimum detectable change in metabolic cost of about 1%.

Human participant

As a point of comparison to the simulations, a single high-functioning human participant with unilateral transfemoral limb loss completed an instrumented gait analysis of walking at their preferred “normal and comfortable” walking speed of 1.30 m/s. Ground reaction forces were measured at 1,000 Hz using piezeoelectric force platforms (Kistler). Positions of retroreflective markers on the body were measured at 200 Hz using 13 optical motion capture cameras (Vicon) and used to calculate joint and segment angles with standard six-degrees-of-freedom methods (Miller et al., 2014). The participant was fit with a chest worn heart-rate monitor (Garmin) and a wearable metabolic system which included a face mask and backpack-worn wireless transceiver (Cosmed K5). Metabolic rate was calculated from the measured steady-state pulmonary gas rates using Brockway’s equation (Brockway, 1987) with the urinary nitrogen term truncated. Participation was approved by the University of Maryland institutional review board, protocol #1826662. The participant gave written informed consent The participant was male, 26 years old, height 1.90 m, mass 96.3 kg, and had undergone limb loss due to traumatic injury. The participant was Medicare K-level 4, used a microprocessor-controlled prosthetic knee, and frequently participated in activities like weightlifting, running, and cycling.

Overview

The models were used to simulate periodic strides of walking at average speed 1.45 m/s and stride duration 1.0 s, values reported by Miller et al. (2014) for healthy adults walking at a subjective “normal and comfortable” pace. The simulations were posed as optimal control problems, finding the time-varying muscle excitations and associated model states that minimized the weighted sum of (i) deviations from average able-bodied gait mechanics and (ii) the gross metabolic cost, in the form of the cost function J: (3)${\displaystyle J={\int }_0^T\left[\frac{1}{NT}\sum _{i=1}^N{\left(\frac{x_i\left(t\right)-{\mu }_i\left(t\right)}{{\sigma }_i}\right)}^2+\frac{w_1}{m\Delta x}\left({\dot{E}}_b+{\dot{E}}_m\left(t\right)\right)\right]dt}$

The first term on the righthand side of Eq. (3) is the “gait deviations tracking error”: $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 variable, averaged over the gait cycle. The N = 37 variables included in the tracking error term were the six translations and rotations of the pelvis, the 25 joint angles, and the three components of the ground reaction force under each foot. The square root of this term gives the average tracking error of the simulation, averaged over all timesteps of all variables, in multiples of the standard deviation. For example, a tracking error of 0.5 means that the simulation’s values for these 37 variables deviated from the average gait mechanics, defined by ${\mu }_i\left(t\right)$, by 0.5 standard deviations on average. For the motions of toe flexion, subtalar inversion, hip internal rotation, and the non-sagittal shoulder motions, we had low confidence in the accuracy of the experimental data, so we increased the value of σ_i to 3σ_i for these variables. This change in weighting allowed for better tracking of variables that can be measured with greater confidence in their accuracy, such as the sagittal joint motions and the ground reaction forces.

The second term on the righthand side of Eq. (3) is the metabolic cost: ${\dot{E}}_b$ is the non-muscular basal metabolic rate, ${\dot{E}}_m\left(t\right)$ is the sum of the muscle metabolic rates, m is the total body mass, and Δx is the horizontal displacement of the center of mass. The weighting constant w₁ defines the relative emphasis the optimizer places on reducing the tracking error vs. the metabolic cost, which are generally in opposition (smaller tracking errors have greater metabolic costs). Values of w₁ = 0.05 to 0.15 produced realistic metabolic costs of roughly 3–4 J/m/kg during model development with the Subject 00 pre-limb loss model, and a middle value of w₁ = 0.10 was chosen. The same value (w₁ = 0.10) was thereafter used in all simulations of all subjects, pre- and post-limb loss. We emphasize that the goal of these simulations was not to produce post-limb loss simulations that necessarily resemble the mechanics, energetics, and control of actual humans walking with a transfemoral prosthesis. Rather, the goal was to define specific clinically relevant objectives (walk with minimum deviations from a target gait and minimum metabolic cost) whose presence and relative weighting were identical for the pre- and post-limb loss simulations, such that differences between simulations could be confidently attributed to independent variable (limb loss) without a confounding effect of differences in the goals of the movement.

The optimal control problems were converted to nonlinear programming problems and solved using the direct collocation method in OpenSim Moco software (Dembia et al., 2020). The model’s state and control variables were discretized with a Hermite-Simpson transcription scheme, with 101 nodal values per variable spaced evenly over the stride duration, i.e., one node every 1% of the gait cycle. Simulation results in model development were invariant to finer grids tested up to 201 nodes per stride, at which point the computational speed of the simulations became impractically slow. A full stride was simulated because the assumption of bilateral symmetry involved in simulating a single step (Anderson & Pandy, 2001) is not valid for asymmetric limb loss locomotion. The nodal state and control variable values were optimized to minimize the value of 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. The IPOPT solver within Moco was used for optimization (Wächter & Biegler, 2006), with constraint tolerance 1•10⁻⁴ and solver tolerance 1•10⁻³.

Validation simulations

Directly validating the model’s accuracy for predicting changes in metabolic cost after limb loss would require impractical longitudinal data on human metabolic cost pre- and post-limb loss. We instead gauged the model’s validity for predicting changes in metabolic cost with changes in lower limb mechanics that conceptually resemble limb loss. Previously, we demonstrated the model’s validity for predicting changes in metabolic cost with changes in foot mass and ankle stiffness (Miller & Russell Esposito, 2021). Here we added evaluations of the predicted changes in metabolic cost with changes in thigh mass and knee stiffness.

Browning et al. (2007) measured net metabolic rate of walking at 1.25 m/s with either +4 or +8 kg of mass added to each thigh (total of +8 or +16 kg). The data were well fit by a simple linear regression where the net metabolic rate increased by an average of 0.075 W/kg per kilogram of added mass. We added the same masses to the thigh segments of the Subject 00 pre-limb loss model and performed simulations of walking at 1.25 m/s using the same cost function and simulation approach described earlier. The model’s net metabolic rate in these simulations increased by 0.077 W/kg per kilogram of added mass.

Lewek, Osborn & Wutzke (2012) used a commercial knee brace that unilaterally limited the knee to about 30°  of flexion, and reported an average increase in net metabolic power of 17% when walking at 1.25 m/s compared to walking with unrestricted knee motion. We simulated knee bracing by adding a “CoordinateLimitForce” actuator to the right knee of the Subject 00 pre-limb loss model that made it essentially impossible for the model to flex its knee beyond 30°, and performed simulations of walking at 1.25 m/s using the same cost function and simulation approach described earlier. We were unable to achieve the magnitudes of metabolic rates reported by Lewek, Osborn & Wutzke (2012). For example, their baseline “control” condition with the knee brace unlocked had an average metabolic rate of 3.8 W/kg, which is unusually high for walking at 1.25 m/s. The model’s net metabolic rate at 1.25 m/s without the CoordinateLimitForce was 2.60 W/kg. The difference here was likely due to unmodeled aspects of the brace that affected the symmetry of the human participant’s walking gaits even with the brace unlocked. However, the model’s increase in metabolic rate with the CoordinateLimitForce on the right knee was +0.74 W/kg, which compares well to the increase in Lewek, Osborn & Wutzke (2012) of 0.66 ± 0.60 W/kg with the brace locked vs. unlocked.

Changing lower limb mass and joint stiffness of an intact biological limb is of course not equivalent to walking with a prosthetic lower limb. However, reasonably accurate model-based predictions of these changes provide confidence in using the model to predict changes in metabolic cost for conditions that cannot easily be validated directly with experiments.

The validity of the pre-limb loss model for predicting realistic muscle excitations when realistic motion and ground reaction forces in walking are imposed on it has been reported elsewhere (Lai, Arnold & Wakeling, 2017). As an additional gauge of the pre-limb loss model’s validity, Fig. 2 shows the on/off timing of excitations for 28 muscles in the baseline model’s pre-limb loss walking simulations, compared to normative timing measured by indwelling electromyograms from the same 28 muscles for walking in healthy adults (Sutherland, 2001). The post-limb loss excitation timings are also shown. Muscles in the model were deemed “on” when the excitation magnitude was above 3% on the range of 0–100%, and “off” otherwise. The main deviations in muscle on/off timing between the pre-limb loss simulation and the referenced electromyograms were the overall timing of the peroneal muscles and the lack of co-activation of all three biarticular hamstrings muscles, but otherwise the timing of the major muscles was reasonably similar to the referenced electromyogram timing.

It is difficult to make comparisons to electromyogram data from human participants with transfemoral limb loss since those datasets are sparse and rather variable across participants and across studies (e.g., Jaegers, Arendzen & de Jongh, 1996; Wetink et al., 2013), but the changes in excitation timing with limb loss were at least theoretically sensible, e.g., prolonged excitation of the prosthetic-side hip flexors. It is important to re-emphasize that the goal of these post-limb loss simulations was to change as little as possible about the model and the simulation workflow, other than replacing one biological limb with a prosthetic limb, such that the effect of the prosthetic limb could be isolated. We did not attempt to create models that include all typical differences between individuals with and without prosthetic limbs, nor to generate simulations that necessarily match the mechanics, energetics, and control of real humans who walk with a prosthetic limb. Other approaches are better suited to address those goals.

Pre-limb loss simulations

Simulations with the pre-limb loss models were performed using the direct collocation method. The same tracking targets ${\mu }_i\left(t\right)$ and σ_i were used in the cost function (Eq. (3)) for all subjects. Each simulation was performed with three different initial guesses using simulation results from a previous similar study (Miller & Russell Esposito, 2021) and the result with the lowest cost function score was retained for analysis. When none of these initial guesses produced a result with a realistic metabolic cost, additional simulations with different initial guesses, specifically using the results of simulations from other “virtual subjects” described earlier as initial guesses, were performed until a result with a metabolic cost within two standard deviations of the mean from Das Gupta, Bobbert & Kistemaker (2019) was found.

Post-limb loss simulations

Simulations with the post-limb loss models were performed identically to the pre-limb loss simulations, with the exception that the standard deviations for the prosthetic limb’s knee, ankle, subtalar, and toe joints in Eq. (3) were increased to 10 σ_i, effectively reducing the weight on accurate tracking of these variables. The weight on these variables was reduced because the tracking targets ${\mu }_i\left(t\right)$ were biological joint motions and it was expected that the model could not track these motions well with prosthetic joints. Motions of these joints in unpowered prostheses also typically do not closely resemble biological joint motions (e.g., Kobayashi et al., 2020). To approximate the patient-specific selection and adjustment process typically done in prosthesis prescription, the k, b, and θ₀ prosthetic joint parameters in Eq. (1) were optimized along with the model states and controls in Moco for each subject’s post-limb loss simulation to minimize the cost function. Stiffness k was bounded on [1, 100] Nm/rad for the knee and [10, 1000] Nm/rad for the ankle. Damping b was bounded on [1, 10] Nm s/rad for both joints. Neutral joint position θ₀ was bounded on [−10, 10] degrees for both joints. The prosthetic subtalar and toe joints were set to k = 400 Nm/rad, b = 1.0 Nm s/rad, and θ₀ = 0°.

Post-limb loss simulations were performed twice for each subject, once with the passive prosthesis in Eq. (1) and once with the prosthesis that included the ActivationCoordinateActuator at the knee, Eq. (2). As noted earlier, the purpose of this actuator was to allow small deviations from the torque profile of the fully passive prosthetic knee, not to simulate a fully powered motorized prosthetic knee. This purpose was achieved by minimizing the control signal of the actuator in the cost function: (4)${\displaystyle J={\int }_0^T\left[\frac{1}{NT}\sum _{i=1}^N{\left(\frac{x_i\left(t\right)-{\mu }_i\left(t\right)}{{\sigma }_i}\right)}^2+\frac{w_1}{m\Delta x}\left({\dot{E}}_b+{\dot{E}}_m\left(t\right)\right)+w_2{u}^2\left(t\right)\right]dt}$ where $u\left(t\right)\in \left[-1,1\right]$ is the prosthetic knee ActivationCoordinateActuator’s control signal. The weighting on this signal was given a high value of w₂ = 100. Given the value of w₁ = 0.1 in combination with w₂ = 100, the cost of using the ActivationCoordinateActuator to produce the torque from the pre-limb loss simulation that tracked the target biological knee angle well equated to a very high increase in metabolic cost of +47 J/m/kg. The cost of using the powered prosthesis to alter the knee torque produced by the passive component of the prosthesis by an average of 1.0 Nm across the gait cycle equated to an increase in metabolic cost of +0.1 J/m/kg. The model with the ActivationCoordinateActuator (Eq. 4) will hereafter be referred to the “non-passive” prosthesis model. For initial guesses in the post-limb loss simulations, the pre-limb loss simulation result of the same subject was used first, then two other pre-limb loss simulations from subjects with heights closest to the present subject. If none of these guesses produced a result with a realistic metabolic cost, other initial guesses from other pre- and post-limb loss subjects were used until a result with a metabolic cost within two standard deviations of the mean from Das Gupta, Bobbert & Kistemaker (2019) was found.

Outcome variables and statistics

The main outcome variable of each simulation was the gross metabolic cost, defined as the metabolic energy expended per unit distance traveled by the model’s center of mass, divided by the body mass (units = J/m/kg). Paired Student’s t-tests were used to test the hypothesis that metabolic cost is greater post-limb loss vs. pre-limb loss, for both the passive and non-passive prosthesis models. The threshold for significance was p = 0.05, with Bonferroni adjustments for multiple comparisons.

Some studies on the limb loss population scale metabolic cost by the “total mass”, which includes the mass of the prosthesis (e.g., Jarvis et al., 2017). Others scale metabolic cost by the “biological mass”, excluding the mass of the prosthesis (e.g., Russell Esposito, Ràbago & Wilken, 2018). Because the prosthesis mass is a substantial fraction of the total body mass for individuals with transfemoral limb loss, the statistical analysis above was performed twice, once with the post-limb loss metabolic costs scaled only by biological mass, and a second time with the post-limb loss metabolic costs scaled by total mass. With two statistical comparisons to the pre-limb loss simulations with the different prosthesis models, and two comparisons with different mass scaling, the threshold for rejecting the null hypothesis was p = 0.05/4 = 0.0125.