## Abstract

Treadmill training is a common intervention to promote healthy walking function for individuals with pathological gait. However, because of the heterogeneity of many patient populations, determining how an individual will respond to new treadmill protocols may require extensive trial and error, causing increased patient fatigue. The purpose of this study was to develop and validate a framework for predictive simulation of treadmill gait, which may be used in the design of treadmill training protocols. This was accomplished through three steps: predict motion of a simple model of a block relative to a treadmill, create a predictive framework to estimate gait with a two-dimensional (2D) lower limb musculoskeletal model on a treadmill, and validate the framework by comparing predicted kinematics, kinetics, and spatiotemporal parameters across three belts speeds and between speed-matched overground and treadmill predictive simulations. Predicted states and ground reaction forces for the block-treadmill model were consistent with rigid body dynamics, and lessons learned regarding ground contact model and treadmill motion definition were applied to the gait model. Treadmill simulations at 0.7, 1.2, and 1.8 m/s belt speeds resulted in predicted sagittal plane joint angles, ground reaction forces, step length, and step time that closely matched experimental data at similar speeds. Predicted speed-matched overground and treadmill simulations resulted in small root-mean-square error (RMSE) values within standard deviations for healthy gait. These results suggest that this predictive simulation framework is valid and can be used to estimate gait adaptations to various treadmill training protocols.

## Introduction

Walking speed is a complex function of body kinematics, propulsive forces, muscle activation, and spatiotemporal parameters [1]. Heterogeneous patient populations often use unique gait compensatory strategies to achieve their self-selected or fast walking speed [2–4]. Treadmill-based gait training at fast speeds can improve walking function and efficiency, increase lower-limb kinematic symmetry, and reduce postural instability for individuals with neuromuscular diseases such as Parkinson's or stroke [5–7]. Although increasing walking speed is a common therapeutic goal, long and fatiguing gait training sessions are often required to study individual-specific speed-related changes in gait. Therefore, the ability to predict speed-related changes in gait mechanics would be beneficial, especially when individual-specific impairments are apparent.

Previous studies have used musculoskeletal models and simulations built from experimental data to provide insight into gait strategies, e.g., how joint contact forces and muscle function change at different walking speeds [8–11]. However, because these simulations rely on existing data at a set walking speed, they cannot predict how an individual may alter their gait at a new walking speed. The ability to predict how gait patterns change with treadmill belt speed may provide additional insight into how individual-specific impairments impact gait mechanics without requiring lengthy, fatiguing experimental studies.

Direct collocation methods that predict movement via trajectory optimization of high-level task goals have recently been employed to predict human gait [12] and offer new options to optimize rehabilitation protocols [13,14]. To our knowledge, this technique has not been used to predict gait adaptations to walking on a treadmill at a range of speeds. The objective of this study was to develop and validate a predictive simulation framework for treadmill gait. Because predictive simulations of treadmill gait have not been generated previously, it was unknown how a moving platform and a body moving relative to the platform would function in this simulation pipeline. Therefore, as a first step toward developing the treadmill gait framework, we began with a simple model of a block on a treadmill. The goal of the block-treadmill model was to predict and substantiate the contact, rotation, and translation of a block on a treadmill to verify that a predictive framework could be developed to estimate motion between a platform and a body moving relative to each other. As a second step in the predictive treadmill gait framework development, we used knowledge gleaned from the block model to build a 2D lower limb treadmill musculoskeletal model and to create a predictive framework to estimate physiologically realistic treadmill gait. Finally, to validate the framework, we compared predicted kinematics and kinetics across three treadmill belt speeds and between results from a predictive simulation of overground versus treadmill walking at the same speed.

## Methods

### Components of a Predictive Simulation Framework for Block Motion and Treadmill Gait.

Predictive simulations of human movement are typically framed as an optimal control problem and assume the central nervous system optimizes certain performance criteria to generate motion [15,16]. Although there are numerous approaches to solving optimal control problems (e.g., shooting methods, direct collocation, etc.), we chose to use direct collocation because it is more computationally efficient than traditional shooting methods [14,17,18]. Direct collocation optimal control methods approximate the system's time-varying states and controls as polynomial splines across a temporal grid and solve the discretized equations while optimizing the cost function [19,20].

