Transtibial limb loss does not increase metabolic cost in threedimensional computer simulations of human walking
 Published
 Accepted
 Received
 Academic Editor
 Antonie van den Bogert
 Subject Areas
 Bioengineering, Kinesiology, Orthopedics
 Keywords
 Limb loss, Prosthesis, Military, Muscle strength, Metabolic cost, Computer simulation, Direct collocation, Optimal control
 Copyright
 © 2021 Miller and Russell Esposito
 Licence
 This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, reproduction and adaptation in any medium and for any purpose provided that it is properly attributed. For attribution, the original author(s), title, publication source (PeerJ) and either DOI or URL of the article must be cited.
 Cite this article
 2021. Transtibial limb loss does not increase metabolic cost in threedimensional computer simulations of human walking. PeerJ 9:e11960 https://doi.org/10.7717/peerj.11960
Abstract
Loss of a lower limb below the knee, i.e., transtibial limb loss, and subsequently walking with a prosthesis, is generally thought to increase the metabolic cost of walking vs. ablebodied controls. However, highfunctioning individuals with limb loss such as military service members often walk with the same metabolic cost as controls. Here we used a 3D computer model and optimal control simulation approach to test the hypothesis that transtibial limb loss in and of itself causes an increase in metabolic cost of walking. We first generated N = 36 simulations of walking at 1.45 m/s using a “prelimb loss” model, with two intact biological legs, that minimized deviations from ablebodied experimental walking mechanics with minimum muscular effort. We then repeated these simulations using a “postlimb loss” model, with the right leg’s ankle muscles and joints replaced with a simple model of a passive transtibial prosthesis. No other changes were made to the postlimb loss model’s remaining muscles or musculoskeletal parameters compared to the prelimb loss case. Postlimb loss, the gait deviations on average increased by only 0.17 standard deviations from the experimental means, and metabolic cost did not increase (3.58 ± 0.10 J/m/kg prelimb loss vs. 3.59 ± 0.12 J/m/kg postlimb loss, p = 0.65). The results suggest that transtibial limb loss does not directly lead to an increase in metabolic cost, even when deviations from ablebodied gait mechanics are minimized. High metabolic costs observed in individuals with transtibial limb loss may be due to secondary changes in strength or general fitness after limb loss, modifiable prosthesis issues, or to prioritization of factors that affect locomotor control other than gait deviations and muscular effort.
Introduction
Numerous studies report on gait deviations and a high metabolic cost of walking in individuals with lower limb amputation (e.g., Pinzur et al., 1992; Sanderson & Martin, 1997; Rábago & Wilken, 2016). These deviations from the mechanics and energetics of ablebodied gait are suspected to affect mobility, quality of life, and risk for chronic musculoskeletal conditions after limb loss (Molen, 1973; Gailey et al., 2008). Minimizing gait deviations and maintaining a reasonably economical gait are therefore important goals of rehabilitation and prosthesis prescription after limb loss. Deviations from ablebodied gait mechanics and energetics are often presumed to result directly from the amputation. However, the studies documenting these deviations are necessarily crosssectional in nature. There are no longitudinal data on actual biomechanical and metabolic changes pretopost limb loss.
Longitudinal studies on gait mechanics pre vs. postlimb loss are impractical but can be approximated in computer simulations (Zmitrewicz, Neptune & Sasaki, 2007; Handford & Srinivasan, 2016). In a previous study, motivated by observations that highfunctioning children and military service members with transtibial limb loss do not have a greater metabolic cost of walking than ablebodied controls (Jeans, Browne & Karol, 2011; Russell Esposito et al., 2014; Jarvis et al., 2017), we used optimal control simulations to perform an in silico longitudinal study of walking mechanics and energetics pre vs. postlimb loss. The results suggested that maintenance of residual limb muscle strength can maintain the prelimb loss metabolic cost of walking, with only small deviations from the prelimb loss gait mechanics (Russell Esposito & Miller, 2018). However, the computer model in this study was limited to only the sagittal plane. Many of the reported gait deviations in individuals with limb loss occur outside the sagittal plane (Rábago & Wilken, 2016), and lateral balance in normal human walking appears to have a modest metabolic cost (Donelan et al., 2004). Neglecting threedimensional motion when making inferences on the energetics of walking with a prosthesis may therefore be a questionable assumption.
Therefore, the purpose of this study was to perform an in silico longitudinal study pre and posttranstibial limb loss using a threedimensional musculoskeletal model to determine the change in metabolic cost, with its typical definition as the energy expended to translate a unit mass by a unit distance (SchmidtNelson, 1972; Tucker, 1975). Based on previous experiments and simulations (Russell Esposito et al., 2014; Russell Esposito & Miller, 2018), we hypothesized that the metabolic cost of walking with a transtibial prosthesis postlimb loss would be similar to the prelimb loss metabolic cost when the prelimb loss muscular properties and body composition were maintained.
Methods
Prelimb loss model
The model was a modified version of the OpenSim (Delp et al., 2007) model described by Rajagopal et al. (2016). The prelimb loss model is visualized in Fig. 1 and consisted of pelvis, trunk, thigh, shank, foot, toes, upper arm, and forearm rigid bodies that articulated with 31 degrees of freedom. The lumbar joint and lower limb joints were actuated by 84 total muscles (42 bilateral pairs). The shoulder and elbow joints were actuated by ideal torque generators with activation dynamics. The model had a total body mass of 75.4 kg and height 1.75 m. Several changes or additions were made to the original model, which lacked some elements necessary or beneficial for optimal control simulations, to make it well suited for direct collocation simulations in Moco software (Dembia et al., 2020):

