Comparison of human gastrocnemius forces predicted by Hill-type muscle models and estimated from ultrasound images

Hill-type muscle models are ubiquitous in biomechanical analyses of human movement and more recently are being applied to animal studies to provide estimates of a muscle's force as a function of its activation state and assumed force–length and force–velocity properties. For example, Hill-type models have been used, together with reconstructions of musculoskeletal geometry, to analyse the moment-generating capacities of animals during terrestrial locomotion (O'Neill et al., 2013; Hutchinson et al., 2015) and to estimate the locomotor capabilities of extinct species based on fossil records (Hutchinson and Garcia, 2002). Muscle models have also been used in simulations of human locomotion to infer the functions of individual muscles during walking (Arnold et al., 2013) and running (Hamner et al., 2010) and to uncover important insights into the biomechanical factors that contribute to gait abnormalities (e.g. Peterson et al., 2010; Steele et al., 2010). However, the accuracy of the muscle forces predicted by Hill-type models, especially during submaximal, dynamic tasks, remains largely unknown.

The predictive accuracy of Hill-type models has been investigated in a small number of animal studies in which tendon forces were directly measured. For example, Sandercock and Heckman (1997) and Perreault et al. (2003) examined whether a Hill-type model could predict in situ forces generated by the cat soleus when they imposed length changes corresponding to those measured during locomotion. Their model reproduced soleus forces measured in vivo, using an implanted tendon buckle transducer, to within 10% of maximal tension during force rise (Sandercock and Heckman, 1997). However, at low motor unit firing rates, differences in the predicted and measured forces were greater than 50% (Perreault et al., 2003). Our own efforts to predict in situ and in vivo forces generated by the gastrocnemius muscles of goats (Wakeling et al., 2012; Lee et al., 2013) confirmed that Hill-type models are sensitive to assumptions about activation state and the force–velocity properties of the model's contractile element. For example, we have shown that Hill-type models reproduce the time-varying in vivo forces with an average r² of 0.40, and errors in force greater than 15% and 28% of maximum isometric force for the medial (MG) and lateral (LG) gastrocnemii of goats, respectively, when averaged across the gait cycle and different locomotor speeds. However, errors in force predictions were reduced when models incorporated the independent activations of slow and fast contractile elements, but only during the fastest locomotor speeds (Lee et al., 2013). Additionally, Millard and co-workers (2013) have demonstrated that Hill-type models are better able to reproduce muscle forces under maximally activated conditions in comparison to submaximally activated conditions. They compared the ability of three different Hill-type formulations to reproduce in vivo forces measured in the cat soleus under maximally active (Krylow and Sandercock, 1997) and submaximally active conditions (Perreault et al., 2003) and found that average errors were 12% of maximum tension for maximal contractions but increased to more than 17% for submaximal contractions. Together, these previous studies and others leave many unanswered questions in regards to the ability of Hill-type models to reproduce time-varying muscle force.

To gain more insight into the strengths and weaknesses of Hill-type models, independent estimates of time-varying muscle forces, for a wide range of movements, are needed. However, muscle forces often cannot be directly measured. We have recently developed an ultrasound-based approach to indirectly estimate the in vivo forces generated by the human triceps surae muscles during dynamic tasks (Dick et al., 2016). In particular, the gastrocnemius–Achilles tendon (AT) complex is advantageous because of its superficial location and short fascicles, which allow us to use ultrasound to measure both tendon and fascicle length changes during movement, and couple these measurements with subject-specific estimates of tendon stiffness and slack lengths to estimate time-varying forces. This is the first study to test Hill-type models against in vivo estimates of time-varying muscle forces from human subjects at a range of mechanical conditions.

