Development of a Subject-Specific Foot-Ground Contact Model for Walking OPEN ACCESS

Musculoskeletal models have been used to analyze the kinematics and muscle function of locomotor activities such as walking and running [1–4] for healthy individuals and those with osteoarthritis [5,6], poststroke hemiparesis [7,8], and cerebral palsy [9–11]. However, rarely have they been used to design an actual clinical treatment for a disorder affecting walking ability [12–16]. Ideally, clinically achievable changes in surgical parameter values, muscle strength, or neural control strategy could be input into a patient-specific computational walking model, and the resulting changes in the patient's joint loads, walking pattern, or walking speed could be accurately predicted. Such capability could facilitate the design of more effective treatments for walking impairments by allowing clinicians to customize treatment based on objective predictions of post-treatment walking function.

A primary reason for limited model use in treatment design is the high level of modeling complexity required to predict how a proposed treatment will change a patient's gait pattern. One of the most challenging aspects of predicting post-treatment gait patterns is modeling of foot-ground contact mechanics. The ability to predict an individual patient's GRFs and moments accurately is critical since model estimates of lower extremity joint motions and loads are dependent on how foot-ground interaction occurs [4]. While generic foot-ground contact models have been used to explore a variety of research questions, patient-specific foot-ground contact models [17] are needed to optimize clinical treatment design for specific patients, especially since a patient's foot-ground contact pattern is likely to change after treatment.

Previous studies have used foot-ground contact models to reproduce experimental GRFs for walking [2–4,18–29], running [3,4,25,30–32], and running jumps [33] (Table 1). Most of these studies modeled foot-ground interactions using a grid of viscoelastic elements placed on the bottom of the foot, with the most notable differences between models being the number and location of elements and the number of foot segments. Researchers have placed viscoelastic elements at a single point [31], along the midline [19,20] or CoP [24] of the foot, at a small number of specific points on the bottom of each foot segment [4,18,21, 22,27,28,30,33], and at 30 or more locations under the foot [2,3,26,32]. Most studies used either a single-segment [3,22, 25,28,30,31,33] or two-segment (HF and FF) [2,4,18–21, 24,27,29] foot model. GRFs produced by single-segment foot models exhibited more discontinuities, especially at heel strike, than did those produced by two-segment foot models, possibly due to their inability to account for rolling motion under a single rigid body. In contrast, new constraint-based single-segment foot models that account for rolling motion have produced realistic GRFs without discontinuities [25]. However, the rolling constraint cannot be used in predictive optimizations because of its dependence on measured CoP data. Some two- and three-segment foot models have also produced discontinuous and/or oscillatory GRFs [18,24,26,27,29]. Results improved when the contact model was modified (e.g., viscoelastic versus volumetric) or when the number of contact elements was increased. Nonetheless, many reported models do not accurately reproduce experimental GRF curves, which is problematic for predictive gait optimization studies that require accurate prediction of GRFs to predict physically realistic walking patterns and corresponding joint moment or muscle force controls.

Despite a large number of foot-ground contact models published in the literature (review Table 1), few studies have reported model predictions of ground reaction moments or CoP location and free moment for walking. Omission of these quantities makes it difficult to predict three-dimensional walking patterns accurately using musculoskeletal models. Two-segment foot models used within a three-dimensional musculoskeletal model can have difficulty matching force plate data because of either (1) the type and number of contact model elements [29] or (2) the lack of subject-specificity in foot segment [34] or shoe stiffness [26] parameter values. While some studies have used optimization methods to determine foot-ground contact model parameter values [19,21,24], either the GRF results did not closely match the experimental data or ground reaction moment results were not reported. In addition, no foot-ground contact model results have been reported for the contralateral foot whose ground contact pattern is discontinuous with that of the ipsilateral foot. Thus, a three-dimensional subject-specific foot-ground contact model capable of reproducing all six ground reaction quantities (expressed as three reaction force components, two CoP coordinates, and a free reaction moment) would be valuable for computational studies that seek to predict new gait motions.

This study describes a novel optimization methodology for constructing subject-specific foot-ground contact models that can reproduce all six quantities of experimentally measured ground reactions for both feet during walking. Similar to previous studies, the method utilizes a two-segment foot model with a grid of viscoelastic elements positioned across the bottom of each segment. Unlike previous studies (see Table 1), it employs a three-stage optimization procedure to calibrate foot-ground contact parameter values using ground reaction data collected from the subject. Each stage either calibrates for the first time or improves calibration of additional parameter values in the model. The three stages involve matching (1) the vGRF, (2) all three components of GRF, and (3) all six ground reaction components (three forces and three moments), which can be expressed equivalently as three reaction force components, two CoP coordinates, and a free reaction moment, for both feet from a single walking trial. Throughout this paper, we will refer to these quantities as “ground reaction quantities”. The ability of the calibrated model to reproduce all six ground reaction quantities is evaluated using four additional walking trials excluded from the calibration process.

