Evaluating cost function criteria in predicting healthy gait

https://doi.org/10.1016/j.jbiomech.2021.110530 0021-9290/ 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). ⇑ Corresponding author. E-mail address: k.veerkamp@amsterdamumc.nl (K. Veerkamp). K. Veerkamp a,b,c,⇑, N.F.J. Waterval , T. Geijtenbeek , C.P. Carty , D.G. Lloyd , J. Harlaar , M.M. van der Krogt a


Introduction
Forward dynamic simulations can be used to predict new patterns of human movement without requiring experimental data. Thus, such simulations have the potential to investigate the mechanistic cause-and-effect relationships in movement impairments, optimise the prescription of assistive devices, and predict outcomes of interventions. A critical aspect in predictive simulations is the cost function, which sets the goal of the movement and aims to replicate human gait optimisation. At present, it is not clearly known how human gait is optimised, and how to specify and weigh the involved criteria.
A number of different physiologically-based criteria have been used in predictive simulations of human gait. Many experimental studies have shown that the energetic cost of transport (CoT) is involved in human gait selection (Abram et al., 2019;Bertram and Ruina, 2001;Donelan et al., 2001;Minetti et al., 2020;Selinger et al., 2015;Zarrugh et al., 1974) and, consequently, it is also the most commonly considered criterion in predictive simulations of gait (Anderson and Pandy, 2001;Dorn et al., 2015;Falisse et al., 2019;Lai et al., 2018;Ong et al., 2019;Song and Geyer, 2015;Wang et al., 2012). Besides energy cost, human gait might also be optimised to prevent muscle fatigue by minimising muscle https://doi.org/10.1016/j.jbiomech.2021.110530 0021-9290/Ó 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). activation, as suggested by the muscle force-endurance relationship (Crowninshield and Brand, 1981;Müller and Grosse-Lordemann, 1937), and as applied in several predictive simulation studies (Ackermann and van den Bogert, 2010;Falisse et al., 2019;Lai et al., 2018). Head stability measures have also been used in simulation studies (Dorn et al., 2015;Nguyen et al., 2019;Ong et al., 2019) supported by the finding that human head accelerations, detected by the otolith organs in the ear (Kandel et al., 2013), remain similar across different walking conditions (Menz et al., 2003). Further, a recent commentary has suggested how the derivative of force, defined as 'yank', is an important factor in sensorimotor systems (Lin et al., 2019). Indeed, mitigating the impact of external forces, or yank, which could be detected by the cutaneous mechanoreceptors in the foot (Kavounoudias and Roll, 1998;Kennedy and Inglis, 2002), may be involved in human walking, as a high rate of loading has been suggested to contribute to injuries by applying high stresses to the leg's tissues (Mündermann et al., 2005;Nigg, 1985). Finally, the musculoskeletal models in previous simulation studies included a joint limit torque, representing ligaments, that are activated during joint hyperflexion or -extension (Dorn et al., 2015;Geyer and Herr, 2010;Ong et al., 2019). Mimicking type III joint receptors firing when reaching these extreme joint angles (Nyland et al., 1994;Zimny and Wink, 1991), those studies minimised using these joint limit tor-ques (Dorn et al., 2015;Geyer and Herr, 2010;Ong et al., 2019), which is particularly relevant for the knee joint that may approach end-of-range of motion during gait.
Although it is likely that human gait is optimised for a variety of these physiologically-based criteria, it is unclear how each of these criteria contribute to gait. To achieve more representative predictions, it is necessary to better understand the individual effects of these criteria, as well as their relative contribution to human gait. Therefore, this study's two aims were to 1) evaluate the effects of minimising the criteria CoT, muscle activity (MusAct), head stability (HeadStab), foot-ground impact (FGImpact), and knee hyperextension (KneeExt) on gait using predictive simulations, and 2) develop and evaluate an optimally weighted, combined cost function tuned to predict healthy gait.

