Predictive Simulations of Gait with Exoskeletons that Alter Energetics

Robotic exoskeletons, designed to augment human locomotion, have the potential to restore function in those with mobility impairments and enhance it in able-bodied individuals. However, optimally controlling these devices, to work in concert with complex and diverse human users, is a challenge. Accurate model simulations of the interaction between exoskeletons and walking humans may expedite the design process and improve control. Here, we use predictive gait simulations to investigate the effect of an exoskeleton that alters the energetic consequences of walking. To validate our approach, we re-created an past experimental paradigm where robotic exoskeletons were used to shift people’s energetically optimal step frequency to frequencies higher and lower than normally preferred. To match the experimental controller, we modelled a knee-worn exoskeleton that applied resistive torques that were either proportional or inversely proportional to step frequency—decreasing or increasing the energy optimal step frequency, respectively. We were able to replicate the experiment, finding higher and lower optimal step frequencies than in natural walking under each respective condition. Our simulated resistive torques and objective landscapes resembled the measured experimental resistive torque and energy landscapes. Individual muscle energetics revealed distinct coordination strategies consistent with each exoskeleton controller condition. Predicted step frequency and energetic outcomes were best achieved by increasing the number of virtual participants (varying whole-body anthropometrics), rather than number of muscle parameter sets (varying muscle anthropometrics). In future, our approach can be used to design controllers in advance of human testing, to help identify reasonable solution spaces or tailor design to individual users.


26
This decade has seen rapid progress in the development of assistive devices designed to augment human 27 locomotion. As early as the 1960s the first robotic exoskeleton prototypes, which are anthropometric in 28 form and are designed to augment the user's movements, were developed [1,2]. However, early devices were 29 often too heavy and cumbersome for practical use. Recent advances in both robotic hardware and software 30 have pushed back technological boundaries, allowing for the development of more sophisticated, light weight, 31 and intelligent designs [1][2][3]. Currently, a wide range of lower-exoskeletons are being tested in research labs 32 around the world and even being sold commercially. Exoskeletons that augment the ankle, knee, hip, or 33 some combination of joints exist. Some exoskeletons are active, or powered, using an actuator to apply 34 torque to a joint, while others are passive, relying on springs or dampers to store and return energy during 35 outcome measures can include optimal device design or control parameters, the resulting individual muscle 90 activations or energy expenditure, and the changes to muscle generated joint torques. For example, Van den 91 Bogert (2003) used this approach to solve for the optimal geometric configuration of a lower-limb passive 92 exotendon system that provided the greatest reduction in muscle generated joint torques [22]. However, 93 an experimental study of this concept showed that while muscle generated joint torques decreased, joint 94 angles changed and metabolic cost of walking was not reduced [25]. Another group has comprehensively 95 investigated various actuated assistive device designs to reduce the metabolic cost of loaded walking and 96 running [23,24]. Their simulations demonstrated that frontal plane hip assistance may lead to the greatest 97 improvements in economy and that ideal device torques can differ significantly from muscle generated torques 98 in natural walking. However, while their simulations predicted large reductions in metabolic cost this was 99 once again not reproduced in accompanying experiments [26]. The current discrepancies between inverse 100 dynamic simulations and real-world human device testing may be because this approach constrains total 101 joint torques and kinematics, failing to account for altered limb dynamics and the adaptive strategies of the 102 user. 103 Forward dynamic simulations allow joint kinematics and muscle generated joint torques to change and 104 may therefore overcome some of the limitations of inverse dynamic simulations when investigating exoskele-105 tons. Furthermore, this approach allows for simulations of movements for which no experimental data is 106 available, so-called predictive simulations. In this class of simulations, a high-level trajectory optimiza-107 tion problem is solved, where a periodic gait cycle is found that minimizes some objective function (or 108 goal) [27][28][29][30]. The objective function often contains a number of weighted terms, such as maximizing smooth-109 ness, minimizing torques, and minimizing an energy term, be it muscular effort [27], metabolic energy [28] 110 or both [29]. While these simulations can in principle be solved without requiring any experimental data, 111 musculoskeletal model inaccuracies currently prevent many of these simulations from producing realistic gait 112 without the use of reference data or expert input. Some researchers include a tracking term that weights 113 kinematics and kinetics that are similar to references gait profiles [30], while others have expertly hand-tuned 114 the objective function to produce a similar result [29]. Importantly though, these approaches do not fully 115 constrain kinematics and kinetics, allowing adaptive strategies to emerge in response to the exoskeleton. 116 Methods to quickly and accurately solve predictive simulations of musculoskeletal models have only been 117 developed in the last ten years [27], meaning their application to exoskeleton design and control has not 118 been fully realized. Predictive simulations have been used to study other aspects of gait including, crouch 119 gait in cerebral palsy [31], joint contact forces in running [32], shoe design on running efficiency [17], and 120 loading asymmetry in prostheses [30]. A few groups have used predictive simulations investigate the optimal 121 stiffness of passive ankle exoskeletons [18,33]. However, these groups greatly simplified the musculoskeletal 122 models, either limiting actuation to hip torques alone [33] or by not modeling tendon behavior [18]. These 123 studies demonstrate that optimal gait kinematics and kinetics are likely to deviate from natural walking 124 patterns in the presence of an assistive aid-both studies found this to be the case. However, accompanying 125 experimental studies were not performed to validate these model predictions.