To develop and evaluate the proposed optimization methodology, we collected motion capture and ground reaction data from one healthy subject (male, age 46, height 1.7 m, weight 69 kg). The study was IRB approved and the subject gave informed consent. A 14-camera Vicon motion capture system (Vicon Motion Systems, Inc., Oxford, UK) measured surface marker positions, and three six-axis Bertec force plates (Type 4060-08, Bertec Corp., Columbus, OH) measured GRFs and moments. The subject was given Adidas Samba Classic sneakers to wear because the bottom is nearly flat (i.e., in a static pose, the vertical distance between ground and toe tip is less than 5 mm) and the insole is neutral with no built-in cushioning, thereby minimizing the complexity of the foot-ground contact model to be developed. Six markers (four dynamic and two static) were placed on each shoe for tracking motion and for defining two foot segments (HF and FF). The four dynamic markers were placed over the tip of the second toe, the heel at the same height as the toe marker, the superior central aspect of the HF over the laces, and the inferior lateral side of the HF just above the sole. The two static markers (removed for dynamic tasks) were placed on each shoe at the base of the first (medial toe marker) and fifth (lateral toe marker) metatarsals to define a single oblique metatarsal phalangeal (toes) axis.

Data collection consisted of static standing trials and dynamic overground walking trials. For the static trials, each foot was located on a separate force plate with the toes pointed forward. One static trial was collected to measure static marker positions on the two feet. Two additional static trials (one per foot) were collected where the shape of the sneaker sole was outlined on the force plate using a marker wand (Fig. 1). For the dynamic trials, we collected five walking trials with clean strikes on all three force plates using the subject's self-selected speed (1.4 m/s). Based on the order of force plate strikes, we defined a gait cycle for both feet as being from heel strike to subsequent heel strike of the right foot. Thus, for one gait cycle, complete force plate data were available from both feet. The walking trial with velocity closest to the mean of all five trials was selected for model calibration purposes. The ability of the calibrated model to reproduce all six ground reaction quantities was evaluated using the four remaining trials withheld from the calibration process.

We developed a parametric two-segment (HF and FF) foot-ground contact model using autolev symbolic manipulation software (Motion Genesis, Palo Alto, CA). The model possessed seven degrees-of-freedom (DoFs) consisting of three translations and three rotations defining the position and orientation of the HF segment in the lab frame (6DoFs) and one rotation defining toes flexion with respect to the HF segment (1DoF). The toes rotation axis was moved to the floor to eliminate any unrealistic gaps that could form when the toes flex during the transition between the flat foot and toe off portions of stance phase. The locations of static and dynamic markers on both segments and of contact elements on the bottom of the foot (see description below) were included in the model. The relevant kinematic and contact element force equations were incorporated into a matlab (The Mathworks, Natick, MA; version 2013 b) program, where all subsequent model development tasks were performed.

Static marker data were used to embed coordinate systems and the toes flexion axis in the HF and FF segments. The origin of the HF segment was defined using the heel marker position and the origin of the FF segment using the toe tip marker position. An anterior direction (x) was defined in both segments using the heel-to-toe tip marker direction projected onto the plane of the floor. A superior direction (z) in both segments was defined to coincide with the superior axis of the lab coordinate system. Medial (right foot) and lateral (left foot) directions (y) were defined as the cross product of the superior direction with the anterior direction. A toes flexion axis was defined in the HF and FF segments of each foot using the vector connecting the medial and lateral toe markers projected onto the floor in the static pose.

Static marker data were also used to construct a uniform rectangular grid of contact elements on the bottom of each foot model. The grid axes were aligned with the AP and ML axes of each model, and the grid boundaries were defined by the positions of the heel, toe tip, medial toe, and lateral toe markers projected onto the floor. Contact elements whose centers were located outside the outline of the sneaker sole were eliminated from the grid. This process reduced the number of grid elements from 55 (5 × 11 grid) down to 38. Remaining contact elements whose centers were anterior to the toes flexion axis was fixed in the FF segment, and those whose centers were posterior to this axis were fixed in the HF segment. To select a grid density, we choose the minimum number of elements (and thus the minimally complex model) that would match experimental ground reactions as well as or better than in previous studies. Based on trial and error with different grid densities, a 5 × 11 grid was the sparsest one that met this criterion.

