## ABSTRACT

Experiments were performed on single-myotome preparations of lamprey muscle, to discover whether force developed by intermittent tetanic stimulation during imposed sinusoidal movement could be predicted by data collected from isometric and constant-velocity experiments. We developed a simple dynamic model consisting of a set of simultaneous ordinary differential equations with unknown parameters. Appropriate values of the parameters were found by fitting numerical solutions of the differential equations to data from the isometric and constant-velocity experiments. Predictions were made of the time course of force developed during imposed sinusoidal movement in which the phase between muscle shortening and tetanic stimulation was varied to cover the whole phase spectrum. The match between the predicted and recorded time courses was very good for all phases, and particularly for those phases that are seen during swimming in the intact animal.

## Introduction

Lampreys are anguilliform swimmers; they generate thrust from the water by passing a wave of lateral curvature down the body from head to tail. The swimming musculature consists of approximately 100 myotomes, and these are activated sequentially from just behind the gills towards the tail, alternating on the left and right sides of the body. The relative timing between activation and curvature during steady swimming is shown in Fig. 1 (Williams *et al*. 1989). The dashed lines along the body of the lamprey show regions of muscle stimulation during swimming. The solid bars on the graph in Fig. 1A represent electromyogram (EMG) activity measured at four different points along the left side of the body, plotted against time. The sinusoidal lines are a measure of the local curvature of the body, which gives a measure of muscle length. Just behind the gills (25 % body length), the muscle is shortening during most of the time that it is stimulated, but this relative timing changes with position on the body. Nearer the tail, more of the activation period occurs during lengthening than during shortening. This relative timing is illustrated in Fig. 1B as a changing phase delay between the beginning of muscle stimulation and the beginning of shortening. This change in phase can be expected to have a large effect on the magnitude and time course of the force developed, as has been shown to be the case in isolated muscle preparations from the bullrout *Myoxocephalus scorpius* L. (Altringham *et al*. 1990), dogfish *Scyliorhinus canicula* (Curtin and Woledge, 1993) and scup *Stenotomus chrysops* (Rome *et al*. 1993).

The motivation for the study reported here was twofold: first, we wished to investigate whether the force developed by muscle stimulated intermittently during continuous oscillatory movement could be accounted for by a simple model in which the tension generated by fully activated muscle depends only upon the instantaneous values of length and velocity of shortening, or whether it would be necessary to postulate another influence not manifest in isometric and constant-velocity experiments, such as dependence on acceleration or on previous stretch or shortening (Abbott and Aubert, 1952). Second, we hoped to develop a predictive set of differential equations to include in our current dynamic model of anguilliform swimming (Bowtell and Williams, 1994; Carling *et al*. 1994; Williams *et al*. 1995).

To approach these questions, we developed a simple kinetic model of force generation by muscle. We found values for the parameters of the model by fitting the model to data collected during isometric and constant-velocity experiments. We then used the model to predict the time course of the force generated during imposed oscillatory movement at different phase values, data not used in developing the model or determining the parameters. Finally, we compared the predicted time courses with the experimental data.

Preliminary results have been published in abstract form (Williams *et al*. 1997). In addition, related work has been performed on dogfish muscle (Curtin *et al*. 1998).

## Materials and methods

Experiments were performed on muscle from the Mississippi River lamprey *Ichthyomyzon unicuspis*. Animals were anaesthetised by immersion in a mixture of ice and tap water containing 3-aminobenzoic acid ethyl ester (MS222, Sigma), then decapitated. The saline contained (mmol l^{−1}): NaCl, 93.5; KCl, 2.0; MgCl_{2}, 1.8; CaCl_{2}, 2.6; Hepes, 20; glucose, 4.0 at pH 7.0 and was equilibrated with 100 % O_{2}. Muscles were kept at temperatures below 10 °C at all times.

Muscle preparations were dissected from the body musculature within 3 cm caudal to the last gill opening (at a position approximately 20 % along the body length). Each preparation consisted of fibres from one myotome, which measured approximately 3 mm from myoseptum to myoseptum. The cross-sectional area of each preparation was approximately 1 mm^{2}. The ends of the preparation were held in foil clips secured with cyanoacrylate glue (gel type).

The arrangement of fast and slow muscle fibres in the lamprey is not as in teleost fish; instead, every myotome consists of a core of fast fibres surrounded by a cortex of slow fibres (Teräväinen, 1971). In the present study, no attempt was made to stimulate the two muscle types differentially.

The preparation was mounted horizontally between a force transducer (custom-made using semi-conductor strain gauges attached to a glass beam) and a motor (Cambridge Technology, Inc., series 300B). Temperature was between 5 and 7.5 °C during the experiments. The preparation was electrically stimulated with 1 ms pulses. Motor position and force were recorded on a digital oscilloscope.