Experimental gait data
Motion capture data was collected using the Plug-In-Gait marker set (Vicon Nexus 2.6.1; Davis, Õunpuu, Tyburski, & Gage, 1991;Kadaba, Ramakrishnan, & Wootten, 1990) for ten healthy adults (4 males, age 26.8 ± 2.6 years, mass 67.2 ± 8.5 kg, height 1.76 ± 0.0 8 m) walking overground at their comfortable walking speed. * indicates thresholds differentiating gait phases that were optimised. C0 was a constant stimulation value. PD was a proportional-derivative reflex, based on the pelvis tilt angle. F was a force-based reflex, L a length-based reflex. + and -were positive and negative feedback, respectively. Reflex loops mostly worked within muscles, but also between muscles (i.e., iliopsoas-hamstrings and soleus-tibialis anterior), to prevent co-contraction.
Ground reaction forces (GRFs) were simultaneously collected from two embedded force plates (AMTI, Watertown, MA, USA). Musculoskeletal modelling was conducted in OpenSim (version 4.0; Seth et al., 2018), scaling the gait2392 musculoskeletal model (Delp et al., 1990) to each participant's proportions using a static T-pose trial. Inverse kinematic and inverse dynamic tools were used to compute joint angles and moments from walking trials. Joint powers were calculated from these angles and moments using custom-written scripts in MATLAB 2016a (The MathWorks Inc., MA). The ensembled averaged GRFs and joint angles, moments and powers were used in further comparisons as normative experimental data for healthy gait. Additionally, normative muscle activation patterns were derived from Bovi et al. (2011), including electromyography data for the gluteus maximus, rectus femoris, biceps femoris, vastus medialis, gastrocnemius medialis, soleus and tibialis anterior muscles using a ZeroWire system (Aurion, Milano, Italy).

Base simulation framework
Our base simulation framework was developed using the opensource software SCONE (version 1.4.0; Geijtenbeek, 2019). A generic musculoskeletal model (gait2392; Delp et al., 1990) was simplified to nine degrees of freedom (sagittal plane only; pelvis and trunk merged) and nine musculotendon units per leg. For each group of muscles working around the same degrees of freedom, the musculotendon pathway and parameters from a representative muscle were selected, and its maximum isometric force was updated by adding the forces from the other muscles to that group. The Millard-equilibrium muscle model was used  with Table 1 providing an overview of its parameters. Of note, tendon strain at maximum isometric force was 4.9% for all muscles Rajagopal et al., 2016), except for the gastrocnemius and soleus that were set to 10% (Arnold et al., 2013;Ong et al., 2019). The ratio between slow and fast twitch fibres for each muscle was similar to Ong et al. (2019). A knee end-of-range of motion limit torque was added that was activated when the knee extended below 5 degrees of flexion (Markolf et al., 1990), with a stiffness of 2 N*m/degree and a dampening of 0.2 N *m/(degrees*second). A Hunt-Crossley foot-ground contact model with two viscoelastic spheres at each foot was added (Hunt and Crossley, 1975) and the geometry of the spheres was initially optimised minimising the error between simulated and experimental GRFs and kinematics (see Appendix A for details).
A reflex-based controller (Geyer and Herr, 2010) was implemented into our base simulation framework, detailed in Table 2. This type of controller results in intrinsically stable neuromuscular simulations, and was therefore favoured over open loop controller Fig. 1. The gait patterns predicted by minimising for different criteria, compared to normative experimental data. The variables from the simulated gait were averaged over time-normalised gait cycles in the ten-second simulation, excluding the first two gait cycles. in this study. The controller consisted of different types of muscle force-and length-based reflex loops working within and between muscles, activating the musculotendon units in different gait phases (i.e., early stance, late stance, pre-swing, swing and late swing). The trunk was controlled by a proportional-derivative (PD) reflex controller based on the pelvis + trunk tilt angle. Reflex delays were also set as in Geyer and Herr (2010). For the initial state, the initial pose rate of change was selected from pilot experiments. The initial muscle activation was set according to initial controller output, while the initial length of the contractile element was found after equilibration of the muscle tendon dynamics, based on initial muscle-tendon length, velocity and activation. The model's initial pose, the thresholds between gait phases, and the reflex gains and offsets of each muscle for each gait phase were optimised. There were 98 design variables: seven for the model's initial pose, four to differentiate gait phases, and 87 for the muscle reflexes. Except for the initial pose, all design variables were symmetric for the left and right legs.
The design variables were optimised by adding different cost functions, specified below, to the base simulation framework. When optimising for a cost function, the score, which is the weighted sum of the involved cost function criteria over the tensecond simulation, was minimised using a Covariance Matrix Adaptation Evolution Strategy (CMA-ES; Hansen 2016). A 24-core Intel Xeon CPU E5-2690 2.60 GHz processor was used. In the optimisation, penalties were applied when, (1) the model walked below a speed of 0.5 m/s, (2) the ankle angle was beyond 60 degrees dorsi-or plantarflexion, (3) simulation time was less than ten seconds, and (4) the centre of mass was below 0.75 m. These penalties were used to avoid unrealistic local minima, with the latter measure, indicating a fall, also prematurely ending simulations to increase optimisation speed. For each cost function, six optimisations with different random seeds were performed, until, averaged over the last 500 generations, the cost had not improved more than 0.01% per generation. The simulation with the lowest score over all six optimisations was selected for subsequent evaluation. From the ten-second simulation, the first two gait cycles were excluded, and waveforms were time-normalised and averaged over all subsequent gait cycles, which were typically all very similar.