Each contact element generated force in three orthogonal directions using a linear spring with nonlinear damping to calculate vertical force and a stick-slip friction model to calculate horizontal force. The vertical position and velocity of the center of any contact element i in the lab coordinate system were used to define penetration depth $yi$ and velocity $y˙i$ for calculating vertical force $Fi$Display Formula

In this equation, $ki$ is spring stiffness unique to each spring, $y0$ is spring resting length (effectively a change in vertical spring position within each foot segment) common to all springs, and $ci$ is a nonlinear damping coefficient unique to each spring. The nonlinear damping term modifies the elastic force and produces physically realistic hysteresis during spring compression [35], in contrast to linear damping which exhibits a force discontinuity at the transitions into and out of contact. Similarly, the horizontal velocity of the center of element i in the lab coordinate system was used to define slip velocity magnitude $vsi$ for calculating the horizontal friction force magnitude $fi$Display Formula

In this equation, $μs$, $μd$, and $μv$ are coefficients of static, dynamic, and viscous friction, respectively, common to all springs. The first term accounts for dynamic friction, the next two terms account for static friction, and the last term accounts for viscous friction. The latching velocity magnitude $vl$ defines the start of a linear transition region through zero velocity between static friction in one direction and static friction in the opposite direction, thereby avoiding a numerical discontinuity at the transition through zero slip velocity. Latching velocity is normally set to an arbitrarily small value, which was selected to be 0.01 m/s in this study. This friction model is derived from the Hollars friction model [36] with modifications to make it continuous and differentiable. The calculated horizontal friction force $fi$ was applied to the contact element in the direction opposite to the slip velocity vector. Once $Fi$ and $fi$ were calculated for each contact element, the net contact force and torque acting on each foot segment about its origin were calculated from the segment's individual contact element forces using the principle of replacement [37].

To simplify the model development process, we assumed that all model parameter values were identical between the two feet and investigated kinematic parameterization methods. This assumption necessitated small adjustments to the contact element grid for one foot so that an identical grid could be used for both feet. For clinical situations where one would expect the two feet to differ, one can decouple the parameter values for the two feet, which will increase the accuracy with which ground reactions can be reproduced by each foot-ground contact model. In addition, we investigated methods for parameterizing kinematic curves to facilitate making small changes to model kinematics during the optimization process. Instead of parameterizing joint position curves directly, we parameterized curves defining deviations away from the experimentally determined joint positions. Parameterization was achieved using B-spline curves with 25 coefficients. The number of coefficients and parameterization method were determined from a sensitivity analysis that used a range of 20–30 coefficients and either a B-spline or polynomial plus Fourier parameterization.

The spring, damper, and friction parameter values in the foot-ground contact model were calibrated for both feet simultaneously using a novel three-stage optimization procedure applied to the calibration walking cycle. All three stages employed Matlab's nonlinear least squares optimization algorithm and used an augmented cost function to account for “soft” constraints (see the Appendix for cost function details for each stage of the calibration process).

In the first stage, we sought to reproduce marker positions and vGRFs for both feet simultaneously by modifying design variables defining model kinematics (5DoFs; horizontal translations excluded), the common spring resting length $y0$, and the 38 spring stiffness $ki$ and damping $ci$ coefficients (i = 1–38). The cost function minimized deviations of individual $ki$ and $ci$ values about the mean values from all contact elements. Soft constraints were used to achieve acceptable marker position (<6 mm) and vGRF (<20 N) errors. The cost function also included terms that minimized errors in the first time derivatives of kinematic (vertical translation, three HF rotations, and toe flexion), marker position, and vGRF curves, which facilitated finding a solution that produced smooth vGRF curves. Model parameter values $ki$, $ci$, and $y0$ were limited to be greater than zero through the use of additional soft constraints in the cost function.

In the second stage, we sought to reproduce marker positions and all three components of GRF for both feet simultaneously by modifying design variables defining model kinematics (all 7DoFs), the three friction coefficients $μs$, $μd$, and $μv$, and the 38 spring stiffness $ki$ and damping $ci$ coefficients. Because of the sensitivity of the vGRF to the common spring resting length $y0$, $y0$ was fixed to the value found in the first stage. However, $ki$ and $ci$ values were allowed to change and their initial values were taken from the previous stage. The cost function was the same as that in the previous stage, with additional terms for matching AP and ML forces and their first derivatives, as well as additional soft constraints encouraging marker errors and friction coefficients to remain close to specified values. These values were 2.5 mm for the marker distance errors (we assume a little marker motion with small changes in kinematics), 0.2 for $μs$, 0.15 for $μd$, and 0.005 for $μv$. A cost function term was also added to enforce $μs$  >  $μd$.

In the third stage, we sought to reproduce marker positions and all six ground reaction quantities for both feet simultaneously by modifying the same design variables as in the second stage. Initial guesses for all design variables were taken from the previous stage. The cost function was the same as in the previous stage, with additional terms for matching CoP location and the free moment, as well as their first time derivatives.

