Sensitivity of musculoskeletal models to variation in muscle architecture parameters

Abstract

Musculoskeletal models, like all theoretical models of physical processes, depend on the assumptions needed to construct the model. For musculoskeletal models, these assumptions include, among other things, the kinematic data, the kinetic data and the muscle parameters. The former (dynamic) data can be acquired relatively easily from living subjects, but the latter are usually based on limited information, frequently determined from cadaver studies performed on elderly individuals. Previously, we determined the sensitivity of forces to dynamic differences among 10 humans walking on a straight path. Here, we assess the sensitivity of the muscle and joint reaction forces developed in human walking to variable muscle parameters obtained from 10 living adults, whose data were recently reported, and compared the results with the values from a standard model that depends on cadaveric data. We found that, while the force patterns across the stance cycle were similar among muscle parameter models, differences of as much as 15% in the force magnitude were produced. Whether or not the variation between the standard model and other muscle parameters is important depends on why the forces are required.

Social media summary: Typical human variation in muscle architecture changes internal forces by 15% in a musculoskeletal model of walking

The nature of the biomechanical requirements and consequences of bipedalism, from its inception in the hominin lineage to its manifestation in modern humans, remains a puzzle that motivates the research of many biological anthropologists, anatomists and functional morphologists. Bipedalism is rare among mammals, but ubiquitous within hominins for at least 4 million years (Harcourt-Smith, 2010), and because modern humans are the only remaining habitual, terrestrial bipedal primate, much work has centred on the bipedalism of modern humans. Insights into how human movement influences the musculoskeletal system has been extrapolated onto the movements of extinct hominins (e.g. Kramer & Eck, 2000; Lovejoy, 1988; Ruff & Higgins, 2013; Sellers, Cain, Wang, & Crompton, 2005). Inherent in this approach, however, is the notion of a ‘generic’, i.e. representative, model of humans, human musculoskeletal systems and human movements. In order to determine the limitations of such a generic approximation, the impact of human variability on the outcome of any particular analysis should be examined.

Our ultimate goal is to assess the impact of the forces generated in walking on the bones of the lower limb. Given that there is no one point in stance where every joint reaction and muscle force is at its maximum (Sylvester, Lautzenheiser, & Kramer, 2021b), we seek to develop a reduced set of generic conditions (each of which represents a point in the stance cycle) that are likely to be relevant to bone morphology. Here, ‘generic’ refers to conditions that are representative of humans, rather than a specific human. This reduced set of conditions that represents many points during stance allows us to create an envelope of cases. We can do this by identifying the points in the stance phase where the applied forces are maximum and whether or not these points are consistent across the stance phases within and among people. A key issue in this approach, however, is that the conditions are representative of humans. Consequently, we are critically concerned about human variability in kinematics and kinetics (motion and forces) and in musculoskeletal anatomy. Models are necessary because empirically measuring muscle forces in living animals depends on the surgical implantation of sensors into the body, which is invasive and expensive, and requires specialised equipment and personnel (e.g. Bey & Derwin, 2012; Pedersen, 1997).

Models are, then, approximations that are designed to represent phenomena that are difficult or impossible to investigate directly. Although models are useful to begin to gauge processes that cannot be empirically investigated, all such approximations have assumptions inherent in their development. Assessing the influence of these assumptions on outcome variables is critically important. Using kinematic and kinetic information from 10 individuals, we have evaluated the muscular forces among individuals (Sylvester et al., 2021b) and found that the functional muscle groups that are critically important to propel (i.e. the ankle plantarflexors) and stabilise (e.g. the hip abductors) the body are consistent among individuals and steps. Inherent in our use of this musculoskeletal model, however, were many assumptions, including those about the values of muscle volumes and optimal fibre lengths. In this work, we seek to determine the degree to which joint reaction and muscle forces of the lower limb vary with different muscle parameters.

