Estimation of the force–velocity properties of individual muscles from measurement of the combined plantarflexor properties

The force-generating properties of individual skeletal muscles have been widely studied to understand how muscles enable a variety of motor functions. Furthermore, the properties of the multiple muscles spanning the same joint synergistically combine to determine the moment-generating properties for a particular movement direction. If the individual muscles within a synergistic group have contractile properties that are identical, differ only in size (e.g. physiological cross-sectional area, PCSA) and share a common tendon, the combined properties of the muscle group are a simple weighted sum of the individual muscles. However, in the case of the ankle plantarflexors, substantial differences have been reported among the architectural properties of the individual muscles (i.e. fiber type, pennation angle, optimum length and shortening velocity) (Brown et al., 1996; Close, 1972; Cui et al., 2008) and mechanical properties of their tendons (Cui et al., 2009; Matson et al., 2012). Moreover, rotation of pennation angle as force changes can decrease a muscle's output force but increase output velocity by allowing the muscle to function at a higher gear ratio (Azizi et al., 2008). This implies that the contributions of individual muscles to the total force and power of the group may change with loading conditions. In general, it is not well understood how differences in the properties of individual muscles influence the properties of a group of muscles acting synergistically at a joint (Biewener and Roberts, 2000; Rajagopal et al., 2015).

This gap in knowledge is especially apparent for dynamic conditions, when the force–velocity (F–V) properties determine the amount of power that can be delivered over a range of velocities and forces. The relationship between the properties of individual muscles and of the muscle group cannot be calculated unless the properties of individual muscles and the muscle group are measured simultaneously. This is because the alternative of physically separating muscles and measuring their properties is complicated by intermuscular force transmission and/or aponeuroses shared by adjacent muscles (Maas and Sandercock, 2010; Maas et al., 2004; Rijkelijkhuizen et al., 2005; Tijs et al., 2014, 2015, 2016). In vivo approaches, such as sonomicrometry or x-ray cinematography, have been used to measure shortening velocity and, in some cases, force in animals but these studies are technically challenging and only provide insight into muscle behavior during specific tasks (Arellano et al., 2016; Biewener and Roberts, 2000; Moo et al., 2017). In human studies, non-invasive in vivo experiments have also been used to characterize the moment–velocity relationship at a joint, but often do not break down the contributions of individual muscles (Hill, 1938; Huijing, 1996). Modeling and optimization methods have been used in combination with in vivo measurements of the torque–angle and torque–angular velocity of the ankle joint to estimate the properties of individual plantarflexor muscles (Hasson and Caldwell, 2012; Hasson et al., 2011); however, these involve multiple assumptions concerning the model parameters of the muscles and their tendons, many of which cannot be validated experimentally. The limitations of these approaches underscore the need for detailed experimental studies that simultaneously characterize the properties of the group and of the individual muscles to elucidate how individual muscles contribute to the group properties (de Brito Fontana et al., 2018; Epstein and Herzog, 1998; Herzog, 2017).

To demonstrate how the force–length (F–L) properties of individual muscles combine to form the group properties, we previously studied the plantarflexors of kangaroo rats under static conditions (Javidi et al., 2019a). Briefly, the kangaroo rat plantarflexors are an excellent model for addressing how muscles work together at a single joint for two reasons. First, they deliver the greatest contribution to the power of a vertical jump (Schwaner et al., 2018), and thus they all are likely synergistically maximally activated during jumping. Second, the soleus (SOL) is the only uniarticular plantarflexor and is very small in kangaroo rats (about 2% of total plantarflexor mass; Biewener and Blickhan, 1988; Biewener et al., 1988; Rankin et al., 2018), so it can be removed without greatly affecting plantarflexor capacity. The remaining plantarflexors, the lateral and medial gastrocnemii and plantaris (LG, MG and PL, respectively), have origins almost at the same place on the femur and are likely to have predominantly fast-type muscle fibers typical of small rodents (Schiaffino and Reggiani, 2011), and their distal tendons attach to approximately the same point on the calcaneus (Rankin et al., 2018). It is also very important to note that in kangaroo rats, the tendons of the LG and MG (i.e. the gastrocnemii, GAS) merge together but the tendon of the PL can be easily separated from the other two; thus, the PL tendon is mechanically independent of the LG and MG tendons (Javidi et al., 2019b).