To address the possibility that the three-stage optimization procedure found a local minimum, we repeated the entire calibration process ten times starting from different initial guesses. We chose the model parameter values from the solution that matched all six ground reactions the closest for subsequent model evaluation.

The calibrated model parameter values $ki$, $ci$, $y0$, $μs$, $μd$, and $μv$ produced by the third stage were evaluated for both feet using the four additional gait trials excluded from calibration. Evaluation involved repeating the third stage optimization with all model parameters fixed and only model kinematics (all 7DoFs) allowed to vary (see the Appendix for cost function details for the evaluation process).

We also investigated the sensitivity of model-predicted ground reactions and foot kinematics to variations in model parameter values. For these analyses, we varied in both directions all stiffness values by 20%, all damping values by a factor of ten, and all friction coefficients (each type separately) by a factor of 2. Sensitivity was quantified by dividing the percent change in RMS error of the output (ground reaction or kinematic curve) by the percent change in the parameter value.

The three-stage calibration process produced subject-specific foot-ground contact models that closely reproduced experimental GRF, CoP location, and free moment curves (Fig. 2) with only minor kinematic changes for both the right (Fig. 3) and left (Fig. 4) foot. The calibrated model parameter values were the same for both feet by design (Table 2). Mean $k$ values were 2596 ± 674 N/m and mean $c$ values were 2.75 × 10⁻⁶ ± 1.63 × 10⁻⁵ s/m. Calibrated friction coefficients were 0.101, 0.088, and 0.013 for static, dynamic, and viscous friction, respectively.

The calibrated models for both feet also achieved close agreement with experimental ground reaction and CoP curves for the four testing trials. RMS errors (Table 3) for GRF curves spanned 6.9–16.3 N for the right foot and 3.8–14.8 N for the left foot (Fig. 2). Predicted CoP locations deviated the most in the AP direction with RMS errors ranging from 21.7 to 27.7 mm for the right foot and from 16.0 to 21.7 mm for the left foot. In contrast, RMS errors in the ML direction ranged from 5.0 to 10.3 mm for the right foot and from 6.0 to 8.5 mm for the left foot. The RMS errors for the free moment curves were between 0.8 and 2.0 N/m for the right foot and between 1.4 and 1.8 N/m for the left foot. The free moment of both feet exhibited moderate sensitivity to the change in stiffness values, while all other ground reaction quantities exhibited little sensitivity to changes in model parameter values (Table 4).

Minimal changes in model kinematics relative to an inverse kinematics solution were needed for both the right (Fig. 3) and left (Fig. 4) foot to match experimental ground reactions well. For both calibration and evaluation trials, RMS errors were under 4.5 mm for all translational DoFs and under 1.5 deg for all rotational DoFs (Table 3 ). The necessary model kinematic changes were comparable for both feet, with the largest mean (of all five gait trials) RMS error being 3.9 mm for the right foot and 4.1 mm for the left foot for translational kinematics and 0.6 deg for both feet for rotational kinematics. The average RMS marker distance error for the HF was below 6 mm for both feet.

The goal of this study was to develop a three-dimensional subject-specific foot-ground contact modeling approach capable of reproducing all six quantities of experimental ground reactions for both feet. Subject-specific models of both feet developed with our approach generated force, CoP, and free moment curves that closely matched experimental curves. Only small variations in foot translational and rotational kinematics relative to inverse kinematic solutions were needed to achieve the desired results. Accurate models of foot-ground interaction may facilitate the development of patient-specific walking optimizations that can predict a patient's post-treatment walking function when large changes in foot-ground contact pattern occur.

It was challenging to develop a foot-ground contact model that could match all six experimental ground reaction curves. Compared to curves reported in the literature [18–20,22,24], the three forces produced by our model for both feet did not exhibit spikes, did not oscillate about the experimental curves, and more closely matched the experimental data. However, it was more difficult to match the CoP and free moment curves. There are at least two possible explanations for this finding. First, the CoP calculation requires dividing by the vGRF. Thus, when the vGRF was “small” at the start and end of stance phase, CoP errors were the largest, and consequently, the optimizer may have been trying to match inaccurate CoP locations. For this reason, we excluded CoP values from the calibration process results (Fig. 2 , Table 3) when vGRF was less than 100 N (about 10% of the maximum vGRF value). Regardless, the influence of CoP prediction errors at the start and end of stance phase on joint moments predicted by a full-body walking model will be minimal since the GRF vector is small at these times. Second, the CoP calculation (and thus the free moment calculation as well) can be influenced by small inaccuracies in force plate measurements. This possibility is supported by two observations. First, the CoP and free moment errors were much larger for the right foot than for the left. Second, when we performed CoP accuracy tests on each force plate using a CalTester stick (C-Motion, Inc., Germantown, MD), we observed differences between actual and measured CoP locations of up to 1 cm for two of the three force plates. While the small magnitudes of the CoP and free moment curves make differences between model and experiment seem large, the model free moment curves exhibited the same trends as the experimental free moment curves. This finding is important because the free moment is critical for counteracting the angular momentum produced by arm swing during walking [25].