126
An existing limitation of predictive simulations is that they are often generated for a single set of mus-127 culoskeletal model parameters, limiting generalizability [18,[27][28][29][30]33]. This is equivalent to generating 128 predictions for a single participant, meaning conclusions may be dependent on the particular choice of 129 anthropometric parameters. Recently, different approaches have been used to try to improve the generaliz-130 ability of predictive simulations. Milller et al. (2013) produced simulations for four different participants, 131 representing younger and older males and females [32]. Dorschky et al. (2019) pioneered the use of a set of 132 'virtual participants', where height and weight are randomly drawn from a reference set of participant data 133 and, for each virtual participant, various sets of muscle parameters are randomly drawn from distributions 134 that capture the variation expected in the population [17]. However, solving a large number of trajectory 135 optimization problems in computationally costly and time consuming, and it is currently unclear how many 136 virtual participants are required to make generalizable conclusions.

137
Our purpose in this study was to use predictive simulations of gait to investigate the effect of an ex-138 oskeleton that alters the energetic consequences of walking, using a virtual participant design. Specifically, 139 we aimed to re-create a past experimental paradigm [4,34], where robotic exoskeletons were used to shift 140 people's energetically optimal step frequency to frequencies higher and lower than normally preferred. In 141 the human experiments, participants adapted their step frequency to converge on the new energetic optima 142 within minutes and in response to relatively small savings in cost. To match the controller in this real-world 143 experiment, we modeled a knee-worn exoskeleton that applied resistive torques that were either proportional or inversely proportional to step frequency-decreasing or increasing the energy optimal step frequency, re-145 spectively. We then performed predictive simulations of human walking with the device for a set of virtual 146 participants, which we created to have similar anthropometric variables to our experimental participants. 147 We compared our predictions to the experimental data to test if we could: i. replicate the exoskeleton torque 148 applied at the knee joint during gait, ii. produce objective landscapes with optima at high and low step 149 frequencies, and iii. solve for optimal gaits through adaptations in step frequency. We also investigated 150 the sensitivity of our results to the sample size of virtual participants used and examined individual muscle 151 energetics, offering insight into distinct coordination strategies adopted in response to the exoskeleton. properties. We calculated individual muscle metabolic rate using the model by Margaria [35]. Ground 162 contact was modelled with a penetration-based model [30]. The multibody dynamics and ground contact

166
We then added an ideal, massless, exoskeleton to our simulations by applying torques at the knee joints that 167 resisted both knee flexion and extension, and thereby added an energetic penalty. We applied two different 168 exoskeleton controller conditions: penalize-high, where higher step frequencies were penalized such that the 169 optimal step frequency was lower than natural, and penalize-low, where lower step frequencies were penalized 170 such that the optimal step frequency was higher than natural. We designed these controllers to replicate our 171 real-world exoskeleton controllers [4].

