Static versus dynamic muscle modelling in extinct species: a biomechanical case study of the Australopithecus afarensis pelvis and lower extremity

Musculoskeletal models

A modern human (henceforth, just ‘human’) musculoskeletal model was used to evaluate the AL 288-1 models (n = 1, specimen ID: Subject03) (Charles et al., 2020). Subject03 was an adult female aged 26, weighing 72.6 kg and 176 cm in height. The Wiseman (Wiseman, 2023) Australopithecus afarensis model was used, based upon the AL 288-1 specimen, with the foot belonging to specimen AL 333-115 (Desilva et al., 2018); see further details in Supplementary Information 1. Both models include 36 muscles of the right side of the pelvis and right lower limb, crossing the hip, knee, ankle and metatarsophalangeal joints (Table 1). Both models each had 16 degrees of freedom (DOF): three rotational DOFs in the pelvis (pitch, tilt and roll), three rotational DOFs in the right hip (flexion-extension, adduction-abduction and long-axis rotation) and one DOF each in the right knee, ankle and metatarsophalangeal joints (flexion-extension). All other DOFs (i.e., in the left limb and translational) were fixed (locked).

Muscular parameter estimation

Seven Au. afarensis (AL 288-1) models were generated, each using different estimation methods for muscle parameters, alongside one human model that was subject-specific. These estimations for AL 288-1 produced either a model composed of muscles that ignored tendon compliance and muscle force-length-velocity relationships (referred to as the ‘static-muscle’ approach), or a model composed of ‘fully’ Hill-type muscles which included an elastic tendon and used force-length-velocity relationships (i.e., the ‘dynamic’ approach; see: Zajac, 1989; Millard et al., 2013). All muscles in the human model followed the dynamic approach. Four different methods were used to estimate architectural properties of the muscles in AL 288-1:

1. Alpha shape centroid estimation (method assumes an isometric ‘static-muscle’ model with no muscle pennation and was thus not suitable for use in the ‘dynamic-muscle’ model which requires pennation inputs)

2. Convex hull centroid estimation (static- and dynamic-muscle models)

3. Arithmetic centroid estimation (static- and dynamic-muscle models)

4. 3D model-based estimation (static- and dynamic-muscle models)

Each variant was based upon a dataset (n = 22) composed of published muscle architecture data from the human lower limb (Charles et al., 2020; Friederich & Brand, 1990; Charles, Moon & Anderst, 2019), which included muscle belly mass (m)_muscle), optimal fibre lengths (ℓ_o) and optimal pennation angle (α_o). This approach follows the estimation methodology developed by Bishop, Cuff & Hutchinson (2021a), but was modified here to include the 3D muscle model approach (Demuth, Wiseman & Hutchinson, 2023b) in which the polygonal muscles reconstructed for AL 288-1 (Wiseman, 2023) provided the muscle mass values for the subsequent calculation of the physiological cross-sectional areas and isometric force (F_max) in the 3D model-based variants. The approach was further modified to include α_o estimates (code provided in Supplementary Information 2). These approaches calculate the normalised muscle mass (m^∗), normalised fibre length (ℓ^∗) and α_o (non-normalised). The dynamic-muscle models ($m_{DynamicMuscleModel}^{\ast }$) defined m^∗as: (1)${\displaystyle m_{DynamicMuscleModel}^{\ast }=\frac{m_{muscle}}{m_{body}}}$

and for the static-muscle models ($m_{StaticMuscleModel}^{\ast }$), was defined as: (2)${\displaystyle m_{StaticMuscleModel}^{\ast }=\frac{m_{muscle\cdot cos\left({\alpha }_o\right)}}{m_{body}}}$

in which m_body was the individual’s body mass. ℓ^∗ was defined in both approaches as: (3)${\displaystyle {\ell }^{\ast }=\frac{{\ell }_o}{m_{body}^{0.33}}.}$

Three values were extracted for each muscle from the dataset (see example in Fig. 1): (1) the average of the input values provided the arithmetic value (i.e., the average of all data points), (2) a 3D convex hull was applied to the data points and the centre point of the hull provided the convex hull value calculated by weighting the centroid of each triangle by the respective face area in comparison to the total surface area of the convex hull—this was mathematically defined by Demuth et al. (2022), and (3) the largest circle which can be fit within the 2D convex hull produced an alpha shape of which the centre point provided the alpha value (Fig. 1). The arithmetic values were used as the input data for the 3D model variant, in which parameters were identical to the arithmetic approach, but the resultant PCSA and F_max calculations were based upon the polygonal muscle mass calculations. Whilst α_o was incorporated into the above equations for the static-muscle model, the value for α_o was left at 0° in OpenSim.

