To develop an anatomically accurate multi-segment foot model, computed tomography (CT) data of the foot of an adult male (age: 42 years, weight: 72 kg, height: 172 cm) were obtained (Ito et al., 2022). Three-dimensional surface models of the foot surface and skeleton were constructed (Figure 1A) using segmentation software (Analyze 9.0, Biomedical Imaging Resource, Mayo Clinic, Rochester, MN, United States). The human foot was represented as a chain comprising three segments (phalanx, forefoot, and hindfoot segments) based on the marker positions attached to anatomical landmarks that are widely used in kinematic studies (Leardini et al., 2007) (Table 1; Figure 1A). The joints between the tibia and hindfoot (ankle joint), hindfoot and forefoot (midtarsal joint), and forefoot and phalanx segments (MTP joint) were defined as the midpoints between the lateral and medial malleoli (ANKL and ANKM), the navicular tubercle (TN) and the fifth metatarsal base (VMB), and the heads of the first and fifth metatarsals (FMH and VMH), respectively. The coordinate systems of the three segments were defined as follows: the y and z axes of the hindfoot segment were defined as the normal vector of the plane defined by the superior (CA) and inferior (HE) points of the calcaneal tuberosity and the midtarsal joint, and the axis connecting the CA and HE, respectively; those of the forefoot segment were defined as the axis connecting the TN and VMB, and the normal vector of the plane defined by the TN, VMB, and MTP joints, respectively; those of the phalanx segment were defined as the axis connecting the FMH and VMH, and the normal vector of the plane defined by FMH, VMH, and the head of the first proximal phalanx (PM), respectively.

We represented the PA using five linear springs (PA1−5, from medial to lateral) connecting the hindfoot and phalanx segments via a point on the metatarsal plantar surface (Figure 1B). The positions of the origin, via point, and insertion of the PA were obtained from the CT data of the foot (Ito et al., 2022). The origin of the PA was defined as the plantar surface of the calcaneal tuberosity. The via points were defined as the plantar surfaces of the sesamoid of the first metatarsal and the second to fifth metatarsal heads. The insertion was defined as the plantar surface of the base of the first to fifth proximal phalanges. The via points were necessary to model the wrapping paths of PA slips around the metatarsal heads, analogous to the pully model in Caravaggi et al. (2009). To create a subject-specific model, the positions of the origin, via point, and insertion extracted from the CT data were described in the segment coordinate systems of the hindfoot, forefoot, and phalanx segments, the origins of which were HE, VMB, and VMH, respectively, and scaled by the length of the line connecting ANKL and VMH.

The length of the ith PA was calculated as the sum of the lengths from the origin to the via point and from the via point to insertion as: L i = | p i n s _ i − p v i a _ i | + | p v i a _ i − p o r i _ i |   (1)