Evaluation cost functions
Various cost functions were added to the base simulation framework, consisting of the following criteria in different compositions: i. CoT, an overall effort measure, was calculated using the muscle metabolic model by Umberger, Gerritsen, & Martin (2003) and implemented based on Uchida, Hicks, Dembia, & Delp (2016). The calculated energy rate was summed for all muscles (m) over the simulation (until t end ), and divided by the mass of the model and distance travelled, i.e.: Muscle fatigue (MusAct) was quantified by the muscle activation squared. The activation of each muscle (m) was squared, summed over the simulation, and divided by the distance travelled, i.e.: iii. Head stability (HeadStab) was represented by the acceleration of the head segment. The sum of the horizontal and vertical accelerations of the head centre of mass (a x and a y , respectively) during the simulation was divided by the distance travelled, i.e.: Fig. 2. Quantified agreement (R 2 and RMSE) of the gait patterns predicted by minimising for different criteria with normative experimental data. RMSE was normalized to the standard deviation of each variable. The right column shows the walking speed for each criterion. HeadStab iv. The foot-ground impact (FGImpact) was assessed by the GRFs derivative or yank, penalizing high increases/decreases in the GRF. The sum of the time-derivative of horizontal and vertical GRFs (GRF x and GRF y ) from both the left and right leg during the simulation was divided by the distance travelled, i.e.: KneeExt involved a penalty to the use of the knee limit torque, representing ligaments, (F limit ) that was summed over the whole simulation for both the left and right knee, i.e.: KneeExt became zero when the knee limit torque was not used (i.e., more than 5 degrees of flexion during the whole simulation) and for this reason, it could not be optimised independently.
Note that t end was ten seconds for each simulation. Besides these five physiologically-based criteria, a cost function criterion was used to examine the base simulation framework's best possi- Fig. 3. The gait patterns predicted by each step in the combined cost function, compared to normative experimental data. The variables from the simulated gait were averaged over time-normalised gait cycles in the ten-second simulation, excluding the first two gait cycles. ble agreement with experimental data (ExpTrack). In this, we minimised the root mean square error (RMSE) between predicted and experimental GRFs, joint angles, and muscle excitations. To give equal weighting to each biomechanical category, the weighting of each variable was divided by its experimental standard deviation, and then divided by the number of variables within each category. Each cost functions' performance was assessed by the agreement between the average predicted and normative experimental waveforms for different categories of biomechanical variables, i.e. GRF, joint angles, moments, powers, and muscle excitations. The performance metrics were coefficients of determination (R 2 ) and RMSE. RMSE of the excitations was excluded since matching magnitudes is not straightforward as these rely on many factors, such as EMG normalisation (Devaprakash et al., 2016). R 2 was considered very weak when smaller than 0.3, weak between 0.3 and 0.5, moderate between 0.5 and 0.7 and strong above 0.7 (Hair Jr et al., 2016). The RMSE for each variable was divided by its average experimental standard deviation, to allow more fair comparison of the RMSE between variables with different magnitudes. R 2 and RMSE were evaluated on three levels: for each variable separately, averaged over each biomechanical category and averaged over all categories. Walking speed was also assessed.
The cost functions were analysed in three stages. In the first stage, criteria i-iv (CoT, MusAct, HeadStab and FGImpact) were optimised independently and their quantified performance was evaluated. Second, these criteria were ranked by their overall average R 2 when optimised independently and then added in a stepwise way to create a combined cost function. The criterion with the highest average R 2 was selected with a set weighting of 1, and the criterion with the second-highest average R 2 was added using five different systematically chosen weightings. Out of these, the weighting providing the best average R 2 and RMSE was selected, after which the next criterion was added with five weightings, and so on until all criteria were combined into an optimal, combined cost function. KneeExt's weighting was tuned when the knee showed more extension than seen in experimental data. For each criterion, its tested weightings were normalised to the criterion's value when independently optimised, except for KneeExt since its minimal value was zero, to be able to evaluate the contribution of each criterion in the combined cost function. When a new criterion was added, the optimised design variables from the previous stepwise combination of cost functions were used as initial values. In stage three, ExpTrack was independently optimised for and the performance of the optimised, combined cost function was compared to this simulation. Here, R 2 and RMSE for the moments and powers were excluded in the evaluation, since they could not be tracked using ExpTrack. Fig. 4. The muscle excitation patterns predicted by each step in the combined cost function, compared to normative experimental data (Bovi et al., 2011). The excitation patterns from the simulated gait were averaged over time-normalised gait cycles in the ten-second simulation, excluding the first two gait cycles.