In each experiment, the relationship between stimulus strength and twitch response was investigated to establish that the stimulus strength was sufficient to activate all the living muscle fibres. Tetanic stimulation at 25 pulses s^{−1} gave a fused tetanus of maximal strength and was used in all other parts of the experiments.

When movement was imposed on the muscle, the total length change was always 0.25 mm, which was approximately 8 % of the slack length of the fibres. This value was chosen to encompass the estimated changes in length due to the maximum and minimum body curvatures that occur during swimming in the lamprey, estimated from the data of McClellan (1984). The relationship between length and isometric force was explored to identify a suitable range of lengths to use in contractions with movement. In general, the resting force of unstimulated lamprey muscle is high at the optimal length and increases steeply with length (N. A. Curtin, unpublished observations). To prevent damage to the muscle, we therefore used lengths on the rising portion of the length–tension curve. Isometric records were made at the shortest, middle and longest lengths used in contractions with movement.

## Sinusoidal movement

Three cycles of sinusoidal movement at 1 Hz were imposed on the muscle, and one tetanus per cycle. The muscle was allowed to recover for 3 min between stimulus trains. The duty cycle (tetanus duration/movement cycle duration) was kept constant at 0.36, the average value for EMG duty cycle in swimming lamprey (Wallén and Williams, 1984). In the lamprey, as in the eel *Anguilla anguilla* (Grillner and Kashin, 1976), there is little or no change in duty cycle over the body length.

The relative timing of the stimulation and the movement, or phase, was defined as the interval between the start of stimulation and the beginning of shortening expressed as a fraction of the cycle duration. This convention is similar to that used by Rome (e.g. Rome *et al*. 1993), but in the units used by Williams *et al*. (1989). This definition was chosen because it expresses the delay between the beginning of the causative factor (muscle activation) and the beginning of the response (muscle shortening) in units directly comparable with duty cycle (fraction of cycle duration).

### Constant-velocity movement

The muscle was stimulated tetanically for 0.36 s, with movement starting 0.04 s before the first stimulus. The fastest velocity was chosen to match the maximum velocity during the sinusoidal movement. Two other velocities, half and one-quarter of the fastest velocity, were also used. Both shortening and stretch were tested.

Records were made of the passive force produced by the unstimulated muscle during both sinusoidal and constant-velocity movement. The active force was found by subtracting this passive record from the total force recorded during movement with stimulation. In the experiment reported here, the passive force was of the order of 10 % of the maximum isometric tetanic force.

Records of isometric tension at the middle length were interspersed between those with sinusoidal movement in order to monitor any change in the ability of the preparation to develop force. The results presented here (see Figs 3–7) are from the experiment in which this value changed least. Over the first 2.7 h of the experiment, the peak isometric value increased by 9 %, and over the remaining 1.3 h of the experiment it decreased to 93 % of the initial value. No normalisation was performed in this experiment. In the worst case of the four, the peak isometric force at the middle length changed by approximately 35 % over the course of the experiment. To compensate, we scaled all the force values in this experiment by the ratio of the most recent isometric peak value to the first isometric peak value before fitting the parameters of the data. The comparison between the predicted and experimental sinusoidal records was qualitatively similar to that shown in Figs 6 and 7 but quantitatively poorer.

### Characterising the series elasticity

The preparation was stimulated tetanically for 0.33 s under isometric conditions and then a rapid step reduction in length was imposed followed by constant-velocity shortening. Ten different velocities of shortening were used. Various sizes of ‘step’ were tested to find the one that reduced force to match that produced during subsequent constant-velocity shortening. The force at the end of the step was plotted against the size of the step. The slope, found by linear regression, was taken as a measure of the stiffness of the series elasticity μ (in mN mm^{−1}).

## Results

Experiments were conducted on preparations from four animals. The results presented in Figs 3–7 are taken from the experiment in which peak isometric force at the middle length changed least over the course of the experiment, as described above. For all experiments, the general shape of the isometric, isovelocity and sinusoidal force traces at all phases were similar to those shown in Figs 3–7.

### Mathematical model

We have devised a simple kinetic model which contains calcium ions, sarcoplasmic reticulum (SR) and contractile filaments (CFs). Both SR and CFs can bind calcium ions. The rates at which calcium ions are bound and released are based on the principle of mass action (see Fig. 2). For example, the rate of binding of calcium ions to the CF is proportional to the product of the concentrations of free calcium ions and unbound filaments, with rate constant *k*_{3}. 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, we found it necessary to introduce some cooperativity in the release of calcium ions from the CF, 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}.

*t*is time. We assumed that the total amount of calcium in the system would be constant, as would the total number of bound and unbound SR and CF sites. If all concentrations are expressed in terms of the total volume of the preparation, these constraints apply to the sums of the concentrations as well. Thus: where

*C, S*and