We chose to report errors in ground reaction quantities using absolute units rather than percentage of expected value, since absolute error is the quantity that will influence simulation accuracy. As an example, consider a situation where the expected AP force is 0 N, the expected vertical force is 800 N, and there is 1 N of error in the AP force. Then the percent error in AP force would be infinite (1/0 × 100), but an infinite percent error would be misleading since it would not reflect the small size of the error relative to the magnitude of the GRF vector. In a simulation with an 800 N ground contact force vector that is primarily in the vertical direction, 1 N of AP force error (regardless of what the actual AP force value is) will be quite acceptable.

We chose to calibrate both feet together since there was no anatomic reason to believe that the two feet should have different parameter values for this subject. However, our model calibration approach can also be applied to each foot individually. When we tried calibrating each foot separately, we obtained somewhat different model parameter values for the two feet but only slightly better ground reaction and marker position matching. This finding suggests that use of a single set of parameter values for both feet may constrain the solution in a physically realistic manner without degrading the results appreciably. It also suggests that other combinations of model parameter values may yield reasonable ground reaction predictions.

To investigate the sensitivity of model-predicted ground reactions and foot kinematics to variations in model parameter values, we performed a series of sensitivity analyses. All ground reaction quantities exhibited little sensitivity to changes in model parameter values, except for the free moment of both feet, which exhibited moderate sensitivity to the change in stiffness values (Table 4 ). Interestingly, the free moment curves were the hardest to match, likely due in part to their sensitivity to stiffness parameter changes, which may be why few studies report them. Overall, the sensitivity study results suggest that while the optimized model parameter values may remain nonunique, they are still able to predict all six ground reaction quantities well for both feet.

We converged on a three-stage optimization methodology through a trial-and-error process. When we initially tried optimizing all six ground reaction quantities directly, the optimizer was unable to find a good solution for any quantity. Since friction force depends on the vGRF and the CoP location and free moment depend on all forces, it is logical that a closely matched vGRF (first stage) helps the optimizer find a solution that matches friction forces (second stage) and ultimately all ground reaction quantities (third stage) well. By breaking the problem into three stages that used the optimized parameter values from the previous stage as initial guesses, we were able to obtain a good solution for all six ground reaction quantities simultaneously for both feet.

We also tried using different optimization algorithms before choosing a nonlinear least squares algorithm. Several of Matlab's optimization algorithms were investigated during development of the calibration process. The constrained nonlinear optimization algorithm fmincon worked well for matching the vGRF. However, when used for the second stage, that algorithm took much longer to find a less optimal solution than did the nonlinear least squares algorithm, and it was unable to find a feasible solution when the moments were included in the cost function for the third stage.

Some parameters were more important than others for finding an optimal solution. To generate force, the spring-damper elements need to penetrate the floor. Penetration could be accomplished by pushing the whole foot into the floor, leading to large marker errors, or by offsetting each contact element within its respective segment frame. For uniformity, we chose to offset all contact elements by the same amount. Similarly, we tried adjusting the height of the elements in the FF segment so that it varied along a parabolic surface to mimic the bottom of sneakers more closely. However, the parabolic adjustment competed with the vertical offset value and often led to a slope near zero. Therefore, we removed this adjustment from the problem formulation, assumed the bottom of the shoe to be flat, and allowed small increases in toes flexion to make up for this geometric omission. The parabolic adjustment could turn out to be valuable for shoes without such a flat bottom.

While our model calibration approach generated good solutions for ground reactions and foot kinematics, our study was not without limitations. First, our method was tested on only one subject with one pair of sneakers, though we demonstrated that the model could generate accurate ground reactions and foot kinematics for both feet over multiple walking trials excluded from calibration. Second, we calibrated and tested our method using only walking data collected at a single walking speed. However, we do not believe that this limitation is important, since our ultimate goal is to predict new walking motions. If we also wanted to predict foot-ground contact patterns for other tasks such as stair climbing or jumping, then use of only walking data for model calibration would be more limiting. Third, our model did not account for soft tissues, ligaments, or muscles that cushion the foot and aid in the coordination of walking. Fourth, we did not model foot motion relative to the shoe, which may have decreased the accuracy of our model. Since we were seeking the simplest model possible that can predict ground reactions well, we modeled the foot and shoe as rigidly connected (though as separate FF and HF segments), similar to the approach used in previous studies [18,32]. Fifth, we evaluated our method using only one size of contact element grid. While our results were good for the chosen 5 × 11 grid of contact elements, further studies should investigate whether other grid densities could lead to improvements in accuracy and/or computation time.

