A new model for force generation by skeletal muscle, incorporating work-dependent deactivation

The major determinants of the amount of force generated by an activated muscle are, in addition to its level of activation, its length and the rate of change of its length (velocity). These properties have been known for many decades from studies measuring the force generated during stimulation with either the length or the velocity of the muscle preparation held constant. More recently, however, it has become clear that these properties are time dependent (for a review, see Josephson, 1999). Some recent attempts have been made to predict the time course of force generated during naturally occurring movements on the basis of assumed time-independent length—tension and force—velocity properties alone (Curtin et al., 1998; Williams et al., 1998; Lichtwark and Wilson, 2005; McMillen et al., 2008). In each of these studies, the predictions worked well during the rising phase of muscle force but tended to fail during relaxation.

This paper describes an expansion of an earlier model of force generation by lamprey muscle (Williams et al., 1998; McMillen et al., 2008), which now incorporates the phenomenon of work-dependent deactivation (WDD), whereby a muscle that has been shortening under load generates less force after the shortening ceases than otherwise expected. This phenomenon was reported in skeletal muscle by Jewell and Wilkie (Jewell and Wilkie, 1960) and later referred to by Edman (Edman, 1975) as shortening-induced deactivation. However, it has been shown in crab muscle (Josephson and Stokes, 1999b) and in whole cat soleus muscle (Herzog et al., 2000) that the observed deactivation is not dependent on the shortening alone but is proportional to the amount of work done by the shortening muscle. Josephson and Stokes thus referred to the phenomenon as WDD (Josephson and Stokes, 1999b). In the present study, this property represents the greatest departure from time independence.

The main goal of this study has been to develop an accurate mathematical model of force generation for incorporation within a comprehensive neural—mechanical—hydrodynamic model of lamprey swimming (Hultmark et al., 2007; McMillen et al., 2008; Tytell and Cohen, 2008; Varkonyi et al., 2008; Hsu et al., 2009).

The data on which this paper is based were collected in collaboration with Professor Nancy Curtin (Imperial College, London, UK), and some of it has been utilised in previous publications (Williams et al., 1998; McMillen et al., 2008). In summary, portions of tissue were dissected from swimming muscle of freshly killed lampreys (Ichthyomyzon unicuspis Hubbs and Trautman 1937). Each preparation was between 2 and 3 mm in length, which constituted the length of one myomere with pieces of the myosepta into which the fibres inserted. The cross-sectional area was approximately 1 mm².

The preparation was mounted horizontally in a tissue bath, with the myoseptal material on one end secured to a force transducer and on the other to a servo-operated motor (Cambridge Technology, Inc., series 300B, Watertown, MA, USA), which controlled the length of the preparation. Electrical stimulation was applied through a pair of bare platinum wires, which ran the length of the tissue bath. The muscle was tetanically stimulated for periods of 0.36 s, representing the activation bursts during swimming at a frequency of one tail-beat cycle per second (Williams et al., 1989).

Three different kinds of experiment were performed. In isometric experiments, the length of the preparation was held constant at 0.96, 1.0 or 1.04 times its in situ length (Fig. 1A); in ramp experiments, the muscle was lengthened or shortened at constant speed during the stimulation (Fig. 1B); during sine experiments, the length of the preparation was sinusoidally varied between 0.96 and 1.04 times the in situ length, at a frequency of 1 cycle per second (Fig. 1C). This range of lengths represents approximately the range estimated to occur during lamprey swimming (McClellan, 1984). The range of velocities in the ramp experiments was chosen to be equal to the range of velocities occurring during the sine experiments. For more details, see Williams et al. (Williams et al., 1998).

A muscle preparation did not always produce reliable data for a full experiment, i.e. a set of isometric, ramp and sine runs. Over time a certain amount of fatigue could occur, more so in some preparations than others. In order to monitor the state of a preparation, isometric tetani were recorded at intervals throughout. If significant change occurred in the peak isometric force, then before fitting the model parameters all the force values in that experiment were scaled by the ratio of the most recent isometric peak value to the first isometric peak value.

A model was devised whose parameters were determined by fitting data collected in the isometric and ramp experiments. The model was then tested by using it to predict the data that had been recorded in the sine experiments, and comparing the model predictions with the recorded experimental data.

