Mechanisms underlying rhythmic locomotion: body–fluid interaction in undulatory swimming

Biological organisms living in an aquatic environment swim by generating a propulsive force through the interaction of rhythmic body movements with the surrounding fluid. There are a variety of swimming styles observed for organisms with different body geometry and internal morphological structure. For instance, carangiform swimmers such as carps constrain the large undulations to the caudal region with small undulations at the thick front part of the body, anguilliform swimmers such as eels have obvious traveling waves from the anterior part down to the posterior and mobuliform swimmers such as rays flap their flexible wings (large pectoral fins) up and down. Studies of locomotion mechanisms aim at understanding, among other questions, why a particular swimming style is chosen by an animal, how the thrust is generated through the dynamic body–fluid interactions and how the central nervous system controls the motion through sensory feedback.

As an essential step toward this understanding, rhythmic body movements during swimming have been analyzed through video recording for various animals including anguilliform (D'Août and Aerts, 1997; Gillis, 1996; Gray, 1933a), subcarangiform (tadpoles) (Wassersug and Hoff, 1985) and carangiform (saithe, mackerel) (Videler and Hess, 1984) swimmers. These studies calculated the kinematic parameters for undulation (frequency, amplitude and speed of traveling waves) and characterized swimming performance using such quantities as the stride length (UT/L₁), propeller efficiency (U/V) and Froude efficiency (E/W), where U is the swim speed, V is the traveling wave speed, T is the cycle period, L₁ is one body wave length in the swimming direction, W is the total work done by the body to the fluid over a cycle and E is the energy returned from the fluid to the body over a cycle to generate the thrust.

The mechanisms underlying propulsive force generation have been studied through exploration of the fluid dynamics surrounding the body. Experimental flow visualization techniques, such as particle image velocimetry (Muller et al., 2001; Tytell and Lauder, 2004; Tytell, 2004) and tracer dyes (Brackenbury, 2002; Brackenbury, 2004), were used to reveal the wake structure and understand, for instance, how vortices contribute to increase propulsive efficiency (Fish and Lauder, 2006; Muller and Leeuwen, 2006; Triantafyllow et al., 2000). Computational fluid dynamics (CFD) simulations complement these experimental studies by providing, in addition to flow visualization, quantitative estimates for the fluid forces acting on the body and the energy exchange between the body and fluid (Borazjani and Sotiropoulos, 2009; Carling et al., 1998; Cortez et al., 2004; Kern and Koumoutsakos, 2006). Although the experimental and computational approaches uncover details of the body–fluid interactions, analytical insights on locomotion mechanisms are hard to develop through these approaches because of the complexity of the experimental data or the computational model.

Classical works by Taylor (Taylor, 1952) and Lighthill (Lighthill, 1960; Lighthill, 1969; Lighthill, 1970; Lighthill, 1971) and their extensions (Cheng et al., 1991) provide analytical models of fluid forces acting on the body during undulatory swimming. The fluid forces were modeled as static functions of kinematic variables: velocity-dependent resistive force (Taylor, 1952) for low Reynolds number (Re) swimming or acceleration-dependent reactive force (Lighthill, 1960) for high Re swimming. These fluid force models have been combined with dynamical models of the body given by a chain of rigid links (Bowtell and Williams, 1991) or a continuum (Bowtell and Williams, 1994) to understand swimming mechanisms from the perspective of body–fluid interactions (Cheng et al., 1998; Hess and Videler, 1984; Jordan, 1996; McMillen and Holmes, 2006; McMillen et al., 2008). The combined body–fluid models enhance our understanding of how muscle activation leads to the observed movements under the influence of the fluid, providing predictions on, for instance, the speed of activation waves in comparison with the traveling body waves. Such analytical models facilitate further studies on swimming mechanisms, including optimal periodic motion for efficient locomotion (Blair and Iwasaki, 2011; Saito et al., 2002) and neuronal control by central pattern generators (Ekeberg, 1993; Ekeberg and Grillner, 1999).

In this paper, we develop a model for body–fluid interactions during undulatory swimming and analyze the mechanisms underlying thrust generation. The animals used in this study are leeches that swim by undulating the body like eels or snakes, except that the body oscillation occurs in a vertical (rather than horizontal) plane. Our model is a simple combination of Taylor and Lighthill's fluid models (Taylor, 1952; Lighthill, 1960) and a link-chain model for the body. Unlike most of the previously developed models, drag coefficients in our fluid model are determined through quantitative analysis of experimental motion data. Our results contribute to a more complete understanding of unsettled issues such as how the reactive and resistive hydrodynamic forces interplay in the anguilliform swimming of elongated animals such as eels or leeches. We have separated the thrust from drag in the total hydrodynamic force and examined the contribution of each body segment to thrust generation. The distributions along the body are estimated for thrust and drag, as well as for the energy supply from the muscle and dissipation into the fluid. The speed of traveling waves in the muscle bending moment is also estimated and compared with the speed of the undulatory body waves.

The results reported here are part of a larger effort to develop an integrated model (Fig. 1), supported by experimental data, for undulatory swimming. We have already developed component models for the leech swimming system, including the central pattern generator (CPG) (Chen et al., 2008; Zheng et al., 2007; Zheng, 2007), impulse adaption in motoneurons (Tian et al., 2010), and passive elastic dynamics (Tian et al., 2007) and motoneuron activation (Tian, 2008) for the longitudinal muscle. In the future, the body–fluid interaction model developed in this paper will be integrated with these models into a closed-loop control system with sensory feedback. Such an analytical integrated model is essential for a systems-level understanding of the neuronal control principles for animal locomotion.

Experiments were carried out on adult leeches, Hirudo verbana Carena 1820 [for several decades misidentified as Hirudo medicinalis (Trontelj et al., 2004; Siddall et al., 2007)] to collect kinematic data during intact swimming. Leeches swam in a transparent Plexiglas® tank (75 cm long, 3 cm wide, 10 cm water depth), which forced leeches to swim within a very narrow distance range from a video camera (which was mounted 1.5 m from the tank) to avoid the parallax error. Swimming behavior was initiated by mechanical tactile or electrical stimulation. The tank was long enough for leech to go into steady-state swimming halfway down the tank, the site of the capturing window for the camera. Swimming profiles were recorded by a video camera (PL-B774 USB color camera; 60 frames s^–1; PixeLINK, Ottawa, ON, Canada).

For kinematic analysis, we used only those video sequences in which the leech swam nearly horizontally at nearly constant speed for at least two cycles of undulation. In each video frame, numerical coordinates of boundary points for the body image were obtained and the midline of the body was calculated to describe the body shape, orientation and location (Fig. 2A). The midline was divided into 18 links defined by 19 points, where the 17 joints between the links correspond to the 17 mid-body ganglia. The coordinates of the 19 points in consecutive video frames give the discrete time courses of kinematic variables (link angles and the location of the center of mass) that were curve-fit by sinusoids or polynomials. The body shape is described by the joint angles (relative angles of adjacent links at each joint), and their harmonic approximations (Fig. 2B) characterize the undulation style in terms of the cycle frequency, amplitude, bias and phase. Fig. 2C shows the phases of joint angles along the body during one swim cycle. The phase lag from head to tail is roughly 360 deg, indicating that the waves traveling down the body have a wavelength that is roughly equal to the body length. The data were obtained for eight swimming episodes from five leeches.

During swimming, the leech uses its dorso-ventral muscles to flatten and elongate the body through internal hydrostatic pressure, which contributes to increasing the resistive/reactive hydrodynamic force on the body. The body dimensions of swimming (rather than resting) leeches were determined from captured video images. The body thickness was measured from those frames in which the lateral axis of the body was perpendicular to the camera screen (Fig. 2A). The thickness varies along the body, and the measured values were averaged from four snapshots for each leech. The body length was determined by averaging the length of the body midline over all frames in the episode. We used the video images in which the leeches swam with the dorsal or ventral side facing the camera to measure the body width. The chances of leeches swimming in this way were small, and only a portion of the body was captured, but we were able to measure the maximum body width of individual leeches around the mid-body. Variations of the width along the body were measured from a single representative medium-sized leech for which swimming episodes were not recorded but a top view image was available (Fig. 3). Because the width variations are similar for different leeches, the measured data were used for all leeches by proportionally scaling the data by the maximum width of each individual leech.

The total mass of the body was measured for the five leeches used in the video recording experiment. The distribution of mass along the body was calculated from the dimension data, assuming uniform density. We measured the body density through the use of sucrose solutions for another set of four leeches. For each leech, the concentration of sucrose dissolved in water was increased until the animal became neutrally buoyant. The density of the sucrose, and hence of the leech, was then determined by weighing a measured volume of the sucrose solution. This measurement was carried out twice on each of the four leeches, whose mass ranged from 0.73 to 3.83 g. The average density of all the leeches was used for the model.

The first equation shows how the muscle bending moment u and fluid torque τ result in the change in the body shape and orientation θ. The second equation describes the net effect of the gravity and buoyancy mg and fluid force h on the motion of the center of mass w.

We determined c_a, c_p and c_t through the best fit of leech swimming video data. The idea is to impose observed undulatory movements on the leech model, calculate the resulting motion with respect to the inertial frame by numerical simulation, compare the swim velocity with the observed velocity, and tune the fluid model coefficients to minimize the error. The equations of motion (Eqns 1a and 1b), given in terms of the variables θ and w, are not directly suitable for this process because the equation for θ contains the unknown muscle bending moment u. It turns out that the direct effect of u is to change the body shape, which in turn induces the inertial motion (body rotation and translation) through the body–fluid interactions. To impose the observed body shape change on the model, we need to decouple the shape dynamics from the orientation dynamics so that u is isolated from the equations of inertial motion.

In addition to quantitative estimates for the fluid coefficients c_p, c_t and c_a, this process predicts the time courses of the fluid force f from Eqns 2a and 2b, and the muscle torque u from Eqn 3a, where B(θ) is a square invertible matrix. The predicted values were used to analyze the swimming mechanisms in terms of the traveling waves along the body and energy flow from the body to the fluid.

The leech body is ribbon-like and is approximately 100 mm long, 10 mm wide and 3 mm thick during swimming. Fig. 5 shows the variations of body dimensions and mass along the body for five leeches. The body tapers towards both ends. The head is thin and narrow, reducing the fluid drag. The tail is also thin, but its width is larger than the head, and so could receive fluid forces for thrust generation. However, leeches do not have a caudal fin like eels or lampreys. Cross-sectional views of the leech body wall reveal that, except for blood in the central crop, longitudinal muscles comprise most of the body mass (Stuart, 1970), the oblique and circular muscles only occupy a small volume relative to the longitudinal muscle. The longitudinal muscles have the same orientation and anatomy along the body, thus the wider and thicker body in the middle portion indicates that the mid-body muscles are stronger and hence are able to supply more mechanical energy for propulsion. Numerical values of the dimension and mass data are summarized for the five leeches in Table 1. Heavier leeches are larger in all three dimensions. The ratios of maximum width and thickness to body length were calculated for the five leeches; mean (±s.d.) ratios were 0.083±0.0035 for width and 0.036±0.0031 for thickness. The standard deviations are relatively small, suggesting that the body shapes in three dimensions are similar for these leeches. The mean leech density was 1.065±0.0022 g cm^–3.

Fig. 6 shows a typical sample of body snapshots during one complete cycle of swimming (taken by a video camera at rate of 60 frames s^–1 and plotted every second frame). The right end is the head. The body exhibits about one quasi-sine wave in each frame. The wave travels backward (crest and trough move to the left) and the body progresses forward (to the right). In addition to the kinematic data, some characteristic parameters for swimming were calculated from the snapshot images. The propeller efficiency is defined as the ratio of the distance traveled in one cycle to the wavelength of the body (0.68 for the leech shown in Fig. 6). The wave number is defined to be the number of waves exhibited by the body at each time instant, and is calculated by (t₁–t₀)/(t₂–t₀), where a crest appears at the head tip at time t₀, travels down the body and arrives at the tail tip at t₁, and the next crest appears at the head tip at t₂. The time intervals were measured by the number of frames and the crest is identified by the maximum body curvature.

Table 1 summarizes some characteristic swimming parameters for leeches with masses ranging from 1 to 3 g. The snapshots in Fig. 6 gave the values for leech 2, episode 3 that are presented in Table 1. The propeller efficiency was ∼0.5–0.7, comparable to values for carp (0.76), mackerel (0.80), trout (0.69) and eel (0.7) (Videler, 1993). The swimming speed ranged between 0.13 and 0.21 m s^–1. The cycle period was between 0.3 and 0.5 s. The Re is defined by Re=ρUL/μ, where U is the body velocity relative to the fluid, and L is the length scale (body length or diameter). For a typical leech of body length L=0.1 m swimming in water at U=0.15 m s^–1, Re is 15,000. If we use the width 0.01 m for L and an average normal velocity of 0.05 m s^–1 for U, this results in an Re of 500.

The body undulations are described by the amplitude, bias and phase of the sinusoidal approximations of joint angles φ_i, which represent the body curvature (see Fig. 2). Variations of these oscillation parameters along the body are plotted in Fig. 7 for the five leeches. Fig. 7A shows that the oscillation amplitude of joint angles (or maximum curvature) gradually increases toward the tail. Fig. 7B shows the biases of the joint angles, which are correlated with the steering direction of locomotion (Saito et al., 2002); the positive (negative) bias at all joints makes the nominal body shape curl down (up) in the middle portion, and results in downward (upward) turning motion. The bias is roughly zero on average in Fig. 7B, indicating that the leeches swam straight forward. The negative bias near the head keeps the head up and may be related to the gravity effect (the leech density is slightly larger than the water density). Fig. 7C shows that the phase lag of the tail joint (joint 1) from the head joint (joint 17) is roughly 360 deg, meaning that it takes approximately one cycle period for a crest of a body wave to propagate from head to tail. The phase lag of x degrees from head to tail approximately corresponds to the wave number x/360. The larger (smaller) the phase lag, the more (less) waves expressed by the body, and the slower (faster) traveling waves for fixed cycle period and body length. Fig. 7C thus indicates that there is ∼0.8 to 1 body wave between the head and tail joints. Note that the wave number in Table 1 is the number of waves between the head and tail tips, and hence is slightly larger than the value read from Fig. 7C.

The fluid coefficients were determined from the swimming episodes 1, 2-1, 3-1 and 4 in Table 1, where the sizes of leeches 1–4 range from small to large. For fixed values of c_t and c_a, c_p was optimized to minimize the RMS error in the simulated swim velocity. Fig. 8A shows the contours of optimal c_p for different values of c_t and c_a. The optimal c_p decreases as c_a increases for a fixed c_t, and increases as c_t increases for a fixed c_a. These observations make sense if we note that the thrust is generated by the normal hydrodynamic force through the terms associated with c_p and c_a, and is balanced with the longitudinal drag due to c_t. To some extent, the roles of the resistive and reactive forces in the normal direction are interchangeable, although their phases are different. Fig. 8B shows the contours of the RMS error in the simulated velocity. The error decreases as c_a approaches zero, which suggests that the added-mass effect is small in undulatory leech swimming. The lowest RMS error happens at c_a=0, c_t=0.6, and the corresponding optimal c_p value is 3.0. The optimal value of c_p is much larger than the drag coefficient for the smooth circular cylinder, which varies from 0.9 to 1.1 when Re ranges in the interval 20<Re<10⁵. The RMS error is not very sensitive to the fluid coefficients near the region of c_a around zero.

The eight swimming episodes were simulated under the optimal fluid coefficients to verify that: (1) Taylor's resistive fluid force formula with the optimal fluid coefficients can reproduce the swimming behavior with a reasonable accuracy for episodes used for optimization (group O), and (2) swimming behaviors predicted by the resulting model are also close to the observations for episodes that were not used for optimization (those used for validation; group V) (see Table 1). The RMS errors in the swim velocity for the group O episodes were 1.7, 3.0, 2.9 and 4.1 cm s^–1 for episodes 1, 2-1, 3-1 and 4, respectively; those for the group V episodes were 3.4, 3.8, 2.8 and 4.4 cm s^–1 for episodes 2-2, 2-3, 3-2 and 5, respectively. The RMS errors for group V are larger than those for group O on average, but are still comparable.

To further evaluate the accuracy of the model, simulated swimming behaviors are plotted in Fig. 9, in comparison with the experimental video-recording data for the best and worst cases of the RMS error: episodes 1 and 4 for group O, and episodes 3-2 and 5 for group V. The simulations of the four swimming episodes show the following properties in common. The angular velocity of the body orientation (column B) was well matched between the simulation and data. The oscillations in v_x (the velocity of the center of mass in the x direction) were basically matched as well. There is a consistent phase advance of the simulated velocity v_y to the data. The simulations with the worst RMS errors are satisfactory, and we conclude that Taylor's resistive force model captures the hydrodynamic force with some accuracy. The optimal fluid coefficients obtained from four swimming episodes are applicable to the other four episodes.

We further calculated the optimal fluid coefficients for each swimming episode to double check the reliability of the optimal fluid coefficients. The results show that the optimal c_a is zero for all eight swimming episodes, and five out of eight swimming episodes (four out of five leeches) have c_p clustered in the range of 2.2–3.2 and c_t in the range of 0.4–1, which are around the optimal values (c_p=3.0, c_t=0.6). Leech 2 was the outlier, with larger values of optimal coefficients; c_p in the range of 4.8–7.6 and c_t in the range of 0.9–1.8 for three swimming episodes.

Fluid forces on the body were estimated by Eqns 2a and 2b through a simulation, using the optimal fluid coefficients, imposing the measured time course of body shape on the model. The fluid forces on the body are visualized in Fig. 10, where the black curves are the body midlines and the blue arrows are the fluid force vectors whose length indicates the magnitude of the force. The leech swims to the right. The sequence of nine frames makes up one cycle of undulation. The f_x and f_y values shown are the net fluid forces in the x and y directions, respectively, at the time instant of the snapshot.

Several observations can be made from the figure. The whole body undulates around the ‘nominal line’, which can be defined as the least square linear approximation of the body midline at each time instant¹. The nominal line has a positive pitch angle (roughly equal to θ_o), and this may be related to the fact that the vertical force f_y is always positive to balance the difference between gravity and buoyancy. The fluid force vectors on each link are pointed almost vertically, suggesting that the fluid is pushed away by the body mostly in the vertical direction in contrast to the rearward zigzagging jet flow observed in carangiform swimmers (Muller et al., 1997; Nauen and Lauder, 2002; Tytell and Lauder, 2004). For each body segment, the fluid force reaches the maximal magnitude when the segment is crossing the nominal line. The maximal force magnitude increases toward the tail. The sign of the net horizontal force f_x varies within the cycle; the leech body accelerates (f_x >0) when the tail is approaching the nominal line, and decelerates (f_x<0) when the tail is moving away from the nominal line. The leech's center of mass has (almost) no net acceleration or deceleration over one undulation cycle in the analyzed sequence as seen from Fig. 9, episode 1.

The longitudinal muscle of leeches is segmented along the body, and a bending moment is generated at each segment through differential contractions of the dorsal and ventral sides. The torque u_i applied at each body joint in our model accounts for this bending moment. The muscle torques were predicted from the body–fluid interaction model and the video-recorded motion data. Fig. 12A shows the time courses of muscle bending moments at several joints. For comparison, Fig. 12B plots the time courses of the measured joint angles (i.e. body curvature). Although the peak locations are spread out over the cycle for the curvature, those for the muscle bending moments are close to each other, indicating that the muscle actuation waves travel faster than the body curvature waves. The amplitude of the muscle torque is larger in the mid-body and becomes smaller toward both ends.

To quantify these observations, the muscle bending moment and joint angle are fit by sinusoids through Fourier analysis to describe the amplitude and phase of these two variables along the body. Because neither muscle bending moment nor joint angle is exactly periodic, each time course datum was fit for one cycle rather than for the whole recorded time duration. We visually checked the sinusoidal fittings of non-sinusoidal muscle bending moments to confirm that the first harmonic of the Fourier series captures the phase and amplitude reasonably well. Fig. 13 shows how the phase and amplitude of the muscle bending moment vary along the body. The phase of the joint angles is also plotted for comparison. The phase curves are almost linear, with positive slopes, indicating traveling waves along the body at almost constant speed. The plot also shows that the phase lag from head to tail for the muscle bending moment is much smaller than that for the joint angle. This means that the speed of actuation waves is two or three times faster than the resulting body waves. The amplitude plot indicates that the muscle bending moments are stronger in the mid-body and are weaker near the head and tail. Thus, the muscle actuation waves travel down the body faster than the body curvature waves, with increasing and then decreasing strength as the waves pass from the anterior to the posterior parts of the body.

Leeches generate swimming undulations by flattening their body by tonic activation of dorso-ventral muscles and phasic activation and relaxation of dorsal and ventral longitudinal muscles (Kristan et al., 2005). The energy, supplied through the bending moment generated by the longitudinal muscles to maintain steady swimming, may be temporarily stored in the form of kinetic energy in the body inertia or potential energy in the body elasticity. However, the energy will eventually be dissipated into the fluid as heat associated with the resistive force or as kinetic energy associated with the reactive force. To understand the energy transmission mechanisms, we have calculated the amount of power supplied and dissipated at each location of the body.

The instantaneous power supplied by the muscle at each joint is given by the product of the bending moment and the rate of change of the joint angle, . The time course of this quantity is shown in Fig. 14 at several joints along the body. The fundamental frequency of the power oscillation is twice as large as the cycle frequency, and has two crests and two troughs within one cycle. This is simply explained by the fact that the product of two sinusoids both at frequency ω has a component oscillating at 2ω. In the anterior half of the body (joints 9, 13 and 17), the supplied power is almost always positive. By contrast, the curves for the posterior body (joints 1 and 5) have large negative portions, indicating that the power is returned to the muscle during part of the cycle and is stored in the passive stiffness. Two troughs in one cycle correspond to the passive elastic energy of ventral and dorsal muscles, respectively. The more negative trough for joint 5 corresponds to the ventral muscle because at this time instant, ventral muscle is being stretched (see Fig. 12). The muscle power supply near the head is always close to zero whereas the power near the tail oscillates with an average close to zero.

We have modeled the leech body by a chain of rigid links and the hydrodynamic forces by static functions of kinematic variables. In particular, the fluid model includes both resistive and reactive forces. The resistive force is modeled as a function of relative velocity between body and fluid, and the reactive force is estimated from acceleration of the fluid pushed by the body in the normal direction. This has led to a set of ordinary differential equations describing the swimming movements resulting from the local muscle bending moments and the effects of the surrounding fluid. Dynamic decoupling of body shape and orientation, as in Eqns 3a–c, separates the unknown muscle torque from the body motion relative to the inertial frame, allowing for direct testing of the fluid force model by experimental motion data. The simple framework for modeling would be applicable to fish swimming in general. The static fluid model saves much time for computation when compared with CFD simulations. More importantly, our body–fluid interaction model is suitable for the analytical study of swimming mechanisms and neuronal control principles.

Most of the analytical and computational models for swimming that are currently available in the literature lack quantitative validation by experimental data. In our study, kinematic data of undulatory leech swimming were obtained via video recording and were used to determine and validate the fluid force model. The values of the drag and added-mass coefficients were chosen to minimize the error in simulated swimming motion (Fig. 8), and the resulting model behavior was reasonably close to the observation even for swim episodes not used for modeling (Fig. 9). Thus we conclude that the simple framework for modeling quantitatively captures the undulatory swimming behavior of leeches with a reasonable accuracy.

The resistive theory by Taylor (Taylor, 1952) describes hydrodynamic forces on a body in a constant flow using static functions of the relative velocity, and is considered applicable for swimming at a low Re. In contrast, the reactive theory by Lighthill (Lighthill, 1960; Lighthill, 1970; Lighthill, 1971) considers the high Re regime and assumes that the thrust for swimming comes from the reactive force of the fluid accelerated by the body motion. The latter theory has proven useful for understanding mechanisms underlying carangiform swimming of fishes with laterally compressed caudal fins and large body undulation constrained to the posterior half or even one-third of body. For anguilliform swimming of elongated animals such as eels or leeches, however, it has been an unresolved issue whether one of the resistive and reactive forces dominates. Lighthill showed that the resistive force is not negligible or is even important when compared with the reactive force in his calculation of energy loss to the wake (i.e. the swimming efficiency) under anguilliform swimming (Lighthill, 1970). A survey of aquatic animal propulsion (Lighthill, 1969; Lighthill, 1971) suggested that undulatory swimming of invertebrates were best studied using the resistive theory. A rationale is that invertebrate swimmers do not have caudal fins like carangiform swimmers, and their cross-sections may not be of the form that would enhance the added-mass effect. However, leeches use their dorso-ventral muscles to flatten the body and this would increase the added-mass effect when compared with a body with a circular cross-section. Thus, leeches, as well as eels, may adopt a combination of resistive and reactive forces for propulsion.

Our experimental data and model-based analysis suggest that the resistive force is dominant in leech swimming. The optimization of the fluid coefficients to fit the motion data by simulated behavior has lead to c_a=0, i.e. the best-fit results when the added-mass effect is removed. Our simulation results using only the resistive force model (Fig. 9) indicate that leech swimming is well described by the resistive theory. However, Fig. 8 shows that the resistive force (represented by the value of c_p) decreases as the reactive force increases, and the net effect may be distributed to the two terms in a somewhat arbitrary manner. Indeed, the error between the simulation and the data is insensitive to the added-mass coefficient c_a, and can be small (in the range 0≤c_a≤0.3) if c_p and c_t are chosen appropriately. We tried another simulation with fluid coefficients c_a=0.3, c_p=3.2 and c_t=0.9, the triple with the largest acceptable added-mass effect, and found that the fluid force distribution along the body were very similar to that shown in Fig. 10. The increasing amplitude of fluid force from the anterior to the posterior is preserved. It shows that the effects of resistive and reactive fluid forces are exchangeable to some extent. In the range of 0≤c_a≤0.3, c_a is quite small compared with the value of c_a=1 in Lighthill's elongated-body theory, and thus the contribution of the reactive force is small. For instance, when c_a=0.2, the reactive force occupies less than 20% of the total resistive force calculated for c_a=0. We conclude that the resistive force is dominant in leech swimming. Our result is consistent with McHenry's results on ascidian larvae (McHenry et al., 2003), suggesting that the acceleration reaction force “does not play a role in the dynamics of steady undulatory swimming at Reynolds number ∼ 10²”.

The pressure drag coefficient (c_p=3.0) turned out to be much larger than the value for a smooth circular cylinder in a constant flow (c_p=0.9–1.1). The difference may be partly accounted for by the geometry of the body and the roughness of body surface. The cross-section of the leech body is elliptical, closer to a flat plate than to a cylinder. The normal drag coefficient for an infinite-length flat plate is c_p=2–2.4 (Crowe et al., 2007). The flatness of the cross-section increases the normal drag coefficient. The roughness of the body surface generates additional normal drag coming from the projection of the circumferential skin friction onto the normal direction. Another factor is the dynamic nature of the fluid flow. The fluid around the body is set in motion by the anterior part of the body and may have a nontrivial effect on the hydrodynamic force. Such an effect is not captured by our model, where static fluid is assumed. The large predicted c_p value suggests that the dynamical flow tends to increase the hydrodynamic force on the body.

For carangiform swimmers that exploit reactive fluid forces, it is well known that the net thrust is determined by the motion of the tail tip (Lighthill, 1960). However, a number of observations have been made for anguilliform swimmers to support the hypothesis that the thrust is generated not only with their tail but also with more anterior body parts (Blickhan et al., 1992; Muller et al., 2001; Tytell and Lauder, 2004). Experimental particle image velocimetry (PIV) measurements of swimming eels suggested that the maximum flow velocities adjacent to the body increase almost linearly from head to tail (Muller et al., 2001) and that the wake velocity is almost in the lateral direction without a prominent downstream jet just behind the tail (Tytell and Lauder, 2004). Both of these results led to the conclusion that eels generate thrust continuously along their body. CFD simulations of anguilliform swimming provided similar results. Solutions to the two-dimensional Navier–Stokes equations supported the view that anguilliform swimmers can use most of their body length to generate vortex structures (Carling et al., 1998). The three-dimensional Navier–Stokes equations also suggested that, in addition to the tail, the middle of the body should contribute significantly to the thrust during efficient steady swimming (Kern and Koumoutsakos, 2006).

Our results on leech swimming, based on the simple fluid force model and experimental kinematic data, are consistent with the previous results on anguilliform swimming. Fig. 11 shows that the thrust is generated continuously along the body, linearly increasing in magnitude toward the tail, as previously predicted (Muller et al., 2001). The head region experiences most of the drag. The vertically pointed fluid force vectors shown in Fig. 10 predict that the wake consists of jets normal to the swimming direction, leading to a pair of vortices. This corresponds to the lateral jets observed in eel swimming (Tytell and Lauder, 2004) because the undulation occurs in a vertical plane for leeches. Moreover, the force vector has a large magnitude at the tail link, consistent with the observation that large lateral jets are generated near (but just anterior to) the tail tip (Tytell and Lauder, 2004). Thus, our simple model appears to capture the essence of the hydrodynamic forces exerted on the body during undulatory swimming. At the same time, the consistency with previous results on eels indicates that the mechanisms underlying leech swimming are similar to anguilliform swimming in spite of the lack of fins in leeches.

Our analysis supports Gray's view of the propulsion mechanisms for anguilliform swimming (Gray, 1933a) and explains continuous thrust generation along the body. Gray found that: (1) the transverse velocity of each body segment is maximal when it crosses the nominal line, and (2) the trajectory of an arbitrary point on the body in the frame of reference moving forward at the average swimming velocity is a ‘figure of eight’ curve, moving backwards at the middle point that lies on the nominal line. We observed that the fluid force on each body segment is almost normal to the body midline, and reaches the maximal magnitude when the segment is crossing the nominal line (Fig. 10). With an appropriate angle of attack set by the figure-of-eight curve, this large force has a component directed forward, resulting in thrust. This mechanism applies to most body segments, explaining continuous generation of thrust over the whole body. The hydrodynamic force on the tail segment becomes prominently large during each cycle. However, this does not lead to markedly large thrust at the tail because the large force acts as thrust when the tail is approaching the nominal line, but as drag when moving away from it. On average, the thrust generated by the tail segment is only slightly larger than those at other segments.

We deduced from Fig. 13A that actuation waves of muscle torque (or tension) travel two or three times faster than the body curvature waves. Similar observations have been made for various fishes. Direct measurements of muscle activity through electromyography (EMG) have shown that EMG signals travel faster than body curvature waves for lamprey, eel, trout, saithe, carp and scup (Wardle et al., 1995; Williams et al., 1989). For anguilliform swimmers, the speed of the EMG signal was approximately 50% faster (Williams et al., 1989), but slower than the speed of muscle torque waves we calculated for leeches. Perhaps it is because the EMG signal, corresponding, though indirectly, to the activation signal from motoneuron to the muscle, does not directly represent the muscle bending moment or tension, which would depend upon both the activation (EMG) and current length. Models with a continuum body subject to reactive hydrodynamic forces predicted, based on kinematic data from saithe, a carangiform swimmer (Videler and Hess, 1984), that the total muscle bending moment is a standing wave (or very fast traveling wave) (Hess and Videler, 1984), but its active component (after subtracting the passive visco-elastic effect) travels at a speed comparable to the EMG signal (Cheng et al., 1998). An analysis of a flexible beam model with resistive hydrodynamic forces has shown that the difference between the speeds of EMG and body curvature waves in anguilliform swimming depends on the passive visco-elasticity and body geometry (McMillen et al., 2008). These studies suggest that the speed of EMG signal is different from that of the total bending moment (or tension) (resulting from both active and passive muscle contractions as our case) and the passive visco-elasticity mediates this speed difference.

A simple analysis of the linearization of Eqn 1a without fluid, , shows that a large part of the speed difference between total bending moment (u) waves and body curvature (φ) waves can be explained by the body geometry, the characteristics of which are embedded in J(θ) and D. In particular, traveling waves of curvature φ with uniform amplitudes at a constant speed would be achieved by waves of torque input u traveling at a greater speed, if the body were floating in space with no gravity or fluid. For instance, the phase lag in φ of 360 deg from head to tail would be achieved if u has a phase lag of approximately 100 deg with a bell-shaped distribution of amplitude over the body. The muscle torque u thus predicted is very much like that shown in Fig. 13, except that the phase lag in Fig. 13 is somewhat larger. The difference may be attributed, in essence, to the hydrodynamic effect.

The mechanical power supplied by the muscle to the body is the product of u_i and the rate of change of the joint angle . In addition to the amplitudes of u_i and φ_i, the phase relationship between them is an important determinant of the amount of energy supplied over a cycle. The phase curve of is obtained by shifting that for φ (shown in Fig. 13A) upward by 90 deg because φ_i(t) is close to a sinusoid (Fig. 12B). We then see that u_i and the curvature derivative are approximately in phase (zero phase difference) in the mid-body (around joint 9). This means the instantaneous power output from the mid-body muscles is always positive over the cycle. However, u_i and are approximately 90 deg different in phase near the head and tail, implying that the integral of the muscle power in one cycle is roughly zero. In addition, Fig. 13B shows that the muscle bending moments have larger amplitude in the mid-body. These observations explain the bell-shaped muscle power distribution shown in Fig. 15.

The supplied energy is used to overcome the fluid drag and maintain a steady speed of swimming. All the energy is eventually dissipated into the fluid at the same rate as the supply during steady swimming. Fig. 15 indicates that the energy dissipation occurs almost uniformly over the body except at the tail. Thus, a large part of the body uses the energy efficiently to generate nontrivial thrust (Fig. 11) without incurring significant energy loss to the fluid; however, a substantial fraction of the energy is released into the fluid near the tail. The distribution of the hydrodynamic force along the body (Fig. 10) suggests that the energy dissipation at the tail is in the form of fast jets of fluid normal to the body surface, upwards and downwards.

Our body–fluid interaction model is based on Taylor's resistive force theory, and the fluid force always does negative work to the body. Hence the energy supplied by the mid-body is transmitted to the tail through the body, not through the fluid. The mechanisms underlying the energy transfer from mid-body to tail can be simply explained by potential energy stored in the passive stiffness of the body wall. Fig. 16 depicts a leech body that is represented by a chain of three links connected by two flexible joints. The active contraction of anterior muscle causes the middle link to rotate clockwise, thereby stretching of posterior muscle and converting the active energy supplied by the anterior muscle to passive elastic energy in the posterior muscle. Indeed, the antiphase relationship between the joint angles and muscle bending moments observed near the tail (Fig. 13A) indicates that the tail oscillation is passive. The energy transfer mechanism depicted in Fig. 16 suggests that mechanical resonance associated with the body inertia and muscle stiffness would be exploited for increased swimming efficiency.

A systematic phase shift between EMG signals and strain along the body has been observed for several species of fish, suggesting that different roles are played by the muscles at different locations of the body. The power generated by muscles at a particular location has been estimated from experiments on isolated muscles under the same activation (EMG) and strain as those observed during intact swimming. The power transmission found here for leeches is similar to the previous result for saithe, a carangiform swimmer: it was found (Altringham et al., 1993) that the power is supplied by anterior muscle and is transmitted to the tail through stiffening of the posterior muscle. However, a review article (Wardle et al., 1995) observed that anguilliform swimmers (eel and lamprey) lack tail blades and that the phase shift of EMG in the strain cycle along the body is not significant; the authors stated that the power generated by the muscle would be passed directly to the water along most of the body length, with the aid of the erect dorsal and ventral fins. Our result indicates that energy transmission mechanisms in leech swimming are different from what is described for anguilliform swimming.

There are several factors that could explain the discrepancy. Rome and coworkers studied another carangiform swimmer, scup, and reached a different conclusion from that for saithe: “most of the power for swimming came from muscle in the posterior region of the fish, and relatively little came from the anterior musculature” (Rome et al., 1993). Thus, the powering mechanism may be sensitive to body form (geometry, mass, stiffness distribution, etc.) and could be different for fishes with similar swimming styles. Another possible factor is that the EMG signal may not be exactly in phase with the actual tension developed because the tension depends on both activation and strain. In this case, the power expenditure estimated from the phase relationship between the EMG signal and strain (Wardle et al., 1995) can be erroneous. Our accompanying paper on muscle activation mechanisms underlying leech swimming (J.C., J. Tian, T.I. and W.O.F., unpublished) indicates that the motoneuron activation signal and the resulting tension have different speeds of traveling waves. In particular, the phase of motoneuron activation in the strain cycle varies little along the body whereas, as shown in Fig. 13A, the phase difference between the torque (or tension) and strain varies greatly, resulting in the energy transmission mechanism illustrated in Fig. 16.

 
• a_ni

  normal acceleration of the fluid

 
• B(θ)

  17×17 input coefficient matrix

 
• c_a

  added-mass coefficient

 
• CFD

  computational fluid dynamics

 
• c_p

  drag coefficient in the normal direction (pressure drag)

 
• CPG

  central pattern generator

 
• c_t

  drag coefficient in the tangential direction

 
• D

  18×17 constant coefficient matrix

 
• d_i

  width of the body link

 
• e

  18-dimensional vector with all its entries equal to one

 
• E

  energy returned from the fluid to the body over a cycle to generate the thrust

 
• E_i

  energy dissipated into the fluid over one cycle at the ith link (i=1,...,18)

 
• EMG

  electromyography

 
• F(θ)

  17×18 coordinate transformation matrix

 
• f_i

  fluid force vector on the ith link

 
• ,

  average thrust and drag over one cycle on the ith link

 
• f_ni

  normal fluid force on the ith link

 
• f_ti

  tangential fluid force on the ith link

 
• f_x, f_y

  net fluid forces in the horizontal and vertical directions, respectively, over the whole body

 
• g

  two-dimensional vector whose first entry is zero and second entry is (1–ρ/ρ_l)g

 
• g

  gravity constant

 
• G(θ)

  18×18 centrifugal matrix

 
• Group O

  group of episodes used for optimization (O) of fluid coefficients

 
• Group V

  group of episodes used for validation (V) of fluid coefficients

 
• h

  net fluid force vector on the center of mass of the body

 
• H(θ)

  17×18 projection of centrifugal matrix

 
• J(θ)

  18×18 moment of inertia matrix

 
• L

  length scale (body length or diameter)

 
• L₁

  straight length of one body wave

 
• l_i

  length of the body link

 
• m

  mass of the body

 
• mg

  net force resulting from gravity and buoyancy

 
• q_i

  fluid torque on the ith link resulting from the normal force on a small segment integrated over the link

 
• Re

  Reynolds number

 
• RMS

  root mean square

 
• T

  cycle period

 
• t₀, t₁, t₂

  time instants when the crest of the body wave appears at the head tip, tail tip and the crest reappears at the head tip, respectively

 
• T₁

  duration of the swimming episode

 
• u

  17-dimensional vector whose ith entry is u_i

 
• U

  swim speed

 
• u_i

  muscle bending moment applied at the ith joint

 
• V

  body wave speed

 
• v_i

  velocity of the center of mass of the ith link

 
• v_ni

  normal component of v_i

 
• v_ti

  tangential component of v_i

 
• v_x, v_y

  horizontal and vertical components, respectively, of the velocity of the center of mass of the body

 
• w

  two-dimensional vector whose entries are x and y

 
• W

  total work done by the body to the fluid over a cycle

 
• 

  simulated velocity of the center of mass of the body

 
• W_i

  work done to the body by muscle over one cycle at the ith joint (i=1,...17)

 
• 

  measured video data of the velocity of the center of mass of the body

 
• x, y

  position of the center of mass of the body

 
• , , ,

  horizontal and vertical components of and , respectively

 
• Δt

  time it takes for the body wave to propagate from joint 5 to joint 13

 
• α_i

  amplitude of φ_i

 
• α_ni

  normal acceleration of the center of mass of the ith link

 
• β_i

  phase of φ_i

 
• γ_i

  bias of φ_i

 
• μ

  viscosity of fluid

 
• θ

  18-dimensional vector whose ith entry is θ_i

 
• θ_i

  the angle between the ith link and the +x-axis, counterclockwise is positive

 
• θ_o

  body orientation, defined as

 
• 

  velocity of θ

 
• 

  18-dimensional vector whose ith entry is

 
• 

  angular velocity of the ith link

 
• 

  velocity of θ_o

 
• 

  acceleration of θ

 
• ρ

  density of fluid

 
• ρ_l

  density of leech body

 
• τ

  18-dimensional vector whose ith entry is τ_i

 
• τ_i

  total fluid torque on the ith link

 
• φ

  17-dimensional vector whose ith entry is φ_i

 
• φ_i

  angle at the ith joint defined as θ_i–θ_i+1 (i=1,...,17)

 
• 

  velocity of φ

 
• 

  velocity of φ_i

 
• 

  acceleration of φ

 
• ψ

  average angular momentum

 
• ω

  angular frequency