Electromyographic (EMG) recordings from the gastrocnemii of cyclists have shown that the activation patterns of different muscle fibre types vary with cadence (Citterio and Agostoni, 1984; Wakeling et al., 2006; Wakeling and Horn, 2009). However, while muscles are composed of a mix of slow and fast fibres with different physiological and biomechanical properties, Hill-type models are typically composed of a single contractile element with force–length and force–velocity properties that have been estimated from in vitro data collected from single fibres under maximally activated conditions. To represent whole muscle, the fibre's force is scaled based on the muscle's level of activation, accounting for the muscle's physiological cross-sectional area, fibre length, pennation angle and parallel elasticity. Here, we tested a traditional Hill-type model, with one contractile element, compared with a differential model, having independent slow and fast contractile elements to account for the contributions of slow and fast muscle fibres. Both models were driven by subject-specific measures of fascicle length, velocity and pennation angle derived from ultrasound images and by activation patterns of slow and fast fibres derived from surface EMG recordings and wavelet-based time–frequency decomposition techniques (von Tscharner, 2000; Wakeling and Horn, 2009). We hypothesized that the two-element model would better reproduce the ultrasound-based estimates of gastrocnemius force at the higher cadences, where preferential recruitment of fast fibres has been reported (Citterio and Agostoni, 1984; Wakeling et al., 2006). This is an extension of our previous studies on the gastrocnemius muscles in goats, where we showed that a two-element model, driven by the independent activation of slow and fast muscle fibres, provided better predictions of time-varying muscle forces than traditional one-element models for some in situ (Wakeling et al., 2012) and in vivo (Lee et al., 2013) conditions.

We collected comprehensive sets of kinematic, kinetic, EMG and ultrasound data from 16 competitive cyclists (8 male, 8 female; Table S1). Subjects pedalled on a stationary cycle ergometer (Indoor Trainer, SRM, Julich, Germany) at cadences between 60 and 140 rpm and crank loads between 14 and 44 N m. In this study, we analysed data from eight different combinations of cadence and load. We predicted the time-varying forces generated by each subject's LG and MG muscles during cycling using a one-element Hill-type model and a two-element Hill-type model (Fig. 1). The models were driven by fascicle lengths, velocities and pennation angles that we determined experimentally from ultrasound images, and by muscle activations calculated from simultaneous recordings of surface EMG (Fig. 1). We compared each subject's predicted LG and MG forces with the forces estimated in vivo from ultrasound-based measures of tendon length changes and stiffness (Fig. 1). Subjects gave informed consent and protocols were approved by the Institutional Ethics Review Boards at Simon Fraser University and Harvard University. The experimental protocol has been described and evaluated in a previous study (Dick et al., 2016), so a brief overview is provided here.

F_max was calculated based on the muscle's scaled volume V_m estimated using the regression equations in Handsfield et al. (2014), and optimal fibre length l_f,opt estimated from a subject-specific scaled musculoskeletal model (OpenSim v3.3; Delp et al., 2007; Arnold et al., 2010; Dick et al., 2016). σ₀ is estimated specific muscle stress (225 kPa; Spector et al., 1980; Roy et al., 1982) (Table S1). c₁ is a scaling factor that accounts for the lower power of the EMG signals that would be expected from action potentials with higher spectral frequencies (Wakeling et al., 2012). c₁ affects the balance between fast and slow fibre contributions and it was given a value of 2 as described below. F_max scales the fibre (or fascicle) force to whole-muscle force and pennation angle β allows the component of fascicle force that is aligned with the tendon to be estimated.

Histological analysis has shown that the human gastrocnemii contain equal cross-sectional areas of slow and fast muscle fibres (Johnson et al., 1973). Therefore, the one-element homogeneous model assumed that the active LG and MG have intrinsic properties that represent the average values for α and between slow and fast fibres (α=0.235; =7.5 s⁻¹) (Wakeling et al., 2012; Lee et al., 2013).

The muscle models (Eqns 1 and 2) were driven directly by the experimentally derived activation, fascicle length and pennation angle, and the predicted muscle forces were subsequently compared with estimates of the tendon force. The models did not involve the balancing of the fascicle force with tendon force (series elastic element force), as is necessary if the whole muscle–tendon unit is incorporated in the model. As such, our approach (which parallels earlier animal studies: Hodson-Tole and Wakeling, 2009; Wakeling et al., 2012; Lee et al., 2013) circumvents the need to apply additional serial damping components to prevent high-frequency oscillations in the solution (Günther et al., 2007; Millard et al., 2013; Haeufle et al., 2014).

