Abstract
Keywords
Introduction
Rehabilitation following a lower limb amputation (LLA) often includes prescription of a prosthesis designed to replace the functionality of the removed limb. For an individual with a below-knee amputation (BKA), a prosthesis system typically consists of a socket, which interfaces with the residual limb, a rigid pylon, and a foot-ankle prosthesis. Use of lower limb prostheses can improve mobility, health, and quality of life. However, abnormal loading of the soft tissues surrounding the truncated shank (e.g. asymmetric pressure distribution and elevated shear forces) can cause tissue deformation and ischemia during load bearing activities. 1 These conditions can lead to cell death, tissue damage, and give rise to ulceration and pain. 2
Dermatological issues are experienced by 75% of individuals using lower limb prostheses2,3 over their lifetime. These conditions lead to prosthesis disuse 4 and can cause an individual with BKA to become wheelchair-bound. Recent estimates suggest that 11–22% of individuals abandon their prosthesis within 1 year of prescription. 5 These data are supported by a report which found that 25% of users abandoned prosthetic limbs, with 29% citing discomfort, 25% citing pain, and 12% claiming poor fit as the determining factor. 6 This represents a substantial reduction in quality of life and increased healthcare-related financial burden.
The socket is a crucial component of mobility and quality of life for individuals with BKA due to its role as the interface between the human, prosthesis system, and gait environment. An improved understanding of biomechanical interaction between the residual limb and prosthetic socket during gait is necessary to attenuate rates of tissue damage and prosthesis disuse. Previous experiments evaluating residual limb-socket interface dynamics have relied on sensors integrated into the prosthetic socket.7–12 However, previous systems have utilized bulky sensors,9,10,12 tethered cables,8–11,13–16 or required modifications to the socket,11,12 thereby compromising the integrity of the socket interface and likely altering gait mechanics of participants.
Simulations based on computational models may be useful for evaluating the relationships between anatomical morphologies, gait mechanics, design of prosthesis systems, and residual limb loading conditions. Previous model-based research has primarily employed finite element (FE) analysis techniques to derive dynamic mathematical models of the residual limb-socket interface.1,11,17–19 While FE models allow for complex characterization of the biological materials and their mechanical properties, they require individual-specific imaging data as inputs, which is a cost and procedure not currently part of typical treatments. They are also computationally costly and thus may be prohibitive when integrated with complex gait models or models of other complex systems (e.g. biomechatronic rehabilitation devices). This computational cost may restrict usability in clinical or personalized medicine scenarios. Other studies have used abstract representations of the interface, such as an idealized joint parameterized with spring and damper force laws. 20 These methods may be appropriate for estimating generalized residual limb kinematics within the socket, but are unable to differentiate limb-socket interaction forces and torques and offer little insight regarding relative load distribution at different anatomical locations around the limb. Thus, there remains a need for a computationally economical model of the biomechanical contact forces arising from the residual limb-socket interaction during gait.
This paper presents the design and development of a spatial contact force model motivated by the material properties of the residual limb and prosthetic socket. The contact model was integrated into a larger computational gait model to simulate kinematics and kinetics at the socket interface during gait with a semi-active variable-stiffness prosthesis. We simulated gait with three stiffness settings of the prosthesis, driven with subject-specific human experiment data, to determine how foot parameters affect limb-socket interface dynamics. It was hypothesized that the lowest limb-socket shear and pressure values would occur in the low stiffness trials. Similarly, it was hypothesized that the lowest axial rotation and vertical displacement of the residual limb inside the socket would occur in the low stiffness trials. This model could assist experimental studies by providing insight into the effects of varied prosthesis design parameters or gait conditions on residual limb-socket interface dynamics.
Materials and methods
Model design
A spatial contact model of the residual limb-socket interface was developed in Simscape Multibody (Mathworks Inc, Natick, MA). The geometry of the residual limb bone element was simplified as a rectangular cuboid with struts to represent the dimensions of the limb inclusive of the soft tissue (Figure 1). Within the residual limb model, soft tissue and bone element mechanics are not differentiated (i.e. the modeled dynamics are considered to be an aggregate of soft tissue and bone mechanics). The prosthetic socket was modeled as a rigid hollow square cone with an aperture of 100 deg. The residual limb interfaces with the socket at the same angle. The residual limb and socket have nine interface frames with attached cuboid structures to model interface dynamics. The shape of this model simplifies the continuous geometry of the residuum and prosthetic socket, while allowing resolution of force distribution among the different interface frames and aspects of the residuum, thus allowing the model to simulate clinically-relevant outcomes (e.g. pressure distribution within the socket) while remaining computationally efficient. Depiction of the rotationally symmetrical residual limb and socket geometries, including the nine interface frames.
The mass of the residual limb was estimated by deriving estimated density of the intact limb, modeled as a conical frustum with an assigned mass estimated per De Leva (1996). 21 The derived density of the intact limb was applied to the residual limb model, which was then truncated at the respective level of amputation for each subject (Figure 1). The mass of the residuum was distributed evenly as point masses about the nine interface frames. The residuum has two contact interfaces (one proximal, one distal) on each of the four sides of the cuboid. The ninth interface is between the distal limb and base of the prosthetic socket. The distance between the distal residual limb and base of the socket was assumed to be 1.5 cm, 22 representative of an air gap, which is common between the socket and liner in prosthetic socket systems. 22 Figure 1 depicts the rotationally symmetrical model.
Contact forces at the interfaces between the limb and socket in the normal plane are implemented as modified Kelvin-Voigt material models with progressive spring and damper characteristics. Shear stresses between the socket and residual limb are considered analogous to the frictional forces arising from these interactions. In total, the residual limb has 4 degrees of freedom (DoF) with respect to the prosthetic socket, including vertical translation (i.e. pistoning) and rotations about three axes.
Model parameterization
A Kelvin-Voigt material model (spring and damper force law) was implemented to estimate residual limb-prosthetic socket interaction forces. The model estimates normal (
Summary of Young’s Modulus values for various anatomical locations on the residual limb.
The static coefficient of friction (μstatic) between the limb and socket was assigned a value of 0.5, based on an in vivo study of the interaction between silicone rubber (a commonly used material for prosthetic socket liners) and the skin of the human leg. 26 Coefficients of friction between 0.5 and 3.0 have been reported for various other socket liner materials.27,28 The dynamic coefficient of friction (μdynamic) was set to 70% of the μstatic based on Cagle et al. (2018). 28 A velocity threshold (μvth) of 0.005 m/s defines the transition between the two values. In the future, subject-specific values for μstatic and μdynamic could be implemented.
The model predicts normal pressure and shear stress at the contact interfaces. Based on these forces, relative kinematics between the residual limb and prosthetic socket are simulated. Model-derived estimates may be compared to the range of experimental values reported in the literature for pressure, shear stress, and residual limb kinematics.7,8,10,14,17,20,29,30 Previously-reported peak values for normal pressure range from 40-160 kPa,8,15,29 and peak values for shear stress range from 3-50 kPa.7,10,11,24 Measuring pressures (Pascals) compared to forces (Newtons) may be more clinically-relevant as it accounts for the variations in surface area between individuals. Previous work has found tissue damage with loads 4–23 kPa. 19 The broad range of values in the literature may be attributed to variation in sensors used, sensor placement, socket materials, individual-specific residual limb tissue properties, and experimental gait protocols. Values should vary based on phase of the gait cycle and anatomical location.1,7,8,24 Nevertheless, values within these ranges may be used as target criteria. Values of 1.0–4.2 cm have been reported for relative vertical translation (i.e. residual limb pistoning).20,30–32 These values may vary based on residual limb shape 33 and type of socket liner used30–32 and whether a vacuum or pin lock is incorporated in the prosthesis attachment. Reported values for axial internal/external rotation (rotation about the residual limb’s long axis) are between 0 and ±20 deg during gait. 20
Gait simulations with experimental data
Participant characteristics.
The experimental motion capture data used to drive the model were insufficient to estimate kinematics of the residual limb with respect to the prosthetic socket. As such, data from the literature were used to constrain motion of the residual limb via a bearing joint. A progressive spring and damper force law was used to constrain motion. Limits of +0.5 cm (upward displacement) and −3.5 cm (downward displacement) were imposed for residual limb vertical translation.20,30–32,37 Constraints of ±10 deg, ±5 deg, and ±5 deg were imposed for axial rotation, anterior-posterior (AP), and medial-lateral (ML) rotations.20,38 These constraints were necessary to account for the differences in contact surface areas between the discretized contact model and the continuous geometry of a biological residual limb and its mechanical interface with a prosthetic socket.
Data processing and statistical methodology
Kinematic and kinetic signals related to the simulated limb-socket interface were low-pass filtered via zero-lag fourth order Butterworth (
We computed discrete outcome metrics from the simulations. Interfacial normal forces, frictional forces, pistoning displacement, and axial rotations were derived from the spatial contact model for the duration of stance phase. From the simulated forces, normalized pressure values were calculated for each of the nine interface frames based on the surface area of the interface and the body weight (BW) of the participant (kPa/BW). Pressure distributions between opposing frames (AP, ML, and proximal-distal (PD)) were calculated as the percent contribution of each contact interface to the aggregate pressure of both opposing interfaces. Peak values for each interface and peak average values (i.e. peak of the average value for the nine interfaces) were reported and used as dependent variable for subsequent analyses.
The relationships between the aforementioned outcome variables and prosthetic foot stiffness condition were addressed using LME regression. For each analysis, the peak of the simulated outcome variable was the dependent variable, stiffness setting was the independent variable, and subject was a random effect. We computed the least squares coefficients to the linear mixed effects model (
Results
Model performance (group data)
MLE statistical parameters.

Group mean data for normal pressure and shear stress (top) and residual limb pistoning and internal/external rotation (bottom) across stance phase for the low, medium, and high stiffness conditions. Mean pressure and shear stress values are the mean pressure and shear stress across the nine interface frames. Kinetic data are normalized to subject body weight. These data are also described in Table 4.

Group mean data for normal pressure (top) and shear stress (bottom) distributions in the anterior-posterior, posterior, medial-lateral, and proximal-distal directions. Pressure and shear stress distributions are the percent contribution of opposing interface frames (AP, ML, and PD) to the sum of the pressure or shear stress of the opposing the interfaces. These data are also described in Tables 5 and 6.
Peak average values for normal and shear stress (average of nine anatomical locations), and peak values for residual limb piston displacement with respect to the prosthetic socket. Data are mean ± SD.

LME regression responses for normal forces (top), shear stresses (middle), and kinematics (bottom).
Average peak values for normal pressure (kPa) by anatomical location, subject, and condition. A zero value for the distal contact interface implies that the distal tibia did not contact the base of the prosthetic socket (i.e., piston displacement < 1.5 cm). Data are mean ± SD.
Average peak values for shear stress (kPa) by anatomical location, subject, and condition. A zero value for the distal contact interface implies that the distal tibia did not contact the base of the prosthetic socket (i.e., piston displacement < 1.5 cm). Data are mean ± SD.
Case study (subject 1)
Mean data for Subject 1 demonstrated less variability compared to the group data (Figures 5 and 6). Average pressure and shear stress values peaked at approximately 20% stance phase, whereas residual limb pistoning peaked and plateaued near 50% stance phase. Maximal pistoning displacement was approximately 1.5 cm for all conditions. A slightly increased rate of pistoning appeared between 15-40% stance phase for the high stiffness compared to the low and medium stiffness conditions (Figures 5 and 6). The subject’s residuum was predominately externally rotated with respect to the prosthetic socket throughout stance phase with maximal external rotation occurring near 50% stance phase. The low stiffness condition resulted in the least external rotation, though high variability was observed late in stance phase for all conditions. Mean data for normal pressure, shear stress, residual limb pistoning, and residual limb internal/external rotation across stance phase for the low, medium, and high stiffness conditions. Kinetic data are normalized to subject body weight. These data are also described in Table 4. Subject 1 mean data for normal pressure (top) and shear stress (bottom) distributions in the anterior-posterior, medial-lateral, and proximal-distal directions. Pressure and shear stress distributions are the percent contribution of opposing interface frames (AP, ML, and PD) to the sum of the pressure or shear stress of the opposing the interfaces. These data are also described in Tables 5 and 6.

The effects of variable prosthesis stiffness on pressure and shear stress distribution appeared primarily during 50–100% of stance phase (Figure 6). However, divergent patterns in the mean data were accompanied by greater variability during this time. From 0-50% stance phase, pressure and shear stress were weighted more heavily toward the anterior aspect of the residual limb for all stiffness conditions. For the low and high stiffness conditions, mean pressure and shear stress trended toward a relatively even AP distribution late in stance phase, whereas the medium stiffness condition resulted in a relatively posterior distribution. Pressure and shear stress distribution outcomes were similar across stiffness conditions for the ML and proximal-distal aspects of the residuum.
Discussion
Model performance (group data)
The objective of this study was to develop a spatial contact model of the residual limb-prosthetic socket interface and evaluate its ability to estimate limb-socket interface dynamics. A secondary objective of this study was to use this model to examine the relationships between prosthetic foot stiffness and limb-socket dynamics. The hypothesis that limb-socket normal pressure and shear stresses would be lowest in the low stiffness condition was supported for the average normal force and shear stress metrics. Specific effects were observed at the anterior proximal and anterior distal interfaces. Stiffness did not significantly affect normal pressure or shear stress at the other interface frames, nor did it affect residual limb kinematics.
Average normal pressure and shear stress increased with prosthesis stiffness during stance phase. This pattern may indicate the need for greater force transfer from the residuum to the prosthetic socket to load a stiffer prosthetic foot. This is supported by the specific increases observed at the anterior interfaces. Participants may have adopted a strategy whereby they increased the anterior forces in the socket to deform the energy storage and return structures of the VSF under stiffer configurations. Increased prosthesis stiffness has been shown to offer potential biomechanical benefits such as increased mechanical efficiency 42 and increased propulsion. 42 However, these data suggest that prosthesis users may in part achieve this through increased loading of the limb-socket interface, which may increase risk of tissue damage via friction or pressure ulcers.
The pressure distribution profiles derived from this model were weighted toward the anterior and lateral aspects of the residual limb. These patterns may be due to gait kinematics of the participants. Future clinical gait evaluations should seek to verify these patterns.
Model-derived values for normal pressure depict spatiotemporal patterns similar to those of a ground reaction force curve during stance phase. The pressure and shear waveforms presented by Sanders et al. (1992) 43 and Laszczak et al. (2016) 10 are similar to those of the present study early in stance phase, but exhibit a brief plateau during mid-stance before values decrease. Comparatively, the present study shows similar loading rates, but a gradual decline in pressure and stress rather than a mid-stance phase plateau (i.e. the waveforms are skewed toward early stance phase) (Figure 2). The lack of a plateau in pressure and shear data in the present results may reflect the different prostheses used for experimental gait trials. Mathematically, it could also be due to inadequate constraining of residual limb motion. Using experimental kinematic data to constrain residual limb motion may help refine the trajectory of the modeled response, or a velocity constraint could be implemented into the present model design. Peak values for normal pressure were similar to those reported in previous sensor-based experiments8,10,29,43 and finite element modeling analyses,1,17,18 which ranged from 40-160 kPa. Values for shear stress were within the 3–50 kPa range reported previously.7,10,11,24
The accuracy of the model to predict pressure and shear stress values specific to different anatomical locations is difficult to discern based on previous experiments. Sensor-based experiments have typically sampled from small, localized areas on the limb or present only resultant data. Nevertheless, broad comparisons can be made with data from Sanders et al. (1992 and 2000)29,43 and Courtney et al. (2016). 8 Results of this study showed peak mean pressure and stress values on the medial side of 60 and 21 kPa (values are the mean of the proximal and distal interfaces under all stiffness conditions). Courtney et al. (2016) reported peak medial pressure of approximately 65–70 kPa, whereas Sanders et al. (2000) present values ranging 40–85 kPa for pressure and 7–12 kPa for shear stress. On the lateral side, results of this study showed peak mean pressure and shear stress values of 117 and 41 kPa. Comparatively, Sanders et al. (2000) present values of 60–140 kPa for pressure and 18–23 kPa for shear stress. Posteriorly, the pressure and shear values of 41 and 13 kPa were lower compared to those presented by Sanders et al. (85–100 and 17–22 kPa). This discrepancy may be due to the increased stiffness of the tissue on the posterior residual limb associated with muscle contraction during gait, which is unaccounted for in this model. Muscular contraction has been shown to increase tissue modulus, for example by 45 kPa 23 in the muscles of the forearm. Muscular contraction would likely have minimal effects on the frictional characteristics of the tissue. In the future, a progressive model of tissue moduli could be implemented into a limb-socket contact model. The model’s predicted anterior pressure and shear stress were 235 and 83 kPa, which were similar to values of 245 and 105 derived from FEA of socket interface dynamics at the patellar tendon. 18 There is a paucity of data from sensor-based experiments related to pressure or shear dynamics along the anterior tibia.
The model predicted peak residual limb displacements between 1.2 and 2.1 cm with respect to the socket. These values are within the range of 1.1–3.6 cm (mean: 2.3 ± 1.0 cm) previously reported in the literature.20,31–33,44–46 It should be noted that the prosthetic socket components (e.g. liner and socket materials) and gait tasks varied among these studies. Data from the present study, among others, support the idea that the amount of residual limb pistoning may be affected by liner and socket type,29,47 residual limb shape, 33 and gait conditions.29,32 Data regarding the prosthetic socket componentry used by participants in this study were not available.
Case study (subject 1)
Across the stiffness conditions, data from Subject 1 showed similar spatiotemporal patterns between 0-50% stance phase (Figure 5), which encompasses the progression from heel strike to foot flat. 48 Divergent responses were observed across the stiffness conditions in the latter half of stance phase for residual limb axial angle and AP pressure and shear stress distribution. Increased variability was also observed for all conditions during this time. The latter half of stance phase is characterized by the progression from foot flat to toe off and involves an anterior shift in the center of pressure. 48 Since the stiffness behavior of the VSF’s heel is unchanged across the conditions, it is logical that the effects of variable forefoot stiffness would be primarily observed in the latter half of stance phase.
The subject presented no discernable effect of variable stiffness on the peak average values for normal pressure, shear stress, or piston motion (Figure 5). Decreased external rotation was observed in the low stiffness condition, and high variability was present in the high stiffness condition. Increased external rotation may direct knee loading out of the sagittal plane and into the frontal plane. Lack of range of motion in the frontal plane compromises the ability of muscles to support load; instead, the joint relies on passive structures (ligaments and cartilage). This effect could overload these structures. 49 Increased external rotation has been associated with high rates of medial knee osteoarthritis, 50 which has been documented for both the amputated and contralateral limb of lower limb prosthesis users. 51 At the limb-socket interface, there were no discernable effects of prosthesis stiffness on distribution of frontal plane pressures or shear stresses. This response is consistent with the mechanical principles of the VSF, which modulates forefoot stiffness primarily in the sagittal plane.
Limitations and future directions
The present model was parameterized using previously reported residual limb tissue mechanical properties and limb-socket kinematics for individuals with BKA. While these methods resulted in pressure and shear stress values within the range reported in the literature, variability in the aforementioned parameters is well documented between individual subjects. Future work should strive toward individualized models by characterizing the unique tissue mechanical properties of subjects. This could be accomplished by measuring tissue stiffness at corresponding sites between the residuum and model and using those measurements to parameterize the model. Similarly, adjusting parameters based on the socket componentry used by subjects would improve the accuracy of the model. For example, coefficients of friction between 0.5-3.0 have been reported for the interaction between human skin and various socket liner materials.27,28 Variation within this range would have a substantial impact on model-derived shear stresses estimates. Further, this model does not account for frictional forces at the liner-socket interface. If frictional coefficients are lower for this interface compared to the skin-liner interface, inaccuracies in modeled shear stresses would arise.
Future work should also seek to quantify kinematics between the residual limb and prosthetic socket through optical motion capture or instrumenting participants with potentiometers. Using these data to constrain residual limb motion during simulations would improve accuracy on an individualized basis. These data could be used to refine the ability of the current model to predict limb-socket kinematics.
The present model assumes oversimplified geometries of the residuum and prosthetic socket. Developing more complex interface geometry could improve model fidelity. For example, using a pentagonal prism shape to model the residual limb geometry would allow for differentiation of the varying moduli of the anterior, anterior lateral, anterior medial, posterior lateral, and posterior medial aspects of the residual limb and would only add two interface frames compared to the present model. Additionally, personalized limb geometry models could be derived using 3D surface geometry and scanning tools and shape analysis software.51,52 This could allow the modeled shape to be more reflective of the in vivo socket-limb interface (e.g. the triangular shape of the tibia).
This study modeled the residual limb as composite of both the bone and soft tissue elements. However, data from X-ray 53 and biplane fluoroscopy 54 studies suggest that residual limb kinematics can be differentiated into the motion of the bone and soft tissue elements. As such, it may be important to distinguish these elements and model the interface between them in future studies. Doing so could lead to improvements when simulating limb-socket dynamics.
Conclusions
Findings from these simulations support the use of reduced order modeling techniques to estimate residual limb-prosthetic socket interfacial pressure and shear stress, as well as residual limb kinematics in a computationally economical manner. Residual limb-prosthetic socket interface dynamics derived from this model were within the range of values reported by previous sensor-based gait experiments. These methods may be useful to aid experimental studies by providing insights into the effects of varied prosthesis design parameters or gait conditions on residual limb-socket interface dynamics.
Simulated data showed increased peak average values for normal pressure and shear stress with a stiffer prosthesis; specific effects were observed on the anterior aspect of the residual limb-socket interface. Data from a case study show promise for evaluating individualized effects of prosthesis stiffness on limb-socket dynamics. Future work could add complexity to the modeled interface geometry in order to better match the shape and variation in tissue material properties of the residual limb. Additionally, the model’s accuracy could be improved by applying subject-specific data for residual limb tissue properties and prosthetic socket componentry when parameterizing the contact interfaces.
Equations
Footnotes
Declaration of conflicting interests
The author(s) declared the following potential conflicts of interest with respect to the research, authorship, and/or publication of this article: Michael McGeehan, Kieran Nichols, and Michael Hahn declare no conflicts of interest for this work. Peter Adamczyk was a partner in Intelligent Prosthetic Systems, LLC, which created the Variable Stiffness Prosthesis used in this project.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Lokey Doctoral Science Fellowship (MAM) and NIH grant HD074424 (PGA).
Guarantor
Dr. Mike Hahn
