A modified Hill muscle model that predicts muscle power output and efficiency during sinusoidal length changes

The relationship between the power output of a muscle and the energetic cost of achieving this power output is critical to the locomotory potential of an animal. Power output of a muscle is modulated by changing activation parameters during cyclical length changes. These include the timing of activation (phase of activation) and the period of activation (duty cycle).

It is generally assumed that, under sub-maximal conditions, muscle activation patterns are optimised to achieve maximum efficiency of work. It has been shown in a range of experiments that both the power output and efficiency of a muscle depend on the frequency of oscillation, length change,duty cycle and phase of activation(Barclay, 1994; Curtin and Woledge, 1996; Ettema, 1996). These studies have demonstrated that a muscle can produce power at a range of efficiencies. For instance, it has been shown that activating the muscle for a longer fraction of the total stretch shortening cycle tends to increase the power output of a muscle, but decrease the efficiency(Curtin and Woledge, 1996). This was due to the excess heat produced during the stretch of muscle. A relatively broad range of activation conditions and length change trajectories would achieve near optimal power output and optimal efficiency, but undertaking sufficient measurements to map these conditions is difficult experimentally.

The reason why the activation conditions for optimum power and optimum efficiency are different is poorly understood. However the series elastic element (SEE) must be accounted for when trying to understand muscle power output and efficiency. The SEE is critical as it can act as an energy storing mechanism, where energy stored during stretching of the SEE can be recovered later in the contraction (Alexander,2002; Biewener and Roberts,2000; Fukunaga et al.,2001; Roberts,2002). This means that the time course of the power output of the contractile element (CE) and of the muscle–tendon unit (MTU) as a whole can differ during a contraction. It has been suggested that this series elasticity makes muscles more versatile under varying locomotor conditions. For instance, when a muscle accelerates an inertial load from rest, early in the movement the CE contraction velocity is higher than that of the MTU because the SEE is stretching; later in the movement the MTU velocity is higher than the CE velocity because the SEE is shortening(Galantis and Woledge, 2003). This should, theoretically, enable the CE to operate at a velocity concomitant with optimum efficiency or optimum power for more of the movement.

Previously the force and power output of muscle have been accurately predicted during contractions with brief tetani during sinusoidal length changes (Curtin et al., 1998; Woledge, 1998). Cost of contraction can also be derived from Hill-type muscle models that incorporate the SEE (Anderson and Pandy,2001; Ettema,2001; Umberger et al.,2003). This is achieved by fitting curves over experimentally derived relationships between energetic cost and power output during contraction. An appropriately validated model of this type makes it possible to explore and map the relationships between power and efficiency of muscle with varying duty cycle, phase of activation and frequency of oscillation,which is difficult to do experimentally. In this paper we: (1) adapt the model used by Curtin et al. (1998) to predict both the power output and the cost of contracting muscles during sinusoidal length changes, (2) validate the model's predictions of muscle energy expenditure (heat + work) by comparing the output of our model to data of force output and heat expenditure during sinusoidal length changes with brief tetani from dogfish Scyliorhinus canicula white muscle and mouse Mus domesticus red muscle (soleus) and (3) determine whether the model could account for the differences between optimum power and optimum efficiency conditions by comparison of the resultant power output and efficiency of these muscle types under experimental conditions to the model predictions.

In the original model of Curtin et al.(1998), a block stimulation was applied, such that during the time period of a train of stimuli pulses the muscle activation level increased to a maximum of 1.0 with an exponential time constant of rise and fall (see Fig. 2A). This stimulation level can be taken to represent the concentration of free calcium (Ca²⁺) available to bind to troponin. This stimulation level is in turn related to the activation level, which represents the relative number of attached crossbridges (Act). This relationship is also shown in Fig. 2A and is described by Curtin et al.(1998).