In the present study, we did not evaluate the ability of our foot-ground contact model to simulate slower walking speeds, as might be expected in a clinical context, or to work well for faster walking speeds when calibrated using a slower walking speed. To address these important issues, we applied our modeling approach to additional walking data collected from a hemiparetic subject walking at 0.5 and 0.8 m/s. We calibrated the parameter values in our foot-ground contact models using walking data collected at 0.5 m/s (same parameter values for both feet) and tested the calibrated models using walking data collected at 0.8 m/s. We found that the RMS errors in predicted ground reaction quantities at both walking speeds were comparable to the values reported in the present study, supporting the use of our foot-ground contact models for slower walking speeds and walking speeds that differ from the calibration speed.

The goal of this study was to develop a subject-specific foot-ground contact modeling methodology that reproduces all six experimental ground reaction curves for both feet. Our modeling approach successfully reproduced experimental ground reactions from a calibration trial and four trials excluded from calibration. For predictive optimizations of three-dimensional walking, all six ground reaction quantities for both feet are required. Such optimization may one day become a valuable tool for developing customized rehabilitation strategies that optimize a patient's post-treatment function. Our foot-ground contact modeling methodology provides valuable functionality for subject-specific predictive gait optimizations that may eventually lead to improvements in rehabilitation medicine.

This study was funded by NSF Grant Nos. CBET 1052754 and CBET 1159735. Additional support was provided by the National Institutes of Health.

Stage 1 both feetDesign variablesCost function termsCost function weights38 Kvals1000*FootError/(initial_marker_errors + tol)R = L = 238 cvalsmarker_slope_diff (coordinates for four foot markers)R = L = (1/0.001)25*5 BSP coefficients (right foot)Q_slope_diff (Q2,5–7)R = L = (1/0.001)25*5 BSP coefficients (left foot)GRFdiff (vGRF only)R = L = 1YvalR (m)GRF_slope_diff (vGRF only)R = L = (1/5)YvalL (m)K_mean_diff1/100c_mean_diff100Virtual springs for 38 cvals: ((cvals − initial_cs)/initial_cs)¹⁰1Virtual springs for 38 Kvals: ((Kvals − initial_Ks)/1)¹⁰1Virtual spring for YvalR: ((YvalR − 0.01)/0.01)¹⁰1Virtual spring for YvalL: ((YvalL − 0.008)/0.008)¹⁰1^*The total cost was divided by 10.

Definitions

Kvals = 38 spring stiffness values (one for each active element; mirrored for second foot)

cvals = 38 spring damping values (one for each active element; mirrored for second foot)

BSP coefficients = used to parameterize deviation curves that were added to original kinematic curves (25 coefficients per curve per foot; total = 250)

Yval = the vertical offset between the ground and the spring-damper elements (only one value per foot for all the elements; total = 2)

FootError = difference between the experimental and model marker coordinates (m); multiplied by 1000 to convert to units of millimeter; 12 values per foot (x, y, z coordinates for each of the four foot markers); total = 24

initial_marker_errors = max coordinate marker errors determined from Soderqvist and Wedin in each direction for each marker (12 per foot; total = 24)

tol = two for hindfoot marker and three for toe marker coordinates

marker_slope_diff = difference between the slopes of the experimental and model marker coordinate curves (12 curves per foot; total = 24)

Q_slope_diff = difference between slopes of kinematic curves (five curves per foot; total = 10)

GRFdiff = difference between the experimental and model GRF curves

GRF_slope_diff = difference between the slopes of the experimental and model GRF curves

K_mean_diff = difference between the Kvals and the mean of all the Kvals (38 total)

c_mean_diff = difference between the cvals and the mean of all the cvals (38 total)

initial_Ks = starting K values obtained from a right foot only optimization yielding a vGRF curve closely matching the experimental vGRF curve

initial_cs = starting c values obtained from a right foot only optimization yielding a vGRF curve closely matching the experimental vGRF curve