In this study, our objective was to estimate the F–V relationships of the individual muscles of the plantarflexors of kangaroo rats from the measured F–V relationship of plantarflexors as a group, which was characterized by the standard in situ technique of applying isotonic loads. To perform this estimation, two different models based on two different assumptions about the contributions of individual muscles to group output were used. In the first model, we assumed that each muscle's percent contribution (based upon cross-sectional area) to the group F–V relationship would be constant across all group forces. The second model assumed that the F–V relationships for the fibers of each individual muscle were identical, which would change each muscle's percent contribution according to the group force (due to differences in architectural parameters amongst muscles). In addition, comparisons of the two model predictions were made to independently derived estimates of force levels of individual muscles based upon the length changes in tendons, which were calculated from sonomicrometry data during the isotonic experiments. By comparing the model predictions with an independent estimate, we could test the hypothesis that inclusion of appropriate fiber properties and muscle architecture is necessary to separate the individual muscle contributions across different force levels.

All procedures completed were approved by the Institutional Animal Care and Use Committee (IACUC) at Washington State University. Experiments were performed with seven wild-caught adult kangaroo rats, Dipodomys deserti Stephens 1887 (4 males, 3 females, mean±s.d. mass 102.9±22.58 g, age unknown).

The surgical procedure and experimental setup were similar to our previous study on the F–L properties of kangaroo rat plantarflexors (Javidi et al., 2019a); thus, they will be briefly described here. Animals were anesthetized with inhalant isoflurane (1%) and an intraperitoneal injection of 50 mg kg⁻¹ ketamine and 0.15 mg kg⁻¹ dexmedetomidine. The plantarflexor muscles of the right leg with their tendons were dissected from surrounding muscles and tissues, leaving the proximal end attached to the femur, major blood vessels intact, and the distal tendons attached to the calcaneus, which was later cut to create a bone chip attached to the tendon. The SOL muscle was removed from the rest of the plantarflexors because it was very small relative to the other plantarflexors (see Introduction). The sciatic nerve was cut proximally (near the hip joint) and distally, the tibial branch was left intact and the other branches (i.e. sural and common peroneal) were cut. To sonometrically measure the muscle belly length changes of the LG, MG and PL, a pair of piezoelectric crystals (1.0 mm diameter; Sonometrics) were implanted near the proximal and distal aponeurosis of each muscle of 5 animals (see Fig. 1A, white circles on LG as an example). The servomotor (309C; Aurora Scientific, Aurora, ON, Canada) measured both force and position at 1000 Hz (Real-Time Muscle Data Acquisition, Aurora Scientific, Aurora, ON, Canada) and was connected to the Achilles tendon via a hook clipped onto the bone chip.

To mechanically ground the origin of the muscle group, a stainless steel intramedullary pin (0.7 mm o.d.) was inserted through the femur and its ends were then secured to a custom-fabricated support arm attached to the animal platform. A wire loop (0.5 mm o.d.) was wrapped around the middle of the femur and support arm to prevent any horizontal displacement of the femur. The muscle group was immersed in a custom-made temperature-controlled (37°C) trough of mineral oil for the duration of testing to prevent muscle dehydration and maintain temperature. The calcaneal bone chip was then attached to the lever-arm of the servomotor and the sciatic nerve was placed onto custom-made bipolar hook electrodes (0.3 mm silver wire). A rectal temperature probe and pulse oximeter (PhysioSuite; Kent Scientific, Torrington, CT, USA) placed on the paw allowed for continuous monitoring of core body temperature, oxygenation level and heart rate.

Servomotor force [equivalent to musculotendon (MT) force, F_GRP], servomotor arm position relative to its initial position (equivalent to the change in musculotendon length, ΔL­_MT) and distance between the pairs of sonometric crystals inserted into individual muscles (L_c) (Fig. 1) were recorded at 1000 Hz.

