The power output of a muscle and its efficiency vary widely under different activation conditions. This is partially due to the complex interaction between the contractile component of a muscle and the serial elasticity. We investigated the relationship between power output and efficiency of muscle by developing a model to predict the power output and efficiency of muscles under varying activation conditions during cyclical length changes. A comparison to experimental data from two different muscle types suggests that the model can effectively predict the time course of force and mechanical energetic output of muscle for a wide range of contraction conditions, particularly during activation of the muscle. With fixed activation properties, discrepancies in the work output between the model and the experimental results were greatest at the faster and slower cycle frequencies than that for which the model was optimised. Further optimisation of the activation properties across each individual cycle frequency examined demonstrated that a change in the relationship between the concentration of the activator (Ca2+) and the activation level could account for these discrepancies. The variation in activation properties with speed provides evidence for the phenomenon of shortening deactivation, whereby at higher speeds of contraction the muscle deactivates at a faster rate. The results of this study demonstrate that predictions about the mechanics and energetics of muscle are possible when sufficient information is known about the specific muscle.

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.

### Force–time predictions

The successful prediction of force production of the Scyliorhinus canicula Linnaeus 1758 white myotomal muscle during sinusoidal movements,based on a two-element Hill-type model, was outlined by Curtin et al.(1998). That model used experimentally determined values for series elastic stiffness and force–velocity properties of the muscle and determined suitable activation parameters to produce an accurate representation of the time course of muscle force production during sinusoidal and ramp length changes. The same relationships were utilised in this study, with the contractile force–velocity relationship and series elastic force–stiffness relationships (and the resultant force–length relationship) illustrated in Fig. 1. A detailed description of how these properties are determined is given in Curtin et al.(1998). The normalised series elastic stiffness is defined by the following equation:
$\ S=\left(\frac{S_{\mathrm{H}}+S_{\mathrm{L}}}{2}\right)+\left[\left(\frac{S_{\mathrm{H}}+S_{\mathrm{L}}}{2}\right){\times}(1-\mathrm{e}^{(-50x|P-X_{\mathrm{o}}|)}){\times}\left(\frac{(P-X_{\mathrm{o}})}{|P-X_{\mathrm{o}}|}\right)\right],$
(1)
where S is the relative SEE stiffness (relative to maximum isometric force divided by optimal muscle fibre length, Po/Lo), SH is the upper limit to the relative stiffness, SL is the lower limit to the relative stiffness and Xo is the value of force relative to the maximum isometric force (Po) where the stiffness switches between the SL and SH in an exponential fashion, and P is the instantaneous force produced by the muscle. The properties of the bundle of Scyliorhinus canicula white myotomal muscle fibres are listed in Table 1.
Table 1.

Properties of the dogfish Scyliorhinus canicula white myotomal muscle and mouse soleus muscle

Cycle frequencyMuscle properties
Optimised activation parameters
Muscle typePo (N)Lo (mm)Vmax (Lo s-1)GSL (Po/Lo)SH (Po/Lo)Xoτ1τ2Kn
Dogfish 0.71 47 7.3 3.8 16 22 0.15 0.035 0.1 0.16 2.2
Dogfish 1.25 60 7.3 3.8 16 22 0.15 0.035 0.1 0.19 2.8
Dogfish 50 7.3 3.8 16 22 0.15 0.035 0.1 0.25
Mouse 48 11.5 16 22 0.15 0.045 0.045 0.22 2.69
Cycle frequencyMuscle properties
Optimised activation parameters
Muscle typePo (N)Lo (mm)Vmax (Lo s-1)GSL (Po/Lo)SH (Po/Lo)Xoτ1τ2Kn
Dogfish 0.71 47 7.3 3.8 16 22 0.15 0.035 0.1 0.16 2.2
Dogfish 1.25 60 7.3 3.8 16 22 0.15 0.035 0.1 0.19 2.8
Dogfish 50 7.3 3.8 16 22 0.15 0.035 0.1 0.25
Mouse 48 11.5 16 22 0.15 0.045 0.045 0.22 2.69

