Mechanisms underlying rhythmic locomotion: dynamics of muscle activation

Muscle contraction is the basis for body movements of animals. Tension is developed as a result of muscle stretch, motoneuron (MN) activation and the simultaneous effects of both. The tension of contractile element T has been modeled as a function of the MN impulse frequency f, muscle strain (or length) x and its time derivative (velocity) . Sliding filament theory (Huxley, 1957; Huxley and Simmons, 1971; Huxley, 1971; Huxley, 1974) suggests that muscle activation scales the tension–length–velocity properties of maximally activated (tetanized) muscle (Hatze, 1977; Zajac, 1989), yielding a multiplicative model structure T(f,x,v)=a(f)h(x,v), where 0<a(f)<1 is the muscle activation level and h(x,v) is the tension of tetanized muscle at length x and velocity v. The tetanic tension is often described by Hill's equation (Hill, 1938) or its variations of the form h(x,v)=p(x)q(v), where the length–tension curve p(x) and velocity–tension curve q(v) are typically determined from tetanic isometric and isovelocity experiments, respectively. The model for the muscle activation a(f) is found from the isometric contraction experiments where the activation level a is proportional to the measured tension because the length and velocity are fixed constants in isometric experiments; the function a(f) can thus be derived using the data at different activation levels (Lloyd and Besier, 2003; Brown et al., 1999).

Although the multiplicative models have been successful in explaining the data from standard isometric and isovelocity experiments under the tetanic condition, their predictive capability remains debatable under conditions different from those used for modeling. Muscle activation dynamics can be too complex to be captured by the simple multiplicative structure; for example, rapid shortening can deactivate the thin filaments, and strongly bound cross bridges can activate skinned fibers in the absence of calcium ions (Gordon et al., 2000). In fact, a number of studies have demonstrated that multiplicative models can fail to predict tensions accurately when the MN stimulation rate is low and/or there are continuous changes in muscle velocity (Perreault et al., 2003; Camilleri and Hull, 2005; Julian and Moss, 1981; Gordon et al., 2000). A common cause of the failure is the assumption that the activation, length and velocity factors in the multiplicative structure are independent of each other. In reality, the activation factor may depend on the velocity (Williams, 2010), and the velocity factor may depend on the activation level (Otten, 1987; Winters and Stark, 1985). The tension developed under such coupling effects may not be well explained by a model with the decoupled multiplicative structure.

In contrast, if a model is developed directly from experimental data under functionally relevant conditions, rather than the standard tetanic isometric/isovelocity conditions, the model may accurately capture the behavior of interest with the simple, decoupled multiplicative terms because the activation level may not have to span the whole range from low to tetanic, and the velocity may be confined in a relatively small range. Indeed, this expectation has been met in the case of periodic movements mimicking lamprey swimming (McMillen et al., 2008). Thus, direct modeling under functionally relevant conditions can be important for accurately characterizing the targeted behavior with least model complexity. However, no systematic methods are currently available for such modeling processes and one has to rely on brute-force parametric optimizations that could suffer from computational difficulties. Moreover, such parametric optimizations have to a priori assume a model structure (multiplicative decoupling, in particular) without verification.

In this paper, we focus on rhythmic movements where MN activation and length change are both periodic, sharing the same frequency with a particular phase relationship between the two. We propose a systematic method for testing the (decoupled) multiplicative structure, i.e. T(f,x,v)=a(f)h(x,v), based on the singular value decomposition of an appropriately constructed tension data matrix. With decoupling confirmed, dynamical models for the activation function a(f) and tension–length–velocity function h(x,v) can then be derived. The method is applied to the longitudinal muscle of leeches that are used for undulatory swimming. The input–output data (f,x,v) and T were collected from novel ‘dual-sinusoid’ experiments on nerve-cord and body-wall preparations where the body wall was activated by rhythmic MN impulse bursts and at the same time subjected to the sinusoidal length changes. The experimental approach and the theoretical method for testing the multiplicative structure are general, and can be employed to model the active tension of the contractile element under periodic (but not necessarily sinusoidal) movements of vertebrate skeletal muscles, provided the series elasticity effect can be neglected as discussed elsewhere (Zajac, 1989).

