Metabolic cost calculations of gait using musculoskeletal energy models, a comparison study

This paper compares predictions of metabolic energy expenditure in gait using seven metabolic energy expenditure models to assess their correlation with experimental data. Ground reaction forces, marker data, and pulmonary gas exchange data were recorded for six walking trials at combinations of two speeds, 0.8 m/s and 1.3 m/s, and three inclines, -8% (downhill), level, and 8% (uphill). The metabolic cost, calculated with the metabolic energy models was compared to the metabolic cost from the pulmonary gas exchange rates. A repeated measures correlation showed that all models correlated well with experimental data, with correlations of at least 0.9. The model by Bhargava et al. (J Biomech, 2004: 81-88) and the model by Lichtwark and Wilson (J Exp Biol, 2005: 2831-3843) had the highest correlation, 0.95. The model by Margaria (Int Z Angew Physiol Einschl Arbeitsphysiol, 1968: 339-351) predicted the increase in metabolic cost following a change in dynamics best in absolute terms.

The raw and processed experimental data are published as a dataset in [30].
Metabolic energy models 90 Seven metabolic energy models were selected for the current study: models BHAR04 [7], 91 HOUD06 [10], UMBE03 [9], LICH05 [21], MINE97 [8], MARG68 [23], and KIMR15 [24]. 92 Six models use muscle states (contractile element length, activation, stimulation) to 93 determine the energy rate of the individual muscles. Model KIMR15 calculates the 94 energy rate for each joint instead of each muscle, using the angular velocity and joint 95 moment. 96 The calculated metabolic cost of walking, C calc , is determined in J/kg/m as follows: 97 where T denotes the motion duration, m the participant's mass, v the speed, N mus the 98 number of muscles, andĖ i the energy rate of muscle i in W.

99
Models BHAR04, HOUD06, UMBE03, and LICH05 calculate the energy rate as a 100 function of work rate,ẇ, and heat rates, due to activation,ḣ a , maintenance of 101 contraction, andḣ m , and muscle shortening and lengthening,ḣ sl [7,9,10,21]: 102 E =ẇ +ḣ a +ḣ m +ḣ sl (2) The implementation of these models is detailed in S1 File. 103 Model MINE97 determines the energy rate for each muscle incorporating an 104 empirical function of the ratio between the contractile element velocity, v CE and the 105 maximum contractile element velocity, v CE(max) [8]: where φ = 0.054 + 0.506v CE + 2.46v 2 where a is the muscle activation, F max the maximum isometric force, andv CE is the shortening, and 120% efficient when lengthening [23]: Model KIMR15 does not use muscle states, but calculates the metabolic rate on the 111 joint level, using the joint moments and angular velocities. The metabolic rate is still 112 the sum of the heat rate and the work [24]: where the power, p i , at joint i is the product of the joint moment M and angular 114 velocityθ [24]: The heat rate is determined as follows [24]: whereḣ M = 0.054 is the heat rate for activation and maintenance,ḣ SL = 0.283 is the 117 shortening-lengthening heat rate for positive power, andḣ SL = 1.423 is the 118 shortening-lengthening heat rate for negative power, andq cc is the cocontraction heat 119 rate. The subscript max indicates the maximum over the gait cycle [24]. 120 Kinetic and Kinematic Data Processing 121 A two step approach was used calculate the joint angles, moments, and muscle states 122 and inputs necessary to determine the metabolic cost of walking level, uphill and 123 downhill using the metabolic energy models. 124 In the first step, the joint angles and moments were determined from marker and 125 ground reaction force data. The data was filtered backwards and forwards with a 126 second order Butterworth filter with a cut-off frequency of 6 Hz. Velocities and 127 accelerations were calculated using finite differences from marker positions. The data 128 was split into gait cycles and resampled to 100 data points per gait cycle.
averaged over all left and right gait cycles to find one average gait cycle.

135
In the second step, the muscle states (activation and contractile element length) and 136 stimulations were determined using the dynamic optimization method introduced in [32] 137 such that the joint moments generated by the muscles matched the moments found 138 using Winter's method. Fig 1 shows    The stimulations u(t), activations a(t), and contractile element lengths l CE (t) were July 25, 2019 6/26 found by solving the following dynamic optimization problem: were used: u(T ) = u(0), a(T ) = a(0), and l CE (T ) = l CE (0).

149
This dynamic optimization problem was solved using direct collocation, with 100 150 nodes for the gait cycle and a Backward Euler formulation. IPOPT 3.11.0 was used to 151 solve the optimization problem [33]. Muscle force patterns of level walking at 1.3 m/s 152 were compared to electromyography (EMG) data by Winter and Yack [34]. The EMG 153 data was digitized from the paper, and averaged when more than one muscle from a 154 muscle group was recorded (e.g. medial and lateral hamstrings). Finally, the muscle 155 state trajectories and stimulations, or the joint angular velocities and moments were 156 inserted in the seven metabolic energy models to find the calculated metabolic cost.
in W/kg was determined as follows for the resting and walking trials [35,36]: The resting trial was subtracted from each walking trial. The metabolic rate was 165 divided by walking speed to find the measured metabolic cost in J/kg/m: Analysis 167 First, the implementation of the metabolic energy models was verified and compared to 168 Miller [26]. To do so, the metabolic rate,Ė, was determined for the soleus for three 169 speeds (shortening at 1 l ce /s, isometric and lengthening at 1 l ce /s), and five activation 170 levels (0.05, 0.25, 0.5, 0.75 and 1). The stimulation was assumed to be the same as the 171 activation. The stimulation time, used in models BHAR04 and LICH05, was set to 1.

