Muscle–tendon mechanics explain unexpected effects of exoskeleton assistance on metabolic rate during walking

The plantarflexor muscle–tendon units seem tuned for near-optimal efficiency and power production during unassisted locomotion. During normal walking, the ankle plantarflexor muscles produce force nearly isometrically throughout mid-stance, while the Achilles tendon lengthens and stores mechanical energy (Fukunaga et al., 2001). This isometric muscle force production is economical because muscles consume relatively little energy to produce force at constant length (Biewener, 1998; Biewener and Roberts, 2000). At the end of stance, the plantarflexor muscles actively shorten and the Achilles tendon simultaneously recoils (Fukunaga et al., 2001; Ishikawa et al., 2005; Rubenson et al., 2012), generating a significant amount of positive power at push-off (Winter, 1990; Meinders et al., 1998). Elastic energy storage and recovery in the Achilles tendon helps to reduce plantarflexor muscle work (Roberts et al., 1997; Ishikawa et al., 2005). Furthermore, the stiffness of the Achilles tendon, in conjunction with the resting length of the plantarflexor muscle fibers, has been shown to maximize plantarflexor muscle efficiency during walking and running by allowing the muscle fibers to operate at favorable lengths and velocities during positive fiber work production (Roberts et al., 1997; Roberts, 2002; Lichtwark and Wilson, 2007; Taylor, 2007; Lichtwark and Wilson, 2008; Lichtwark and Barclay, 2010; Arnold et al., 2013). Any change to the stiffness of the Achilles tendon can affect the mechanics of the plantarflexor muscle fibers and consequently alter muscle energy consumption (Lichtwark and Wilson, 2007). The architecture of the plantarflexor muscles, the compliance of the Achilles tendon, and the interaction between these mechanisms enables economical operation.

The complexity of these plantarflexor muscle–tendon mechanics poses a challenge for the design of exoskeletons intended to operate in concert with the musculoskeletal system. Previous experiments and simulations of a musculoskeletal model have shown that elastic exoskeletons worn during bilateral hopping significantly reduce plantarflexor muscle force, but not muscle work (Farris and Sawicki, 2012; Farris et al., 2013, 2014; Robertson et al., 2014). Although large reductions were observed in whole-body metabolic rate, estimated metabolic energy consumed by the plantarflexor muscles was not significantly reduced, likely due to unfavorable changes in the operating lengths and velocities of the muscle fibers (Farris et al., 2014). Simulations of a simplified, lumped model of the plantarflexor muscle–tendon units acting in parallel with a passive exoskeleton during walking, with fixed joint kinematics, similarly suggest a disruption to the normal operation of the plantarflexor muscle–tendon units (Sawicki and Khan, 2016). We were curious to see if similar mechanisms could explain the effect of different types of exoskeleton assistance on locomotor coordination and metabolic rate that we observed in a prior study.

We previously conducted an experiment in which subjects walked in eight conditions with different amounts of net work and average plantarflexion torque provided by an exoskeleton worn on one ankle (Jackson and Collins, 2015). We expected that providing net positive exoskeleton work at the ankle joint would replace or augment positive work performed by the plantarflexor muscles and reduce the associated metabolic cost (Donelan et al., 2002; Gordon et al., 2009). We expected that providing plantarflexion torque about the ankle joint, without providing any net work, would offload plantarflexor muscle forces and reduce the metabolic cost associated with force production (Grabowski et al., 2005). Providing increasing amounts of net exoskeleton work decreased metabolic rate as expected. In contrast with our predictions, providing increasing amounts of average exoskeleton torque increased metabolic rate. We thought these surprising results might be explained by changes in the dynamic interactions between muscles and tendons at the assisted joint.

We were unable to explore changes at the muscle–tendon level during assisted walking using direct measurement in our previous study. Although muscle fiber length changes can be measured using ultrasound imaging, the number of muscles that can be imaged is limited. Furthermore, it is not yet feasible to directly measure individual muscle force and metabolic rate during locomotor tasks in humans. An alternative approach for investigating how plantarflexor muscle–tendon mechanics are affected by different exoskeleton behaviors is to conduct simulations with a musculoskeletal model. Driving a musculoskeletal model with experimentally measured electromyography and joint kinematics is one promising simulation technique for generating realistic estimates of muscle–tendon mechanics (Lloyd and Besier, 2003; Arnold et al., 2013; Farris et al., 2014; Markowitz and Herr, 2016). Simulated muscle–tendon mechanics can be fed into models of muscle energy consumption to obtain estimates of muscle-level energetics (Umberger et al., 2003; Bhargava et al., 2004; Umberger and Rubenson, 2011; Uchida et al., 2016). Such estimates could potentially provide an explanation for the observed changes in whole-body energy consumption.

The purpose of this study was to explore how the mechanics and energetics of the plantarflexor muscle–tendon units change when subjected to different perturbations applied by an ankle exoskeleton. We used muscle activity and joint kinematics data to drive simulations of a musculoskeletal model and obtain estimates of muscle-level mechanics and energetics. We focused our musculoskeletal analyses on the soleus because observed changes were most pronounced in this muscle–tendon unit, and it is the muscle–tendon unit most analogous to the exoskeleton. We hypothesized that providing exoskeleton torque without providing any net work detuned the soleus muscle–tendon unit, leading to reduced elastic recoil of the tendon and increased work by the muscle fibers. We hypothesized that providing net positive exoskeleton work, focused at the end of stance, more fully replaced the role of the soleus muscle–tendon unit, thereby reducing energy consumed at the ankle joint and elsewhere. We expected the results from this study to shed light on how exoskeletons should interact with the muscles and tendons to achieve the greatest benefits.