Stage 2: This stage matches all three GRF curves and tweaks all the seven kinematic curves. Stage 2 both feetDesign variablesCost function termsCost function weights38 Kvals1000*FootError/(initial_marker_errors + tol)R = L = 138 cvalsxGRFdiffR = 3; L = 525*7 BSP coefficients (right)yGRFdiffR = 2; L = 525*7 BSP coefficients (left)zGRFdiffR = 2; L = 5mu_sxGRF_slope_diffR = L = (1/3)mu_dyGRF_slope_diffR = L = (1/5)mu_vzGRF_slope_diffR = L = 2marker_slope_diff (coordinates for four foot markers)R = L = (1/0.001)Q_slope_diff (all seven Q’s)R = L = (1/0.0001)K_mean_diff(1/10)K_err(1/100)c_mean_diff100Virtual springs for marker dist errors: ((FootDistRMSErrors − 2.5)/2.5)²⁰R = L = 1Virtual springs for 38 cvals: ((cvals − initial_cs)/1)¹⁰1Virtual springs for 38 Kvals: ((Kvals − initial_Ks)/initial_Ks)¹⁰1Virtual spring for mu_s: ((mu_s − 0.25)/0.05)¹⁰1Virtual spring for mu_d: ((mu_d − 0.2)/0.05)¹⁰138 virtual spring for mu_v: ((mu_v − 0.005)/0.005)²⁰1Virtual spring for mu’s: ((mu_s − mu_d) − 0.06)/0.06)¹⁰1^*The total cost was divided by 50.

Definitions

Kvals = 38 spring stiffness values (one for each active element; mirrored for second foot)

cvals = 38 spring damping values (one for each active element; mirrored for second foot)

BSP coefficients = used to parameterize deviation curves that were added to original kinematic curves (25 coefficients per curve per foot; total = 350)

mu_s = static coefficient of friction

mu_d = dynamic coefficient of friction

mu_v = viscous coefficient of friction

FootError = difference between the experimental and model marker coordinates (m); multiplied by 1000 to convert to units of millimeter; 12 values per foot (x, y, and z coordinates for each of the four foot markers); total = 24

initial_marker_errors = max coordinate marker errors determined from Soderqvist and Wedin in each direction for each marker (12 per foot; total = 24)

tol = two for hindfoot marker and three for toe marker coordinates

GRFdiff = difference between the experimental and model GRF curves in the x, y, and z directions

GRF_slope_diff = difference between the slope of the experimental and model GRF curves in the x, y, and z directions

marker_slope_diff = difference between the slopes of the experimental and model marker coordinate curves (12 curves per foot; total = 24)

Q_slope_diff = difference between the slopes of kinematic curves (seven curves per foot; total = 14)

K_mean_diff = difference between the Kvals and the mean of all the Kvals (38 total)

K_err = difference between each K value and its starting value

c_mean_diff = difference between cvals and the mean of all the cvals (38 total)

FootDistRMSErrors = RMS foot marker distance errors (mm); one per marker, four markers per foot (total = 8)

initial_cs = starting c values obtained from stage 1 calibration

initial_Ks = starting K values obtained from stage 1 calibration

Stage 3: This stage matches all three GRF curves, the center of pressure location, and the free moment curve. Stage 3 both feetDesign variablesCost function termsCost function weights38 Kvals1000*FootError/
(initial_marker_errors + tol)R = L = 138 cvalsxGRFdiff (right)3 (tf = 1)25*7 BSP coefficients (right)3 (tf = 2:63 (stance))25*7 BSP coefficients (left)30 (tf = 64:101)mu_sxGRFdiff (left)3mu_dyGRFdiff (right)3 (tf = 1)mu_v6 (tf = 2:63 (stance))30 (tf = 64:101)yGRFdiff (left)3zGRFdiff (right)3 (tf = 1)3 (tf = 2:63 (stance))30 (tf = 64:101)zGRFdiff (left)3xGRF_slope_diffR = 2; L = (2/3)yGRF_slope_diffR = (2/5); L = (3/5)zGRF_slope_diffR = 2; L = 4marker_slope_diff (coordinates for four foot markers)R = L = (1/0.001)Q_slope_diff (all seven Q’s)R = L = (1/0.0001) for Q1–3 and (1/(0.0001*0.0873)) for Q4–7K_mean_diff1/50C_mean_diff10,0001000*xCOPdiffR = L = 81000*zCOPdiffR = L = 8FTdiffR = L = 8xCOP_slope_diffR = L = (2/0.01)zCOP_slope_diffR = L = (2/0.01)FT_slope_diffR = L = (2/0.01)K_err1cvalneg1Virtual springs for marker dist errors: ((FootDistRMSErrors−2.5)/2.5)¹⁰R = L = 1Virtual springs for 38 cvals: ((cvals−initial_cs)/1)¹⁰1Virtual springs for 38 Kvals: ((Kvals−initial_Ks)/initial_Ks)¹⁰1Virtual spring for mu_s: ((mu_s−optval)/0.02)A¹⁰1Virtual spring for mu_d: ((mu_d−optval)/0.02)A¹⁰138 virtual springs for mu_v: ((mu_v−optval)/0.003)²⁰1Virtual spring for mu's: ((mu_s−mu_d)−0.05)/0.05)¹⁰1000^*The total cost was divided by 250.