In total, our workflow produced seven different estimates of architecture per muscle (n = 36) in the AL 288-1 pelvis and lower limb. There were four static-muscle models based upon (1) an alpha shape, (2) a 2D convex hull, (3) the arithmetic centre, and (4) a 3D model; and three dynamic-muscle models based upon (1) a 3D convex hull, (2) the arithmetic centre, and (3) a 3D model. See Fig. 1 for an example of how different these parameters are between variants in one muscle example (data available: Wiseman, Charles & Hutchinson, 2024).

Each muscle’s physiological cross-sectional area (PCSA) for all seven models was calculated as: (4)${\displaystyle \text{PCSA}=\frac{cos\left({\alpha }_o\right)\cdot m_{muscle}}{\rho \cdot {\ell }_o}}$

where p is tissue density of 1,060 kg m⁻³ (Mendez & Keys, 1960), for the 3D muscle models, volumes were equivalent to m_musclep⁻¹. α_o was ignored in the dynamic-muscle model’s Eq. (4) because it is separately accounted for in the OpenSim pipeline (Seth et al., 2011; Millard et al., 2013; Cox et al., 2019). Each muscle’s maximal isometric force (F_max) subsequently was calculated as: (5)${\displaystyle F_{max}=PCSA\cdot \sigma }$

where σ is the muscle’s maximal isometric stress, with a representative value of 300 kN m⁻² used across all muscles (Medler, 2002; Wells, 1965; James, Altringham & Goldspink, 1995; Zajac, 1989; Maganaris et al., 2001; Saito et al., 2021; also Hutchinson, 2004a; Hutchinson, 2004b and references therein). The specific stress values used do not matter for our simulations, which do not attempt to estimate true in vivo maximal performance, but rather focus on data exploration. Tendon slack length (ℓ_S) was calculated for the dynamic-muscle model variants only, with the assumption that muscle fibres had the capability to range in length from 0.5 to 1.5 × ℓ_o throughout the muscle–tendon unit’s length across the pertinent limb joints’ ranges of motion (Manal & Buchanan, 2004). For the static-muscle models, tendon slack length was left at a default value of 1 m. All muscles were modelled in OpenSim using the Millard muscle class (Millard et al., 2013). For the dynamic-muscle models only, force-length-velocity properties were included through the Static Optimisation procedure. For the static-muscle models, this relationship was ignored.

For the human, muscle architecture from ‘Subject03’ was previously determined via diffusion tensor imaging and muscle paths and masses were established via muscle resonance imaging (Charles et al., 2020), producing a subject-specific model rather than a generic model. The human model is provided here as a comparative model to the australopith, rather than as an additional model for examining the influences of static versus dynamic muscles or different architecture estimates. The impacts of variations in muscle force-generating properties on inferences of muscle function in humans have been previously described (Charles et al., 2022; Kramer et al., 2022).

Inverse simulations

A previously captured walking cycle was used (Wiseman et al., 2022), tracked and scaled to the musculoskeletal model of the human model using the ‘model scale’ and ‘inverse kinematics’ tools in OpenSim 4.3 (Delp et al., 2007; Seth et al., 2018). In this previous study, one participant walked with a typical walking posture (five trials along a 12 m long trackway) at a speed of 1.0 m/s and a 14-camera optoelectronic 3D motion capture system (250 Hz, Oqus Cameras, Qualysis AB, Gothenburg, Sweden) was used to capture kinematics via a skin-attached reflective marker-set. A mid-stance (30% of the gait cycle), posture was extracted from each right limb’s stride, averaged and applied to the human model in OpenSim (Supplementary Information 3). The left limb was modelled in a mid-swing pose at 30% of the gait cycle, with toe clearance ensured. Each left limb segment had inertial properties, but all DOFs were locked.

These joint angles were used as a guideline to position the AL 288-1 model in a similar mid-stance pose. However, it was not possible to directly apply the joint angles from the human to AL 288-1 due to geometric differences (e.g., Gatesy & Pollard, 2011) which resulted in the foot mis-aligning with the ground (i.e., a flat foot on the ground was a requirement of this study). This possibly was due to differences in iliac flaring of the pelvis and lower limb proportions (Lovejoy, 2007; Lovejoy, 2005; Jungers, 1982; Stern & Susman, 1983; Ward, 2002; Wiseman, 2023) influencing foot placement. Slight manual modifications were made to the hip joint angle to position the limb with the foot flat on the ground (i.e., greater hip abduction and knee adduction were required; knee adduction angle was subsequently locked). These joint angles were applied to all seven AL 288-1 models.

