and key terms Introduction: Loading of a residual limb within a prosthetic socket can cause tissue 50 damage such as ulceration. Computational simulations may be useful tools for 51 estimating tissue loading within the socket, and thus provide insights into how 52 prosthesis designs affect residual limb-socket interface dynamics. The purpose of this 53 study was to model and simulate residual limb-socket interface dynamics and evaluate 54 the effects of varied prosthesis stiffness on interface dynamics during gait. Methods: A 55 spatial contact model of a residual limb-socket interface was developed and integrated 56 into a gait model with a below-knee amputation. Gait trials were simulated for four 57 subjects walking with low, medium, and high prosthesis stiffness settings. The effects of 58 prosthesis stiffness on interface kinematics, normal pressure, and shear stresses were 59 evaluated. Results: Model-predicted values were similar to those reported previously in 60 sensor-based experiments; increased stiffness resulted in greater average normal 61 pressure and shear stress (p<0.05). Conclusions: These methods may be useful to aid 62 experimental studies by providing insights into the effects of varied prosthesis design 63 parameters or gait conditions on residual limb-socket interface dynamics. The current 64 results suggest that these effects may be subject-specific. 65 66 Key terms: Simulation, Computational model, prosthetic socket, residual limb, 67 rehabilitation68 Glossary of terms Nomenclature 69 Acronyms widely used in text AP Anterior-Posterior BKA Below-Knee Amputation DoF Degrees of Freedom LLA Lower Limb Amputation ML Medial-Lateral Ode23t Ordinary Differential Equation 23 trapezoidal solver PD Proximal-Distal SD Standard Deviation VSF Variable Stiffness Foot Abbreviations F Force, N k Linear stiffness, N/mm µ Coefficient of friction v Linear velocity m/s ! Penetration depth, m !̇ Penetration velocity, m/s Superscripts and subscripts Ff Frictional force, N Fn Normal force, N vpoc Linear velocity at point of contact, m/s vthreshold Linear velocity threshold, m/s µkinetic Coefficient of kinetic friction µstatic Coefficient of static friction 70 Manuscript Body 1.0 Introduction 71 Rehabilitation following a lower limb amputation (LLA) often includes prescription of a 72 prosthesis designed to replace the functionality of the removed limb. For an individual with a 73 below-knee amputation (BKA), a prosthesis system typically consists of a socket, which 74 interfaces with the residual limb, a rigid pylon, and a foot-ankle prosthesis. Use of lower limb 75 prostheses can improve mobility, health, and quality of life. However, abnormal loading of the 76 soft tissues surrounding the truncated shank (e.g., asymmetric pressure distribution and 77 elevated shear forces) can cause tissue deformation and ischemia during load bearing 78 activities1. These conditions can lead to cell death, tissue damage, and give rise to ulceration 79 and pain2. 80 Dermatological issues are experienced by 75% of individuals using lower limb 81 prostheses2,3 over their lifetime. These conditions lead to prosthesis disuse4 and can cause an 82 individual with BKA to become wheelchair-bound. Recent estimates suggest that 11-22% of 83 individuals abandon their prosthesis within one year of prescription5. These data are supported 84 by a report which found that 25% of users abandoned prosthetic limbs, with 29% citing 85 discomfort, 25% citing pain, and 12% claiming poor fit as the determining factor6. This 86 represents a substantial reduction in quality of life and increased healthcare-related financial 87 burden. 88 The socket is a crucial component of mobility and quality of life for individuals with BKA 89 due to its role as the interface between the human, prosthesis system, and gait environment. 90 An improved understanding of biomechanical interaction between the residual limb and 91 prosthetic socket during gait is necessary to attenuate rates of tissue damage and prosthesis 92 disuse. Previous experiments evaluating residual limb-socket interface dynamics have relied on 93 sensors integrated into the prosthetic socket7–12. However, previous systems have utilized bulky 94 sensors9,10,12, tethered cables8–11,13–16, or required modifications to the socket11,12, thereby 95 compromising the integrity of the socket interface and likely altering gait mechanics of 96 participants. 97 Simulations based on computational models may be useful for evaluating the 98 relationships between anatomical morphologies, gait mechanics, design of prosthesis systems, 99 and residual limb loading conditions. Previous model-based research has primarily employed 100 finite element (FE) analysis techniques to derive dynamic mathematical models of the residual 101 limb-socket interface1,11,17,18. While FE models allow for complex characterization of the 102 biological materials and their mechanical properties, they require individual-specific imaging 103 data as inputs, which is a cost and procedure not currently part of typical treatments. They are 104 also computationally costly and thus may be prohibitive when integrated with complex gait 105 models or models of other complex systems (e.g., biomechatronic rehabilitation devices). This 106 computational cost may restrict usability in clinical or personalized medicine scenarios. Other 107 studies have used abstract representations of the interface, such as an idealized joint 108 parameterized with spring and damper force laws19. These methods may be appropriate for 109 estimating generalized residual limb kinematics within the socket, but are unable to 110 differentiate limb-socket interaction forces and torques and offer little insight regarding relative 111 load distribution at different anatomical locations around the limb. Thus, there remains a need 112 for a computationally economical model of the biomechanical contact forces arising from the 113 residual limb-socket interaction during gait. 114 Manuscript Body This paper presents the design and development of a spatial contact force model 115 motivated by the material properties of the residual limb and prosthetic socket. The contact 116 model was integrated into a larger computational gait model to simulate kinematics and 117 kinetics at the socket interface during gait with a semi-active variable-stiffness prosthesis. We 118 simulated gait with three stiffness settings of the prosthesis, driven with subject-specific human 119 experiment data, to determine how foot parameters affect limb-socket interface dynamics. It 120 was hypothesized that the lowest limb-socket shear and pressure values would occur in the low 121 stiffness trials. Similarly, it was hypothesized that the lowest axial rotation and vertical 122 displacement of the residual limb inside the socket would occur in the low stiffness trials. This 123 model could assist experimental studies by providing insight into the effects of varied 124 prosthesis design parameters or gait conditions on residual limb-socket interface dynamics. 125 126 2.0 Materials and Methods 127 2.1 Model Design 128 A spatial contact model of the residual limb-socket interface was developed in Simscape 129 Multibody (Mathworks Inc, Natick, MA). The geometry of the residual limb bone element was 130 simplified as a rectangular cuboid with struts to represent the dimensions of the limb inclusive 131 of the soft tissue (Figure 1). Within the residual limb model, soft tissue and bone element 132 mechanics are not differentiated (i.e., the modeled dynamics are considered to be an aggregate 133 of soft tissue and bone mechanics). The prosthetic socket was modeled as a rigid hollow square 134 cone with an aperture of 100 deg. The residual limb interfaces with the socket at the same 135 angle. The residual limb and socket have nine interface frames with attached cuboid structures 136 to model interface dynamics. The shape of this model simplifies the continuous geometry of the 137 residuum and prosthetic socket, while allowing resolution of force distribution among the 138 different interface frames and aspects of the residuum. 139 The mass of the residual limb was estimated by deriving estimated density of the intact 140 limb, modeled as a conical frustum with an assigned mass estimated per De Leva (1996)20. The 141 derived density of the intact limb was applied to the residual limb model, which was then 142 truncated at the respective level of amputation for each subject (Figure 1). The mass of the 143 residuum was distributed evenly as point masses about the nine interface frames. The 144 residuum has two contact interfaces (one proximal, one distal) on each of the four sides of the 145 cuboid. The ninth interface is between the distal limb and base of the prosthetic socket. The 146 distance between the distal residual limb and base of the socket was assumed to be 1.5 cm21, 147 representative of an air gap, which is common between the socket and liner in prosthetic 148 socket systems21. Figure 1 depicts the rotationally symmetrical model. 149 Contact forces at the interfaces between the limb and socket in the normal plane are 150 implemented as modified Kelvin-Voigt material models with progressive spring and damper 151 characteristics. Shear stresses between the socket and residual limb are considered analogous 152 to the frictional forces arising from these interactions. In total, the residual limb has 4 degrees 153 of freedom (DoF) with respect to the prosthetic socket, including vertical translation (i.e., 154 pistoning) and rotations about three axes. 155 156 2.2 Model Parameterization 157 Manuscript Body A Kelvin-Voigt material model (spring and damper force law) was implemented to 158 estimate residual limb-prosthetic socket interaction forces. The model estimates normal (Fn) 159 and frictional forces (Ff) associated with the contact between a viscoelastic residual limb (spring 160 and damper system) and rigid prosthetic socket (eq. 1). The present model does not include a 161 socket-liner interface, but one could be implemented in the future. The spring force (k) 162 increases as a function of penetration depth (!), whereas damping force (b) increases with 163 penetration velocity (!)̇. Damping force is only applied when ! ̇> 0. Frictional forces are 164 calculated as the product of normal force and a user-defined coefficient of friction (µ) (eq. 2). A 165 stick-slip friction law defines the transition between static (µstatic) and kinetic (µkinetic) 166 coefficients of friction based on a velocity threshold (vthresh). Forces are applied along a common 167 contact plane and conform to Newton’s Third Law of Motion. 168 Values for spring stiffness in the normal plane (kn) were formulated according to 169 Hooke’s Law (eq. 3), as described in Zheng et al. (1999)22 and Noll et al. (2017)23. The effective 170 tissue moduli for individuals with a below-knee LLA (Table 1) were previously described in 171 Zheng et al. (1999)22 and Mak et al. (1994)24. In both studies, Poisson’s ratio was assumed to be 172 0.45. The stiffness values were parameterized independently for the anterior, posterior, medial, 173 lateral, and distal contact interfaces to best represent the varying moduli at corresponding 174 anatomical locations. Due to a lack of information reported in previous literature, damping 175 coefficients (N×s/mm) were set to half the numerical value of stiffness (N/mm) in an effort to 176 reduce high frequency oscillations at the contact interfaces, which can lead to rapidly-evolving 177 ordinary differential equations, and thus computational instability when simulating interface 178 dynamics. 179 The static coefficient of friction (µstatic) between the limb and socket was assigned a 180 value of 0.5, based on an in vivo study of the interaction between silicone rubber (a commonly 181 used material for prosthetic socket liners) and the skin of the human leg26. Coefficients of 182 friction between 0.5 and 3.0 have been reported for various other socket liner materials27,28. 183 The dynamic coefficient of friction (µdynamic) was set to 70% of the µstatic based on Cagle et al. 184 (2018)29. A velocity threshold (µvth) of 0.005 m/s defines the transition between the two values. 185 In the future, subject-specific values for µstatic and µdynamic could be implemented. 186 The model predicts normal pressure and shear stress at the contact interfaces. Based on 187 these forces, relative kinematics between the residual limb and prosthetic socket are simulated. 188 Model-derived estimates may be compared to the range of experimental values reported in the 189 literature for pressure, shear stress, and residual limb kinematics7,8,10,14,17,19,30,31. Previously-190 reported peak values for normal pressure range from 40-160 kPa8,15,30, and peak values for 191 shear stress range from 3-50 kPa7,10,11,23. The broad range of values in the literature may be 192 attributed to variation in sensors used, sensor placement, socket materials, individual-specific 193 residual limb tissue properties, and experimental gait protocols. Values should vary based on 194 phase of the gait cycle and anatomical location1,7,8,23. Nevertheless, values within these ranges 195 may be used as target criteria. Values of 1.0 to 4.2 cm have been reported for relative vertical 196 translation (i.e., residual limb pistoning)19,31–33. These values may vary based on residual limb 197 shape 34 and type of socket liner used31–33 and whether a vacuum or pin lock is incorporated in 198 the prosthesis attachment. Reported values for axial internal/external rotation (rotation about 199 the residual limb’s long axis) are between 0 and ±20 deg during gait19. 200 Manuscript Body 201 2.3 Gait Simulations with Experimental Data 202 The spatial contact model of the residual limb-socket interface was integrated into a gait 203 model with a BKA and a semi-active variable stiffness foot (VSF) prosthesis. This model was 204 previously described in McGeehan et al. (2021a and 2021b)35,36. Briefly, gait simulations were 205 driven by experimental data from four individuals (Table 2) walking with the VSF37 configured to 206 “high”, “medium”, and “low” stiffness settings, corresponding to forefoot stiffness values of 10, 207 19, and 32 N/mm. Forward kinematics simulations were computed for three trials per setting, 208 as in McGeehan et al. (2021b)36. Subject 2 did not complete one medium stiffness trial, and 209 Subject 4 did not complete one high stiffness trial. In total, 34 simulations were computed. All 210 simulations were computed in Simscape Multibody using the ode23t solver profile with variable 211 step sizes for numerical integration. 212 The experimental motion capture data used to drive the model were insufficient to 213 estimate kinematics of the residual limb with respect to the prosthetic socket. As such, data 214 from the literature were used to constrain motion of the residual limb via a bearing joint. A 215 progressive spring and damper force law was used to constrain motion. Limits of +0.5 cm 216 (upward displacement) and -3.5 cm (downward displacement) were imposed for residual limb 217 vertical translation19,31–33,38. Constraints of ±10 deg, ±5 deg, and ±5 deg were imposed for axial 218 rotation, anterior-posterior (AP), and medial-lateral (ML) rotations19,39. These constraints were 219 necessary to account for the differences in contact surface areas between the discretized 220 contact model and the continuous geometry of a biological residual limb and its mechanical 221 interface with a prosthetic socket. 222 223 2.4 Data Processing and Statistical Methodology 224 Kinematic and kinetic signals related to the simulated limb-socket interface were low-225 pass filtered via zero-lag a 4th order Butterworth (fc: 6Hz). Data were indexed from heel strike to 226 toe-off and resampled to 101 data points via cubic spline interpolation. These methods allow 227 for comparison of stance phases of different lengths. Model-derived values were compared to 228 those previously reported in the literature and the effects of stiffness setting on limb-socket 229 dynamics were evaluated using linear mixed effects (LME) regression analysis. It was 230 hypothesized that these effects would be subject-dependent, and as such, an exemplary case 231 study for Subject 1 is presented along with group level data. 232 We computed discrete outcome metrics from the simulations. Interfacial normal forces, 233 frictional forces, pistoning displacement, and axial rotations were derived from the spatial 234 contact model for the duration of stance phase. From the simulated forces, normalized 235 pressure values were calculated for each of the nine interface frames based on the surface area 236 of the interface and the body weight (BW) of the participant (kPa/BW). Pressure distributions 237 between opposing frames (AP, ML, and proximal-distal (PD)) were calculated as the percent 238 contribution of each contact interface to the sum of the pressure of both opposing interfaces. 239 Peak values for each interface and peak average values (i.e., peak of the average value for the 240 nine interfaces) were reported and used as dependent variable for subsequent analyses. 241 The relationships between the aforementioned outcome variables and prosthetic foot 242 stiffness condition were addressed using LME regression. For each analysis, the peak of the 243 simulated outcome variable was the dependent variable, stiffness setting was the independent 244 Manuscript Body variable, and subject was a random effect. We computed the least squares coefficients to the 245 linear mixed effects model (metric = ß0*stiffness + Oi), along with confidence intervals and 246 p-values for each coefficient. The random effect appears as a unique intercept (Oi) for each 247 subject. The overall effect of stiffness on each outcome measure (coefficient ß0) was 248 evaluated (α = 0.05). We report the slope, mean intercept, adjusted R2, and p-value for each 249 model. The slope (ß0) represents the sensitivity of the outcome to changes in stiffness, 250 whereas the mean intercept (mean of all Oi) allows for characterization of the group 251 average values. The adjusted R2 value represents the strength of the independent variable 252 to explain variations in the dependent variable, excluding the effect of individual intercepts. 253 P-value quantifies the statistical confidence that the slope (sensitivity) is different from 254 zero. These methods were adapted from previous work evaluating the effects of prosthesis 255 design on gait outcomes40–42. Analyses were performed using MATLAB 2020b (Mathworks 256 Inc., Natick, MA). 257 258 3.0 Results 259 3.1 Model Performance (Group Data) 260 Contact model-derived values for normal pressure and shear stress were dependent 261 upon anatomical location (Table 3) and progression of stance phase (Figures 2 and 3). Peak 262 average normal pressure across stance phase was 70.4 ± 4.28, 75.9 ± 4.44, and 85.0 ± 13.0 kPa 263 for the low, medium, and high stiffness conditions, whereas peak average shear stress values 264 were 25.0 ± 1.52, 26.6 ±1.55, and 29.9 ± 4.61 kPa for the same conditions (Figure 2, Table 5). 265 Increased prosthesis stiffness was associated with increased peak average normal pressure and 266 shear stress (p < 0.05) (Figure 4, Table 3); specific increases were observed at the anterior 267 proximal and anterior distal interfaces (p < 0.05). Stiffness did not significantly affect normal 268 pressure or shear stress at the other interface frames, nor did it affect residual limb pistoning or 269 axial rotation (p > 0.05), though effects were subject-dependent (Table 3). 270 Spatiotemporal patterns for pressure and shear stress distribution were variable between 271 participants, but show similar variability across stiffness conditions (Figure 3, Tables 4 and 5). On 272 average, participants displayed predominantly anterior pressure and shear distributions early in 273 stance phase, but trended toward a more even distribution later in stance phase. Pressure 274 trended slightly toward the lateral and proximal aspects of the residual limb thoughout stance 275 phase. High variability was observed among participants for the AP and ML distributions 276 throughout stance phase. 277 278 3.2 Case Study (Subject 1) 279 280 Mean data for Subject 1 demonstrated less variability compared to the group data 281 (Figures 5 and 6). Average pressure and shear stress values peaked at approximately 20% 282 stance phase, whereas residual limb pistoning peaked and plateaued near 50% stance phase. 283 Maximal pistoning displacement was approximately 1.5 cm for all conditions. A slightly 284 increased rate of pistoning appeared between 15-40% stance phase for the high stiffness 285 compared to the low and medium stiffness conditions (Figures 5 and 6). The subject’s residuum 286 was predominately externally rotated with respect to the prosthetic socket throughout stance 287 phase with maximal external rotation occurring near 50% stance phase. The low stiffness 288 Manuscript Body condition resulted in the least external rotation, though high variability was observed late in 289 stance phase for all conditions. 290 The effects of variable prosthesis stiffness on pressure and shear stress distribution 291 appeared primarily during 50-100% of stance phase (Figure 6). However, divergent patterns in 292 the mean data were accompanied by greater variability during this time. From 0-50% stance 293 phase, pressure and shear stress were weighted more heavily toward the anterior aspect of the 294 residual limb for all stiffness conditions. For the low and high stiffness conditions, mean 295 pressure and shear stress trended toward a relatively even AP distribution late in stance phase, 296 whereas the medium stiffness condition resulted in a relatively posterior distribution. Pressure 297 and shear stress distribution outcomes were similar across stiffness conditions for the ML and 298 proximal-distal aspects of the residuum. 299 300 4.0 Discussion 301 4.1 Model Performance (Group Data) 302 The objective of this study was to develop a spatial contact model of the residual limb-303 prosthetic socket interface and evaluate its ability to estimate limb-socket interface dynamics. 304 A secondary objective of this study was to use this model to examine the relationships between 305 prosthetic foot stiffness and limb-socket dynamics. The hypothesis that limb-socket normal 306 pressure and shear stresses would be lowest in the low stiffness condition was supported for 307 the average normal force and shear stress metrics. Specific effects were observed at the 308 anterior proximal and anterior distal interfaces. Stiffness did not significantly affect normal 309 pressure or shear stress at the other interface frames, nor did it affect residual limb kinematics. 310 Average normal pressure and shear stress increased with prosthesis stiffness during 311 stance phase. This pattern may indicate the need for greater force transfer from the residuum 312 to the prosthetic socket to load a stiffer prosthetic foot. This is supported by the specific 313 increases observed at the anterior interfaces. Participants may have adopted a strategy 314 whereby they increased the anterior forces in the socket to deform the energy storage and 315 return structures of the VSF under stiffer configurations. Increased prosthesis stiffness has been 316 shown to offer potential biomechanical benefits such as increased mechanical efficiency43 and 317 increased propulsion43. However, these data suggest that prosthesis users may in part achieve 318 this through increased loading of the limb-socket interface, which may increase risk of tissue 319 damage via friction or pressure ulcers. 320 The pressure distribution profiles derived from this model were weighted toward the 321 anterior and lateral aspects of the residual limb. These patterns may be due to gait kinematics 322 of the participants. Future clinical gait evaluations should seek to verify these patterns. 323 Model-derived values for normal pressure depict spatiotemporal patterns similar to 324 those of a ground reaction force (GRF) curve during stance phase. The pressure and shear 325 waveforms presented by Sanders et al. (1992)45 and Laszczak et al. (2016)46 are similar to those 326 of the present study early in stance phase, but exhibit a brief plateau during mid-stance before 327 values decrease. Comparatively, the present study shows similar loading rates, but a gradual 328 decline in pressure and stress rather than a mid-stance phase plateau (i.e., the waveforms are 329 skewed toward early stance phase) (Figure 2). The lack of a plateau in pressure and shear data 330 in the present results may reflect the different prostheses used for experimental gait trials. 331 Mathematically, it could also be due to inadequate constraining of residual limb motion. Using 332 Manuscript Body experimental kinematic data to constrain residual limb motion may help refine the trajectory of 333 the modeled response, or a velocity constraint could be implemented into the present model 334 design. Peak values for normal pressure were similar to those reported in previous sensor-335 based experiments8,10,30,44 and finite element modeling analyses1,17,18, which ranged from 40-336 160 kPa. Values for shear stress were within the 3-50 kPa range reported previously7,10,11,23. 337 The accuracy of the model to predict pressure and shear stress values specific to 338 different anatomical locations is difficult to discern based on previous experiments. Sensor-339 based experiments have typically sampled from small localized areas on the limb or present 340 only resultant data. Nevertheless, broad comparisons can be made with data from Sanders et 341 al. (1992 and 2000)30,45 and Courtney et al. (2016)8. Results of this study showed peak mean 342 pressure and stress values on the medial side of 60 and 21 kPa (values are the mean of the 343 proximal and distal interfaces under all stiffness conditions). Courtney et al. (2016) reported 344 peak medial pressure of approximately 65-70 kPa, whereas Sanders et al. (2000) present values 345 ranging 40-85 kPa for pressure and 7-12 kPa for shear stress. On the lateral side, results of this 346 study showed peak mean pressure and shear stress values of 117 and 41 kPa. Comparatively, 347 Sanders et al. (2000) present values of 60-140 kPa for pressure and 18-23 kPa for shear stress. 348 Posteriorly, the pressure and shear values of 41 and 13 kPa were lower compared to those 349 presented by Sanders et al. (85-100 and 17-22 kPa). This discrepancy may be due to the 350 increased stiffness of the tissue on the posterior residual limb associated with muscle 351 contraction during gait, which is unaccounted for in this model. Muscular contraction has been 352 shown to increase tissue modulus, for example by 45 kPa22 in the muscles of the forearm. 353 Muscular contraction would likely have minimal effects on the frictional characteristics of the 354 tissue. In the future, a progressive model of tissue moduli could be implemented into a limb-355 socket contact model. The model’s predicted anterior pressure and shear stress were 235 and 356 83 kPa, which were similar to values of 245 and 105 derived from FEA of socket interface 357 dynamics at the patellar tendon18. There is a paucity of data from sensor-based experiments 358 related to pressure or shear dynamics along the anterior tibia. 359 The model predicted peak residual limb displacements between 1.2 and 2.1 cm with 360 respect to the socket. These values are within the range of 1.1-3.6 cm (mean: 2.3 ± 1.0 cm) 361 previously reported in the literature19,32–34,47–49. It should be noted that the prosthetic socket 362 components (e.g., liner and socket materials) and gait tasks varied among these studies. Data 363 from the present study, among others, support the idea that the amount of residual limb 364 pistoning may be affected by liner and socket type 30,50, residual limb shape 34, and gait 365 conditions30,33. Data regarding the prosthetic socket componentry used by participants in this 366 study were not available. 367 368 4.2 Case Study (Subject 1) 369 Across the stiffness conditions, data from Subject 1 showed similar spatiotemporal 370 patterns between 0-50% stance phase (Figure 5), which encompasses the progression from heel 371 strike to foot flat51. Divergent responses were observed across the stiffness conditions in the 372 latter half of stance phase for residual limb axial angle and AP pressure and shear stress 373 distribution. Increased variability was also observed for all conditions during this time. The 374 latter half of stance phase is characterized by the progression from foot flat to toe off and 375 involves an anterior shift in the center of pressure (CoP)51. Since the stiffness behavior of the 376 Manuscript Body VSF’s heel is unchanged across the conditions, it is logical that the effects of variable forefoot 377 stiffness would be primarily observed in the latter half of stance phase. 378 The subject presented no discernable effect of variable stiffness on the peak average 379 values for normal pressure, shear stress, or piston motion (Figure 5). Decreased external 380 rotation was observed in the low stiffness condition, and high variability was present in the high 381 stiffness condition. Increased external rotation may direct knee loading out of the sagittal plane 382 and into the frontal plane. Lack of range of motion in the frontal plane compromises the ability 383 of muscles to support load; instead the joint relies on passive structures (ligaments and 384 cartilage). This effect could overload these structures52. Increased external rotation has been 385 associated with high rates of medial knee osteoarthritis53, which has been documented for both 386 the amputated and contralateral limb of lower limb prosthesis users52. At the limb-socket 387 interface, there were no discernable effects of prosthesis stiffness on distribution of frontal 388 plane pressures or shear stresses. This response is consistent with the mechanical principles of 389 the VSF, which modulates forefoot stiffness primarily in the sagittal plane. 390 391 4.4 Limitations and Future Directions 392 The present model was parameterized using previously reported residual limb tissue 393 mechanical properties and limb-socket kinematics for individuals with BKA. While these 394 methods resulted in pressure and shear stress values within the range reported in the 395 literature, variability in the aforementioned parameters is well documented between individual 396 subjects. Future work should strive toward individualized models by characterizing the unique 397 tissue mechanical properties of subjects. This could be accomplished by measuring tissue 398 stiffness at corresponding sites between the residuum and model and using those 399 measurements to parameterize the model. Similarly, adjusting parameters based on the socket 400 componentry used by subjects would improve the accuracy of the model. For example, 401 coefficients of friction between 0.5-3.0 have been reported for the interaction between human 402 skin and various socket liner materials27,28. Variation within this range would have a substantial 403 impact on model-derived shear stresses estimates. Further, this model does not account for 404 frictional forces at the liner-socket interface. If frictional coefficients are lower for this interface 405 compared to the skin-liner interface, inaccuracies in modeled shear stresses would arise. 406 Future work should also seek to quantify kinematics between the residual limb and 407 prosthetic socket through optical motion capture or instrumenting participants with 408 potentiometers. Using these data to constrain residual limb motion during simulations would 409 improve accuracy on an individualized basis. These data could be used to refine the ability of 410 the current model to predict limb-socket kinematics. 411 The present model assumes oversimplified geometries of the residuum and prosthetic 412 socket. Developing more complex interface geometry could improve model fidelity. For 413 example, using a pentagonal prism shape to model the residual limb geometry would allow for 414 differentiation of the varying moduli of the anterior, anterior lateral, anterior medial, posterior 415 lateral, and posterior medial aspects of the residual limb and would only add two interface 416 frames compared to the present model. 417 This study modeled the residual limb as composite of both the bone and soft tissue 418 elements. However, data from x-ray54 and biplane fluoroscopy55 studies suggest that residual 419 limb kinematics can be differentiated into the motion of the bone and soft tissue elements. As 420 Manuscript Body such, it may be important to distinguish these elements and model the interface between them 421 in future studies. Doing so could lead to improvements when simulating limb-socket dynamics. 422 423 5.0 Conclusions 424 Findings from these simulations support the use of reduced order modeling techniques 425 to estimate residual limb-prosthetic socket interfacial pressure and shear stress, as well as 426 residual limb kinematics in a computationally economical manner. Residual limb-prosthetic 427 socket interface dynamics derived from this model were within the range of values reported by 428 previous sensor-based gait experiments. These methods may be useful to aid experimental 429 studies by providing insights into the effects of varied prosthesis design parameters or gait 430 conditions on residual limb-socket interface dynamics. 431 Simulated data showed increased peak average values for normal pressure and shear 432 stress with a stiffer prosthesis; specific effects were observed on the anterior aspect of the 433 residual limb-socket interface. Data from a case study show promise for evaluating 434 individualized effects of prosthesis stiffness on limb-socket dynamics. Future work could add 435 complexity to the modeled interface geometry in order to better match the shape and variation 436 in tissue material properties of the residual limb. Additionally, the model’s accuracy could be 437 improved by applying subject-specific data for residual limb tissue properties and prosthetic 438 socket componentry when parameterizing the contact interfaces. 439