172
After the verification of the metabolic energy models, the metabolic cost was 173 calculated for every subject and trial (Eq. 1) and compared to the measured metabolic 174 cost (Eq. 2). Additionally, the calculated metabolic cost was determined for all joints 175 individually. The metabolic cost of biarticular muscles was split between the joints on 176 which they act using the ratios of the moment arms. The difference between the 177 calculated and measured metabolic cost, averaged over all subjects, is also presented.

178
Finally, we investigated the ability of the metabolic energy models to predict a 179 change in energy following a change in dynamics -due to altered walking speed and/or 180  Metabolic rate of soleus, calculated with the metabolic energy models for different conditions. Metabolic rate as calculated by the metabolic energy expenditure models that are based on individual muscles of the soleus muscle at four activation levels, for an isometric condition, and a shortening and lengthening velocity of one optimal fiber length per second. muscle force graph shows EMG data for comparison [34], scaled to the maximum muscle 204

Joint Kinetics and Kinematics
force. The muscle force pattern matches the EMG data. In the gastrocnemius, the EMG activation increases earlier in the gait cycle than the muscle force, while in the 206 rectus femoris, there is an extra burst in the muscle force at about 60-80% of the gait 207 cycle, and in the tibialis anterior, the activity of the EMG is generally higher than the 208 muscle force.   One measurement was missing for the downhill trial at 1.3 m/s due to a malfunction of 224    while model BHAR04 and model LICH05 correlated best with the experimental data 259 (see Table 2). The regression model of model MARG68 had a slope closest to unity, 260 which indicates that model MARG68 calculates absolute differences between two trials 261 most accurately. All models had a slope larger than one, meaning that the trials with 262 higher metabolic demands were underestimated to a larger extent than those with a 263 smaller metabolic demand, which is also visible in Table 1 cost. A possible reason is the approach used to find the muscle activations by minimizing muscular effort. Then, muscular co-contraction was disregarded. However, 278 the muscle forces matched EMG patterns (see Fig 3), which supports that the muscle activation pattern is accurate, but the quality of the activation level cannot be assessed. 280 Additionally, the model used in the current study was a two-dimensional model of the 281 trunk and lower leg, meaning that metabolic cost due to lateral and rotational motion, 282 as well as arm swing, were not taken into account, which could cause an 283 underestimation of metabolic cost. Finally, the calculated metabolic cost could be lower 284 since negative work was subtracted in the current study. However, the high correlation 285 indicates that the mechanical analysis still yielded a correct prediction of increases or 286 decreases in metabolic cost, which was the main aim of the current study.

287
The calculated metabolic cost in the current study is lower than calculated by resting metabolic rate, and used a three dimensional model [7,9], while the current 291 study used a sagittal plane model with eight muscles and no resting metabolic rate.

292
When accounting for a resting metabolic rate of 1. and much lower than reported by Bhargava et al. [7] and Umberger et al. [9]. Therefore, 299 the underestimation was not unexpected. While the dimensionality of the model played 300 a role, the difference in metabolic cost between the original and more recent work could 301 also be caused by the limited number of parameters that were used for the muscle 302 stimulation or due to differences in kinetics or kinematics, which were not reported by 303 either [7,9] both of which are discussed later in this section. Therefore, the absolute performance of 318 the metabolic energy models should be investigated further.

319
The handling of muscular energy expenditure during lengthening is still debated [26]. 320 During lengthening, the energy rate can be negative, which is physically impossible.

321
However, the negative work should be subtracted from the metabolic cost in models 322 BHAR04, HOUD06, UMBE03, and LICH05 [7,9,10,21]. We aimed to see if predictions 323 improved without subtracting negative work, which is physically more sensible. between 27% and 34%, the hip between 23% and 39%, and the ankle between 29% and 362 50%, while the absolute metabolic cost varied greatly between models [26]. Therefore, 363 when creating simulation of walking, some energy minimization is required to distribute 364 energy between joints more accurately.

365
Similar trends existed between the metabolic energy models in the breakdown of the 366 metabolic cost per joint, though differences exist in the breakdown for the specific trials 367 (see Fig 5). Since no measurements were taken, it cannot be said which model is more 368 correct. However, the trends with speed and slope were very similar between models.

369
The hip was mainly responsible for the change of energy expended with the slopes.

370
When walking uphill, the energy expended in the hip increased, while it decreased when 371 walking downhill. Between the speeds, the most obvious difference was an increase in energy expended in the ankle for the larger speed, while less energy is expended in the 373