Before the active muscle measurements, the motor arm was fixed vertically (ΔL­_MT=0), and the motor was moved until the passive force of the muscle group was equal to a very small positive value. The distance between the servomotor center and the origin of the muscle group (L_ser in Fig. 1A) was then measured by calipers and defined to be the slack length of the muscle–tendon unit (L­_MT,s). Distances between the sonometric crystal pairs of individual muscles also were collected at this time. We used the slack length as the initial guess of optimal MT length (L­_MT,o*).

We determined the actual optimum MT length of the group (L­_MT,o) so that the isotonic measurement could be made near optimal length to minimize any F–L effects. To do so, the isometric force of the muscle group was measured at a minimum of seven different random MT lengths within a range of L­_MT=L­_MT,o*+ΔL­_MT, where −4≤ΔL­_MT≤+3 mm, using supramaximal isometric tetanic contractions generated by 300 ms trains of 0.2 ms biphasic pulses at 150 Hz. We fitted a polynomial curve to the active force and length data, and the optimum MT length (L­_MT,o) and the maximum isometric force of the group (F_GRP,max) were estimated from the values associated with the peak of the fitted curve (more detailed methods are found in Javidi et al., 2019a; Rehwaldt et al., 2017). Subsequently, all isotonic trials were started at L­_MT,o (Fig. 1).

Isotonic quick-releases were performed with a minimum of 6 isotonic levels, with 9 levels being typical and not exceeded. These levels were relative to maximal tetanic force (F_GRP,max) and within a range of 0.1–1.0 F_GRP,max. We could not use a level of zero force because overshoot in the servomotor caused an artifact as the MT went slack. For each level, the MT was stretched to L­_MT,o and we elicited an isometric supramaximal tetanic contraction with a 250 ms train of 0.2 ms biphasic pulses at 150 Hz (Fig. 1). The servomotor then switched from length control to force control and isotonic tetanic shortening occurred for a 150 ms interval. Isotonic levels were randomized and muscles were rested for 3 min between each level. At the end of the experiment, animals were terminally anesthetized, muscles were dissected from their tendon and each muscle mass was measured.

The objective of the data analysis was to estimate the F–V curve by calculating the velocity and force following the application of the isotonic load, when it is assumed that velocity and force have reached a steady state. For each isotonic level, ΔL­_MT was fitted with a line from 10 to 25 ms following the force drop (Fig. 1C) to estimate shortening velocity of the musculotendon (V_MT). The isotonic force of the muscle group (F_GRP,isot) was also calculated for the same time interval (Fig. 1B).

To pool the muscle group F–V and power–velocity (P–V) relationships across animals, V_MT was normalized by V_MT,max and F_GRP,isot was normalized by F_GRP,max for each animal's data. The pooled normalized F–V data were fitted to the rectangular hyperbolic curve (Eqn 1). Note that in Eqn 1, F_GRP,max and V_MT were replaced by 1 and normalized velocities, respectively. The normalized P–V relationship was calculated as the product of the normalized forces and velocities.

The normalized F_GRP–V_MT relationship for the data pooled across animals was then used to estimate the F–V curves of the LG, MG and PL. Briefly, we used two sets of assumptions that resulted in two different models, which we called Mdl1 and Mdl2. An overview of how we used the two models is shown in Fig. 2 and described in detail below.

where F_M and V_M are the isotonic load and shortening velocity of the individual muscles (i.e. PL, MG or LG).