Results
When optimising for each physiologically-based criterion independently, different criteria performed best for distinct biomechanical categories (Figs. 1 and 2; more detailed in Appendix B). Minimising CoT was best for predicting kinematics as it scored highest R 2 and lowest RMSE for all angles (strong average R 2 : 0.80; RMSE: 1.54 SD), although, it was worst for predicting GRFs compared to the other criteria (moderate average R 2 : 0.53; RMSE: 9.09 SD). It resulted in the highest walking speed of 1.36 m/s. Minimising for MusAct resulted in the lowest average RMSE for joint powers (3.66 SD), and resulted in the slowest walking speed (0.51 m/s). Over all criteria, minimising HeadStab was best for pre-dicting GRFs (strong average R 2 : 0.92; RMSE: 3.22 SD) and moments (moderate average R 2 : 0.68; RMSE: 2.48 SD), performing best for the knee moment (moderate R 2 : 0.54: RMSE: 2.33 SD). It resulted in a walking speed of 0.92 m/s. FGImpact minimisation was best for R 2 of the joint powers (weak average R 2 : 0.46), for which it resulted in the best ankle push-off power (moderate R 2 : 0.67; RMSE: 3.30 SD). Its walking speed was 1.13 m/s. Overall, average R 2 was highest for FGImpact (0.56), followed by HeadStab (0.51), CoT (0.49) and MusAct (0.37). Hence, the cost functions were combined in this stepwise order.
When combining the criteria, the average R 2 and RMSE improved in each step, indicating that each criterion contributed positively to the combined cost function (Figs. 3, 4 and 5; more Fig. 5. Quantified agreement (R 2 and RMSE) with normative experimental data of the gait patterns predicted in each step of combining the cost functions. RMSE was normalized to the standard deviation of each variable. The right column shows the walking speed in each step. Table 3 Overview of composition of the cost function combining all criteria, and the minimised cost function score.

Criteria
Criterion's score when optimised independently Criterion's optimal cost function weighting* Criterion's weighted score in best simulation optimal combined cost function detailed in Appendix C). When combining all cost function criteria, the average RMSE was 2.10 SD and the average R 2 was 0.72, performing strong for GRFs, kinematics and moments (average R 2 : 0.93, 0.80 and 0.83, respectively), moderate for powers (average R 2 : 0.67) and weak for excitations (average R 2 : 0.36). In this cost function, CoT had the biggest contribution (Table 3) and it also had the highest normalised weighting. The performance of the combined cost function was close to the maximal quantified agreement obtained from the ExpTrack simulation, (Figs. 6 and 7; comparison of excitation patterns in Appendix D). When comparing GRFs, joint angles and excitation patterns, the average R 2 (0.70) was only slightly smaller than for the Exp-Track simulation (0.73), while the average RMSE (2.15 SD) was slightly higher than for ExpTrack (1.15 SD).

Discussion
In this study we demonstrated that optimising gait for FGImpact, HeadStab, CoT, and MusAct independently predicted different gait patterns which showed good agreement with experimental data for distinct biomechanical categories. Combining physiologically-based criteria with tuned weightings in a stepwise approach resulted in an overall improved agreement with experimental data. Compared to ExpTrack simulations, the combined cost function performed very well, which indicates that its agreement with experimental data is close to its maximum given the underlying constraints of the framework.
When optimising for the criteria separately, FGImpact overall showed the best R 2 compared to experimental data and its predicted walking speed was closest to 1.2 m/s, which in general cor- Fig. 6. The gait pattern and muscle excitation patterns predicted by the cost function combining all criteria, compared to a simulation tracking normative experimental GRFs, joint angles and muscle excitations. The variables from the simulated gait were averaged over time-normalised gait cycles in the ten-second simulation, excluding the first two gait cycles. K. Veerkamp, N.F.J. Waterval, T. Geijtenbeek et al. Journal of Biomechanics 123 (2021) 110530 responds to the preferred walking speed of humans (Ralston, 1958). Its main benefit was the improved ankle push-off, with a better timing of the ankle power compared to the other criteria, which showed a peak in ankle power generation either too early (CoT and HeadStab) or too late (MusAct). Remarkably, even though CoT had only a weak R 2 when minimised individually, it had the highest contribution in the optimal combined cost function (Table 3), underlying its importance in human gait as also shown before in predictive simulations and human experiments. Optimising for CoT by itself led to relatively poor GRFs with a high peak directly after initial contact for both the anterior-posterior and vertical components. Ackermann and van den Bogert (2010) showed comparable peaks in the GRFs when optimising for comparable effort and fatigue criteria. Contrarily, they showed that minimising for a MusAct criterion contributed to a physiologically consistent early stance knee flexion, whereas we did not see this when optimising for MusAct. This difference could possibly be explained by a difference in walking speed, which was set to 1.1 m/s in the study of Ackermann and van den Bogert, but was self-selected in our study and went down to 0.5 m/s for MusAct. Our combined, optimal cost function improved upon previous work using a reflex-based controller in SCONE by showing a better agreement with experimental data (current study vs Ong et al., 2019: R 2 : 0.72 vs 0.63; RMSE: 2.10 vs 2.47 SD). The main improvement was obtained for the knee moment, which lacked an external knee flexion moment in early stance in their results, even though their framework included minimisation of the head acceleration (HeadStab) and of using the knee limit torque (KneeExt), which we found to be important for obtaining knee flexion in early stance. Other recent work, using open-loop control with a more complex model and a manually-tuned five-term cost function, also found limited knee flexion in early stance after including metabolic energy rate in the cost function (Falisse et al., 2019). These comparisons underline the importance of not only including multiple cri-teria but also optimising their relative contributions to get an optimal agreement with experimental data.
Overall, performance of our framework was worst for tracking the experimental muscle excitations, even in the ExpTrack simulation. The reflex-based controller may be constraining the match with experimental EMG, not allowing a wide enough range of muscle excitation pattern for better tracking. It has been suggested that reflexes, and their spinal loops, are important during walking (Duysens and Forner-Cordero, 2019), but it is hard to validate which type of reflexes are predominant for each muscle during walking, as well as their gain values. We would suggest for future studies to further tune the controller, by broadly testing and implementing different sets of reflexes. For example, anticipatory activation patterns that act to stabilize the knee in the frontal plane (Besier et al., 2001;Lloyd and Buchanan, 2001) may be underpinned by ligamento-muscular reflexes (Kim et al., 1995). Additionally, the musculotendon parameters of the musculoskeletal model in our framework were based on generic values and may not be accurate for our population, possibly explaining the overall bad performance for muscle excitations while performing moderate and strong for other biomechanical categories. Previous studies have shown the importance of personalising these musculotendon parameters (Lloyd and Besier, 2003) and we would suggest for future studies to apply this personalisation, especially when working with pathological gait (Davico et al., 2020;Veerkamp et al., 2019). Additionally, the match with the experimental moments and powers was generally not as good as for the GRFs and kinematics. This could have been affected by the impact artifacts as seen in the experimental data, which are too rapid to be produced by the muscle model. Since this concerns only a small part of the gait cycle (~3%), effects are expected to be minimal. Likely, it is also caused by only having two foot-ground contact points, limiting the path of the centre of pressure, directly affecting the joint moments and thus also joint powers. Fig. 7. Quantified agreement of the gait predicted by the cost function combining all criteria and a simulation tracking normative experimental data with normative experimental data. RMSE was normalized to the standard deviation of each variable. K. Veerkamp, N.F.J. Waterval, T. Geijtenbeek et al. Journal of Biomechanics 123 (2021) 110530 Several factors need to be considered to ascertain whether our weighted, combined cost function is indeed optimal. A more ideal approach to obtain an optimal cost function matching experimental data is to perform a bilevel optimisation, in which the weightings are optimised such that an optimal agreement with experimental data is obtained (Mombaur and Truong, 2010;Nguyen et al., 2019). We were limited by computational power to do so, as solving for a cost function took around 14 hours using a shooting approach. The use of open loop control with trajectory optimization methods, such as direct collocation, is generally much faster, enabling the use of such bilevel optimisation. However, those simulations are not intrinsically stable, meaning they cannot handle uncertainty (such as perturbations and noise) like reflex-based controllers (Groote and Falisse, 2021), limiting their application. Nevertheless, both bilevel optimisation and the stepwise approach we used could be considered a fitting process, with the risk of overfitting, especially given the low dimensionality of human gait (Schwartz and Rozumalski, 2008). In our stepwise design, when adding a criterion to the cost function, the match with experimental data could not get worse than the previous result that lacked this criterion (assuming the global minimum was found each time), and it got close to optimal for our simulation framework. However, this does not mean that this cost function is necessarily best in predicting gait in all new conditions. Hence, to be able to conclude what the best cost function in terms of predictive ability is, cross-validation on other gait conditions is needed. For example, in a validation test, a person could walk at different speeds and on different slopes, or with assistive/resistive devices, and the best cost function should be able to capture all observations in the different experimental conditions. Also, before clinical application, it is important to establish whether the used cost function is valid for a range of pathological conditions. A shortcoming of our framework is that it is only twodimensional and will therefore only be applicable to answer questions in which gait deviations occur mostly in the sagittal plane. A three-dimensional analysis would become far more complex, and other factors, such as medio-lateral stability, as discussed above, will come into play. Also, more muscles and more control parameters will have to be added to the framework, which will considerably slow down optimisation times even further. Another challenge of our framework is that within predictive simulations the risk of ending up in local instead of global minima depends on chance. Therefore, we optimised each scenario for six random seeds enabling broad exploration. Even though we found differences between the outcomes of these seeds to be minimal (Appendix E), we cannot be sure that when the number of seeds would have been increased, possibly lower and slightly different minima would have been found. However, given the fact that we observed clear trends in the stepwise combinations, and both the RMSE and R 2 clearly improved when combining criteria, we are confident that it would not have affected the interpretation of the results. Lastly, in this study, the speed and step length differed between simulations with different cost functions, which allows the evaluation of how these variables change in different conditions in future work (i.e., in pathology). However, when interpreting the mechanical effects of the different cost functions, it needs to be realized that the speed and step length can be confounding factors.
In conclusion, this study showed how various optimisation criteria can contribute to different gait features and that careful weighting of them is essential in predicting healthy gait. Remarkably, a high contribution of CoT to gait only comes to expression when combined with other criteria. After developing our framework incorporating an optimal cost function that is tuned to perform well in predicting healthy gait, a next step is to validate its predictive ability. This step is important to be able to answer clin-ical what-if questions and to get closer to the ultimate goal of predicting treatment outcomes.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  composed of spheres of specified size and force parameters, which are positioned at specified locations on the foot. Previous studies used different number of spheres, location of the spheres and the force parameters (i.e. stiffness and dampening), which are often determined by manual tuning and might therefore be suboptimal and not consistent for the experimental data. Given recent research suggesting that the predicted gait is highly sensitive to changes in the foot contact model (Millard and Mombaur, 2019), we decided to tune the foot-ground contact model in such a way that it can fit with both experimental kinematics and GRFs as good as possible. In order to do this, we setup a tracking optimisation within SCONE, in which the root mean square error (RMSE) between the predicted and experimental kinematics and GRFs was minimised while the geometries of the spheres (i.e. x-position, y-position, and radius) were added to the design variables. The stiffness, dissi-pation, and friction of both spheres were set to pre-determined values, based on the values used in previous studies. We performed the optimisation with three random initial guesses and the results from the best simulation (i.e. lowest RMSE) were selected, which are shown in Fig. A1. Table A1 provides all parameters for the foot-ground contact model used in the predictive simulations in this study.
Ground reaction forces joint angles.

Appendix B
There was no EMG available for the iliopsoas and the biceps femoris short head (see Fig. B1).  . The variation between the simulation outcomes of the optimal combined cost function, in which six different random seeds were used. Blue indicates the best score while red is the worse score.