The model consists of two linked parts. The first part is a simple model of activation of the muscle protein filaments by Ca²⁺ ions released from the sarcoplasmic reticulum (SR) of muscle, and is similar to that used in previous studies (Williams et al., 1998; McMillen et al., 2008). The second part is a modification of the classic model of force development in muscle based on the mechanical model of Hill (Hill, 1938) and incorporating the length—tension and force—velocity properties of muscle (Jewell and Wilkie, 1958). The model presented in this paper differs from our previous versions in two ways: first, the stiffness of the series elastic component is assumed not to be constant but to depend on the level of activation, and second (and more importantly) the model incorporates the phenomenon of WDD.

The simple kinetic model for activation of the protein filaments assumes that Ca²⁺ is released from the SR at a rate proportional to its bound concentration and is then free to bind to protein filament sites at a rate proportional to the concentration of both free Ca²⁺ and free filament sites. Similar mass-action processes govern the release of Ca²⁺ from the filaments and its rebinding in the SR. However, there were two departures from simple mass-action kinetics. First, k₁, the rate constant for release of calcium from the SR, was assumed to be negligible while the stimulus was off, and likewise k₂, the rate constant for binding of calcium ions to the SR, became negligible while the muscle was being stimulated. Second, it was found necessary to introduce some cooperativity in the release of Ca²⁺ from the contractile filaments, so the rate of release was made proportional to the product of the concentrations of bound and unbound filament sites, with rate constant k₄. Without this modification, relaxation began too abruptly.

The details of this model are given in an earlier publication (McMillen et al., 2008); only the resulting simplified equations are given below: where Ca and Caf are the (non-dimensional) concentrations of free Ca²⁺ and Ca-bound filaments, respectively; during tetanus the value of Caf approaches 1. C and S represent (non-dimensional) total concentrations of Ca²⁺ and sarcoplasmic-reticular binding sites, respectively, and k₁—k₄ are positive rate constants [for derivation of the above equations, see McMillen et al. (McMillen et al., 2008)]. In that study, all four binding parameters were assumed to be constant for a given preparation, but in the version reported here, k₃ and k₄ will be made to vary in response to loaded shortening of the muscle and hence to be responsible for WDD.

In the classic Hill model (Hill, 1938) the SE is purely passive and independent of muscle activation, abiding primarily in its attachments. It has long been known, however, that much of the elasticity attributed to the SE resides in the activated muscle filaments (Huxley and Simmons, 1971; Huxley and Simmons, 1973). Measurements of muscle stiffness (the rate of change of strain with respect to stress) have traditionally been made on fully activated, tetanised muscle [e.g. in frog (Jewell and Wilkie, 1958); in dogfish (Curtin et al., 1998)]. In such studies the value of stiffness has been found to be independent of muscle force except at quite low values of force. However, when the measurements are made during the initial development of tension, stiffness is found to increase as activation rises and muscle force increases [in rabbit cardiac muscle (Edman and Nilsson, 1968); in frog skeletal muscle (Edman and Josephson, 2007)]. In particular, Edman and Josephson (Edman and Josephson, 2007) found that the stiffness of the SE increased linearly with developed force from an intercept value of approximately 4% of the maximal stiffness.

The parameters of these equations (α_m, α_p, P₀, L₀ and λ₂) are determined in the first instance by plotting peak force values in experiments as in Fig. 1 against the appropriate variables. These functions are additionally restricted such that 0≤α(ν_c)≤α_max and 0≤λ(l_c)≤1.

Using Eqn 11 to calculate ν_c from Eqn 6 for use in Eqn 9 can lead to instability in the numerical computation. This can be prevented by deriving another formulation for dP/dt, as follows: Substituting ν_c from Eqn 6 into Eqn 9 and thence into Eqn 10 gives: where and μ is given by Eqn 5.

Substituting Eqn 12 into Eqn 11 and solving for dP/dt gives: where α₁ and μ are as in Eqn 12.

Although the full model (at this point) is given by Eqns 1, 2 and 3, 5, 6 and Eqns 8, 9, 10 and 11, using Eqn 14 instead of Eqn 11 in Eqn 6 aids the numerical computation by providing stability.

It was suggested by Edman (Edman, 1975) that WDD may be caused by a structural change in the myofilament system, and Moréchal and Plaghki (Moréchal and Plaghki, 1979) postulated that inhibition of the formation of cross bridges could be the mechanism. Corr and Herzog, in a study on recovery from WDD in cat soleus muscle, have provided some evidence for this hypothesis (Corr and Herzog, 2005).