Subjects pedalled at eight different combinations of cadence and crank torque: 60 rpm at 44 N m, 80 rpm at 14 N m, 80 rpm at 26 N m, 80 rpm at 35 N m, 80 rpm at 44 N m, 100 rpm at 26 N m, 120 rpm at 13 N m and 140 rpm at 13 N m, corresponding to crank powers of 275, 115, 220, 290, 370, 270, 165 and 200 W, respectively. Each trial lasted 15 s. Each block of eight experimental conditions was repeated, in a random order, four times so that we could separately image both the LG and MG muscle bellies and the LG and MG muscle–tendon junctions (MTJs). Subjects rested for 5 min in between blocks of conditions. ‘Maximum effort’ sprint trials (high power and cadence) were collected at the beginning and end of each session in an effort to elicit maximum muscle activity, and we used these data as a reference when normalizing muscle EMG intensity (e.g. Rouffet and Hautier, 2008). We used an optical motion capture system (Certus Optotrak, NDI, Waterloo, Canada) sampling at 100 Hz to track the 3D locations of 32 LED markers placed bilaterally on the lower extremities, the pedals and the ultrasound probe. Pedal reaction forces effective and ineffective (normal and radial) to the right and left cranks were recorded at 2000 Hz (Powerforce, Radlabor, Freiburg, Germany). Previously, we have verified that the plantarflexion moments estimated from ultrasound-based measures of tendon strain are consistent with the net ankle moments determined from inverse dynamics (Dick et al., 2016).

A linear-array B-mode ultrasound probe (7 MHz, 60 mm field-of-view; Echoblaster, Telemed, Vilnius, Lithuania) recording at 40 Hz was used to alternately image LG and MG fascicles from the right leg during pedalling. The order in which muscles were imaged was randomized; at the end of each block of conditions, the probe was re-positioned over the other muscle. The probe was secured in a custom-made foam support and positioned to image muscle fascicles that were in the middle of the muscle belly (Fig. S1), where the fascicle architecture is homogeneous (Maganaris et al., 1998), and in plane with the scanning image (Kawakami et al., 1993). An ultrasound gel pad (Parker Laboratories, NJ, USA) was placed at the probe–skin interface to enhance image quality and allow muscle bulging.

Ultrasound images were manually digitized (ImageJ, NIH, Bethesda, MD, USA) to determine each subject's time-varying fascicle lengths and pennation angles for the eight pedalling conditions. Eight points in each ultrasound image were digitized – three points on each of the superficial and deep aponeuroses and two points on a representative muscle fascicle (Fig. S2). Pennation angle was calculated as the mean of the angles made by the fascicle with the deep and superficial aponeuroses. Fascicle lengths and pennation angles from 10 complete crank revolutions were fitted using a least-squares minimization of a Fourier series: c₂+a₁sin(Ø₁+ω)+a₂sin(Ø₂+2ω), where c₂, a₁, a₂, Ø₁ and Ø₂ are fitted constants and ω is the crank angle. These fascicle lengths and pennation angles were used as inputs to the muscle models for each muscle, each pedalling condition and each subject.

Bipolar Ag/AgCl surface electrodes (10 mm diameter and 21 mm spacing; Norotrode; Myotronics, Kent, WA, USA) were used to record muscle excitations from the LG, MG, soleus and seven other muscles on the contralateral (left) limb to which the ultrasound transducer was placed. EMG signals were preamplified (gain 1000), band-pass filtered (bandwidth 10–500 Hz; Biovision, Wehrheim, Germany) and sampled at 2000 Hz as described elsewhere (e.g. Wakeling and Horn, 2009; Blake and Wakeling, 2014). Electrodes were placed in the mid-region of the muscle belly after the skin had been shaved and cleaned with isopropyl alcohol solution. We secured electrodes with stretchable adhesive bandages and tubular net bandages to reduce movement artefacts during pedalling.

To compare predicted forces derived from the total intensities (one-element model) with the forces derived from the optimized wavelets (two-element model), the time resolutions from the two approaches were defined to be the same; to do this, we set the Gauss filters at the end of the wavelet analysis to a 50 ms time resolution for all wavelets, which is within the range of physiological time to peak twitch for skeletal muscle (10 and 100 ms; Levin et al., 2000; Spägele et al., 1999). We used the square-root of the EMG intensity as a measure of muscle excitation e, as muscle force is linearly related to the EMG amplitude and not its power (Milner-Brown and Stein, 1975). A lower EMG intensity would be expected from action potentials with higher spectral frequencies (Wakeling et al., 2012). Here, we selected a value of c₁=2 (in Eqn 2) to scale the excitations (based on the ratio of their centre frequencies: Table 2) to account for this effect.