We performed electromyography-driven simulations of a musculoskeletal model to explore changes in plantarflexor muscle-tendon mechanics under a variety of systematically chosen ankle exoskeleton perturbations (Fig. 1). Electromyography and kinematic data were fed into a musculoskeletal simulation that generated estimates of muscle-level mechanics. The simulated muscle–tendon mechanics were input into a muscle-level metabolics model to obtain estimates of individual muscle metabolic rate. We analysed these simulations to gain insight into how the ankle plantarflexor muscle–tendon units are impacted during walking with an ankle exoskeleton.

We previously conducted an experiment exploring the independent effects of a particular form of ankle exoskeleton torque support and work input on human coordination and metabolic energy consumption during walking (Jackson and Collins, 2015). Eight healthy, able-bodied subjects (seven men and one woman; age=25.1±5.1 years; body mass=77.5±5.6 kg; leg length=0.89±0.03 m) wore a tethered, unilateral ankle exoskeleton, capable of providing up to 120 N m of plantarflexion torque (Witte et al., 2015), while walking on a treadmill at 1.25 m s^–1. We ran two separate parameter sweeps. In the first parameter sweep, average exoskeleton torque was increased across conditions, while net exoskeleton work was held constant at approximately zero. In the second parameter sweep, net exoskeleton work rate was increased across conditions, while average exoskeleton torque was held constant. Metabolic rate, ground reaction forces, motion capture marker positions and muscle activities were measured across all conditions. Data from these experiments were used to drive the musculoskeletal simulations of the current study.

A small fraction of electromyographical signals could not be properly analysed due to poor electrode connectivity and a faulty sensor. Raw electromyographical signals were identified as erroneous if they crossed a threshold of 2 mV. An erroneous signal for a specific subject, muscle and condition led us to exclude that signal across all exoskeleton torque conditions or all exoskeleton work conditions, ensuring that averages were always computed across the same subjects for all conditions in the relevant sweep. Approximately 10% of electromyographical data were thereby excluded from the current study.

We drove a generic lower-body musculoskeletal model adapted from a previously published model (Arnold et al., 2010). The model included the pelvis and both legs, with segments and degrees of freedom as defined by Arnold et al. (2010). The bones in this model were created by digitizing bones of an average-height male (Delp et al., 1990; Arnold et al., 2010). We chose this model because it has previously been used to examine muscle fiber dynamics during human walking and running at different speeds (Arnold et al., 2013) and to understand the effects of elastic ankle exoskeletons on the mechanics and energetics of muscles during hopping (Farris et al., 2014).

Of the original 35 lower-limb muscles in the model, we only included the muscles for which we had electromyographic data: lateral gastrocnemius, medial gastrocnemius, soleus, tibialis anterior, vastus medialis, rectus femoris and biceps femoris long head. Each muscle was modelled as a Hill-type muscle, with a single fiber contractile element and series tendon. Muscle-specific parameters included in the model were optimal fiber length, pennation angle at optimal fiber length, tendon slack length and maximum isometric force. These parameters were based on measurements of 21 cadavers (Ward et al., 2009) and the values used for tendon slack length and maximum isometric force were further based on those computed by Arnold et al. (2010). The exact muscle-specific parameters of the generic model used in this study are provided in Table 1. We used 10% tendon strain at maximum isometric force for the lateral gastrocnemius, medial gastrocnemius and soleus based on results from another study that found these values to result in strains that more closely matched ultrasound measurements (Arnold et al., 2013). Maximum fiber contraction velocity was set to 10 optimal fiber lengths per second for all muscles (Arnold et al., 2010).

We performed electromyography-driven simulations of muscle–tendon dynamics during walking with an ankle exoskeleton using the OpenSim musculoskeletal modelling software (version 3.1; Delp et al., 2007). As we were interested in understanding muscle-level mechanics, it was important that individual muscle activation patterns were estimated appropriately. Driving simulations with electromyographical data, rather than estimating muscle excitations by solving a constrained optimization problem, helps to ensure proper estimates of muscle activations (Lloyd and Besier, 2003; Zajac et al., 2003; Buchanan et al., 2005; Arnold and Delp, 2011). Prescribing joint kinematics helps to ensure that total muscle–tendon unit length changes are simulated accurately and that experimentally measured motions are obeyed even when muscles are omitted (Lloyd and Besier, 2003; Arnold et al., 2013; Farris et al., 2014). Muscle-level activations and muscle–tendon unit length changes provide sufficient information to obtain estimates of muscle–tendon mechanics, specifically muscle fiber force, muscle fiber length, muscle fiber velocity and tendon length (Arnold et al., 2010). Given measured muscle activity and joint kinematics, we were able to generate estimates of plantarflexor muscle–tendon mechanics. Muscle–tendon mechanics generated using electromyography-driven simulations, with prescribed joint kinematics, have shown reasonable qualitative agreement with ultrasound measurements (Farris et al., 2014). Such methods have also been shown to successfully match net joint moments measured via dynamometers (Manal et al., 2012).