The parameters are explained in Eq. 1-3 and 7 and in the List of symbols and abbreviations. Activation parameters (τ12, K and n) were optimised to produce a best fit for the force—time data across a range of activation conditions at each cycle frequency. Po represents the maximum isometric force of the muscle and differs between cycle frequencies due to muscle fatigue in the experimental protocol (for details of this procedure, see Curtin and Woledge, 1996). Lo is the optimal length to achieve Po.

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 (Ca2+) 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).

We have found that this model is not very accurate at predicting the time course of force rise and decline during twitch contractions or during deactivation of the muscle (after cessation of stimuli whilst force is still produced). Therefore a new model of stimulation was developed that was designed to respond to each individual stimulus rather than trains of stimuli. This stimulation model was essentially designed to model the influx and efflux of the activator (Ca2+ concentration) associated with an individual stimulus (twitch). This can be explained with the following equation:
$\ \frac{\mathrm{d}a}{\mathrm{d}t}=\frac{(1-a)}{{\tau}_{1}}{\\ }(\mathrm{during\ individual\ stimulus})=\frac{-a}{{\tau}_{2}}{\\ }(\mathrm{otherwise}),$
(2)
where a is the activator (Ca2+) concentration and theτ 1 and τ2, respectively, are the time constants dictating the time course of rise and fall of this activator. The individual stimulus can be thought of as a gate opening and allowing an influx of calcium. This gate opens only for a finite period of time, depending on the amplitude and time of the individual stimulus. This model is supported experimentally from measurements of free Ca2+ transients(Baylor and Hollingworth,1998), whereby the activator concentration rises faster than it falls. This relationship can be seen in Fig. 2B where the distinct peaks of activator concentration can be seen in the activator concentration during the twitch.
The crossbridge activation level was modelled in the same way as in Curtin et al. (1998), where the activation level depends on the free concentration of the activator(a) according to the following equation:
$\ Act=\frac{a^{\mathrm{n}}}{(a^{\mathrm{n}}+K^{\mathrm{n}})},$
(3)
where Act is the crossbridge activation level, n is the Hill coefficient for expressing the cooperativity of binding and K is the value of a at which 50% of the crossbridge activation sites are occupied. The time constants for rise and fall of the activator and the binding coefficients used to calculate the relationship of the activator to the crossbridge activation level were optimised to fit the response to a single stimulus trial from the raw data used in Curtin and Woledge(1996)(Fig. 3).
Fig. 1.

Properties of the muscle and definition of terms used in the model. These Scyliorhinus canicula white myotomal muscle properties were determined experimentally in the work of Curtin et al.(1998). (A)Force–velocity relationship of the contractile component at different levels of activation. (B) Variation in the SEE stiffness with force and (C)the resulting force–length relationship of the SEE. (D) Phase of activation is defined as the time between the start of stimulation and the start of shortening expressed as a percentage of cycle duration and is demonstrated with respect to one cycle of length change (a negative value corresponds to activating the muscle whilst the muscle is stretching). Duty cycle is expressed as the fraction of the cycle that the muscle/model is stimulated.

Fig. 1.

Properties of the muscle and definition of terms used in the model. These Scyliorhinus canicula white myotomal muscle properties were determined experimentally in the work of Curtin et al.(1998). (A)Force–velocity relationship of the contractile component at different levels of activation. (B) Variation in the SEE stiffness with force and (C)the resulting force–length relationship of the SEE. (D) Phase of activation is defined as the time between the start of stimulation and the start of shortening expressed as a percentage of cycle duration and is demonstrated with respect to one cycle of length change (a negative value corresponds to activating the muscle whilst the muscle is stretching). Duty cycle is expressed as the fraction of the cycle that the muscle/model is stimulated.

Fig. 2.

