## SUMMARY

A model is developed to predict the force generated by active skeletal muscle when subjected to imposed patterns of lengthening and shortening, such as those that occur during normal movements. The model is based on data from isolated lamprey muscle and can predict the forces developed during swimming. The model consists of a set of ordinary differential equations, which are solved numerically. The model's first part is a simplified description of the kinetics of Ca^{2+} release from sarcoplasmic reticulum and binding to muscle protein filaments, in response to neural activation. The second part is based on A. V. Hill's mechanical model of muscle, consisting of elastic and contractile elements in series, the latter obeying known physiological properties. The parameters of the model are determined by fitting the appropriate mathematical solutions to data recorded from isolated lamprey muscle activated under conditions of constant length or rate of change of length. The model is then used to predict the forces developed under conditions of applied sinusoidal length changes, and the results compared with corresponding data. The most significant advance of this model is the incorporation of work-dependent deactivation, whereby a muscle that has been shortening under load generates less force after the shortening ceases than otherwise expected. In addition, the stiffness in this model is not constant but increases with increasing activation. The model yields a closer prediction to data than has been obtained before, and can thus prove an important component of investigations of the neural—mechanical—environmental interactions that occur during natural movements.

## INTRODUCTION

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).

## MATERIALS AND METHODS

### Experimental data

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^{2}.

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.

### Mathematical model

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^{2+} 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.

### Activation by Ca^{2+}

The simple kinetic model for activation of the protein filaments assumes that Ca^{2+} 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^{2+} and free filament sites. Similar mass-action processes govern the release of Ca^{2+} from the filaments and its rebinding in the SR. However, there were two departures from simple mass-action kinetics. First, *k*_{1}, the rate constant for release of calcium from the SR, was assumed to be negligible while the stimulus was off, and likewise *k*_{2}, 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^{2+} 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*_{4}. 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^{2+} and Ca-bound filaments, respectively; during tetanus the value of *Caf* approaches 1. *C* and *S* represent (non-dimensional) total concentrations of Ca^{2+} and sarcoplasmic-reticular binding sites, respectively, and *k*_{1}—*k*_{4} 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*_{3} and *k*_{4} will be made to vary in response to loaded shortening of the muscle and hence to be responsible for WDD.

### Mechanical model

*L*of the muscle is given by the sum of the length

*l*

_{c}of the CE and the length of the SE. When the SE is stretched, it develops a tension force,

*P*, given by the product of its stretched length and its stiffness, μ (its unstretched length assumed not to be significantly different from zero). This gives:

*ν*

_{c}and

*V*are the rates of change of

*l*

_{c}and

*L*, respectively,

*P*the force developed, and μ the stiffness of the SE. Traditionally,

*V*is taken as positive during muscle shortening but in these equations

*V*is defined as the time derivative of length and is thus negative for shortening, positive for lengthening.

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.

_{0}and μ

_{1}are the intercept and gradient of the dependence of SE stiffness on

*Caf*.

#### Physiological properties

*Caf*and the effects of the length and velocity of the CE:

*P*

_{0}is the peak force in an isometric tetanic contraction (when

*Caf*=1) at the optimum length

*L*

_{0}. In all four muscle preparations, the value of

*L*

_{0}was somewhat larger than the length

*in situ, L*

_{is}(see Table 1). The functions λ(

*l*

_{c}) and α(

*ν*

_{c}) are estimated from force measurements (as in Fig. 1). Over the range of lengths and velocities used, the data are fit well enough for the model by a parabolic function for λ(

*l*

_{c}) and a piecewise linear function (approximating a hyperbola) for α(

*ν*

_{c}):

*l*

_{c0}is the length of the CE at the optimum length,

*L*

_{0}, from Eqn 3.

The parameters of these equations (α* _{m}*, α

_{p}, P_{0},

*L*

_{0}and λ

_{2}) 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.

*P*

_{0}and values of length (and hence velocity) by

*L*

_{0}. The only equation whose form is changed by this manoeuvre is Eqn 7, which becomes:

#### The combined model

*P*(in the SE) and

*P*

_{c}(in the CE) must be equal in the steady state, because the CE and SE are in series. However, setting these variables as exactly equal makes the calculation of

*ν*

_{c}by Eqn 6 extremely difficult, because it requires differentiation of Eqns 8, 9 and 10. Furthermore, in reality the stretch of the SE by the CE surely cannot be instantaneous. For these reasons the transfer of force is modelled here by simple linear kinetics:

*k*

_{5}is chosen large enough (see Table 1) that

*P*is nearly identical to

*P*

_{c}. The full model thus consists of Eqns 1, 2 and 3, 5, 6 and Eqns 8, 9, 10 and 11.

#### Computational considerations

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 d*P*/d*t*, 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 d*P*/d*t* gives: where α_{1} and μ are as in Eqn 12.

#### WDD

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^{2+} 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 _{m}*

_{1}and

*k*

_{m}_{2}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*_{3} and *k*_{4}, the rate constants of the binding of Ca^{2+} 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*_{3} and *k*_{4} are no longer constant but are given by: where *k*_{30} and *k*_{40} are the respective values when *m*=1. The ratio *k*_{4}/*k*_{3} is thus proportional to *m*, and Eqns 16 must be incorporated into Eqns 1, 2.

The changing value of the ratio *k*_{4}/*k*_{3} 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_{m}_{1}and

*k*

_{m}_{2}) are used to fit the full model to the force trajectories of the ramp experiments. Provisional values for

*k*

_{1},

*k*

_{2},

*k*

_{30}and

*k*

_{40}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.

### RESULTS

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.

### Effects of non-constant SE

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.

### WDD

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.

### Sinusoidal predictions

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).

### Work done by muscle during swimming

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.

## DISCUSSION

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*_{1}—*k*_{5}) 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*_{0}. The values of both *P*_{0} and *L*_{0} 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^{2+} 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.

## Acknowledgements

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.

## LIST OF SYMBOLS AND ABBREVIATIONS

*C*total Ca concentration

*Ca*free Ca

^{2+}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*_{0}optimum length

*m*variable responsible for WDD

*P*force developed by SE

*P*_{c}force developed by CE

*P*_{0}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

- μ
_{0}intercept of the dependence of SE stiffness on

*Caf* - μ
_{1}gradient of the dependence of SE stiffness on

*Caf*