The original “Millard2012Equilibrium”type muscles were replaced with “DeGrooteFregly2016”type muscles, which have activation and contractile dynamics suitable for gradientbased optimization in Moco (De Groote et al., 2016).

Wrapping objects that defined muscle paths on the skeleton were replaced with “Conditional Path Points” that gave similar relationships between the skeletal pose and muscle lengths and moment arms. Moco does not currently support wrapping objects.

The studies that motivated the present study examined samples of young men (Russell Esposito et al., 2014; Jarvis et al., 2017). Therefore, muscle optimal fiber lengths and unloaded tendon lengths were adjusted so that the model produced maximum isometric joint torques at average joint angles for young men from Anderson, Madigan & Nussbaum (2007). Specifically, optimal fiber length was adjusted unless the adjustment required shortening it to a length that prevented the muscle from producing force across a realistic range of motion, in which case unloaded tendon length was adjusted instead. The adjustments made from the original model’s parameters were small, e.g., under three cm for the longest muscles in the model.

The original model did not include musclespecific time constants for activation dynamics. These muscle activation dynamics time constants were defined as functions of muscle mass and fasttwitch fiber ratio using fibertype data from Miller (2018) and time constant functions from Winters & Stark (1985).

The original model did not include passive forces limiting joint ranges of motion beyond the parallel elastic components of muscles. Therefore, coordinate forces in the form of passive torsional springdampers were added to the joints representing lumped effects of ligaments and other nonmuscular joint structures on range of motion (Anderson, 1999).
Additionally, ground contact was modeled by 11 sphereshaped HuntCrossley contact elements (Sherman, Seth & Delp, 2011) located on the plantar surface of each foot (Fig. 1). The modulus of the contact elements (3.06 MPa) and damping coefficient (2.0 s/m) were set so that the heel contact sphere deformation and energy return were similar to the heel region of a human foot in an athletic shoe (Aerts & De Clercq, 1993). The frictional contact forces were a Stribeck friction model (Sherman, Seth & Delp, 2011) with static, dynamic, and viscous friction coefficients of 0.8, 0.8, and 0.5, respectively, and transition velocity 0.2 m/s.
The model’s muscle volumes were determined as functions of the model’s mass and height using average MRI measurements of muscle volume for healthy young adults from (Handsfield et al., 2014). Maximum isometric forces were the determined by dividing the volume by the muscle’s optimal fiber length to determine the physiological crosssectional area, and multiplying this area by a specific tension of 60 N/cm^{2}, the same value assumed in the original model (Rajagopal et al., 2016). This specific tension is larger than the values typically determined from in vitro singlefiber experiments, which average ∼13 ± 6 N/cm^{2} in human muscle (Rospars & MeyerVernet, 2016), but is similar to the in vivo result on human quadriceps of 55 ± 11 N/cm^{2} for healthy adult men from O’Brien et al. (2010). The model’s maximum isometric knee extension torque when emulating typical conditions of human dynamometry experiments on hip position, voluntary activation, and antagonism (Kubo et al., 2006), was 3.92 Nm/kg, which compares reasonably well to dynamometry experiments on human subjects that reported torque–angle data (Fig. 2).
Muscle metabolic rates were calculated as functions of muscle activation and fiber velocity using the model described by Umberger, Gerritsen & Martin (2003). The muscle rates were summed and added to an assumed basal metabolic rate of 1.0 W/kg to determine the model’s wholebody metabolic rate. The metabolic rate was then averaged over time and divided by average walking speed to determine the metabolic cost, i.e., the energy expended per unit distance traveled.
Postlimb loss model
The model of limb loss and a generic prosthesis was based on our previous 2D simulations (Russell Esposito & Miller, 2018). The 13 muscles spanning the right ankle were removed and the doubleexponential joint torque–angle relationships representing the biological ankle, subtalar, and toe joints were replaced by linear relationships of the form: (1)$\tau =k\theta b\stackrel{\u0307}{\theta}$
where τ is the passive moment at the joint, θ is the joint angular position, and k and b are the stiffness and damping constants. The mass and moment of inertia of the prosthetic leg’s bodies below the knee were reduced to 65% and 40%, respectively, of the intact leg values, representing the typical mass and mass distribution in a transtibial prosthesis (Smith & Martin, 2013). To compute the postlimb loss model’s basal metabolic rate, the minimum heat rates of all 84 prelimb loss muscles (1.0 W per kg muscle mass; (Umberger, Gerritsen & Martin, 2003) were subtracted from the total basal metabolic rate (1.0 W per kg biological mass) to determine the nonmuscular basal metabolic rate (54 W average). The nonmuscular rate was assumed to be unchanged postlimb loss, and the minimum heat rates of the remaining 71 postlimb muscles were added to the nonmuscular rate to determine the postlimb loss basal metabolic rate. The postlimb loss model was otherwise identical to the prelimb loss model.
A limitation of this prosthesis model is that the prosthesis was rigidly attached to the residual limb, which is more similar to an osseointegrated prosthesis than a traditional prosthesis with deformation at the socketresiduum interface. Currently, osseointegration is primarily used in patients with transfemoral limb loss who have negative outcomes using a traditional prosthesis, although transtibial osseointegration is feasible (Frossard, Leech & Pitkin, 2020). In model development, a model that allowed pistoning between the prosthesis and residual limb was tested. The force–deformation relationship of the distal residuum was a doubleexponential: (2)$F={a}_{1}exp\left({a}_{2}y\right)+{a}_{1}exp\left(2{a}_{2}y\right)b\stackrel{\u0307}{y}$
where F is the soft tissue load at the socketresiduum interface and y is the soft tissue deformation along the long axis of the residuum. Values for a_{1} = 7.8 and a_{2} = 350 were set to approximate the average deformation profile from 0–100% bodyweight loading from fluoroscopy data (Darter, Sinitski & Wilken, 2016). Damping b = 100 N/(m/s) was set to suppress highfrequency oscillations of the prosthesis during the swing phase. Adding pistoning changed the model’s metabolic cost of walking by only 0.5% and did not affect its accuracy in tracking experimental gait data, so the rigid socketresiduum interface without pistoning was used for simplicity.
Simulations
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 ‘PreLimb Loss Model’–‘PostLimb Loss Model’ with different parameter values. “Subject #1” used the default values of the prelimb loss model’s parameters. To generate additional prelimb 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 prelimb loss metabolic costs of roughly 3.4–3.8 J/m/kg. These models were then converted to postlimb loss models (‘PostLimb Loss Model’). For each subject’s postlimb 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 subjectspecific 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 timevarying muscle excitations and associated model states that minimized a cost function J that was the weighted sum of gait deviations plus muscular effort: (3)$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 righthand 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 betweensubjects of the experimental data, averaged over the stride cycle. The N_{1} = 37 variables included in this term were the six translations and rotations of the pelvis, the 25 joint angles, and the 3D 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 nonsagittal shoulder motions did not track experimental data due to low confidence in accurately measuring these motions with skin markerbased 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 righthand of Eq. (3) is the “muscular effort per unit distance”: ${u}_{i}\left(t\right)$ ranges from 01 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_{2} = 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 minimumeffort 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$\u2215\stackrel{\u0304}{m}$, where $\stackrel{\u0304}{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$\u2215\stackrel{\u0304}{m}$ was therefore used in all simulations. For the upper limb torque generators, m_{i} was set equal to $\stackrel{\u0304}{m}$. The same value of w = 100$\u2215\stackrel{\u0304}{m}$ was used in all simulations for all subjects, both pre and postlimb 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 predetermined result. Keeping the same value of w in both the pre and postlimb 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 highresolution muscle excitations and relatively fast computational speeds (Van den Bogert, Blana & Heinrich, 2011). The constraint tolerance for convergence of the optimizations was 10^{−4}.
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. postlimb 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 prelimb 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. postlimb loss.
Prelimb loss simulations
Simulations with each subject’s prelimb loss model described in ‘PreLimb 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 prelimb loss model with different values of w = 10$\u2215\stackrel{\u0304}{m}$, 100$\u2215\stackrel{\u0304}{m}$, and 1000$\u2215\stackrel{\u0304}{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 101node 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 timevarying 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 prelimb loss simulation that converged on a local minimum with a high metabolic cost could bias the postlimb loss simulation towards having a comparatively low metabolic cost.
Postlimb loss simulations
Simulations with the postlimb loss models were performed identically to the prelimb 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_{1} = 34 in Eq. (3)) because the tracking targets were ablebodied data, and (ii) the 13 muscles of the right ankle were removed (N_{2} = 79 in Eq. (3)). To approximate the patientspecific 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 postlimb 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 postlimb loss simulations, the same subject’s prelimb loss simulation result was used as the initial guess. The cost function in the postlimb loss simulations used the values of ${\mu}_{i}\left(t\right)$ from the ablebodied subjects of Miller et al. (2014) as tracking targets. The postlimb loss model was therefore not explicitly trying to track gait data from individuals with limb loss. The goal of these postlimb loss simulations was to walk with minimal deviations from a typical ablebodied 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 prelimb loss metabolic cost was scaled by the prelimb loss biological mass, and the postlimb loss metabolic cost was scaled by the postlimb loss biological mass. The biological mass was 75.4 kg prelimb loss and 72.8 ± 0.3 kg postlimb loss. Scaling by biological body mass (excluding the prosthesis mass in the postlimb loss case) was done to avoid favoring the hypothesis of similar metabolic costs pre vs. postlimb loss, which would have been favored by scaling the postlimb 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 ttest to determine if the null hypothesis of no change could be rejected. Equivalence of the pre vs. postlimb loss metabolic costs was then tested using two onesided tests (Lakens, Schell & Isager, 2018): $d\le {d}_{min}d\ge {d}_{min}$
where d is the effect size in the data for the postlimb 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. postlimb loss metabolic costs are practically equivalent.
Results
Gait mechanics
The prelimb loss model tracked the mean experimental 3D kinematics and GRFs with an average deviation of 0.51 SD. The postlimb loss model had an average deviation of 0.65 SD, excepting the joint angles of the prosthetic limb which were not tracked. Small nonzero GRFs appear in the “swing” phase from the foot briefly clipping the ground as the leg is swung forward.
Figures 4–7 present the time series of pre and postlimb loss GRFs, lumbopelvic kinematics, lower limb kinematics, and upper limb kinematics, respectively, averaged over subjects. Compared to prelimb loss, the peak anterior GRF on the prosthetic limb was reduced after limbloss (26.7 ± 0.5 vs. 21.2 ± 0.6%bodyweight, p < 0.001; Fig. 4), and the second peaks of the vertical GRF became bilaterally asymmetric, with on average 26% less vertical force under the prosthetic limb vs. the intact limb (106 ± 7 vs. 132 ±1%bodyweight, p < 0.001; Fig. 4). Similar differences in GRFs have been reported in highfunctioning individuals with transtibial limb loss vs. controls (Sanderson & Martin, 1997). Although kinematics were not evaluated statistically, these changes in GRFs were accompanied by numerous small/subtle changes in the lumbopelvic and lower limb kinematics postlimb loss of ∼5° or less, but few larger changes (Figs. 5 and 6). Limb loss did not affect upper limb kinematics (Fig. 7).
Metabolic cost
There was no significant change in the metabolic cost of walking prelimb loss vs. postlimb loss (3.58 ± 0.10 J/m/kg prelimb loss vs. 3.59 ± 0.12 J/m/kg postlimb loss, p = 0.65). The change in cost on average (+0.006 J/m/kg) fell within the equivalence bounds of ± 0.068 J/m/kg at the 95% confidence level (p < 0.001; Fig. 8). Twenty of the 36 subjects (56%) decreased metabolic cost postlimb loss. The changes in metabolic cost were generally small: of the 36 subjects, only seven (19%) had a change outside the equivalence bounds of ± 0.068 J/m/kg (Fig. 9).
When the postlimb loss metabolic cost was scaled by the total system mass (including the prosthesis mass) instead of the biological body mass, the change in metabolic cost postlimb loss averaged −0.028 J/m/kg. This change was significant (p = 0.030) but still fell within the equivalence bounds (Fig. 10). When no mass scaling was performed, metabolic cost decreased postlimb loss (average −8.8 J/m). This change was significant (p < 0.001) and also outside the equivalence bounds (Fig. 10).
Table 1 presents the pre and postlimb loss metabolic costs and their work and heatrelated components. The right ankle muscles on average accounted for 65 J of metabolic energy, or 17% of the prelimb loss metabolic cost. These muscles were absent and consumed no energy postlimb loss, but these savings were coupled with increases in energy expenditure by most other muscles postlimb loss (Fig. 11): the energy expended by the remaining 71 muscles postlimb loss exceeded their prelimb loss energy expenditure by 52 J on average. The postlimb loss model still expended less total metabolic energy than the prelimb loss model (378 J vs. 391 J on average), but when divided by the slightly lower biological body mass postlimb loss (72.8 kg vs. 75.4 kg), the metabolic cost of transporting a unit biological mass by a unit distance did not differ prelimb loss vs. postlimb loss.
Prelimb loss  Postlimb loss  

Quantity (J/m/kg)  Whole body  Right ankle  Whole body  Prosthesis  
METABOLIC COST:  3.58 ± 0.10  0.59 ± 0.02  3.59 ± 0.12  –  
NET WORK:  0.55 ± 0.04  0.28 ± 0.01  0.61 ± 0.03  −0.01 ± 0.00  
Concentric:  1.36 ± 0.04  0.39 ± 0.01  1.28 ± 0.05  0.36 ± 0.02  
Eccentric:  −0.81 ± 0.03  −0.11 ± 0.01  −0.67 ± 0.03  −0.37 ± 0.02  
NET HEAT:  3.03 ± 0.08  0.31 ± 0.01  2.98 ± 0.09  –  
Activation:  2.14 ± 0.09  0.24 ± 0.01  2.13 ± 0.09  –  
Shortening:  0.40 ± 0.01  0.07 ± 0.01  0.33 ± 0.02  –  
*Basal:  0.49 ± 0.01  –  0.51 ± 0.01  –  
Quantity (J)  
METABOLIC COST:  391 ± 11  65 ± 2  378 ± 12  –  
NET WORK:  60 ± 5  31 ± 2  64 ± 3  −0.5 ± 0.2  
Concentric:  149 ± 5  42 ± 2  135 ± 5  26.6 ± 1.6  
Eccentric:  −89 ± 3  −12 ± 1  −71 ± 3  −27.1 ± 1.6  
NET HEAT:  331 ± 9  34 ± 1  314 ± 10  –  
Activation:  234 ± 10  26 ± 1  225 ± 10  –  
Shortening:  44 ± 2  8 ± 1  35 ± 2  –  
*Basal:  54 ± 1  –  54 ± 1  – 
Prosthesis parameters
The optimized prosthesis stiffness parameters (k in Eq. (1)) averaged 426 ± 51, 540 ± 63, and 306 ± 90 Nm/rad for the ankle, subtalar, and toe joints, respectively. The optimized damping parameters (b in Eq. (1)) averaged 0.06 ± 0.06, 1.27 ± 0.52, and 0.73 ± 0.91 Nm/(rad/s) for the ankle, subtalar, and toe joints, respectively. Owing to the low damping, the prosthesis had highefficiency energy return, storing an average of 27.1 J of mechanical strain energy and returning an average of 26.6 J (Table 1).
Figure 11 includes the energy expended by each muscle for Subject 06, who had the largest absolute increase in metabolic cost (+0.27 J/m/kg, or about +8%). No particular muscle or small group of muscles explained this subject’s increase in cost. This subject did not have outlier values for their optimized prosthesis parameters: ankle, subtalar, and toe stiffness parameters were 409, 544, and 233 Nm/rad, and the associated damping parameters were 0.20, 1.20, and 0.35 Nm/(rad/s).
Discussion
The purpose of this study was to determine the change in metabolic cost when walking with a transtibial prosthesis in threedimensional optimal control simulations. With a computer modeling approach, this study performed a pre vs. postlimb loss comparison that is difficult to make in experiments on human subjects. Our hypothesis of similar pre and postlimb loss metabolic costs using this threedimensional model was supported in walking simulations. While we did not test for the effect of muscle strength as we did previously with 2D models of limb loss (Russell Esposito & Miller, 2018), the postlimb loss model was identical to the prelimb loss model in all parameter values including muscle strength, body mass, and mass distribution. The only difference between the pre and postlimb loss models was the loss of the ankle muscles on the amputated side, and the replacement of the biological foot/ankle complex with a prosthetic one. This result joins a growing body of evidence suggesting that belowknee limb loss per se does not directly or inevitably cause the relatively high metabolic cost of walking seen in the general limb loss population (Jeans, Browne & Karol, 2011; Russell Esposito et al., 2014; Russell Esposito & Miller, 2018). An increase in metabolic cost may instead be due to changes that tend to accompany limb loss such as loss of muscle strength or other fitness and mobilityrelated factors such as balance or endurance (Hewson, Dent & Sawers, 2020).
This study was hypothesisdriven and addressed a research question with a yes/no answer, “Are metabolic costs equivalent pre vs. postlimb loss?” The suggested mechanism for this result is that limb loss removes the metabolic cost of the lost ankle muscles but increases metabolic cost of other muscles to compensate for the lost function (Table 1; Fig. 11), with the caveat that these compensations do not necessarily need to exceed the prelimb loss metabolic cost. The natural followup questions of why metabolic cost did not change on average (Fig. 8) or why some of the 36 “subjects” had different magnitudes and directions of change (Fig. 9) are interesting but are difficult to address with confidence from the present data. For example, the subject with the largest increase in metabolic cost postlimb loss did not have an obvious muscle or muscle group that appeared to explain this increase (Fig. 11), and their model parameters, prosthesis parameters, and tracking accuracy in either simulation were not unusual. Isolating the cause of this subjectspecific response would require additional simulations to discover why metabolic cost sometimes increased or decreased by relatively large amounts for some of the subjects, and to determine if this effect can be induced in a predictable way.
The model’s tracking accuracy of human gait data (Figs. 4–7) was similar to other studies that have used this simulation method with simpler models (Van den Bogert, Blana & Heinrich, 2011; Nguyen et al., 2020) and is generally better than what can be achieved with the traditional “single shooting” approach to optimal control simulations in biomechanics (e.g., Miller, Esterson & Shim, 2015), owing to the ability to use highresolution muscle excitations that more closely resemble human neuromuscular control. The differences in ground reaction forces pre vs. postlimb loss were similar to those reported in highfunctioning individuals with transtibial limb loss vs. controls (Sanderson & Martin, 1997) but the postlimb loss simulations did not exhibit the magnitude of kinematic deviations typically seen in highfunctioning individuals with transtibial limb loss (Rábago & Wilken, 2016), most of which could be reasonably expected to influence metabolic cost. As noted earlier, nonzero GRF appeared in the “swing” phase from the foot briefly clipping the ground as the leg is swung forward. These GRF were brief and small (∼50 N), are not unusual in these types of simulations (e.g., Anderson & Pandy, 2001), and appeared in both the pre and postlimb loss simulations, and therefore did not bias the pre vs. postlimb loss changes in metabolic cost. Within the limitations of the present 3D model, it was possible to walk after unilateral limb loss with a passive prosthesis without increasing metabolic cost and without major deviations from prelimb loss gait mechanics. This result suggests that large gait deviations, even in highfunctioning individuals with unilateral transtibial limb loss, may not necessarily be due to limb loss directly. Increasing the weight on muscular effort in the cost function (parameter w in Eq. (3)) will generally increase gait deviations and decrease metabolic cost within a subject, but the magnitude of gait deviations did not correlate strongly with metabolic cost between subjects (r^{2} < 0.01 for both pre and postlimb loss). It is difficult to comment further on possible relationships between specific gait mechanics like external work and metabolic cost (e.g., Houdijk et al., 2009), as this study was not designed to determine these relationships. Relatedly, the cost function used here explicitly considered only two factors, gait deviations and muscular effort. Other unmodeled factors such as balance, comfort, and agility are also important to prosthesis users (Bell et al., 2021) and could involve tradeoffs with gait deviations.
The main finding of this study of no increase in metabolic cost postlimb loss agrees with our previous conceptually similar simulation study that used a 2D sagittal plane model (Russell Esposito & Miller, 2018). The 3D model expends muscular effort moving and controlling itself in the nonsagittal planes, vs. the 2D model where these planes did not exist and all effort was expended in the sagittal plane. The 2D model was not simply a 3D model that only moved in the sagittal plane, it was a model for which no other planes were defined. The present 3D result confirms that the previous 2D result was not due to the neglected nonsagittal mechanics and gives more confidence in the generalizability of the previous result to real human walking since it was replicated here with a more complex 3D model. We caution though that the similarity in results between these two studies with 2D and 3D models should not be taken to imply that the mechanics and energetics of limb loss can necessarily be modeled realistically in 2D as a general rule. A more conservative conclusion would be that the energy cost of frontal plane balance is similar with and without transtibial limb loss. There is limited information in the limb loss literature on the energy cost of frontal plane balance. Ijmker et al. (2014) reported that adding lateral stabilization via springcords at the waist caused small reductions in metabolic cost in individuals with transtibial limb loss (−5%) and in controls without limb loss (−3%), but an increase of 6.5% in individuals with transfemoral limb loss. However, these data do not necessarily imply that the cost/task of controlling balance in the frontal plane is unaffected by transtibial limb loss.
The muscle energy model includes the fraction of fasttwitch fibers as an input parameter, which affects the computed metabolic cost primarily through muscle activation/maintenance heat rate, which is a linear function of the fasttwitch fiber fraction (Umberger, Gerritsen & Martin, 2003). Average values from the literature were assumed for the model’s fibertypes and the same values were used for all subjects. This assumption would not be expected to have a major influence on the present results given the repeatedmeasures design. Our previous 2D study used random perturbations to define fibertypes for each subject and produces the same overall finding as the present study (Russell Esposito & Miller, 2018). However, extreme cases, such as a subject with an unusually low fraction of fasttwitch fibers in the right ankle muscles and unusually high fraction elsewhere, could feasibly produce different results than the average result reported here.
Further concerning model complexity and assumptions, the research question here was whether limb loss per se increases the metabolic cost of walking. We therefore developed a modeling approach for testing the effect of limb loss in isolation of other factors that may affect metabolic cost, such as secondary changes in fitness, a poorly adjusted prosthesis, pain during walking, etc. While the present model is on the upper end of the complexity spectrum for models used previously in these types of simulations, the model and simulation approach required simplifying assumptions compared to real humans that should be considered when generalizing the results to real human walking. For example, as noted earlier we modeled a rigid prosthesis attachment, which was justified for the sake of simplicity by the small effect of a nonrigid attachment on the primary outcome variable (metabolic cost), but the present results would not necessarily generalize to individuals who have a loose prosthesis fit or substantial discomfort at the socketresiduum interface. For another example, we did not directly minimize metabolic cost in the cost function of the simulations because this is not currently possible in the Moco software (Dembia et al., 2020) and instead used a surrogate “muscular effort” cost. This assumption is justified by the fact that the approach predicted realistic magnitudes of metabolic cost and realistic changes in metabolic cost in comparison to human experiments (Fig. 3), but we cannot rule out the possibility that other plausible cost functions would produce different results. Lastly, the perturbations to model parameters defining different “subjects” were done only to a subset of parameters. These perturbations are intended to increase confidence in the effect of the independent variable (limb loss). Other unperturbed parameters would surely affect metabolic cost and the variance between subjects in Fig. 8, but whether they would systematically affect the change in metabolic cost with limb loss in this repeatedmeasures design is unclear and speculative, excluding extreme cases such a “subject” with ankle muscle and joint properties that are unusually economical for walking.
The prosthesis stiffness and damping parameters were optimized in the postlimb loss simulations to represent the clinical process of fitting/adjusting the prosthesis to the user, similar to our previous 2D study (Russell Esposito et al., 2014). From a technical perspective, this process also ensures that the simulated gait deviations and metabolic costs are not inflated due to inadequacy in generic, unadjusted prosthesis parameters for a given subject. For all subjects the optimized stiffnesses fell within the range of stiffnesses for typical commercial prosthetic feet of ∼300–900 Nm/rad (Lehmann et al., 1993) even though bounds on these parameters allowed values outside this range, lending some confidence that the prosthesis model, although fairly simple, was reasonably realistic. The optimized damping parameters were low, as would be expected for modern prosthetic feet with efficient mechanical energy return, although some damping can theoretically be beneficial for reducing damage from highfrequency vibrations (Burnett et al., 2021). The modeling approach here could feasibly be adapted to subjectspecific simulations for determining ideal prosthesis characteristics to accomplish particular patientspecific goals, similar to proposed modelbased approaches on personalized gait modification and joint loading in osteoarthritis (Fregly et al., 2007).
The prosthesis model here with its low damping had an average energy return of 98% (Table 1), which is rather high compared to most conventional prosthetic feet (e.g., Prince et al., 1998) and likely contributed to the postlimb loss metabolic costs, but is not necessarily unrealistic. Recent “ProFlex” (Össur, Reykjavík, Iceland) prosthetic feet with carbon fiber leaf springs offer considerably greater energy return than conventional feet (Childers & Takahashi, 2018), in some cases appearing to approach 100% return (Heitzmann et al., 2018), and springbased “energyrecycling” feet with a small (0.8 W) battery can also return energy at very high rates (Zelik et al., 2011). However, prosthetic mechanical energy is not a direct component of metabolic cost and this study was not designed for causal inferences on this rather complicated relationship. For example, powering a prosthesis to act as a net generator of mechanical energy does not necessarily reduce metabolic cost (Quesada, Caputo & Collins, 2016), removing the ability of the ankle to store and return energy does not necessarily increase metabolic cost (Vanderpool, Collins & Kuo, 2008), and increasing energy return does not necessarily reduce metabolic cost (Zelik et al., 2011).
Metabolic cost data in biomechanics are often scaled by body mass for analysis. With limb loss subjects there is a choice of whether this scaling includes or excludes the mass of the prosthesis. Some studies have scaled by biological body mass excluding the prosthesis mass (Russell Esposito et al., 2014), others by total mass including the prosthesis mass (Jarvis et al., 2017), and others still by an unspecified mass (Jeans, Browne & Karol, 2011). This choice is nontrivial for research questions like the present study: for a typical transtibial prosthesis mass of 1–2 kg (Mattes, Martin & Royer, 2000), the scaling choice affects metabolic cost by ∼1–3% for body masses ranging from 70–100 kg, which can feasibly change conclusions when the minimum effect of interest is ∼2%. Here regardless of the choice of scaling, limb loss did not increase metabolic cost, but different conclusions could be reached depending on the scaling choice (Fig. 10). With no scaling, limb loss decreased metabolic cost, and the decrease was greater than the 2% minimum change of interest. With total mass scaling, inclusive of both the biological mass and the prosthesis mass, limb loss also decreased metabolic cost, but the decrease was less than the defined meaningful difference of 2%. With biological mass scaling only, limb loss did not change the metabolic cost; the average result was a small increase in cost but this increase was too small to be meaningful. We contend that to address our specific research question on whether limb loss per se changes metabolic cost, scaling by biological mass is the most appropriate and conservative choice since it accounts for change in biologic tissue mass that occurs with limb loss, and divides the energy consumed by a smaller mass than the total mass to avoid favoring the hypothesis of no increase in cost postlimb loss. Metabolic cost scaled by biological mass could be interpreted as a metric of efficiency, the energy expended to walk by the available biological mass that is capable of consuming energy, while scaling by total mass could be interested as a metric of economy, the energy expended to transport a unit mass. Regardless of this choice, many studies on limb loss that involved mass scaling do not report the specific mass that outcome variables were scaled by. There are certainly cases where scaling by total mass rather than biological mass is a sensible approach. We are not advocating that limb loss data should always be scaled by biological mass regardless of the research question. However, reporting the specific mass used for scaling is critical for interpreting metabolic energy expenditure data and should be mandatory.
The choice of a 2% change in metabolic cost to define the minimum “meaningful” effect of interest was based on a recent study on the reliability of modern portable metabolic units (Guidetti et al., 2018). The smallest “clinically relevant” change for walking with limb loss or for walking in general is not well defined to our knowledge, but is presumably at least as large as the minimum reliably measurable difference used here. The small value used here avoided favoring our hypothesis of no necessary increase in metabolic cost with limb loss. In addition, “small” changes in gross metabolic cost (under 5–10%) are often of interest in prosthesis research (e.g., Zelik et al., 2011; Nolan, 2012).
A final limitation of note is that these results pertain to unilateral transtibial limb loss only, despite any ambiguous references in this report to general “limb loss”. We would expect similar results for simulations of lowerlevel limb loss, e.g., midfoot, but not for higherlevel limb loss e.g., aboveknee or bilateral limb loss, but these expectations are speculative. Relatedly, these results do not necessarily generalize to other more demanding movements such as running or walking with substantial load carriage.
In conclusion, walking with a transtibial prosthesis after belowknee limb loss did not increase the metabolic cost in optimal control simulations with a 3D musculoskeletal model. This finding supports our previous similar finding from a simpler 2D model. The results further suggest that major deviations in the body’s intact joints and segments from ablebodied gait are also not an inevitable consequence of limb loss and may be due to secondary changes in fitness, or to priorities in the motor control of gait other than muscular effort. Lastly, studies on limb loss should report the mass used in scaling outcome variables.
Supplemental Information
Raw data and processing/plotting code
Each subject’s simulation output and Matlab code for generating plots presenting the data.