JEB desktop wallpaper calendar 2016

Journal of Experimental Biology partnership with Dryad

Mechanisms underlying rhythmic locomotion: body–fluid interaction in undulatory swimming
J. Chen, W. O. Friesen, T. Iwasaki


Swimming of fish and other animals results from interactions of rhythmic body movements with the surrounding fluid. This paper develops a model for the body–fluid interaction in undulatory swimming of leeches, where the body is represented by a chain of rigid links and the hydrodynamic force model is based on resistive and reactive force theories. The drag and added-mass coefficients for the fluid force model were determined from experimental data of kinematic variables during intact swimming, measured through video recording and image processing. Parameter optimizations to minimize errors in simulated model behaviors revealed that the resistive force is dominant, and a simple static function of relative velocity captures the essence of hydrodynamic forces acting on the body. The model thus developed, together with the experimental kinematic data, allows us to investigate temporal and spatial (along the body) distributions of muscle actuation, body curvature, hydrodynamic thrust and drag, muscle power supply and energy dissipation into the fluid. We have found that: (1) thrust is generated continuously along the body with increasing magnitude toward the tail, (2) drag is nearly constant along the body, (3) muscle actuation waves travel two or three times faster than the body curvature waves and (4) energy for swimming is supplied primarily by the mid-body muscles, transmitted through the body in the form of elastic energy, and dissipated into the water near the tail.


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/L1), 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, L1 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.


Collection of swimming motion data

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

Fig. 1.

Diagram of the leech swimming system. Commanded by the central pattern generator (CPG), motoneurons send impulse trains to the muscles and cause them to contract, where the tension is determined by both motoneuron activation and muscle strain. The body, subjected to the muscle actuation and fluid forces, forms a quasi-sine shape, sending traveling waves rearward for thrust generation. Sensory feedback of the muscle tension through stretch receptors modifies the cycle period and phases of the CPG to adjust the undulatory swimming style appropriately.

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.

Measurements of body dimension, mass and density

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.

Fig. 2.

Video-captured data of leech swimming. (A) Body contour of a swimming leech and discretization of the body midline into 18 links. The right end is the head and the joint angles are denoted by φ1,...,φ17 from tail to head. The joint angle is positive when the body is in a configuration such that ventral muscle contracts and dorsal muscle extends. (B) The time courses of joint angles φ5 and φ13 in swimming cycles. Joint angles change nearly sinusoidally [compare data (circles) with fitted sinusoids (solid lines)]. The body wave propagates from right to left in A, and φ5 lags φ13 by 177 deg in B, which is calculated by (Δt/T)×360. (C) Phases of all joint angles. The abscissa is the body joint index. The joint angle at the tail end (φ1) lags that at the head end (φ17) by 360 deg, which indicates that one full body wave occurs between the two end joints.

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.

Body–fluid interaction model

The leech body is modeled as a chain of 18 rigid links connected by 17 frictionless rotational joints, constrained to move within a vertical plane. We take the x- and y-axes as horizontal and vertical, respectively. The counterclockwise angular displacement of the ith link from the x-axis is denoted by θi, and the 18-dimensional vector obtained by stacking θi in a column is θ∈R18. The coordinates of the center of mass for the whole body are denoted by a two-dimensional vector wR2 with components (x, y). The longitudinal muscle generates bending moment uiR at each joint, which is defined as negative when it tends to make the body bend into a V-shape. In addition to the muscle tension caused by motoneuron activation, the passive visco-elastic tension caused by the length change in the body wall is also contained in ui. Each link is subject to force fiR2 and torque τiR due to the fluid; τ∈R18 is defined by stacking τi in a column vector. The equations of motion are given by (Saito et al., 2002): Embedded Image (1a) Embedded Image (1b) where J(θ),G(θ)∈R18×18 are the matrices for the moment of inertia and the centrifugal term, respectively, DR18×17 is a constant coefficient matrix, uR17 is the vector obtained by stacking ui in a column, Embedded Image is the vector whose ith entry is Embedded Image, mR is the total mass of the body and hR2 is the net fluid force on the center of mass. The term mg is the net force resulting from gravity and buoyancy; gR2 is the vector whose first entry is zero and the second entry is (1–ρ/ρl)g, where ρ and ρl are the densities of the fluid and leech body, respectively, and g is the gravity constant.