The rise and fall of the calcium transients (green) and the resulting activation level (blue) are shown for the original block model (A) and the twitch model (B). Twitches were applied to the model at a frequency of 22 Hz,the same as in the experimental design, and each individual stimulus (twitch)lasted 7.5 ms. The constants τ1 (0.03), τ2(0.10), K (0.19) and n (2.8) were optimised to produce fits to respond to a single stimulus condition(Fig. 3). A time constant,(d=0.015 s), was also required to delay the onset of activation, as was seen in the experimental results. This time constant is thought to be due to the time taken for ATP turnover to proceed and for the calcium signalling to cause a calcium influx.

Fig. 2.

The rise and fall of the calcium transients (green) and the resulting activation level (blue) are shown for the original block model (A) and the twitch model (B). Twitches were applied to the model at a frequency of 22 Hz,the same as in the experimental design, and each individual stimulus (twitch)lasted 7.5 ms. The constants τ1 (0.03), τ2(0.10), K (0.19) and n (2.8) were optimised to produce fits to respond to a single stimulus condition(Fig. 3). A time constant,(d=0.015 s), was also required to delay the onset of activation, as was seen in the experimental results. This time constant is thought to be due to the time taken for ATP turnover to proceed and for the calcium signalling to cause a calcium influx.

If the force–velocity properties of the CE and the force–length properties of the SEE are known, it is also possible to determine the activation level of a muscle fibre bundle from its force and length changes in time. The activation level basically scales the force–velocity curve(Fig. 1) and therefore,providing one knows the force (and hence the stretch of the series elastic component) and also the velocity of the contractile component, it is possible to estimate the activation level; i.e. the percentage of the total maximum number of crossbridges bound. This is shown numerically below:
$\ Act=P{\\ }{/}P^{{^\prime}},$
(4)
where P is the instantaneous force produced by the muscle and P′ is maximum isometric force scaled for the instantaneous muscle velocity. Therefore an estimated activation level for each experimental condition could be calculated from raw experimental data.

### Energetic model

Efficiency is defined as the work produced by a mechanical system divided by the energetic cost of doing that work; this represents the mechanical efficiency (Ettema, 2001). Efficiency is therefore defined as:
$\ \mathrm{Efficiency}=\mathrm{Work}{/}(\mathrm{Heat}+\mathrm{Work}).$
(5)
Here we define work as the integral of the force over the change in length of the whole muscle–tendon complex, and heat as the heat produced by the muscle. Positive work is defined as mechanical work produced in shortening the muscle complex, while stretching the muscle complex with an external force is seen as negative work, or work absorbed by the muscle as a whole. If the elongation occurs in the SEE, work is elastic deformation and the energy can be recovered with 100% return (no hysteretic damping in this model). If the elongation occurs in the CE, this work is converted to heat; this has a small metabolic cost in the model.
Fig. 3.

Comparison of output of model with that of the experimental results for the single optimised condition (frequency=1.25, duty cycle=0.121, stimulus phase=5). The experimental results are taken from raw data used in Curtin and Woledge (1996). (A) Length trajectory (relative to Lo) of MTU for the model (dotted line) and experimental results (solid line); (B) Force (relative to Po) output for the model (dotted line) and experimental results (solid blue line); (C) Energy output (heat + work) for the model (blue dotted line) and experimental results (blue solid line). The experimental results show small amplitude fluctuations as a result of heat measurement from thermopile. The experimental force recordings (solid blue line in B) and the estimated activation level (red line in D) were used as inputs into the energetic model to approximate energetic output using the experimental results(red line in C). (D) Activation level (the relative number of attached crossbridges) predicted by the model (dotted blue line), estimated from the experimental results (red line) and the stimulation pattern from the experiment (solid blue line). (E) CE velocity (Los–1) as predicted by the model (dotted blue line), and as estimated from the experimental results (solid blue line). This is shown in reference to the MTU velocity (red line).

Fig. 3.

