A simulation-based analysis of the effects of variable prosthesis stiffness on interface dynamics between the prosthetic socket and residual limb.

Introduction: Loading of a residual limb within a prosthetic socket can cause tissue damage such as ulceration. Computational simulations may be useful tools for estimating tissue loading within the socket, and thus provide insights into how prosthesis designs affect residual limb-socket interface dynamics. The purpose of this study was to model and simulate residual limb-socket interface dynamics and evaluate the effects of varied prosthesis stiffness on interface dynamics during gait.
Methods: A spatial contact model of a residual limb-socket interface was developed and integrated into a gait model with a below-knee amputation. Gait trials were simulated for four subjects walking with low, medium, and high prosthesis stiffness settings. The effects of prosthesis stiffness on interface kinematics, normal pressure, and shear stresses were evaluated.
Results: Model-predicted values were similar to those reported previously in sensor-based experiments; increased stiffness resulted in greater average normal pressure and shear stress (p < 0.05). Conclusions: 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. The current results suggest that these effects may be subject-specific.

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. These conditions can lead to cell death, tissue damage, and give rise to ulceration and pain.

Dermatological issues are experienced by 75% of individuals using lower limb prostheses[2,3] over their lifetime. These conditions lead to prosthesis disuse 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. 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. 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. 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.

Figure 1.

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). 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, representative of an air gap, which is common between the socket and liner in prosthetic socket systems. 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 (F ) and frictional forces (F ) associated with the contact between a viscoelastic residual limb (spring and damper system) and rigid prosthetic socket (eq. 1). The present model does not include a socket-liner interface, but one could be implemented in the future. The spring force (k) increases as a function of penetration depth , whereas damping force (b) increases with penetration velocity . Damping force is only applied when > 0. Frictional forces are calculated as the product of normal force and a user-defined coefficient of friction (μ) (eq. 2). A stick-slip friction law defines the transition between static (μstatic) and kinetic (μkinetic) coefficients of friction based on a velocity threshold (vthresh). Forces are applied along a common contact plane and conform to Newton’s Third Law of Motion.

Values for spring stiffness in the normal plane (kn) were formulated according to Hooke’s Law (eq. 3), as described in Zheng et al. (1999) and Noll et al. (2017). The effective tissue moduli for individuals with a below-knee LLA (Table 1) were previously described in Zheng et al. (1999) and Mak et al. (1994). In both studies, Poisson’s ratio was assumed to be 0.45. The stiffness values were parameterized independently for the anterior, posterior, medial, lateral, and distal contact interfaces to best represent the varying moduli at corresponding anatomical locations. Due to a lack of information reported in previous literature, damping coefficients (N⋅s/mm) were set to half the numerical value of stiffness (N/mm) to reduce high frequency oscillations at the contact interfaces, which can lead to rapidly evolving ordinary differential equations, and thus computational instability when simulating interface dynamics.

Table 1.

Summary of Young’s Modulus values for various anatomical locations on the residual limb.