We estimated in vivo muscle forces from subject-specific tendon length changes and stiffness to evaluate the predictions of our Hill-type models. A linear-array B-mode ultrasound probe was secured over the LG or MG distal MTJ. A rigid-body cluster of LED markers (Prager et al., 1998) firmly attached to the ultrasound probe, allowed the 3D coordinates of the LG or MG MTJ to be tracked (ImageJ). We calculated the length changes of each tendon l_AT−LG or l_AT−MG during cycling as the distance from the AT insertion on the calcaneus (identified using an LED marker) to the 3D position of the LG or the MG MTJ. We estimated each muscle's contribution to tendon force, F_LG and F_MG (Eqn 11; shown for MG), during cycling based on the measured length changes l_AT−LG or l_AT−MG, the subject-specific values of toe region stiffness k_SEE–LG,T, k_SEE–MG,T and linear region stiffness k_SEE–LG, k_SEE–MG for either the LG or MG (Table S1), and the tendon slack length l_0,AT–LG or l_0,AT–MG:

Tendon stiffness was estimated from ramped isometric tests (Morrison et al., 2015; Dick et al., 2016), where we quantified the total stiffness of the AT and also determined the proportion of the total AT stiffness contributed by the LG (k_SEE–LG) and MG (k_SEE–MG) separately. These muscle-specific stiffnesses were calculated based on the ratio of the maximum force-generating capacity F_max of either the LG or MG to the combined triceps surae maximum force (Table S1). Each tendon force–length curve had a toe region and a linear region. Linear stiffnesses, k_SEE–LG and k_SEE–MG, were subject specific but due to the difficulty in measuring the toe region properties of the force–length curve in all subjects during the isometric protocol (Dick et al., 2016), the same toe region stiffness, k_SEE–LG,T and k_SEE–MG,T, was used for the LG and the MG in all subjects. F_max,T–MG corresponds to the estimated force at the end of the toe region. Tendon slack lengths for the LG and MG portion of the AT (l_0,AT−LG and l_0,AT−MG) were measured at 310 deg of the crank cycle, averaged over all cycles; this choice was motivated by published tendon buckle data (Gregor et al., 1987) showing force rise to occur after 310 deg. The time-varying tendon force from 10 cycles was used to compute a mean force trace for each condition.

We used the one- and two-element models to predict time-varying LG and MG forces for the 16 subjects and 8 pedalling conditions. Additionally, because the forces predicted by Hill-type models are sensitive to the assumed maximum intrinsic speed ν₀ and curvature α of the force–velocity relationship (Shue et al., 1995; Perreault et al., 2003; Wakeling et al., 2012; Lee et al., 2013), we performed a sensitivity analysis and ran the models, varying the values of α (0.155–0.315) and v₀ (4–10 s⁻¹) to determine the extent to which the modelled forces were sensitive to these force–velocity properties. We compared the predicted forces with the forces estimated from ultrasound-based measures of tendon length changes and stiffness, and characterized differences across the crank cycle using two measures: the coefficient of determination r² and the root mean square error (RMSE). A general linear model ANOVA was conducted to determine whether differences in the r² and RMSE (dependent variable) existed between models, muscles, cycling conditions, subject (random), and values of α and v₀. We used a Chi-square test at each pedalling condition to determine whether the two-element model reproduced the estimated forces in the LG and MG better than the one-element model (higher r² and lower RMSE), over the entire crank cycle, more often than the one-element model. Additionally, an ANCOVA identified changes in the mean muscle activity for total, fast and slow motor units with cadence, crank load, and fascicle length as covariates and subject as a random factor. was included as a covariate because muscle strain has been shown to affect the frequency content of the EMG signal (Doud and Walsh, 1995). Differences were considered significant at the P<0.05 level. Data are reported as means±s.e.

The Hill-type models tested here, driven by EMG-derived activations and ultrasound-based measures of fascicle length, velocity and pennation angle, were able to capture the general features of the ultrasound-based estimates of force (Fig. 3). On average, the models performed with an r² of 0.54 and an RMSE of 13.52% F_max across muscles, conditions and models (Table 4; Fig. S3). However, models performed even better under low cadence–high load conditions, with r² values as high as 0.62 and 0.79 and errors as low as 8.98% F_max and 9.86% F_max for the LG and MG, respectively. A comparison of the modelled forces between the pedalling conditions, between muscles and between the one- and two-element models is shown in Table 4 and Fig. S3.