Muscle forces are difficult to determine in musculoskeletal models because they exist in functional groups that are multiply redundant with respect to the joint of interest, meaning that multiple muscles can produce forces that create moments that act at particular joints and about particular axes. For example, soleus, plantaris and medial and lateral gastrocnemius all produce ankle plantarflexion. Consequently, the only way to resolve individual muscle forces (i.e. to solve the muscle redundancy problem) is to assume how the body allocates muscle force or, in other words, how the body determines the degree of active contraction to produce in each muscle at any given point. This physical system is difficult to study empirically, so modern musculoskeletal models use muscle strength (derived from muscle parameters such as fibre length, muscle volume, and pennation angle) to determine activation level via an allocation or cost algorithm. This allocation algorithm is the critical assumption for determining the force produced by a particular muscle because muscles in a functional group are all capable of producing forces that produce similar joint moments. Given that our goal is to use muscle forces to interpret bone morphology, the force that a particular muscle produces is critical because some muscles in a functional group originate or insert on different bones than other muscles in the same functional group. For example, soleus connects the tibia to the calcaneus across the ankle and subtalar joints, but medial and lateral gastrocnemius and plantaris originate on the femur, crossing subtalar, ankle and knee joints. Even though the primary action of all the muscles of the triceps surae complex is to cause the ankle to plantarflex, for examining femoral evolution, the distribution of force among the muscles is important because soleus does not directly interact with the femur.

Our musculoskeletal model produces muscular forces that are consistent among steps, trials and individuals for functional muscle groups (Sylvester et al., 2021b). We have, therefore, some understanding of the influence of variability in the kinematics and kinetics that is due to inter- and intra-individual differences on joint reaction and muscle forces. However, musculoskeletal models, including ours, are typically constructed using muscle architecture and geometry parameters drawn from either a single individual or small and morphologically homogeneous samples (e.g. Klein Horsman, Koopman, van der Helm, Prosé, & Veeger, 2007). Nonetheless, even for models of muscle function that are simpler than the Hill type, such as the muscle formation that we use, muscle architectures (e.g. optimal fibre length, muscle volume) are foundational inputs as they determine strength. The majority of reference datasets are drawn from elderly cadaveric specimens (e.g. Klein Horsman et al., 2007), on whom muscle dissection and examination are used to obtain the geometric parameters. These reference datasets are then used to construct generic musculoskeletal models which are assumed to be representative of the general population. How representative such limited samples are of the ‘average’ healthy adult remains an open question, however, as muscle architecture and geometry, including volume, fibre length and pennation angle, have been shown to change with age (Handsfield, Meyer, Hart, Abel, & Blemker, 2014; Narici, Maganaris, Reeves, & Capodaglio, 2003; Tate, Williams, Barrance, & Buchanan Thomas, 2006), among other individual characteristics.

Recently, Charles et al. (2019) used diffusion tensor magnetic resonance imaging and fibre tractography to visualise with a high degree of resolution the in vivo anatomy of a variety of muscle groups and to create one of the most comprehensive in vivo lower limb muscle architecture datasets available to date (Charles, Suntaxi, & Anderst, 2019). This novel dataset includes muscle architecture measurements from 20 lower limb muscles in 10 young, healthy individuals, with muscle volume, length, optimal fibre length, pennation angle, muscle physiological cross-sectional area and maximum isometric force provided for each subject. Continuing this work, Charles et al. (2020) compared the accuracy of both subject-specific lower limb models and generic models built using data from both elderly and young individuals with muscle values obtained experimentally using an isokinetic dynamometer (Charles, Grant, Août, & Bates, 2020). Their results demonstrated that, while the generic models were not as accurate as the subject-specific models in predicting the experimental muscle values of particular individuals, generic models still have utility in musculoskeletal modelling studies. One of the major limitations of this study was that the modelling predictions were only validated in comparison with muscle forces obtained in the sagittal plane involving one degree of freedom. It remains unclear how different subject specific and generic musculoskeletal models are for making functional predictions when simulating more complex movements such as gait.