To maintain dynamic congruity between the model and the experimentally generated ground reaction forces (GRFs) and kinematic data (e.g., the hypothetical positioning of AL 288-1 into a mid-stance single-limb support pose) during static optimization (i.e., consistency between the simulated model’s GRFs and kinematic data), residual actuators—i.e., additional forces that help to achieve dynamic congruity by adjusting the model to match the input data more closely during static optimisation—were applied (M_x,y,z and F_x,y,z, in which the x-axis is anterior-posterior, the y-axis is perpendicular to the ground, and the z-axis is medio-lateral) to the pelvis-body COM (see: Hicks et al., 2015). The activations of these actuators during the simulation are heavily penalised by the optimisation algorithm and were thus only recruited if necessary. Results from these residual actuators as outputs from the simulations can be found in Supplementary Information 4, showing that they were well below the tolerances recommended by Hicks et al. (2015).

Following Bishop, Cuff & Hutchinson (2021a), gravity was increased accordingly relative to m_body during the static simulations to produce net force balance and to help achieve dynamic consistency (see Fig. 2). For example, in the australopith, the vertical GRF (vGRF; see below) applied at 1.6 * BW was 544.171 N and the corresponding gravitational force was changed to 15.691 m/s². By increasing gravity relative to m_body, the gravitational force acting on the body is artificially enhanced. This increase in gravitational force helps achieve balance with the applied vGRF (Bishop, Cuff & Hutchinson, 2021a). This ensured that there was force equilibrium (i.e., all forces were balanced, assuring no net acceleration) when the applied vGRF was smaller or greater than 1 * BW.

Vertical and medio-lateral GRFs were applied to the static inverse simulations (all joints had zero accelerations and velocities). No antero-posterior GRFs were appended. The centre of pressure of the foot (point of GRF origin) was assumed to be at the midpoint between the COMs of the pes (defined here cumulatively as the hindfoot, midfoot and metatarsals) and digit segments. The medio-lateral GRF vectors were calculated as the ratio between the vertical and mediolateral vectors from a large dataset (GAITREC) of published GRF profiles (n = 211 healthy participants; n = 2,382 trials in total) extracted at 30% of the gait cycle at a controlled walking speed of 0.98 m/s (Horsak et al., 2020b; Horsak et al., 2020a). Because the data were not evenly distributed, the median ratio of 0.011265 was used rather than the mean. The mediolateral GRFs were thus appended in AL 288-1 as a fraction of the vGRF. The vGRFs were input as a factor of body weight (BW) starting at 0.2 * BW and increased by increments of 0.2 * BW (i.e., see Supplementary Information 1) until the static optimization algorithm used for the inverse simulations could no longer find a solution to achieve static equilibrium for any combination of muscle activations (i.e., simulation failure), while minimising the sum of squared activations for each DOF (n = 16) (Rowninshield & Brand, 1981).

Evaluation of simulations

Simulated muscle activations from the human were qualitatively compared to published electromyography (EMG) studies that report on walking in healthy, non-pathological individuals (Cappellini et al., 2006; Van Criekinge et al., 2018) to evaluate the reliability of the human simulation (Hicks et al., 2015). We inspected the mid-stance on-off timings of the EMG signals, but ignored magnitudes because this study only modelled a static, mid-stance pose (See: Supplementary Information 5). Subsequently, we set the performance criteria of the AL 288-1 simulations (e.g., Demuth, Wiseman & Hutchinson, 2023b) as follows: First, we used the maximal amount of vGRF exertion in the human as a benchmark for the expected vGRF exertion (vs. BW) in the AL 288-1 simulations, assuming that two phylogenetically closely related and morphologically similar individuals employing similar locomotor behaviours should be able to sustain comparable forces on an extended limb. If a simulation was unable to support the body below this benchmark, the underlying model was assumed to be ‘weak’ in light of the current setup in which a hypothetical limb pose was modelled to examine the range of outputs that can be generated if muscle architectural parameters are changed. A different limb posture would likely generate different results in both the human and AL 288-1 individuals and would lead to the selection of a different benchmark value.

Second, the operating fibre lengths, muscle specialisation (quantified as the ratio of ℓ_o to PCSA) and simulated activations of the AL 288-1 models were compared to those of the human model and simulations to provide preliminary insights into functional differences between species.

How do the different estimation methods for muscle parameters influence maximal performance (sustainable ground reaction force)?

