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

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

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

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 Ca2+ 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 Ca2+

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

Mechanical model

The mechanical portion of the model begins with Hill's (Hill, 1938) formulation of muscle as consisting of a passive elastic element (SE) in series with a contractile element (CE). The CE is responsible for the active generation of force, subject to the length—tension and force—velocity properties of the muscle (Hill, 1938; Hill, 1949; Jewell and Wilkie, 1958). In this model, the length L of the muscle is given by the sum of the length lc 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:
(3)
(4)
where νc and V are the rates of change of lc 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.

Because of these findings, in our model we have made stiffness dependent on the level of activation. With this formulation, fully activated muscle has stiffness independent of force, as in the previous model (McMillen et al., 2008) whereas during force development and relaxation, stiffness is linearly dependent on the level of activation. Thus:
(5)
where μ0 and μ1 are the intercept and gradient of the dependence of SE stiffness on Caf.
In consequence, Eqn 4 becomes:
(6)

Physiological properties

As in the previous model (McMillen et al., 2008), the time course of force development by the CE is given by the product of the activation Caf and the effects of the length and velocity of the CE:
(7)
where the constant P0 is the peak force in an isometric tetanic contraction (when Caf=1) at the optimum length L0. In all four muscle preparations, the value of L0 was somewhat larger than the length in situ, Lis (see Table 1). The functions λ(lc) 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 λ(lc) and a piecewise linear function (approximating a hyperbola) for α(νc):
(8)
where lc0 is the length of the CE at the optimum length, L0, from Eqn 3.

The parameters of these equations (αm, αp, P0, L0 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≤λ(lc)≤1.

All calculations are made non-dimensional (except for time) by dividing values of force by P0 and values of length (and hence velocity) by L0. The only equation whose form is changed by this manoeuvre is Eqn 7, which becomes:
(10)

The combined model

The forces P (in the SE) and Pc (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:
(11)
where k5 is chosen large enough (see Table 1) that P is nearly identical to Pc. 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 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 α1 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.

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 Ca2+ 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 Pc 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 km1 and km2 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 k3 and k4, the rate constants of the binding of Ca2+ to the filaments and its release. This results in smaller values of Caf and, by Eqn 10, smaller force than would otherwise occur. Hence, k3 and k4 are no longer constant but are given by: where k30 and k40 are the respective values when m=1. The ratio k4/k3 is thus proportional to m, and Eqns 16 must be incorporated into Eqns 1, 2.

The changing value of the ratio k4/k3 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, km1 and km2) are used to fit the full model to the force trajectories of the ramp experiments. Provisional values for k1, k2, k30 and k40 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 lc 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 Pc 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, k1k5) 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 P0. The values of both P0 and L0 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 Ca2+ 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 Ca2+ concentration

  •  
  • Caf

    Ca-bound filament concentration

  •  
  • CE

    contractile component

  •  
  • F

    total filament concentration

  •  
  • lc

    length of CE

  •  
  • lc0

    length of CE at its optimum

  •  
  • L

    length of preparation

  •  
  • Lis

    in situ length

  •  
  • L0

    optimum length

  •  
  • m

    variable responsible for WDD

  •  
  • P

    force developed by SE

  •  
  • Pc

    force developed by CE

  •  
  • P0

    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 lc

  •  
  • V

    time derivative of L

  •  
  • WDD

    work-dependent deactivation

  •  
  • α(νc)

    dependence of Pc on νc

  •  
  • λ(lc)

    dependence of Pc on lc

  •  
  • μ

    stiffness of SE

  •  
  • μ0

    intercept of the dependence of SE stiffness on Caf

  •  
  • μ1

    gradient of the dependence of SE stiffness on Caf