Fig. 3.

The top view of a medium-sized leech during swimming, from which the body width variation along the body was measured. The right end is the head.

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.

The model for the hydrodynamic force on the link is based on the resistive theory of Taylor (Taylor, 1952) and the reactive theory of Lighthill (Lighthill, 1960). For each link, the velocity of its center of mass vi is split into the normal and tangential components vni and vti, respectively, as shown in Fig. 4. The hydrodynamic force in the normal and tangential directions on the ith link (fni and fti, respectively) are then modeled as: Embedded Image (2a) Embedded Image (2b) where ρ and μ are the density and viscosity of the fluid; di and li are the width and length of the body link, respectively; cp and ct are drag coefficients in the normal (pressure drag) and tangential directions, respectively; ca is the added-mass coefficient; and ani is the normal acceleration of the fluid, which is assumed to be pushed by the link in the normal direction only and to slide in the tangential direction (Embedded Image, where αni is the normal acceleration of the center of mass of the link). The function sgn extracts the sign of a real number. The tangential force fti and the first term in the normal force fni are static functions of the relative velocity between the body link and the presumably static fluid, representing the resistive forces (‘static functions’ mean that the fluid force at any time instant is determined by the current velocity only, and is independent of the velocity history). These force components are estimated from the drag on an infinite-length circular cylinder placed obliquely in a constant flow (Taylor, 1952). The second term in Eqn 2b is a static function of the acceleration of the body, representing the added-mass effect due to the fluid acceleration resulting from the body motion (Lighthill, 1960). The fluid torque qi results from the normal force on a small segment integrated over the link. The first and second entries of the fluid force vector h in Eqn 1b are obtained by summing the x and y components of fi over the body, respectively. The total fluid torque τi on the ith link (the ith entry of vector τ in Eqn 1a) comprises the direct effect of the local qi and the indirect effect of the fluid forces acting on the other links through the mechanical linkage.

Fig. 4.

Hydrodynamic force on a link. The normal and tangential components of the link velocity are used to model the force components. See ‘List of symbols and abbreviations’ for parameter definitions.

Determination of fluid drag coefficients and prediction of muscle bending moment and fluid force

We determined ca, cp and ct 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.

We use the joint angle φi (=θiθi+1, where i=1,...,17) to describe the body shape and the average angular momentum ψ [Embedded Image] to represent the change in the body orientation, where eR18 is the vector with all its entries equal to one. By the coordinate transformation to replace the variable Embedded Image by Embedded Image and ψ∈R, Eqn 1a is decoupled into the body shape equation and the body orientation equation, respectively (Saito et al., 2002). As a result, we have: Embedded Image (3a) Embedded Image (3b) Embedded Image (3c) where H(θ), B(θ) and F(θ) are appropriately defined coefficient matrices depending on θ. The unknown muscle bending moment u is now decoupled from the equations of inertial motion (rotation ψ and translation w).

The inertial motion can be predicted by numerically simulating the equations for ψ and w using the observed swimming data for the shape change φ as the input, where the fluid torque τ and force h depend on φ and Embedded Image. The polynomial fitting of φ is used to describe the body shape change in simulations because the body undulation from cycle to cycle is not exactly periodic. The predicted inertial motion depends on the fluid coefficients in Eqns 2a and 2b. The values of cp, ct and ca were determined by minimizing the root mean square (RMS) of the difference between the simulated velocity of the center of gravity Embedded Image and the measured video data Embedded Image, i.e.: Embedded Image (4) where Embedded Image and Embedded Image are the horizontal and vertical components of Embedded Image and Embedded Image, respectively, and T1 is the time length of the video data (roughly two cycles of undulation). The summation is taken over four swimming episodes of four different leeches whose sizes range from small to large. In the numerical optimization, the integral was approximated by gridding the time axis. The body orientation θo: Embedded Image (5) could have been used as part of the error function being minimized, but we chose to use it for testing the fluid coefficients resulting from this optimization. Four additional swimming episodes were used to validate the model.