Consequently, in this work we evaluate the impact of variability in muscle parameters and the choice of exponent in the cost algorithm on muscle functional groups of the lower limb that have both uniarticular and multiarticular muscles. Muscle functional groups are characterised by their primary function, although individual muscles within the group can have secondary functions. We identify the following functional muscle groups: the triceps surae (plantaris, soleus, and medial and lateral gastrocnemius), the hamstring muscles (semitendinosus, semimembranosus, biceps femoris caput breve, and biceps femoris caput longum) and the quadriceps (vastus lateralis, vastus intermedius, vastus medialis, and rectus femoris). Changing muscle parameters within a model does not change the external kinematic and kinetic conditions, so joint moments do not change. However, because muscle parameters determine muscle strength and the cost algorithm distributes force to muscles based on their strength, the distribution of force among the various muscles in a muscle functional group or among groups could change and this could change joint reaction forces. Consequently, although our approach is predominantly data exploration, we hypothesise that:

1. (1) The shape of the joint reaction and muscle force profiles with respect to stance time will be similar among muscle parameter models in the same individual and among all individuals.

2. (2) The distribution of maximum force within a functional group will vary based on strength of the individual muscles relative to the other muscles in that functional group.

Materials and methods

Dynamic cases: participants and protocol

Ten healthy participants (five female, five male; age 21–55 years; body weight 53.1–98.5 kg; stature 1.59–1.82 m) were recruited from the student population at the University of Washington. Exclusion criteria included gait abnormalities and a recent history of injury. Biometric data were collected from each participant, including stature, mass, limb segment length and limb segment circumference. The University of Washington's Institutional Review Board approved all aspects of this study (IRB no. STUDY00001125). The experimental protocol for the kinematic data collection process has been described in detail elsewhere (Sylvester, Lautzenheiser, & Kramer, 2021a, 2021b). We refer to the motion and force data of these participants as dynamic cases to distinguish them from subjects whose muscle parameters were used in the simulations described below.

The participants walked at their self-selected normal velocity across four force plates (Kistler, Switzerland) while the motions of 30 infrared-reflective markers placed on anatomical landmarks were measured with a 10 camera motion capture system (Qualisys, Sweden) in the Amplifying Movement and Performance Laboratory at the University of Washington. Each participant completed five trials, but the data from two trials for two participants were unusable, resulting in 46 total dynamic cases with two consecutive stance periods (92 total steps).

Baseline Musculoskeletal Model

The baseline model used to calculate lower limb joint reaction and muscle forces in this study was the MoCap model included as part of the commercially available AnyBody Modelling System (v. 7.3, AnyBody Technology, Denmark) and hosted on the AnyBody musculoskeletal model repository (AnyBody Managed Model Repository AMMR v. 2.3.0). This is a validated, multi-trial, full-body, motion-capture-driven human gait model with head, trunk, and left and right lower limb components (TLEM 2.0; De Pieri et al., 2018). The muscle architecture and geometry data for the baseline model were drawn from a single cadaveric individual, described in detail in Klein Horsman et al. (2007). We refer to this baseline MoCap model as the KH model.

Each lower limb of the KH model comprises thigh, patella, shank, talus and foot segments, and has six total degrees of freedom including all three rotations at the hip and one each at the knee (flexion/extension), ankle (plantarflexion/dorsiflexion) and subtalar (inversion/eversion) joints. The pelvis (relative to the ground) has six degrees of freedom, three translational and three rotational. Thirty-seven anatomical muscles are represented by 169 muscle elements (force actuators) in each lower limb (e.g. gluteus medius is composed of 12 separate muscle element actuators).

As detailed in our previous reports using these dynamic cases (Sylvester et al., 2021a, b), after initial scaling to the participant's body mass and stature, segment lengths are further adjusted to improve the match of the model with the motions described by the markers. This results in an individualised model for each participant based on their anthropometrics and motions. Then, an inverse dynamics solution of the individualised model is obtained for each trial and the muscle forces are allocated. Each trial produces a set of joint reaction and muscle forces. Variation within an individual is due to slight differences in motions or ground reaction forces, while variation among individuals includes differences in anthropometrics as well as the kinetic and kinematic variation. Thus, each dynamic model is ‘subject-specific’ for those input parameters.