This study is part of our larger effort to develop an integrated analytical model of the leech swimming system. We have developed mathematical models of the central pattern generator (CPG) (Zheng et al., 2007; Zheng, 2007), impulse adaptation in MNs (Tian et al., 2010), passive properties of longitudinal muscles of the body wall (Tian et al., 2007; Tian, 2008) and body-fluid interactions during undulatory swimming (Chen et al., 2011). The muscle model with MN activation developed here will be integrated with these component models later to study the feedback control mechanisms underlying animal locomotion. Our focus is therefore more on systems-level modeling for accurately describing the rhythmic muscle contraction behavior than on uncovering physiological mechanisms underlying tension development.

Medicinal leeches (Hirudo verbana Carena 1820) were obtained from Leeches USA (Westbury, NY, USA) and maintained in aquaria at ∼20°C in a light-controlled room on a 12 h:12 h light:dark cycle. Experiments were performed on preparations of nerve cord and body wall excised from medium-sized leeches, ∼11 cm long when gently stretched. The dissected body wall consisted of a 6-annuli-long (∼6 mm) rectangle, the width of which was half of dorsal sector excised from the mid-body. The length was one annulus longer than one segment, to reflect the innervation pattern of the MN, and xtended from the dorsal to the lateral midline. This dorsal half of the body wall is innervated by the dorsal posterior (DP) and the posterior posterior nerves (Stuart, 1970). Two dorsal excitatory MNs, DE-3 and DE-5, were present as bilateral homologs. Their axons cross the midline to synapse with muscle fibers contralateral to their somata. The nerve cord and body wall were pinned dorsal-side up in a glass-bottomed recording dish. The body wall was suspended between a tension transducer (Harvard Apparatus, Holliston, MA, USA) to record body wall tension and a servomotor actuator (model 94757, Airtronics Inc., Anaheim, CA, USA) to rhythmically change body length (Fig. 1). Serotonin is an important neuromodulator in leeches, which modulates muscle tension (Tian et al., 2007). Hence, in this investigation, we utilized a 10 μmol l^–1 concentration of serotonin to emulate physiological hormonal conditions.

To mimic the dual inputs of MN impulse activation and muscle length changes contributing to muscle tension during swimming, we injected sinusoidal current into identified excitatory MNs and simultaneously imposed sinusoidal length changes to the body wall (‘dual sinusoid’ experiments; Fig. 1). Because of MN innervation patterns (Stuart, 1970), preparations were activated by stimulation of either DE-3 or DE-5. Sinusoidal current was delivered to the soma of DE-3 or DE-5 through sharp glass capillary microelectrodes. The evoked impulse frequencies of the stimulated MNs were estimated from spikes recorded by extracellular suction electrodes in peripheral nerves that carry the axons of contralateral MNs. It had previously been shown that, because of electrical coupling between homologs, impulse frequency in one MN can be estimated by multiplying the impulse frequency in its homolog by a factor of 3 (Tian, 2008; Tian et al., 2010). The amplitude of the injected current ranged from 2 to 3 nA to evoke MN impulse bursts with impulse frequencies similar to those recorded in swimming leeches (Yu et al., 1999). The sinusoidal length change imposed on the body wall was ±8% of the nominal body length during swimming in accordance with measured kinematic data (Chen et al., 2011; Kristan et al., 1974). Nominal length was measured prior to dissection at the lateral mid-line of the segment in video frames of the swimming leech. The cycle frequencies of the two sinusoidal inputs (at 0.5, 1, 2 and 3 Hz) were identical and spanned the range observed in swimming animals. We varied the phase of sinusoidal length change with respect to the current injected into MN soma from 0 deg to 330 deg in 30 deg increments; thus, the phase differences between two inputs spanned one full cycle of 360 deg. In addition to the dual-sinusoid experiments described above, we also conducted experiments to measure passive muscle tension. In these experiments, the length of isolated flaps of body wall was varied sinusoidally to span ±8% of the nominal length; there was no activation of MNs. Body wall tension was detected by the tension transducer. All data collected were amplified, displayed and recorded with Powerlab hardware and Chart 5 software (ADInstruments Inc., Mountain View, CA, USA). Signal sampling rate was 4 kHz.