172
To accomplish this, we made the applied torque (M exo ) proportional to step frequency (s) throughout the stride and proportional to knee angular velocity (q knee ) within the stride. To ensure that the exoskeleton torque was penalizing, we determined its absolute value and multiplied this by the opposite sign of the knee angular acceleration (q knee ). The absolute peak torque was limited to 12 Nm. This yields the following equation: where t is time, c T = 3.36x10 −2 Nm/A is the motor's torque constant and c q = 60 C/rad is a constant that 173 represents the relationship between motor current (A) and knee angular velocity (rad/s). The zero-torque 174 step frequency (s 0 ) is set to -15% of the natural step frequency (s nat ) in the penalize-high condition and 175 +15% of the natural step frequency in the penalize-low condition.

176
To allow the optimization to be solved with a gradient-based algorithm, we converted Equation 1 to a 177 twice-differentiable function. This involves three steps.

178
First, the absolute value of the applied torque (M exo,abs ) was determined using the following soft absolute  Second, the applied torque was limited to an absolute peak of 12 Nm (M exo,lim ) using an exponential 181 soft-max function: Third, the sign function was approximated using the inverse tangent: where we used the factor 10 to increase the nonlinearity and similarity to the sign function. 184

185
We generated muscle-driven simulations of walking at 1.3 m/s by solving trajectory optimization problems. 186 The objective was to minimize: i) a weighted sum of muscular effort, which was calculated as the cubed muscle 187 stimulation; ii) a tracking error, which was calculated as the difference between predicted and literature joint 188 angles and ground reaction forces during normal walking [36]; and iii) a regularization term [28,37]. The 189 weighting ratio between effort and tracking error was 1000:1, a factor one hundred times higher than previous 190 work with the same model [30]. We chose this high ratio to allow simulations to deviate from tracked data as 191 much as possible, and therefore adapt step frequency, while still maintaining realistic gaits. We included a 192 small regularization term to minimize the derivatives of the states and controls, which enhances convergence 193 of the optimization without affecting simulation accuracy [28,37]. We assumed symmetry and simulated 194 only half a gait cycle. This yields the following optimization problem: where T is the duration of the gait cycle, u i is the stimulation of muscle i, q j is the angle for joint j, q j,trk  We created walking simulations for three conditions: natural, where no added knee torques were applied 207 by the exoskeleton, as well as penalize-high and penalize-low, where each respective exoskeleton controller 208 applied torques to the knees. We added the applied exoskeleton torques to the torques generated by the 209 muscle forces. For each condition, we first optimized step frequency to investigate if we could predict the 210 energy optimal gait adaptation. This optimal step frequency is the step frequency where the full objective 211 (Equation 5) is minimized. Next, we created simulations at a range of fixed step frequencies (-15% to +15% 212 of the natural optimal step frequency, at increments of 1%) to generate the landscapes for all three conditions 213 (natural, penalize-high and penalize-low). We created a full objective landscape, which is the full objective 214 function evaluated for the range of fixed step frequencies, an effort term landscape, which is the effort term 215 of the objective function alone evaluated for the range of fixed step frequencies, and the metabolic rate 216 landscape, which is a sum of the individual muscle metabolic rates (to estimate whole-body metabolic rate) 217 for the range of fixed step frequencies. The step frequency was fixed by constraining the duration of the gait To create all simulations, we solved trajectory optimization problems with direct collocation and a back-220 ward Euler discretization. We used 30 nodes per half gait cycle. We coded the problems in MATLAB 221 R2018a (Mathworks, NA, USA) and used IPOPT [38] to solve the resulting large-scale constrained opti-222 mization problem. We first created a simulation of standing by constraining the derivatives of the degrees 223 of freedom to be equal to zero:ẋ = f (x, u) = 0 while minimizing muscular effort, as described in [37]. 224 We repeated this problem with 50 different random initial guesses and used the solution with the lowest 225 objective. Then, we used standing as an initial guess to find walking simulations with free, or optimal, step 226 frequencies in the natural, penalize-high, and penalize-low conditions. Next, to create the landscapes, for 227 each condition we first solve for a simulation where the step frequency is equal to the optimal step frequency 228 in the natural condition. To find these simulations for each condition, we use the previously solved walking 229 simulations with free, or optimal, step frequencies for the same condition as an initial guess. Next, for each 230 condition, we used the simulation with step frequency fixed to the natural optimum as an initial guess for 231 simulations fixed at +1% and -1% of the natural optimum. For each condition, we then used these solved 232 simulations as initial guesses to find the simulations fixed at +2% and -2% of the natural optimum, and 233 we repeated this process until we reached +15% and -15% of the natural optimum, the full range of step We first identified and removed all simulations of virtual participant and muscle parameter set combinations 255 that produced at least one unrealistic simulation. We identified simulations as unrealistic when the difference 256 in metabolic rate between at least one simulation with fixed step frequency and free step frequency in the 257 same condition was larger than 2 W/kg. We further verified these outliers by visual inspection of the gait 258 cycle, to confirm that the gait was unrealistic and likely the result of a local minimum. For all simulated 259 outcome measures, we then averaged across the remaining muscle parameter variations within each virtual 260 participant. We then calculated the mean and standard deviation of the outcome measures across all 50 261 virtual participants. We compared the average height and body mass of our virtual and experimental 262 participants to ensure they were similar.