We use a muscle recruitment algorithm which seeks to minimise the sum of the ratios of each muscle force to its strength raised to a power:

(1)$${\rm \Sigma }( F/N) ^p$$

where F is the muscle force, N is the muscle strength and p is the power. Both N and p are model assumptions. We used the strength algorithm with an exponent of 3 as the baseline, consistent with our previous work (Sylvester et al., 2021a, b) and others (e.g. Simonsen et al., 2021). We also completed a set of simulations using a strength algorithm with an exponent of 2. We distinguish these versions of the same musculoskeletal model with KHP3 and KHP2, respectively.

New muscle parameter model

The data from the 10 participants (Table 1) from Charles et al. (2019) were used to create MoCap models with revised muscle parameters. We refer to these models as CS1, CS2, and so on, where the number indicates the Charles et al. (2019) subject ID. Two main parameters that influence muscle force generation were identified as of interest for our simulations: optimal fibre length (L _opt; Table 2) and muscle volume (muscvol; Table 3). To create the 10 new MoCap models modified to include muscle architecture parameters drawn from Charles et al. (2019), we adjusted the parameters of interest to the size of the KH model. Specifically, to adjust muscle volume, we summed the total volume of muscle for the KH model based on the published values provided for the KH cadaveric specimen (Klein Horsman et al., 2007). We then calculated the ratio of muscle volume for each muscle and muscle element (proximal, mid and distal; inferior and superior; medial and lateral depending on the muscle) to the total muscle volume in the Charles et al. (2019) subjects and multiplied by the total muscle volume of the KH baseline model. CS muscle volume in model = (CS muscle volume/CS total muscle volume) × KH total muscle volume.

Table 1. Characteristics of subjects with muscle parameter data (Charles et al., 2019)

Table 2. Optimal fibre lengths (mm)

Table 3. Muscle volumes (ml)

To scale optimal fibre length (L _opt), a similar process was undertaken using stature: the stature of the KH model was divided by the stature of the Charles et al. subject and multiplied by the optimal fibre length values provided in Charles et al. (2019) for that subject. As Charles et al. (2019) does not distinguish among the muscle regions for all muscles included in the MoCap model (e.g. soleus is divided into medial and lateral portions in the MoCap model, but not in the Charles et al., 2019 data), we elected to use the same L _opt for each region.

The adjustment of these two muscle parameters generated 10 new, modified models – one model based on each of the 10 subjects from the Charles et al. (2019) study. Table 2 lists the optimal fibre length and Table 3 details the muscle volumes used for the baseline KH model and each of the modified CS models in this study.

Model simulation

Each CS model was simulated using the gait kinematic data drawn from each of the 10 participants in the experimental walking protocol described by Sylvester et al. (2021b), resulting in 460 simulations. To this dataset were added the baseline muscle parameter models (KHP3 and KHP2), resulting in a total of 552 simulations and 1104 steps. All CS simulations used an exponent of 3 in the strength algorithm.

Sensitivity analysis

We created three muscle functional groups that contain both uniarticular and multiarticular muscles that cross the knee (Table 4). We visually assessed the patterns of the joint reaction forces and muscle forces within dynamic cases for all muscle parameter models and among individual averages for all muscle parameter models for the functional muscle groups. We determined the strength of the functional muscle groups and the uni- and multiarticular portions by summing the strengths of the muscle elements of the group. The strength ratio of each functional muscle group is the strength of the uniarticular muscles divided by the multiarticular muscles. We determined the point in the stance of the maximum total muscle force of a functional group and determined the ratio of the force of the uniarticular muscle(s) to that of the multiarticular muscle(s) at that point.

Table 4. Definition of muscle functional groups

^a Plantaris is not traditionally included in triceps surae, but we include it here as it has similar function to soleus and gastrocnemius.

Although our focus here is on data exploration – e.g. whether the joint reaction and muscle force profiles are consistent – to address hypothesis 2, we predicted the maximum force ratio from the strength ratio in linear and quadratic forms, using linear regression with correction for repeated measures (i.e. the dynamic cases; Stata, StataCorp, College Station, TX USA), to gain insight into how the muscle force allocation algorithm influences the joint reaction and muscle forces.