The modeling method described below may be applicable to contractile elements in both vertebrate and invertebrate muscles. To avoid confusion between studies of leech, invertebrate, muscle and vertebrate skeletal muscle, it should be noted, first, that parallel and in-series elastic elements are absent in leech muscle. Leech longitudinal muscle, whose rhythmic relaxation and contractions cause the body undulations, extends the length of the animal, with overlapping muscle fibers several millimeters in length (Sawyer, 1986). There are no tendons to generate in-series elasticity. The leech epidermis is slack within the operating length during swimming (Miller, 1975; Tian et al., 2007), hence there is no parallel elasticity either. Consequently, the tension we measured in dual-sinusoid experiments corresponds to the tension generated by the contractile element in Hill-type models. Second, leech muscle has intrinsic tonus tension at rest, in the absence of activation (Miller, 1975; Mason and Kristan, 1982; Tian et al., 2007). Because of the similarities of the ultrastructure of leech longitudinal muscles and those of other animals (Lanzavecchia et al., 1985; Paniagua et al., 1996), and the parallel between resting tonus in the leech and unstimulated tension in smooth muscle (Fung, 1993), it is reasonable to assume that the tonus tension results from the resting concentration of calcium ions. The activation of excitatory MNs would then increase the calcium ion concentration and thus increase the muscle tension from this resting tonus value. It appears reasonable to hypothesize that the tonus tension shares the same mechanisms as the active tension induced by the MN activations. Hence, the activation variable described in this paper combines the intrinsic tonus value and the activation by MNs. The intrinsic tonus tension is also called the ‘passive tension’ because it is generated without MN activation. This should not be confused with the passive tension in vertebrate muscles arising from parallel elastic elements. When applied to contractile elements of skeletal muscles, for which series elasticity can be neglected, our method would simply yield a model with no tonus effect.

We will consider a general situation where the time course data of tension T, strain x, and stimulation f have been generated from N experiments as follows. In the first experiment, x(t) and f(t) are varied periodically (with cycle period p) and the resulting tension T₁(t) is measured. In the second experiment, the tension T₂(t) is measured when the stimulation f(t) remains the same but the strain is time-shifted from x(t) to x(t+τ) where τ=p/N. Repeating this procedure N times using x(t+kτ) for the (k+1)th experiment, N sets of tension data T_k+1(t) (k=0,...,N–1) are collected. In our dual-sinusoid experiments, x(t) is a sinusoid and f(t) is the MN frequency resulting from sinusoidal current injection into the identified MN DE-3. The number of experiments is N=12, corresponding to the experimental procedure wherein the phase of strain with respect to the current injection was increased by k×30 deg in the (k+1)th experiment.

The muscle is viewed as a dynamical system with strain x and stimulation f as inputs and tension T as output. If the effects of x and f on T are multiplicative, the tension can be expressed as T(x,f)=a(f)h(x), where h(x) is the length factor, and a(f) is the activation factor due to stimulation. The functions h(x) and a(f) are possibly dynamic; this means that, if the outputs of functions h(x) and a(f) are denoted by h(t) and a(t) with a slight abuse of notation, h(t) may depend on x(t) as well as its derivatives (velocity, acceleration, etc.) and its past history, and similarly for the relationship between a(t) and f(t). Note that this h(x) represents the same quantity as h(x,v) in the introduction section, without indicating its dependence on v explicitly.

A method is described below to test whether the measured tension data can be explained as the product of length and activation factors. The basic idea is to divide the time courses of h(t) and a(t) for one cycle of the first experiment (i.e. k=0) into N equal intervals (corresponding to N shifts of strain with respect to current injection to MN). Then, if the multiplicative structure is correct, tension data from the remaining N–1 experiments can be reproduced by circulating N intervals of h(t) with respect to a(t), one interval at a time, and multiplying the circulated sets of h(t) by a(t). The validation of the multiplicative model structure is thereby transformed into checking whether the computed values of tension for the remaining N–1 experiments replicate the measured data.

With the multiplicative structure, tension data recorded in the (k+1)th (k=0,...,N–1) experiment is described by T_k+1(t)=h(t+kτ)a(t).

The matrix T(t) admits a factorization of the form h(t)a(t)^T if, and only if, T(t) has rank one. Therefore, the tension T(t) can be expressed as the product of strain factor h(t) and activation factor a(t) if, and only if, the tension data matrix T(t), defined by Eqn 2, has rank one for all time t. By periodicity, it suffices to check the rank of T(t) for the interval 0≤t≤τ. Suppose the tension measurements are made at discrete time instants t_i for i=1,...,m in this interval, and are repeated with the same time spacing in the subsequent intervals kτ≤t≤(k+1)τ for k≥1. Then, the multiplicative effect can be tested by calculating the rank of m tension data matrices T(t_i). Because the measured data always contain some noise, the rank will never be exactly one. In practice, one can compute the singular values of T(t_i), and if the second largest singular value is substantially smaller than the first, T(t_i) can be regarded as effectively rank one. In that case, we conclude that the effects of the strain and MN impulse frequency on the tension generation are multiplicative.