*F*are the total number of calcium ions, sarcoplasmic reticulum sites and filament sites, respectively, divided by the volume of the preparation. (This use of the preparation volume is implicit; its value is not required.)

*P*

_{c}) proportional to the concentration of calcium-bound filaments, with a proportionality constant dependent on

*v*

_{c}, the rate of change of sarcomere length: where α(

*v*

_{c}) represents the force–velocity function. For simplicity, we assumed the function to be piecewise linear, with different slopes for

*v*

_{c}<0 and

*v*

_{c}>0, to be determined during the course of the analysis. Hence: where α

_{1}has different values for

*v*

_{c}<0 and

*v*

_{c}>0.

*P*

_{c}) to the force measured between the tendons of the muscle preparation (

*P*) and to relate the length and velocity of filament overlap to the length and velocity of the whole preparation, we used Hill’s (1938) mechanical model of muscle containing a contractile element (CE) in series with an elastic element. This gives: where

*l*

_{c}and

*v*

_{c}are the length and velocity of the contractile element,

*L*(

*t*) and

*V*(

*t*) are the predetermined length and velocity of the whole preparation, and μ is the stiffness of the series elastic element. In this paper, we take shortening velocities to be negative and lengthening to be positive, so that

*V*and

*v*

_{c}are the time derivatives of

*L*and

*l*

_{c}, respectively.

*F*with length. We know of no simple molecular model for the effect of shortening below the resting length of the thin filaments, so we chose to model the effect of length as an effect on the rate constants for calcium binding or release by the filaments,

*k*

_{3}and

*k*

_{4}. Our most successful attempt was to make their ratio dependent on length, as follows: where β(

*l*

_{c}) is a function to be determined.

*P*), six parameters [

*k*

_{1},

*k*

_{2},

*k*

_{3},

*k*

_{5},

*F*and α

_{1}], one unknown function [β(

*l*)] and two known functions [

*L*(

*t*) and

*V*(

*t*)]:

### Estimation of functions and parameters from isometric and constant-velocity data

*F*was chosen by assuming that the filament binding sites were approximately 99 % saturated at the maximum isometric tension. If α

_{0}of equation 12 is set equal to 1.0, so that α(

*v*

_{c}) in equation 11 is equal to 1.0 for zero velocity (with appropriate implied units), then: where

*P*

_{max}is the maximum force measured in a 0.36 s isometric tetanus at the longest length.

The function β(*l*_{c}) was derived from the force records during isometric stimulation, by assuming that by the end of the tetanic stimulation period (i) the binding of calcium ions to the filament sites had come to equilibrium (d[cf]/d*t*=0) and (ii) the velocity of filament sliding (*v*_{c}) had become zero.

*F*at all times gives: Setting

*v*

_{c}=0 gives, from equation 11: For a large enough

*k*

_{5}in equation 16, by the end of the tetanus: Combining equations 6–10, 21, 22 and 23 gives: For each of the three isometric lengths, the value of β(

*l*

_{c}) was determined from equation 24, where

*P*is the force obtained by the end of the tetanus (Fig. 3A) and

*F*is as given in equation 20. These values of β(

*l*

_{c}) were plotted against

*l*

_{c}obtained from equation 13. This is shown in Fig. 3B. Linear regression on these data gave the following function for this particular experiment:

*v*

_{c}) was determined from the data obtained with ramp stretches and releases. It can be seen in the experimental data (see Fig. 4A) that, at the end of the stimulation during a ramp or (for the fastest ramps) at the end of the ramp, the rate of change of force (d

*P*/d

*t*) had become constant and could easily be measured. For a large enough

*k*

_{1}, the amount of calcium remaining bound to SR ([cs]; see Fig. 2) will be negligible by this time. Solving equation 11 for α(

*v*

_{c}) with

*P*

_{c}from equation 16 and [cf] from equations 6, 10 and 21 gives:

*v*

_{c},

*l*

_{c},

*F*and β(

*l*

_{c}) are as in equations 14, 13, 20 and 25, respectively, evaluated at the time of measuring

*P*and d

*P*/d

*t*. At this point of the analysis,

*k*

_{5}was yet to be determined. Therefore, it was assumed in the first instance to be large enough that the term (d

*P*/d

*t*)/

*k*

_{5}in equation 26 was negligible compared with

*P*. This allowed values of α(

*v*

_{c}) at the end of each ramp to be determined from equation 26 and plotted against

*v*

_{c}, as in Fig. 4B. The initial estimate of the function α(

*v*

_{c}) was used to find a value for

*k*

_{5}(see below), after which new values of α(

*v*

_{c}) were determined from equation 26, and the final form of the function was determined.

*v*

_{c}) (Fig. 4B) reflects the underlying length–tension relationship, bearing in mind that positive values of velocity correspond to lengthening. The data were fitted by a piecewise linear regression with a break-point at α(0)=1. For the data from this experiment, the function was: The remaining parameters of the model,