In Mdl1, we assumed that each muscle's percent contribution (based on PCSA; see below) to the group F–V relationship would be constant across all group forces (Fig. 2). In other words, the normalized F–V relationships of individual muscles were assumed to be the same as the normalized F–V relationship of the muscle group. In Mdl2, we assumed that the normalized F–V relationships of the fibers within each individual muscle were identical, and this would result in a change in each muscle's percent contribution according to the group force (as a result of differences in architectural parameters among muscles). For example, fiber length scaled the normalized F–V relationship, so that differences in fiber length between muscles resulted in differences in predicted force for a given absolute velocity. Because the kangaroo rat plantarflexors are likely all fast fibers (Graffe et al., 2017), we used a value of 0.29 for the curvature of the F–V curve (a_f), which is appropriate for fast fibers based upon a literature survey of 59 species (Dick et al., 2017). To find the normalized (by initial fiber length, l_f) maximum fiber shortening velocity (V_f,n-max), we optimized V_f,n-max, such that the summation of forces of the individual muscles, which were calculated by using the F–V equation (Eqn 1) scaled by fiber length and maximum isometric force (see below), were equal to the group force for all velocities (Fig. 2). The guess of V_f,n-max was then optimized until model calculation of the group force matched the experimental measurement of group force.

In both models, the velocity of individual muscles was assumed to be equal to that of the muscle group. The force of individual muscles at each level of isotonic load or constant shortening velocity of the muscle group also was estimated by models. These estimations of the force with the velocity can simply be used to calculate percent contributions of individual muscles to the muscle group force and power.

The general basis for the calculation of F–V properties based on sonomicrometry data was using the amount of tendon stretch and recoil during the isotonic protocol to estimate force in each muscle. Specifically, similar to the models explained above, we assumed that the shortening velocities of individual muscles were equal to the shortening velocity of the muscle group. Then, for each isotonic trial, the length of MT and distances between crystal pairs inside individual muscles were calculated in three conditions (see Fig. 1). The initial condition was when the muscle group was at its slack length, and the distances between the pairs of sonometric crystals inserted into the individual muscles were calculated as L_c,s. The maximum isometric contraction condition included two different parts: the passive stretch of the muscle group to its optimum length (Fig. 1, brown solid lines) and then the isometric contraction (Fig. 1, green solid lines). The optimum musculotendon length (L_MT,o), and the individual crystal distances after passive stretch (L_c,p) and after isometric contraction (L_c,o) were measured in the isometric condition. The last condition was immediately after the tendon recoil phase (Fig. 1C,D, red solid lines), when the MT length (L_MT,isot) and crystal distance of individual pairs (L_c,isot) were measured.

where ΔL_MT,s-o and ΔL_MT,o-isot are the changes in MT length from the slack length to the optimum length of the group (Fig. 1C, green solid line) and from the optimum length to the isotonic contraction (ΔL_MT,o-isot was equal to the amount of change in the red solid line, Fig. 1C), respectively. ΔL_c,s-p, ΔL_c,p-o and ΔL_c,o-isot are the changes in distance between the pair of sonometric crystals of individual muscles from slack length to passive stretch, shortening during isometric contraction and from isometric to the start of isotonic contractions, respectively. Note that because this method calculates the isotonic force as force normalized by the isometric force, it does not require knowledge of the specific value of tendon stiffness (i.e. tendon stiffness cancels in the numerator and denominator of the formulation of Eqn 8). To separate the contributions of the LG and the MG, we also had to assume that their tendons were effectively independent, and as a consequence, the muscle length changes (as measured by sonomicrometry) were influenced by the apparent stiffness of the tendon attached to each muscle.

Sonomicrometry was performed on 5 animals. The normalized forces of individual muscles (F_M,n) and the normalized shortening velocities of the muscle group (V_MT,n) for 5 animals were pooled and fitted to a rectangular hyperbolic curve (Eqn 1).

The error metrics were calculated as the absolute values of the difference between each point from sonomicrometry data, as well as between each point from sonomicrometry data and the F–V relationships predicted by Mdl1 and Mdl2. Then, the mean errors were used in an ANOVA table to test whether the method to estimate muscle force (Mdl1, Mdl2 or sonomicrometry) had a significant effect on the predictions of individual muscle properties. If the effect was significant (P<0.05), the methods were compared with the Tukey–Kramer test using the multicompare function of MATLAB to find which paired comparison had a significant difference at a level of 0.05.