Raw data from the results reported by Curtin and Woledge(1996) were compared to the predicted force and energy outputs (heat + work) with respect to time. A subset of varying duty cycles, stimulus phases and frequencies (0.71, 1.25 and 5 Hz) of sinusoidal MTU length changes were chosen to compare to the model. The activation parameters (τ₁, τ₂, Kand n) were first optimised to minimise the sum of the force differences between the model and the experimental results at each time point for one individual condition (frequency=1.25, duty cycle=0.121, stimulus phase=3; from the raw data of Curtin and Woledge, 1996) using the Nelder–Mead simplex (direct search)method (Matlab, Mathworks Inc, Natick, MA, USA). The length change, force,energetic output and activation level were then compared between the model and the experimental results. An estimate of the activation level was made from the experimental data using Eq. 4, and an estimate of the CE velocity was made from Eq. 12. These were input into the energetic model along with the experimentally determined force output to approximate energetic output.

The activation parameters that control the uptake of the activator(K and n) were then optimised to fit force output for each of the individual cycle frequencies (Table 1). This was done to account for possible variation in these activation parameters due to shortening speed – a phenomenon known as shortening deactivation. The activation parameters (K and n)were optimised to minimise the sum of the force differences between the model and the experimental results at each time point for three individual conditions at each cycle frequency, and the average of these taken as activation parameters for each cycle frequency. The resultant relationship between the activator concentration and activation level for the optimised activation parameters for each condition is shown in Fig. 5. A comparison between the model with the optimised activation parameters and the experimental results for force output, activation level, energetic output and muscle fibre velocity is shown in Fig. 6.

Similarly, the model was fitted with the appropriate values for the mouse soleus muscle and the protocol employed experimentally by Barclay(1994) was simulated with the model. This was to ensure that the model could emulate muscle from other species and had not merely been fitted to reproduce the dogfish muscle data. The only parameters that were changed were V_max (scaled to L_o s^–1) and the activation constants(τ₁, τ₂, K and n) and the force data was scaled to F_max(Barclay, 1994); these parameters are muscle specific. The appropriate optimal cycle phases, duty cycles and stimulations were applied to the model across a range of frequencies. In this model, it was assumed that the series elastic component of mouse soleus muscle had a similar relative stiffness to that of the dogfish muscle.

The activation parameters (τ₁, τ₂, Kand n) were first optimised to minimise the sum of the force differences between the model and the experimental results at each time point for the following condition: frequency=1.25, duty cycle=0.121, stimulus phase=5. A comparison of the model's response to this stimulus is compared to the experimental results in Fig. 3. The values of the activation parameters that achieved this fit are given in Table 1. It is apparent from the time course of the predicted activation level of the muscle fibre bundle that the activation parameters used in the model provide a good representation of the activation level.

The force–time and energy–time outputs from the model were then compared to experimental results from a single muscle fibre bundle preparation(as used in the results of Curtin and Woledge, 1996) across a range of duty cycles, stimulus phases and oscillating frequencies. The results of the simulations are compared to the experimental results in Fig. 4.

Force development within the model was shown to match that of the experimental examples at the frequency for which the activation parameters had been optimised (1.25 Hz). At the fastest frequency (5 Hz), the model develops force at a faster rate than the experimental data suggests. The model also predicts a significant increase in force output as the muscle begins to lengthen, the effect of which is less evident in the experimental results. The major discrepancy between the model and the experimental data seems to be due to the time course of activation level. The model seems to activate at a higher rate than the experimental prediction initially (causing a higher velocity of the CE as the SEE is stretched) and then falls at too slow a rate during deactivation. In contrast, at the slowest frequency of 0.71 Hz, the predicted activation level has a slower rate of decline during deactivation,which results in a maintenance of force.