263
To investigate the quality of predictions from our simulations, we compared the simulation results to 264 the experimental results of   [4]. Simon Fraser University's Office of Research Ethics 265 approved the protocol, and participants gave their written, informed consent before experimentation. We 266 determined the coefficient of determination (R 2 ), between simulation predictions and experimental data, for 267 the across-stride average and within-stride torque profiles in the penalize-high and penalize-low conditions, 268 as well as for the objective, effort, and energy landscapes in all three conditions. To compare the optimal, 269 or freely chosen, step frequencies and the accompanying full objective, effort term, or metabolic rate values 270 between the experimental and simulated data, we calculated the percent changes from the natural condition 271 to the penalize-high and penalize-low conditions. We then compared the means and standard deviations of 272 theses percent changes between the experimental and simulated data. 273 We also investigated the changes in simulated metabolic rate for individual muscles across the three 274 conditions. For all three conditions, we generated and compared individual muscle metabolic rate landscapes, 275 which were the individual muscle metabolic rates computed for the range of fixed step frequencies. For 276 simulations with an optimal, or freely chosen step frequency, we also determined the percent change in 277 individual muscle metabolic rate between the natural and the penalize-high and the penalize-low conditions, 278 during both stance and swing phase. To split individual muscle metabolic rates throughout the stride into 279 stance and swing, we used the vertical ground reaction force to determine when the foot was in contact with 280 the ground. To calculate the relative contribution of stance and swing to metabolic rate, we divided by the 281 full duration of the gait cycle, such that the sum of both was equal to the total metabolic rate.

282
Finally, we investigated how different virtual participant and muscle parameter set sample sizes affected 283 our simulation outcomes. We repeated our previously described analysis to determine the percent change in 284 optimal, or freely chosen, step frequencies between the natural condition and the penalize-high and penalize-285 low condition for virtual participant sample sizes between 1 and 50 and muscle parameter sample sizes 286 between 1 and 10. We also calculated the accompanying percent change in whole-body metabolic rate. To 287 accomplish this, we start with the full set of 50 virtual participants by 10 muscle parameters sets. For each of 288 the desired number of virtual participants and muscle parameter sets, we randomly select samples from the 289 full set. We then repeat this sampling 10 times, to account for randomness in the data, and then calculate 290 the percent difference between the outcome measures of the sample and the full dataset. For each outcome 291 measure, we then calculate the root mean square (RMS) error across these 10 repeats for each combination 292 of sample sizes. 293 Figure 2: Comparison of simulated and experimental within-stride and across-stride torques. Average torques throughout the gait cycle for the penalize-high condition (A) and penalize-low condition (B) at -10%, 0% and +10% of natural step frequency. C. Average torques across the stride for the penalize-high (blue) and penalize-low (red) conditions.