Optimal fiber length and tendon slack length were scaled to each participant's anthropometry using marker data collected from a static trial, such that they maintained the same ratio as in the generic model. For each participant, the same muscle–tendon parameters were used across all experimental conditions. Data from a single averaged stride, for each participant for each condition, was provided as the input to the simulation. Marker data were fed into OpenSim's inverse kinematics tool, which generated joint angles. A processed version of electromyographical data was used as the control, i.e. excitation, signal in OpenSim's forward simulation tool. The raw electromyographical data were high-pass filtered (fourth order, Butterworth, 20 Hz cut-off frequency, recursive) to remove movement artifact, full-wave rectified and low-pass filtered (fourth order, Butterworth, 6 Hz cut-off frequency, recursive) to smooth the signal (Ferris et al., 2006; De Luca et al., 2010). It was then normalized to maximum muscle activity measured during normal walking, scaled and delayed. The process of selecting the scaling and delay factors is discussed in the next subsection.

To improve the accuracy of our simulations, we optimized the electromyography scaling and delay factors such that the error between muscle-generated ankle joint mechanics and those derived through inverse dynamics was minimized for each subject across the conditions with increasing average exoskeleton torque. We chose to optimize these parameters because they had a large impact on muscle-generated ankle joint mechanics. Muscle-generated ankle joint moments were calculated by summing the joint moments, defined as the product of the tendon force and moment arm, of the medial gastrocnemius, lateral gastrocnemius, soleus and tibialis anterior. Inverse dynamics-derived ankle joint moments were obtained using OpenSim's inverse dynamics tool, which required joint angles from inverse kinematics, measured ground reaction forces, and exoskeleton torques as inputs. Exoskeleton torques were modelled as equal and opposite external torques applied to the shank and the foot. Both computed ankle joint moments were multiplied by ankle joint velocity to obtain muscle-generated and inverse dynamics-derived ankle joint powers. Other studies have reported that the combination of the soleus, medial gastrocnemius and lateral gastrocnemius contribute about 90% of the total ankle plantarflexion moment, and the tibialis anterior contributes more than 50% of the total ankle dorsiflexion moment in the model (Arnold et al., 2013), suggesting that these muscles are sufficient for generating realistic ankle joint mechanics. We did not, however, expect a perfect match between the two methods (Hicks et al., 2015).

To obtain the optimal values of the scaling and delay factors for the muscles acting about the ankle joint, we performed gradient descent optimization. In order to address differences across subjects, we used scaling and delay factors that were subject specific. For a given subject, the same delay was used for all muscles, while a different scaling factor was used for each muscle. Peak muscle activation during walking, relative to maximum voluntary contraction of that muscle, varies significantly across muscles, therefore suggesting the importance of muscle-specific scaling factors (Perry and Burnfield, 2010). Differences in electromechanical delay across muscles is a more complicated issue (Corcos et al., 1992; Hug et al., 2011). While studies have shown that the delay may be muscle dependent (Conchola et al., 2013), we were able to achieve sufficiently accurate timing of joint moments and powers without such added complexity. Furthermore, previous studies have used a single electromechanical delay across muscles and subjects and obtained reasonable results (Lloyd and Besier, 2003; Arnold et al., 2013).

In total, there were five optimization parameters for each subject: the delay, the medial gastrocnemius scaling factor, the lateral gastrocnemius scaling factor, the soleus scaling factor, and the tibialis anterior scaling factor. The root mean square errors (RMSEs) between the muscle-generated ankle joint moments and powers and the inverse dynamics-derived ankle joint moments and powers were used to quantify the quality of fit. The norm of the RMSEs across the five increasing average exoskeleton torque conditions was chosen as the objective function. The optimized parameters for each subject are provided in Table 2. Because our simulations only include three muscles that cross the knee and hip joints, muscle-generated knee and hip joint mechanics should not be expected to match inverse dynamics-derived knee and hip joint mechanics (Arnold et al., 2013). Therefore we did not optimize the scaling factors for these three muscles, but estimated them as the percentage of the maximum voluntary contraction produced during normal walking observed in other experiments (Perry and Burnfield, 2010).

The optimized parameters produced reasonable ankle joint moments and powers (Fig. 2). The average RMSE, over subjects and conditions, between muscle-generated and inverse dynamics-derived ankle joint moments was 0.13 N m kg^–1, which was 11% of the average peak of the inverse dynamics-derived ankle joint moment (1.2 N m kg^–1). The average RMSE between muscle-generated and inverse dynamics-derived ankle joint powers was 0.19 W kg^–1, which was 9% of the average peak of the inverse dynamics-derived ankle joint power (2.2 W kg^–1). Muscle-generated ankle joint moments were found to be within two standard deviations of inverse dynamics-derived ankle joint moments, on average, which has been considered acceptable by other researchers (Hicks et al., 2015). The error in the timing of peak subject-averaged joint moments and powers had a maximum value of 1.6% of the gait cycle across all conditions. We were most interested in trends in ankle joint moments and powers with increasing average exoskeleton torque and net exoskeleton work, so an exact match in the absolute values of muscle-generated and inverse dynamics-derived ankle joint mechanics was not necessary.