Upon completion of the factorization process for each cycle frequency, we had a data set of time course (f,a) for muscle activation and (x,h) for the length–tension relationship. The next step was to find system models that relate the input–output data.

For the pair (x,h), one can first determine if the length–tension relationship is static by plotting the data on the (x,h) plane. If the data points lie on a single curve, then the relationship is static (i.e. the output h at time t is determined by the input x at time t, and is independent of the past history of x). In this case, a model can be simply obtained by expressing h as a linear combination of some basis functions of x (such as polynomials and exponentials) and calculating the coefficients through least square optimization. If the data plot on the (x,h) plane makes a loop, the relationship is dynamic. In this case, one may plot h as a function of (x,v), which can be visualized as a surface in the three-dimensional space (x,v,h), where is the velocity. The surface can then be modeled by a least square fit to some basis functions, such as Hill's equation. The length and velocity factors would be decoupled if the cross sections of the surface along each axis direction have the same shape after a normalization by some scalar.

The relationship between a and f is typically a dynamic mapping. One approach to finding a model is to assume an analytical structure for the relationship such as the calcium ion kinetic model (Williams et al., 1998) and to determine the parameter values by data fit. Another approach, taken here, is to assume linearity and fit the data by a transfer function. The transfer function is a description of the linear dynamical system that converts the time-domain integration and differentiation operations of the system on the input signal (e.g. f) to the division and multiplication of the Laplace transform of the input signal by the Laplace variable s. The output of the system is obtained by the inverse Laplace transform of the product of the transfer function and the Laplace transform of the input signal. By replacing the Laplace variable s with jω, where j is the imaginary unit and ω is the frequency of the sinusoidal input signal, the transfer function becomes a complex number, the magnitude of which is the amplification (gain) from input to output and the angle is the phase lag of the output from the input.

As mentioned earlier, the tonus tension is assumed to share the same mechanisms as the active tension induced by MN activation. Thus the activation factor in our model contains the combined effect of the tonus and MN activation. Assuming linearity, this gives the activation model structure a=P(s)f+b, where P(s) is a transfer function and b is a constant. Here, P(s)f should be interpreted as the output signal of the system P(s) in response to the input f(t), that is, P(s)f is the time-domain signal obtained by the inverse Laplace transform of the product of P(s) and the Laplace transform of f(t). Note that P(s)f captures the direct effect of MN impulses and equals zero when the MN activation is off. We chose b=1 so that h represents the passive tonus tension. For the contractile element of striated muscles having no passive tension, b should be set to zero.

Four data sets of a(t) and f(t) were obtained at ω=0.5, 1, 2 and 3 Hz, where f(t) at each cycle frequency was the average of the measured MN impulse frequencies over the 12 phase-shifted experiments. Using these data, we estimated the gain and phase of P(jω) (denoted by |P(jω)| and ∠P(jω), respectively) at the specified cycle frequencies. Neither a(t) nor f(t) were sinusoids, and we expanded them into Fourier series. Ideally, each harmonic pair of the Fourier series of a(t) and f(t) gives data points |P(jω)| and ∠P(jω) at every integer-multiple of the cycle frequency, generating a discretized frequency response plot. However, the Fourier analysis showed that the amplitude of the fundamental harmonic was much larger than those of the higher order harmonics. This suggested that the |P(jω)| and ∠P(jω) determined by the fundamental harmonics were more reliable than the ones by the higher order harmonics when the recording noise in a(t) and f(t) were considered. Hence, we only used the fundamental harmonics of a(t) and f(t) to calculate |P(jω)| and ∠P(jω).

The data sets of a(t) and f(t) at the four cycle frequencies gave rise to four points on each of the gain and phase frequency response plots. An additional datum point is the static gain P(j0)=(a₀–b)/f₀ at zero frequency, where a₀ and f₀ are the average values of a(t) and f(t) over the cycle. The static gain value should be the same for the four data sets (f,a) from different cycle frequencies. The four scaling parameters ρ in the factorization process were chosen to satisfy this constraint while minimizing the difference between the length factor h(t) and the observed passive tension data. The structure and parameters of P(s) were determined from the data points P(jω) at ω=0, 0.5, 1, 2 and 3 Hz. In addition to the frequency response plot of P(s), the shape of the time courses of a(t) and f(t) assisted in determining the structure of P(s), as explained later in the Results section.