Results

Our first hypothesis is that muscle parameter models, i.e. the CS models and KH, would all exhibit similar shapes of their joint reaction and muscle functional group force curves. Figures 1–3 (muscle forces for the triceps surae, quadriceps and hamstring muscle functional groups, respectively) and Figure 4 (magnitudes of the knee joint reaction force; similar figures for the hip and ankle are provided in the Appendix, Figures A1 and A2) demonstrate that our expectations are broadly confirmed for each dynamic case. The muscle forces of the triceps surae group (Figure 1) and the joint reaction forces (Figure 4, A1 and A2) are consistent in pattern within a dynamic case and among them. The triceps surae force peaks in late stance (~80% of stance) and the joint reaction forces also exhibit peaks around 80% of stance. Of note, however, is that the shapes of the curves, while similar among muscle parameter models within a dynamic case, can vary considerably among dynamic cases. The muscle forces of the quadriceps group peak in early stance (~20% of stance) but dynamic case 7 (< 250N) exhibits much lower force peak than do cases 2 or 9 (> 1000N; Figure 2). The shapes of hamstring muscle functional group (Figure 3) demonstrate this variability: in dynamic cases 8 and 9 a distinct peak in hamstring force at 70–80% of stance is clear, while no such peak is apparent in dynamic case 1. Additionally, while the shapes of the curves are consistent, the difference in joint reaction force magnitudes within a dynamic case for the various muscle parameter models can exceed 700N or greater than 15% (e.g. hip forces for dynamic case 9; Figure A1).

Figure 1. Force in triceps surae. For the CS subjects, males are indicated with a solid line while females are indicated with a dashed line. The colour of the line indicates the CS subject (Table 1): CS1 and CS7 are shown in red; CS2 and CS6 are shown in green; CS3 and CS8 are shown in blue; CS4 and CS10 are shown in magenta; CS6 and CS9 are shown in cyan. KHP3 is shown with a solid black line, while KHP2 is indicated with a dotted black line.

Figure 2. Force in quadriceps. For the CS subjects, males are indicated with a solid line while females are indicated with a dashed line. The colour of the line indicates CS subject (Table 1): CS1 and CS7 are shown in red; CS2 and CS6 are in green; CS3 and CS8 are in blue; CS4 and CS10 are in magenta; CS6 and CS9 are in cyan. KHP3 is shown with a solid black line, while KHP2 is indicated with a dotted black line.

Figure 3. Force in hamstrings. For the CS subjects, males are indicated with a solid line while females are indicated with a dashed line. The colour of the line indicates the CS subject (Table 1): CS1 and CS7 are shown in red; CS2 and CS6 are shown in green; CS3 and CS8 are shown in blue; CS4 and CS10 are shown in magenta; CS6 and CS9 are shown in cyan. KHP3 is shown with a solid black line, while KHP2 is indicated with a dotted black line.

Figure 4. Knee joint forces. For the CS subjects, males are indicated with a solid line while females are indicated with a dashed line. The colour of the line indicates the CS subject (Table 1): CS1 and CS7 are shown in red; CS2 and CS6 are shown in green; CS3 and CS8 are shown in blue; CS4 and CS10 are shown in magenta; CS6 and CS9 are shown in cyan. KHP3 is shown with a solid black line, while KHP2 is indicated with a dotted black line.

The second hypothesis, that variation in muscle strength predicts variation in the distribution of the forces to bi- and uniarticular muscles, was not rejected. Figure 5 demonstrates the relationship between peak muscle force ratio and strength ratio for the three muscle functional groups. In triceps surae, CS and KHP3 follow a strong trend (r ² = 0.91, p < 0.001). In the hamstring muscles, CS and KHP3 also follow a strong trend (r ² = 0.71, p < 0.001). In the quadriceps muscles, CS and KHP3 demonstrate the same trend (r ² = 0.12, p = 0.002; Figure 5).

Figure 5. Peak muscle forces vs. strength in triceps surae: (a), quadriceps (b) and hamstrings (c).