In both the LG and MG, there were minor differences between the muscle forces predicted by the two-element model when compared with the ultrasound-based estimates of tendon force and those predicted by the traditional one-element model. In particular, when averaged across subjects and conditions, the two-element model predicted forces with r² values of 0.47±0.09 and 0.63±0.1 for the LG and MG, respectively, compared with 0.46±0.1 and 0.61±0.12 for the one-element model LG and MG, respectively. RMSE errors were similar between the two models. The two-element model did perform significantly better than the one-element model at high cadences (100–140 rpm; Table S2), although the differences were small (Table 4; Fig. S3).

Consistent with previous literature (Wakeling et al., 2006; Blake and Wakeling, 2014), the conditions tested in our study elicited a range of activation patterns between slow and fast muscle fibres, providing the rationale to test the two-element model (Fig. 4). When averaged over the whole pedal cycle, the total activation increased with both crank load and cadence in both the LG and MG across all subjects (P<0.05). Activation of the slow motor units decreased with cadence (LG: P=0.02, MG: P=0.045) but increased with crank load (LG and MG: P<0.01). Fast motor unit activation increased with both cadence and load in the LG and MG (P<0.01). Therefore, increased cadence resulted in a reduction of slow motor unit activity coupled with an increase in fast motor unit activity.

The Hill-type models tested in this study reproduced the estimated forces for the MG more accurately than for the LG. Both the one- and two-element models predicted forces with a higher r² over the crank cycle for the MG as compared with the LG across all pedalling conditions (P<0.05) (Table 4; Fig. S3). The average r² across subjects, conditions and models decreased from 0.62 (range: 0.45–0.79) for the MG to 0.47 (range: 0.29–0.62) for the LG. However, the average RMSE over the crank cycle was similar between muscles: 13.55% F_max (range: 9.86–18.25) for the MG and 13.51% F_max (range: 9.98–21.91) for the LG.

Time-varying patterns of the predicted and estimated forces, as assessed by r² and RMSE, showed greater differences at high cadences in comparison with low cadences (Table 4; Fig. S3); this was a common feature for both the one-element and two-element models. In particular, when averaged across the two different models, the RMSE over the crank cycle was 11.1±1.2% F_max at cadences of 60–80 rpm and increased to 17.6±2.1% F_max at cadences of 100–140 rpm when averaged across muscles, subjects and conditions.

The forces predicted by the Hill-type models tested here were sensitive to the force–velocity curvature α used (Fig. 5). In both the LG and MG, the models predicted force with a higher r² and lower RMSE (data not shown) when using α=0.235 or α=0.315 for the one-element model, as compared to α=0.155 when averaged across all pedalling conditions (Fig. 5A). The most notable differences in model performance with α were evident at cadences of ≥100 rpm, where models predicted force with a higher r² and lower RMSE when using greater α values (Fig. 5A). In contrast, the predicted forces were not sensitive to the maximum intrinsic speed v₀ implemented within the models. These trends were similar for the one-element and two-element models tested (data shown for the one-element model; Fig. 5).

Our results reveal that Hill-type muscle models, driven by EMG-derived activation and ultrasound-based measurements of fascicle length, velocity and pennation angle, predicted, on average, 54% the gastrocnemius forces generated by human subjects during cycling, as estimated via ultrasound-based estimates of tendon length changes and stiffness. In general, the models performed best when the fascicles shortened at lower velocities and the activation levels were high, with models under these conditions able to predict nearly 80% of the medial gastrocnemius force over a complete pedal cycle. This work provides the first comparison of forces predicted by Hill-type models – driven with in vivo human experimental data under submaximal dynamic tasks – to ultrasound-based estimates of muscle–tendon force. The two-element model with independent slow and fast contractile elements showed no significant improvement for predicting time-varying and maximum muscle forces, when compared with a traditional one-element model for the slower conditions tested; however, small but significant improvements were noted at the fastest cadences. This suggests that other factors play an important role in determining time-varying patterns and magnitudes of muscle force in addition to the contractile properties of the constituent motor units. Identifying the underlying sources of error in force predictions from Hill-type models remains challenging, particularly when having to rely on non-invasive estimates of muscle–tendon force. However, the approach described here provides an important step and an experimental framework for investigating how Hill-type models may be improved.