294
Of the 500 combinations of virtual participants and muscle parameter sets (50 x 10), we removed 6 that  Our simulated exoskeleton torques within the stride and averaged across the stride, resembled that from 304 human experiments (Fig. 2). The coefficients of determination between the simulated and experimental 305 within-stride torques ranged between 0.33 and 0.49. For both the penalize-high and penalize-low controllers, 306 our simulated within-stride torques better matched experimental when angular velocities were greater, dur-307 ing late stance and swing (>50% of gait cycle, Fig. 2AB). During early stance within-stride estimates were 308 poorer, though still exhibited two moderate peaks. Although within-stride differences existed in the torque 309 profiles, the average simulated torques across strides were very consistent with that applied in human ex-310 periments (Fig. 2C). The coefficients of determination between the simulated and experimental across-stride 311 torques were 0.97 and 0.98 for the penalize-high and penalize-low conditions, respectively. The penalize-low 312 and penalize-high conditions exhibited the desired relationships with step frequency: being negatively and 313 positively sloped, respectfully, and delivering the same magnitude torque at the 0% step frequency. The torques applied by our simulated exoskeleton controllers produced full objective, effort term, and 315 metabolic rate landscapes similar in shape to human experimental energy landscapes (Fig. 3), though differ-316 ences existed in the non-normalized natural optimal step frequencies and metabolic rates. In all landscapes, 317 the simulated penalize-low controller produced a negatively sloped gradient about the initial preferred step 318 frequency (0%) and shifted the optimal gaits to higher step frequencies, as was the case in the experimental 319 energy landscapes. The simulated penalize-high controller produced a positively sloped gradient about the 320 initial preferred step frequency (0%) and shifted the simulated energy optimal gaits to lower step frequen-321 cies, once again in a manner consistent with the experimental landscapes. However, in the full objective 322 landscape the optimal step frequency occurred at 2.24 Hz, which is higher than the step frequency of the 323 tracked data and the experimental optima, which were both 1.8 Hz [4,36]. The effort term and metabolic rate 324 landscapes in the natural condition had optimal step frequencies that were both lower than that for the full 325 objective landscape (2.14 Hz and 2.16 Hz, respectively), but were still high. Moreover, the simulated energy 326 landscapes predicted metabolic rates that were notably higher than the experimental metabolic rates (+0.5 327 to 1.5 W/kg), as well as smaller relative changes in metabolic rate under the penalize-high and penalize-low 328 conditions. Overall, the coefficients of determination, between the simulated and experimental landscapes, 329 ranged between 0.62 and 0.98 and coefficients tended to be higher for the penalize-high and penalize-low 330 conditions than the natural condition.

331
Like our human participants, our virtual participants responded to the new landscapes, adapting toward 332 the energy optimal step frequencies (Fig. 3). Under the penalize-high and penalize-low conditions virtual 333 participants adapted their preferred step frequency by -5.7% ± 1.5% and 8.1% ± 3.1%, respectively (these 334 step frequencies are those that minimize the full objective). These adaptation magnitudes were both indis-335 tinguishable from those displayed by experimental participants in each respective condition (-5.7% ± 3.9%, p 336 =0.79; -6.9% ± 4.3%, p =0.16). In the penalize-high condition, virtual participants' full objective magnitude 337 at the new preferred step frequency was 1.9% ± 1.7% lower than the full objective magnitude at the initial 338 preferred step frequency under the penalize-high control function (Fig. 3A). For the effort term the equiv-339 alent reduction was 3.6% ± 2.4% (Fig. 3B), while for the simulated metabolic rate the reduction was 2.5% 340 ± 0.91% (Fig. 3C). These reductions are smaller than the 8.1% ± 7.0% reductions in metabolic rate seen in 341 the experiment under the penalize-high condition (Fig. 3C, [4]). In the penalize-low condition, virtual par-342 ticipants' full objective magnitude at the new preferred step frequency was 2.2% ± 2.1% lower than the full 343 objective magnitude at the initial preferred step frequency under the penalize-low control function (Fig. 3A). 344 For the effort term the equivalent reduction was 2.5% ± 2.5% (Fig. 3B), while for the simulated metabolic 345 rate the reduction was 1.4% ± 1.0% (Fig. 3C). These reductions are again smaller than the 4.0% ± 3.8% 346 reductions in metabolic rate seen in the experiment under the penalize-low condition (Fig. 3C, [4]). Similar 347 to the experiment, the simulation predicts that the changes in effort and metabolic rate in the penalize-low 348 condition are roughly double those changes in the penalize-high condition, but this was not observed for the 349 full objective.