Construction and evaluation of the tension model is based on the data set derived from one leech and illustrated in Fig. 2. This figure shows the body wall strain, the current injected into the MN soma, the MN impulse trains and the MN impulse frequency and tension data, from one dual-sinusoid experiment. Muscle strain oscillated between plus and minus 8% of nominal swim length (top trace). The injected sinusoidal current had an amplitude of 2 nA (second trace from top). Impulse trains in a DE-3 MN, electrically coupled with the penetrated DE-3 MN, were recorded from the DP nerve contralateral to the body wall nerve (middle trace). Impulse frequency was calculated as the reciprocal of interspike intervals. The impulse frequency of the target MN DE-3 was estimated to be three times that computed from extracellular recordings of impulses generated by the contralateral MN (Tian et al., 2010). The bottom trace is the tension resulting from the strain and MN input to the longitudinal muscles; maximum tension is ∼3 gram force (0.02958 N). In additional experiments on the same leech preparation, the sinusoidal strain advanced current injection by phases spanning from 0 deg to 330 deg in 30 deg increments, generating 12 time series of tension data. The complete data set for all 12 experiments, performed at four cycle frequencies (ω=0.5, 1, 2 and 3 Hz), are plotted in Fig. 3, left column. One conclusion is that peak tension varies with the phase difference between muscle strain and MN activation impulses. This relationship is most obvious at low cycle frequencies.

For each cycle frequency, the tension data were used to construct the tension data matrices T(t_i) (i=1,...,m) defined at each sampling time instant t_i according to the left hand side of Eqn 2. The largest singular value of T(t_i) was substantially larger than the second largest (the ratio of the two was more than 10 at all time instants t_i for all four cycle frequencies). This means, roughly speaking, that the modeling error due to the multiplicative structure assumption is less than 10% of the measured tension. Thus, the effects of muscle strain and MN impulse frequency on tension are experimentally verified to be multiplicative.

Following the modeling procedure, we factored the length factor h(t) and activation factor a(t) from T(t_i) (i=1,...,m). The products of a(t) and the time-shifted h(t+kτ) (k=0,...,11) are plotted in Fig. 3, middle column. The measured tension data (Fig. 3, left column) are reproduced accurately by the product of h(t+kτ) and a(t). Thus, the multiplicative model structure T(x,f)=h(x)a(f) is validated.

The activation factor a(t) and the length factor h(t) for the cycle frequency ω=1 Hz, as an example, are plotted in Fig. 4A,B (black curves). We see that h(t) has the magnitude of ∼1 gram force, which is in the range of recorded passive tension value, and a(t) oscillates above 1 with peak value near 3 gram force, indicating that the passive tension h(t) can be increased up to threefold through MN activation. Twelve time courses of the MN impulse frequency at 1 Hz are shown in Fig. 4C, and those at the other three cycle frequencies are similar. The maximum firing frequency is ∼100 Hz, which is consistent with the MN activity during intact swimming (Yu et al., 1999). The data f(t) are defined as the average of the measured twelve time courses.

The data for x(t) and h(t) are plotted parametrically on the strain–tension plane (x,h) in Fig. 5B. Although time is not shown explicitly, the data points progress around the narrow loops in the clockwise direction, once each cycle. Note that the tension curve obtained during body wall stretching is not exactly replicated during the contraction portion of the cycle. The implication is that the tension h is a dynamic function of the strain x, that is, h(t) depends not only on x(t) but also possibly on its derivative v(t). We could thus attempt to characterize h as a static function of x and v as in the Hill's equation h(x,v)=p(x)q(v). However, the area enclosed by each of the four loops is relatively small, with loops generated by the different cycle frequencies overlapping. Hence, it is reasonable to model the tension h as a static function of x, ignoring the small dynamic effects. Under this simplification, the tension can be modeled as h(x)=μe^x/ζ, a simple exponential curve, where the constant parameters μ and ζ are tuned to provide a good fit to the loops of all four cycle frequencies. The resulting values are μ=0.135 gram force and ζ=1/29; the model tension–strain curve is plotted in Fig. 5B (thick blue curve). Finally, the data from the dual-sinusoid experiments (Fig. 5B) are compared with the tension data obtained from the passive sinusoidal stretching experiments (Fig. 5C). The passive (tonus) tension h(t) extracted from the active tension data is remarkably close to the directly measured passive tension, consistently with our hypothesis that the tonus tension shares the same mechanisms as the tension induced by the MN activation, and the active tension is obtained by scaling the passive tonus tension.