The performance of the Hill-type models depended on which muscle was being modelled. The human LG and MG vary in architecture (Maganaris et al., 1998), activation patterns (Wakeling et al., 2011) and motor-unit twitch profiles (Vandervoort and McComas, 1983). Specifically, the MG has shorter, more pennate fascicles than the LG (Wakeling et al., 2006; Ward et al., 2009), while the LG has a more complex architecture than the MG, with multiple heads that have different fascicle arrangements (Wolf et al., 1998). This complexity may have diminished our ability to accurately estimate representative fascicle lengths, velocities and pennation angles in the LG using ultrasound – and suggests that models are likely to be more accurate for muscles with homogeneous architecture. Consequently, these differences may help to explain the variation in model performance shown here, particularly the finding that models reproduced the estimated forces more closely across the pedal cycle for the MG than for the LG, which is consistent with the results from the gastrocnemii of goats (Lee et al., 2013).

The muscle models were compared with estimates of tendon force that were derived from ultrasound measures of tendon length and subject-specific tendon stiffness (methods described in Dick et al., 2016). The tendon stiffness for the LG and MG portions of the AT was previously determined on an isometric frame in which the knee was held at 130 deg and the subjects performed maximal effort plantarflexions (Morrison et al., 2015; Dick et al., 2016). However, it is possible that at different knee angles, where the relative lengths of the LG, MG and soleus are affected differently as a result of the biarticular nature of the LG and MG versus the uniarticular nature of the soleus, we may get a different value of stiffness and thus a different value of force. During cycling, in the middle of the pedal upstroke (270 deg), the AT is unloaded and muscle activation is zero; consequently, both estimated tendon force and predicted muscle force will be zero. Towards the top of the pedal cycle, the AT starts to become loaded; however, as the knee is still relatively flexed, the contribution of the gastrocnemii to the combined AT force is likely to be less relative to the soleus than occurred during our isometric measures of force and tendon strain. This may lead to an over-estimation of the tendon force during this phase, although the forces predicted by the muscle models are independent of tendon stretch and are thus unaffected by this. However, when the knee angle reaches 130 deg at a crank angle of 90 deg, this matches the isometric test conditions and coincides with the period when the muscle forces are the greatest. Additionally, the ankle plantarflexors (LG, MG and soleus) have a similar activity during the low cadence–high load conditions (Dick et al., 2016), which most closely matches the isometric conditions for which tendon stiffness was determined. There was also a greater resolution of ultrasound frames per pedal cycle at the lowest cadences as a result of the finite sampling rate of the ultrasound scanner. Thus, we would expect the best match between the modelled muscle force and the estimated tendon stiffness to occur for the low cadence–high load conditions. Indeed, the r² reached 0.79 with a corresponding RMSE of 10% F_max for these conditions in the MG (Table 4), and this is likely to be the limit to the accuracy that can be achieved using these non-invasive experimental techniques.

The accuracy of the models tested here is partly limited by the measures of fascicle length, pennation angle and tendon length as well as estimates of fascicle and tendon slack length. We have previously shown that small increases in tendon slack length delayed the onset of AT force during cycling and decreased the magnitude of the predicted forces (Dick et al., 2016). Additionally, fascicle slack length is set at the time when the tendon is at its slack length. It is possible that the accuracy of these models could be further improved if the fascicle and tendon slack lengths were independently determined by optimization (e.g. Gerus et al., 2012; Gerus et al., 2015), rather than being experimentally guided; however, this was not the purpose of this study. Additionally, the overall accuracy of the muscle model is scaled by F_max (Eqn 3), that was determined from both an optimal fibre length and an estimated muscle volume that were derived from regression equations from other studies (Delp et al., 2007; Handsfield et al., 2014).

Previous studies have assessed the accuracy of Hill-type models against direct measures of force from tendon buckle recordings in animals (Sandercock and Heckman, 1997; Perreault et al., 2003; Wakeling et al., 2012; Lee et al., 2013; Millard et al., 2013; Kim et al., 2015) or against measures of heat and work (Umberger et al., 2003; Lichtwark and Wilson, 2005). In particular, modelling in situ forces yielded higher r² values (Sandercock and Heckman, 1997; Perreault et al., 2003; Wakeling et al., 2012) than in vivo forces (Lee et al., 2013), which probably relates to the more controlled contractions studied during in situ experiments in comparison to in vivo experiments. The errors achieved in this study for the slower cadences (RMSE of 10% F_max with r² of 0.79) match the best in vivo models (RMSE of 10%; Lee et al., 2013) that have been achieved in animal studies where the data were collected from more invasive techniques (fine-wire EMG for activation, sonomicrometry for fascicle length, and tendon buckles for force).