Because direct collocation is a gradient-based optimization method, all model elements used were continuous and differentiable. We used a twice continuously differentiable Hunt–Crossley contact model to predict the interaction between the block or feet and the treadmill belt [12,21]. Additionally, for the gait model, we used a previously published muscle model that is compatible with direct collocation [14,18].

Finally, to complete the predictive framework we chose our cost function performance criteria [22]. Various performance criteria have been used to predict gait including minimizing cost of transport [23], minimizing fatigue [17], and multi-objective cost functions that minimize metabolic energy squared, muscle activations squared, and joint accelerations squared [12]. We chose to minimize effort for both the block controls and gait model muscles, as this is a very simple and commonly used cost function to predict overground gait [17].

### Development and Validation of a Block Model.

To demonstrate the feasibility of a predictive framework of motion of an object relative to a moving treadmill belt we constructed a simple model in opensim [24,25] consisting of a block on a treadmill connected via a planar joint. We assumed the treadmill to be massless, set the block's mass to 62 kg, and restricted the treadmill to only translate along the anterior-posterior axis. We prescribed the treadmill motion as a linear function with an intercept of 0 and a slope of –1.2. The intercept represented the starting position of the treadmill center in the anterior–posterior direction relative to the global origin in meters and the slope represented the treadmill velocity in meters per second. The negative sign defined motion in the posterior direction. Three coordinate actuators were applied to facilitate movement of the block in the sagittal plane, one for each degree-of-freedom: translation in the anterior-posterior and superior– inferior directions and rotation about the medial-lateral axis.

*J*, minimized the controls cubed, where

*u*(

*t*) is the time-varying force applied by each coordinate actuator and

*t*is the movement duration (Eq. (1)). We applied this cost function to generate all three motions. Kinematic constraints (Table 1) and initial and final conditions for the states and velocities (Table 2) were applied to define each motion separately. To validate the predicted block motions, root-mean-square error (RMSE) values were computed between the predicted and tracked states.

_{f}Kinematic constraints on block model | |||
---|---|---|---|

tx | ty | rz | |

Fall | $[\u22120.5\u2009m,\u20090\u2009m]$ | $[0\u2009m,\u20090.71\u2009m]$ | $[0\u2009deg,\u20090\u2009deg]$ |

Rotate | $\u2009[\u22120.5\u2009m,\u20090\u2009m]$ | $\u2009[0\u2009m,\u20090.8\u2009m]$ | $[0\u2009deg,\u200910\u2009deg]$ |

Translate | $\u2009[\u22120.5\u2009m,\u20090\u2009m]$ | $\u2009[0.08\u2009m\u2009,\u20090.09\u2009m]$ | $\u2009[0\u2009deg,\u20090\u2009deg]$ |

Kinematic constraints on block model | |||
---|---|---|---|

tx | ty | rz | |

Fall | $[\u22120.5\u2009m,\u20090\u2009m]$ | $[0\u2009m,\u20090.71\u2009m]$ | $[0\u2009deg,\u20090\u2009deg]$ |

Rotate | $\u2009[\u22120.5\u2009m,\u20090\u2009m]$ | $\u2009[0\u2009m,\u20090.8\u2009m]$ | $[0\u2009deg,\u200910\u2009deg]$ |

Translate | $\u2009[\u22120.5\u2009m,\u20090\u2009m]$ | $\u2009[0.08\u2009m\u2009,\u20090.09\u2009m]$ | $\u2009[0\u2009deg,\u20090\u2009deg]$ |

Degrees of freedom include translation in the anterior-posterior direction (*tx*), translation in the superior-inferior direction (*ty*), and rotation about the medial-lateral axis (*rz*). Distances measured from the global coordinate frame origin to the block center.

Initial and final conditions for states and velocities of block model | |||
---|---|---|---|

tx | ty | rz | |

Fall | (I)$\u20090\u2009m$ | (I)$\u20090.71\u2009m$ | (I)$\u20090\u2009deg$ |

(F)$\u20090\u2009m$ | (F) $[0.58,\u20090.6\u2009m]$ | (F)$\u20090\u2009deg$ | |