We used the results of the electromyography-driven simulations to estimate the energy consumed by each muscle using a modified version of Umberger's muscle metabolics model (Umberger et al., 2003; Umberger, 2010; Uchida et al., 2016). The metabolics model contains three different heat rates: the activation/maintenance heat rate, the fiber shortening/lengthening heat rate, and the fiber mechanical work rate. These heat rates depend in part on the muscle's excitation, activation, fiber length, fiber velocity and fiber force. We explored the effect of each heat rate on the total metabolic rate for the muscles under consideration. Additionally, we summed the metabolic rates for each of the muscles simulated in our study and investigated how well estimated trends in individual and summed muscle metabolic rates explained changes in whole-body metabolic rate.

A lack of comparative experimental data makes it difficult to validate metabolics models. Other studies have validated their metabolics estimates by comparing simulated whole-body metabolic rate, defined as the sum of the individually simulated muscle metabolic rates, and indirect calorimetry (Umberger et al., 2003; Markowitz and Herr, 2016). As our study only includes a subset of potentially costly muscles, we did not expect the sum of metabolic rates of these muscles to accurately represent absolute changes in whole-body metabolic rate. Furthermore, we were most interested in trends across the different experimental conditions, as opposed to absolute differences. For these reasons, we limit our validation to the percentage change in the sum of the individually simulated muscle metabolic rates.

The version of Umberger's metabolics model that is implemented in OpenSim is configurable, and we chose to use the original version of Umberger's metabolics model (Umberger et al., 2003) with two modifications introduced by Uchida et al. (2016), as these modifications provided more accurate estimates compared with indirect calorimetry in similar studies. The first modification was the addition of a model of orderly fiber recruitment. Umberger's model assumes that the ratio of slow- to fast-twitch fibers that are recruited is equal to the ratio of slow- to fast-twitch fibers comprising the muscle. In the modified model, the ratio of slow- to fast-twitch fibers that are recruited instead varies with excitation so that fast-twitch fibers are increasingly used as excitation increases (Bhargava et al., 2004). The second modification was that the total metabolic rate at any time could not be negative (theoretically, total metabolic rate could be negative if the fiber mechanical work rate was negative and exceeded the total heat rate in magnitude). This change was consistent with the argument that eccentric work cannot cause a net synthesis of ATP (Miller, 2014).

We compared changes in trajectories of muscle fiber forces, muscle fiber lengths, muscle fiber velocities, muscle fiber powers and tendon lengths across all experimental conditions. Average values of outcomes of interest were computed by integrating stride-averaged trajectories over the period of interest and dividing by average stride time of the corresponding trial. Instantaneous values of outcomes of interest were computed by taking the values of the stride-averaged trajectories at the defined times for each subject for each condition. Timing of these values varied across conditions; this was taken into account when calculations were performed. Muscle fiber force was normalized to maximum isometric force as defined in the model, tendon length was normalized to tendon slack length, muscle fiber length was normalized to optimal muscle fiber length, and fiber velocity was normalized to the maximum shortening velocity. Muscle fiber power was calculated as the product of muscle fiber force and muscle fiber velocity at each instant in time. Muscle fiber power and work, as well as all measures of metabolic rate, were normalized to body mass. All outcomes were averaged across subjects. Correlations between estimated percentage changes in muscle-level metabolic rate and measured percentage changes in whole-body metabolic rate were performed on both subject-specific data and data averaged across all subjects. Percentage changes that were calculated on data averaged across subjects are referred to in the text as average percentage changes. Standard deviations represent inter-subject variability.

We first performed a linear mixed-model ANOVA (random effect: subject; fixed effect: average torque or net work) to test for trend significance across experimental conditions in the different measured outcomes. We applied the Jarque–Bera test of normality to ensure samples being compared were normally distributed. For measures that showed trend significance and were normally distributed, we performed paired t-tests to compare two conditions. For measures that showed trend significance but were not normally distributed, we used the Wilcoxon signed rank test to compare two conditions. Across those experimental conditions for which average exoskeleton torque was systematically altered, pairwise statistical comparisons were made with respect to the condition that provided zero average exoskeleton torque. Across those experimental conditions for which net exoskeleton work rate was systematically altered, pairwise statistical comparisons were made with respect to the condition that provided zero net exoskeleton work with a controlled non-zero amount of average exoskeleton torque. After performing pairwise comparisons, we applied the Holm–Šídák step-down correction for multiple comparisons (Glantz, 2012) and used a significance level of α=0.05. The data used to produce our results are publicly available in Dryad (Jackson et al., 2017).

To test the sensitivity of simulated muscle mechanics to model parameters, we conducted a sensitivity analysis. We varied soleus activation and deactivation time constants by ±10% relative to the initial value, maximum fiber contraction velocity by ±20% relative to the initial value, maximum isometric force by ±10% relative to the initial value, tendon slack length by ±5% relative to the initial value, and tendon strain at maximum isometric muscle force by an absolute ±1%. Varying these model parameters as described produced human-like values of muscle mechanics and did not significantly affect trends observed in the outcomes of interest (Figs S1–S6). The figures presented in the supplementary materials on the sensitivity analysis are representative of the changes observed in muscle mechanics and metabolic rates when the reported model parameters were varied.