Comparison of output of model with that of the experimental results for the single optimised condition (frequency=1.25, duty cycle=0.121, stimulus phase=5). The experimental results are taken from raw data used in Curtin and Woledge (1996). (A) Length trajectory (relative to Lo) of MTU for the model (dotted line) and experimental results (solid line); (B) Force (relative to Po) output for the model (dotted line) and experimental results (solid blue line); (C) Energy output (heat + work) for the model (blue dotted line) and experimental results (blue solid line). The experimental results show small amplitude fluctuations as a result of heat measurement from thermopile. The experimental force recordings (solid blue line in B) and the estimated activation level (red line in D) were used as inputs into the energetic model to approximate energetic output using the experimental results(red line in C). (D) Activation level (the relative number of attached crossbridges) predicted by the model (dotted blue line), estimated from the experimental results (red line) and the stimulation pattern from the experiment (solid blue line). (E) CE velocity (Los–1) as predicted by the model (dotted blue line), and as estimated from the experimental results (solid blue line). This is shown in reference to the MTU velocity (red line).

The rate of heat production from a muscle is a function of crossbridge activation level (Act), velocity of the contractile component(VCE), the time relative to the start of the train of stimulation (t) and the relative force produced (P). For the purpose of this study, where the length of the contractile element remains within the plateau of the force–length relation, length need not be taken into account. The rate of heat production can be further divided into four distinct functions (f) of heat production, which sum to give the overall heat rate:
$\ \frac{\mathrm{d}H}{\mathrm{d}t}=\frac{\mathrm{d}H_{\mathrm{M}}}{\mathrm{d}t}+\frac{\mathrm{d}H_{\mathrm{L}}}{\mathrm{d}t}+\frac{\mathrm{d}H_{\mathrm{S}}}{\mathrm{d}t}+\frac{\mathrm{d}H_{\mathrm{T}}}{\mathrm{d}t}=f(Act)+f(Act,t)+f(Act,V_{\mathrm{CE}})+f(P),$
(6)
where HM has been termed the stable' heat, HL the labile' heat, HS theshortening' heat and HT is the thermoelastic' heat(Aubert, 1956; Hill, 1938; Woledge, 1961).
The stable heat rate can be thought of as the minimum heat rate required to produce an isometric force at any given activation state. This includes the heat produced to activate the muscle (transportation of Ca2+ to activate muscle) and heat produced to maintain force production at the level of the crossbridge. Numerous investigators working on a variety of skeletal muscles have found that this stable heat rate can be approximated by a constant in the range of (a×b), the product of Hill's force–velocity constants (Woledge et al., 1985). When normalised for PoLo units and scaled for activation level:
$\ \frac{\mathrm{d}H_{\mathrm{M}}}{\mathrm{d}t}=Act\left(\frac{a}{P_{\mathrm{o}}}{\times}\frac{b}{V_{\mathrm{max}}}{\times}\frac{V_{\mathrm{max}}}{L_{\mathrm{o}}}\right)=Act\left(\frac{V_{\mathrm{max}}}{G^{2}}\right),$
(7)
where Vmax is the maximum shortening velocity and G=Po/a=Vmax/b.
Over the time course of a contraction the heat rate is not completely stable. Aubert (1956) described a phenomenon he termed labile heat production, where if a muscle is contracted over a period of time the maintenance heat rate could fall from a rate of 2–3 times that of the stable heat rate in an exponentially decaying manner. He termed this extra heat the labile heat'. Assuming that the stable heat rate is as in Eq. 5 and using constants to control the rate of decay of the labile heat rate (adapted from data of Linari et al., 2003) we get the equation:
$\ \frac{\mathrm{d}H_{\mathrm{L}}}{\mathrm{d}t}=\left\{\left[0.8\left(\frac{\mathrm{d}H_{\mathrm{M}}}{\mathrm{d}t}\right)e^{(-0.72^{{\ast}}t)}\right]+\left[0.175\left(\frac{\mathrm{d}H_{\mathrm{M}}}{\mathrm{d}t}\right)e^{(-0.022^{{\ast}}t)}\right]\right\}.$
(8)
Shortening' heat rate can be thought of as the extra energetic cost associated with shortening muscle at any given activation level. Once again,numerous investigators have found a relationship between velocity of the contractile component and the heat rate and it has been approximated to a linear relationship with respect to velocity with a gradient of a(Woledge et al., 1985). Normalising for PoLo:
$\ \frac{\mathrm{d}H_{\mathrm{S}}}{\mathrm{d}t}=Act\left(\frac{a}{P_{\mathrm{o}}}{\times}\frac{V_{\mathrm{CE}}}{L_{\mathrm{o}}}\right)=Act\left(\frac{V_{\mathrm{CE}}}{G}\right)(\mathrm{when}{\\ }V_{\mathrm{CE}}{>}0,\mathrm{shortening}),$
(9)
where VCE is the contractile element velocity relative to Lo (i.e. Lo s–1). Like the maintenance heat rate, the effective shortening heat rate must also be scaled for activation as it is dependent on the number of bound crossbridges.
Energy output is reduced as a result of active lengthening. Studies by Lou et al. (1998) revealed that during an isometric contraction, at least 30% of the heat produced by muscle was the result of activating the muscle (i.e. calcium turnover). Therefore,during active lengthening, the minimum heat rate must be at least 30% of the stable heat rate. Studies by Linari et al.(2003) also revealed that there is an exponential decay of the rate of energy production as the lengthening velocity increases. This heat rate must also be scaled for activation. During stretch, work done on the contractile component also becomes heat within a short period of time(Linari et al., 2003). This model accounts for this energy. However, it ignores the small time delay. Therefore the heat rate during lengthening can be approximated with the following equation:
$\ \frac{\mathrm{d}H}{\mathrm{d}t}=0.3\left(\frac{\mathrm{d}H_{\mathrm{M}}}{\mathrm{d}t}\right)+0.7\left\{\left(\frac{\mathrm{d}H_{\mathrm{M}}}{\mathrm{d}t}\right)e^{-8\left[\left(\frac{P}{Act}\right)-1\right]}\right\}+P{\times}V_{\mathrm{CE}}(\mathrm{when}{\\ }V_{\mathrm{CE}}{<}0,\mathrm{lengthening}).$
(10)
One further factor that is not associated with the contractile component also contributes to the rate of heat production. This is the thermoelastic effect, which causes heat to be absorbed by muscle. This is proportional to the rate of force production and has been characterised(Woledge, 1961) with the following relationship:
$\ \frac{\mathrm{d}H_{\mathrm{T}}}{\mathrm{d}t}=-0.014\left(\frac{\mathrm{d}P}{\mathrm{d}t}\right).$
(11)
The above model of heat expenditure during dynamic contraction can therefore be applied providing the following parameters are known: force,length of contractile component, activation level, Vmaxand G. Using experimental measurement of the force output and assuming that the stiffness of the SEE is known, it is possible to calculate the length (and velocity) of the SEE and the CE as follows:
$\ L_{\mathrm{CE}}=L_{\mathrm{MTU}}-L_{\mathrm{SEE}}=L_{\mathrm{MTU}}-\left(\frac{P}{S}\right),$
(12)
where LCE is the length of the CE, LMTU is the length of the MTU, LSEE is the length of the SEE, P is the instantaneous force produced by the muscle and S is the relative stiffness of the SEE at that force. Therefore, by extracting these data from the experimental results and by predicting the activation level as in Eq. 4, it is possible to apply the model to experimental results to estimate energetic output.