In addition to quantitative estimates for the fluid coefficients cp, ct and ca, 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.


Basic data for leech body and swimming motion

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.

Table 1.

Data summary for leech body and swimming episodes

Fig. 5.

Variation of body width, thickness and mass along the body of five leeches. The leech body is divided into 18 segments of equal length, and the abscissa is the segment index. The measurements for body width and thickness are taken at the middle of each segment. The body width tapers significantly toward the head. Body thickness is reduced at both ends. The mass is concentrated in the thick and wide middle portion of the body.

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 (t1t0)/(t2t0), where a crest appears at the head tip at time t0, travels down the body and arrives at the tail tip at t1, and the next crest appears at the head tip at t2. 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 ReUL/μ, 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.

Fig. 6.

Snapshot sequence of a leech body during swimming. One complete cycle of undulation is shown. The body contours are extracted from video images recorded by a high-resolution camera at 60 frames s–1. The snapshots are plotted every second frame. The leech swims from left to right. The body presents about one quasi-sine wave and the cycle period is 0.38 s.

Optimal fluid coefficients

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 ct and ca, cp was optimized to minimize the RMS error in the simulated swim velocity. Fig. 8A shows the contours of optimal cp for different values of ct and ca. The optimal cp decreases as ca increases for a fixed ct, and increases as ct increases for a fixed ca. These observations make sense if we note that the thrust is generated by the normal hydrodynamic force through the terms associated with cp and ca, and is balanced with the longitudinal drag due to ct. 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 ca approaches zero, which suggests that the added-mass effect is small in undulatory leech swimming. The lowest RMS error happens at ca=0, ct=0.6, and the corresponding optimal cp value is 3.0. The optimal value of cp 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<105. The RMS error is not very sensitive to the fluid coefficients near the region of ca around zero.

Fig. 7.