Perturbing the biological ankle joint with an active exoskeleton altered plantarflexor muscle–tendon mechanics and energetics as well as whole-body coordination patterns. Applying exoskeleton torques in parallel with the biological ankle muscles, without providing any net work, reduced soleus activation and force, but increased muscle fiber excursion, contraction velocity, and consequently, positive muscle fiber work. Increased positive muscle fiber work offset the observed decrease in activation heat rate of the exoskeleton-side soleus. Providing net work with an ankle exoskeleton reduced soleus activation and force during push-off, without significantly altering muscle fiber excursion and velocity, leading to an overall decrease in metabolic rate. Trends in estimated individual and combined muscle metabolic rates correlated well with experimentally observed trends in whole-body metabolic rate.

As exoskeleton average torque was independently increased, the mechanics of the soleus muscle–tendon unit at the assisted ankle joint were disrupted. Average exoskeleton-side soleus muscle activation decreased by 69% during mid-stance and by 21% during late stance across exoskeleton torque conditions (P=8×10^–3 and P=0.02, respectively; Fig. 3A). Average exoskeleton-side soleus muscle fiber force decreased by 65% during mid-stance and by 45% during late stance across exoskeleton torque conditions (P=8×10^–3 and P=2×10^–5, respectively; Fig. 3B). Change in tendon length, from the instant the soleus muscle fiber started lengthening to the instant it transitioned from lengthening to shortening, decreased by 74% across exoskeleton torque conditions (P=1×10^–3; Fig. 3C). Soleus muscle fiber length, at the instant the soleus muscle transitioned from lengthening to shortening, increased by 12% and muscle fiber contraction velocity, at the instant of peak muscle fiber power, increased by 155% across exoskeleton torque conditions (P=1×10^–3 and P=0.02, respectively; Fig. 3D,E). Positive muscle fiber work during late stance increased by 232% across exoskeleton torque conditions (P=0.01; Fig. 3F). Similar trends were observed in the medial and lateral gastrocnemius muscles for a majority of these outcomes, but to a lesser extent (see Appendix, Figs A1, A2).

The positive work rate of the exoskeleton plus tendon decreased with increasing average exoskeleton torque (ANOVA, P=2×10^–3; Fig. 4) while the positive work rate of the soleus muscle increased by 142% across exoskeleton torque conditions (P=0.02). The positive work rate of the combined exoskeleton, tendon and soleus muscle remained relatively unchanged as exoskeleton torque was increased (ANOVA, P=0.9).

The trend in estimated metabolic rate of the exoskeleton-side soleus as average exoskeleton torque was increased appeared to be similar to the trend in measured whole-body metabolic rate (Fig. 5). Average activation/maintenance heat rate decreased by 28% across exoskeleton torque conditions (P=2×10^–3). Average shortening/lengthening heat rate appeared to increase with increasing average exoskeleton torque, but the trend was not significant (ANOVA, P=0.1). Positive mechanical work rate increased by 144% from the condition with no exoskeleton torque to the condition with the second-highest exoskeleton torque (P=8×10^–3). Correlating the estimated percentage change in soleus metabolic rate (⁠⁠) with the experimentally observed percentage change in whole-body metabolic rate (⁠⁠), the best fit line was found to be (R²=0.3, P=2×10^–3). Correlating the average estimated percentage change in soleus metabolic rate (⁠⁠) with the average experimentally observed percentage change in whole-body metabolic rate (⁠⁠), the best fit line was found to be (R²=0.8, P=0.1).

Estimated metabolic rate of the contralateral limb vastus increased with increasing average exoskeleton torque (ANOVA, P=0.02; Fig. 5) and matched trends in measured whole-body metabolic rate. Correlating the estimated percentage change in contralateral limb vastus metabolic rate (⁠⁠) with the experimentally observed percentage change in whole-body metabolic rate (⁠⁠), the best fit line was found to be (R²=0.8, P=2×10^–8). Correlating the average estimated percentage change in contralateral limb vastus metabolic rate (⁠⁠) with the average experimentally observed percentage change in whole-body metabolic rate (⁠⁠), the best fit line was found to be (R²=0.9, P=0.05).

The trend in the sum of the metabolic rates of simulated muscles with increasing average exoskeleton torque was similar to the trend observed in measured whole-body metabolic rate (Fig. 5). Correlating the estimated percentage change in the sum of the metabolic rates of the simulated muscles (⁠⁠) with the experimentally observed percentage change in whole-body metabolic rate (⁠⁠), the best fit line was found to be (R²=0.6, P=7×10^–8). Correlating the average estimated percentage change in the sum of the metabolic rates of the simulated muscles (⁠⁠) with the average experimentally observed percentage change in whole-body metabolic rate (⁠⁠), the best fit line was found to be (R²=0.9, P=0.02).

Effort-related measures of the assisted soleus decreased with increasing net exoskeleton work. Average exoskeleton-side soleus muscle activation and fiber force during mid-stance increased as net exoskeleton work was increased (ANOVA, P=8×10^–3 and P=3×10^–3, respectively; Fig. 3A,B). Average exoskeleton-side soleus muscle activation and fiber force during late stance decreased by 66 and 73%, respectively, across exoskeleton work conditions (P=5×10^–6 and P=2×10^–6, respectively). Change in tendon length, from the instant the soleus muscle started lengthening to the instant it transitioned from lengthening to shortening, remained relatively unchanged as net exoskeleton work was increased (ANOVA, P=0.2; Fig. 3C). Soleus muscle fiber length, at the instant the soleus muscle transitioned from lengthening to shortening, and fiber contraction velocity at the instant of peak muscle fiber power, remained relatively constant across exoskeleton work conditions (ANOVA, P=0.06 and P=0.06, respectively; Fig. 3D,E). Muscle fiber work during late stance decreased by 77% across exoskeleton work conditions (P=8×10^–3; Fig. 3E). Similar trends were observed in the medial and lateral gastrocnemius muscles for a majority of these outcomes, but to a lesser extent (see Appendix, Figs A1, A2).