351
Individual muscle changes in metabolic rate, across the energy landscapes, offer insight how the exoskeleton 352 controllers produce changes in whole-body metabolic rate. Across the full landscape (-15% to +15% change 353 in step frequency) the range of metabolic rates (difference between the minimum and maximum rate across 354 penalize-high, penalize-low, and natural) are similar, about 0.3 W/kg, for all muscles except for the vastus, 355 where the range is 0.6 W/kg (Fig. 4). This indicates that the vastus has the largest influence on the whole-356 body energy landscape, while all other muscles have a similar, but lower, influence. We also found that 357 individual muscle changes in metabolic rate across the energy landscape differed between muscles that cross 358 the knee and those that do not. Under natural conditions, muscles that cross that knee (hamstrings, rectus 359 femoris, vastus, and gastrocnemius) tend to display an optimum (minimum metabolic rate) within the range 360 of step frequencies in the landscape, while those muscles that do not cross the knee (iliopsoas, gluteals, 361 soleus, and tibialis anterior) tend to have a consistently increasing or decreasing slope across the landscape. 362 Although the steepness of these slopes is altered by the exoskeleton controller conditions, the direction 363 remains unchanged (i.e iliopsoas has a positively sloped gradient for natural, penalize-high and penalize-low 364 conditions). Conversely, for muscles that cross the knee, the two exoskeleton controllers create slopes of 365 opposite direction-creating a positively sloped gradient under the penalize-high condition (decreasing the 366 optimal step frequency) and a negatively sloped gradient under the penalize-low condition (increasing the 367 optimal step frequency). For these muscles, the penalize-high condition tends to create steeper slopes than 368 the penalize-low condition, the effect of which is evident in the whole-body metabolic rate landscapes.

369
An examination of simulated metabolic rate at the individual muscle level during the stance and swing 370 phase reveals distinct coordination strategies consistent with each exoskeleton controller condition (Fig. 5). 371 Our simulated optimal gaits show changes in muscle metabolic rate not only for muscles spanning the knee, 372 but also muscles crossing the hip or ankle. In the penalize-high condition, muscle metabolic rate decreased 373 across more swing muscles, but showed large increases for muscles crossing the ankle during stance (i.e. 374 Figure 3: Comparison of simulated full objective, effort term, and metabolic rate landscapes with experimental metabolic rate landscapes. Experimental metabolic rate landscapes as well as simulated full objective (A), effort term (B), and metabolic rate (C) landscapes for the natural (grey), penalize-high (blue), and penalize-low (red) conditions. In all plots, solid lines are simulated data, while dashed lines are experimental data [4]. Squares indicate optimal gaits, that minimize the full objective, for each simulated condition, with error bars representing the standard error, while the circles indicate experimental optimal gaits for each condition. In C, shading represents standard error.
iliopsoas, swing: -9%; tibialis anterior, stance: +15%). This is consistent with an adaptation toward lower 375 step frequencies, requiring more work during push off at the ankle to lengthen steps and relatively less work 376 Figure 4: Simulated individual muscle contributions to energy landscapes. Simulated muscle metabolic rate across step frequency for the penalize-high (blue), penalize-low (red) and natural (black) conditions. In all plots, solid lines show the average metabolic rate, across 50 virtual participants, while the faded fills show the standard error. All muscle y-axes span the same magnitude range (0.03 W/kg), accept the vastus, which spans 0.6 W/kg. to swing the limb. In the penalize-low condition, the opposite was the case; muscle metabolic rate increased 377 across nearly all muscles during swing, but showed decreases for many muscles during stance, particularly 378 the tibialis anterior that crosses the ankle (i.e. gluteals, swing: +35%; tibialis anterior, stance: -21%). This 379 is consistent with an adaptation toward higher step frequencies, requiring more work at the hip during swing 380 and relatively less work during a shorter stance phase. 381 Figure 5: Effects of exoskeleton on simulated individual muscle metabolic rate. Average muscle metabolic rate, during stance (left, standing leg) and swing phase (right, swinging leg), in W/kg for the natural condition (A), and relative to natural for the penalize-high (B) and penalize-low (C) conditions.