The tensions predicted by Eqn 6 are plotted in Fig. 3 (right column), where the inputs are the recorded data of MN impulse frequency f(t) and muscle strain x(t). The predicted tensions are reasonably close to the measured data, capturing the peak amplitudes and their modulations via the phase shift between f(t) and x(t).

Many models assume a multiplicative structure T(f,x,v)=a(f)p(x)q(v) for the independent contributions of activation, length and velocity factors to tension development. The structure and parameters of each component function are then determined based on the observed muscle properties and the tension data from isometric and isovelocity experiments. Because of the complexity of muscle dynamics (Gordon et al., 2000; Julian and Moss, 1981; Balnave and Allen, 1996) such models can fail to predict the tension under functionally relevant conditions that are remote from the condition under which the model is developed. Perreault and colleagues modeled the cat soleus muscle and found that the error between the model prediction and data was consistently large for a variety of length change and stimulation patterns; sometimes greater than 50% during large length changes (Perreault et al., 2003). Based on their analysis, they concluded that the error could be attributed to the absence of coupling between the activation factor a(f) and the velocity–tension function q(v) in Hill-type models. Camilleri and Hull determined the parameter values of a Hill-type model to fit the data of rapid front-arm release movements at low muscle activation levels, and found that the model parameter v_max (shortening velocity at zero load) in function q(v) deviated considerably from the value determined from the tetanic isovelocity experiments (Camilleri and Hull, 2005). The dependence of v_max on the calcium ion concentration was also discussed in the earlier literature (Julian and Moss, 1981; Gordon et al., 2000). To compensate for this model error, v_max was parameterized by the muscle activation level in some applications (Otten, 1987; Winters and Stark, 1985), and the resulting models no longer had the decoupled multiplicative structure.

Our results indicate that the multiplicative structure holds for the longitudinal muscle of leeches under rhythmic movements that emulate the undulatory swimming condition. The reason that the length and activation factors are decoupled can be understood if we hypothesize that the rate of change of muscle length v during normal swimming is much smaller than v_max under no load. Recall that Hill's equation for the velocity–tension relationship is given by q(v)=(v_max–v)/(v_max+rv), where r is the ‘shape parameter’. The dependence of v_max on the activation level has been identified as a major factor that destroys the decoupled multiplicative structure, as discussed above. However, if the magnitude of v is much smaller than that of v_max, the velocity–tension function can be approximated as q(v)≅1, regardless of the value of v_max, which may depend on the activation level. In this case, the tension can be modeled as the product of the independent activation and length factors. In our study, the dependence of the tension on the velocity was indeed found to be small; the relationship between the length factor h and strain x was almost static (Fig. 5B). Further isovelocity experiments on leech muscle are needed to test the hypothesis that the muscle length change during swimming is much slower than the maximum shortening velocity.

Our study, like those of others, was motivated by the desire to understand the mechanisms that underlie animal locomotion through modeling of muscle properties under periodic stimulation and rhythmic movements. In pursuit of such understanding, Williams and colleagues addressed the fundamental question of whether a muscle model developed from isometric/isovelocity experiments could predict the forces generated during swimming (Williams et al., 1998). They answered the question in the affirmative by using an activation model that describes the kinetics of calcium ions with ordinary differential equations (ODEs). Some coefficients of the ODEs were dependent on the length, and the tension was modeled as the product of two factors: a length-dependent activation factor and a velocity factor.