The positive work rate of the ankle exoskeleton plus the tendon increased by 137%, while the positive work rate of the soleus muscle decreased by 73% across exoskeleton work conditions (P=1×10^–5 and P=2×10^–4, respectively; Fig. 4). The positive work rate of the combined system increased by 72% across exoskeleton work conditions (P=1×10^–4).

Estimated metabolic rate of the exoskeleton-side soleus decreased with increasing net exoskeleton work (Fig. 5). Average activation/maintenance heat rate and positive mechanical work rate decreased by 30 and 72%, respectively, across exoskeleton work conditions (P=5×10^–4 and P=8×10^–3, respectively). Average shortening/lengthening heat rate remained relatively unchanged, while negative mechanical work rate decreased as net exoskeleton work was increased (ANOVA, P=0.3 and P=1×10^–3, respectively). Total estimated soleus metabolic rate decreased by 66% across exoskeleton work conditions (P=8×10^–3). Correlating the estimated percentage change in soleus metabolic rate (⁠⁠) with the experimentally observed percentage change in whole-body metabolic rate (⁠⁠), the best fit line was found to be (R²=0.4, P=8×10^–6). Correlating the average estimated percentage change in soleus metabolic rate (⁠⁠) with the average experimentally observed percentage change in whole-body metabolic rate (⁠⁠), the best fit line was found to be (R²=0.8, P=0.03).

Estimated total metabolic rate of the contralateral limb vastus decreased with increasing net exoskeleton work (ANOVA, P=9×10^–8; Fig. 5). Correlating the estimated percentage change in contralateral limb vastus metabolic rate (⁠⁠) with the experimentally observed percentage change in whole-body metabolic rate (⁠⁠), the best fit line was found to be (R²=0.5, P=2×10^–4). Correlating the average estimated percentage change in contralateral limb vastus metabolic rate (⁠⁠) with the average experimentally observed percentage change in whole-body metabolic rate (⁠⁠), the best fit line was found to be (R²=0.9, P=8×10^–3).

The sum of the metabolic rates of simulated muscles decreased with increasing net exoskeleton work (ANOVA, P=3×10^–8; Fig. 5). Correlating the estimated percentage change in the sum of the metabolic rates of the simulated muscles (⁠⁠) with the experimentally observed percentage change in whole-body metabolic energy consumption (⁠⁠), the best fit line was found to be (R²=0.5, P=1×10^–7). Correlating the average estimated percentage change in the sum of the metabolic rates of the simulated muscles (⁠⁠) with the average experimentally observed percentage change in whole-body metabolic rate (⁠⁠), the best fit line was found to be (R²=0.9, P=6×10^–3).

Providing increasing amounts of average exoskeleton torque, while maintaining zero net exoskeleton work, had both detrimental and beneficial effects on soleus muscle–tendon interactions. Normally, the soleus produces large forces throughout the dorsiflexion phase of stance, allowing the tendon to lengthen substantially and store mechanical energy. In this study, however, the exoskeleton displaced and reduced force in the soleus during early and mid-stance. This caused less stretch in the tendon and greater excursion of the muscle fibers than observed during unassisted walking. The decrease in tendon stretch had the detrimental effect of shifting work from the tendon to the muscle fibers. Reduced tendon stretch meant reduced elastic recoil during push-off, which was not adequately compensated for by the exoskeleton. The muscle fibers, therefore, did more work to maintain normal levels of total ankle positive work, but doing positive work with muscles is costly. The increase in muscle fiber excursion had complicated effects on the muscle's force-generating capacity. The muscle fibers operated closer to their optimal length at the time of peak power in late stance, which was beneficial to the muscle's ability to generate force. However, the muscle fibers also had to shorten a greater distance during push-off, thereby significantly increasing shortening velocity. Although the increase in fiber velocity helped to increase fiber power, force-generating capacity of muscle drops sharply with increased contraction velocity, thereby explaining the substantial reductions in soleus fiber force during late stance, despite much smaller reductions in activation.

Energy consumed by a muscle can be approximated by a combination of different heat and work rates (Hill, 1938; Mommaerts, 1969). As average exoskeleton torque increased, exoskeleton-side soleus activation/maintenance heat rate decreased, due to the forces applied in parallel with the soleus by the exoskeleton (Fig. 5). The shortening/lengthening heat rate, however, appeared to increase due to increased lengthening of muscle fibers during mid-stance and increased shortening during push-off. Exoskeleton-side soleus positive mechanical work rate increased with increasing average exoskeleton torque. Summing the activation/maintenance heat rate, shortening/lengthening heat rate, and net mechanical work rate together, soleus metabolic rate did not change significantly across conditions, but seemed to follow a similar trend to experimentally measured whole-body metabolic rate.