Rotate | (I)$\u20090\u2009m$ | N/A | (I)$\u200910\u2009deg$ |

(F)$\u20090\u2009deg$ | |||

(F)$\u20090\u2009m$ | |||

Translate | (I)$\u20090\u2009m$ | N/A | (I)$\u20090\u2009deg$ |

(F)$\u20090\u2009m$ | (F)$\u20090\u2009deg$ | ||

All | (I)$\u20090\u2009m/s$ | (I)$\u20090\u2009m/s$ | (I)$\u20090\u2009m/s$ |

(F)$\u20090\u2009m/s$ | (F)$\u20090\u2009m/s$ | (F)$\u20090\u2009m/s$ |

Initial and final conditions for states and velocities of block model | |||
---|---|---|---|

tx | ty | rz | |

Fall | (I)$\u20090\u2009m$ | (I)$\u20090.71\u2009m$ | (I)$\u20090\u2009deg$ |

(F)$\u20090\u2009m$ | (F) $[0.58,\u20090.6\u2009m]$ | (F)$\u20090\u2009deg$ | |

Rotate | (I)$\u20090\u2009m$ | N/A | (I)$\u200910\u2009deg$ |

(F)$\u20090\u2009deg$ | |||

(F)$\u20090\u2009m$ | |||

Translate | (I)$\u20090\u2009m$ | N/A | (I)$\u20090\u2009deg$ |

(F)$\u20090\u2009m$ | (F)$\u20090\u2009deg$ | ||

All | (I)$\u20090\u2009m/s$ | (I)$\u20090\u2009m/s$ | (I)$\u20090\u2009m/s$ |

(F)$\u20090\u2009m/s$ | (F)$\u20090\u2009m/s$ | (F)$\u20090\u2009m/s$ |

Distances measured from the global coordinate frame origin to the block center.

### Development of a Treadmill Gait Model and Predictive Cost Function.

We modified a 2D lower limb musculoskeletal model with a mass of 62 kg, ten degrees-of-freedom, 18 musculotendon actuators, and a coordinate actuator controlling lumbar flexion–extension to include a treadmill body with a prescribed linear motion [14]. Additionally, we placed two contact spheres on each foot—one on the heel and one on the forefoot—and the contact half space on the top surface of the treadmill (Fig. 2(a)). The slope of the prescribed treadmill motion was defined as one of three values: –0.7, –1.2, or –1.8, to represent the velocity (m/s) of the treadmill belt in the posterior direction at slow, comfortable, and fast speeds.

*a*(

_{i}*t*) of the