Anatomical location	Effective modulus (kPa)	Corresponding interface(s) on the model
Tibial tuberosity 23	105	Anterior
Posterior tibia 54	30	Posterior
Distal tibia 23	60	Distal
Medial proximal tibia 25	56	Medial
Lateral proximal tibia 23	78	Lateral

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. 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). 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. 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 and type of socket liner used[30-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.

Gait simulations with experimental data

The spatial contact model of the residual limb-socket interface was integrated into a gait model with a BKA and a semi-active variable stiffness foot (VSF) prosthesis. This model was previously described in McGeehan et al. (2021a and 2021b).[34,35] Briefly, gait simulations were driven by experimental data from four individuals (Table 2) walking with the VSF configured to “low”, “medium”, and “high” stiffness settings, corresponding to forefoot stiffness values of 10, 19, and 32 N/mm. Forward kinematics simulations were computed for three trials per setting, Briefly, these methods involve using experimental motion capture data (retroreflective markers) collected during in vivo gait trials to control the model’s kinematics in corresponding in silico gait trials. The kinetics are derived using a 25-point ground contact model, which allows contact forces under the prosthesis to evolve according to the model’s kinematics and the mechanical properties of the prosthesis. These methods are described in greater detail in McGeehan et al. (2021a and 2021b).[34,35] Subject 2 did not complete one medium stiffness trial, and Subject 4 did not complete one high stiffness trial. In total, 34 simulations were computed. All simulations were computed in Simscape Multibody using the ode23t solver profile with variable step sizes for numerical integration. Prior written informed consent was provided by all subjects as approved by the Health Sciences Institutional Review Board at the University of Wisconsin-Madison.

Table 2.

Participant characteristics.

Subject	Sex	Age (y)	Height (cm)	Mass (kg)	K-level	Amputation side	Years postamputation
1	Male	34	181	77.3	4	Right	15
2	Male	51	175	111	3	Right	8
3	Male	70	180	83.8	3	Left	14
4	Female	61	163	63.8	3	Right	8
Mean ± SD	–	54 ± 15	175 ± 19.9	84.0 ± 19.9	–	–	11 ± 3.8

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 (fc: 6Hz). Data were indexed from heel strike to toe-off and resampled to 101 data points via cubic spline interpolation. These methods allow for comparison of stance phases of different lengths. Model-derived values were compared to those previously reported in the literature and the effects of stiffness setting on limb-socket dynamics were evaluated using linear mixed effects (LME) regression analysis. It was hypothesized that these effects would be subject-dependent, and as such, an exemplary case study for Subject 1 is presented along with group level data.

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 (metric = ß *stiffness + O ), along with confidence intervals and p-values for each coefficient. The random effect appears as a unique intercept (O ) for each subject. The overall effect of stiffness on each outcome measure (coefficient ß ) was evaluated (α = 0.05). We report the slope, mean intercept, adjusted R2, and p-value for each model. The slope (ß ) represents the sensitivity of the outcome to changes in stiffness, whereas the mean intercept (mean of all O ) allows for characterization of the group average values. The adjusted R2 value represents the strength of the independent variable to explain variations in the dependent variable, excluding the effect of individual intercepts. p-value quantifies the statistical confidence that the slope (sensitivity) is different from zero. These methods were adapted from previous work evaluating the effects of prosthesis design on gait outcomes.[39-41] Analyses were performed using MATLAB 2020b (Mathworks Inc., Natick, MA).

Results

Model performance (group data)

Contact model-derived values for normal pressure and shear stress were dependent upon anatomical location (Table 3) and progression of stance phase (Figures 2 and 3). Peak average normal pressure across stance phase was 70.4 ± 4.28, 75.9 ± 4.44, and 85.0 ± 13.0 kPa for the low, medium, and high stiffness conditions, whereas peak average shear stress values were 25.0 ± 1.52, 26.6 ±1.55, and 29.9 ± 4.61 kPa for the same conditions (Figure 2, Table 4). Increased prosthesis stiffness was associated with increased peak average normal pressure and shear stress (p < 0.05) (Figure 4, Table 4); specific increases were observed at the anterior proximal and anterior distal interfaces (p < 0.05). Stiffness did not significantly affect normal pressure or shear stress at the other interface frames, nor did it affect residual limb pistoning or axial rotation (p > 0.05), though effects were subject-dependent (Table 3).

Table 3.

MLE statistical parameters.

Outcome metric	Slope	Mean intercept	Adjusted R2	p-value
Normal pressure (kPa/BW)
Peak Average	0.009	0.078	0.742	0.003*
Anterior Proximal	0.040	0.213	0.524	0.027*
Anterior Distal	0.018	0.115	0.552	0.022*
Posterior Proximal	0.005	0.046	0.115	0.371
Posterior Distal	0.002	0.028	0.191	0.240
Lateral Proximal	−0.009	0.196	0.092	0.428
Lateral Distal	0.002	0.100	0.052	0.554
Medial Proximal	−0.001	0.077	0.025	0.686
Medial Distal	0.004	0.047	0.130	0.340
Shear stress (kPa/BW)
Peak Average	0.003	0.028	0.678	0.006*
Anterior Proximal	0.015	0.073	0.555	0.021*
Anterior Distal	0.006	0.040	0.552	0.022*
Posterior Proximal	0.002	0.016	0.116	0.371
Posterior Distal	0.001	0.010	0.151	0.302
Lateral Proximal	−0.003	0.069	0.097	0.416
Lateral Distal	<0.001	0.036	0.031	0.651
Medial Proximal	−0.001	0.028	0.025	0.685
Medial Distal	0.001	0.017	0.106	0.393
Kinematics (varying units)
Axial Angle (deg)	−0.078	0.426	0.052	0.557
Pistoning (cm)	0.025	1.446	0.334	0.103

Figure 2.

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.

Figure 3.

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.

Table 4.

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.

Subject (condition)	Normal pressure (kPa)	Shear stress (kPa)	Piston displacement (cm)
1 (low)	75.7 ± 3.71	26.8 ± 1.26	1.48 ± 0.02
1 (medium)	72.7 ± 6.42	25.5 ± 2.22	1.51 ± 0.01
1 (high)	81.54 ± 4.11	28.5 ± 1.44	1.50 ± 0.03
2 (low)	91.2 ± 10.93	32.0 ± 3.84	1.95 ± 0.10
2 (medium)	103 ± 6.47	36.1 ± 2.30	2.07 ± 0.15
2 (high)	119.2 ± 32.6	42.1 ± 11.7	2.09 ± 0.19
3 (low)	50.5 ± 1.08	18.6 ± 0.47	1.22 ± 0.01
3 (medium)	52.1 ± 0.61	18.2 ± 0.19	1.26 ± 0.02
3 (high)	61.0 ± 13.9	21.5 ± 4.89	1.26 ± 0.02
4 (low)	64.1 ± 1.41	22.5 ± 22.5	1.19 ± 0.00
4 (medium)	76.1 ± 4.24	36.6 ± 1.50	1.22 ± 0.03
4 (high)	78.3 ± 1.26	27.4 ± 0.42	1.19 ± 0.03
Group (low)	70.4 ± 4.28	25.0 ± 1.52	1.46 ± 0.03
Group (medium)	78.1 ± 4.44	26.6 ± 1.55	1.52 ± 0.05
Group (high)	85.0 ± 13.0	29.9 ± 4.61	1.51 ± 0.07

Figure 4.

LME regression responses for normal forces (top), shear stresses (middle), and kinematics (bottom).

Table 5.

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.

Subject (condition)	Anterior proximal	Anterior distal	Posterior proximal	Posterior distal	Medial proximal	Medial distal	Lateral proximal	Lateral distal	Distal
1 (low)	259 ± 30.6	132 ± 6.52	41.5 ± 24.9	20.6 ± 4.86	72.9 ± 29.9	47.4 ± 1.61	84.1 ± 2.25	68.6 ± 3.28	0.08 ± 0.06
1 (medium)	238 ± 41.9	123 ± 12.0	53.4 ± 19.4	22.2 ± 4.23	78.1 ± 44.4	43.8 ± 7.41	128 ± 17.6	73.7 ± 4.35	4.91 ± 5.34
1 (high)	292 ± 14.2	141 ± 6.99	39.2 ± 15.2	20.9 ± 2.74	58.3 ± 15.9	40.0 ± 2.19	99.7 ± 10.4	73.3 ± 4.95	4.91 ± 5.34
2 (low)	199 ± 19.2	118 ± 13.9	34.5 ± 6.84	24.5 ± 0.95	56.0 ± 11.3	40.5 ± 2.72	275 ± 32.4	128 ± 11.9	89.5 ± 33.8
2 (medium)	358 ± 44.3	138 ± 12.3	29.9 ± 17.9	25.0 ± 1.03	49.9 ± 11.2	41.0 ± 3.13	283 ± 52.4	132 ± 23.3	145 ± 39.7
2 (high)	320 ± 21.4	196 ± 49.4	27.5 ± 3.37	36.6 ± 11.7	53.4 ± 15.7	71.7 ± 26.4	182 ± 31.8	131 ± 38.4	233 ± 104
3 (low)	131 ± 14.5	81.8 ± 4.62	80.1 ± 5.52	42.5 ± 3.11	76.5 ± 4.08	43.7 ± 3.01	111 ± 16.6	74.5 ± 7.04	0.00 ± 0.00
3 (medium)	151 ± 7.19	87.5 ± 2.46	57.7 ± 8.35	32.4 ± 3.90	69.4 ± 12.5	43.1 ± 5.46	90.1 ± 13.1	55.4 ± 3.07	0.00 ± 0.00
3 (high)	160 ± 3.20	87.7 ± 3.86	88.2 ± 51.4	45.7 ± 26.0	87.0 ± 13.3	54.1 ± 7.38	103 ± 23.6	75.7 ± 25.1	0.00 ± 0.00
4 (low)	188 ± 12.5	37.2 ± 5.50	17.1 ± 0.72	15.7 ± 0.72	38.1 ± 2.53	32.4 ± 2.02	138 ± 6.97	71.8 ± 2.53	0.00 ± 0.00
4 (medium)	256 ± 25.0	122 ± 9.51	26.1 ± 6.90	16.4 ± 0.57	43.3 ± 3.56	37.5 ± 0.39	145 ± 5.49	75.2 ± 1.53	0.00 ± 0.00
4 (high)	271 ± 0.24	128 ± 0.92	42.3 ± 2.83	17.8 ± 0.78	35.9 ± 4.11	30.3 ± 2.20	139 ± 7.78	73.6 ± 3.19	0.00 ± 0.00
Group (low)	194 ± 19.2	107 ± 7.65	43.3 ± 9.50	25.8 ± 2.41	60.7 ± 12.0	41.0 ± 2.34	152 ± 14.5	85.8 ± 6.19	22.4 ± 8.46
Group (medium)	251 ± 29.6	118 ± 9.05	41.8 ± 13.1	24.0 ± 2.44	60.2 ± 17.9	41.3 ± 4.10	162 ± 22.2	84.0 ± 8.05	36.9 ± 10.0
Group (high)	261 ± 9.77	138 ± 15.3	49.3 ± 18.2	30.3 ± 10.3	58.6 ± 12.3	49.0 ± 9.55	131 ± 18.4	88.5 ± 17.9	59.5 ± 27.4

Table 6.

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.

Subject (condition)	Anterior proximal	Anterior distal	Posterior proximal	Posterior distal	Medial proximal	Medial distal	Lateral proximal	Lateral distal	Distal
1 (low)	90.5 ± 10.7	46.3 ± 2.33	15.4 ± 8.13	8.08 ± 2.10	26.9 ± 9.58	16.7 ± 0.76	29.7 ± 1.01	24.7 ± 1.14	0.03 ± 0.02
1 (medium)	83.1 ± 14.7	43.3 ± 4.17	19.0 ± 6.27	8.61 ± 1.78	27.9 ± 15.1	16.5 ± 3.47	44.7 ± 6.25	25.8 ± 1.61	0.98 ± 0.23
1 (high)	102 ± 4.98	49.5 ± 2.44	13.7 ± 5.35	7.51 ± 0.96	20.8 ± 5.05	14.1 ± 0.82	35.4 ± 3.97	26.1 ± 1.74	1.90 ± 1.99
2 (low)	69.9 ± 6.71	41.3 ± 4.89	12.1 ± 2.34	8.62 ± 0.35	19.6 ± 3.95	14.4 ± 1.23	96.5 ± 11.4	45.0 ± 4.32	35.8 ± 12.9
2 (medium)	127 ± 18.1	48.2 ± 4.30	10.5 ± 6.25	8.89 ± 0.23	17.5 ± 3.92	14.3 ± 1.10	99.3 ± 18.5	46.5 ± 8.35	51.7 ± 13.1
2 (high)	120 ± 10.2	68.9 ± 17.5	10.2 ± 1.90	12.9 ±4.18	19.4 ± 5.94	25.7 ± 9.65	64.7 ± 11.5	46.7 ± 13.9	81.9 ± 36.4
3 (low)	45.8 ± 5.10	28.6 ± 1.63	28.3 ± 1.97	15.4 ± 1.10	26.8 ± 1.43	15.3 ± 1.06	40.1 ± 5.32	27.3 ± 2.38	0.00 ± 0.00
3 (medium)	52.8 ± 2.51	30.6 ± 0.87	20.4 ± 2.91	11.6 ± 1.42	24.7 ± 4.76	15.5 ± 2.17	31.5 ± 4.58	19.5 ± 1.01	0.00 ± 0.00
3 (high)	55.9 ± 1.13	30.7 ± 1.34	31.0 ± 18.0	16.3 ± 8.99	30.4 ± 4.66	19.0 ± 2.69	36.3 ± 8.41	26.8 ± 8.95	0.00 ± 0.00
4 (low)	65.7 ± 4.37	34.0 ± 1.93	6.01 ± 0.27	5.54 ± 0.31	13.8 ± 0.69	12.2 ± 1.21	48.3 ± 2.43	25.1 ± 0.85	0.00 ± 0.00
4 (medium)	89.8 ± 8.74	42.8 ± 3.33	9.11 ± 2.42	5.73 ± 0.19	15.3 ± 1.37	13.5 ± 0.38	50.9 ± 1.93	16.3 ± 0.54	0.00 ± 0.00
4 (high)	95.0 ± 0.08	44.9 ± 0.32	15.9 ± 0.08	6.71 ± 0.73	13.7 ± 2.53	11.0 ± 0.37	48.5 ± 2.71	25.8 ± 1.16	0.00 ± 0.00
Group (low)	68.0 ± 6.72	37.6 ± 2.70	15.4 ± 3.18	9.41 ± 0.97	21.8 ± 3.91	14.7 ± 1.06	53.7 ± 5.03	30.5 ± 2.17	8.95 ± 3.23
Group (medium)	88.2 ± 11.0	41.2 ± 3.17	14.7 ± 4.46	8.70 ± 0.91	21.4 ± 6.28	14.9 ± 1.78	56.6 ± 7.82	29.5 ± 2.88	13.2 ± 3.33
Group (high)	93.3 ± 4.09	48.5 ± 5.40	17.7 ± 6.33	10.9 ± 3.72	21.1 ± 4.54	17.5 ± 3.38	46.2 ± 6.59	31.3 ± 6.43	21.0 ± 9.58

Spatiotemporal patterns for pressure and shear stress distribution were variable between participants, but show similar variability across stiffness conditions (Figure 3, Tables 5 and 6). On average, participants displayed predominantly anterior pressure and shear distributions early in stance phase, but trended toward a more even distribution later in stance phase. Pressure trended slightly toward the lateral and proximal aspects of the residual limb thoughout stance phase. High variability was observed among participants for the AP and ML distributions throughout stance phase.

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.

Figure 5.

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.

Figure 6.

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

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 and increased propulsion. 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) and Laszczak et al. (2016) 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 experiments[8,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). 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 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. 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, and gait conditions.[29,32] Data regarding the prosthetic socket componentry used by participants in this study were not available.

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. 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. 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. Increased external rotation has been associated with high rates of medial knee osteoarthritis, which has been documented for both the amputated and contralateral limb of lower limb prosthesis users. 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 and biplane fluoroscopy 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

: normal force (N)

: contact stiffness (N/mm)

: penetration depth (mm)

: contact damping coefficient (N⋅s/mm)

: frictional force (N)

: coefficient of friction (unitless)

: velocity at point of contact (mm/s)

: velocity threshold (mm/s)

: Stiffness in the normal plane (N/mm)

: Young’s modulus of the tissue (N/mm2)

: Area of the contact point (mm2)

: Width of the residual limb (mm)