Activation level changes with shortening speed (a phenomenon known as shortening deactivation). To account for this we optimised the activation parameters that control the uptake of the activator (K and n) to determine whether this effect could be accounted for at each individual cycle frequency (Table 1). The activation parameters (K and n) were optimised to minimise the sum of the force differences between the model and the experimental results at each time point. This was done for three individual conditions at each cycle frequency and the average of these taken as the activation parameters for each cycle frequency. The resultant relationship between the activator concentration and activation level for the optimised activation parameters for each condition is shown in Fig. 5. A comparison between the model with the optimised activation parameters and the experimental results for force output, activation level, energetic output and muscle fibre velocity is shown in Fig. 6.

A comparison of the energetic output of the muscle and the model (with optimised activation parameters for cycle frequency) suggests that the model makes a reasonable prediction of the experimental energetic output (consisting of work + heat) at all speeds during activation, but does less well during deactivation (the period when the muscle is still producing force while no stimulation is being applied) (Fig. 6). This is true whether the experimental work output along with the predicted heat output (using the force, CE length and the activation level) is used to calculate the energy, or whether the model's work output is used. The biggest discrepancy between the model and the experimental results occurs during deactivation at 0.71 Hz. It is apparent from the experimental results that there is a continuation of energy output (in the form of heat)shortly after the cessation of force production. In comparison, the model ceases to produce heat when the force reaches zero. Therefore the final energy expenditure predicted by the model is smaller than that of the experimental findings at this speed. There are also discrepancies between the time course of the experimental energetic output and the model during the 1.25 Hz trials(Fig. 4), where there seems to be some delay between the traces. However the rate of energetic output and the final energetic output compares favourably.

Power output (work/cycle time) and efficiency was also estimated by the model across a range of duty cycles and compared to the average data reported by Curtin and Woledge (1996)(Fig. 7). These simulations were performed at the optimal stimulus phase for each duty cycle as reported by Curtin and Woledge (1996). The optimised activation constants (K and n) for each frequency were used to generate these data(Table 1). The model reproduced the experimental relationship between power and duty cycle and also efficiency and duty cycle, and can be used to predict the duty cycles where optimum power and efficiency occur for all cycle frequencies. The magnitude of power output and efficiency calculated by the model were also accurate for these conditions.

The force and length data reported by Barclay(1994) were scaled according to the maximum isometric force and the optimal muscle fibre length. The parameters of the model were kept essentially the same, however; the stimulation rate constants used were τ₁=0.045 andτ ₂=0.045 and V_max=4.0 L_o s^–1. A comparison of the time course of changes in muscle length, force production and energy output is shown in Fig. 8A. It is apparent from this figure that there is a reasonable prediction of force and work output across the two cycles although, as for the previous comparisons, the model is less reliable during deactivation. Further discrepancies occurred in the time course of heat output, where the model shows a higher rate of heat output during the force production, while the experimental data suggests that this heat output continues to rise as the muscle deactivates and force declines. Despite this discrepancy, the absolute value of energy (heat + work) output across whole cycles seems to be very similar.

Fig. 8B shows the power output and heat rate output (heat/cycle time) of the mouse soleus across a range of frequencies for both experimental conditions and model predictions. The results indicate that the relationship between cycle frequency and energy output is predicted well by the model. Again, the heat rate output from the model has a slightly lower magnitude than the experimental data; however, it is always within 25% of the measured values and it is apparent that the model also accurately predicts the optimal frequency of length change for maximising total energy output.

A comparison of the time course of force and energy changes between the model and the experimental condition yielded positive results that are valuable in our understanding of muscle mechanics. It was possible to use a modified Hill-type muscle model, with an energetic component, to reproduce the optimum activation conditions for achieving maximum power output and maximum efficiency of muscle, as per the results of Curtin and Woledge(1996).

The model makes robust predictions for the time course of the force, energy(heat + work), activation level and contractile element velocity(Fig. 2) when the activation parameters are optimised for force alone. Although the predicted activation level is based partly on the model itself, it does demonstrate small peaks in the activation level, which correspond to each individual muscle twitch (with a time delay of approximately 0.05 s). Providing enough is known about the properties of the muscle in question (force–velocity, force–length and series elastic properties), this technique could be used to estimate the activation level of a muscle in a range of activities where force and length of the muscle–tendon unit is directly measured.

