Introduction: Simulations based on computational musculoskeletal models are powerful tools for evaluating the 24 effects of potential biomechanical interventions, such as implementing a novel prosthesis. However, the utility of 25 simulations to evaluate the effects of varied prosthesis design parameters on gait mechanics has not been fully 26 realized due to lack of a readily-available limb loss-specific gait model and methods for efficiently modeling the 27 energy storage and return dynamics of passive foot prostheses. The purpose of this study was to develop and 28 validate a forward simulation-capable gait model with lower limb loss and a semi-active variable-stiffness foot 29 (VSF) prosthesis. Methods: A seven-segment 28-DoF gait model was developed and forward kinematics 30 simulations, in which experimentally-observed joint kinematics were applied and the resulting contact forces under 31 the prosthesis evolved accordingly, were computed for four subjects with unilateral below-knee amputation walking 32 with a VSF. Results: Model-predicted resultant ground reaction force (GRFR) matched well under trial-specific 33 optimized parameter conditions (mean R2: 0.97, RMSE: 7.7% body weight (BW)) and unoptimized (subject-34 specific, but not trial-specific) parameter conditions (mean R2: 0.93, RMSE: 12% BW). Simulated anterior-posterior 35 center of pressure demonstrated a mean R2 = 0.64 and RMSE = 14% foot length. Simulated kinematics remained 36 consistent with input data (0.23 deg RMSE, R2 > 0.99) for all conditions. Conclusions: These methods may be 37 useful for simulating gait among individuals with lower limb loss and predicting GRFR arising from gait with novel 38 VSF prostheses. Such data are useful to optimize prosthesis design parameters on a user-specific basis.39 Manuscript body 2 1. Introduction 40 Individuals with below-knee amputation who use prostheses for ambulation exhibit distinct gait characteristics, 41 which may limit mobility and decrease quality of life. They may display gait asymmetry [1,2], elevated metabolic 42 cost during locomotion [3], increased muscle activity [4], and a variety of psychological disorders including anxiety 43 and depression [5]. Sustained prosthesis use may also lead to musculoskeletal injury through overloading of intact 44 joints [6] and residual limb tissue damage due to pressure and shear forces at the prosthetic socket interface [7–9]. 45 Each of these issues may be attenuated by improving user specificity in the design characteristics of foot prostheses. 46 However, the effects of foot prosthesis design parameters (e.g. stiffness) are not well characterized, and thus 47 achieving meaningful improvements in gait through prosthesis design has proven elusive [10,11]. In order to achieve 48 improvements, a robust understanding of the relationships among anatomical morphologies associated with lower 49 limb amputation (LLA), gait mechanics, and prosthesis design is necessary. 50 Gait simulations based on computational musculoskeletal models can be powerful tools to investigate human 51 movement. Simulations have been used previously to evaluate musculoskeletal dynamics under a variety of 52 conditions, including walking [12], running [13], and jumping [14]. However, the models used in these studies are 53 not suitable for simulating gait among individuals with LLA, due to the altered lower limb mass, inertial properties 54 and mechanical properties associated with a prosthetic foot and socket. LaPré et al. (2018) [15] and Willson et al. 55 (2020) [16] developed limb loss-specific models in OpenSim [12], a popular musculoskeletal modeling and 56 simulation platform. These studies used simulations to investigate the effects of prosthesis alignment and a 57 biarticular clutched spring mechanism on gait mechanics, respectively. However, the models used in these studies 58 do not account for the energy storage and return (ESR) properties of the prosthetic foot, thus limiting their 59 ecological validity. Other studies that did incorporate the force and torque contributions of ESR feet into gait models 60 focused on characterizing biomechanical and muscle activation responses with prosthesis use, but did not report 61 systematic validation of the underlying gait and prosthesis models [17,18]. While these studies made important 62 progress toward investigating the relationship among anthropometry, gait mechanics, and prosthetic foot design, 63 they had limited ability to verify simulation results through comparison with experimental data. Due to these 64 limitations, the use of simulations to inform the design of ESR foot prostheses has not been fully realized. 65 Previously, a reduced order computational model of a semi-active variable-stiffness foot (VSF) prosthesis [19] 66 was developed and validated in [20]. However, there remains a need to integrate this model into a scalable forward 67 kinematics and forward dynamics simulation capable gait model. Such a model would be useful for simulating gait 68 and elucidating the relationships between ESR prosthesis design and gait mechanics on an individualized basis. The 69 purpose of this study was to develop and validate a forward simulation-capable gait model with lower limb loss and 70 a semi-active VSF. 71 72 2. Methods 73 74 2.1 Gait model design 75 A seven-segment, 28 degree-of-freedom (DoF) gait model (Figure 1) was developed in Simscape Multibody 76 (Mathworks, Inc., Natick, MA). Simscape was chosen for its application to model-based design, customizability in 77 simulation computing (e.g. differential equation solver profiles, numerical integration settings, and tolerance 78 controls), and rich libraries of contact models, mechatronic components, and signal processing modules. Further, its 79 integration with the rest of the Mathworks computational suite provides a link by which complex algorithms and 80 data processing tools can be easily incorporated into the model. 81 The anatomical structures for the generic model were derived from open-source cadaveric skeletal 3D surface 82 geometry data from Mitsuhashi et al. (2009) [21]. A massless head and torso (HAT) segment is included in the 83 model for future uses; however, the present study utilized only kinematics of the pelvis and lower limbs. The pelvis 84 interacts with the world frame via a 6-DoF bushing joint, which is the model’s parent joint. For this study, pelvis 85 translational movement in the frontal plane was prescribed using the experimental gait data described below in order 86 to account for subjects that walked diagonally or performed minor turning during the gait trial. Other translational 87 movements were unconstrained. Pelvis rotations in three planes with respect to the world frame were prescribed 88 using experimental gait data. Each joint of the lower limb is modeled as a 3-DoF gimbal joint, with the coronal and 89 axial rotations of the knee constrained to have no movement. Rotational movements in each plane may be 90 constrained through user-defined inputs. To represent a below-knee amputation (BKA), the tibia and fibula were 91 transected via a planar cut and encapsulated within a prosthetic socket. The residual limb and prosthetic socket were 92 connected via a compliant bushing joint to assess 6 DoF socket-residual limb interface dynamics, similar to previous 93 Manuscript body 3 work by LaPrè et al. (2018) [15]. The rotational and translational stiffness as well as displacement and velocity 94 constraints were designed according to previous gait experiments [15] and finite element analysis [22]. The 95 translational aspect of the bushing joint has a stiffness of 20 N/mm and damping of 10 N·s/mm. Vertical 96 displacement of the residual limb with respect to the prosthetic socket (i.e. socket pitoning) is constrained to 35 mm 97 maximum displacement and 100 mm/s maximum displacement velocity. The rotational aspect of the joint has a 98 stiffness of 10 N·mm/deg and damping of 5 (N·mm)·s/deg The VSF model [20] is rigidly attached to the socket via a 99 pylon and pyramid adapter using weld joints. Separate generic models with right and left side amputations were 100 developed for this study and are available in the supplementary materials. 101 Methods for modeling the mechanical stiffness and ground contact dynamics of the VSF were previously 102 described in [20]. Briefly, the VSF model was designed using the lumped parameter approach for approximating 103 flexible body dynamics. The keel of the VSF is discretized into 16 segments connected with alternating revolute 104 joints and weld joints, yielding eight DoF. The most posterior segment is 66 mm in length, which matches the 105 minimum possible fulcrum position. The rest of the keel consists of 11.64-mm segments for a total beam length of 106 229 mm. The stiffness and damping values for the revolute joints were parameterized to represent the material 107 properties of the VSF’s G10/FR4 composite keel. A custom MATLAB script controls continuous fulcrum position 108 (i.e. variable stiffness) between the posterior position at 66 mm and a maximum forward position of 151 mm. 109 Stiffness of the VSF model was characterized by simulating static compression tests and optimizing the load-110 displacement response based on that of the physical VSF across the full range of fulcrum positions [19]. 111 Segment masses were estimated according to DeLeva (1996) [23]. Mass distribution and inertial properties of 112 the lower limbs and pelvis are modeled as conical frusta and an ellipsoid, respectively. The geometry of the conical 113 frustum for each limb is defined using the proximal and distal joint radii for each segment, as derived from marker 114 coordinates (see “Model scaling and parameterization”). The ellipsoid’s radii are defined based on the distances 115 between the right and left anterior and posterior superior iliac spines. The estimated mass of the HAT segment was 116 added to the mass of the pelvis. 117 The model is designed with variant subsystems, which are Simscape coding structures that allow for multiple 118 implementations of code where only one implementation is active during a simulation. In this design, gait 119 simulations can be computed within a forward kinematics framework (i.e. joint kinetics and end effector positions 120 estimated given joint kinematics as inputs) or forward dynamics framework (i.e. joint kinematics and end effector 121 positions estimated given joint torques as inputs). Implementation of the variant subsystems (i.e. whether the model 122 computes a forward kinematics or forward dynamics simulation) is controlled via user-defined inputs. 123 124 2.2 Model scaling and parameterization 125 Model scaling and parameterization are completed through a pipeline of custom MATLAB scripts (Figure 2). 126 The pipeline uses raw marker coordinates from optical motion capture trials of subjects walking with the VSF to 127 programmatically generate a scaled subject-specific gait model, parameterize the model, and drive forward gait 128 simulations. First, data from a static capture trial are used to scale and assemble the generic gait model according to 129 subject-specific anthropometrics. Estimates for segment dimensions, mass, and inertial properties are derived from 130 the static data [23]. Then, the subject-specific gait model is used for subsequent forward kinematics or forward 131 dynamics simulations to determine the movements or forces of the model based on inputs of joint angles or joint 132 torques, respectively. Custom methods for signal processing and derivation of joint angles (e.g. Cardan-Euler 133 sequences or quaternions) may be defined by the user. A variety of MATLAB-based approaches can be used to 134 solve the model. These include both continuous and discrete options, an extensive set of ordinary differential 135 equation (ODE) solvers, and a variety of numerical integration settings (e.g. fixed or variable time steps, step size 136 constraints, and tolerance values). Each of these parameters may be strategically selected based on their suitability 137 for a given simulation scenario, such as system dynamics, solution stability, computation speed, or solver 138 robustness. The specific data processing and computational configurations used for simulations in this study are 139 described in Section 2.4. 140 141 2.3 Contact model design and parameterization 142 Foot-ground contact models were designed for the VSF and the intact foot. A model of VSF-ground contact 143 consists of 24 sphere-to-plane contact models [24] parameterized to represent the geometry and dynamics of the 144 VSF’s foam base. On the intact side, contact spheres were placed on anatomical locations that experience high 145 localized pressure during stance phase (Figure 1) [25]. Methods for the design and parameterization of the contact 146 model were previously described in [20]. Briefly, each contact sphere estimates normal (Fn) (eq. 1) and frictional 147 (Ff) (eq. 2) forces associated with the collision of a viscoelastic sphere (a massless spring and damper system) and a 148 Manuscript body 4 rigid plane [26,27]. For the intact foot, each contact sphere is parameterized independently according to Equations 1-149 2. 150 151 !!= $(&∗(!)++(,∗()̇(>0,(̇>0&∗((>0,(̇<00(<0 (1) 152 Fn: normal force 153 k: contact stiffness 154 2: penetration depth 155 n: penetration exponent 156 y: damping force scaling factor 157 3: contact damping coefficient 158 159 !"= 4!!∗5#$%$&'6()'<6$*+,#*)-.!!∗5/&!,$&'6()'>6$*+,#*)-. (2) 160 Ff: frictional force 161 7: coefficient of friction 162 8012: velocity at point of contact 163 8345674189: velocity threshold 164 165 For the VSF, the overall foot contact model was divided into five zones; the sphere-to-plane models were 166 parameterized by zone (Figure 1). The heel of the VSF model is comprised of three zones; this choice was motivated 167 by the sensitivity of contact parameters when few spheres are in contact with the walking plane (e.g. the heel of the 168 foot early in stance phase). Contact parameters for each zone were optimized to match the experimental ground 169 reaction force data from a sample trial using procedures described below. The foam base of the physical VSF 170 undergoes compression throughout stance phase. To account for these effects, a modified Kelvin-Voigt nonlinear 171 spring and damper force law (eq. (1)) was implemented to represent contact between the VSF and walking plane 172 [26,27]. 173 174 2.4 Gait simulations 175 The data presented in this study were derived via the processing and simulation pipeline described above. 176 Forward kinematics gait simulations were computed for four subjects (one female) with a unilateral BKA (Table 1) 177 walking with the VSF configured to “high”, “medium”, and “low” stiffness settings for three trials each (Table 2). 178 To be included, participants must have been at least 2 years post-amputation and able to safely complete nine over-179 ground gait trials walking at 1.2 ± 0.1 m/s. Prior written informed consent was provided by all subjects as approved 180 by the Health Sciences Institutional Review Board at the University of Wisconsin-Madison. All simulations were 181 computed in Simscape Multibody using the ode15s solver profile with variable step sizes for numerical integration. 182 Motion capture (12-camera Optitrack Prime13 system, Natural Point, Inc. Corvallis, OR) and ground reaction 183 force (GRF) data (Bertec Inc. Columbus, OH) collected from subjects during nine gait trials (Table 2), a static 184 neutral pose, and functional joint movements [28] that were used to generate scaled subject-specific models and 185 compute forward kinematics gait simulations. Marker coordinates and GRF data were sampled at 200 and 1000 Hz, 186 respectively. Marker data were low-pass filtered (4th order Butterworth, ƒc: 6 Hz). Force data were down-sampled 187 (polynomial interpolation, ƒs: 200 Hz) and low-pass filtered (4th order Butterworth, ƒc: 40 Hz). Segment and joint 188 kinematics were estimated using Cardan Euler rotation calculations in accordance with International Society of 189 Biomechanics recommendations [29–31]. These data were then used to drive the corresponding joints in the model 190 during simulations. Localized joint center coordinates were estimated for the hips, knees, and intact ankle joint using 191 functional joint movements and methods described by Schwartz and Rozumalski (2005) [28]. The location of the 192 pylon-pyramid adapter interface (i.e. prosthetic “ankle”) was calculated based on a measured offset from markers 193 placed on the VSF. This interface was assumed to be rigid. Motion at the socket-limb interface was considered to be 194 passive based on the aforementioned velocity and displacement constraints. Residual limb length was estimated as 195 the distance from the knee joint center to the base of the prosthetic socket and scaled accordingly. The mass and 196 inertial properties of the residual limb were estimated first by deriving estimated density of the intact limb modeled 197 as a conical frustum with an assigned mass estimated per De Leva (1996). The derived density of the intact limb was 198 applied to the residual limb model, which was also modeled as a conical frustum. This limb was then truncated at 199 the respective level of amputation for each subject. 200 Manuscript body 5 Contact model parameters were optimized on a subject-specific basis. Contact model parameterization was 201 formulated as a least-squares optimization problem with the objective of minimizing the sum of squared errors 202 between model-predicted and experimentally measured resultant ground reaction force (GRFR). For a description of 203 the comparison between model-predicted and experimentally-measured data, see Section 2.5. A Latin Hypercube 204 Sampling-based optimization algorithm was used to programmatically derive the stiffness, damping, and friction 205 terms for each sphere. For each subject, contact model-derived GRFR prediction was optimized using data from one 206 medium stiffness trial. The GRFR error resulting from this trial represents the theoretical optimal performance of the 207 comprehensive VSF-ground contact model. The transferability of the optimized parameter values was determined by 208 simulating the two remaining medium stiffness trials and the three remaining trials each for the low and high 209 stiffness configurations. No GRF data were available to optimize or verify contact model parameters on the intact 210 side. As such, stiffness and damping for contact spheres on the intact side were parameterized to simplify the model 211 of the intact side, so that the contact spheres were able to support the model’s mass throughout stance phase. 212 213 2.5 Gait model evaluation 214 The gait model’s performance was evaluated under two scenarios: static and dynamic gait conditions. Of the 215 three trials per condition (Table 2), Subject 2 did not complete one medium stiffness trial, and Subject 4 did not 216 complete one high stiffness trial. In total, forward kinematics simulations were computed for 38 trials (4 static, and 217 34 dynamic). For the static condition, the model was simulated with anatomically neutral joint angles for ten 218 seconds. Model-predicted GRFR was measured at the end of the trial and compared to the mass of the subject. This 219 comparison represents the accuracy of the contact parameters to estimate the ground contact force imparted by the 220 subject with no dynamic component. For the dynamic conditions, forward kinematics gait simulations were 221 computed for the nine gait trials per subject with high, medium, and low VSF stiffness configurations. 222 The Simscape solver computes the model’s initial conditions by finding state values at the initial time step that 223 exactly satisfy the model’s system of equations. For all simulations in this study, the model conditions were assumed 224 to start from steady state (i.e. zero derivative). The model was initialized in a posture taken from the experimental 225 motion capture data at the first time step. The model’s initial position in the gait environment was derived based on 226 the initial coordinates of a virtual motion capture marker placed at the pelvis’ origin frame. Then, joint kinematics, 227 also from experimental data, were applied at each joint and the model was simulated forward in time, allowing 228 kinematic movements and contact with the ground to evolve in time according to the model’s dynamics. Continuous 229 differential equations were integrated with respect to time to compute all model variables as a function of time. A 230 forward kinematics-based approach was selected for this study because the goal was to apply experimentally-231 observed joint kinematics and allow the resulting contact forces under the prosthesis to evolve accordingly. 232 Joint kinematics and GRFR data derived from simulations were low-pass filtered (4th order Butterworth: fc: 6 Hz 233 and 40 Hz, respectively). Simulation and experimental GRFR were time locked and indexed to 0.25 s before and 234 0.25 s after stance phase. Including the brief period before and after stance phase provides insights regarding how 235 the contact model behaves outside of stance phase and whether or not key gait events (e.g. heel strike and toe off) 236 occur at similar time points in the simulated and experimental data. Resultant ground reaction force time series were 237 re-sampled to 101 data points via cubic spline interpolation to allow for comparison between stance phases of 238 differing lengths. Ensemble curves (mean ± SD) were generated for each condition. The cumulative impulse of 239 GRFR was calculated to verify the dynamic compatibility of the simulation with the measured GRF. The time 240 integral of GRFR has been shown previously to be indicative of whole-body kinetics [32]. 241 Anterior-posterior center of pressure (CoPAP) position was calculated as the weighted sum of each contact 242 sphere’s predicted normal force multiplied by its anterior-posterior position (x) (eq. (3)). Raw normal forces arising 243 from each sphere during stance phase were low-pass filtered (4th order Butterworth: fc: 40 Hz) and summed. 244 Anterior-posterior CoP position was calculated throughout the stance phase. The CoPAP time series data were low- 245 246 CoP:;= ∑=!>"!#!$%∑>" (3) 247 CoPAP: Anterior-posterior center of pressure position 248 xi: Anterior posterior coordinate of contact sphere relative to the prosthesis heel (x = 0). 249 250 pass filtered (4th order Butterworth: fc: 6 Hz) and re-sampled to 101 data points via cubic spline interpolation to 251 allow for comparison between stance phases of differing lengths. Joint kinematics, GRFR, and CoPAP derived from 252 the simulations were compared to the corresponding data measured during experimental gait trials using coefficient 253 of determination and RMSE. 254 255 Manuscript body 6 3. Results 256 In the static condition, model-predicted subject mass (from ground reaction force) matched the measured mass 257 of each participant within 0.1%. In the dynamic conditions, simulated joint angles matched experimental joint angles 258 well (mean RMSE: 0.14 ± 0.01 deg, mean R2: 1.00 ± 0.00), indicating a successful forward-kinematics simulation. 259 All simulated joint angles matched experimental values with less than 0.23 deg RMSE. 260 Optimization of the medium stiffness trial for each subject resulted in a mean GRFR RMSE value of 7.7 ± 1.5% 261 body weight (BW) and a mean R2 of 0.97 ± 0.01 across stance phase (Figure 3). Cumulative impulse also matched 262 well (RMSE: 0.76 ± 0.26% BW·s, R2: 1.00 ± 0.00). In the time domain, model-predicted stance phases were, on 263 average, 0.03 ± 0.03 s shorter compared to experimental data for this optimized trial. 264 In the unoptimized trials, experimental GRFR and GRFR impulse responses were similar in the time and 265 amplitude domains across the three stiffness conditions (Figure 4). On average, stance phase times were 0.01 ± 0.01 266 s shorter in the simulations across the stiffness conditions. Time errors were least for the medium stiffness and 267 greatest for the high stiffness simulations. Variability for GRFR was greatest during the first 15% of stance phase for 268 simulations and lowest for the experimental data during this time. The contact parameters optimized for the medium 269 stiffness condition transferred well across the low and high stiffness conditions, which is evident by the similar 270 RMSE values for GRFR and GRFR impulse (Table 3). Model-predicted GRFR values were better in the high stiffness 271 configuration, whereas predictions for GRFR impulse were best in the medium and high stiffness conditions (Table 272 3). 273 Anterior-posterior CoP trajectory during stance phase was similar between simulated and experimental data 274 (Figure 5, Table 3). However, divergent trajectories were observed, primarily between 15–40% of stance phase. 275 Overall, RMSE values were 14 ± 4.4 percent foot length across all conditions, with the low stiffness condition 276 performing the best. Simulated data exhibited a strong correlation (R2 > 0.70) with experimental data for all stiffness 277 conditions. Simulation-predicted CoPAP showed similar variability compared to experimental data. 278 279 4. Discussion 280 The objective of this study was to develop and evaluate a forward simulation-capable gait model with below-281 knee amputation and a semi-active VSF prosthesis. The gait model and VSF-ground contact model were evaluated 282 under static and dynamic gait conditions. Under static conditions, model-predicted mass exhibited an error < 0.01% 283 with respect to the mass of the subjects. This value represents the performance of the contact model with no 284 dynamic component. Under dynamic gait conditions, the ability of the model to predict GRFR, GRFR impulse, and 285 COPAP within a forward kinematics framework was evaluated. For all trials, simulated joint angles were strongly 286 correlated with experimental angles (R2 > 0.999 for all joints) with RMSE values of less than 0.23 deg for all joints. 287 These values are logical for forward kinematics simulations, but importantly indicate that the model is numerically 288 stable when actuated by joint kinematics measured during gait with the VSF. Further, concomitant agreement 289 between the model’s kinematic and kinetic response is critical for forward dynamics simulations of human gait. 290 Optimization of contact model parameters in the medium stiffness trials achieved a mean RMSE value of 7.5 ± 291 1.5 percent BW for GRFR and 0.76 ± 0.26% BW·s for GRFR impulse. These values represent the theoretical optimal 292 performance of the comprehensive contact model within the present simulation framework. This contact model 293 performance is similar to the values previously reported in biomechanical contact modeling work [33–35]. However, 294 the focus of those studies was to predict foot-ground contact dynamics during gait for individuals with intact limbs. 295 Direct comparison of these data was limited to work in intact limb biomechanical modeling due to a lack of studies 296 reporting validation data for prosthesis-ground contact modeling in gait biomechanics. The strong correlation and 297 low error for GRFR impulse indicates that the contact model is able to predict the shape and trajectory of the GRFR 298 waveform arising from gait kinematics. Accurate prediction of GRFR impulse is important for capturing whole-body 299 energetics throughout gait [32]. The concomitant agreement for both simulated kinematics and kinetics further 300 suggests that these methods are viable for simulating whole-body movement during gait. 301 The transferability of the optimized contact model parameters from the medium stiffness condition was assessed 302 by simulating the two remaining medium stiffness trials and the six total trials with low and high stiffness 303 configurations per subject. Compared to the optimized trials, error for simulation-derived GRFR predictions 304 increased by an average of 4% across the unoptimized trials. Mean GRFR RMSE and R2 were 12.7 ± 6.9% BW and 305 0.91 ± 0.02 respectively for the remaining medium stiffness trials. These errors were slightly greater than those for 306 the low and high stiffness trials (Table 3). There was one outlier trial for the medium stiffness condition, classified 307 as a RMSE response greater than three standard deviations from the mean. This trial was not removed from the 308 mean values presented in Figures 3 and 5 and Table 3, since it could not be discerned what the source of the 309 increased error was. With this trial removed, the model’s GRFR response for the medium stiffness was 9.6 ± 2.6% 310 BW RMSE and 0.96 ± 0.01 R2. The same subject had a single high error trial in the low and high stiffness 311 Manuscript body 7 conditions. However, these trials did not meet the criteria for outlier classification and hence were included. Data 312 from these trials show a brief but large overshoot (>1 BW) for model-predicted GRFR at initial heel contact. Each of 313 these trials likely contributed to the higher mean error and variability in model-predicted GRFR early in stance phase 314 (Figure 4). Given that the outlier trials were all from the same subject, these errors may be due to errors or noise in 315 the kinematic data used to drive the model. Errors early in stance phase could be compounded by the relative 316 sensitivity of contact model parameters when a single sphere is in contact with the walking surface. 317 Model-predicted stance phase times agreed well with experimental values. On average, simulated stance phases 318 were 0.01 ± 0.01 s shorter than experimental times. Errors were similar between stiffness conditions. Accurate 319 prediction of stance phase length is important because it contributes to the model’s ability to quantify metrics such 320 as the time integral of GRFs, which may be indicative of whole-body kinetics [32]. Presented by increasing stiffness 321 condition, experimentally-measured stance phase times were 0.68 ± 0.06, 0.72 ± 0.06, and 0.71 ± 0.07 s. 322 Simulated CoPAP values were strongly correlated (R2 > 0.70) with experimental data in all stiffness conditions. 323 The RMSE values achieved using this model were similar to those reported in a previous study [36], which used a 324 subject-specific 38 contact point model to predict CoPAP for individuals with intact limbs. Accurate mapping of 325 CoPAP throughout stance phase is vital for simulating the effects of variable prosthesis stiffness on joint forces and 326 moments during gait. Errors in model-predicted CoPAP may be reduced by increasing the density of contact spheres 327 distributed on the plantar surface of the foot, which would improve the resolution of CoPAP predictions. However, 328 this would likely result in increased execution time for simulations due to increased model complexity and also 329 increase complexity of the contact parameter optimization problem. The low density of contact spheres in the heel of 330 the VSF model is a likely source of the steep inflection in the CoPAP trajectory early in stance phase, which contrasts 331 the more gradual progression depicted in the experimental data (Figure 5). 332 Computation times are an important consideration in simulation-based approaches. Previous work has shown 333 the lumped parameter approach for modeling ESR prosthesis dynamics to be more computationally efficient 334 compared to more detailed finite element models [20]. Total execution times for the gait model were 13.7 ± 2.48, 335 16.8 ± 4.66, and 64.3 ± 71.2 times slower than real time for the low, medium, and high stiffness conditions. 336 Increased execution times for the stiff conditions are a reflection of the need for small time-steps in solving a 337 rapidly-evolving, stiff differential equation. Computation times were also increased for participants with more mass. 338 These times could be reduced by utilizing “Accelerator” or “Rapid Accelerator” modes in Simscape Multibody, 339 which improve simulation times by generating a C-code executable of the model. For example, one subject’s 340 simulations computed 60% faster in Accelerator mode compared to Normal mode across the nine trials. Rapid 341 Accelerator mode was not tested in this study. All simulations were performed on a quad-core computer with a 4.0 342 GHz processor. 343 The present data show promise for predicting GRFR arising from gait with a semi-active VSF prosthesis. These 344 methods may be applied to the design and prescription of lower limb prostheses and forward dynamics simulations 345 in robotics and biomechanics. For example, simulations could be used to evaluate the potential effects of varied 346 prosthesis design parameters on the gait mechanics of a user. Simulations could also be computed within an optimal 347 control framework to identify optimal device configurations and manufacture customized prostheses. Evaluating 348 these effects within a simulation-based framework rather than traditional in vivo experimentation minimizes risk and 349 time spent by the user. Further, a broad spectrum of prosthesis design parameters could be modeled and simulated 350 without the need to manufacture multiple devices or the costs associated with doing so. Within biomechanics, 351 simulated stiffness could be varied to minimize a biomechanical cost function such as peak or average joint moment 352 or metabolic cost. Further optimization of the VSF-ground contact model may be necessary for simulation scenarios 353 with error tolerances less than 12% BW (Table 3). Similar improvements may be required if the mean difference 354 between simulation conditions is less than the error of the model, as was the case for the variation in GRFR by VSF 355 stiffness condition presented in this study. Reducing error in model-predicted GRFR may be accomplished through 356 improved methods in optimizing contact model parameters. For example, the objective function could be evaluated 357 under a variety of conditions, thereby improving the generalizability of the contact model. A deformable contact 358 model, such as presented in Jackson, Hass, and Fregly (2016), may also be a viable means of representing foam 359 deformation throughout stance phase and thus reducing error. 360 Previous work in gait simulations with biomechatronic devices has often relied on abstract representations of 361 components such as motors, electronics, and control systems [16,37,38]. In contrast, Simscape Multibody and 362 Simulink offer a large library of these components, which can be readily integrated into the model for more realistic 363 representations of biomechatronic systems. Further, control system parameters derived from simulations can be 364 readily deployed to prototype devices. As such, this model may be advantageous for simulating gait with 365 biomechatronic devices. 366 Manuscript body 8 These methods assume accurate estimation of segment length, joint centers, and joint angles which were 367 derived from marker-based motion capture data. Each of these metrics likely suffers from small errors due to marker 368 placement, localization, and coordinate system design. Such errors would contribute to decrements in contact model 369 performance. Characterizing the model’s sensitivity to varied joint kinematics is possible within the present 370 simulation paradigm. However, this analysis was not performed due to lack of experimental gait data necessary to 371 quantify accuracy of those simulations. 372 The components and joints of the prosthetic limb were also modeled as rigid, which may not be completely 373 accurate to represent the physical limb. This discrepancy would manifest as small differences in kinematics and 374 energy transfer between the components of the prosthetic limb. Nevertheless, simulated motions were consistent 375 with experimental data of subjects walking with the VSF and other previously reported data of spatiotemporal gait 376 patterns among persons with lower limb loss [39,40]. Another limitation is inherent to the reduced order design of 377 the lumped parameter VSF keel, which constrains keel motion to the sagittal plane. While this design is 378 computationally efficient compared to more detailed finite element models, it fails to account for small torsional 379 keel motions that would be possible under ecological gait conditions with the physical VSF. Future work could 380 develop a lumped parameter model that allows for torsional keel movement in the frontal plane. Within the gait 381 model, future work should incorporate biologically-inspired muscle models (e.g. Hill [41] or Thelen [42] models) to 382 actuate joints for forward dynamics simulations. This would improve the biomechanical and physiological validity 383 of the model and also allow for neuromuscular evaluation of simulated gait. 384 385 5. Conclusions 386 The present study demonstrates that the ESR properties of a semi-active VSF can be modeled and integrated 387 into a scalable gait model that accounts for altered lower limb mass and inertial properties associated with LLA. The 388 model captured whole-body kinetics associated with gait with varied prosthesis stiffnesses. Foot-ground contact 389 models were used to estimate GRFR with 7.7% BW mean RMSE in optimized gait trials, which translated to a mean 390 RMSE of 12% across unoptimized trials. The contact models also predicted COPAP with RMS error of 14% foot 391 length. This model performance may be sufficient for gait simulations among persons with lower limb loss. Such 392 simulations may be used to aid in the prosthesis design and prescription process in order to improve user mobility. 393 These methods may also be helpful to identify other important prosthesis design parameters, which can be modified 394 to optimize gait. Further contact model optimization and error reduction may be required for simulation-based 395 comparisons of varied prosthesis stiffness, where differences in GRFR magnitude may be nuanced. 396 Statement of acknowledgements Statement of acknowledgements: The authors would like to thank Evan Glanzer for his work in developing and 397 testing the variable-stiffness foot and Terry Denery for his intellectual contributions in optimizing the contact model. 398 399 Funding sources: This work was funded in part by the Lokey Doctoral Science Fellowship (MAM) and NIH grant 400 HD074424 (PGA).401 402 Nomenclature 403 Acronyms widely used in text BKA Below-Knee Amputation BW Body Weight; M*g CoP Center of Pressure DoF Degrees of Freedom ESR Energy Storage and Return GRFR Resultant Ground Reaction Force, N; =>?!?@+ >?!A@+ >?!B@ LHS Latin Hypercube Sampling LLA Lower Limb Amputation ode15s Ordinary differential equation 15 solver SD Standard Deviation RMSE Root Mean Square Error; =∑(D=(,+&E,!$%-!F H&EI-%$&)!!)&#!$%K VSF Variable Stiffness Foot Abbreviations a Overhung length, mm b Damping coefficient, N×s/mm B Support fulcrum position, mm D Displacement, mm F Force, N k Linear stiffness, N/mm L Total beam length, mm l Supported length, mm n Penetration exponent R2 Coefficient of determination µ Coefficient of friction v Linear velocity, mm/s y Scaling factor ( Penetration depth, mm (̇ Penetration velocity, mm/s w Angular velocity, rad/s Superscripts and subscripts CoPAP Anterior-posterior (Center of Pressure) Ff Frictional force, N Fn Normal force, N GRFR Resultant ground reaction force, N vpoc Linear velocity at point of contact, mm/s vthreshold Linear velocity threshold, mm/s µkinetic Coefficient of kinetic friction µstatic Coefficient of static friction 404