382
We found that both simulated step frequency adaptations and metabolic rate were primarily affected by the 383 number of virtual participants, although increasing the number of muscle parameters sets had an effect at 384 larger virtual participant numbers (Fig. 6). The RMS difference, between outcomes with a limited dataset 385 and that with the full dataset (50 virtual participants and 10 muscle parameter sets), clearly increases when 386 the number of virtual participants is decreased, while there is a smaller effect when the number of muscle 387 parameter sets is decreased, for both step frequency (Fig. 6A) and metabolic rate (Fig. 6B). For example, 388 when using only one virtual participant and one set of muscle parameters, the RMS difference in step 389 frequency is 1.9% in the penalize-high condition and 3.8% in the penalize-low condition. These differences 390 are high, being roughly half the predicted step frequency adaptation magnitudes under the full dataset.

391
However, when using all 50 virtual participants and one set of muscle parameters, the RMS difference 392 decreases to less than 0.2% in the penalize-high condition and less than 0.4% in the penalize-low condition 393 (a ten-fold reduction). Conversely, if only a single virtual participant is used and all 10 sets of muscle 394 parameters are included, the RMS difference remains high, at 1.5% in the penalize-high condition and 3.0%  We used predictive simulations of gait to investigate the effect of an exoskeleton that alters the energetic 404 consequences of walking using a virtual participant design. We were able to generate predictive simulations 405 that well-matched experimental results [4]. While simulated within-stride torques tended to differ from ex-406 perimental, particularly during stance, our simulated torques averaged across the stride were very similar to 407 experimental. This resulted in full objective, effort term, and metabolic rate landscapes with optima shifted 408 to lower step frequencies under the penalize-high controller condition and higher step frequencies under the 409 penalize-low controller condition, as desired. Our simulated optimal gaits, under each condition, displayed 410 step frequency adaptations consistent with these shifts in optima and indistinguishable from our past exper-411 imental results. Simulated individual muscle metabolic rates provided insight, beyond that available from 412 experimental data, into what drives whole-body changes metabolic rate and the new optima. In particular, 413 the slope of the individual muscle energy landscapes change direction under the differing controllers only for 414 muscles that cross the knee. Furthermore, individual muscle stance and swing costs at the optimal gaits re-415 veal distinct coordination strategies consistent with adaptations under each exoskeleton controller condition. 416 Finally, our virtual participant study showed that increasing the number of virtual participants improved 417 simulated outcomes much more than increasing the number of muscle parameter sets.

418
Our use of predictive simulations is not without its limitations. While trends in metabolic rate and 419 step frequency adaptations were similar between simulation and experiment, magnitudes differed. Across 420 all landscapes, our model predicted metabolic rate was often 0.5-1.5 W/kg higher than that measured 421 experimentally. This is not unexpected. Various metabolic models often produce reliable and reproducible 422 relative changes in metabolic rate, but magnitudes can vary widely between models [40,41]. In particular, 423 the commonly used Margaria metabolic model we implemented [35] has been shown to produce estimates 424 that tend to be higher than other models [40]. We also found that the step frequency that minimized the 425 full objective under the natural condition was higher than that observed experimentally (2.24 Hz vs. 1.8 Hz, 426 respectively). Interestingly, the step frequencies that minimized effort (2.14 Hz) and metabolic rate (2.16 Hz) 427 alone were more closely aligned with the experimental, indicating that the tracking term in the full objective 428 favored higher step frequencies. This was surprising because the step frequency of the tracked data (1.8 Hz) 429 matched the experimental. Although removing the tracking term entirely would cause gait simulations to 430 become more inaccurate [28], we set the weight of the tracking term to be as small as possible to minimize 431 this effect. We also explored removing the hip angle from the tracked variables (which we expect to be most 432 correlated with step frequency) and including a term to minimize metabolic rate instead of effort in the 433 objective using the model described in [28]. However, in both cases we found that the step frequency that 434 minimized the full objective remained higher than the step frequency that minimized the effort or metabolic 435 rate term.