Results (presented as means±s.d.) are the data from 7 animals. The mass of the PL, MG and LG muscles was 0.35±0.09, 0.48±0.12 and 0.48±0.10 g, respectively. The calculated PCSA of the PL, MG and LG muscles was 34.30±9.90, 30.54±8.01 and 27.14±5.88 mm², respectively.

The muscle group F–V and the P–V curves of each animal were fitted to Hill's equation to visualize the variability across all animals (Fig. 3A,B). The means of F_GRP,max and V_MT,max were 25.63±5.83 N and 146.54±41.97 mm s⁻¹, respectively. By using the maximum isometric force of the group and its PCSA, the maximum stress of the kangaroo rat plantarflexors combined was 283±37.8 kPa. The mean maximum power (P_GRP,max) and optimum shortening velocity of the MT (V_MT,o) were 335.02±89.49 mW and 43.54±8.95 mm s⁻¹, respectively. Normalized F–V and P–V data from all animals were pooled and fitted to Hill's equation (Fig. 3C,D). The maximum power of the muscle group normalized to the weight of plantarflexors was 255±24.86 W kg⁻¹. The curvature (a/F_GRP,max) of the fitted normalized F–V curve was 0.24.

The normalized force and fiber F–V relationships and the normalized P–V relationships of individual muscles were estimated by the two models (Fig. 4). The differences between the estimates of the individual muscle properties by the different models, specifically for PL and LG are clear (Fig. 4A,C). The properties of individual muscles estimated by the two different models are given in Table 1. By using the calculated PCSAs and Eqn 6, the maximum isometric force of PL, MG and LG was estimated as 37.3%, 33.2% and 29.5% of F_GRP,max.

Noticeable differences were observed between the contributions of individual muscles to the group force and power at different isotonic levels predicted by the two models (Fig. 5). Specifically, at lower isotonic levels (or higher MT velocities), Mdl2 predicted a decrease in the contribution of PL to the muscle group force (or power), while for Mdl1, the contribution remained constant (Fig. 5E,F). Similarly, there was an increase in contributions by the MG and LG as isotonic level decreased for Mdl2, while they remained constant for Mdl1.

The forces of individual muscles normalized by their maximum isometric forces at different isotonic levels of the muscle group were estimated from the sonomicrometry data and fitted to Hill's equation (Fig. 6). The statistical analysis indicated that the models have significant effects on the predictions of the isotonic load of individual muscles (P=0.001). The individual group comparison (Table 2) indicated that the prediction of Mdl1 was significantly different from the fitted curve to sonomicrometry-based data (P=0.0007) while there was no significant difference between the Mdl2 prediction and sonomicrometry-based calculation (P=0.4690).

The motivation of this study was to address the challenges involved with estimating the contributions of individual muscles within a synergistic group, especially under in vivo conditions when it is only possible to measure the mechanical output from the group. By comparing model predictions with independent estimates of the F–V curves, we tested the hypothesis that the measured properties of the kangaroo rat plantarflexor muscle group could be used to estimate the F–V and P–V relationships of individual muscles, if muscle-specific architecture is incorporated into a model.

There was variation in the absolute maximum shortening velocity and maximum power of the plantarflexor muscles as a group for 7 animals (Fig. 3A,B). The animals used in this study were different in sex, size and likely age, which was unknown because they were wild caught. Any of these factors could have effects on the performance of the muscles and the absolute values of the maximum isometric force, maximum shortening velocity and maximum power. However, the variation decreased by normalizing the values (Fig. 3C,D).

Moreover, variation in the F–V curves likely arose from the difficulty in collecting data at low isotonic levels without artefacts due to oscillations from the servomotor (see Materials and Methods). By not including data from the lowest isotonic forces, we could have underestimated the maximum shortening velocities obtained from curve fitting (Marsh and Bennett, 1986). Practically speaking, we preferred to estimate the curve shape and its maximum based on the more reliable forces that were collected at a higher level of the forces. Furthermore, to minimize the effects of the transition from isometric to isotonic load and to maintain consistency in our analysis, we measured the velocity within the period of 10–25 ms for all trials, which may have led to an underestimation of the isotonic velocity.