A Hill-type multiplicative model was developed (McMillen et al., 2008) by removing the length dependence of the ODE coefficients and multiplying a separate length factor in the previous model (Williams et al., 1998). With the parameters of this model determined from isometric/isovelocity experiments, the multiplicative model was able to predict the tension under rhythmic activation or movements with a reasonable accuracy. This model was further refined to improve its accuracy (Williams, 2010) by incorporating the effect of work-dependent deactivation, which reduces tension following shortening more than would otherwise be expected. By having some of the calcium-ion-binding constants dependent on the velocity of muscle shortening, the refined model successfully predicted the tension observed during sinusoidal changes in length. The refined model had independent length and velocity factors, but the activation factor was dependent on the work done by the contractile element (the integral of the force–velocity product). It is interesting to note that the accuracy of the multiplicative model in McMillen et al. was greatly improved when model parameters were fit directly to experimental data during sinusoidal movements, without increasing the model complexity (McMillen et al., 2008). We encountered a similar situation when modeling the passive (tonus) tension of leech muscles (Tian, 2008; Tian et al., 2007). In those studies, the tonus tension model developed from step-stretch experiments was not able to accurately predict the tension induced by sinusoidal length change. Directly fitting the tension data derived from sinusoidal length changes to a simpler, static model gave a much better fit and greater predictive capability. This simple model predicted periodic tension changes for a range of cycle frequencies much better than the model based on step-stretch experiments. These previous findings demonstrated the importance of modeling data derived from functionally relevant conditions, and motivated the work reported here.

There are no systematic methods for modeling dynamics of tension development that are specifically tuned for rhythmic movements such as those observed in animal locomotion. The modeling work described above (Williams et al., 1998; Williams, 2010) was based on tetanic isometric/isovelocity experiments, and experimental tension data under sinusoidal movements and intermittent tetanic activation were used for validation purpose only. McMillen et al. are among the few researchers who have performed modeling directly based on tension data under sinusoidal movements (although this part of their modeling was not the focus of this reference) (McMillen et al., 2008). A model structure was fixed a priori, and least-square curve-fitting software was used to determine the model parameters. Although this particular case was successful, such a gray-box approach, when it fails, suffers from the inability to identify the crucial elements that limit model accuracy. For instance, if a muscle model fails to explain the data, it may be difficult to determine whether the failure is due to the multiplicative assumption or oversimplified calcium kinetics.

In this paper, we have proposed a systematic method for modeling muscle tension resulting from periodic activation and length change. The method first validates the multiplicative model structure, and then determines the time courses of activation and length factors from tension data for developing component models. The novel dual-sinusoid experiments were conducted on nerve-cord and body-wall preparations, where the muscle was activated through sinusoidal current injection into a MN to induce rhythmic bursts at realistic impulse frequencies, rather than tetanic electrical stimulation directly applied to muscle intermittently (Williams et al., 1998; Williams, 2010). The experiments were designed for the physiological conditions during swimming – rhythmicity and the phase relationship between muscle strain and activation – thus were specifically tuned for locomotion study. Periodic time-course data for the tension, length and MN impulse frequency were obtained for 48 cases (four cycle frequencies and 12 phase shifts between the length and MN frequency). Given such data, we have shown that the singular value decomposition of appropriately constructed tension data matrices can determine whether the tension is generated through multiplicative effects of the length and MN frequency (and their derivatives), without making any a priori assumptions. This method applies to the general problem of determining if a given output resulted from two independent inputs in a multiplicative manner, and is analogous to the well-established theory of linear system identification through factorization of Hankel matrices (Akaike, 1974; Van Overschee and De Moor, 1993).

Our study predicts that the tension developed within leech longitudinal muscle during swimming is highly dependent on muscle length, but much less on shortening and lengthening velocity. In contrast, for vertebrate muscles, dependence of tension on velocity tends to be more significant than on length as muscles normally operate within the plateau of the length–tension curve. The multiplicative structure could be destroyed under large dependence on velocity as discussed above, especially when the activation level varies over a wide range. In such cases it is important to test the accuracy of the structural assumption on the model, and our modeling method may be found useful in the context of rhythmic movements. However, care must be taken when the method is applied to muscle fibers with significant series elasticity, for then measured fiber length does not directly correspond to the length of the contractile element. Although Zajac made reasonable arguments for neglecting the series elasticity in muscle fibers, the problem remains of how to isolate the effect of series elasticity, when significant, on muscle tension developed under rhythmic movements (Zajac, 1989).

