SUMMARY
Undulatory animal locomotion arises from three closely related propagating waves that sweep rostrocaudally along the body: activation of segmental muscles by motoneurons (MNs), strain of the body wall, and muscle tension induced by activation and strain. Neuromechanical models that predict the relative propagation speeds of neural/muscle activation, muscle tension and body curvature can reveal crucial underlying control features of the central nervous system and the powergenerating mechanisms of the muscle. We provide an analytical explanation of the relative speeds of these three waves based on a model of neuromuscular activation and a model of the body–fluid interactions for leech anguilliformlike swimming. First, we deduced the motoneuron spike frequencies that activate the muscle and the resulting muscle tension during swimming in intact leeches from muscle bending moments. Muscle bending moments were derived from our videorecorded kinematic motion data by our body–fluid interaction model. The phase relationships of neural activation and muscle tension in the strain cycle were then calculated. Our study predicts that the MN activation and body curvature waves have roughly the same speed (the ratio of curvature to MN activation speed ≈0.84), whereas the tension wave travels about twice as fast. The high speed of the tension wave resulting from slow MN activation is explained by the multiplicative effects of MN activation and muscle strain on tension development. That is, the product of two slower waves (activation and strain) with appropriate amplitude, bias and phase can generate a tension wave with twice the propagation speed of the factors. Our study predicts that (1) the bending moment required for swimming is achieved by minimal MN spike frequency, rather than by minimal muscle tension; (2) MN activity is greater in the midbody than in the head and tail regions; (3) inhibitory MNs not only accelerate the muscle relaxation but also reduce the intrinsic tonus tension during one sector of the swim cycle; and (4) movements of the caudal end are passive during swimming. These predictions await verification or rejection through further experiments on swimming animals.
INTRODUCTION
Neuromuscular dynamics, linking the central neural networks and the mechanical body, plays a crucial role in controlling rhythmic movements during undulatory swimming. The timing of the neural activation of muscle and the resulting muscle tension in the strain cycle have been investigated to uncover the mechanisms for activation of and energy supply by muscle (powering mechanisms). It has been found that the phase relationship between tension and strain varies along the body for many species of fish so that midbody muscles generate most of the power required to maintain swimming, and the energy, transmitted through body elasticity, is dissipated into the water at the tail region (Hess and Videler, 1984; Cheng et al., 1998; Chen et al., 2011b; Altringham et al., 1993). This variation in the phase relationship between tension and strain results, in part, from the specific speed of the neural activation wave. Electromyographic (EMG) recordings have shown that activation waves travel along the body of a fish at a much slower speed than modelpredicted tension waves, but slightly faster than observed body curvature waves (Grillner and Kashin, 1976; Williams et al., 1989; Hess and Videler, 1984; Wardle and Videler, 1993; Van Leeuwen et al., 1990; Rome et al., 1993; Wardle et al., 1995).
The activation mechanisms leading to the different speeds of the traveling waves are not completely understood and only a few studies have been directed towards an understanding of the underlying mechanisms. To gain insight into the relative timing of EMG and body curvature waves, Bowtell and Williams simulated a mechanical model for lamprey swimming, containing inertia and viscoelastic properties but no fluid force effects; they imposed a forcing function to generate active tension (Bowtell and Williams, 1991; Bowtell and Williams, 1994). These researchers found that the speed difference between the active tension and body curvature waves was a consequence of body viscoelasticity. In particular, high damping and low stiffness resulted in slower body curvature waves. Because the hydrodynamic forces were absent in the model, these were not responsible for the differences in wave speeds. McMillen and colleagues refined the model by incorporating the dynamical properties of active tension development and hydrodynamic forces, and identified three primary factors responsible for the wave speed difference: viscoelastic damping, body taper and the nonlinear dependence of muscle force on length and shortening velocity (McMillen et al., 2008). In those studies, the timing of total tension (sum of the passive and active components) in the strain cycle was not mentioned; however, timing is important for determining the powering mechanisms of muscle along the body. It also remains unclear how relatively slow activation waves along the body generate much faster tension waves.
In this paper, we give an analytical explanation for the differences between the propagation speeds of neural activation, muscle tension and body curvature waves, based on two previous modeling studies of swimming movements in the leech Hirudo verbana Carena 1820: one explored a model for neuromuscular activation (Chen et al., 2011b) and the other provided a model for body–fluid interactions (Chen et al., 2011a). In the latter study, experimental measurements of body kinematics were used to predict the muscle bending moment during swimming. We used the neuromuscular model to predict the motoneuron (MN) activation spike frequencies and the resulting muscle tensions that generated the predicted bending moment and observed body curvature waves. The MN activation and body curvature waves were found to have roughly the same speed, whereas the tension waves travel much faster (about twice as fast). We argue, based on the premise that the speed of the curvature wave has evolved for efficient swimming, that the body–fluid interaction dynamics require that the tension wave propagates much faster to achieve this efficient curvature (strain) wave. Consequently, the speed of the MN activation wave arose to yield the tension needed under the influence of the strain change. A simple mathematical calculation of the product of two slow waves (MN activation and muscle strain) is shown to generate a tension wave that propagates at twice the speed of the strain and activation waves.
MATERIALS AND METHODS
Control of muscle activation by MNs
The shape of the leech body during swimming is ribbon like, about 10 cm long, 1 cm wide and 0.3 cm thick. The body undulates in the vertical plane with rearward traveling waves. Central interneuronal circuits (the central pattern generator, CPG) in midbody ganglia drive the MNs that, in turn, activate the dorsal and ventral segmental longitudinal muscles, whose rhythmic contractions lead to dorsal and ventral bending. Tension in dorsal or ventral segmental muscle is controlled by two sets of bilaterally symmetrical MNs: three dorsal excitatory MNs innervate the contralateral dorsal longitudinal muscle, and three ventral excitatory MNs activate the ventral longitudinal muscle – one innervates ipsilateral muscle, the other two innervate contralateral muscle (Kristan et al., 2005). There are corresponding inhibitory MNs that actively reduce the muscle tension commanded by the excitatory MNs (Mason and Kristan, 1982). Dorsal inhibitory MNs are active in antiphase to the dorsal excitors; similarly, ventral inhibitory MNs fire out of phase with the ventral excitors (Kristan et al., 2005). The mechanisms by which the inhibitory MNs regulate tension during swimming are unknown. We simply attributed the negative value of predicted MN spike frequency to the effect of the inhibitory MNs. Our experimental results show that when two dorsal excitatory MNs, DE3 and DE5, activate the muscle simultaneously, the resulting tension is greater than the sum of tensions generated by these neurons individually; more specifically, the tension generated by simultaneous activation was about three times that evoked by the individual MNs (Tian, 2008). We account for this interactive effect by multiplying the sum of spike frequency from all excitatory MNs by a factor of 1.5, and use this value as the input to the muscle model. This linear summation of MN frequencies captures the interactive effect, under the assumption that the excitatory (or inhibitory) MNs on one (dorsal or ventral) aspect of the body always burst together during swimming.
Relationship between MN spike frequency and muscle bending moment
We derive MN spike frequency from muscle bending moments by utilizing the tension model developed in our previous work (Chen et al., 2011b). The muscle bending moment was predicted from the measured kinematic motion data from video recordings by utilizing a body–fluid interaction model (Chen et al., 2011a). In our model, segmental units of the CPG are distributed in 17 midbody ganglia; overall, we modeled the body by a chain of 18 rigid links connected by frictionless rotational joints. At each joint (or body segment), the muscle bending moment was related to muscle tension, strain and MN spike frequency by: (1) where u(t) is the muscle bending moment at the ith joint (the subscript i is omitted for simplicity), x(t) is the muscle strain, T(t) is the muscle tension, f(t) is the comprehensive representation of spike frequencies of MNs in the dorsal or ventral body wall, subscripts d and v represent the dorsal and ventral muscle, respectively, and r is the moment arm of the tension and is approximated by half of the body thickness.
The thickness varies along the body (thicker in the midbody), and the measured thickness data for Leech 11 in Chen et al. (Chen et al., 2011a) was used for the analysis in this paper. The length of the moment arm may change over time because contraction of dorsal muscle would deform the body shape to make the arm length shorter/longer on the ventral/dorsal side. This detail is ignored in our model and the moment arms on the dorsal and ventral sides are assumed to have the same constant length. Our sensitivity analyses (not reported here in detail) indicated that the arm length varies within 10% during swimming if the muscle volume is assumed to be constant, resulting in small perturbations (within 10%) in the magnitude of strain, tension and MN spike frequency, and almost no visible change in the phase when compared with the predictions by the simpler model reported here. Thus, the general conclusions of this paper are not sensitive to these variations in moment arm.
The tension model is described by Chen et al. (Chen et al., 2011b): (2) where μ is the passive tension at nominal swim length (0.27 gram force), η is the length constant of the tension model (1/29), s is the Laplace variable, α is the gain of muscle activation (190 Hz^{–1}), τ_{1} is the time constant of muscle activation (0.23 s) and τ_{2} is the time delay in muscle activation (0.13 s). The parameter μ was doubled [μ=0.135 gram force in our previous paper (Chen et al., 2011b) described the tension of half of the dorsal or ventral body wall] such that T is the tension of the entire dorsal or ventral body wall within one segment, P(s) is the transfer function of muscle activation dynamics, and h(x) describes the model for passive tonus tension. We used the notation P(s)f in Eqn 2 to represent the timedomain signal L^{–1}{P(s)f(s)}, where L^{–1}{·} is the inverse Laplace transform.
The dorsal and ventral tensions and MN spike frequencies were calculated by solving Eqns 1 and 2 using the predicted bending moment and measured strains. There is a redundancy in T_{d}(t) and T_{v}(t) because increasing both by the same amount leads to the same value of the bending moment. This redundancy was eliminated by assuming that the animal achieves a desired bending moment with a minimal activation cost measured by the magnitude of the MN spike frequency. Once T_{d}(t) and T_{v}(t) are obtained, the MN spike frequency can be solved from Eqn 2 given that the only other unknown variable x (muscle strain) is available from video recordings of the swimming leech.
Optimization problem formulation
We now give the mathematical formulation for solving the MN spike frequencies f_{v}(t) and f_{d}(t) from the available data for u(t), x_{v}(t) and x_{d}(t) by utilizing Eqns 1 and 2. The following procedure was applied to each of the 17 segments of the leech body. We rewrite Eqn 2 as: (3) where we used a(f) to replace P(s)f+1, and * stands for d or v, representing the dorsal or ventral muscle. We will refer to h(x) and a(f) as the strain factor and the muscle activation factor, respectively. With muscle strain x_{*} fixed from the measured joint angle, T_{*} is a linear function of a_{*}. We express the discretized time courses T_{d}(t), T_{v}(t), u(t), a_{d}(t), a_{v}(t), h_{d}(t) and h_{v}(t) as column vectors of dimension n, the number of discretized time instants, and denote the vectors by T_{d}, T_{v}, u, a_{d}, a_{v}, h_{d} and h_{v}. Then Eqn 3 can be written in matrix form as T_{*}=H_{*}a_{*}, where H_{*} is a diagonal matrix with the diagonal being the entries of h_{*}. Replacing T_{v}(t) and T_{d}(t) in the first equation of Eqn 1 by their matrix form expression T_{*}=H_{*}a_{*} and u(t) by the vector u, we get a system of n linear equations with 2n unknowns (i.e. the entries of vectors a_{d} and a_{v}). To eliminate the redundancy in the solution, we look for the solution a_{d} and a_{v} such that the quantity: (4) is minimized with the additional constraints: (5) where t_{1} is the duration of the recorded swimming episode. In skeletal muscle, the activation heat can be approximated by a linear function of MN spike frequency (Hatze and Buys, 1977; Davy and Audu, 1987). Thus, the quantity defined by Eqn 4 is related to heat dissipation. The MN spike frequencies f_{d}(t) and f_{v}(t) in Eqn 4 are related to a_{d}(t) and a_{v}(t) by the dynamic function a(f)=P(s)f+1, where P(s) is a firstorder lowpass filter with a time delay, as given in Eqn 2. To express Eqn 4 in terms of a_{d}(t) and a_{v}(t), we need the derivatives of a_{d}(t) and a_{v}(t), which are approximated by [a(t+Δt)–a(t)]/Δt. The discretized Eqns 1–5 lead to a minimization of a quadratic cost function over Euclidean vectors a_{d} and a_{v}, subject to linear equality and inequality constraints, which can be solved efficiently using, for instance, the MATLAB linear matrix inequality (LMI) control toolbox. The signals f_{d}(t), f_{v}(t), T_{d}(t) and T_{v}(t) are then obtained from a_{d}(t) and a_{v}(t) through the relation a(f)=P(s)f+1 and Eqn 3. The spike frequency of each MN is estimated by dividing the result of f_{d} or f_{v} by 4.5 (3 MNs times the interactive summation factor 1.5).
RESULTS
Predicted MN spike frequency and muscle tension
Muscle strain and bending moment data for a typical swim episode of a midsize leech were taken from episode 11 of Chen et al. (Chen et al., 2011a) for our analysis, and are plotted in Fig. 1, where the strain is proportional to the joint angle of the linkchain model under the linear approximation. The colored curves correspond to five different body segments (or link joints of the body model), as indicated in the legend. The MN spike frequency and muscle tension on the dorsal and ventral aspects were predicted from Eqns 1 and 2, and are plotted in Figs 2 and 3 with the same color scheme as in Fig. 1. Figs 2 and 3 both show more than one full cycle of MN spike frequency and tension, and they do not repeat exactly from cycle to cycle. One reason is that the leech swimming movements are not repeated precisely from cycle to cycle either, as seen in the time course of the measured strains in Fig. 1A. The difference in the magnitude of MN spike frequency in consecutive cycles (Fig. 2B) may not be small but is in the reasonable range. The magnitude of spike frequency on the dorsal aspect approximates the observed values of 80–100 Hz recorded from the swimming of nearly intact leeches (Yu et al., 1999), and the values on the ventral aspect (∼100–200 Hz) seem too large. We get the spike frequency for single MNs based on the data that there are three pairs of excitatory MNs each on the dorsal and ventral aspects (Kristan et al., 2005); however, there may exist additional unidentified excitatory MNs. Although somewhat high, we view the predicted values for the ventral excitatory MNs as acceptable. We extract the amplitude and phase information from the time course data in the next sections.
The negative value of spike frequency implies an accelerated tension reduction. We attribute the negative spike frequency to the action of inhibitory MNs, firing out of phase with the excitatory MNs and reducing muscle tension (Mason and Kristan, 1982). Fig. 3 shows that tensions are near zero during a portion of the cycle. The predicted time courses of the MN spike frequencies and tensions suggest that the inhibitory MNs not only accelerate the muscle relaxation but also reduce the intrinsic tonus tension to nearly zero during a portion of the cycle.
Amplitude distribution along the body
MN activities at the head (J17) and tail (J1) joints are very low, as shown in Fig. 2. The amplitudes of the oscillations in the corresponding CPG needed to drive the MNs would need to be small as well. To compare the MN activation strength along the body quantitatively, we calculated the root mean square (r.m.s.) values of the excitatory MN spike frequencies (Fig. 2) over the firing interval, i.e. the interval with positive frequencies, and plotted the result in Fig. 4A. The r.m.s. values at the midbody (joint 8 to joint 13) are ∼100–120 Hz for those exciting the ventral muscles and ∼50–70 Hz for those exciting the dorsal muscles. The magnitude of excitatory MN spike frequency decreases from the midbody towards the head and tail, and is greater on the ventral aspect than on the dorsal aspect. The greater activation of the ventral muscle is perhaps due to the gravity effect; that is, the body bends more strongly towards the ventral aspect to generate lift force. The density of the leech body was measured for four leeches whose mass ranged from 0.73 to 3.83 g, and was found to be 1.065±0.0022 g cm^{–3} (Chen et al., 2011a). Simulations of the body–fluid interaction model for leech swimming (Chen et al., 2011a) indicated that a biased muscle bending moment generates the hydrodynamic lift force to compensate for the negative buoyancy, consistent with our observation of the MN frequency. The mean tension over one cycle was calculated and is plotted in Fig. 4A. Tension is consistent with MN activation data; that is, the mean tension over one cycle decreases from the midbody towards the two ends, and the tension on the ventral aspect is larger than that on the dorsal aspect.
Phases of MN activation and tension in the strain cycle
The peak sequence of midbody MN spike frequency (arrows in Fig. 2) indicates a rostrocaudally propagating wave. Similarly, Fig. 3 shows traveling waves for dorsal and ventral tensions. As the peak locations for the anterior and posterior segments are further apart in time for MN frequency than for tension, the propagation speed of MN activation is slower than that of tension. For quantitative analysis, we calculated the phases of MN spike frequency and tension in the strain cycle relative to the time instant at which the dorsal strain at the head segment takes a maximum value (t_{0}=0.28 s from Fig. 1), using the data of predicted MN spike frequency, muscle tension and the measured strain. The phase of a periodic signal is defined to be 2π(t_{0}–t_{1})/T rad, where T is the cycle period and t_{1} is the time instant defined by the following property: for the region under the time course and above the time axis within one cycle, the area on the left of the vertical line at time t_{1} is equal to the area on the right (the equal area rule). The results are plotted in Fig. 4B, where the spike frequency of MNs near the head and tail ends are too small and noisy to quantify the phases accurately. The slopes of the phase curves are inversely proportional to the propagation speed of the corresponding waves. This figure shows that body curvature waves travel at a speed close to (about 84% of) that of MN activation waves traveling along the body. The tension wave speed is roughly twice that of the other two. The greater speed of the tension wave resulting from the MN activation wave is explained in the Discussion (see ‘Mechanisms underlying different wave speeds’).
DISCUSSION
Predicted MN activation and muscle tension in comparison to existing results
In fish as well as leeches, undulatory swimming is characterized by body waves traveling from head to tail. The origin of the waves is in the rhythmic neural signals sent from the CPG in the form of MN action potentials to activate muscle contractions. MN activation (spike frequency) and the resulting muscle tension are both viewed as traveling waves as well. The intriguing fact is that these waves do not necessarily travel at the same speed as the body wave. Neuromuscular dynamics convert the MN activation into tension and, at the same time, vary the phase of tension in the strain cycle along the body. The interactions between activation and muscle strain in tension development, as well as the varied muscle properties along the body in some species, lead to the observed different propagation speeds among MN activation, muscle tension and body curvature waves. As a result, the phase of tension in the strain cycle can vary along the body, leading to different roles of muscles for power generation and storage during swimming.
The relationship of MN output to the longitudinal muscles and swimming undulations was investigated in previous experiments on medicinal leeches. The occurrence of MN spikes (generated by dorsal excitatory MN cell DE3) was recorded from electrodes implanted at two spatially separated points in the midbody of nearly intact swimming leeches (Pearce and Friesen, 1984; Yu et al., 1999). In these experiments, the leeches were tethered by threads attached to the denervated head and tail suckers and suspended in a water tank. Body shape was monitored with video recordings. The intersegmental phase lags of MN activity found in these independent experiments were 11.8 deg per segment (Pearce and Friesen, 1984) and 13.8 deg per segment (Yu et al., 1999). The ratio of wave propagation speed between body curvature and MN activation was about 0.77 (0.74=11.8/16 in the former experiment and 0.80=13.8/17.3 in the latter). Our model predicted a speed ratio of 0.84, a relatively slower MN activation wave compared with these experiment measurements. The results for leeches are comparable to the EMG measurements in fish. For anguilliform swimming, like eels and lampreys, EMG activity propagates a little faster than the body curvature wave (Grillner and Kashin, 1976; Williams et al., 1989; Wardle et al., 1995). Although leeches are not closely related to anguilliform swimming fish, the swimming gait and the MN activation patterns are analogous in these species.
In our study, the tension wave was predicted to be about twice as fast as the activation wave. The higher speed of the tension wave results from the interference between muscle strain and MN activation, as discussed in the next section. Wardle and coworkers (Wardle et al., 1995) carefully discuss the use of EMG signals as a first approximation to estimate the timing of muscle force generation in a strain cycle, and our result suggests that such approximation can be inaccurate. There currently are no data on tension, measured or predicted, for anguilliform swimming to compare with our prediction. However, models of carangiform swimmers also predict a much faster (almost standing) wave of bending moment than the body curvature and activation waves (Cheng et al., 1998; Hess and Videler, 1984). The bending moment results from the difference between tensions of antagonistic muscles (left/right in these studies), thus the tension and bending moment waves travel at essentially the same speed.
Mechanisms underlying different wave speeds
Muscle tension is generated through multiplicative effects of MN activation and strain factor as in Eqn 3 and, hence, depends strongly on the phase relationship of the two. To understand how a fast tension wave can result from slower MN activation along the body, we did a simple calculation of the product of two sinusoidal waves traveling at the same speed. Let the activation signal at the ith body segment be approximated by a_{i}(t)≈α_{i}sin(ωt+β_{i})+γ_{i}, where ω is the oscillation frequency, α_{i} and β_{i} are the amplitude and phase of the fundamental harmonic component, and γ_{i} is the bias. We note that the slope of phase curve (β_{i}, i=1,..., 17) along the body is inversely proportional to the wave speed. In particular, it takes (β_{i}–β_{j})/ω seconds for the wave to travel from the ith segment to the jth segment. Thus, the (average) wave speed is ω/Δβ body lengths per second, or 2π/Δβ body lengths per cycle, which is clearly faster if the phase lag from head to tail, Δβ=β_{17}–β_{1}, is smaller. Similarly, let the time course of the strain factor h(x) in Eqn 3 be described by . Their product gives the tension: (6) where the amplitude r_{i} and phase θ_{i} of the fundamental harmonic component are given by: (7) with j being the imaginary unit, and where c_{i} contains the constant and higher order harmonic (2ω) terms. A simple calculation then shows that the phase of the resulting tension θ_{i} is given by: (8) where ∠(·) denotes the angle of a complex number in the polar coordinates, ρ_{i} is the ratio of normalized amplitudes of a_{i}(t) and h_{i}(t) and δ_{i} is the phase difference between a_{i}(t) and h_{i}(t). This equation indicates how the phase (and hence the speed) of the tension wave depends on the phase (or speed) of the activation and strain waves through the amplitude and bias parameters.
Eqn 8 explains how two slow MN spike frequency and strain waves generate a faster tension wave as shown in Fig. 4B. To see this, first note that the muscle activation dynamics induces the phase lag ∠P(jω)≈–217 deg at ω=3 Hz from the MN spike frequency to activation if the signals were perfectly sinusoidal. Hence, the phases of a_{i}(t) on the ventral side, obtained by applying the equal area rule to the solutions of a_{i}(t) from the minimization problem, lag the phase curve of MN spike frequency (red circle curve in Fig. 4B) by roughly 200 deg and are plotted in Fig. 5 (red curve). The phases of strain factor h_{i}(t) are the same as that of muscle strain (red dot curve in Fig. 4B) and are plotted in Fig. 5 (black curve). The parallel phase curves show that a(t) and h(t) share the same propagation speed. The phase of the tension wave θ_{i}, computed from Eqn 8, is plotted in Fig. 5 (blue curve). This curve is close to the phase curve of the predicted tension calculated by the equal area rule (i.e. the red cross curve in Fig. 4B, reproduced as the cyan curve in Fig. 5). The slopes of the tension phase curves are roughly half of the slopes of the phase curves for a(t) and h(t), implying the tension wave speed was twice as fast.
The underlying mechanism for faster wave generation can be seen in Eqn 8. Roughly speaking, the phases of the two waves a(t) and h(t) linearly decrease by 360 deg from head to tail with the same slope, keeping a constant distance δ_{i} of about 180 deg. If δ_{i}=π–ε with a small perturbation parameter ε>0 as in Fig. 5, the phase lag from head to tail is given by: (9) provided ρ_{17} and ρ_{1} are smaller and larger than one, respectively. We now see that the phase lag of activation (and strain) Δβ≈360 deg reduces by approximately 180 deg to yield the phase lag of tension Δθ≈180 deg or, equivalently, the two slow waves of speed 1 body length per cycle multiply to give the faster wave of speed 2 body lengths per cycle. Thus, it is possible to generate a tension wave that is twice as fast from slower MN activation and strain waves of the same speed if they are roughly antiphase to each other and the amplitude and bias of both are set appropriately. In general, muscle dynamics allow for the generation of a tension wave that propagates at a different speed from the activation or strain wave.
Our results show that multiplicative effects of MN activation and muscle strain can generate a different wave speed for tension. It should be mentioned that multiplication between activation and muscle strain is not the only mechanism that generates a faster tension wave. Addition of passive and active tensions can also lead to a faster total tension wave (Cheng et al., 1998). A harmonic analysis similar to the one given above indicates that a faster tension wave can result when the passive and active tensions are roughly antiphase, and the ratio of the two amplitudes is distributed over the body appropriately. Our result on the leech muscle can also be explained from this point of view when the tension in Eqn 2 is interpreted as the sum of active tension h(x)P(s)f and tonus tension h(x). Computations show that the active tension h(x)P(s)f is (nearly) antiphase with the passive tonus tension h(x)P at all body segments. Thus, both active and tonus tension have the same traveling speed as the body curvature waves. The active tension tends to cancel the tonus tension through inhibitory MN activation when the muscle is stretched, and the remainder (total tension) turns out to have faster traveling waves.
Energy consumption along the body
The energy consumed by muscle has two main components: the cost of activation and mechanical work output. The former is mainly associated with heat dissipation resulting from biochemical processes, while the latter is the time integral of the product of tension and shortening velocity. The work output is positive/negative when the muscle in tension contracts/extends. Hence, the timing of tension generation during the strain cycle is important for understanding the different roles played by muscles distributed over the body. The functionality of muscles varies along the body of swimming animals. In the saithe, a carangiform swimmer, for instance, the midbody muscles are the power generator, and the muscles near the tail behave like a passive spring to transmit this power to the caudal fin (Altringham et al., 1993). The phase relationship between muscle tension and strain in our study of leeches predicts the powering mechanisms for anguilliform swimming.
The tension and strain in the tail region of the swimming leech were found to be nearly in phase (Fig. 4B). Together with the prediction that MN activity in the tail region is very small (Fig. 4A), this suggests that leech tail movements are almost passive. Shifting the phase curve of strain in Fig. 4B downward by 90 deg gives the phase of muscle shortening velocity. In the midbody, muscle shortening velocity is in phase with tension, thus generating positive instantaneous power over the whole cycle. This inphase property may indicate a resonance because the force and velocity of a mechanical system are in phase when it oscillates at a resonance. In contrast, the shortening velocity of muscle in the tail (head) region lags (precedes) the tension by 90 deg, making the average power output over one cycle from tail and head muscles close to zero. Therefore, swimming is powered mainly by the midbody muscles. This result is consistent with our previous prediction based on the body–fluid interaction model (Chen et al., 2011a). In addition to the favorable phase relationship between the shortening velocity and tension in midbody muscles, the greater tension amplitude (Fig. 4A) further assists in generating more power for swimming.
The quantity defined by Eqn 4 that we minimized to eliminate the redundancy in the antagonistic tension pair is related to the heat dissipation associated with active muscle contraction. In skeletal muscle, heat dissipation is composed of activation heat, maintenance heat and shortening heat (Davy and Audu, 1987). There are positive correlations between activation heat and MN spike frequency (Hatze and Buys, 1977), between maintenance heat and muscle tension (Gibbs and Gibson, 1972), and between shortening heat and shortening velocity (Hill, 1938). The fact that the minimization of an energyrelated quantity has led to reasonable predictions indicates that the energy cost may be important for leeches when swimming around in search of food. Our result predicts that the activation heat is larger in the midbody than in the head or tail because of the larger amplitude of the MN frequencies (Fig. 4A). This property is shared with the mechanical work output, but is not automatically expected because the activation heat depends on the amplitude alone, whereas the work output depends not only on the amplitude but also on the phase timing, which varies along the body.
To assess sensitivity with respect to the choice of cost functions, we also tried minimizations of the activation variable a_{*}(t) and the muscle tension T_{*}(t). However, these optimizations did not give reasonable predictions for MN spike frequencies – the peak magnitude of the negative spike frequency (i.e. inhibitory MNs spike frequency) turned out to be over 400 Hz for the midbody segment J9, for example. Thus, leeches appear to generate the bending moment required for swimming with minimal MN spike frequency, rather than with minimal activation or tension. If leeches indeed optimize energy efficiency during swimming, the result may imply that the activation heat is much larger than the maintenance heat of leech muscle and the MN spike frequency is more closely related to the total heat dissipation than the muscle tension or the activation variable. In addition, the phase of MN spike frequency was found to be insensitive to the quantities we minimized. Thus, overall conclusions on the wave speeds and their implications are robust against the choice of the cost function for minimization.
Perspectives on locomotion control mechanisms from the systems view point
The research reported here is part of a larger effort towards understanding the control mechanisms underlying undulatory swimming. Fig. 6 is a diagram of the feedback control system for leech swimming. We have developed a mathematical model for each dynamical component shown in this figure based on experimental data, including the CPG (Zheng et al., 2007), MN spike adaptation (Tian et al., 2010), passive tension in the body wall (Tian et al., 2007), muscle tension development by MN activation (Chen et al., 2011b) and body–fluid interactions (Chen et al., 2011a). The current study estimated the MN spike frequency and muscle tension during swimming, based on the previously developed muscle model, predicted bending moment and measured strain. In addition to the mechanisms for energy consumption and wave speed setting discussed earlier, the results can further be interpreted in the context of the whole swimming system as follows.
The body of swimming leeches exhibits approximately one spatial cycle of a quasisinusoidal body wave at each time instant, with the peak and trough traveling down the body at the speed of one body length per cycle. We have found that the MN activation wave has about the same propagation speed as the strain wave, but that the tension wave travels roughly twice as fast. How the leech chooses a particular speed for each wave has not been fully answered and therefore remains an open question. Approaches based on integrated models (Ekeberg and Grillner, 1999; McMillen et al., 2008; Tytell et al., 2010) would allow for sensitivity analysis of wave speed with respect to system parameters, but the mechanisms of speed determination may be hidden behind the model complexity. However, the analysis of each component in Fig. 6, including the results reported here, suggest the following. The body shape, or the speed of strain waves, is first determined by the fluid dynamics around the body of a particular geometry to achieve efficient swimming. This statement requires further justification, but appears reasonable (Taylor, 1952; Lighthill, 1960). Assuming that this is true, then body–fluid interaction dynamics require the tension (or bending moment) wave be much faster to achieve the prescribed curvature (or strain) wave (Hess and Videler, 1984; Cheng et al., 1998; Chen et al., 2011a). Consequently, the speed of the MN activation wave is set (through animal evolution) to yield the required tension under the influence of the strain change in accordance with the muscle activation dynamics.
The leech CPG is influenced by muscle tension through sensory feedback signals generated by tension receptors. This feedback modulates the oscillation pattern (frequency, amplitude and phase) of the membrane potentials to generate the appropriate MN activation pattern (Yu et al., 1999; Kristan et al., 2005). Without sensory feedback, the wave of membrane potentials in the CPG propagates down the isolated nerve cord at the same rate as the fast tension wave. In the intact animal, with sensory feedback, segmental CPG membrane potentials propagate at the speed of the much slower body curvature wave. These properties suggest that the leech CPG adaptively modifies the oscillation pattern through sensory feedback. In contrast, this adaptive pattern formation property does not seem to be shared by the lamprey CPG for anguilliform swimming under nominal conditions. In lamprey, previous studies (Wallén and Williams, 1984; Williams et al., 1989) have shown that the oscillation pattern of the isolated spinal cord is about the same as that of intact swimming animals, suggesting that the sensory feedback does not alter the CPG oscillation pattern. However, a more recent study (Guan et al., 2001) found that the phase lag was smaller for a semiintact preparation than for an isolated spinal cord preparation, and the studies in Tytell and Cohen (Tytell and Cohen, 2008) have clearly shown entrainment of CPG activities to mechanical bending. Thus, the lamprey CPG would also adapt its oscillation pattern through sensory feedback from edge receptors (Grillner et al., 1984) under perturbed conditions, as suggested by Ekeberg and Grillner (Ekeberg and Grillner, 1999). Detailed comparative studies of the seemingly different CPG control mechanisms underlying similar undulatory swimming of leeches and lamprey (and other fishes) would help us gain deeper understanding of the principles underlying pattern formation.
FOOTNOTES