436
Simulated individual muscle energetics can offer insight into the full-body energetics and the resulting 437 gait adaptations. First, individual muscle energy landscapes revealed that the slope of the energetic gradient 438 changed direction under the two controller conditions only for muscles that cross the knee. This indicates 439 that the optimal step-frequencies in the penalize-high and penalize-low conditions are largely driven by the 440 energetics of these knee-crossing muscles. These muscles favor a change in step frequency towards a lower 441 resistive torque, while the hip-and ankle-crossing muscles favor a higher or low step frequency independent 442 of the controller condition. However, this is not to say that muscles crossing the hip or ankle do not adapt 443 activity and display changes in metabolic rate. While their slope direction remains largely unchanged, their 444 slope steepness, magnitude and timing of expenditure during gait phases do change under the controller 445 conditions. This indicates that exoskeleton torques are not simply counteracted or offset, but rather the 446 predictive simulations solve for complex and adaptive lower limb coordination strategies. In particular, 447 we observed changes in individual muscle energetics consistent with: i) an adaptation toward higher step 448 frequencies, requiring more work at the hip during swing and relatively less work during a shorter stance 449 phase and ii) an adaptation toward lower step frequencies, requiring more work during push off at the ankle to 450 lengthen steps and relatively less work to swing the limb [42]. Finally, individual muscle metabolic rates may 451 explain why we find a steeper whole-body energetic gradient in the penalize-high condition compared to the 452 penalize-low condition, in both simulation and experiment. It is muscles that cross the knee that drive these 453 steeper slopes; the individual muscle energetic gradients are steeper under the penalize-high than penalize-454 low condition. This appears to be because under the penalize-high condition, the highest resistive torques 455 are applied during high step frequencies, when knee angular velocity is relatively higher. Conversely, under 456 the penalize-low condition, the highest resistive torques are applied during low step frequencies, when knee 457 angular velocity is relatively lower. High resistive torques applied to muscles moving at higher velocities in 458 turn produce greater muscle work and therefore greater changes in metabolic rate. This leads to a higher cost 459 at +10% step frequency under the penalize-high condition than at -10% step frequency under the penalize-460 low condition. Another possible contributor is that exoskeleton resistive torques are proportional to knee 461 angular velocity. This occurs in the physical exoskeleton because we used the motor as a generator, where 462 rotational motion induces a voltage in the motor's windings and in turn a current that generates a magnetic 463 field that resists the motion of the knee [43,44]. At low velocities, current and therefore resistance, cannot be 464 generated. Our simulated controller replicated these effects. Therefore, in both simulation and experiment, 465 lower torques are applied at -10% in the penalize-low condition than at +10% in the penalize-high condition. 466 Our virtual participant study revealed that adding virtual participants tends to improve simulation 467 outcomes more so than adding muscle parameter. However, at larger numbers of virtual participants, 468 additional muscle parameter sets can meaningfully improve accuracy. If starting with one virtual participant 469 and one muscle parameter set, a 50% reduction in outcome measure RMS error (for example from 4% to 2%) 470 can be achieved by adding roughly 5-10 virtual participants. It is not possible to consistently achieve this 471 same reduction by adding additional muscle parameter sets alone. For most predictive simulation studies, 472 using one muscle parameter set per virtual participant appears to create sufficient variation in the dataset, 473 without adding unnecessary computational cost. However, accuracy can still be improved with a larger 474 muscle parameter set size, especially when the number of virtual participants is large. For example, if 475 50 virtual participants are included, a 50% reduction in outcome measure RMS error (for example from 476 0.3% to 0.15%) to can be achieved by adding roughly 1-5 muscle parameter sets. Therefore, in studies 477 where expected differences in outcome measures are on the order of 1% or less, for example when detecting 478 metabolic rate differences between running shoe designs [17], using a larger number of muscle parameter sets 479 could meaningfully improve accuracy.

480
Our predictive simulation approach has a number of fundamental and applied uses. Fundamentally, it 481 can be used to further investigate aspects of human optimization that are inherently difficult to test exper-482 imentally. For example, we can systematically alter the weighting of objectives to predict and potentially 483 understand the effect on gait adaptation. We can also implement physiologically realistic optimization al-484 gorithms, such as reinforcement learning, to understand how humans may adapt and learn gaits over time. 485 In a more applied sense, it can be used to explore strategies for improved exoskeleton design and control. 486 While here we used experimental results to validate our predictive simulations, in future we can do the 487 opposite, designing and iterating controllers in advance of human testing to help identify reasonable solution 488 spaces. Our ability to simulate diverse participant pools may also allow us to tailor controller design based 489 on individual users' abilities or disabilities.