The KH muscle parameter model with an exponent of 2 for the cost equation can change the hip and knee joint reaction forces relative to using an exponent of 3 with slightly higher knee, but slightly lower hip forces (Figures 4 and A1). The effect of the exponent on the ankle is minimal (Figure A2). An exponent of 2 produces lower triceps and hamstring muscle forces, but higher quadriceps forces (Figures 1–3). An exponent of 2 produces a different relationship between uni- and multiarticular muscles than has an exponent of 3 for triceps and hamstring muscles, but not for quadriceps (Figure 5).

Discussion

Our ultimate goal is to produce a generic set of joint reaction and muscle forces that can be used to interrogate the effect of variation in bone morphology that is seen in the fossil record. We had previously found consistency within and among individuals in the forces exerted by the propulsive and stability-focused muscle functional groups when we used the same muscle parameters (Sylvester et al., 2021b). Consequently, we conducted the current study to evaluate the impact of varying muscle parameters on the joint reaction and muscle forces and found that the shapes of the curves, especially within a dynamic subject, are similar, although the forces within a subject can vary among muscle parameter models.

The triceps surae muscle force for the KHP3 muscle parameter model is consistently greater than that for all of the CS muscle parameter models. Interestingly, the KH muscle parameter model has very large muscle volume values for soleus compared with the CS muscle parameter models: medial segment = 414.8 ml; lateral segment = 378 ml in KH vs. medial segment = 214.9–281 ml; lateral segment = 195.8–256.1 ml in the CS models. Gastrocnemius muscle volumes are more similar between KH and CS muscle parameter models. Because muscle volume is a contributor to muscle strength and because muscle strength determines the distribution of force among the available muscles, the larger soleus muscle volume in the KH model potentially draws more force to the soleus in the muscle redundancy algorithm. Similarly, the smaller soleus muscle volumes of the CS models cause the algorithm to allocate more muscle force to the medial and lateral gastrocnemius. In the KH models where soleus is larger (higher strength) and, therefore, is allocated more force, the soleus's somewhat shorter moment arm requires a higher muscle force, which increases joint reaction forces of the ankle (Figure A2).

The muscle forces in the quadriceps group in the KHP3 model are similar to, or somewhat higher than, the CS muscle model parameter models in the earlier parts of the stance phase, but smaller in the later portion of stance. The hamstring group demonstrates a more complex picture that depends on the portion of stance and the particular CS muscle parameter model. The knee joint forces also depend on the muscle parameter model to some degree with the KH model producing lower and higher forces than the CS model in some dynamic cases. As with all models, assumptions are necessary and whether or not the differences in magnitudes that we found are important depends on the question of interest. To demonstrate this point we offer two possible analyses that might utilise these data.

First, our proximate goal is to create an envelope of boundary conditions in a finite element model in order to interrogate the impact of localised morphology, e.g. anteroposterior flattening of the femoral neck in Australopithecus on femoral neck strain, during bipedal walking. Each muscle insertion or origin and the joint contact surfaces of the bone actively load the bone across the stance phase. Each boundary condition constitutes a non-trivial exercise in pre- and post-processing the finite element model. Our initial focus is the hominin femur; our goal, therefore, is to identify those points in the stance phase likely to produce the highest strains in the femur. Consequently, we aim to reduce the number of conditions to an envelope that contains as few conditions as possible to produce the maximum strains in the femur. We specifically seek to identify the points in the stance cycle when hip, knee or patellar joint reaction or femoral muscle forces are maximum within and among people given inherent kinematic, kinetic and muscular variability. While the creation of this envelope is beyond the scope of this paper, it is clear that conditions in early (~20% of stance phase) and late (~80%) stance are needed due to the peaks in the hip and knee joint contact and muscle forces that occur during those gait phases and that conditions near midstance might be less useful. Most biomechanical analyses of the evolution of the morphological features associated with bipedalism have focused on a single stance (Lovejoy, 1988; Ruff & Higgins, 2013), but our current work indicates that more attention needs to be directed towards the loading response (immediately after foot flat) and terminal stance (as the heel loses contact with the substrate) phases of gait. Additionally, our results suggest that it is possible to construct a ‘generic human’ model. The triceps surae functional group, the main source of propulsion in human bipedalism, demonstrates a consistent muscle force pattern across muscle models (as demonstrated here; Figure 1) and among individuals.