Our estimate of muscle group maximum stress σ_GRP,max=283±37.8 kPa is very similar to the value we measured in the study of the F–L properties of kangaroo rat plantarflexor (Javidi et al., 2019a). This value was slightly larger than specific tensions from in situ experiments on the SOL, PL and GAS of several species, such as mouse, rat, rabbit, cat and hopping mouse, which were measured in the range from 150 to 250 kPa (with a few exceptions) (Rospars and Meyer-Vernet, 2016). This discrepancy may be due to the calculation of specific tension by muscle cross-sectional area (without the correction for pennation angle, Eqn 7). The specific tension values without accounting for pennation angle were σ_GRP,max=263±35.1 kPa.

The values reported for maximum intrinsic speed (V_f,n-max) of mouse, cat and rat muscle measured at physiological temperatures have a large amount of variability in that they range from 4.8 to 7.3 s⁻¹ for slow fibers and from 9.2 to 24.2 s⁻¹ for fast fibers (Dick et al., 2017; Medler, 2002). The values of the maximum shortening velocity of the individual muscles, estimated by either Mdl1 or Mdl2 (Table 1), were consistent with the fast fiber type data of similarly sized species. This conclusion is also supported by a preliminary electrophoresis study which showed the LG was composed entirely of fast myosin heavy chain isoforms (Graffe et al., 2017).

Pennation angle has been shown to change with length and load (Azizi et al., 2008), and this affects the calculation of fiber velocity, specifically when converting muscle velocity to fiber velocity. We did not factor in this effect in the current study, and we used a constant pennation angle at all levels of force. An increasing pennation angle at smaller isotonic loads would have increased the estimation of fiber velocity (for a given muscle velocity; see Eqn 3) and produced a larger estimate of maximum velocity, V_max.

The maximum power output of vertebrate skeletal muscles has been reported to be in the range from 150 to 480 W kg⁻¹ (Askew and Marsh, 1998; Nelson et al., 2004; Rehwaldt et al., 2017; Sawicki et al., 2015; Swoap et al., 1997; Weis-Fogh and Alexander, 1977). Previous studies of the jumping performance by kangaroo rats (Biewener and Blickhan, 1988; Schwaner et al., 2018) estimated maximum power in individual muscles to be within this range of reported values and concluded that kangaroo rat muscles are not comparatively more powerful. The maximum power that was calculated here by the two different models either for the muscle group or for the PL, LG and MG individually (Table 1) also did not exceed this range.

We estimated the F–V properties of each muscle from a measurement of the group F–V properties using two models, with each model having a different set of assumptions. Mdl1 assumed that the percent contribution of each muscle to the group F–V does not change with isotonic level. For example, this assumption has been used to measure the contribution of vastus lateralis muscle to moments generated by isokinetic contraction of knee extensors (Ichinose et al., 2000). Although this approach is not theoretically correct (de Brito Fontana et al., 2018; Herzog, 2017), we wanted to see how much the predication would be different from Mdl2 and the estimation based on sonomicrometry measurements. Mdl2 assumed that the fiber F–V properties were identical, but because of muscle architectural differences, the muscles would operate on different portions of the fiber F–V curves for different isotonic levels. As a result, the percent contribution of each muscle to the group force would change as the isotonic level (and shortening velocity) changed. In fact, the PL does not contribute any force at the lowest isotonic levels because the high group muscle velocity results in a fiber velocity greater than the PL's V_max. The assumptions of Mdl2 are commonly used in musculoskeletal models to estimate F–V relationships of individual muscles from the torque angular velocity relationship at the joint (Anderson et al., 2007; Caldwell and Chapman, 1991; Garner and Pandy, 2003). Although these methods also can be combined with imaging techniques to estimate more realistic model parameters (Fukashiro et al., 1995; Hasson and Caldwell, 2012; Rugg et al., 1990), it is not possible to measure all musculoskeletal parameters experimentally. Like these other studies, we could not measure all the muscle parameters directly in that we needed to optimize the parameter of fiber maximal shortening velocity, and this produced uncertainty in the model predictions.