Reductions in the elastic recoil in the tendon, as well as increased lengthening and shortening of the soleus muscle fibers, negated the reductions in muscle activation and fiber force afforded by the exoskeleton. The human-robot system proved less efficient than the human system alone. Similar changes in soleus muscle–tendon mechanics were also observed during hopping with passive ankle exoskeletons (Farris et al., 2014); plantarflexor muscle fiber force decreased when passive assistance was provided, but fiber shortening velocity increased, resulting in no significant change in positive muscle fiber work. Such trade-offs in the observed changes in plantarflexor muscle–tendon mechanics led the metabolic rate of the plantarflexor muscles to remain relatively constant when hopping with and without assistance, which is comparable to our results for walking. Hopping with passive ankle exoskeleton assistance, however, still led to a reduction in whole-body energy consumption, likely due to off-loading of other muscle forces, particularly about the knee joint (Farris et al., 2014). Simulations of a simple, lumped model of the plantarflexor muscle–tendon units during walking with an elastic ankle exoskeleton also showed similar results: increasing exoskeleton stiffness decreased activation and force of the plantarflexor muscle fibers, but increased muscle fiber length changes and led to no change in work done by the muscle fibers (Sawicki and Khan, 2016). In contrast with our results, this simulation study showed that plantarflexor muscle metabolic rate decreased with increasing exoskeleton stiffness. This could be a result of differences in the way in which exoskeleton torque was applied in our study compared with Sawicki (2016), a result of different changes elsewhere in the body, or a result of different constraints on joint kinematics. In general, it seems the soleus muscle–tendon unit is sensitive to changes in operation. Slight alterations to the nominal system can have significant effects on coordination, which can be beneficial or detrimental to individual muscle and whole-body metabolic energy consumption, depending on the specific task.

As average exoskeleton torque was increased, changes in contralateral limb vastus mechanics and energetics were observed, which helps further explain the increase in experimentally measured whole-body metabolic rate. Changes in estimated contralateral limb vastus metabolic rate correlated well with experimentally observed changes in whole-body metabolic rate. Summing the metabolic rate of each muscle for which we had electromyographic data, we found that trends matched experimentally observed trends in whole-body metabolic rate well (Fig. 5).

Joint work is not necessarily a good predictor of muscle work and, consequently, energy consumed by a muscle. Positive exoskeleton-side soleus muscle fiber work increased with increasing average exoskeleton torque, but the biological ankle joint work remained relatively unchanged according to the muscle-generated ankle joint work computations, and actually decreased according to the inverse dynamics-derived ankle joint work computations.

Changing the amount of net work the exoskeleton provided also impacted exoskeleton-side soleus muscle mechanics and energetics, but in ways that were more expected. With increasing net exoskeleton work, peak exoskeleton-applied torque occurred later in stance, leading to significant changes in muscle–tendon dynamics. Reduced activation, in addition to reduced positive power during late stance, resulted in reduced effort of the soleus (Fig. 3A). Soleus muscle fiber force (Fig. 3B) and work (Fig. 3F) were reduced as net exoskeleton work was increased, thereby compromising the normal capabilities of the biological ankle. Positive work provided by the exoskeleton more than compensated for the reduced performance of the biological mechanisms, leading to an improved human-robot cooperative system. Metabolic rate of the exoskeleton-side soleus muscle significantly decreased with increasing net exoskeleton work, which accounted for a portion of the reduction in whole-body energy expenditure (Fig. 5).

As net exoskeleton work was increased, changes in contralateral limb vastus mechanics and energetics were observed, helping to further explain reductions in experimentally measured whole-body metabolic rate. Decreases in exoskeleton-side soleus metabolic rate were greater than those observed in the contralateral limb vastus, but both contributed to reductions in whole-body metabolic rate. Summing the metabolic rate of each muscle for which we had electromyographic data, trends fit experimentally observed reductions in whole-body metabolic rate well (Fig. 5).

Tendon stiffness and other muscle–tendon properties seem to be tuned such that the biological ankle joint operates efficiently. The results of this study support the idea that the physiological value of the Achilles tendon stiffness is optimal for muscle efficiency during walking and running (Lichtwark and Wilson, 2007). The lengthening and shortening of the Achilles tendon, instead of the muscle fibers, allows for energy to be stored and returned passively throughout stance. Positive work done by elastic elements can reduce the amount of positive work done by muscles.

Usefully interacting with biological muscles and tendons, via an external device, is complicated. Muscle–tendon mechanics are important and should be taken into account when designing devices to assist human motion. Adding an external device to the human body may affect muscle-level mechanics and energetics in unexpected ways. Disrupted muscle–tendon interactions were observed in this study and have similarly been observed in human hopping with ankle exoskeletons (Farris et al., 2013, 2014). Assistive devices should be designed and controlled to compensate for any compromised performance or functioning of muscles and tendons. Analyses similar to those discussed above can be used to help understand how different exoskeleton behaviors affect muscle-level mechanics, and provide insights into why certain device behaviors are more effective than others at assisting locomotion. For instance, torque support with a device can be an effective assistance strategy (Collins et al., 2015), but subtleties of how the external torques are applied and how the device interacts with the human musculoskeletal system greatly impact coordination patterns and overall effectiveness.