*i*th muscle in the system with m-muscles (m = 18) and the error between the reference and simulated states ($etrack$) from the initial time (

*t*

_{0}) to stride time ((

*t*, Eq. (2)). Additionally, the tracking problem was subject to physiologically realistic boundary constraints for each joint, and bilateral periodicity was imposed for coordinate values (except pelvis anterior– posterior motion), coordinate actuator controls, and muscle activations to output a full gait cycle while only tracking a half gait cycle ((Eq. (3)). Bilateral periodicity enforces periodic constraints on states (

_{f}*x*(

*t*)) and muscle activations so the right and left values, and vice versa, are identical at the beginning and half of the gait cycle [28]. Imposing bilateral symmetry via periodicity is a commonly used method to decrease computation time [15].

*v*

_{COM}, equal 0 m/s.

### Validation of Treadmill Predictive Framework.

It is well established for young adults that walking speed affects joint kinematics, ground reaction forces, and spatiotemporal parameters [29]. To validate our predictive framework of treadmill gait we computed the differences between predicted peak joint angles (hip, knee, and ankle flexion), peak ground reaction forces (vertical and horizontal), step length, and step time between the different belt speeds. We then compared those differences against changes reported in the literature for healthy adults walking at slow, comfortable, and fast speeds.

As a secondary means of validating the usefulness of the framework, we compared the predicted joint angles, ground reaction forces, and spatiotemporal parameters from the comfortable speed treadmill simulation against results from a speed-matched predictive simulation of overground gait. We used the same cost function, constraints, and boundary conditions for both the overground and treadmill predictive simulations. There were two main differences between the overground and treadmill simulations. First, in the overground simulation the contact model was defined between the feet and the ground whereas in the treadmill model, contact was between the foot and the top surface of the treadmill belt. Second, for the overground tracking simulation the unmodified predicted states, with anterior translation of the pelvis, from a previous predictive simulation of healthy overground gait [12] were set as the reference coordinates. Previous experimental gait analyses showed that sagittal plane joint angles, ground reaction forces, and spatiotemporal parameters are not significantly different when walking on a treadmill and overground at the same speed [30,31]. Root mean square errors were computed between the predicted treadmill and overground joint angles and ground reaction forces for validation. The block and treadmill gait models and corresponding MATLAB simulation code are available to access on simtk.^{2}

## Results

### Block Model.

Tracked and predicted states matched closely for all block motions, with RMSE values of 0.2 m, 4 deg, and 0.04 m for falling, rotating, and translating, respectively (Fig. 3, top). Peak vertical ground reaction force (vGRF) for the falling block occurred at the instant of initial contact of the block with the treadmill. Once an equilibrium position was reached the vGRF equaled one block weight (bw). During rotation the initial vGRF was slightly below one bw since most, but not all the weight was on top of the contact sphere. As the block rotated less of the weight was over the sphere and the vGRF decreased. Finally, in translation a constant vertical position of the block center corresponded to a constant predicted vGRF of one bw, suggesting the solution converged on the equilibrium position of the contact model (Fig. 3, bottom).

### Validation of the Treadmill Predictive Framework.

Predicted and tracked sagittal plane hip, knee, and ankle angles, and vertical and horizontal ground reaction force trajectories align well for the treadmill gait model simulated in the comfortable speed condition. RMSE values between tracked and predicted joint angles were 5.50 deg for hip flexion, 6.46 deg for knee flexion, and 8.49 deg for ankle plantar and dorsiflexion. Ground reaction force RMSE values were 0.08 body weight (BW) for vertical and 0.3 BW for horizontal ground reaction forces.

The treadmill framework successfully predicted many trends commonly observed experimentally with changes in walking speed including changes in sagittal plane joint angles, ground reaction forces, and spatiotemporal parameters (Fig. 4). Peak hip extension increased from 5.11 deg to 14.90 deg and ankle plantarflexion increased from 0.62 deg to 15.01 deg from the slow to the fast condition. Peak vertical and anterior ground reaction forces also increased with the increasing belt speed. From the slow to fast belt speed, peak vGRF ranged from 1.08 BW to 1.95 BW, and peak anterior ground reaction force ranged from 0.1 BW to 0.23 BW. Finally, step length increased from 0.39 m to 0.56 m while step time decreased from 0.59 s to 0.40 s from the slow to the fast condition.

Root-mean-square error values comparing speed-matched (1.2 m/s) overground and treadmill gait predictions were 2.22 deg for hip angle, 2.94 deg for knee angle, and 3.16 deg for ankle angle (Fig. 5, top). Vertical and horizontal ground reaction forces were also comparable between the comfortable belt speed and overground conditions with RMSE values of 0.37 BW and 0.09 BW (Fig. 5, bottom). Spatiotemporal parameters between the overground and comfortable belt speed predicted simulations with differences of 0.02 m in step length and 0.03 s in step time.

## Discussion

We successfully developed the first predictive simulation framework for treadmill gait and used it to solve two predictive problems: motion of a block and motion of a 2D gait model relative to a moving treadmill belt. Because this was the first simulation framework of gait relative to a moving ground, it was important to begin the development with the simple model of a block to mimic foot-treadmill contact and then apply an understanding of contact model definition and kinematic constraints to the more complex gait model.

Predicted block states for all motions closely matched tracked states. Although predicted peak vGRF during falling is high at 8.1 bw compared to average peak vGRF recorded during human gait of 1 to 1.5 BW [32], this result may still be reasonable. Block impact force should be higher given that the block is free-falling while the foot is in a controlled descent into heel-strike. Additionally, the block weight is entirely concentrated over the contact sphere while body weight is distributed over both feet at heel strike. Values of predicted vGRF are also consistent with basic rigid body dynamics for the rotation and translation simulations.

Once we applied the knowledge from the block model to construct the 2D treadmill gait model and predictive simulation framework, we were able to capture different mechanisms to achieve each gait speed. As an initial validation for the framework we first compared the predicted results from the comfortable speed treadmill simulation to experimental treadmill data for healthy young adults. Predicted kinetics and kinematics including peak anterior and vertical ground reaction forces and ankle plantarflexion and hip extension angles were all within two standard deviations of the experimental mean [31], meeting the validation criteria established by Hicks et al. [33]. Predicted step length [34] and step time [35] were also within two standard deviations of experimental gait data at comfortable speeds.

Predicted speed-related changes in gait kinematics, kinetics, and spatiotemporal parameters were also consistent with experimental results. The increase in peak ankle plantarflexion by about 15 deg and hip extension by about 10 deg from the slow to fast belt speed is consistent with experimental data of young adults walking on a treadmill with a similar speed range between slow and fast conditions [36]. The predicted increase in peak vGRF by about 0.9 BW and peak anterior force by about two times from the 0.7 to 1.8 m/s belt speed conditions, matches experimental results from a study with a similar walking speed range of 1 to 2 m/s [37]. Finally, the predictive simulation framework estimated an increase in step length by 43% and a decrease in step time by 32% with an increase in speed consistent with previous experimental studies (+56% step length, –32% step time) [35,38].

In addition to the speed comparisons, the predictive treadmill framework also estimated sagittal plane joint angles and ground reaction forces similar to a model of overground gait with matched speed. RMSE values for the hip, knee, and ankle were all within typical standard deviation ranges for healthy individuals during gait [36]. Errors between the predicted vGRF and horizontal ground reaction force curves for overground and treadmill gait were also similar to differences observed in a previous experimental study comparing treadmill and overground gait at matched speeds [31]. The small differences between the simulated predicted joint angles and ground reaction forces in the speed-matched overground and treadmill simulations are likely due to the differences in the contact model definition and the states in the initial guess. However, overall the close match between the predicted simulation results and experimental treadmill gait data at a range of belt speeds and between the predicted overground and treadmill simulations at matched speeds helps to validate this novel framework for predictive simulation of treadmill gait. Future work will build off these results and use this predictive framework with a musculoskeletal model representing individual-specific impairments (e.g., range of plantarflexor weakness) to estimate compensatory gait strategies employed at various belt speeds.

As with all modeling studies, there are limitations to this predictive treadmill gait framework. The application of a 2D model limited reproduction of out of plane movements and this, along with the small number of muscles in the model for simplicity may explain why the predicted results did not capture increased knee flexion during swing phase typically observed with increased belt speed [39]. Predicted vertical and horizontal ground reaction force curves also slightly deviated from the characteristic experimental trajectories [37]. To improve this result, more contact spheres could be added to each foot and the properties of the spheres could be calibrated to be individual-specific [13]. Adjusting cost function terms and applying a different set of reference coordinates to the tracking simulation may also result in more accurate solutions. However, despite these limitations, the framework captured many kinematic, kinetic, and spatiotemporal patterns commonly observed in experimental treadmill gait and may be a valid tool to predict treadmill gait mechanics.

This workflow established the first predictive framework to estimate several common gait adaptations on a treadmill at a range of belt speeds. Future work will explore the feasibility of this framework to predict gait adaptations on a treadmill when modeling impairments like unilateral muscle weakness following a stroke and to predict responses to training protocols that combine treadmill training with other interventions like functional electrical stimulation or robotic exoskeletons.

## Acknowledgment

The authors would like to thank Nick Bianco, Christopher Dembia, Carmichael Ong, and Jon Stingel for their assistance with this work.

## Funding Data

Delaware Space Grant Consortium (Grant No. 80NSSC20M0045; Funder ID: 10.13039/100005680).

National Science Foundation (Grant No. NSFGRFP 1940700; Funder ID: 10.13039/100000001).

University of Delaware (Grant No. Department of Mechanical Engineering Helwig Fellowship; Funder ID: 10.13039/100006094).