*k*

_{1}–

*k*

_{3}and

*k*

_{5}, were found by a least-squares fit to the middle-length isometric data, as shown in Fig. 5. As described above, the fitting was carried out in two steps, using first the preliminary form of α(

*v*

_{c}) to find values for all four parameters (

*k*

_{1}–

*k*

_{3}and

*k*

_{5}). Using this value of

*k*

_{5}, the final form of α(

*v*

_{c}) was determined. A final least-squares fit, using this form of α(

*v*

_{c}) and the previous value of

*k*

_{5}, was used to give final values to the parameters

*k*

_{1}–

*k*

_{4}.

This determination of the parameters of the model by finding a least-squares fit of a set of ordinary differential equations to experimental data was accomplished using the computer software *xpp*, developed by G. Bard Ermentrout and available free of charge from web site http://www.pitt.edu/.

### Predicting the sinusoidal data

The differential equations 17–19 with best-fit parameters were then used to predict the time course of force development during experiments with imposed sinusoidal movement. The range of phase values used spanned the entire cycle. Fig. 6 shows the time courses obtained at four of the eight phase values used. The solid lines represent the predictions, the dashed lines the data. Fig. 7 summarises the results for all phases used, comparing the predicted and observed values of peak force (Fig. 7A) as well as the time from the beginning of stimulation to the time of peak force (Fig. 7B). In Fig. 7B, the muscle was stimulated between the time represented by the origin of the *y*-axis and the horizontal dashed line.

## Discussion

In the present study, we have shown that a dynamic model based on simplified calcium kinetics and muscle filament properties can be used to predict the time course of force generation when muscle is stimulated during imposed movement. The data in Fig. 7 show that the peak force during sinusoidal movements and the time of its occurrence can be reasonably well predicted by this model. The time course of the development of force was also well predicted, including inflection points (as in Fig. 6, phases 0.1 and 0.7). The fit is especially close over the phase range 0–0.3 (see Figs 6, 7), a range that encompasses the phase values seen during swimming (see Fig. 1B).

Our model contains 11 parameters that have been fitted from the isometric and isovelocity data. It seems likely that there could be other models that could fit the data as well or better, and equally likely that such a model would predict the sinusoidal data as well or better. Without an exploration of the possible models, it is not possible to judge how well this particular model reflects reality. The conclusions of this study, however, do not depend on how well the model reflects reality, but only on its predictive power, which has been shown to be high.

Clearly, the calcium kinetics have been over-simplified, with only one intracellular calcium store in the model and the implicit assumption that diffusion is instantaneous. However important these features may be for understanding force development at the biochemical or molecular level, we believe that this work has shown that they are not necessary for understanding the physiology of swimming muscle at the systems level.

One of the most significant findings of this study is that it was not necessary to include any effects of acceleration on force development. The experiments used for developing the model and fitting its parameters were performed under isometric or constant-velocity conditions, when acceleration is zero. The model was then tested by predicting the force developed while the preparation was subjected to sinusoidal lengthening and shortening, in which case there is significant acceleration. Yet the experimental data were well predicted by the model. Thus, any effects of acceleration on force development must be insignificant compared with the effects of length and velocity changes. Furthermore, it was not necessary to include any ‘memory’ of previous position or velocity.

In experiments such as these, movements are imposed by a servomotor. In contrast, the changes in length and velocity of a given myotome in intact muscle during swimming come about in response to forces developed by the water, the inert tissues of the body and neighbouring myotomes, as well as the force generated by the muscle itself. By imposing movements on isolated preparations which mimic those during swimming, we are able to deduce the approximate time course of the sum of all the forces that oppose the action of the muscle. To understand swimming, we need to know more about these forces, and it is for this reason that we are developing a dynamic model of the system, including both body mechanics and fluid dynamics (Williams *et al*. 1995).

The next step in this part of the study will be to investigate muscle from several positions along the length of the lamprey body, to see whether the muscle properties change with position and, if so, to encode this variation in a position dependence of the parameters of the model.

## Acknowledgement

This work was funded by the Wellcome Trust and the BBSRC.

## References

*J. Physiol., Lond*

*J. exp. Biol*

*J. math. Biol*

*Mechanics and Physiology of Animal Swimming*

*J. exp. Biol*

*Scyliorhinus canicula*

*J. exp. Biol*

*Neural Control of Locomotion*

*Proc. R. Soc. Lond. B*

*in vitro*preparation of the lamprey brainstem/spinal cord

*Brain Res*

*Science*

*J. Neurophysiol*

*in vitro*compared with swimming in the intact and spinal lamprey

*J. Physiol., Lond*

*Biological Fluid Dynamics*

*SEB Symposium*

*J. Physiol., Lond*. (in press)

*J. exp. Biol*