FUNDING
This work was supported by the National Science Foundation [grant no. 0654070 to T.I.], the Office of Naval Research [MURI grant no. N000140810642 to T.I.] and the National Institutes of Health (NIH)/National Institute of Neurological Disorders and Stroke (NINDS) [grant no. 1 R01 NS4605701] as part of the National Science Foundation (NSF)/NIH Collaborative Research in Computational Neuroscience Program. We gratefully acknowledge ‘Interim Funding’ from the University of Virginia. Deposited in PMC for release after 12 months.
LIST OF SYMBOLS AND ABBREVIATIONS
 a(f)
 activation factor for leech muscle
 a_{d}(t)
 time course of activation factor of dorsal muscle
 a_{i}(t)
 activation factor at the ith body segment
 a_{v}(t)
 time course of activation factor of ventral muscle
 a_{d}
 column vector of discretized a_{d}(t)
 a_{v}
 column vector of discretized a_{v}(t)
 a_{*}
 column vector of discretized a_{d}(t) or a_{v}(t)
 CPG
 central pattern generator
 DE3
 dorsal excitatory motoneuron 3
 DE5
 dorsal excitatory motoneuron 5
 EMG
 electromyographic
 f
 motoneuron spike frequency
 f_{d}
 net spike frequency of MN activating the dorsal muscle
 f_{v}
 net spike frequency of MN activating the ventral muscle
 h(x)
 strain factor in tension model
 h_{d}(t)
 time course of strain factor of dorsal muscle tension
 h_{i}(t)
 strain factor at the ith body segment
 h_{v}(t)
 time course of strain factor of ventral muscle tension
 h_{d}
 column vector of discretized h_{d}(t)
 h_{v}
 column vector of discretized h_{v}(t)
 h_{*}
 column vector of discretized h_{d}(t) or h_{v}(t)
 H_{*}
 diagonal matrix with diagonal being h_{*}
 LMI
 linear matrix inequality
 MN
 motoneuron
 P(s)
 transfer function in muscle activation model
 r
 half of body thickness
 r_{i}, θ_{i}
 amplitude and phase of the fundamental harmonic component of the product of h_{i}(t)a_{i}(t)
 s
 Laplace variable
 T
 tension of leech body wall
 T_{d}(t)
 time course of tension of dorsal muscle
 T_{v}(t)
 time course of tension of ventral muscle
 T_{d}
 column vector of discretized T_{d}(t)
 T_{v}
 column vector of discretized T_{v}(t)
 T_{*}
 column vector of discretized T_{d}(t) or T_{v}(t)
 u(t)
 time course of muscle torque
 u
 column vector of discretized u(t)
 x
 muscle strain
 x_{d}
 strain of dorsal muscle
 x_{v}
 strain of ventral muscle
 α
 static gain of the transfer function for muscle activation
 αi, βi, γi
 amplitude, phase and bias of activation factor in the tension model at the ith joint
 amplitude, phase and bias of strain factor in the tension model at the ith joint
 δ
 phase difference between a_{i}(t) and h_{i}(t)
 ε
 small perturbation parameter
 η
 length constant of tension model
 μ
 passive (tonus) tension at nominal swim length (i.e. x=0)
 ρ
 ratio of normalized amplitudes of a_{i}(t) and h_{i}(t)
 τ1
 time constant in muscle activation
 τ2
 time delay in muscle activation
 ω
 oscillation frequency
 *
 representing dorsal or ventral muscle
 © 2012.