In the current study, this property was modelled by introducing a new variable m, which affects the rate of Ca²⁺ binding and release to and from the muscle filaments. This variable increases during muscle shortening under load, at a rate proportional to the product of P_c and —ν_c, which equals the rate at which work is doing done by the activated filaments. When the muscle is not shortening, m decays to its initial value of 1. Thus: where k_m1 and k_m2 are positive rate constants to be determined by fitting the force trajectories of the ramp experiments to the full model.

The effect of this variable m is to change the ratio of k₃ and k₄, the rate constants of the binding of Ca²⁺ to the filaments and its release. This results in smaller values of Caf and, by Eqn 10, smaller force than would otherwise occur. Hence, k₃ and k₄ are no longer constant but are given by: where k₃₀ and k₄₀ are the respective values when m=1. The ratio k₄/k₃ is thus proportional to m, and Eqns 16 must be incorporated into Eqns 1, 2.

The changing value of the ratio k₄/k₃ affects the force obtained during shortening, as well as after shortening ceases. For this reason, the peak values obtained during the ramp experiments are no longer adequate for the determination of the parameter α_m. Instead, all three parameters (α_m, k_m1 and k_m2) are used to fit the full model to the force trajectories of the ramp experiments. Provisional values for k₁, k₂, k₃₀ and k₄₀ are found by fitting the model to the isometric data. Reiterative fitting of the isometric and ramp traces is then used to obtain the final values of all the adjustable parameters.

The behaviour of the variable m during both isometric and sine movements is shown in Fig. 2. During the isometric experiment, ν_c becomes negative during the rising phase of the force (as SE is being stretched by the shortening of CE) and positive during relaxation. During the sine experiment, such changes in ν_c are superimposed upon those caused by the changing length L of the whole preparation.

The full model for force generation in response to stimulation is thus given by Eqns 1, 2 and 3, 5 and 6 and 8, 9, 10, 11, 12, 13, 14, 15 and 16. These equations are then solved, with parameters fitted to the isometric and ramp data, using the least-squares curve-fitting facilities in the software XPPAUT devised by G. Bard Ermentrout and available free of charge at www.pitt.edu/~phase/. The parameters so obtained (listed in Table 1) are then used to predict the force trajectories of the sine experiments.

Data were analysed from four complete experiments, which gave qualitatively similar results. All the data for the figures and in Table 1 have been taken from the one experiment during which the shape and height of isometric tetanus changed least over time.

The classic Hill model of SE in parallel with CE behaves differently when μ is dependent on muscle activation. In the left-hand panel (A) of Fig. 3 are shown the time courses of the variables of Eqns 3 and 4 during an isometric contraction in which μ is kept constant. The effects of non-constant μ (Eqns 3, 5 and 6) are revealed in the right-hand panel (Fig. 3B). In this case, the low value of stiffness at the onset of contraction allows the SE to stretch more quickly and thus both l_c and ν_c plateau before the stimulus ends. Furthermore, because μ falls as the value of the force falls, these plateaus are maintained longer.

As can be seen in the upper panels of Fig. 3, the isometric force data were well fitted with the model in both cases. In general, the model with variable μ (as in Fig. 3B) gave slightly better fits to the force data of both the isometric and ramp experiments than the model with constant μ.

In both panels of Fig. 3 the variable m has been included in the model.

In the isometric case, the variable m makes little difference to the data fit, because its value does not increase much above 1 (see Fig. 2A). However, in both the ramp and sine cases, including the effects of WDD makes a significant difference to the model's behaviour.

Without m, the model could fit only the rising phases of the force traces during shortening ramps, as can be seen in Fig. 4. Panels A and B show the results of the original model, without the variable m, at two different shortening ramp speeds. The isometric traces are included for comparison. Although the rising phase of the force is well fitted for both speeds, the model force falls too slowly during the relaxation phase, especially at the faster shortening speed (Fig. 4B). For Fig. 4C,D, WDD has been included in the model, and it can be seen that the force traces are quite well fitted over the full time course, for both shortening speeds illustrated.

The importance of including WDD is shown further when the model devised from the isometric and ramp data is used to predict the sine data. Fig. 5A shows the predictions when the model does not include the variable m. Although the model captures the main features of the effect of movement phase on force production, the greatest errors in the time courses occur in the smaller phase values, which are the only phases that occur during steady swimming. By contrast, when the variable m is included, the model produces a much closer prediction to the behaviour seen at all phases, as shown in Fig. 5B. The poorest fit here is seen at a phase of 40%, which does not normally occur during swimming (Williams et al., 1989; Altringham and Ellerby, 1999).