The modelling approaches used in this study can be applied to a wide array of human motions. The results suggest that, given a coordination pattern, via measured muscle activity and joint kinematics, it is possible to generate reasonable estimates of individual muscle mechanics and metabolic rate. In the future it may be possible to invert the process. Based on what we know about the mechanics and energetics of individual muscles, we can try to generate a set of desirable coordination patterns. It may even be possible to prescribe exoskeleton behaviors that elicit desirable changes in coordination.

Our modelling and simulation approach required making a number of assumptions and choices that need to be considered when evaluating the generated results. If the parameters used in the model were inaccurate, this could have led to invalid estimates of muscle mechanics and energetics. The parameters we used are, however, comparable to previously published work (Arnold et al., 2000, 2010) which is based on cadaver studies (Ward et al., 2009). Furthermore, to validate our approach, we compared muscle-generated ankle joint moments and powers to inverse dynamics-derived ankle joint moments and powers (Fig. 2). We optimized parameters to reduce the RMSE between the two, and performed an in-depth sensitivity analysis (Figs S1–S6) which shows the qualitative trends are robust to model parameters.

The soleus, lateral gastrocnemius and medial gastrocnemius were each modelled with a separate tendon, as opposed to one shared tendon. It is unclear which modelling choice is more appropriate for our study, but our fiber and tendon excursions were consistent with experimental ultrasound studies (Fukunaga et al., 2001; Lichtwark and Wilson, 2006; Cronin et al., 2010; Rubenson et al., 2012; Cronin et al., 2013). Moreover, qualitative trends in elastic element negative, positive and net work (Fig. 4) held for the combination of all plantarflexor tendons. Combined with the results of our sensitivity analysis, we are confident that this modelling choice does not affect our conclusions.

We were limited by the number of muscles we could measure experimentally. In particular, we did not measure muscle activity from the glutei, or other muscles acting about the hip, which are thought to consume a substantial amount of energy during walking. Nonetheless, the change in the sum of metabolic energy consumption from simulated muscles showed a similar trend to the change in whole-body metabolic energy consumption measured via indirect calorimetry; this independent validation increases our confidence in the primary findings of the study. Including more muscles in future experiments would make these analyses more complete.

Muscle-generated ankle joint mechanics did not perfectly match inverse dynamics-derived ankle joint mechanics, but most trends were consistent across the two methods. Results from inverse dynamics suggested that total exoskeleton-side positive ankle joint work decreased as average exoskeleton torque increased, while results from the electromyography-driven simulations suggested that total exoskeleton-side positive ankle joint work remained relatively unchanged. This inconsistency could have implications for our understanding of why contralateral limb knee mechanics and vastus metabolic rate were affected by torque applied at the exoskeleton-side ankle joint. We only optimized across those conditions with increasing average exoskeleton torque, but do not expect a better match would be obtained if we optimized across those conditions with increasing net exoskeleton work as well. There are inherent trade-offs that prevent errors across all conditions from simultaneously improving. These results illustrate the importance of knowing the limitations and assumptions inherent in a model and taking these into consideration when analysing and interpreting its outputs. To account for these limitations, we conducted sensitivity analyses and minimized inconsistencies between inverse dynamics-derived and muscle-generated joint mechanics by optimizing those model parameters in which we had the least confidence.

We simulated plantarflexor muscle–tendon mechanics and individual muscle energetics during walking with an ankle exoskeleton to gain a deeper understanding of how different exoskeleton assistance strategies affect the operation of the plantarflexor muscles and tendons. Providing increasing amounts of average plantarflexion torque with an ankle exoskeleton while providing no net work, disrupted soleus muscle–tendon interactions. Reduced tendon recoil was not sufficiently compensated for by the exoskeleton and this led to an increase in positive work done by the soleus muscle, which is costly. Providing increasing amounts of net exoskeleton work more than compensated for reduced work done by the soleus muscle–tendon unit, leading to a reduction in soleus force, work and total metabolic rate. Trends in the sum of the metabolic rates of the simulated muscles correlated well with trends in experimentally observed whole-body metabolic rate, suggesting that the mechanical and metabolic changes observed in the simulated muscles contributed to the measured changes in whole-body metabolic rate.

By performing these analyses we were able to explain experimentally observed changes in coordination patterns and metabolic energy consumption. Models without muscles and tendons would not have been able to capture these effects. Due to the sensitivity of muscle–tendon units to external disturbances, assisting locomotion by placing a device in parallel with muscles is challenging. When designing assistive devices, it is therefore important to consider how muscle–tendon mechanics might change due to interactions with the device and to ensure that the device sufficiently replaces any compromised function of the human musculoskeletal system.

The medial and lateral gastrocnemius muscles, which operate about the ankle joint, were also directly impacted by exoskeleton-applied assistance. These muscles, however, are biarticular, causing both plantarflexion of the ankle joint and flexion of the knee joint, and exhibit behaviors slightly different from the soleus muscle during normal walking. To deepen our understanding of how ankle exoskeletons affect those muscles involved in plantarflexion during assisted walking, we analyzed the changes in muscle-level mechanics of the medial and lateral gastrocnemius muscles under varying levels of exoskeleton assistance. Trends observed in the muscle–tendon mechanics of the medial and lateral gastrocnemius muscles were similar to those observed in the soleus muscle–tendon unit for most outcomes, but to a lesser extent (Figs A1, A2).

The authors thank Thomas Uchida for assistance with the OpenSim metabolics model.