The PA tension force ( f P A i ) can be calculated as: f P A i = max ( 0 ,   k i ( L i − L 0 i ) ) (2) where k i is the spring constant and L 0 i   is the natural length of the ith PA. The natural length of the PA was estimated to be 0.98 times the PA length during quiet standing ( L q s _ i ), since the foot (hence, the PA) should be slightly stretched from its natural state due to flattening of the foot arch during quiet standing. The spring constants of the PAs were calculated based on studies reporting that the strain (ε) of the PA was approximately 0.07 (the range of 0.03–0.12 based on Caravaggi et al. (2009, 2010) and Gefen (2003), and the force generated by the PA was approximately 1.5 times the body weight (Caravaggi et al., 2009) at the time of push-off during walking. The spring constant can be calculated as: k i = ( 9.8 ) ( 1.5 B W ) 5 ( 0.07 ) ( 0.98 L q s _ i ) (3) where B W is the individual-specific body weight. Therefore, the spring constant varies by individual depending on the body weight, but this is reasonable since heavier individuals have thicker (Pascual Huerta and Alarcón García, 2007) and hence harder PAs (the spring constant is proportional to the cross-sectional area of the PA if the material property is the same; see Ito et al., 2022). The tension force vectors ( f i n s P A i , f o r i P A i , f v i a P A i ) acting on the insertion, origin, and via point of the PA can be defined as: f i n s P A i = f P A i p v i a _ i − p i n s _ i | p v i a _ i − p i n s _ i |    (4) f o r i P A i = f P A i p v i a _ i − p o r i _ i | p v i a _ i − p o r i _ i |   (5) f v i a P A i = ( − f i n s P A i ) + ( − f o r i P A i ) (6)

Based on the free body diagrams in Figure 2, the Newton−Euler equations of the motion of the multi-segment foot model consisting of the phalanx, forefoot, and hindfoot segments can be written as: m P r ¨ P G = f P J + f P E 1 + ∑ i = 1 5 f i n s P A i + m P g (7) M P { ω P × ( I P P ω ) + I P P ω ˙ } = n P J + ( r P J − r P G ) × f P J + ( r P E 1 − r P G ) × f P E 1 + n P E 1 + ∑ i = 1 5 ( r i n s P A i − r P G ) × f i n s P A i (8) m F r ¨ F G = f F J + ( − f P J ) + f F E 1 + f F E 2 + ∑ i = 1 5 f v i a P A i + m F g (9) M F { ω F × ( I F F ω ) F + I F ω ˙ } = n F J + ( − n P J ) + ( r F J − r F G ) × f F J + ( r P J − r F G ) × ( − f P J ) + ( r F E 1 − r F G ) × f F E 1 + n F E 1 + ( r F E 2 − r F G ) × f F E 2 + n F E 2 + ∑ i = 1 5 ( r v i a   P A i − r F G ) × f v i a   P A i (10) m H r ¨ H G = f H J + ( − f F J ) + f H E 2 + ∑ i = 1 5 f o r i P A i + m H g (11) M H { ω H × ( I H H ω ) + I H H ω ˙ } = n H J + ( − n F J ) + ( r H J − r H G ) × f H J + ( r F J − r H G ) × ( − f F J ) + ( r H E 2 − r H G ) × f H E 2 + n H E 2 + ∑ i = 1 5 ( r o r i P A i − r H G ) × f o r i P A i (12)

The variables in the equations are presented in Table 2. The left superscripts P, F and H denote the coordinate systems in which the corresponding vector or matrix was represented.

To calculate the inertial properties (i.e., mass, inertial tensor around the center of mass (COM), and position of the center of mass in each segment), the foot surface was divided into segments by planes passing through ANKL and ANKM, TN and VMB, and FMH and VMH (Figure 1C). Therefore, the hindfoot, forefoot, and phalanx segments correspond to the calcaneus, talus, cuboid, and navicular; cuneiforms and metatarsals; and five phalanges, respectively. The inertial parameters of each segment were calculated using a computer-aided design software (Autodesk Inventor Professional, Autodesk, United States), assuming a homogeneous segment composition and a density of 1.1 g/cm³ (Winter, 1990). To create an individual-specific foot model, the relative mass of the foot segments and the relative inertial tensor around the segment COM were computed based on the inertial parameters from the CT data. Relative foot segment mass was defined as the mass of the foot segment as a percentage of the total foot mass. The relative inertial tensor around the segment COM was defined as the inertial tensor around the segment COM normalized by 5/3 power of each segment mass. Relative COM position was defined as the location of the COM, expressed as a percentage of the segment length from the proximal point of the segment, assuming that the COM is located on the line connecting the centers of the proximal and distal joints.

Ten males (age: 23.9 ± 3.0 years, height: 171.8 ± 5.1 cm, weight: 62.8 ±8.2 kg) without any deformity of the foot and lower extremity and with no history of orthopedic, neurological, and musculoskeletal disorders that are likely to affect gait were recruited for data collection. The number of participants was determined by referring to previous studies (Leardini et al., 2007; Portinaro et al., 2014). All participants provided written informed consent following a detailed explanation of the study’s purpose and the risks involved. The experimental procedures used in this study complied with the Declaration of Helsinki and were approved by the Ethics Committee on Human Experimentation at Saitama Prefectural University (No. 29508).

Infrared-reflecting markers (diameters: 9.5 and 14 mm) were attached to 65 landmarks on the foot and the whole body, according to Table 1 and Figure 1A, and the Plug-in-Gait Full-body Ai model (Davis et al., 1991; Vicon Motion Systems Limited, 2016) (Supplementary Figure S1; Supplementary Table S1), respectively. The participants walked at a self-selected speed along a walkway with four force plates (two per side) placed in series close together. The three attempts in which the fore- and hindfoot contacted the front and rear force plates, respectively (Supplementary Figure S2), were selected for the inverse dynamic analysis.

Marker trajectories were collected using the Vicon Nexus 2.10.2, a three-dimensional motion analysis system (Vicon, Oxford, United Kingdom) with 20 infrared cameras at 100 Hz. Ground reaction force was collected from four force plates (Kistler Instrumente AG, Winterthur, Switzerland) at 1,000 Hz. All data were synchronized using the Vicon Workstation v4.5 software and saved for offline analysis.

Marker trajectories and ground reaction forces were filtered using a zero-phased lag and fourth-order Butterworth filter with a cutoff frequency of 10 Hz. The 3D rotation angles of the ankle joint were described by the y-x-z Euler angle, and those of the midtarsal and MTP joints were described by z-x-y Euler angles. The rotational angles around the y-, x-, z-axes represent plantarflexion–dorsiflexion, inversion–eversion and adduction−abduction, respectively. The MTP joint in the present model is regarded as a hinge joint allowing only plantarflexion–dorsiflexion; hence the rotational angles around the other two axes are not presented.

For inverse dynamics, the total mass of each individual’s foot was estimated to be 0.0145 times the body weight of the individual (Winter, 1990). The mass, inertial tensor, and COM position of each foot segment were then calculated based on the relative mass, inertial tensor, and COM position obtained from the CT scan of the foot used to create the present multi-segment foot model (See Section 2.1). The ground reaction forces were applied to the foot segments based on the force plate data. The positional relationship between the center of pressure and the markers specified the segments on which the two force vectors were applied. The tension forces due to the PAs were applied to the insertions, via points, and origins of the PAs (Figure 1D). The joint moments were calculated by solving the equations for motions (Eqs 7–12) consecutively from the distal phalanx segment to the proximal hindfoot segment. The calculated joint moment vectors were transformed to the corresponding proximal segment coordinate systems to match them with the joint angles. The joint powers were calculated by multiplying the joint angular velocities and joint moments. The joint power curves were integrated into the positive and negative portions during one gait cycle, which represented energy generation and energy absorption, respectively. All data were analyzed using MATLAB 2018a (MathWorks, Natick, MA, United States).

In this study, the inverse dynamic analysis was performed with and without the PA to investigate how inclusion of the PA in the model alters the estimated joint moments and work. For this, the maximum and minimum values of the calculated joint moments during one gait cycle were compared between the two conditions. In addition, the PA contribution rate (%PA contribution) was calculated as the percentage of the joint moment generated by the PA with respect to the net joint moment, to quantify the contribution of the PA on joint moment generation at the time of push-off during walking. The changes in the positive and negative work generated or absorbed during walking were also quantified and compared.

We estimated the relative mass of the segment and the inertia tensor of each segment from the CT data of a single male participant. Further, due to the uncertainty regarding the stiffness and the natural length of each PA spring, a sensitivity analysis was conducted to evaluate how the estimated joint moments were altered by the changes in these parameters within a reasonable range. For the mass and inertia of the segment, we made 0.5- and 1.5-fold changes in the mass and inertial parameters of the foot model to assess how changes in the mass and inertial parameters affected the joint moments, which was calculated using inverse dynamics. We made the 0.5- and 1.5-fold changes in the mass and inertial parameters because the standard deviations of the mass and inertial parameters were reportedly about a quarter of the respective mean values (Zatsiorsky, 2002) and approximately 95% of the population lies within two standard deviations. For the stiffness, we decreased and increased the natural PA length (default L 0 i = 0.98 L q s _ i ) by 0.02 L q s _ i , and the PA strain at toe-off (default ε = 0.07) by 0.03 to assess how changes in the stiffness altered the results. The ranges were determined by referring to previous studies (Gefen, 2003; Caravaggi et al., 2009).

To test for possible statistical differences in the joint moment and power profiles between the models with and without PA in a continuous manner, we performed a one-dimensional Statistical Parametric Mapping (SPM) paired t-test (Pataky, 2012). To compare the statistical differences in the estimated peak joint moments, joint works, and %PA contributions, we performed a paired t-test if the normality was confirmed using the Kolmogorov–Smirnov test. If the normality was violated, we used a signed-rank test for statistical comparisons. All statistical tests were performed with a significance level of 5% using MATLAB 2018a (MathWorks, Natick, MA, United States).

The relative mass of the segment, relative COM position, and relative inertial tensor around the COM of the phalanx, forefoot, and hindfoot segments are shown in Table 3. The size-normalized position vectors of the insertions, via points, and origins of the PAs are presented in Table 4. The mean and standard deviations of the masses of the phalanx, forefoot, and hindfoot segments of the 10 participants were calculated to be 0.131 ± 0.017 kg, 0.386 ± 0.051 kg, and 0.393 ± 0.052 kg, respectively.

The mean and standard deviations of the speed, cycle duration, and stance phase duration of the measured walking were 1.33 ± 0.17 m/s, 1.03 ± 0.05 s and 0.63 ± 0.05 s, respectively. The 3D joint angle profiles of the MTP, midtarsal, and ankle joints during walking are presented in Figure 3. The joint angles during quiet standing are also presented as a reference point. The joints angles were positive for eversion, dorsiflexion, and abduction. The ankle joint plantarflexed after heel-contact and dorsiflexed during single support phase so that the shank vaulted over the foot. The joint then plantarflexed while approaching toe-off and returned to its natural angle during the swing phase. Following heel contact, the midtarsal and MTP joints dorsiflexed and plantarflexed, respectively. During the single support phase, both joints remained largely unchanged; however, in the late stance phase, the midtarsal and MTP joints plantarflexed and dorsiflexed, respectively. Eversion−inversion and abduction−adduction of the two joints were relatively small, though slight eversion of the ankle and midtarsal joints and abduction of the ankle joint were observed during the stance phase.

Figure 4 displays the changes in the PA length and force during walking. The length profiles were generally consistent with each other for the four lateral PAs. They were stretched after the heel-contact and remained stretched until 50% of the gait cycle. During push-off, they were sharply shortened and remained shortened during the swing phase. PA1 was more largely stretched at 50% of the gait cycle than the other four PAs, mainly due to the elongation of the distal portion of the PA connecting the via point and insertion. The maximum tension forces exerted by the five PAs were 0.23 ± 0.06 N/kg, 0.15 ± 0.04 N/kg, 0.16 ± 0.04 N/kg, 0.17 ± 0.04 N/kg, 0.16 ± 0.04 N/kg, respectively. The resultant maximum net PA force was 0.87 ± 0.20 N/kg.

The joint moment and joint power profiles calculated with and without incorporating the PA are presented in Figure 5. Plantarflexion moments were generated by all three joints during walking, particularly in the late stance phase. The magnitudes of the peak plantarflexion moments were 0.18 ± 0.04 Nm/kg, 0.98 ± 0.18 Nm/kg, and 1.58 ± 0.18 Nm/kg for the MTP, midtarsal, and ankle joints, respectively. Slight inversion and abduction moments were generated in phase with the plantarflexion moment by the midtarsal joint. The power generated by the midtarsal joint was calculated to be positive, as the joint plantarflexed while generating the plantarflexion moment. In contrast, the power generated by the MTP joint was negative because the joint dorsiflexed while generating the plantarflexion moment. Therefore, positive and negative work was performed by the midtarsal and MTP joints, respectively, during walking. The peak positive power generated by the midtarsal joint and peak negative power absorbed by the MTP joint were estimated as 1.81 ± 0.73 W/kg and −1.01 ± 0.48 W/kg, respectively, during walking. The power in the eversion−inversion and abduction−adduction directions was almost zero.

Due to the presence of the PA model, the plantarflexion moments and power of the MTP and midtarsal joints were significantly smaller in magnitude in the model with incorporating the PA than that without incorporating the PA (p < 0.01). The abduction and inversion moments of the midtarsal joint during the stance phase were also significantly smaller in the model with incorporating the PA than that without incorporating the PA (p < 0.01). The mean peak plantar flexion moments of the MTP were 0.07 Nm/kg smaller (p = 0.01), and the mean peak inversion, plantar flexion, and abduction moments of the midtarsal joint were 0.02 Nm/kg, 0.2 Nm/kg, and 0.04 Nm/kg smaller, respectively (p < 0.01, p < 0.01, and p = 0.01, respectively) when incorporating the PA model (Figure 6A). The mean negative work done by the MTP joint and the positive work done by the midtarsal joint were significantly reduced by 3.0 and 2.2 J, respectively (p < 0.01 and p < 0.01, respectively) (Figure 6B).

Figure 7 compares the mean joint moment profiles of the three joints when the inertial parameters were altered. The change in the inertial parameters had virtually no effect on the calculated joint moments. Figure 8A compares the mean peak moments at toe-off of the midtarsal and MTP joints when the PA stiffness parameters were altered. The change in the PA stiffness altered the peak plantarflexion moments of the two joints but not the moments in the other directions. It was estimated that 13%–45% of the plantarflexion moment of the MTP joint and 8%–29% of plantarflexion in the midtarsal joints were generated by the PA at the time of push-off during walking (Figure 8B).