### Comparison to experimental data and analysis

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 (τ1, τ2, 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 Vmax (scaled to Lo s–1) and the activation constants(τ1, τ2, K and n) and the force data was scaled to Fmax(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.

### Experimental comparisons

The activation parameters (τ1, τ2, 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.

Fig. 4.

Comparison of the model output to the experimental results across a range of stimulation conditions with activation constants optimised to fit a single stimulus condition (Fig. 3). Length, force, energetic output, activation and contractile component velocity are compared between the model output (dotted blue lines) and the experimental results (solid blue lines). The convention for these figures is the same as for the single condition in Fig. 3, with the red lines indicating the modelled energetic output for the experimental data (energy plot), the estimated experimental activation level (activation plot) and the MTU velocity [contractile component (CE)velocity plot]. The comparison is made at three different cycle frequencies:(A) 5 Hz, (B) 1.25 Hz and (C) 0.71 Hz, with varying duty cycle and phase.

Fig. 4.

Comparison of the model output to the experimental results across a range of stimulation conditions with activation constants optimised to fit a single stimulus condition (Fig. 3). Length, force, energetic output, activation and contractile component velocity are compared between the model output (dotted blue lines) and the experimental results (solid blue lines). The convention for these figures is the same as for the single condition in Fig. 3, with the red lines indicating the modelled energetic output for the experimental data (energy plot), the estimated experimental activation level (activation plot) and the MTU velocity [contractile component (CE)velocity plot]. The comparison is made at three different cycle frequencies:(A) 5 Hz, (B) 1.25 Hz and (C) 0.71 Hz, with varying duty cycle and phase.

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.

Fig. 5.

The relationship between the activator (Ca2+) concentration and the activation level (the relative number of attached crossbridges) for each of the optimised values of K and n that produced best fits to the force data at each of the three cycle frequencies (5, 1.25 and 0.71 Hz).

Fig. 5.

The relationship between the activator (Ca2+) concentration and the activation level (the relative number of attached crossbridges) for each of the optimised values of K and n that produced best fits to the force data at each of the three cycle frequencies (5, 1.25 and 0.71 Hz).

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 τ1=0.045 andτ 2=0.045 and Vmax=4.0 Lo 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).

Fig. 6.

Comparison of the model output to the experimental results across a range of stimulation conditions with activation constants optimised to best fit each individual cycle frequency (Fig. 3). Length, force, energetic output, activation and contractile component velocity are compared between the model output (dotted blue lines)and the experimental results (solid blue lines). The convention for these figures is the same as for the single condition in Fig. 3, with the red lines indicating the modelled energetic output for the experimental data (energy plot), the estimated experimental activation level (activation plot) and the MTU velocity [contractile component (CE) velocity plot]. The comparison is made at the frequencies of (A) 5 Hz and (B) 0.71 Hz.

Fig. 6.

Comparison of the model output to the experimental results across a range of stimulation conditions with activation constants optimised to best fit each individual cycle frequency (Fig. 3). Length, force, energetic output, activation and contractile component velocity are compared between the model output (dotted blue lines)and the experimental results (solid blue lines). The convention for these figures is the same as for the single condition in Fig. 3, with the red lines indicating the modelled energetic output for the experimental data (energy plot), the estimated experimental activation level (activation plot) and the MTU velocity [contractile component (CE) velocity plot]. The comparison is made at the frequencies of (A) 5 Hz and (B) 0.71 Hz.

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.

Fig. 7.

Comparison of the experimental (blue) power output (A) and efficiency (B)to the model predictions (green) for a range of cycle frequencies and increasing duty cycle. The same phase of activation was used in each condition. Model results use the optimal activation parameters for each cycle frequency, as shown in Table 1.Power is defined as the work per cycle time and is scaled up from PoLo/cycle time units based on the properties of the muscle reported by Curtin et al. (1996).

Fig. 7.

Comparison of the experimental (blue) power output (A) and efficiency (B)to the model predictions (green) for a range of cycle frequencies and increasing duty cycle. The same phase of activation was used in each condition. Model results use the optimal activation parameters for each cycle frequency, as shown in Table 1.Power is defined as the work per cycle time and is scaled up from PoLo/cycle time units based on the properties of the muscle reported by Curtin et al. (1996).

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 (τ1, τ2, 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 τ1 and τ2 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 (Ca2+) 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 Ca2+ 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 Ca2+ 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).

Fig. 8.

(A) Comparison of the model (dotted lines) to individual records of length,force, energetic output and stimulation/activation for an example soleus muscle of a mouse (solid lines) as reported in Barclay(1994). Values are scaled according to average data reported by Barclay and personal communication. (B)Comparison of model (dotted lines) to experimental (solid lines) power and rate of heat and total energy output (power + heat rate) across a range of frequencies. The stimulation rate constants used were τ1=0.045 and τ2=0.045 and Vmax=4.0 Lo s–1.

Fig. 8.

(A) Comparison of the model (dotted lines) to individual records of length,force, energetic output and stimulation/activation for an example soleus muscle of a mouse (solid lines) as reported in Barclay(1994). Values are scaled according to average data reported by Barclay and personal communication. (B)Comparison of model (dotted lines) to experimental (solid lines) power and rate of heat and total energy output (power + heat rate) across a range of frequencies. The stimulation rate constants used were τ1=0.045 and τ2=0.045 and Vmax=4.0 Lo s–1.

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 τ1 and τ2 and Vmax). 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 (Ca2+) 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.

 a activator (Ca2+) concentration a,b constants Act crossbridge activation level CE contractile element f function of... G Po/a = Vmax/b HL labile' heat HM stable' heat HS shortening' heat HT thermoelastic' heat K value of a at which 50% of the crossbridge activation sites are occupied LCE length of the CE LMTU length of the MTU Lo optimal muscle fibre length LSEE length of the SEE MTU muscle–tendon unit n Hill coefficient P instantaneous force produced by muscle P′ maximum isometric force scaled by muscle velocity Po normalised maximum isometric force S relative SEE stiffness SEE series elastic element SH upper limit to the relative stiffness SL lower limit to the relative stiffness t time VCE contractile element velocity Vmax maximum shortening velocity Xo force relative to Po where stiffness changes from SH to SL τ1, τ2 time constants
 a activator (Ca2+) concentration a,b constants Act crossbridge activation level CE contractile element f function of... G Po/a = Vmax/b HL labile' heat HM stable' heat HS shortening' heat HT `thermoelastic' heat K value of a at which 50% of the crossbridge activation sites are occupied LCE length of the CE LMTU length of the MTU Lo optimal muscle fibre length LSEE length of the SEE MTU muscle–tendon unit n Hill coefficient P instantaneous force produced by muscle P′ maximum isometric force scaled by muscle velocity Po normalised maximum isometric force S relative SEE stiffness SEE series elastic element SH upper limit to the relative stiffness SL lower limit to the relative stiffness t time VCE contractile element velocity Vmax maximum shortening velocity Xo force relative to Po where stiffness changes from SH to SL τ1, τ2 time constants

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.

Alexander, R. M. (
2002
). Tendon elasticity and muscle function.
Comp. Biochem. Physiol. A
133
,
1001
-1011.
Anderson, F. C. and Pandy, M. G. (
2001
). Dynamic optimization of human walking.
J. Biomech. Eng.
123
,
381
-390.
Askew, G. N. and Marsh, R. L. (
2001
). The mechanical power output of the pectoralis muscle of blue-breasted quail(Coturnix chinensis): the in vivo length cycle and its implications for muscle performance.
J. Exp. Biol.
204
,
3587
-3600.
Aubert, X. (
1956
).
Le Couplage Energetique de la Contraction Musculaire
. Brussels:Arscia.
Barclay, C. J. (
1994
). Efficiency of fast- and slow-twitch muscles of the mouse performing cyclic contractions.
J. Exp. Biol.
193
,
65
-78.
Baylor, S. M. and Hollingworth, S. (
1998
). Model of sarcomeric Ca2+ movements, including ATP Ca2+binding and diffusion, during activation of frog skeletal muscle.
J. Gen. Physiol.
112
,
297
-316.
Biewener, A. A. and Roberts, T. J. (
2000
). Muscle and tendon contributions to force, work, and elastic energy savings: a comparative perspective.
Exerc. Sport Sci. Rev.
28
,
99
-107.
Curtin, N. and Woledge, R. (
1996
). Power at the expense of efficiency in contraction of white muscle fibres from dogfish Scyliorhinus canicula.
J. Exp. Biol.
199
,
593
-601.
Curtin, N. A., Gardner-Medwin, A. R. and Woledge, R. C.(
1998
). Predictions of the time course of force and power output by dogfish white muscle fibres during brief tetani.
J. Exp. Biol.
201
,
103
-114.
Ettema, G. J. (
1996
). Mechanical efficiency and efficiency of storage and release of series elastic energy in skeletal muscle during stretch–shorten cycles.
J. Exp. Biol.
199
,
1983
-1997.
Ettema, G. J. (
2001
). Muscle efficiency: the controversial role of elasticity and mechanical energy conversion in stretch-shortening cycles.
Eur. J. Appl. Physiol.
85
,
457
-465.
Fukunaga, T., Kubo, K., Kawakami, Y., Fukashiro, S., Kanehisa,H. and Maganaris, C. N. (
2001
). In vivo behaviour of human muscle tendon during walking.
Proc. R. Soc. Lond. B
268
,
229
-233.
Galantis, A. and Woledge, R. C. (
2003
). The theoretical limits to the power output of a muscle-tendon complex with inertial and gravitational loads.
Proc. R. Soc. Lond. B
270
,
1493
-1498.
Hill, A. V. (
1938
). The heat of shortening and the dynamic constants of muscle.
Proc. R. Soc. Lond. B
126
,
136
-195.
Josephson, R. K. (
1999
). Dissecting muscle power output.
J. Exp. Biol.
202, 23
,
3369
-3375.
Leach, J. K., Priola, D. V., Grimes, L. A. and Skipper, B. J. (
1999
). Shortening deactivation of cardiac muscle:physiological mechanisms and clinical implications.
J. Invest. Med.
47
,
369
-377.
Linari, M., Woledge, R. C. and Curtin, N. A.(
2003
). Energy storage during stretch of active single fibres from frog skeletal muscle.
J. Physiol.
548
,
461
-474.
Lou, F., Curtin, N. A. and Woledge, R. C.(
1998
). Contraction with shortening during stimulation or during relaxation: how do the energetic costs compare?
J. Mus. Res. Cell Motil.
19
,
797
-802.
Maganaris, C. N. and Paul, J. P. (
2000
). Hysteresis measurements in intact human tendon.
J. Biomech.
33
,
1723
-1727.
Pollock, C. M. and Shadwick, R. E. (
1994
). Relationship between body mass and biomechanical properties of limb tendons in adult mammals.
Am. J. Physiol.
266
,
R1016
-R1021.
Roberts, T. J. (
2002
). The integrated function of muscles and tendons during locomotion.
Comp Biochem. Physiol.
133A
,
1087
-1099.
Umberger, B. R., Gerritsen, K. G. and Martin, P. E.(
2003
). A model of human muscle energy expenditure.
Comput. Methods Biomech. Biomed. Engin.
6
,
99
-111.
Wang, Y. and Kerrick, W. G. (
2002
). The off rate of Ca(2+) from troponin C is regulated by force-generating cross bridges in skeletal muscle.
J. Appl. Physiol.
92
,
2409
-2418.
Wilson, A. M. and Goodship, A. E. (
1994
). Exercise-induced hyperthermia as a possible mechanism for tendon degeneration.
J. Biomech.
27
,
899
-905.
Woledge, R. C. (
1961
). The thermoelastic effect of change of tension in active muscle.
J. Physiol.
155
,
187
-208.
Woledge, R. C. (
1998
). Muscle energetics during unfused tetanic contractions. Modelling the effects of series elasticity.