Differences between the forces predicted by the one- and two-element models in this study were smaller than the differences predicted by similar one- and two-element models in the goat gastrocnemii (Lee et al., 2013). However, both the goat study and the current study did find the only significant improvements in the two-element model to occur at the fastest muscle strain rates. Surprisingly, the greater benefits from the two-element model occurred in the goat study, which did not show such fast strains rates as seen here. However, the differences in model performance may reflect the additional challenges of extracting EMG information from surface EMG signals in humans, or from the effect whereby contractile differences between muscle fibre types may be mitigated in larger muscles that are submaximally activated (Holt et al., 2014; Ross and Wakeling, 2016).

Previous studies have shown that the errors in forces predicted by Hill-type models are often related to the assumed maximum intrinsic speed v₀ and curvature α of the force–velocity relationship (Shue et al., 1995; Perreault et al., 2003; Wakeling et al., 2012; Lee et al., 2013). In contrast to our previous study (Lee et al., 2013), predicted forces of the models here were not sensitive to maximum intrinsic speed v₀ (Fig. 5). However, predicted forces were sensitive, in a cadence-specific manner, to the curvature α of the force–velocity relationship. Most notably, models reproduced the estimated forces more closely (higher r² and lower RMSE over the pedal cycle) using low α values at low cadences, but models performed better using higher α values at high cadences (Fig. 5A). For a given shortening velocity, a higher α would result in a larger force. As mentioned previously, slow fibres are characterized by a higher curvature (lower α) whereas fast fibres are characterized by a lower curvature (higher α) (Otten, 1987; Umberger et al., 2003). These results suggest that the forces predicted by traditional Hill-type models may be improved by incorporating a force–velocity curvature α that varies appropriately with the type of muscle fibres recruited for the specific mechanical demands of the task – a larger α for high-speed tasks where faster fibres are recruited, and a smaller α for low-speed tasks where slower fibres are recruited. Umberger and colleagues (2003) have previously alluded to this, and have shown that Hill-type models with force–velocity parameters, to include v₀ and α, that are matched to the proportions of slow and fast muscle fibres, generate reasonable predictions of whole-muscle energetics during isolated muscle contractions, single joint motion and whole-body movement.

Consistent with previous cycling studies (Wakeling et al., 2006; Wakeling and Horn, 2009; Blake and Wakeling, 2014), higher cadences elicited differential recruitment patterns that displayed an increase in the activation of the faster motor units coupled with a decrease in the activation of the slower motor units (Fig. 4). However, the differences between predicted forces for the one- and two-element models at even these highest cadences were minor (average 3%). This could be attributed to the fact that (i) the one-element model was not sensitive to ⁠, (ii) the differences in recruitment across the range of loads and cadences tested here may not have been large enough to elicit substantial differences in model performance, or (iii) the forces from human-sized muscle at submaximal activation may be less sensitive to differences in fibre-type recruitment than previously expected. To date, our understanding of whole-muscle behaviour has largely come from studies in which muscles are maximally activated, and from single-fibre data or data from smaller muscles (<10 g). However, recent evidence suggests that muscle behaviour also depends on the muscle's size (Ross and Wakeling, 2016) and level of activation (Rack and Westbury, 1969; Josephson and Edman, 1988; Rassier et al., 1999; Holt and Azizi, 2014, 2016; Holt et al., 2014; Ross and Wakeling, 2016). Understanding the role of different fibre types within large mixed muscles during voluntary submaximal contractions is clearly an avenue that warrants further investigation.

In this study, we compared gastrocnemius forces predicted by Hill-type models with the forces estimated from ultrasound-based measures of tendon length change and stiffness. Similar forces were predicted for the traditional one-element model and the two-element model that accounted for the independent contributions of different muscle fibre types across most of the speeds tested, and the models performed best when fascicle velocities were low and muscle activation was high. What additional factors should we consider in order to yield more accurate muscle models? Future work should aim to additionally consider the mechanical effects of a muscle's physical properties (e.g. mass, viscosity), as well as the potential contributions of elastic tissue (e.g. aponeurosis) and surrounding structures on muscle contractile function. This is likely to be critical for improving the reliability of muscle models as well as for understanding the mechanisms underlying whole-muscle function during dynamic sub-maximal contractions – which we currently know little about. Benchmark experimental data sets, like the data we have analysed here, are necessary and useful for testing and improving models.

The authors thank Dr Allison Arnold for insightful discussions.