We have found that the length–tension relationship can be approximated by a static function (Fig. 5B), leading to a constant velocity factor in our model when viewed as a Hill-type model. This result is consistent with our previous study of passive (tonus) tension in leech longitudinal muscle (Tian et al., 2007); step-stretch experiments revealed that the relaxation time constants were larger than the swim cycle period and hence the muscle would act as a nonlinear spring with almost no damping (velocity-dependent) effect during swimming. The length–tension relationship in our model is an exponential function h(x)=μe^x/ζ. It should be noted that this h(x) represents the tonus tension, rather than the tetanic length–tension relationship in the standard Hill-type models for vertebrate muscles. The isometric tetanic length-tension relationship and the passive properties of leech longitudinal muscle were studied previously (Miller, 1975). Similar to skeletal muscle, the active length–tension curve had a ∩-shape, and the passive length–tension relationship was exponential. The difference between the shapes of the tetanic length–tension curve and h(x) identified in this paper suggests that the leech muscle property during swimming cannot be modeled accurately by scaling the tetanic tension from isometric/isovelocity experiments by the activation level.

Muscle activation in response to MN impulses was modeled by a first-order low pass filter (time constant 230 ms) and a time delay structure (time constant 130 ms). This time decay and delay should account for the processes of conduction of MN impulses to the end plate of muscle cells, calcium diffusion and the attachment and detachment of cross bridges (Huxley and Simmons, 1971). The large time decay and delay are not explained by impulse travel time in MN axons, which extend ∼5–10 mm from the soma to the longitudinal muscle fibers. Based on impulse-conduction velocities in leech interneurons, the impulse travel time should be no more than 30 ms, most likely considerably less. Synaptic delays in the leech are ∼5 ms (Granzow et al., 1985). So the time decay and delay mainly come from the muscle activation process, which includes calcium diffusion and the attachment and detachment of cross bridges.

 
• a(f)

  activation factor in muscle tension

 
• a(t)

  time course of activation factor in muscle tension

 
• a(t)

  vector related to a(t) defined by Eqn 2

 
• a₀

  constant term of Fourier series of a(t)

 
• c

  vector with its ith entry c_i

 
• c_i

  scaling parameter for relative magnitude

 
• CPG

  central pattern generator

 
• DE-3

  dorsal excitatory motoneuron 3

 
• DE-5

  dorsal excitatory motoneuron 5

 
• f

  motoneuron impulse frequency

 
• f₀

  constant term of Fourier series of f(t)

 
• h(t)

  time course of length factor in muscle tension

 
• h(t)

  vector related to h(t) defined by Eqn 2

 
• h(x)

  abbreviated notation for h(x,v) when the velocity effect is negligible

 
• h(x,v)

  length–velocity-dependent factor of muscle tension

 
• j

  imaginary unit

 
• J(c)

  quadratic cost function in terms of c to optimize the continuity of h(t)

 
• M

  positive semi-definite matrix such that J(c)=c^TMc

 
• MN

  motoneuron

 
• N

  number of experiments

 
• p

  period of dual-sinusoid experiments

 
• P(s)

  transfer function of muscle activation dynamics

 
• p(x)

  length–tension curve of tetanic isometric muscle contraction

 
• q(v)

  velocity–tension curve of tetanic isovelocity muscle contraction

 
• r

  shape parameter of q(v)

 
• s

  Laplace variable

 
• T(f,x,v)

  tension function

 
• T(t)

  tension data matrix defined by Eqn 2

 
• T_k(t)

  time course of tension data in the kth recording

 
• u_i

  left singular vector of T(t_i) corresponding to σ_i

 
• u_ij

  left singular vector of T(t_i) corresponding to σ_ij

 
• v

  muscle strain rate (time derivative of x)

 
• v_i

  right singular vector of T(t_i) corresponding to σ_i

 
• v_ij

  right singular vector of T(t_i) corresponding to σ_ij

 
• v_max

  shortening velocity at zero load of tetanized muscle

 
• x

  muscle strain

 
• α

  static gain of the transfer function for muscle activation

 
• ζ

  length constant of tension model

 
• μ

  passive (tonus) tension at nominal swim length (i.e. x=0)

 
• ρ

  magnitude scaling parameter

 
• σi

  largest singular value of T(t_i)

 
• σij

  jth singular value of T(t_i)

 
• τ

  time shifting of sinusoidal strain change in the consecutive experiments

 
• τc

  time constant in muscle activation

 
• τd

  time delay in muscle activation

 
• ω

  cycle frequency of dual-sinusoid experiments

 
• ∠P(jω)

  phase response of P(s) at frequency ω

 
• |P(jω)|

  magnitude response of P(s) at frequency ω

We gratefully acknowledge ‘Interim Funding’ from the University of Virginia. Deposited in PMC for release after 12 months.