Relationship between Joint Stiffness, Limb Stiffness and Whole–Body Center of Mass Mechanical Work across Running Speeds

: The lower–extremity system acts like a spring in the running stance phase. Vertical stiffness ( K vert ) and leg stiffness ( K leg ) reﬂect the whole–body center of mass (COM) and leg–spring system loading and response in running, while joint stiffness ( K joint ) represents joint–level dynamic loading and response. This study aimed to investigate whether K joint is associated with K vert and K leg across different running speeds. Twenty healthy subjects were recruited into a treadmill running study (1.8 to 3.8 m/s, with 0.4 m/s intervals). We found that K joint accounted for 38.4% of the variance in K vert ( p = 0.046) and 42.4% of the variance in K leg ( p = 0.028) at 1.8 m/s; K joint also accounted for 49.8% of the variance in K vert ( p = 0.014) and 79.3% of the variance in K leg ( p < 0.0001) at 2.2 m/s. K knee had the strongest unique association with K vert and K leg at 1.8 and 2.2 m/s. K joint was associated with K leg at a wider range of speeds. These ﬁndings built a connection between joint stiffness and limb stiffness within a certain range of running speeds. K knee may need to be considered as an important factor in future limb stiffness and general running performance enhancement.


Introduction
The lower extremity is compliant during the stance phase of running [1], with the joints going through a flexion and then an extension movement [1]. These motions suggest that in response to external force, the lower extremity musculoskeletal system acts like a spring, absorbing energy in the first half of the stance and returning a portion of elastic energy in the second half of the stance [2][3][4]. This results in the whole-body center of mass (COM) position reaching its minimum height at mid-stance'; the movement trajectory is similar to a bouncing ball [1,3]. Using this analogy, a simplified spring-mass model has been proposed, and is widely used in the analysis of human running gait [1,[5][6][7][8][9].
The loading and unloading characteristics of the leg spring system under external moment and force in the running stance phase can be regarded as stiffness patterns. Vertical stiffness (K vert ), leg stiffness (K leg ) and joint stiffness (K joint ) can be directly calculated from running activities [6]. Moreover, K vert and K leg can be calculated via the spring-mass model mentioned previously. K vert is the peak ground reaction force (GRF) divided by the vertical COM displacement, and it reflects COM vertical movement and oscillation characteristics in the stance phase [6,10,11]; it has been reported to increase with running speeds [6,[12][13][14]. This may be attributed to an increase in the peak vertical GRF whilst COM displacement decreases when running speeds increase [6]. K leg is the peak GRF divided by the maximum leg vertical displacement during ground contact [5][6][7], and K leg

Experimental Protocol and Data Collection
We measured the participants' body mass, height and leg length (L 0 ) before the running test. Leg length (L 0 ) was measured as the vertical distance from the greater trochanter to the floor during static standing [9]. Then, fifty-five retro-reflective markers were placed on the skin surface of the participants, based on a previously published wholebody marker set [36]. The participants were asked to run on a force-instrumented treadmill (Bertec, Inc., Columbus, OH, USA) at six different speeds, from 1.8 to 3.8 m/s (0.4 m/s intervals), for 75 s per stage. Data were extracted from the middle strides (20 strides on average) of each stage. Segmental kinematic data were collected at 120 Hz using an 8-camera motion capture system (Motion Analysis Corp., Santa Rosa, CA, USA). Ground reaction force data were collected at 1200 Hz using the force-instrumented treadmill. Kinematic and kinetic data were filtered with a low-pass fourth-order Butterworth filter at 6 Hz and 50 Hz, respectively.

Data Analysis
The whole-body COM position (X com ) was calculated from the weighted sum of a 15-segment (head, trunk, pelvis, upper arms, lower arms, hands, thighs, shanks, and feet) full-body model [37] for each subject in Visual 3D (C-Motion, Inc., Germantown, MD, USA). Specifically, it was calculated as follows: where n is the number of the segment, m seg i is each individual segment's mass, X seg com i is each individual segment's center of mass coordinate, and m b is the whole-body mass.
The spring-mass model vertical stiffness (K vert ) was calculated from the peak vertical ground reaction force (vGRF peak ) divided by the vertical displacement of the COM from ground contact until mid-stance (∆y) ( Figure 1) [2,[5][6][7][8][9], expressed as iomechanics 2022, 3, FOR PEER REVIEW stiffness ( ) was calculated as the change in the sagittal plane joint mom divided by the sagittal plane joint angular displacement (∆ ) in the first h contact, based on the anterior-posterior ground reaction force value [22,39], Figure 1. Schematic representative of a spring-mass model in the running stance ph consists of a point mass (COM) equivalent to the body mass and the leg as a massles The leg spring is compressed, and reaches maximum compression (∆ ) at mid-sta displacement in the vertical direction is denoted as ∆ . The half swept angle of th denoted as . IC: initial-contact. MS: mid-stance. TO: toe-off. Point O is the ground c The COM gravitational potential energy ( ) was calculated as the p body mass ( ), gravitational constant ( = 9.81 m/s 2 ), and instantaneous (ℎ ) [28], expressed as The leg spring is compressed, and reaches maximum compression (∆L) at mid-stance. The COM displacement in the vertical direction is denoted as ∆y. The half swept angle of the leg spring is denoted as θ. IC: initial-contact. MS: mid-stance. TO: toe-off. Point O is the ground contact location.
The half swept angle (θ) was defined as the angle between the leg-spring at ground contact and mid-stance (Figure 1), and it was calculated from the running speed (µ), ground contact time (t c ) and initial leg length (L 0 ) [2,8], expressed as The leg-spring maximum displacement (∆L) can be calculated via the expression of changes in the vertical COM displacement (∆y), half swept angle (θ) and initial leg length (L 0 ) ( Figure 1) [2,5,7,9], expressed as Leg stiffness (K leg ) was calculated as the peak vertical ground reaction force (vGRF peak ) divided by the leg-spring maximum displacement (∆L) (Figure 1) [2,[5][6][7][8][9], expressed as Lower-extremity joint moments were calculated using a standard inverse dynamics model [38] coded in Visual 3D. In this study, each joint neutral position was defined as the zero-degree reference angle in the sagittal plane, joint extension was defined as positive, and flexion was defined as negative, in comparison with the neutral position. Joint stiffness (K joint ) was calculated as the change in the sagittal plane joint moment (∆M joint ) divided by the sagittal plane joint angular displacement (∆θ joint ) in the first half of ground contact, based on the anterior-posterior ground reaction force value [22,39], expressed as The COM gravitational potential energy (E pot ) was calculated as the product of the body mass (m b ), gravitational constant (g = 9.81 m/s 2 ), and instantaneous COM height (h i ) [28], expressed as The value of the COM kinetic energy (E kin ) was calculated from the sum of E kin in both the horizontal and vertical direction [28], expressed as where v h and v v are the COM velocity in the horizontal and vertical directions, respectively. We also calculated the COM instantaneous power in the horizontal (P comh ) and vertical (P comv ) directions, and the sagittal plane (P coms ), based on the definition of a previous study [28], expressed as where a h and a v are the COM acceleration in the horizontal and vertical directions, respectively. Moreover, the COM positive (W + com ) and negative external mechanical work (W − com ) in the horizontal and vertical directions, as well as in sagittal plane were calculated as the instantaneous positive (P + com ) or negative power (P − com ) in each direction or plane integrated over time, respectively [28].
Ground reaction force (GRF) and virtual leg length (instantaneous leg length/L 0 ) force-length relationships were plotted for the average of the twenty participants, for further interpretation (Figure 2). The curve slope was estimated via the tangent function between the curve ascending phase starting point and the ascending phase ending point horizontal and vertical axis coordinate values, respectively ( Figure 2). The group mean COM potential energy (E pot ), kinetic energy (E kin ) and sagittal plane COM instantaneous power (P coms ) were plotted from three representative speeds (1.8, 2.6, 3.8 m/s) as well (Figure 3a,b). All of the graphs were plotted in the MATLAB program (R2018a, Mathworks, Natick, MA, USA). tive work ( ) and negative work ( ), and COM peak positive ( ) and negative power ( ) in the sagittal plane were examined using a one-way ANOVA to compare among six speeds. The initial alpha level was set to 0.05. When the main effect was detected, Bonferroni adjustments were used for pairwise comparison, such that the alpha level was divided by the number of comparisons (adjusted α = 0.0033 for all pairwise comparisons in this study). Additionally, multiple linear regression analysis was conducted to develop models to build potential associations between (ankle, knee and hip joint stiffness), and within each running speed. Lastly, simple linear regression analysis was used to examine the relationships between the sagittal plane COM positive work ( ) and and across the speeds. All of the statistical analyses were performed using SPSS (V22.0, IBM, Armonk, NY, USA). All of the outcome variables were calculated and averaged from both limbs, and were averaged across three selected gait cycles among the 20 strides for each stage. In order to make better comparisons with previous studies, only K joint was normalized to body weight. The vertical stiffness (K vert ), leg stiffness (K leg ), joint stiffness (K joint ), COM positive work (W + com ) and negative work (W − com ), and COM peak positive (P + coms ) and negative power (P − coms ) in the sagittal plane were examined using a one-way ANOVA to compare among six speeds. The initial alpha level was set to 0.05. When the main effect was detected, Bonferroni adjustments were used for pairwise comparison, such that the alpha level was divided by the number of comparisons (adjusted α = 0.0033 for all pairwise comparisons in this study). Additionally, multiple linear regression analysis was conducted to develop models to build potential associations between K joint (ankle, knee and hip joint stiffness), K vert and K leg within each running speed. Lastly, simple linear regression analysis was used to examine the relationships between the sagittal plane COM positive work (W + coms ) and K vert and K leg across the speeds. All of the statistical analyses were performed using SPSS (V22.0, IBM, Armonk, NY, USA).

Stiffness
The comparison of among all of the running speeds was not significant (p = 0.413). Speed's main effect for was significant (p < 0.0001), therefore, pairwise comparison was conducted (Table 1): at 1.8 m/s was significantly lower than all speeds between 2.6 and 3.8 m/s (p < 0.0001), at 2.2 m/s was lower than all speeds between 3.0 and 3.8 m/s (p ≤ 0.0001), at 2.6 m/s was lower than at 3.4 m/s (p = 0.001) and 3.8 m/s (p = 0.0002), and at 3.0 m/s was lower than at 3.8 m/s (p = 0.0032). For comparison, speed's main effect was significant in (p < 0.0001), and pairwise

Stiffness
The comparison of K leg among all of the running speeds was not significant (p = 0.413). Speed's main effect for K vert was significant (p < 0.0001), therefore, pairwise comparison was conducted (Table 1): K vert at 1.8 m/s was significantly lower than all speeds between 2.6 and 3.8 m/s (p < 0.0001), K vert at 2.2 m/s was lower than all speeds between 3.0 and 3.8 m/s (p ≤ 0.0001), K vert at 2.6 m/s was lower than at 3.4 m/s (p = 0.001) and 3.8 m/s (p = 0.0002), and K vert at 3.0 m/s was lower than at 3.8 m/s (p = 0.0032). For K joint comparison, speed's main effect was significant in K knee (p < 0.0001), and pairwise comparison was conducted: K knee at 1.8 m/s was lower than that at 3.0 m/s (p = 0.002) and 3.8 m/s (p = 0.001), and K knee at 2.2 m/s was lower than that at 3.0 m/s (p = 0.001) and 3.8 m/s (p = 0.003).

Mechanical Work and Power
Speed's main effects were significant in both W + coms (p < 0.0001) and W − coms (p = 0.002); therefore, pairwise comparison was conducted ( Table 2)   Additionally, speed's main effects were significant in both P + coms and P − coms (p < 0.001), and pairwise comparison was conducted ( Table 2): P + coms at 1.8 m/s was lower than at all speeds between 2.2 and 3.8 m/s (p < 0.001); P + coms at 2.2 m/s was lower than at 2.6, 3.4 and 3.8 m/s, respectively (p < 0.001); P + coms at 2.6 m/s was lower than at 3.4 and 3.8 m/s, respectively (p ≤ 0.001). For P − coms at 1.8 m/s, it was lower than at 2.2, 3.0, 3.4 and 3.8 m/s, respectively (p ≤ 0.002); additionally, P − coms at 2.2 m/s was lower than at 3.0 and 3.4 m/s, respectively (p < 0.001).

Multiple and Simple Linear Regression
The results from the multiple linear regression analysis showed that K joint was associated with K vert at 1.8 m/s and 2.2 m/s (Table 3). At 1.8 m/s, the model accounted for 38.4% of the variance in K vert (R 2 = 0.384, p = 0.046), and K knee had the strongest unique association with K vert at this speed (β = 0.509, p = 0.022). At 2.2 m/s, the model accounted for 49.8% of the variance in K vert (R 2 = 0.498, p = 0.014), and K knee again had the strongest unique association with K vert at this speed (β = 0.553, p = 0.011). Additionally, the multiple linear regression analysis revealed that K joint was associated with K leg among most speeds, except at 3.0 m/s and 3.8 m/s (Table 3). At 1.8 m/s, the model accounted for 42.4% of the variance in K leg (R 2 = 0.424, p = 0.028), and K knee had the strongest unique association with K leg (β = 0.532, p = 0.014). At 2.2 m/s, the model accounted for 79.3% of the variance in K leg (R 2 = 0.793, p < 0.0001). For this speed, however, K knee (β = 0.553, p = 0.0004) and K hip (β = 0.526, p = 0.001) both had a strong unique association with K leg . At 2.6 m/s, the model accounted for 39.9% of the variance in K leg (R 2 = 0.399, p = 0.039), and K knee had a unique association with K leg (β = 0.456, p = 0.04). At 3.4 m/s, the model accounted for 47.4% of the variance in K leg (R 2 = 0.474, p = 0.026), and K hip had a strong unique association with K leg (β = 0.721, p = 0.009).
Simple linear regression analysis showed that K leg was not associated with W + coms across the speeds (R 2 = 0.133, p = 0.477). However, K vert was positively associated with W + coms across the speeds (R 2 = 0.902, r = 0.95, p = 0.004) ( Table 4).

Interpretation of Graph Patterns
Based on the stance phase ground reaction force and the virtual leg length relationship for three representative speeds, we found that the slope of the curve increased as running speeds increased (estimated curve slope value: 30 at 1.8 m/s, 38 at 2.6 m/s, 56 at 3.8 m/s), and the virtual leg length magnitude at both initial ground contact and in the take-off phase tended to decrease (virtual leg length value: 0.97 at 1.8 m/s, 0.96 at 2.6 m/s, 0.94 at 3.8 m/s; Figure 2). The COM E pot slightly decreased as the running speeds increased (initial contact phase: 670-650 J from 1.8 to 3.8 m/s; mid-stance phase: 630-620 J; take-off phase: 680-660 J), while the magnitude of E kin increased dramatically when speeds increased (initial contact phase: 125-510 J from 1.8 to 3.8 m/s; mid-stance phase: 105-470 J; take-off phase: 120-513 J; Figure 3a).

Discussion
The primary goal of this study was to investigate whether K ankle , K knee and K hip were associated with K vert and K leg using multiple linear regression models for each running speed. Additionally, we investigated whether W + coms was associated with K vert or K leg across the running speeds. The initial hypothesis that K joint would be associated with K vert and K leg was supported. The hypothesis that W + coms was associated with K vert and K leg was partially supported.
K joint was associated with both K vert and K leg in the multiple linear regression models at slow speeds (1.8 and 2.2 m/s) (Table 3). Furthermore, K knee had a significant unique association with K vert and K leg at these speeds. However, K joint was not associated with K vert among speeds from 2.6 to 3.8 m/s. One reason may be that K vert tended to increase as running speeds increased, due to the increased vertical GRF and decreased COM displacement [6]. However, the change of running speeds had mixed effects on K joint . Specifically, when the running speed increased, K knee tended to increase, while K ankle and K hip remained almost constant, and they did not have a linear relationship with the change of running speeds (Table 1). Another reason may be that K vert is more related to whole-body COM bouncing and oscillation patterns [6,10,11], and K joint likely has a closer relationship with leg-spring stiffness than with COM oscillation characteristics.
For multiple linear regression analysis between K joint and K leg , the values of K knee and K hip were more associated with K leg (Table 3). Both K knee and K hip were associated with K leg at 2.2 m/s. Interestingly, K ankle did not have much association with K leg across all of the running speeds in this study. However, K knee was associated with K leg among most speeds. This may be attributed to the idea that the human leg is a system comprised of multiple springs, and the sub-springs can be coordinated with each other during ground contact in running. Under similar loading conditions, the spring with the smallest stiffness will undergo the largest displacement, and this would have the most influence on the overall leg-spring system stiffness [23]. In this study, K knee tended to be lower than K ankle and K hip across all the running speeds (Table 1). Besides having more association with K leg among speeds, knee joint flexion (indicating relatively lower stiffness) could also be beneficial for elastic energy storage in the first half of the stance phase and the subsequent energy return in the second half of the stance [17,22]. The joint level stiffness is influenced by both tendon stiffness and active control of the knee muscle activation [22].
In the simple linear regression analysis, K vert and W + coms had a strong positive association across the running speeds. This may be due to the observation that both K vert and W + coms tended to increase with the running speed. In response to greater GRF impacts, decreasing COM sagittal plane displacement and oscillation may reduce the amount of mechanical energy being absorbed via the spring-mass system in the first half of the stance phase; this would allow for more positive mechanical work to be generated through the whole-body COM.
The other goal of the study was to examine whether a change of running speeds would affect K vert , K leg , W com and P com . The initial hypothesis was partially supported. The results showed that K vert increased with the running speeds, while K leg remained unchanged from 1.8 to 3.8 m/s. These findings agreed with previous findings [5][6][7][12][13][14][15][16].
Changes of speed influenced both positive and negative W com in the sagittal plane and in the horizontal direction, as well as the sagittal plane peak positive and negative P com (Table 2). However, changes of speed did not have significant effects on either positive or negative W com in the vertical direction. This finding can be explained by the COM E pot and E kin curve patterns (Figure 3a). Among the three representative running speeds, both the maximum and minimum E pot values slightly decreased around 3%, from 1.8 to 3.8 m/s, while the magnitude dramatically increased around 124% for E kin as the running speeds increased (Figure 3a). This indicates that changes of running speed had more effect on E kin than E pot . Additionally, there was a greater change of COM velocity in the horizontal direction than in the vertical direction in this speed increment running protocol, and E kin was affected more by a speed change in the horizontal direction than in the vertical direction. Furthermore, GRF increased in both the vertical and horizontal directions as speeds increased, indicating that COM energy absorption was greater in the first half of the stance, and that higher speeds required more positive mechanical work generated on the COM to assist the body to move forward in the following propulsive phase. This helps explain why W + comh and W − comh increased as speeds increased. We also investigated the vertical GRF and virtual leg length relationship in three representative speeds (Figure 2). The curve consisted of an ascending and a descending phase. The ascending phase represents the loading period, and the descending phase represents the unloading period. Within the ascending phase, the "yielding" pattern became more obvious as speeds increased. Additionally, the virtual leg length at initial contact decreased as speed increased, indicating that the leg-spring compressed more with increased speed. This would be beneficial for energy absorption as external impact forces increase, and it could also be beneficial for the reduction of COM height and E pot as speed increases. Moreover, the magnitude of the virtual leg length change tended to decrease as the speed increased ( Figure 2). This indicates that the leg-spring became stiffer as the running speed increased.
One limitation of this study is that the leg spring was assumed not to be compressed at initial ground contact in the spring-mass model. As speed increased, the initial leg length was less than the static standing leg length (L 0 ), which was used in the K leg calculation. This likely affected the K leg results at relatively higher speeds. Furthermore, the model used to calculate the leg-spring displacement [5] underestimated the real leg-spring displacement, and this may also affect the K leg values [21,40]. However, we checked the K leg and K vert results with a sine-wave model [14], and they both derived similar results. Additionally, a treadmill running protocol was used in this study, with controlled locomotion speeds, and thus some individual variations may have been constrained. Another limitation is that we investigated a slow-to-medium range of running speeds. Whether the COM dynamic patterns would be different in a wider range of speeds requires further investigation.
Future studies should compare the accuracy of different models in the estimation of COM dynamic patterns during locomotion. In this study, we calculated the COM instantaneous mechanical power from kinematic variables of COM movement (COM velocity and acceleration). The method was previously shown to be reliable in the estimation of COM displacement compared with the method derived from GRF [28]. Other studies have used the dot product of GRF and COM velocity to estimate COM external mechanical power, with the COM velocity being derived from the integration of GRF in these studies [29,30,41]. Further comparison between these two methods in both walking and running across different speeds is needed.

Conclusions
In conclusion, when running at slow-to-medium speeds, whole-body COM positive and negative mechanical work tended to increase in the sagittal plane and in the horizontal direction. Lower-extremity joint stiffness was associated with both leg stiffness and vertical stiffness using multiple linear regression models at 1.8 m/s and 2.2 m/s. Joint stiffness was associated with leg stiffness at a wider range of running speeds compared with vertical stiffness. The knee joint was more associated with vertical stiffness and leg stiffness. Sagittal-plane COM positive work and vertical stiffness had a strong positive association when running speed increased. These findings suggest that leg-spring system stiffness was associated with subsystem joint-level stiffness characteristics. Lastly, whole-body COM mechanical work had a strong positive association with COM oscillation patterns in the stance phase of running across different speeds. These findings build a connection between joint stiffness and limb stiffness, as well as whole-body COM mechanical work and oscillation patterns across different running speeds. The outcomes of this study may be applicable to improved designs for running-specific prostheses, and to enhancing our understanding of general running performance based on joint and limb stiffness.  Informed Consent Statement: Informed consent was obtained from all of the subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.