Definitions

Kvals = 38 spring stiffness values (one for each active element; mirrored for second foot)

cvals = 38 spring damping values (one for each active element; mirrored for second foot)

BSP coefficients = used to parameterize deviation curves that were added to original kinematic curves (25 coefficients per curve per foot; total = 350)

mu_s = static coefficient of friction

mu_d = dynamic coefficient of friction

mu_v = viscous coefficient of friction

FootError = difference between experimental and model marker coordinates (m); multiplied by 1000 to convert to units of mm; 12 values per foot (x, y, and z coordinates for each of the four foot markers); total = 24

initial_marker_errors = max coordinate marker errors determined from Soderqvist and Wedin in each direction for each marker (12 per foot; total = 24)

tol = two for hindfoot marker and three for toe marker coordinates

GRFdiff = difference between the experimental and model GRF curves in the x, y, and z directions

GRF_slope_diff = difference between the slope of the experimental and model GRF curves in the x, y, and z directions

marker_slope_diff = difference between the slopes of the experimental and model marker coordinate curves (12 curves per foot; total = 24)

Q_slope_diff = difference between the slopes of kinematic curves (seven curves per foot; total = 14)

K_mean_diff = difference between the Kvals and the mean of all the Kvals (38 total)

c_mean_diff = difference between the cvals and the mean of all the cvals (38 total)

COPdiff = difference between the experimental and model COP curves in the x and z directions

FTdiff = difference between the experimental and model free moment curves

COP_slope_diff = difference between the slopes of the experimental and model COP curves in the x and z directions

FT slope diff = difference between the slopes of the experimental and model free moment curves

Testing: This stage tests the values for the design variables that were optimized with the three-stage calibration procedure. Testing for both feetDesign variablesCost function termsCost function weights25*7 BSP coefficients (right)1000*FootError/(initial_marker_errors + tol)R = L = 125*7 BSP coefficients (left)xGRFdiff (right)5 (tf = 1)5 (tf = 2:63 (stance))50 (tf = 64:101)xGRFdiff (left)5yGRFdiff (right)2 (tf = 1)4 (tf = 2:63 (stance))20 (tf = 64:101)yGRFdiff (left)2zGRFdiff (right)5 (tf = 1)5 (tf = 2:63 (stance))50 (tf = 64:101)zGRFdiff (left)5xGRF_slope_diffR = (1/3); L = 1yGRF_slope_diffR = (1/5); L = (1/5)zGRF_slope_diffR = 1; L = 2marker_slope_diff (coordinates for four foot markers)R = L = (1/0.001)Q_slope_diff (all seven Q’s)R = L = (1/0.0001)1000*xCOPdiffR = L = 51000*zCOPdiffR = L = 5FTdiffR = L = 20xCOP_slope_diffR = L = (1/0.01)zCOP_slope_diffR = (1/0.001); L = (1/0.01)FT_slope_diffR = (1/0.001); L = (1/0.01)Virtual springs for marker dist errors: ((FootDistRMSErrors−2.5)/2)¹⁰R = L = 1^*The total cost was divided by 100.

Definitions

BSP coefficients = used to parameterize deviation curves that were added to original kinematic curves (25 coefficients per curve per foot; total = 350)

FootError = difference between the experimental and model marker coordinates (m); multiplied by 1000 to convert to units of millimeter; 12 values per foot (x, y, and z coordinates for each of the four foot markers); total = 24

initial_marker_errors = max coordinate marker errors determined from Soderqvist and Wedin in each direction for each marker (12 per foot; total = 24)

tol = two for hindfoot marker and three for toe marker coordinates

GRFdiff = difference between the experimental and model GRF curves in the x, y, and z directions

GRF_slope_diff = difference between the slope of the experimental and model GRF curves in the x, y, and z directions

marker_slope_diff = difference between the slopes of the experimental and model marker coordinate curves (12 curves per foot; total = 24)

Q_slope_diff = difference between the slopes of kinematic curves (seven curves per foot; total = 14)

COPdiff = difference between the experimental and model COP curves in the x and z directions

FTdiff = difference between the experimental and model free moment curves

COP_slope_diff = difference between the slopes of the experimental and model COP curves in the x and z directions

FT_slope_diff = difference between the slopes of the experimental and model free moment curves

FootDistRMSErrors = RMS foot marker distance errors (mm); one per marker, four markers per foot (total = 8)