A second example of how these data might be used involves the variations in knee joint contact forces (Figure 4) that arise owing to the activation level of the muscle of triceps surae. Our results suggest that researchers interested in knee function in modern humans should focus attention on the influence of muscle parameters on the force distribution among the muscles of the triceps surae complex. Because triceps surae inserts on the calcaneus via the Achilles tendon, the cost algorithm solution for the ankle plantarflexion moment ‘sees’ all of the anatomical muscles of the complex as similar. Soleus does not, however, cross the knee, while gastrocnemius and plantaris do. The impact of triceps surae on lower limb biomechanics is well understood, as is the importance of gastrocnemius activation on knee contact force (e.g. Killen et al., 2021; Mengarelli et al., 2018). Nonetheless, our results demonstrate that knee joint forces increased when the muscle volume of soleus decreased because more propulsive force was distributed to medial and lateral gastrocnemius (which cross the knee), even though the patterns are consistent. Consequently, muscle parameter choices may be important determinants of joint forces and have relevance for such questions as the development of the femoral bicondylar angle in bipeds (Shefelbine, Tardieu, & Carter, 2002) or the differences in shapes of the medial and lateral condyles between australopiths and modern humans that are assumed to be caused by the distribution of forces (e.g. Preuschoft & Tardieu, 1996).

Beyond the utility of this analysis for specific questions, it seems clear more generally that the pattern among and within individuals for functional muscle groups with large magnitudes (e.g. for walking, triceps surae) is consistent when the muscle parameters (muscle volume and optimal fibre length) are varied, which is consistent with our previous work examining kinematic and kinetic variation (Sylvester et al., 2021b). Nonetheless, the magnitude of the forces may vary considerably (15%; Figure 1 and 4). Consequently, comparisons among morphological forms using the same external kinematics and kinetics and muscle model assumptions are preferable to analyses that focus on absolute magnitudes without control of kinematic, kinetic or muscle assumptions.

This simulation uses an algorithm based on muscle strength to resolve the joint moments (that are determined by the external kinematic and kinetic data) into the muscle forces, but this seemingly straightforward relationship is variably reflected in multiarticular muscles. While the relationship between relative maximum force and strength exhibits the shape expected from the algorithm for the triceps surae and hamstring functional groups (Figure 5), that of the quadriceps group is more variable.

A principal limitation of these analyses is that the kinematic and kinetic data derived from different individuals than the muscle parameter model data. Consequently, we do not know if the movements of the dynamic individuals could be achieved by each individual whose data were used to create the muscle parameter models. Most musculoskeletal model simulations use generic muscle parameter models, typically derived from elderly cadavers (e.g. Klein Horsman et al., 2007), so this limitation is not unique to this study. Additionally, we only varied muscle volume and optimal fibre length to create the muscle parameter models. While the Charles et al. (2019) data include pennation angle, we lacked the data to relate the orientations of the AnyBody muscle elements to the Charles et al. full muscles. We also did not address other areas of interindividual variation that could impact joint reaction and muscle forces, including variability in bone shape (e.g. pelvic shape (Cox, 2021), femoral torsion (De Pieri et al., 2021)), muscle origin and insertion locations (e.g. tibialis anterior (Zielinska et al., 2021), popliteus (Berthaume, Barnes, Athwal, & Willinger, 2020)), or even the presence or absence of anatomical elements (e.g. fabella (Berthaume & Bull, 2020), plantaris (Olewnik et al., 2020)). Finally, we chose to use a simplified muscle formulation rather than the more sophisticated Hill-type muscle because the Hill-type muscle formulation requires more detailed data than were available. Given that our approach is more focused on data exploration with the goal of developing an envelope of boundary conditions, this simpler model formulation is reasonable.