We made the assumption that in both Mdl1 and Mdl2, GAS and PL muscle–tendon units are two independent actuators with the same mechanical actions in the sagittal plane, thus the total force of the group is equal to the sum of the force generated by individual muscles and transferred through their tendons (see Eqn 4). However, force transmission between muscle–tendon units through the connective tissue due to the relative displacement of the muscles (Maas and Sandercock, 2010; Tijs et al., 2016) or tendon interactions (Tijs et al., 2014) can occur. In this study, the impact of these two factors is likely to be minimal because the PL and GAS muscles both have minimal shortening during the recoil (see Fig. 1D) and then shorten with similar relative displacements during the isotonic phase, and in kangaroo rats, the PL and GAS tendons are completely mechanically separate (Javidi et al., 2019b) and thus do not interact (see Introduction). To separate the contributions of the LG and the MG, we also had to assume that their tendons were effectively independent, and as a consequence, the muscle length changes (as measured by sonomicrometry) were influenced by the apparent stiffness of the tendon attached to each muscle.

Because the two models provided different predictions for how the percent contribution of each muscle varies with isotonic level, we could test the validity of each model if we had an independent estimate of individual muscle force. The sonomicrometry measurements provided that independent estimate, but this methodology also had several key assumptions that were simplifications of the physiology involved. The main assumption was that the tendon had a constant stiffness throughout the range of forces. The non-linearity of tendon stiffness, specifically in the toe region of the stress–strain curve of the tendon properties, could affect the method of estimation of individual muscle contribution to the group muscle force. Specifically, the average stiffness value of the isometric phase (from passive to maximal isometric force) is likely to be lower than the stiffness value during the isotonic recoil phase (from maximal isometric force to a lower isotonic force), leading to an underestimation of the normalized isotonic force (see Eqns A1–A3). This would potentially shift the F–V curves (e.g. as shown in Fig. 3) upwards slightly, causing a larger estimate of V_max.

The statistical analysis indicated that the estimations of the properties by Mdl1 and Mdl2 were significantly different, and there were no significant differences between Mdl2 estimations of individual muscle properties and the sonomicrometry-based estimations (Table 2). The fact that the Mdl2 predictions agreed with the independently derived estimates of individual muscle F–V properties does support the use of MT models with the appropriate architecture parameters to estimate the individual muscle properties from measured group properties.

The second model (Mdl2), which was the more accurate predictor of the individual muscle properties, predicted that an individual muscle's contribution changes with loading condition (Fig. 5F). The change in contribution is consistent with the effect of differences in muscle architecture, namely fiber length and pennation angle, which determined fiber velocity for a given muscle group velocity (Bodine et al., 1982; Segal et al., 1986; Spector et al., 1980; Wakahara et al., 2013).

The results presented here clearly indicate that the assumption of a constant contribution of a muscle to the F–V or P–V relationship throughout the entire range of isotonic loads or shortening velocities is not correct. Consistent with our results, in vivo measurements combined with optimization methods to estimate the individual muscle properties of the human plantarflexor (Hasson and Caldwell, 2012) indicated changes in the contribution of GAS and SOL muscles with the level of muscle group isotonic load. However, that study could not validate the results by independently obtained data as we did in this study, with the caveat that the methodology used to estimate muscle forces from the sonomicrometry data had several simplifying assumptions.

More generally, our results indicated that the contribution of the individual muscles to the group muscle force is not distributed entirely based on the PCSA of the muscles. This is consistent with the findings in the human wrist extensors, the extensor carpi radialis brevis and longus (ECRB and ECRL, respectively). The ECRB produces more extensor moment than the ECRL under isometric conditions, but the reverse is true at velocities nearer to maximal shortening velocity of the combined muscle group (Lieber et al., 1997). In kangaroo rats, PCSA of the PL is higher than that of the LG and MG muscles (Rankin et al., 2018). Our previous study (Javidi et al., 2019a) and this study show that the percent contribution of the PL can be lower than that of the MG and LG muscles at certain muscle group lengths or isotonic levels, respectively.