The human model was able to exert a maximum of 2.2 * BW vGRF on a single limb at 30% of the gait cycle. The maximally exerted vGRF in the AL 288-1 simulations differed according to model variant (Table 2). The static-muscle models typically weakly supported the capability of this representative mid-stance pose compared with the human results as a baseline, with some AL 288-1 variants unable to withstand >1.8 * BW vGRF on a single limb. Contrastingly, the static-muscle model 3D variant withstood 3.6 * BW vGRF, indicating that the muscle mass estimations following the alpha and 2D convex hull approaches possibly were underestimated and thus under-performed. However, the 3D muscle modelling approach produced muscles with greater force capacities and greater sustained vGRFs. The dynamic-muscle models were able to withstand greater vGRF than the static-muscle models, excluding the 3D variant. All dynamic-muscle models indicated that the limb was capable of supporting greater vGRF than the human, ranging between 2.6 * BW (the convex hull variant) to 2.8 * BW (the 3D variant) vGRF. Notably, the dynamic-muscle 3D variant was markedly weaker than the static-muscle 3D variant. Overall, the average amount of vGRF that the limb could withstand in the static-muscle models was 2.2 * BW and in the dynamic-muscle models was 2.6 * BW. Only the static-muscle models (excluding the 3D variant) were below the 2.2 * BW benchmark and are thus assumed to be ‘weak’.

Do the different estimation methods for muscle parameters produce muscles with different inferred architectural specialisations (force versus velocity?)

We used the ratio of ℓ_o to the PCSA to investigate how the different variants influenced whether a muscle was specialised more for force (∼small ratio value) or for velocity (∼high ratio value) (Fig. 2). Whilst is impossible to dictate a ‘cut-off’ value to distinguish between low versus high ratios, the general trend across muscles can instead be examined to elucidate differences (i.e., Charles et al., 2022). Small changes in ratio are expected because architectural parameters are not identical between each estimation method or variant, but larger changes between variants indicate major functional changes of a muscle—for example, if a variant’s ratio is twice as large, then we can interpret this as a difference in predicted muscle specialisation. Generally, most muscles across the variant spectrum appeared similarly specialised (small changes such as those found in the GRA, GMed, SM, ST, EDL, PL, PB, and SOL muscles, for example, were negligible). Some muscle specialisations switched depending on the variant used. For example, the static-muscle model’s variant estimates for the MG (Fig. 2D) predicted that this muscle was velocity-specialised (∼higher ratios), but the dynamic-muscle model’s variants for the same muscle instead predicted that the muscle was suited for force-generation due to lower ratios that were 3.2× smaller than the ratios of the static-muscle variants. The static-muscle 3D-based model of the TP (Fig. 2E) predicted a velocity-specialised muscle, yet all other variants predicted a force-specialised TP (i.e., the 3D-based variant was almost 3× the amount of the other ratios, a profound difference). In the VL (Fig. 2C), the dynamic-muscle variants suggested that this muscle was a velocity-specialised muscle, but the ratios of the static-muscle model variants were 2.2× greater; more force-specialised. In the GMax, both static-muscle and dynamic-muscle models of the 3D variant were found to have higher ratios (velocity-specialised) than all other variants. Overall, there was broad similarity among variants with only a few aforementioned outlying differences.

Do any of the methods produce a muscle which could generate relatively greater forces or operate on a different part of its force-length curve?

We sought to determine if any of the variants produced a muscle which generated relatively greater forces. Using the plot function in OpenSim, we plotted each muscle’s ℓ^∗ from the three dynamic-muscle models of AL 288-1 and also for the human for comparison of ℓ^∗ only (Fig. 3). In the AL 288-1 models, the convex hull and arithmetic variants had comparably lower muscle force (maximum: ∼330 N and ∼220 N, respectively; Figs. 3C, 3D) than that of the 3D variant (∼580 N; Fig. 3B).

Next, we assessed if a muscle’s force-length curve was the same or different between each of the variants based on their normalised fibre length values from the chosen single time point of a walking stride A value >1 was on the descending limb, and a <1 was on the ascending limb. The patterns of where fibre lengths lay on their force-length curves were broadly similar between the convex-hull and arithmetic variants, with some differences in the 3D model-based variant. For example, the AB operated on the same region of its force-length curve in all variants, but had a muscle force of around 10 N in the 3D model-based variant (Fig. 3B) in comparison to greater force in the other variants, ranging between 75–250 N (Figs. 3C, 3D). The AM was on the descending limb (i.e., actively lengthened) for the 3D model-based variant (Fig. 3B), but around the plateau of the curve (near-isometric) for the other variants (Figs. 3C, 3D). The AL was at the plateau of the curve in the 3D model-based variant (Fig. 3B), but instead was on the descending limb in the other variants (Figs. 3C, 3D).

