Comparison of the validity of Hill and Huxley muscle–tendon complex models using experimental data obtained from rat m. soleus in situ

All human locomotor tasks (e.g. walking, cycling, rowing and swimming) involve periodic movements. Both mechanical behaviour and metabolic energy expenditure are important variables in the study of such movements. A musculoskeletal model that can accurately predict both is therefore an essential tool in biomechanics research. Currently, the Hill muscle model (Hill, 1938) is most used in musculoskeletal modelling (e.g. Delp et al., 2007; Lee et al., 2013; Lai et al., 2014; Biewener et al., 2014). Musculoskeletal models using the Hill model have been successfully applied to study phenomena in which only mechanical behaviour is considered (e.g. Van Soest et al., 1993; Bobbert, 2012; Azizi and Roberts, 2010; Lai et al., 2014). However, in the Hill model, a direct relationship between metabolic and mechanical behaviour is absent. Attempts have been made to establish such a relationship (Anderson and Pandy, 2001; Bhargava et al., 2004; Lichtwark and Wilson, 2005a), but these involve adding extra, phenomenological relationships to the model, which introduces extra parameters, the values of which are not straightforward to determine. Moreover, it is unclear whether metabolic cost models based on heat terms obtained in stylized experiments can be generalized to real-life stretch–shortening contractions (SSCs). This problem may (partially) explain the lack of success in predicting whole-body metabolic energy expenditure in human walking using such musculoskeletal models (Anderson and Pandy, 2001). An alternative to the Hill muscle model is the Huxley model (Huxley, 1957), which models the dynamics of cross-bridge cycling and, as such, does entail a direct relationship between metabolic and mechanical behaviour. The original Huxley model and multi-state variations thereof have been widely studied, but mainly to address fast time-scale phenomena on the single fibre level (e.g. muscle's transient response to a quick release or stretch), which could not be explained by the Hill model (Huxley and Simmons, 1971; Eisenberg and Hill, 1978; Tözeren, 1985; Campbell, 2009; Stoecker et al., 2009). In explaining the latter phenomena, series elasticity and/or submaximal activation are not modelled. However, both tendons (Alexander, 2002; Roberts, 2002; Lichtwark and Barclay, 2010; Lai et al., 2014) and the time course of muscle activation (Van Soest and Casius, 2000; Lichtwark and Wilson, 2005b) play a crucial role in the functioning of human muscle–tendon complexes (MTC) in vivo. Thus, if the Huxley model is to be used in musculoskeletal modelling, it should be implemented as an MTC model, which includes modelling of series elastic structures and activation dynamics. Although several attempts have been made (Comincioli et al., 1984; Julian, 1969; Williams, 2011; Stoecker et al., 2009), the latter is mathematically complicated, and difficult to formulate in a generic way, such that the MTC model can be incorporated in a larger musculoskeletal model. Moreover, the Huxley model was previously considered to be computationally too demanding for use in musculoskeletal modelling (Van den Bogert et al., 1998; Winters and Stark, 1987; Zahalak, 1981; Cadova et al., 2012). Zahalak (1981) mitigated the computational issues by introducing the simplifying assumption that the distribution of the fraction of attached cross-bridges always follows a Gaussian curve: the distribution moment (DM) approximation. Many authors aiming to describe macroscopic muscle behaviour using the Huxley model have adopted the DM approach (Cadova et al., 2012; Winters and Stark, 1987; Cholewicki and McGill, 1995; Gielen et al., 2000; El Makssoud et al., 2011). Notwithstanding the important contribution the DM approach has made to the usage of cross-bridge models in MTC modelling, it is problematic that, in many applications, it is impossible to establish to what extent the assumption on the distribution affects the predicted muscle behaviour (Wu et al., 1997; Tözeren, 1985). Moreover, the computational issues mentioned above are no longer relevant at present. In the present study, the original two-state Huxley model was extended to include both activation dynamics and parallel/series elasticity, and was formulated as a generic state-space model. As an important step towards implementation of the Huxley model in musculoskeletal modelling, the aim of the present study was to assess the validity of the latter model, with respect to mechanical behaviour. To evaluate the relative merit of the Huxley MTC model in the context of large-scale musculoskeletal modelling, its validity was compared with that of the Hill model. Although widely used in musculoskeletal modelling, there are only a few studies in which the Hill model is validated for realistic muscle contractions and in which parameters were estimated from data independent of the data used to validate the model (Krylow and Sandercock, 1997; Sandercock and Heckman, 1997; Williams et al., 1998; Lee et al., 2013). Therefore, as recently pointed out by Biewener et al. (2014), there is an ongoing need for validation of the Hill model. As such, a secondary aim of the present study was to contribute to previous work in which the Hill model was validated for realistic muscle contractions. To these ends, data were collected from rat m. soleus (SOL) in situ. The experimental data were divided into two parts. Data from the first part were used for the estimation of parameter values for both models. Data from the second part, consisting of realistic SSCs, were used to evaluate model performance, by comparing the experimental data with model predictions thereof.

where τ_act and τ_deact are the activation and deactivation time constants, respectively, and the superimposed dot represents differentiation with respect to time.

if q≤q_crit, with the minimal value of the scale factor at q=q_min. To ensure maximum CE contraction velocity is not affected by CE length at CE lengths above optimum, a^rel was scaled by the isometric force from Eqn 4 if ⁠.

The eccentric part of the Hill force–velocity curve is modelled as a hyperbola with an oblique asymptote. The hyperbola is defined by the ratio between the concentric and eccentric slope in the isometric point (r_slope), and the intercept of the asymptote with the ordinate (F_asymp). The slope of the asymptote is equal to the slope of the concentric part at q=q_min and ⁠.

where is the first-order moment of the steady-state solution of Eqn 10, evaluated at u=0.

where g₂ and g₃ are the detachment rate parameters. Note that g(x)=0 in the region where cross-bridges can attach [i.e. Huxley's g₁=0; our notation thus agrees with Zahalak (1981)], which is normally not the case because this would result in zero energy expenditure during steady-state isometric contraction. However, in the present study only mechanical behaviour is considered and it was found that the force–velocity curve can be adequately modelled while constraining g₁ to be zero (see Fig. 5D). Allowing the value of g₁ to be determined in the optimization aimed at describing the force–velocity data was found not to improve the quality of the fit. We emphasize that, if data on energy consumption were available, g₁ would have been made a tunable parameter, used to fit the additional data regarding metabolic energy expenditure. This would have resulted in values for the other rate parameters that are different from those reported here.

The stepwise nature of Eqns 13 and 14 required a small integration step size during simulations. Therefore, to improve numerical tractability, the rate functions were smoothed using a sine function in an interval of 0.1 h around the transition points, resulting in rate functions with a continuous first derivative, as depicted in Fig. 2. It was found that the smoothening did not materially change model behaviour.

Note that with [n l_CE γ]^T as a quasi state, combinations of n and l_CE can occur that violate Eqn 2. Thus, the initial conditions must be chosen such that Eqn 2 is satisfied at t=0, and the violation of Eqn 2 must be monitored throughout the simulation. In all simulation results presented here, the violation of Eqn 2 remained below 10⁻⁸ N.

All experiments, animal handling and procedures were approved by the animal ethics committee of the Vrije Universiteit Amsterdam (protocol number FBW 12-03). Data were collected from male Wistar rat SOL (n=3, see Table 1 for rat and muscle mass). We chose SOL for this study because it is a parallel-fibred, homogeneous, fatigue-resistant muscle, which means ignoring these respective properties in our modelling efforts will have minimal influence on our results. After anaesthetization, which consisted of two doses of 6 ml urethane per 1 kg body weight, administered 15 min apart, the hind body of the rat was shaved bare and the skin was opened and partly removed. Subsequently, SOL was exposed by dissecting surrounding muscles. Care was taken to maximally separate SOL from all its surrounding tissue, except for the origin and the neurovascular tract, which passes through the m. gastrocnemius lateralis. Specifically, the n. plantaris was removed, and all connections to the m. plantaris and m. gastrocnemius medialis were severed. The complete Achilles tendon was dissected free and at its insertion a small part of the calcaneus was cut and left attached to the tendon, leaving the complete SOL MTC intact. The operation technique proved crucial for the success of the experiments, as exemplified in Fig. S1. After dissection of the SOL, the distal end of the Achilles tendon was tightly knotted to a steel rod, using the leftover part of the calcaneus. The femur was made accessible (blunt dissection) and attached to a screw terminal. Finally, a cuff electrode was placed over the n. ischiadicus and connected to a direct current stimulus isolator (Digitimer Limited, Hertfordshire, UK, model DS3). Hereafter, the rat was fixated on the experiment table by securing the screw terminal and fixating the foot. The steel rod was connected to a servo-controlled, linear motor (custom made, VU University, Amsterdam, The Netherlands), via a force transducer (miniature strain gauge transducer, Honeywell International, Morrice Township, NJ, USA). In Fig. 3, an overview of the setup is depicted. The resolution of the entire setup was 10 mN and 1.0 μm for force measurement and length manipulation, respectively. All data were collected at 1000 Hz. The laboratory was climate-controlled and kept at a constant temperature of 24°C and humidity of 80%.

At the start of each trial, the MTC was brought to its desired initial length by a manually operated micromanipulator (analogue resolution of position control 0.1 mm). After each trial, the muscle was shortened below slack length, by the same means. During all measurements, the muscle was kept moist by applying saline after each trial. To prevent fatigue, time between trials was at least 2 min. First, supramaximal current amplitude was determined at a length close to MTC optimum length, by increasing the amplitude until peak force during a tetanic stimulation no longer increased. During the experiment, maximal isometric force was monitored to ensure the chosen current amplitude still resulted in supramaximal stimulation. Supramaximal current amplitude was approximately 0.5 mA in all experiments. Stimulation frequency and pulse width were 100 Hz and 100 μs, respectively, in all experiments. MTC optimum length was found iteratively in a series of trials, by applying tetanic stimulation at increasing MTC length (1 mm per step) until maximum force decreased. Next, the muscle was alternately stimulated isometrically at optimum length and at short length, until forces at both lengths stabilized.

For part 1 of the experiment, the applied contraction protocols were based on assumptions following the Hill MTC model. To collect data pertaining to the force–length relationships of SEE, PEE and CE, 8 to 10 contractions with a quick-release protocol were applied, similar to Blümel et al. (2012), at different MTC lengths ranging from active slack length to 1 mm over optimum MTC length, in 1 mm steps. The quick-release protocol started with an isometric phase until force levelled off, followed by a rapid 0.1 mm MTC length decrease. Stimulation ceased shortly thereafter (Fig. 4A). Next, to collect data pertaining to the concentric part of the CE force–velocity relationship, a series (15–20) of step-ramp contractions were applied, similar to Blümel et al. (2012) and Curtin et al. (1998). Each step-ramp started with an isometric phase at MTC optimum length, followed by a step decrease in length and immediately thereafter a ramp shortening (Fig. 4B), at a velocity intended to result in constant force during the ramp. The protocol was applied at different combinations of step sizes and ramp slopes, such that the forces measured during the ramp covered the force range from close to 0 to maximal isometric force. To collect data pertaining to the eccentric part of the CE force–velocity relationship, a series (four to six) of isokinetic, eccentric contractions were applied. Each contraction started isometrically at MTC optimum length followed by an eccentric ramp lengthening (Fig. 4C). Note that the eccentric contraction protocol was applied (contrary to the order in which it is listed here) at the very end of the experiment, to minimize the influence of any adverse effects of the lengthening contractions. As we did not have access to biochemical measures in the present study, the parameters pertaining to the activation dynamics were estimated by fitting them to behavioural data, similar to Curtin et al. (1998). These data included several trials of the previous protocols and a single isokinetic release trial, which consisted of an isometric force build-up, followed by isovelocity shortening, at a velocity close to the maximum MTC velocity encountered during the SSCs of part 2 of the experiment (see below). Stimulation was ceased at onset of the shortening phase. In part 2 of the experiment, SSCs were imposed on the muscle. Each SSC consisted of five quasi-sinusoidal movement cycles (frequency 1 Hz, amplitude 1 mm) in which stimulation was applied in the middle three cycles, during shortening (Fig. 6, top inset; note that in the figure only the final two cycles during which stimulation was present are shown). The stimulation timing was chosen such that net positive mechanical work performed by the muscle was large. After the measurements, SOL was dissected, weighed and stored in fixation fluid (75% tridest, 4% formaldehyde, 15% ethanol, water and 1.7 g l⁻¹ thymol). After removal of SOL, the animal was euthanised by injection of an overdose of pentobarbital.

Data analysis was based on assumptions pertaining to the Hill version of the model depicted in Fig. 1. For both models, parameters were intentionally estimated as much as possible from behavioural data, to allow for a comparison between models that was not compromised by different levels of uncertainty in each model’s parameter values. Table 2 shows the values of all parameters that were not fitted to the data, and lists their respective literature sources.

After was found, SEE slack length was determined as: ⁠. Subsequently, PEE slack length was determined as: ⁠. Parameters pertaining to the SEE, PEE and CE force–length relationships were assumed to be equal for both models.

Data from the step-ramp and isokinetic eccentric trials were treated similarly. First, l_CE(t) and F_CE(t) were determined by solving Eqns 1–3 with the now known parameters on the experimental force and length traces. Next, the 20 consecutive samples of the force trace for which the sum of squared difference from their average was smallest were selected from the set of samples (after the step) for which stimulation was present and MTC velocity was constant. After normalization, the average CE force and the slope of the best-fit straight line through l_CE(t) of these 20 samples provided one data point in the force–velocity curve (Fig. 5D). For the Hill model, the inverse of Eqn 7, with ⁠, was fitted to these data (i.e. the difference between model and data CE force was minimized in the fitting procedure).

For the treatment of the data from the quick-release trials it is assumed that the step happens instantaneously. In practice, this is obviously not the case. The violation of this assumption was corrected for by first computing the CE displacement during the step as the product of the velocity corresponding to the average force level (ΔF) and the duration of the step (∼6 ms). Next, the estimated CE displacement was subtracted from the length change (Δl) and the analysis of the quick-release trials was repeated, yielding corrected values for the parameters pertaining to the force–length relationships. The latter parameters were also used in the treatment of the trials pertaining to the force–velocity relationship, which could thus also be re-analysed. The above procedure was iterated until the difference between the parameter estimates in consecutive loops was smaller than some arbitrary small value. For all rats, the iteration procedure required no more than 15 iterations, and yielded marginally changed parameter values.

The Huxley model rate parameters were estimated using the previously found force–velocity data points. Huxley model CE force was computed from the steady-state solution to Eqn 12, evaluated with at each CE velocity. The rate parameters were then found by minimizing the difference between CE model and data force.

For both models separately, the time constants in Eqn 6 were estimated by fitting simulation results to four different, experimentally obtained force traces, simultaneously. The four trials were: (1) a quick-release trial around optimum length (Fig. 4A); (2) a step-ramp trial with the constant shortening velocity similar to the maximum MTC shortening velocity encountered during part 2 of the experiment (Fig. 4B); (3) an eccentric isokinetic trial, at a velocity where both the Hill and the Huxley model accurately fitted the corresponding data point in the force–velocity relationship (Fig. 4C); and (4) an isokinetic release trial, which consisted of an isometric phase followed by isokinetic shortening with stimulation ending at the onset of the shortening phase, at a shortening velocity at or around maximum MTC shortening velocity encountered during part 2 of the experiment (Fig. 4D). The latter four trials represent a wide range of contraction conditions. We therefore expected the parameter set resulting from this fit to generalize well to other contraction types.

The statements made in this section hold for all three rats, unless specified otherwise. The muscles remained metabolically stable throughout the experiment. Maximum isometric force at optimum MTC length at the end of the experiment did not differ by more than 10% from that at the beginning of the experiment.

The contraction protocols from part 1 of the experiment resulted in measured force traces (Fig. 4A–C) from which data pertaining to both the force–length relationships and the concentric and eccentric force–velocity relationships could be accurately extracted. Fig. 5A and B show typical examples (rat 3) of the fit of Eqn 3 to, respectively, the SEE and PEE data obtained from the quick-release trials. For both relationships, the model described the data well (PEE: R²>0.97; SEE: R²>0.99). The data from the CE force–length relationship show a steep ascending limb, which was followed by a plateau and a decrease in force at CE lengths larger than optimum CE length (Fig. 5C). The proposed model for the CE force–length relationship (Eqn 4) described these data well (R²>0.99). Data pertaining to the concentric part of the force–velocity curve followed a typical hyperbolic shape (Fig. 5D), which was described well by both the Hill and Huxley models. Fewer data points were collected on the eccentric part of the force–velocity curve, and these data showed a more irregular pattern for each rat, which was described similarly well by both models. For the force–velocity relationship, we found R²>0.98 for both models.

Results of simulations of data trials used to estimate the force–length relationships of CE, SEE and PEE matched the observed data closely (Fig. 4A), particularly those parts from which the data in Fig. 5A–C were derived. This is also the case for the force–velocity relationship, for both the Hill and the Huxley models in equal measure (Fig. 4B,C). For both models, the fit of the activation dynamics parameters yielded parameter values with which all four contractions could be adequately described (Fig. 4), resulting in low root mean squared (RMS) differences between model predictions of force and experimental data (Table 3, part 1).

For both models, all parameter values estimated from the experimental data obtained in part 1 of the experiment are provided in Table 4. Note that the parameter values of rat 1 and 3 are similar, whilst the values of rat 2 deviate from the latter two rats. Note also that the only between-model difference in parameter values concerns the activation dynamics parameters, as the parameters pertaining to the CE, SEE and PEE force–length characteristics are assumed to be equal for both models.

The stretch–shortening protocol applied during part 2 of the experiment resulted in considerable force production during shortening (Fig. 6). The observed force traces were similar for the three SSCs within a trial (difference between lowest and highest observed peak force, averaged over animals <0.03 N). As can be appreciated from Fig. 6, both models yielded similar predictions of the experimental data. During simulations with both models, CE started contracting eccentrically while significant force was still present, as indicated by the crosses in the force traces in Fig. 6. Predicted peak force was lower for both models in equal measure. Forces predicted by both models deviated from experimental data during relaxation, both during the concentric and the eccentric phase, in equal measure. Table 3, part 2, lists the RMS differences between the measured force traces and force traces pertaining to simulations with the Hill and Huxley models, corresponding to the data in Fig. 6. Absolute RMS differences between model predictions and experimental data were similar to those found for the trials used to fit the activation dynamics time constants (Table 3, part 1). As can be appreciated from Table 3, the difference between simulated and measured force traces was similar for both models. Furthermore, the difference between the two models was smaller than the difference between each of the models and the data. In all simulation results of SSCs, CE length varied between and ⁠, and CE velocity varied between and ⁠. Following the similarity in force traces (given equal SEE characteristics), simulated CE length and velocity traces were similar for both models.

The aim of this study was to evaluate and compare the validity of the Huxley and Hill MTC models, with respect to mechanical behaviour, during SSCs in which the muscle was primarily active during MTC shortening. Data were collected from rat SOL in situ. Importantly, the data used to fit the parameters of both models (part 1) were not used to evaluate model performance (part 2). Furthermore, to minimize bias in the comparison between the Hill and the Huxley MTC models, parameters were estimated from the same behavioural data as much as possible, for both models. For both models, the fitting of the model parameters on data from part 1 of the experiment was exceptionally good (Figs 4, 5). Comparison of simulation results with data from part 2 of the experiment (Fig. 6) shows that both models predict experimental data reasonably well during realistic SSCs. Moreover, comparison between the simulation results of the Hill and the Huxley models (Figs 4, 6) indicates that both models behave similarly under a wide range of contraction conditions.

The number of animals included in this study is modest (n=3). However, given the exploratory nature of this work, and the similarity of the results between rats, we deemed the number of animals included sufficient.

The between-animal variability in the parameter values warrants attention. Specifically, rat 2 differs from rats 1 and 3 (Table 4). These differences reflect the progressive interdependency in the parameter estimation procedure, i.e. the parameters estimated after the force–length relationships are influenced by the previously estimated parameter values. Therefore, the relatively small differences between rats in parameters pertaining to the force–length relationships are amplified in the estimates of the force–velocity parameter values for both models. Any differences between model behaviour and data are subsequently ‘corrected’ for in the estimates of the activation dynamics parameters. Therefore, although the behaviour during the SSCs was similar for all three rats, estimated parameter values differed. Note, however, that it is not uncommon to find large variety in animal-specific parameter estimates, as illustrated by Blümel et al. (2013). Thus, the large variability in parameter values found here probably reflects a combination of (small) between-animal biological variance and amplification of differences arising from the estimation procedure. In short, we feel that the methodological issues discussed above do not jeopardize the validity of the results and do not compromise the possibility to draw conclusions regarding the key questions of this study.

The results of the present experiments regarding the model parameter estimation (part 1) were similar to those previously reported by Blümel et al. (2012) (stick insect), Curtin et al. (1998), Williams et al. (1998) (lamprey) and Krylow and Sandercock (1997). The experimental data obtained in the present study agree with data obtained in the latter studies. The fitting of parameters to these data was successful for both models in equal measure. For the Hill model, the quality of the fits was similar to that presented by Blümel et al. (2012). The Hill model parameters presented in the present study agree with previously reported values in rat SOL, near optimum length (see Table 1 in Krylow and Sandercock, 1997). This suggests that the measurements and procedures for estimating the parameters pertaining to the force–length and force–velocity curves employed in the present study are reliable and valid.

For the activation dynamics, we obtained the parameter values pertaining to the steady-state relationship between free calcium concentration and active state (Eqn 5) from literature. The time constants of the (de)activation dynamics were fitted on four representative trials (Fig. 4), an approach similar to that of Curtin et al. (1998). For reasons mentioned previously, these values differed between rats. Note that for all rats, both the activation and deactivation time constants were smaller for the Huxley model than for the Hill model (Table 4). This is not surprising: the Huxley model contains an additional dynamical step in the transformation of stimulation to force as compared with the Hill model. During an isometric contraction of the CE, a stepwise change in active state would result in a stepwise change in force in the Hill model, but lead to a dynamic response in the Huxley model, as it takes time for the cross-bridges to attach or detach and for the distribution n(x,t) to settle on its new steady-state value (see also Materials and methods section). Thus to arrive at a given force level within a given time limit, the activation dynamics time constants must be smaller in the Huxley model compared with the Hill model.

An important point to note here is that the estimated parameter values pertaining to the force–length and force–velocity relationships were consistent for both models. The parameter estimation procedure applied to simulations of the experimental data yields similar parameter values, particularly as those parts of the force traces from which the data were extracted to form the force–length and force–velocity relationships were adequately simulated (Fig. 4). For the Hill model, this was not surprising, as the assumptions on which the data analysis was based were tailored to the Hill model. However, for the Huxley model this was not obvious, as cross-bridge force at any time depends not only on the state at that time, but also on its history. Yet, a prominent feature of the simulation results presented in Figs 4 and 6 is the similarity between the model predictions of the Hill and Huxley models, over a wide range of contraction conditions. It might be that this similarity arises because the dynamics of cross-bridge cycling are so fast that the distribution n(x,t) is similar to the steady-state distribution at the instantaneous velocity, at all time points during a simulation. If this were the case, the macroscopic behaviour of the Huxley model would effectively be the same as that of the Hill model, given the similarity in the shape of the Hill and Huxley (steady-state) force–velocity curves (Fig. 5D). To investigate this possibility, the distribution n(x) during a simulation of an SSC was compared with the steady-state distribution at the same instantaneous velocity, activation level and relative CE length, at four different time points (Fig. 7). As can be appreciated from Fig. 7, the two distributions differed at t=0.1 s and t=0.34 s and were similar at t=0.22 s and t=0.46 s. Thus, the similarity between the behaviour of the Hill and Huxley models is not trivially caused by a large difference in time scales between microscopic and macroscopic behaviour of the Huxley model, which would render modelling the cross-bridge distribution superfluous.

The experimental data of the SSCs obtained in the present study agree well with those from similar experiments performed previously on mammalian muscle (Lichtwark and Barclay, 2010). The level of agreement between Hill model predictions and experimental data found in the present study agrees with previous studies in which a similar approach was adopted (Biewener et al., 2014; Williams et al., 1998; Krylow and Sandercock, 1997; Sandercock and Heckman, 1997). As model predictions were similar for both models (see Table 3, Fig. 6), we will discuss the results for both models simultaneously. For both models, predictions of the SSCs deviated from the data both during relaxation and, to a lesser degree, during stimulation (Fig. 6). Note that in Table 3, the absolute RMS differences between part 1 and part 2 are not large, but the RMS differences relative to the peak force attained during the trial are larger for the simulations of part 2. Below, we discuss possible causes for this discrepancy between model predictions and experimental data.

With respect to the activation phase, the main difference between the trials on which the activation parameters were fitted and the SSCs was that during the former, activation always occurred isometrically (on the MTC level), whilst during the latter it occurred during shortening of the MTC. This observation suggests that the concentric part of the force–velocity relationship is dependent on the time course of activation, in a way that is not modelled here.

The discrepancy during relaxation can possibly be attributed to inadequate extrapolation of the eccentric part of the force–velocity curve, of which we have collected only limited data. During relaxation, MTC velocity was close to zero whilst muscle force was decreasing, thus SEE was shortening and CE was probably contracting eccentrically. Although the latter was indeed the case, the force was already deviating well before the CE started contracting eccentrically (see crosses in Fig. 6). Moreover, during the portion of the experiment in which significant force was present, eccentric velocities did not exceed the range of values for which experimental data were collected. Therefore, in our view, inadequate extrapolation of the eccentric force–velocity curve cannot explain the discrepancy between experimental data and model predictions during relaxation.

There are two well-known phenomena in muscle physiology that were not modelled in the present study and that should also be considered as possible explanations of the discrepancy mentioned above. These are: (1) muscle inhomogeneities, stemming from either inhomogeneous fibre type distributions or serial sarcomere strain inhomogeneity (Palmer et al., 2011) and (2) contraction history effects [e.g. potentiation (Krarup, 1981) and shortening-induced force depression (Meijer et al., 1997; Rassier and Herzog, 2004)]. With respect to the first option, Lee et al. (2013) have recently shown that explicit modelling of fibre type distribution can improve Hill model predictions of goat m. gastrocnemius mechanical behaviour in vivo. However, rat SOL consists for the most part of slow-twitch, type I fibres (∼85%, Delp and Duan, 1996; Pousson et al., 1991). It is therefore unlikely that lumping all fibre types has had a large influence on our results. Serial sarcomere strain inhomogeneity has been shown experimentally (Palmer et al., 2011), and was also pointed out by Curtin et al. (1998) to explain the discrepancy between model predictions and experimental data during the relaxation phase of SSCs. In a Huxley modelling framework, it has been shown that modelling sarcomere strain inhomogeneity can account for the stretch response of muscle (Campbell, 2009) and can result in emergent behaviour that cannot be explained by a lumped fibre model (Stoecker et al., 2009). Whilst not ruling out the possibility of such a mechanism being at work here, we have no grounds to assume this mechanism would work differently during the contractions in part 1 as compared with part 2, and consequently we have no grounds to assume that this mechanism explains the discrepancy between experimental data and simulation results during relaxation. With respect to the second option, it is noted that the force traces pertaining to each sinusoid were nearly equal, which indicates that effects of history dependence were minimal on the time scale of the sinusoid's period. However, within one sinusoid, the muscle is passively stretched prior to active shortening. This stretch may have altered the excitability of calcium-gated potassium channels (Mallouk and Allard, 2000), caused the release of extra calcium (Ji et al., 2002) or altered cross-bridge contraction dynamics (Haugen, 1991); these processes were not modelled in the present study. However, if these processes were to result in increased force during the relaxation phase of the shortening contraction, an increased force at maximum length should also be expected. As the latter was not observed, we consider it unlikely that any of these processes underlies the observed discrepancies.

Another possible explanation may be related to the duration of stimulation: >0.5 s during part 1 versus <0.25 s during part 2. The duration of stimulation may affect the calcium saturation characteristic (Eqn 5). This could lead to over-saturation of the model during part 2, which may explain the overestimation of the force during the last part of relaxation. Although, in this study we did not fit the parameters in Eqn 5 on the data in part 1, but based them on Curtin et al. (1998), it is still possible that these parameters resulted in more adequate prediction of the data in part 1, as compared with part 2. In short, we believe that the discrepancy between model predictions and experimental data is caused by inadequate modelling of physiological effects that arise due to the differences in the contraction protocols between part 1 as compared with part 2.

For the Hill model, the level of agreement between predicted and experimental data is not perfect, which is in line with previous work (Lee et al., 2013; Krylow and Sandercock, 1997; Williams et al., 1998; Sandercock and Heckman, 1997). This suggests that caution is necessary when applying the model in situations in which accurate quantification of muscle force is required (e.g. in patient-specific modelling). However, Hill models have been successfully applied in musculoskeletal models aimed at addressing conceptual questions regarding design and control of the motor system (Van Soest et al., 1993; Bobbert, 2012; Azizi and Roberts, 2010; Richards and Clemente, 2012; Kistemaker et al., 2007; Lichtwark and Wilson, 2007). Given the similarity in predicted behaviour of the Hill and Huxley MTC models in the present study, it stands to reason that the Huxley model could at least be applied in a similar fashion. The level of agreement between predicted and experimental data displayed in the present study is deemed adequate for the latter type of applications.

The present study is the first to rigorously validate a Huxley MTC model and compare it with the Hill MTC model in the context of physiologically relevant SSCs. It is concluded that, with respect to mechanical behaviour, the current Huxley MTC model predicts the experimental data as well as the Hill MTC model. As such, large-scale musculoskeletal modelling using the current Huxley MTC model is a viable option. It remains to be determined whether the Huxley model can adequately describe both mechanical behaviour and metabolic energy consumption from a unified framework. In future research we will investigate the validity of the present Huxley model in terms of metabolic energy consumption and its relationship with mechanical behaviour.

The authors thank Dr Alistair Vardy for helpful discussions regarding the formulation of the Huxley model.