Despite the model's ability accurately to predict the time course of force production at the 1.25 Hz frequency, the results from Fig. 3 demonstrate that the model is less accurate at 0.71 Hz and 5 Hz; this is most apparent during deactivation. At the fastest frequency, it is apparent that the model maintains a high force level once the real muscle–tendon complex begins to lengthen. The experimental results show that the muscle force is low during this period. Analysis of the predicted contractile element velocity from the experimental results suggests that the contractile element needs to be lengthening during the deactivation, rather than shortening, as predicted by the model. To resolve this problem, the activation parameters need to be changed so that the muscle can deactivate at a faster rate. The opposite effect is required at low frequencies, with a reduction in the deactivation rate required.

Numerous investigators have described a phenomenon termed `shortening deactivation', whereby at high velocities of muscle shortening, the muscle tends to deactivate and the force trace is depressed(Askew and Marsh, 2001; Josephson, 1999; Leach et al., 1999). The mechanism behind shortening deactivation is not well known. However the results of this study both support its existence and also provide some information as to how the cycle frequency may influence the activation level. Optimising the activation constants (τ₁, τ₂, K and n) to minimise the sum of the force differences between the model and the experimental results for each of the nine individual conditions (Fig. 4) revealed that the constants τ₁ and τ₂ could remain relatively constant and still provide the best fit. By varying just the parameters K and n, it was possible to get good fits between the model and the experimental force–time data for each individual frequency.

The constants K and n can be thought of as representing the rate of binding of the activator (Ca²⁺) to the troponin, which allows for binding and dissociating of the crossbridges and hence force production. Recent experimental evidence suggests that the off-rate of calcium from troponin increases with the dissociation of the force-generating crossbridges (which occurs with increasing speed of contraction; Wang and Kerrick, 2002). Therefore the mechanism behind the phenomenon of shortening deactivation may be the change in affinity of Ca²⁺ to troponin. The predicted change in the relationship between the activator and the activation level demonstrated here (Fig. 5)provides further evidence that shortening deactivation results from a change in the affinity of the Ca²⁺ to troponin. However, although the optimisation procedure showed that the optimal values of K and n could be characterised across a range of contraction conditions at any given cycle frequency, optimisation under different contraction conditions within the same cycle frequency did show some variation in the activation constants. Therefore the instantaneous speed of contraction is likely to be important, not just the cycle frequency. Further investigation into this area is beyond the scope of this paper and would require a vigorous experimental protocol on live muscle bundles.

The energetic model has been shown to perform relatively well at all frequencies, which is reflected by the ability of the model to predict the duty cycle that produces optimal efficiency. The rate of energetic output(heat + work) during activation in the model provides consistently good agreement with the model. Discrepancies in the onset of the energetic output at 1.25 Hz may be due to the experimental setup, where the muscle may have shifted across the thermopile during contraction. However, it is apparent that the same total energy is measured during the period of one cycle. This is not the case in the 0.71 Hz contractions, where although the energetic outputs of the experiment and the model match during activation, they do not agree during deactivation and as a result the total energetic output during the cycle is underestimated by the model.

The discrepancies in both force and energetic output during deactivation highlight some possible processes that need further investigation within contraction dynamics of muscle. A common finding in the experimental data is that during the longer periods of activation (>0.2 s), the decline in force is associated with a delayed rise in the rate of energetic cost. This is not simulated in the model, which instead predicts a fall in rate of energetic cost once force has declined. This delayed onset of heat production has been cited elsewhere and can partly be explained by the release of heat due to conversion of work by the CE and partly by ATP turnover due to crossbridge cycling (Linari et al., 2003; Curtin and Woledge, 1996). Another possible source for some of the energy liberation during the fall of the force is hysteresis of the elastic tissues. During shortening of the elastic tissues, some of the energy stored in them is lost as heat(Wilson and Goodship, 1994). In biological tissues the range of energy liberated as heat could be as much as 7–30% of total energy stored(Maganaris and Paul, 2000; Pollock and Shadwick,1994).

The experimental muscle continues to produce force for a significant time after the cessation of stimulation at 0.71 Hz compared to the model, even with the optimised activation constants (Fig. 6B). This result suggests that crossbridges are still attached,either due to continuation of ATP turnover, or perhaps some other passive process. The experimental observation that the rate of energetic cost actually plateaus during this period of force maintenance (before increasing again during unloading; see Fig. 6B,duty factor=0.6) suggests that ATP turnover is not responsible for this force maintenance, and instead some other process is involved. The plateau in the force record may also be due to an experimental artefact; however, inspection of experimental trials with similar activation conditions suggests that this phenomenon is consistent for a range of conditions. Therefore perhaps some parallel structure at the fibre level (possibly elastic) is being engaged to produce this force as the force maintenance occurs during muscle lengthening.

The model was highly successful at predicting the various conditions under which the optimal power output and efficiency could occur across two different muscle types. The comparisons to a second set of muscle data, the mouse soleus muscle of Barclay (1994),yielded very positive results for the extension of the model to other muscle types. As with the dogfish muscle, the model was particularly successful at mapping the optima for power output and the rate of energy output, despite changing only three parameters from those used in the dogfish model (the activation constants τ₁ and τ₂ and V_max). The good results may also be assisted by the faster relaxation rate of the mouse muscle, which may hide some of the strange phenomena that occur in the dogfish muscle during relaxation.

Accurate modelling of muscle can effectively allow investigators to simulate large amounts of muscle experiments where the conditions of muscle activation and length changes are changed. Experimentation with muscle fibres,bundles or whole muscles is limited by the life of the muscle. Hence, changing the conditions under which contractions are performed, such as duty cycle,phase of activation and frequency, is difficult without fatiguing/damaging the muscle. Instead, a thorough modelling approach such as that presented here is very useful for determining why muscles function the way they do. More accurate muscle models can also improve simulation of movement with forward dynamics and allow us to determine the effect that varying muscle properties has on muscle mechanics and energetics. Caution should, however, be used when applying this model of energetics across a broad range of muscle types. Knowledge of the properties of individual muscle types (both of the CE and the SEE) is essential in applying this model. These properties are known to vary greatly across the biological spectrum and care should be taken in determining these properties before applying the model.

Although the model predicts the optimal power output and efficiency conditions, further refinement to the model may improve its robustness under varying conditions. For instance, the current model neglects the force–length relationship of muscle because the amplitude of length change is not thought to be large enough to exceed the plateau of this relationship. During animal movement, however, muscles are often subject to length changes that exceed the plateau and some muscles routinely operate in the ascending limb of the force–length relationship. Therefore,application of the energetic model to biological cases should include a scaling of the energy consumed by this relationship.

In conclusion, it has been demonstrated that a Hill-type muscle model can effectively predict the energetics of muscle contraction (heat + work) for two different muscle types using experimentally determined muscle properties. Using the model, it was demonstrated that the activation parameters for achieving optimal power output and optimal efficiency can be predicted and are in line with experimental data for most conditions. With increases in cycle frequency, it was necessary to vary the activation parameters that control the affinity of the activator (Ca²⁺) to the force generator (troponin)in such a way that the off-rate of the activator was increased. This provides further evidence for the phenomenon known as shortening deactivation. The validated model is useful for exploring how activation conditions affect power output and efficiency of a muscle, and how properties of the muscle affect these relationships.

The authors would like to sincerely thank Chris Barclay (Griffith University, Australia), Nancy Curtin (Imperial College London, UK) and Roger Woledge (Kings College London, UK) for graciously providing data and contributing valuable information and advice for the preparation of this paper. The authors would also like to thank the Royal National Orthopaedic Hospital Trust, the BBSRC and the British Council for supporting this work.