Introduction: Passive energy storage and return (ESR) feet are the current performance standard in lower limb 24 prostheses. A recently developed semi-active variable-stiffness foot (VSF) prosthesis balances the simplicity of a 25 passive ESR device with the adaptability of a powered design. The purpose of this study was to model and simulate 26 the ESR properties of the VSF prosthesis. Methods: The ESR properties of the VSF were modeled as a lumped 27 parameter overhung beam. The overhung length is variable, allowing the model to exhibit variable ESR stiffness. 28 Foot-ground contact was modeled using sphere-to-plane contact models. Contact parameters were optimized to 29 represent the geometry and dynamics of the VSF and its foam base. Static compression tests and gait were 30 simulated. Simulation outcomes were compared to corresponding experimental data. Results: Stiffness of the model 31 matched that of the physical VSF (R2: 0.98, RMSE: 1.37 N/mm). Model-predicted resultant ground reaction force 32 (GRFR) matched well under optimized parameter conditions (R2: 0.98, RMSE: 5.3% body weight,) and unoptimized 33 parameter conditions (R2: 0.90, mean RMSE: 13% body weight). Anterior-posterior center of pressure matched well 34 with R2 > 0.94 and RMSE < 9.5% foot length in all conditions. Conclusions: The ESR properties of the VSF were 35 accurately simulated under benchtop testing and dynamic gait conditions. These methods may be useful for 36 predicting GRFR arising from gait with novel prostheses. Such data are useful to optimize prosthesis design 37 parameters on a user-specific basis.38 Manuscript body 2 1. Introduction 39 Individuals with lower limb loss exhibit distinct gait characteristics, which may limit mobility and decrease 40 quality of life. Those using lower limb prostheses may display gait asymmetry [1,2], elevated metabolic cost during 41 locomotion [3], and a variety of psychological disorders including anxiety and depression [4]. Sustained prosthesis 42 use may also induce overloading of intact joints and ultimately, musculoskeletal ailments [5]. Each of these issues 43 may be attenuated by improving user specificity in the design characteristics of foot prostheses. However, the effects 44 of foot prosthesis design parameters (e.g. stiffness) are not well characterized, and thus achieving meaningful 45 improvements in gait has proven arduous [6,7]. In order to achieve improvements, a robust understanding of the 46 relationships between anthropometry, gait mechanics, and prosthesis design are necessary. 47 One of the primary design goals of a lower limb prosthesis is to replace the coordinated energy absorption and 48 generation properties of a lost limb. Passive energy storage and return (ESR) foot prostheses are the current standard 49 for mimicking this functionality. However, the fixed stiffness behavior of these devices contrasts that of the healthy 50 foot-ankle complex, which modulates its behavior in response to varied gait conditions (e.g. velocity and terrain) 51 [8,9]. Glanzer and Adamczyk (2018) [10] recently developed a variable-stiffness foot (VSF) prosthesis designed 52 with an actuated keel support fulcrum to semi-actively control sagittal forefoot stiffness and thereby adapt to 53 different gait conditions with low power (Fig. 1). The ESR keel of the VSF is a composite leaf spring designed as an 54 overhung beam, which modulates the supported length (l) via an actuated keel support fulcrum (B). The total beam 55 length (L) is 229 mm, whereas the overhung length (a) is variable between 66–151 mm. By modulating overhung 56 length, the VSF’s forefoot is capable of exhibiting roughly a three-fold range of forefoot stiffness values (10–32 57 N/mm). The heel component of the VSF has a consistent linear stiffness of 65 N/mm. The VSF’s fulcrum position is 58 designed to be adjusted during swing phase, thus minimizing the power necessary for actuation. As such, the VSF 59 behaves primarily as a passive ESR prosthesis, which can adapt stiffness in response to variable gait conditions. 60 Simulations based on computational models can be powerful tools for evaluating potential biomechanical 61 interventions, such as the implementation of a novel ESR prosthesis. Recently, simulations have been used to aid in 62 the iterative design process and improve user-specificity [11–13]. Inverse simulations provide the ability to estimate 63 values that cannot be measured in vivo (e.g. socket-residual limb interface dynamics), whereas predictive 64 simulations suggest hypotheses regarding how humans may interact with and adapt to new prosthetic devices. 65 Computational modeling has been used to investigate the effects of prosthesis alignment [14] and a biarticular 66 clutched spring mechanism [15] on gait mechanics among persons with lower limb loss. However, these models do 67 not account for the ESR properties of the prosthetic foot, thus limiting their ecological validity. Other studies, which 68 did incorporate the force and torque contributions of ESR feet into gait models focused on characterizing 69 biomechanical and myophysiologic responses with prosthesis use, rather than validation of the prosthesis model 70 [16,17]. While these studies made important progress toward investigating the relationship between anthropometry, 71 gait mechanics, and prosthetic foot design, they had limited ability to verify simulation results in the context of 72 experimental values. Due to these limitations, the use of simulations to inform the design of ESR foot prostheses has 73 not been fully realized. The purpose of this study was to further couple experimental and simulation prosthesis data 74 by modeling and validating the mechanical stiffness properties and resulting ground reaction forces of a semi-active 75 VSF. 76 77 2. Methods 78 79 2.1 Model design 80 A computational model of the VSF was developed in Simscape Multibody (Mathworks, Inc., Natick, MA). 81 The assembly, geometry, mass, and inertial properties were derived from SolidWorks (Dassault Systemes Inc., 82 Waltham, MA). A reduced order model of the VSF’s variable-stiffness elastic keel was designed using the lumped 83 parameter approach for approximating flexible body dynamics. This approach involved discretizing the continuous 84 geometry of the keel into finite rigid segments coupled via revolute joints, springs, and dampers (Fig. 2). This 85 simplification of the original state space of the continuous elastic keel system to finite dimensions allows the partial 86 differential equations of the infinite-dimensional time-space states of the physical VSF to be represented by ordinary 87 differential equations with a finite number of parameters. 88 The keel of the VSF model was discretized into 16 segments (eight DoF). The most posterior segment is 66 89 mm in length, which matches the minimum possible fulcrum position. The rest of the keel consists of 11.64-mm 90 segments for a total beam length of 229 mm (Fig. 2). The stiffness and damping values for the revolute joints were 91 parameterized to represent the material properties of the VSF’s G10/FR4 Garolite keel (flexural elastic modulus: 92 Manuscript body 3 18.6 GPa, Poisson’s ratio: 0.136). A MATLAB script controls continuous fulcrum position (i.e. variable stiffness). 93 The VSF model was rigidly attached to a prosthetic pylon and socket via a pyramid adapter, as the device would be 94 used in vivo. These connections were modeled as weld joints. Each segment is independently scalable, allowing the 95 model to be integrated into an anatomically scaled computational gait model. 96 Foot-ground contact consists of 24 sphere-to-plane contact models [18] parameterized to represent the 97 geometry and dynamics of the VSF’s foam base. Each of these models estimates normal (Fn) and frictional (Ff) 98 forces associated with the collision of a viscoelastic sphere (a massless spring and damper system) and a rigid plane 99 (Fig. 3). The overall foot contact model was divided into five zones; the sphere-to-plane models were parameterized 100 by zone (Fig. 3, Table 1). The heel of the VSF model is comprised of three zones; this choice was motivated by the 101 sensitivity of contact parameters when few spheres are in contact with the walking plane (e.g. the heel of the foot 102 early in stance phase). Contact parameters are less sensitive when many spheres are in contact with the walking 103 plane (e.g. the midfoot and forefoot late in stance phase). The foam base of the physical VSF undergoes 104 compression throughout stance phase. To account for these effects, a modified Kelvin-Voigt nonlinear spring and 105 damper force law (eq. (1)) was implemented to represent contact between the VSF and walking plane: 106 107 !!= $(&×(!)++(,×()̇(>0,(̇>0&×((>0,(̇<00(<0 (1) 108 Fn: normal force 109 k: contact stiffness 110 2: penetration depth 111 n: penetration exponent 112 y: damping force scaling factor 113 3: contact damping coefficient 114 115 The spring force increases exponentially as the sphere penetrates the contact plane. The damping force is multiplied 116 by a scaling factor (y), which increases from zero to one as a polynomial as it approaches a user-defined value for 117 full damping. Frictional force (eq. (2)) is the product of the normal force and coefficient of friction (µ). A stick-slip 118 friction law defines the transition between static (µstatic) and kinetic (µkinetic) coefficients of friction based on a 119 velocity threshold (vthresh): 120 121 !"= 4!! × 5#$%$&'6()'<6$*+,#*)-.!! × 5/&!,$&'6()'>6$*+,#*)-. (2) 122 Ff: frictional force 123 7: coefficient of friction 124 8012: velocity at point of contact 125 8345674189: velocity threshold 126 127 Static and kinetic coefficients of friction were set to 0.5 and 0.3 with a velocity threshold of 0.1 m/s. Resultant 128 ground reaction force (GRFR) was derived by summing and low-pass filtering (4th order Butterworth, ƒc: 40 Hz) the 129 normal and frictional forces arising from each contact sphere. 130 In order to improve GRFR predictions, contact model parameterization was formulated as a least-squares 131 optimization problem with the objective of minimizing the sum of squared errors between model-predicted and 132 experimentally measured GRFR (see “Model Validation”). Initial parameter settings at the outset of the optimization 133 were derived by increasing stiffness until the contact spheres were able to support the weight of the model. Initial 134 damping coefficients (N×s/mm) were set to half the numerical value of stiffness (N/mm). Penetration exponents and 135 penetration for full damping values were initialized at 1 and 1 mm, respectively. These initial values were used as 136 inputs to the problem. Latin hypercube sampling (LHS) was applied to generate simulation scenarios with pseudo-137 random sets of parameters. The LHS approach is a method of stratified sampling, which divides parameter values 138 into equal strata based on an assumed normal distribution and constrained by user-defined bounds. Random 139 parameter values are sampled from within these strata to generate a simulation scenario with a pseudo-random set of 140 parameters. The LHS technique effectively samples the search space, while providing the randomness required to 141 explore the efficacy of a range of variable values to minimize the objective function. The objective function value of 142 each iteration is compared to the previous iteration; the parameter scenario which best minimizes the objective is 143 passed to the next iteration of the algorithm. The optimization algorithm proceeds for 100 iterations or until an 144 Manuscript body 4 objective function tolerance of 0.1 N is reached (i.e. convergence). If the optimization algorithm did not meet any of 145 the termination criteria, the initial parameter values were updated using the results of the first run, and an additional 146 run was initiated. Parameter tolerances were set to 0.001 (varying units) in order to avoid false minima. 147 148 2.2 Model validation 149 150 2.2.1 Static compression testing 151 The operational stiffness range of the physical VSF was determined through static compression testing 152 (TestResources, Shakopee, MN) (Glanzer and Adamczyk 2018). Load was applied at a constant speed of 50 153 mm/min to a point 30 mm proximal to the anterior tip of the VSF (i.e. supported beam length = 199 mm). To 154 validate the ESR properties of the VSF model, a simulated materials testing system (MTS) was developed in 155 Simscape Multibody. The MTS simulator consists of a massless body, which translates vertically according to a 156 user-defined time-position vector (Fig. 2). Simulated static compression tests were performed as in Glanzer and 157 Adamczyk (2018). Contact was maintained throughout VSF deflection. Contact dynamics between the VSF and 158 MTS were estimated using a sphere-to-sphere contact model. Stiffness (k) (eq. (3)) was computed as the average 159 slope of the load-displacement data for loads above 200 N. 160 161 &= ∆;!"#$∆;%&'()#*+,+-. for loads 200 N to !<%= (3) 162 163 Deformation for loads under 200 N was considered to arise primarily from foam compression, rather than keel 164 displacement. Mid-range keel displacement was also calculated for the VSF model as the displacement of the keel at 165 50 percent of the maximal load applied during the static compression test. 166 Static compression tests were simulated at five discrete fulcrum positions (66, 87, 108, 129, and 151 mm), 167 which span the full continuous range of possible positions. These ascending fulcrum positions represent decreases in 168 overhung length (a) depicted in Figure 1, and therefore yield increases in endpoint stiffness. Simulation-derived 169 values were compared to those from static compression tests of the physical VSF via coefficient of determination 170 and root mean squared error (RMSE). Simulations were calculated in Simscape Multibody using the ode15s solver 171 profile with variable step size. 172 173 2.2.2 Gait conditions 174 Model-predicted GRFR was validated under two scenarios: static and dynamic gait conditions. For both 175 validations, the VSF model was integrated into a seven-segment, 28-DoF anatomically-scaled gait model of a 176 subject with a unilateral transtibial amputation. Three-dimensional optical motion capture data (Optitrack, Natural 177 Point, Inc. Corvallis, OR) of a male subject (181 cm, 78.0 kg) with a right side transtibial amputation walking with 178 the physical VSF were used as inputs to the model. Retroreflective marker coordinates from a static motion capture 179 trial were used to estimate and scale limb dimensions for the pelvis, leg, intact shank, residual shank, and intact foot. 180 Within the gait model, the residual shank was encapsulated in a prosthetic socket and welded to the pyramid adapter 181 of the VSF model (Fig. 2). The interface between the prosthetic socket and residual limb was modeled as a high-182 stiffness 6-DoF bushing joint, similar to previous work by LaPrè et al. (2018). The rotational and translational 183 stiffness as well as displacement and velocity constraints were designed according to previous gait experiments [19] 184 and finite element analysis [20]. The mass and inertial properties of the lower limbs and pelvis were modeled as 185 conical frusta and an ellipsoid, respectively. Segment masses were estimated according to De Leva (1996). 186 For the static condition, the model was simulated with anatomically neutral joint angles for ten seconds. Model-187 predicted GRFR was averaged over the course of the trial and compared to the mass of the subject. Dynamic gait 188 simulations were calculated based on experimental motion capture trials of the subject walking over ground between 189 1.0 and 1.2 m/s with the VSF under low, medium and high stiffness configurations (fulcrum positions: 66, 108, and 190 151 mm). Three trials were collected for each stiffness configuration for a total of nine trials. Three-axis pelvis, hip, 191 knee, and ankle angles were calculated from three-dimensional marker coordinate data [22,23] and used as inputs to 192 drive the corresponding joints of the model. Motion at the socket-limb interface was considered to be passive based 193 on the aforementioned velocity and displacement constraints. The pyramid adapter-pylon interface was assumed to 194 be rigid. 195 Contact model-derived GRFR prediction was optimized for a single trial at the 66-mm fulcrum position. 196 The GRFR error resulting from this trial represents the theoretical optimal performance of the comprehensive VSF-197 ground contact model. The transferability of the optimized parameter values was determined by simulating the two 198 remaining low stiffness trials and the three remaining trials each for the medium and high stiffness configurations. 199 Manuscript body 5 Joint kinematics and GRFR data were low-pass filtered (4th order Butterworth: fc: 6 Hz and 40 Hz, 200 respectively). Simulation and experimental GRFR were time locked and indexed to 0.25 s before and 0.25 s after 201 stance phase. Including the brief period before and after stance phase provides insights regarding how the contact 202 model behaves outside of stance phase and whether or not key gait events (e.g. heel strike and toe off) occur at 203 similar time points in the simulated and experimental data. Resultant ground reaction force time series were re-204 sampled to 101 data points via cubic spline interpolation to allow for comparison between stance phases of differing 205 lengths. Ensemble curves (mean ± SD) were generated for each condition. The impulse of GRFR was calculated to 206 assess the simulation’s ability to predict GRFR trajectory. 207 Anterior-posterior center of pressure (CoPAP) position was calculated as the weighted sum of each contact 208 sphere’s predicted force multiplied by its anterior-posterior position (x). Raw normal forces arising from each sphere 209 during stance phase were low-pass filtered (4th order Butterworth: fc: 40 Hz) and summed. Anterior-posterior CoP 210 position was calculated across stance phase (eq. (4)). 211 212 CoP>?= ∑=/A0/1/23∑A0 (4) 213 CoPAP: Anterior-posterior center of pressure position 214 xi: Anterior posterior coordinate of contact sphere 215 216 The CoPAP time series data were low-pass filtered (4th order Butterworth: fc: 6 Hz) and re-sampled to 101 data points 217 via cubic spline interpolation to allow for comparison between stance phases of differing lengths. Joint kinematics, 218 GRFR, and CoPAP data measured during experimental gait trials were compared to those derived from the 219 simulations using coefficient of determination and RMSE. 220 221 3. Results 222 223 3.1 Static compression tests 224 Simulated VSF stiffness effectively reproduced experimental stiffness across the five fulcrum 225 configurations (R2 > 0.98, RMSE = 1.37 N/mm) (Fig. 4, Table 2). Simulated mid-range displacement also matched 226 well (R2 > 0.99) with small offset from experimental displacement in each condition (RMSE = 0.45 mm). 227 Experimental load-displacement relationships were most linear in the 66 and 87 mm fulcrum configurations, as 228 indicated by variance in the slope of the relationship. The stiffest three conditions exhibited curvilinear relationships. 229 Simulated load-displacement data were linear in all conditions due to the linear spring and damper force parameters 230 for the revolute joints in the lumped parameter keel model. 231 232 3.2 Resultant ground reaction force predictions 233 In the static condition, model-predicted subject mass was 2.6 ± 0.0% less than measured mass. In the 234 dynamic conditions, simulated joint angles matched experimental joint angles well, but exhibited a small phase lag 235 (mean RMSE: 1.9 ± 1.0 deg, mean R2: 0.98 ± 0.02). Simulated and experimental GRFR data agreed well in the time 236 domain (Fig. 6). Amplitude discrepancies, quantified via RMSE, were least in low stiffness configuration and 237 greatest in the high stiffness configuration. Coefficient of determination values were similar for the low and medium 238 stiffness conditions and lower for the high stiffness condition. Impulse was similar in the low and high stiffness 239 conditions and lower for the medium stiffness condition (Table 3). 240 Optimization of the single low stiffness trial resulted in a GRFR RMSE of 5.3% body weight (BW) and R2 241 of 0.98 across stance phase. Impulse also matched well (RMSE: 0.01 BW×s, R2 > 0.99) (Fig. 5). In the time domain, 242 model-predicted heel contact preceded experimental heel contact by 0.02 s, resulting in a 0.02-s longer stance phase. 243 Simulating the two additional low stiffness trials with the optimized contact parameters resulted in average RMSE 244 and R2 values of 0.10 ± 0.05 BW and 0.93 ± 0.05 for GRFR and 0.02 ± 0.01 BW×s and > 0.99 ± 0.01 for GRFR 245 impulse (Fig. 6, Table 3). 246 Experimental GRFR and GRFR impulse responses were similar in the time and amplitude domains across 247 the three stiffness conditions (Fig. 6). On average, stance phase time was 0.05 ± 0.03 s longer in the simulations 248 across the stiffness conditions. Time errors were least in the low stiffness condition and greatest in the high stiffness. 249 Variability for GRFR was greatest during the first 25% of stance phase for all conditions. Variability for GRFR 250 impulse was greatest near the end of stance phase. The ability of the contact parameters optimized for the low 251 stiffness condition transferred well across the other two conditions, which is evident by the similar RMSE values for 252 GRFR (Table 3). Resultant ground reaction force RMSE and R2 values were better in the medium stiffness 253 configuration, whereas RMSE and R2 were better in the high stiffness condition for GRFR impulse. The medium 254 Manuscript body 6 stiffness condition demonstrated the least variability for the GRFR response, whereas the low and high stiffness 255 conditions showed similarly low variability for GRFR impulse (Table 3). 256 Anterior-posterior CoP trajectory during stance phase was similar between simulated and experimental data 257 (Fig. 7). Root mean squared errors were 8.9 ± 1.0, 9.5 ± 0.9, and 5.7 ± 1.4 percent foot length for the low, medium, 258 and high stiffness conditions, respectively (Table 3). Simulated data correlated well with experimental data across 259 all conditions. Coefficient of determination values were 0.95 ± 0.01, 0.94 ± 0.01, and 0.97 ± 0.01 for the low, 260 medium, and high stiffness conditions. 261 262 4. Discussion 263 264 The goal of this study was to develop a reduced-order computational model of a semi-active variable-265 stiffness foot prosthesis. Results from simulated static compression tests showed good agreement with experimental 266 data. These outcomes suggest that the variable-stiffness ESR properties of the VSF were modeled with high fidelity 267 using a reduced order lumped parameter approach for approximating flexible body dynamics. One of the goals of 268 reduced order modeling is to capture a structure’s dynamic behavior in a computationally inexpensive way. A 269 common benchmark for reduced-order models is the ability to simulate at or near real-time [24,25], which contrasts 270 with more computationally expensive methods such as mesh-based finite element modeling. Including initialization 271 time, static compression simulations computed 3.3 ± 0.8 times faster than real-time (i.e. the length of time required 272 to complete the experimental static compression test) on computer with a four core 4.0 GHz processor. Initialization 273 time, which includes model compiling and building, can be minimized using “Accelerator” and “Fast Restart” 274 modes in Simscape Multibody. Using these tools, simulations computed 39 ± 16 times faster than real time. This 275 computational efficiency is useful if the model is to be simulated iteratively, for example in parameter optimization 276 or machine learning frameworks. 277 The VSF is an ideal device for studying the effects of prosthesis stiffness on gait mechanics because it can 278 readily exhibit a range of forefoot stiffness values, thereby eliminating the need to purchase or manufacture multiple 279 prostheses as in [26–28]. In doing so, this also eliminates confounding variables that accompany a foot-switching 280 experimental design, such as mismatched or out-of-order stiffness from foot prostheses of different categories. The 281 VSF can also modulate stiffness along a continuous scale, which provides improved resolution compared to the 282 typical discrete stiffness options available for fully passive designs. The range of forefoot stiffness values exhibited 283 by the physical VSF and captured by the VSF model represent a range of stiffness values available in many 284 commercially-available prosthetic feet [29,30]. Accurate characterization of this range is important, should this 285 model be used to inform the design and/or prescription of prosthetic feet. Further, this model can be easily re-286 parameterized to exhibit a different range of stiffness values, which could aid in the selection of keel dimensions or 287 material properties to meet design goals. Two primary limitations are present for the static compression testing 288 simulations. Experimental load-displacement data were only available for positive loading conditions, and thus a 289 comparison of the model’s hysteresis behavior was not possible. Similarly, experimental data were only available 290 for the 50 mm/min loading rate. A robust characterization of the VSF’s stiffness behavior under a range of higher 291 loading rates would likely improve the model’s behavior under dynamic conditions. Experimental load-displacement 292 data could also be influenced by imperfections in maintaining a constant contact point with the prosthesis. Results of 293 the present study are difficult to compare to previous work, as there is a paucity of previous research that evaluated 294 simulated prosthesis dynamics compared to mechanical testing data of a physical prosthesis under multiple 295 conditions. However, errors exhibited by this model are similar to those reported in Tryggvason et al. (2020), who 296 compared the angular stiffness response of a finite element foot prosthesis model to data from mechanical tests [13]. 297 Under dynamic gait conditions, simulated joint angles agreed well with experimental values, indicating that 298 the model is numerically stable when actuated by joint kinematics measured during gait with the VSF. Joint angles 299 were strongly correlated, but exhibited a small phase lag, possibly due to ODE solver settings and numerical 300 integration. This phase lag may be also be present in the kinetic data, but masked by the larger inherent variability of 301 the simulated GRFR. Total simulation times were 8.95 ± 3.92, 12.7 ± 0.67, and 46.2 ± 1.19 times slower than real 302 time for the low, medium, and high stiffness configurations respectively. Execution times were 3.12 ± 0.10, 3.40 ± 303 0.67, and 38.4 ± 1.19 times slower than real time. Increased execution times for the stiff conditions may reflect the 304 need for small time-steps in solving a rapidly-evolving, stiff differential equation. 305 Optimization of the GRFR for the low stiffness configuration achieved a RMSE of 5.3% BW and R2 of 306 0.98. These values are similar to those reported in previous biomechanical contact modeling work [31–33]. 307 However, those studies focused on quantification of foot-ground contact during gait for individuals with intact 308 limbs. Direct comparison of these data was limited to work in intact limb biomechanical modeling due to a lack of 309 studies reporting validation data for prosthesis-ground contact modeling in gait biomechanics. The strong correlation 310 Manuscript body 7 and low error for GRFR impulse indicates that the contact model is able to predict the shape and trajectory of the 311 GRFR arising from gait kinematics. Accurate predictions of GRFR impulse is important for capturing whole-body 312 energetics throughout gait. The concomitant agreement for both kinematics and kinetics further suggests that these 313 methods are viable for simulating whole-body energetics during gait. 314 The transferability of the optimized contact model parameters from the low stiffness condition was 315 assessed by simulating two additional low stiffness trials and three trials each with medium and high stiffness 316 configurations. Compared to the optimized trial, simulation-derived GRFR predictions did not perform as well in the 317 unoptimized trials. Mean GRFR RMSE and R2 were 12.7 ± 1.44% BW and 0.91 ± 0.02 for the remaining low 318 stiffness trials. These values were similar for the medium and low stiffness trials (Table 3). The impulse of these 319 data matched well across the unoptimized trials (RMSE: 0.03 ± 0.02 BW×s, R2: 0.98 ± 0.01). Variability of the 320 model’s performance was similar across the unoptimized conditions for all outcome measures. It is possible that the 321 contact model parameters were over-fitted to the specific conditions of a single trial, resulting in decreased 322 generalizability. Future work should assess the balance between optimization specificity and generalizability. 323 The amplitude and shape of experimental GRFR waveforms were similar across the three stiffness 324 conditions. However, stance phase times did vary by condition for the subject tested. The medium stiffness 325 condition resulted in the longest stance phase time (0.79 ± 0.01 s), high stiffness resulted in the shortest (0.71 ± 0.02 326 s), and low stiffness (0.73 ± 0.02 s) was in the middle. The same pattern was present in the simulated data, although 327 stance phase times were 0.05 ± 0.03 s longer on average compared to the experimental data. Stance phase times 328 derived from simulations were correlated with experimental times (R2 = 0.65). More data are necessary to discern 329 the strength, repeatability, and significance of these relationships. 330 Simulated CoPAP values agreed well with experimental values. The RMSE values achieved using this 331 model were similar to those reported in previous work involving subject-specific biomechanical contact modeling 332 for individuals with intact limbs [34]. Accurate mapping of CoPAP throughout stance phase is vital for simulating the 333 effects of variable prosthesis stiffness on joint forces and moments during gait. Errors in model-predicted CoPAP 334 may be reduced by increasing the density of contact spheres distributed on the plantar surface of the foot, which 335 would improve the resolution of CoPAP predictions. However, this would likely result in increased execution time 336 for simulations and also increase complexity of the contact parameter optimization problem. 337 The present data show promise for predicting GRFR arising from a semi-active VSF prosthesis. These 338 methods may be applied to the design and prescription of lower limb prostheses and forward dynamics simulations 339 in robotics and biomechanics. Within biomechanics, future work could integrate the VSF model into a gait model of 340 an individual with lower limb loss. Gait simulations could be formulated as an optimal control problem in which 341 prosthesis stiffness is tuned to minimize a biomechanical cost function such as joint loading or metabolic cost. 342 Evaluating these effects within a simulation-based framework rather than traditional in vivo experimentation 343 minimizes risk and time spent by the user. Further, a broad spectrum of prosthesis design parameters could be 344 modeled and simulated without the need to manufacture multiple devices or the costs associated with doing so. 345 Further optimization of the VSF-ground contact model may be necessary for simulation scenarios with error 346 tolerances less than 12% BW. Similar improvements may be required if the mean difference between simulation 347 conditions is less than the error of the model. Reducing error in model-predicted GRFR may be accomplished by 348 evaluating the objective function under a variety of conditions and choosing the parameter set that achieves the best 349 minimization across several conditions. A deformable contact model, such as presented in Jackson, Hass, and Fregly 350 (2016), may also be a viable means of representing foam deformation throughout stance phase and thus reducing 351 error. 352 These methods assume accurate estimation of segment length, joint centers, and joint angles which were 353 derived from marker-based motion capture data. Each of these metrics likely suffers from small errors due to marker 354 placement, localization, and coordinate system design. Such errors would contribute to decrements in contact model 355 performance. The components and joints of the prosthetic limb were also modeled as rigid, which may not be 356 completely accurate to represent the physical limb. This discrepancy would manifest as small differences in 357 kinematics and energy transfer between the components of the prosthetic limb. Nevertheless, simulated motions 358 were consistent with experimental data of subjects walking with the VSF and other previously reported data of 359 spatiotemporal gait patterns among persons with lower limb loss [35,36]. Another limitation is inherent to the 360 reduced order design of the lumped parameter VSF keel, which constrains keel motion to the sagittal plane. While 361 this design is computationally efficient compared to more robust finite element models, it fails to account for small 362 torsional keel motions that would be possible under ecological gait conditions with the physical VSF. 363 364 5. Conclusions 365 366 Manuscript body 8 The present study demonstrates that the ESR properties of a semi-active VSF can be modeled with high 367 fidelity. Foot-ground contact models were used to estimate GRFR with 5.3% BW error in an optimized gait trial, 368 which translated to mean errors of 13% for unoptimized trials. The contact models also predicted COPAP with mean 369 error of 9.3% foot length. This model performance may be sufficient for gait simulations among persons with lower 370 limb loss. Such simulations may be used to aid in the prosthesis design and prescription process in order to improve 371 user mobility. These methods may also be helpful to identify other important prosthesis design parameters, which 372 can be modified to optimize gait. Further contact model optimization and error reduction may be required for 373 simulation-based comparisons of varied prosthesis stiffness, where differences in GRFR magnitude may be nuanced. 374 Statement of acknowledgements Statement of acknowledgements: The authors would like to thank Evan Glanzer for his work in developing and 375 testing the variable-stiffness foot and Terry Denery for his intellectual contributions in optimizing the contact model. 376 377 Funding sources: This work was funded in part by the Lokey Doctoral Science Fellowship (MAM) and NIH grant 378 HD074424 (PGA).379 380 Nomenclature 381 Acronyms widely used in text BW Body weight; M*g CoP Center of pressure DoF Degrees of Freedom ESR Energy storage and return GRFR Resultant Ground Reaction Force, N; =>?!BC+ >?!DC+ >?!EC LHS Latin Hypercube Sampling MTS Material Testing System ode15s Ordinary differential equation 15 solver SD Standard Deviation RMSE Root Mean Square Error; =∑(G=(,+&<,!$%-/H J&<K-%$&)!/)41/23M 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 y Scaling factor ( Penetration depth, mm (̇ Penetration velocity, mm/s w Angular velocity, rad/s Superscripts and subscripts CoPAP Anterior-posterior (Center of Pressure) Dsim Simulation (Displacement) Dexp Experimental (Displacement) Ff Frictional force, N Fn Normal force, N GRFR Resultant ground reaction force, N ksim Simulation (stiffness), N/mm kexp Experimental (stiffness), N/mm vpoc Linear velocity at point of contact, mm/s vthreshold Linear velocity threshold, m/s µkinetic Coefficient of kinetic friction µstatic Coefficient of static friction 382