Altringham
J. D.
,
Ellerby
D. J.
(
1999
).
Fish swimming: patterns in muscle function
.
J. Exp. Biol.
202
,
3397
-
3403
.
Corr
D. T.
,
Herzog
W.
(
2005
).
Force recovery after activated shortening in whole skeletal muscle: transient and steady-state aspects of force depression
.
J. Appl. Physiol.
99
,
252
-
260
.
Curtin
N. A.
,
Gardner-Medwin
A. R.
,
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
.
Edman
K. A. P.
(
1975
).
Fatigue vs shortening-induced deactivation in striated muscle
.
J. Physiol.
246
,
255
-
275
.
Edman
K. A. P.
,
Josephson
R. K.
(
2007
).
Determinants of rise time during isometric contraction of frog muscle fibres
.
J. Physiol.
580
,
1007
-
1019
.
Edman
K. A. P.
,
Nilsson
E.
(
1968
).
The mechanical parameters of myocardial contraction studied at a constant length of the contractile element
.
Acta Physiol. Scand.
72
,
205
-
219
.
Herzog
W.
(
1998
).
History dependence of force production in skeletal muscle: a proposal for mechanisms
.
J. Electromyog. Kines.
8
,
111
-
117
.
Herzog
W.
,
Leonard
T. R.
,
Wu
J. Z.
(
2000
).
The relationship between force depression following shortening and mechanical work in skeletal muscle
.
J. Biomech.
33
,
659
-
668
.
Hill
A. V.
(
1938
).
The heat of shortening and the dynamic constants of muscle
.
Phil. Trans. Roy. Soc. Lond. B
,
126
,
136
-
195
.
Hill
A. V.
(
1949
).
The abrupt transition from rest to activity in muscle
.
Phil. Trans. Roy. Soc. Lond. B
136
,
399
-
420
.
Hsu
C.-Y.
,
Tytell
E.
,
Fauci
L.
(
2009
).
An integrated muscle mechanic-fluid dynamic model of lamprey swimming
. http://meetings.aps.org/Meeting/DFD09/Event/112355
Hultmark
M.
,
Leftwich
M.
,
Smits
A. J.
(
2007
).
Flowfield measurements in the wake of a robotic lamprey
.
Exp. Fluids
43
,
683
-
690
.
Huxley
A. F.
,
Simmons
R. M.
(
1971
).
Proposed mechanism of force generation in striated muscle
.
Nature
233
,
533
-
538
.
Huxley
A. F.
,
Simmons
R. M.
(
1973
).
Mechanical transients and the origin of muscle force
.
Cold Spring Harb. Symp. Quant. Biol.
37
,
669
-
680
.
Jewell
B. R.
,
Wilkie
D. R.
(
1958
).
An analysis of the mechanical components in frogs striated muscle
.
J. Physiol.
143
,
515
-
540
.
Jewell
B. R.
,
Wilkie
D. R.
(
1960
).
The mechanical properties of relaxing muscle
.
J. Physiol.
152
,
30
-
47
.
Josephson
R. K.
(
1999
).
Dissecting muscle power output
.
J. Exp. Biol.
202
,
3369
-
3375
.
Josephson
R. K.
,
Stokes
D. R.
(
1999a
).
The force—velocity properties of a crustacean muscle during lengthening
.
J. Exp. Biol.
202
,
593
-
607
.
Josephson
R. K.
,
Stokes
D. R.
(
1999b
).
Work-dependent deactivation of a crustacean muscle
.
J. Exp. Biol.
202
,
2551
-
2565
.
Lichtwark
G. A.
,
Wilson
A. M.
(
2005
).
A modified hill muscle model that predicts muscle power output and efficiency during sinusoidal length changes
.
J. Exp. Biol.
208
,
2831
-
2843
.
McClellan
A. D.
(
1984
).
Descending control and sensory gating of fictive swimming and turning responses elicited in an in vitro preparation of the lamprey brainstem/spinal cord
.
Brain Res.
302
,
151
-
162
.
McMillen
T.
,
Williams
T. L.
,
Holmes
P.
(
2008
).
Nonlinear muscles, passive viscoelasticity and body taper conspire to create neuro-mechanical phase lags in anguilliform swimmers
.
PLoS Computational Biology
.
4
,
e1000157
.
Moréchal
G.
,
Plaghki
L.
(
1979
).
The deficit of the isometric tetanic tension redeveloped after a release of frog muscle at a constant velocity
.
J. Gen. Physiol.
73
,
453
-
467
.
Tytell
E. D.
,
Cohen
A. H.
(
2008
).
Rostral versus caudal differences in mechanical entrainment of the lamprey central pattern generator for locomotion
.
J. Neurophysiol.
99
,
2408
-
2419
.
Varkonyi
P. L.
,
Kiemel
T.
,
Hoffman
K.
,
Cohen
A.
,
Holmes
P.
(
2008
).
On the derivation and tuning of phase oscillator models for lamprey central pattern generators
.
J. Comput. Neurosci.
25
,
245
-
261
.
Williams
T. L.
,
Grillner
S.
,
Smoljaninov
V. V.
,
Wallén
P.
,
Kashin
S.
,
Rossignol
S.
(
1989
).
Locomotion in lamprey and trout: the relative timing of activation and movement
.
J. Exp. Biol.
143
,
559
-
566
.
Williams
T. L.
,
Bowtell
G.
,
Curtin
N. A.
(
1998
).
Predicting force generation by lamprey muscle during applied sinusoidal movement using a simple dynamic model
.
J. Exp. Biol.
201
,
869
-
875
.