Oscillation profile of the joint angles φi (body curvature) along the body, measured during steady swimming in a fixed, roughly horizontal direction. The oscillation of joint angles can be well fitted by sinusoids (see Fig. 2), i.e. φi(t)≈αisin(ωt+βi)+γi, where αi, γi and βi are defined as the amplitude, bias and phase of φi, respectively. (A) The oscillation amplitude of joint angles. The colored lines represent five leeches (e.g. 2-1 is leech 2, episode 1). The abscissa is the body joint index. The oscillation amplitude increases toward the tail. (B) The bias of joint angles. The biases are almost zero all along the body except near the head. Hence, the nominal body shape, around which undulations occur, is nearly straight with the head curling slightly upward. (C) The phase of joint angles. The phase lag from head to tail is roughly 360 deg for leeches 1, 2 and 5, expressing one full wave along the body. The phase lag is a little less than 360 deg for leeches 3 and 4, and the body exhibits approximately 80 to 90% of a full wave between the tail joint φ1 and the head joint φ17.

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 vx (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 vy 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 ca is zero for all eight swimming episodes, and five out of eight swimming episodes (four out of five leeches) have cp clustered in the range of 2.2–3.2 and ct in the range of 0.4–1, which are around the optimal values (cp=3.0, ct=0.6). Leech 2 was the outlier, with larger values of optimal coefficients; cp in the range of 4.8–7.6 and ct in the range of 0.9–1.8 for three swimming episodes.

Fig. 8.

Contours of optimal pressure drag coefficient cp (A) and RMS error (cm s–1) in the simulated velocity (B) in the cact plane, where ca is the added-mass coefficient and ct is the tangential drag coefficient.

Hydrodynamic force on the body

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 fx and fy 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 instant1. The nominal line has a positive pitch angle (roughly equal to θo), and this may be related to the fact that the vertical force fy 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 fx varies within the cycle; the leech body accelerates (fx >0) when the tail is approaching the nominal line, and decelerates (fx<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.

We define the thrust and drag according to Gray (Gray, 1933b) and Taylor (Taylor, 1952), who studied propulsive mechanisms of anguilliform swimming. The body waves travel backward at a speed faster than the forward swimming, and this makes the body push the fluid at some angle of attack and generates the normal and longitudinal components of the hydrodynamic force. The normal component is responsible for the thrust and the longitudinal component always drags the body backward. Hence, for each link, the thrust is calculated by projecting the normal fluid force onto the x-axis, whereas the drag is calculated by projecting the tangential fluid force onto the x-axis. The average thrust and drag over one cycle are thus given by: Embedded Image (6a) Embedded Image (6b) for i=1,...18 (see Fig. 4). The thrust and drag calculated in this way were normalized by the total thrust of the whole body and are plotted in Fig. 11, where the negative of Embedded Image is shown to indicate that the drag is in the direction of the negative x-axis. Four swimming episodes are indicated by curves of different colors. The abscissa is the body location; 0 is the tail tip and 1 is the head tip. The net thrust (i.e. integral over the body) is roughly equal to the net drag because they balance during steady swimming at a constant speed. Fig. 11 confirms that the tangential force is always negative (drag), whereas the normal force is mostly positive (thrust), except near the head. The thrust increases roughly linearly towards the tail, but the drag is more nearly constant.

Body actuation by muscle bending moment

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

Fig. 9.

Simulations of four swimming episodes (blue) and their comparisons with the video data (red). (A) Velocity of the center of mass; vx and vy are the horizontal and vertical components of Graphic, respectively. (B) Angular velocity of the body orientation Graphic. The observed time course of the body shape is imposed on the model during each simulation. The fluid model is Taylor's resistive force (i.e. ca=0) with the optimal drag coefficients cp=3 and ct=0.6. The plots for episodes 1 and 4 show a dynamical curve fit of the data by the model, whereas those for episodes 3-2 and 5 are model predictions compared with the data. This is a full comparison in the sense that there are no other hidden variables to compare; exact matching of vx, vy and Graphic implies that the swimming behaviors are identical.

Energy transmission from muscle to fluid

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, Embedded Image. 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.

The work done on the body by the muscle over a cycle was calculated at each joint by: Embedded Image (7) for i=1,...,17. The energy dissipated from the body to the fluid in a cycle was calculated at each link by: Embedded Image (8) for i=1,...,18 (see Fig. 4). These quantities were normalized by the total energy supplied by the muscle in one cycle, and their variations along the body are shown in Fig. 15 for four swimming episodes, where the negative of Ei is plotted to indicate that the energy is lost. The sum of all Wi is the total work done by the muscle, and is equal to the sum of all Ei, the total energy dissipation into the fluid. The energy supplied by the muscle is bell-shaped – large in the mid-body and smaller toward both ends. The energy dissipated into the fluid is small and almost constant along the body except at the tail, where the velocity of the tail link is large because this link is unconstrained at its posterior terminus. Clearly, the muscle power is supplied around the mid-body and is released to the fluid at the tail.


Simple modeling framework is broadly applicable to fish swimming

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.

Fig. 10.

Fluid force on the body in one swim cycle. The leech swims to the right. The blue arrows are the fluid force vectors, each of which comprises the normal and tangential components. The parameters fx and fy indicate the net fluid forces in the horizontal and vertical directions, respectively; the former oscillates around zero, but the latter is always positive to account for gravity. A large fluid force is observed at each body location when it crosses the nominal line. The tail tends to generate large fluid forces most of the time. The anterior body experiences smaller fluid force than the posterior body. The whole body has a positive pitch angle, which may be related to the gravity effect.

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.

Resistive force dominates

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.

Fig. 11.

Hydrodynamic thrust and drag generated along the body over one swim cycle, normalized by the total thrust. The abscissa indicates the location along the body; 0 is the tail tip and 1 is the head tip. The colored curves (as in Fig. 13B) represent four swimming episodes (1, 2-1, 3-1 and 4). Thrust and drag are defined as projections of the normal and tangential forces onto the x-axis, averaged over one cycle. Negative drag indicates that the drag is in the direction opposite to the swim direction (the leeches swim toward the positive x-axis). The length of the head and tail links are twice as long as the middle links, and hence the force on each of these links is divided in two and is plotted as two data points. The thrust increases roughly linearly toward the tail; the drag is nearly constant.

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 ca=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 cp) 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 ca, and can be small (in the range 0≤ca≤0.3) if cp and ct are chosen appropriately. We tried another simulation with fluid coefficients ca=0.3, cp=3.2 and ct=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≤ca≤0.3, ca is quite small compared with the value of ca=1 in Lighthill's elongated-body theory, and thus the contribution of the reactive force is small. For instance, when ca=0.2, the reactive force occupies less than 20% of the total resistive force calculated for ca=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 ∼ 102”.

Fig. 12.

Time courses of body actuation waves and the resulting body curvature waves. (A) The muscle bending moments during swimming, predicted by the body–fluid interaction model with experimental kinematic data. (B) The joint angles (body curvature) measured through video recording. The corresponding swimming episode is 3-1.

The pressure drag coefficient (cp=3.0) turned out to be much larger than the value for a smooth circular cylinder in a constant flow (cp=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 cp=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 cp value suggests that the dynamical flow tends to increase the hydrodynamic force on the body.

Thrust generated continuously along 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).

Fig. 13.

Phase and amplitude variations along the body for the muscle bending moment at 17 joints (solid lines). For comparison, the phase of joint angles (body curvature) is also shown by dashed lines. The colored lines represent four swimming episodes (1, 2-1, 3-1 and 4). The abscissa is the joint index indicating the location along the body. The amplitude and phase were calculated from sinusoidal fittings of the time course. (A) Phase variation. The phase lag from head to tail is proportional to the time it takes for the waves to travel from the head joint to the tail joint. The speed of the muscle bending moment wave is two or three times faster than the body curvature wave. (B) Amplitude of muscle bending moment at 17 joints. Muscle actuation is bell-shaped: large in the mid-body and small near the head and tail.

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.

Muscle actuation waves travel faster than body curvature waves

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.

Fig. 14.

Time course of the instantaneous power Graphic by the muscle at several joints along the body. The power supply is almost always positive in the anterior part of the body, but takes both positive and negative values in the posterior. Negative power implies energy storage in passive elasticity. The head has almost no energy supply or storage, but the tail exhibits energy exchange with an average of ∼0. The corresponding swimming episode is 3-1.

Fig. 15.

Mechanical work done to the body by the muscle (Wi, for i=1,...,17) and by the fluid (–Ei, for i=1,...,18) in one cycle along the body, normalized by the total muscle work (sum of all Wi). Because the head and tail links are twice as long as the middle links, each of E1 and E18 was divided by two and is plotted as two data points. The abscissa is the body location; 0 is the tail tip and 1 is the head tip. The colored curves (as in Fig. 13B) represent four swimming episodes (1, 2-1, 3-1 and 4). The power is supplied to the mid-body and is dissipated to the fluid at the tail.

A simple analysis of the linearization of Eqn 1a without fluid, Embedded Image, 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.

Energy supplied by muscle at mid-body, dissipated into fluid at tail

The mechanical power supplied by the muscle to the body is the product of ui and the rate of change of the joint angle Embedded Image. In addition to the amplitudes of ui and φi, the phase relationship between them is an important determinant of the amount of energy supplied over a cycle. The phase curve of Embedded Image 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 ui and the curvature derivative Embedded Image 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, ui and Embedded Image 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.

Fig. 16.

Schematic for energy transmission along the body. The dashed line indicates the initial configuration of the body links. Muscle power is supplied to the anterior body through active contraction of the dorsal muscle. The contraction causes the middle link to rotate clockwise and induces stretching of the posterior muscle, storing the supplied energy in passive muscle stiffness.

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.


  • 1 The nominal line would coincide with the so-called ‘axis of forward movement’ (see Gray, 1933a) if there were no gravity and the nominal line became horizontal.

  • This work is supported by the National Science Foundation under grant no. 0654070; the Office of Naval Research, under MURI grant no. N00014-08-1-0642; and NIH/NINDS 1 R01 NS46057-01 as part of the 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.


normal acceleration of the fluid
17×17 input coefficient matrix
added-mass coefficient
computational fluid dynamics
drag coefficient in the normal direction (pressure drag)
central pattern generator
drag coefficient in the tangential direction
18×17 constant coefficient matrix
width of the body link
18-dimensional vector with all its entries equal to one
energy returned from the fluid to the body over a cycle to generate the thrust
energy dissipated into the fluid over one cycle at the ith link (i=1,...,18)
17×18 coordinate transformation matrix
fluid force vector on the ith link
Embedded Image, Embedded Image
average thrust and drag over one cycle on the ith link
normal fluid force on the ith link
tangential fluid force on the ith link
fx, fy
net fluid forces in the horizontal and vertical directions, respectively, over the whole body
two-dimensional vector whose first entry is zero and second entry is (1–ρ/ρl)g
gravity constant
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
net fluid force vector on the center of mass of the body
17×18 projection of centrifugal matrix
18×18 moment of inertia matrix
length scale (body length or diameter)
straight length of one body wave
length of the body link
mass of the body
net force resulting from gravity and buoyancy
fluid torque on the ith link resulting from the normal force on a small segment integrated over the link
Reynolds number
root mean square
cycle period
t0, t1, t2
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
duration of the swimming episode
17-dimensional vector whose ith entry is ui
swim speed
muscle bending moment applied at the ith joint
body wave speed
velocity of the center of mass of the ith link
normal component of vi
tangential component of vi
vx, vy
horizontal and vertical components, respectively, of the velocity of the center of mass of the body
two-dimensional vector whose entries are x and y
total work done by the body to the fluid over a cycle
Embedded Image
simulated velocity of the center of mass of the body
work done to the body by muscle over one cycle at the ith joint (i=1,...17)
Embedded Image
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
Embedded Image, Embedded Image, Embedded Image, Embedded Image
horizontal and vertical components of Embedded Image and Embedded Image, respectively
time it takes for the body wave to propagate from joint 5 to joint 13
amplitude of φi
normal acceleration of the center of mass of the ith link
phase of φi
bias of φi
viscosity of fluid
18-dimensional vector whose ith entry is θi
the angle between the ith link and the +x-axis, counterclockwise is positive
body orientation, defined as Embedded Image
Embedded Image
velocity of θ
Embedded Image
18-dimensional vector whose ith entry is Embedded Image
Embedded Image
angular velocity of the ith link
Embedded Image
velocity of θo
Embedded Image
acceleration of θ
density of fluid
density of leech body
18-dimensional vector whose ith entry is τi
total fluid torque on the ith link
17-dimensional vector whose ith entry is φi
angle at the ith joint defined as θi–θi+1 (i=1,...,17)
Embedded Image
velocity of φ
Embedded Image
velocity of φi
Embedded Image
acceleration of φ
average angular momentum
angular frequency


View Abstract