How do these different methods affect simulated influence muscle activation patterns (i.e., different groups activated, or the same muscle to different activation levels)?

Simulated muscle activations during single-limb stance in the AL 288-1 models are reported in Fig. 4. Muscles are heuristically simplified into primary functional groups. The static-muscle models (excluding the 3D variant) were broadly comparable with heavily activated ankle plantarflexors (LG, PB and SOL) and moderately activated hip adductors, hip extensors and knee extensors. In the 3D static-muscle variant (Fig. 4A), there was increased activation in comparison to the other static-muscle variants (Figs. 4B–4D) in the hip adductors and knee extensors. The 3D static-muscle variant also had different patterns in the plantarflexors in which maximal activation only occurred at maximal vGRF (i.e., just before simulation failure), rather than being activated earlier in the simulations at ∼0.8 * BW vGRF as per the other static-muscle variants.

The muscle activations in the dynamic-muscle variants were mostly similar to the static-muscle models with few differences. In the dynamic-muscle convex (Fig. 4G) and 3D variants (Fig. 4E), the TP and FHL muscles were maximally activated, but had low activation in the dynamic-muscle arithmetic variant (Fig. 4F). The dynamic-muscle arithmetic variant (Fig. 4F) was the only variant to not have maximal activation of the PB. Activations patterns overall were comparable. The activation patterns of the SM and ST were switched between the dynamic-muscle convex (Fig. 4G) and arithmetic variants (Fig. 4F), in which the SM had greater activation in the arithmetic variant, but in the convex hull the ST muscle activation was greater, and the convex and arithmetic variants also were the only ones to have highly activated hip adductors. Both 3D variants were the only variants to not activate the SM and ST muscles. The dynamic-muscle 3D variant (Fig. 4E) had reduced hip extensor activation in comparison to the other dynamic-muscle variants.

In Fig. 5, the maximum activation of each muscle (i.e., from maximal sustainable vGRF) from each of the variants is presented, alongside the average activation and standard deviation across all variants. Whilst activations of most muscles across each of the variants are broadly comparable, such as the SOL which was fully activated (full activation = 1) across all variants, other muscles had activations that depended more on the variant. The TP was highly activated in the dynamic-muscle 3D and convex hull variants, but weakly activated in all other variants at the point of maximal vGRF exertion. The PL, TA and PB had similar patterns, with activations that were variant-dependent.

Next, we calculated the mean activations in each major muscle group (see Supplementary Information 6). All variants mostly produced comparable mean activations per muscle group with only minor deviations, such as increased activations of hip extensors in the dynamic-muscle arithmetic variant in comparison to the other variants. Only the 3D (static- and dynamic-muscles) and dynamic-muscle convex hull variants had a minimally activated ankle dorsiflexor group, which was inactive (activation = 0) across all other variants.

What are the differences in muscle activations and force length states between species?

The simulated muscle activations in the human model were consistent with published EMG studies of muscle excitation at mid-stance in healthy individuals (Perry & Burnfield, 2010; Van Criekinge et al., 2018; Cappellini et al., 2006; Wall-Scheffler et al., 2010) and are reported in Supplementary Information 5. As such, the human simulation was considered sufficiently representative.

In both the human and the AL 288-1 simulations, it was the ankle plantarflexor group that limited maximal vGRF exertion (Fig. 4). The hip adductors were only activated in the australopith, not in the human, alongside more moderate activations in the hip extensor and knee extensor groups. The human was the only model in which the PIRI muscle was fully activated. Overall, more muscles in the AL 288-1 limb were required to maintain single-limb stance than in the human, but these muscles were only mildly to moderately activated.

Similar patterns in the simulated operation of ℓ^∗ on its force-length curve were observed between the human and AL 288-1 variants for the mid-stance pose (Fig. 3). For example, the MG in all variants generated high force in comparison to the other muscles and was on the descending limb of its force-length curve. The FDL and EDL were the same (on the descending limb) in all AL 288-1 variants, but different in comparison to the human, in which the FDL and EDL instead were on the ascending limb (actively shortened) or were at the plateau of the curve (isometric) at 30% of the gait cycle. The vastus muscle group was on the ascending limb in the human, but around the plateau of the curve for all AL 288-1 variants in which this muscle group was weakly activated. The SM was on the ascending limb for the human (but inactive) whereas it on the descending limb for all AL 288-1 variants. The SOL was activated in all AL 288-1 variants, but only weakly activated in the human model (descending limb), and heavily activated in the AL 288-1 variants (force-length plateau).