All the energy to power swimming must be manifest through work done by the swimming muscles on their attachments. This quantity can be calculated as the time integral of the product of the force P_c generated by the CE and its shortening velocity —ν_c. In Fig. 6A is shown the time course over one cycle of these two quantities, along with their cumulative product, for three phase values. These phase values encompass the range of phases seen during swimming in the lamprey. Although there are periods during which negative work is performed, net positive work is performed over a full cycle at all three phase values.

In Fig. 6B, the net work per cycle at all phases tested is shown for both the model output and the experimental data. It can be seen that the model reproduces the general pattern of work done at all phases.

This model of force development can predict with reasonable accuracy the forces and power observed in response to stimulation of muscle during imposed movements simulating those occurring naturally. It thus represents a contribution to efforts to understand the complex neuro-mechanical interactions occurring during natural movements. In addition, this demonstrates that an adequate description of muscle behaviour under these conditions can be achieved by adding an analysis of WDD to the force—velocity and length—tension properties of muscle.

It should be noted that the particular set of values found for the parameters of the chemical portion of the model (C, S, k₁—k₅) is not a unique set, and the goodness of the fit cannot be taken as evidence that the details of the model provide a quantitative description of reality. Nonetheless, any set of parameters that provided a good fit for the isometric and ramp data also provided a good prediction for the sine data; this is the unique power of the model. Furthermore, it may be possible to adapt the model to other kinds of muscle, by an appropriate choice of parameters.

The data on which the model were based were collected under conditions of full tetanic activation whereas muscle in situ is not generally tetanically activated. This can be accommodated when using the model to simulate normal behaviour by scaling the value of P₀. The values of both P₀ and L₀ need to be altered at any rate, to match a particular modelling task, for example to approximate the changing body geometry along the length of the lamprey.

The observation that the model works somewhat better when the stiffness of the SE depends on activation cannot be taken as evidence that this is the case in reality, because the model of Ca²⁺ activation is oversimplified and ad hoc. Further experiments such as those by Edman and Josephson (Edman and Josephson, 2007) are necessary to investigate the properties of the SE in non-tetanised muscle. Likewise, implementing WDD by changing the ratio of the Ca-binding rate constants of our model cannot be taken as evidence that this reflects the actual mechanism for this effect, although it may be a factor. That such accurate predictions can be made on the basis of such simple assumptions is very heartening to a modeller.

A further property of muscle that has not been included in the current model is a time-dependent augmentation of force by muscle stretch (for reviews, see Herzog, 1998; Josephson and Stokes, 1999a). The inadequacy of the fit for the sine experiment at 40% in Fig. 5B can be attributed to this property. Because the muscles of swimming fish do not normally operate at this phase (see Altringham and Ellerby, 1999), no attempt has been made to include this property in the model. However, further work needs to be done on this phenomenon.

I am grateful to Nancy Curtin for consenting to my use of our published data, and to Philip Holmes and Tyler McMillen for very helpful comments on the manuscript. Travel support was provided by NIH CNRS Grant 1RO1NS054271. Deposited in PMC for release after 12 months.

 
• C

  total Ca concentration

 
• Ca

  free Ca²⁺ concentration

 
• Caf

  Ca-bound filament concentration

 
• CE

  contractile component

 
• F

  total filament concentration

 
• l_c

  length of CE

 
• l_c0

  length of CE at its optimum

 
• L

  length of preparation

 
• L_is

  in situ length

 
• L₀

  optimum length

 
• m

  variable responsible for WDD

 
• P

  force developed by SE

 
• P_c

  force developed by CE

 
• P₀

  peak force in an isometric tetanic contraction

 
• S

  total concentrations of sarcoplasmic-reticular binding sites

 
• SE

  series elastic component

 
• SR

  sarcoplasmic reticulum

 
• ν_c

  time derivative of l_c

 
• V

  time derivative of L

 
• WDD

  work-dependent deactivation

 
• α(ν_c)

  dependence of P_c on ν_c

 
• λ(l_c)

  dependence of P_c on l_c

 
• μ

  stiffness of SE

 
• μ₀

  intercept of the dependence of SE stiffness on Caf

 
• μ₁

  gradient of the dependence of SE stiffness on Caf