Sink and swim: kinematic evidence for lifting-body mechanisms in negatively buoyant electric rays Narcine brasiliensis

Swimming in negatively buoyant fish is like flight in birds: gravity will bring you down. Even though the physical properties of air and water differ dramatically, soaring birds, gliding mammals (Bishop, 2006; Bishop, 2007), gliding underwater vehicles (Rudnick et al., 2004) and gliding negatively buoyant fishes use the same mechanism of controlled falling: potential gravitational energy is converted into the kinetic energy of descent and forward motion. To maximize the forward component of the glide, the body's lift-to-drag ratio must be maximized. In negatively buoyant fishes, lift is generated by paired fins or by the body itself (Nowroozi et al., 2009; Wilga and Lauder, 1999; Wilga and Lauder, 2000). Because body-generated lift will be maximized by wing-like shapes, we hypothesized that dorso-ventrally flattened negatively buoyant fishes swimming and gliding in the water column not only use their body to generate lift, but also actively adjust the body's orientation, shape and motion to control locomotor performance.

Dorso-ventrally flattened negatively buoyant fishes are exemplified by batoids, i.e. the skates and rays (Superorder Batoidii; Class Elasmobranchii). Their bodies, compared with those of the closely related sharks, are dorso-ventrally flattened, with pectoral fins that have been enlarged and integrated into an expanded neurocranium, forming the large body disc characteristic of the group (Rosenberger and Westneat, 1999). With this evolutionarily derived morphology, most batoids swim by undulating or oscillating the body disc, a modality known as median–paired-fin swimming (Webb, 1984). However, electric rays (Order Torpediniformes), even though they possess the body disc, use only the body–caudal–fin (BCF) swimming mode, undulating the caudal fin and the portion of the body posterior to the disc (Roberts, 1969; Rosenberger, 2001). Because previous research on batoid locomotion has focused primarily on species using the median–paired-fin mode (Rosenberger and Westneat, 1999), relatively little is known about the swimming of batoid species using the BCF mode.

For BCF-swimming electric rays, we propose the following functional model (Fig. 1): (1) the body disc functions as an underwater wing, a hydrofoil that generates lift during both BCF swimming and gliding; (2) the electric ray's tail, the portion of the body from the pelvic fins to the tip of the caudal fin, functions as an undulating propeller, generating forward thrust during BCF swimming. If this model is correct, then the following predictions will hold true: during horizontal BCF swimming, any changes in the pitch of the body disc will correlate with changes in swimming speed; and during gliding, glide angle and speed on the glide path will correlate with changes in the pitch angle of the body's disc and the feathering of the pectoral fins. We test these functional predictions by measuring the correlation between the motions of the tail, pelvic fins and body disc and changes in locomotor performance – measured by horizontal speed, Strouhal number, speed on the glide path and glide angle – in neonates of the lesser electric ray, Narcine brasiliensis.

Ten neonate Narcine brasiliensis (Olfers 1831) were born in captivity at Vassar College, Poughkeepsie, NY, USA (Table 1). All rays were maintained in a 10 gallon aquarium for approximately 1 month, at which point half the cohort was transferred to another 10 gallon tank. Water chemistry and animal health were monitored daily by trained staff. The rays were fed cubes of bloodworms three times a week. Five days after the first set of trials, rays 1 through 4, which were housed together in one of the two tanks, died suddenly, when the tank's filtration system failed. The remaining rays survived past the end of the trials and died individually, from a variety of causes. On the days that the rays swam for kinematic trials, they were fed in the evenings after videotaping. Animals were cared for under IACUC protocol no. 08-10B in the vivarium at Vassar College.

Video data were collected in one session from rays 1 to 4 and in five sessions from rays 5 to 10 (Table 2). Filming trials began during the fourth week of life and continued for 5 weeks thereafter, during which the animals remained in the neonate stage. Filming was conducted with a video camera (Sony HDR HC-5 Handycam) set at a 30 Hz framing rate and 1080 horizontal lines of resolution. The filming rate was deemed appropriate for two reasons: (1) pilot studies showed tailbeat frequencies always well below 3.0 Hz, giving more than 10 images per cycle; and (2) variables were means calculated over several tailbeats (see below). The camera was placed 3.35 m from the tank. The swimming tank was 25 cm wide and 53 cm in length, and filled with seawater to a depth of 15 cm; given the small size (see Table 1) and slow speed (see Results) of the rays, this tank provided room, on average, for 9 s and nine to 18 tailbeat cycles of horizontal swimming. A mirror, angled at 45 deg relative to the horizontal, was placed under the tank to allow both lateral and ventral views of a swimming ray to be captured in a single video frame. The overall field of view, which was scaled using a calibrated grid that appeared in both lateral and ventral views, was 60×80 cm. Rays were identified, weighed and photographed prior to each swimming session. Each ray swam for 5 min on camera and was not provoked to swim.

To control for the possibility that the bodies of the electric rays might glide passively (without neuromuscular control), we allowed post mortem neonates to sink in a tank filled with water to a depth of 30 cm. Prior to release, the post mortem neonates were submerged and cleared of any pockets of air in their spiracles or pharynx. We then released each neonate three times and observed the direction of their movement (up or down in the water column) and the trajectory of their movement [glide angle, θ (deg)], relative to the horizontal.

The trials on video were first screened, and three different behaviors were cataloged: (1) in BCF mode, swimming up in the water column from the bottom (upward BCF swimming); (2) in BCF mode, swimming horizontally (horizontal BCF swimming); and (3) gliding, sinking at an angle relative to the horizontal, with only intermittent undulations of the tail (gliding). Given the behavioral variability in upward BCF swimming, we chose to analyze only the horizontal BCF swimming and gliding. Within each of theses two swimming modes, video segments were further screened to make sure that the ray was clearly visible in both lateral and ventral aspects.

Horizontal BCF swimming and gliding data were analyzed from 10 rays (Table 2). The motion of the ray was manually digitized using Logger Pro 3.6.0 (Vernier Software and Technology, Beaverton, OR, USA). From the lateral view, we digitized seven body points (Fig. 2), three on the body disc, one on the pelvic fin and three on the tail. From the ventral view, we digitized six points (Fig. 2), four along the disc midline and two at the juncture of tail and pelvic fin.

Four performance variables were calculated (Fig. 3). For horizontal BCF swimming, we calculated the horizontal speed, U_h (m s^–1), as the mean, over the length of each trial, of the finite interframe differences of the displacement of the mid-disc body point, m₂, from the ventral view (mean U_h from ventral and lateral views did not differ statistically). The mid-disc body point was chosen to represent speed because, compared with the other points on the disc, m₂ yields a more conservative estimate of U_h because it varies less than m₁ or m₃. The U_h over all trials was significantly lower for m₂ than for the two other mid-points (F_2,270=24.16, P<0.0001). The mean variance of U_h over all trials, as measured by standard deviation of speed over the course of each trial, was also lower at m₂ (F_2,268=65.71, P<0.0001). During gliding, the lowest variance also occurred at m₂, as well as m₃ (F_2,39=7.31, P=0.002).

For gliding, we calculated speed along the glide path (Fig. 1C), U_p (m s^–1), as the mean, over the length of the trial, of the finite interframe differences of m₂ from the lateral view. For gliding, we also calculated the mean glide angle by trigonometry using U_h and U_p (Fig. 1B,C).

We considered a number of variables as possible candidates to correlate with changes in swimming performance (Table 3). These candidate control variables were of two types: (1) BCF undulatory thrust power and (2) lift in a static glider. To identify candidate control variables for BCF undulatory thrust power, we used Lighthill's (Lighthill, 1971) large-amplitude elongated-body theory, which has been shown, by careful experiments in swimming eels, to serve as a reasonable first approximation for calculating force and power (Tytell and Lauder, 2004). Because thrust power, which is dependent on the rate at which the body is producing mechanical work, is proportional to the square of f and the square of A (Wu, 1977), these variables were chosen. In addition, the `propulsive wavelength' is also a variable in the thrust-power calculations (Lighthill, 1971); however, because the propulsive wavelength varies along the body (Long and Nipper, 1996), we measured a more accurate local curvature, the maximum lateral flexion of the anterior and posterior portions of the tail, Ψ_A (deg) and Ψ_P (deg), respectively. Following the method of Jayne and Lauder (Jayne and Lauder, 1995), lateral flexion was measured as the deviation of three consecutive points along the body's midline from a straight line. For rays we calculated Ψ_A and Ψ_P as the difference between the angles formed by points m₃, m₄, m₅ and m₄, m₅, m₆, respectively, where a flexion of 0 deg indicates that segment is straight (Fig. 2 and Table 3).

Lifting-body theory (Harris, 1936; Vogel, 1994) inspired seven more candidate control variables (Fig. 3 and Table 3). Crucial to the lift-and-drag ratio of rigid and fixed airfoils is the angle of attack, α_A (deg), i.e. the orientation of the wing relative to the oncoming flow (Vogel, 1994). For our rays swimming horizontally, α_A is equivalent to the disc's pitch angle, α_D (deg), which is the angle relative to the horizontal. This equality is violated during gliding, when a non-zero θ alters the direction of the oncoming flow such that α_A equals the difference between θ and α_D. Thus, because α_A has an obvious mathematical relationship with θ, and is hence not independent of our performance variable θ, we chose to use only α_D as a candidate control variable. α_D was calculated as the angle from the line segment between points d₁ and d₃ and the horizontal.

The surface area of the pelvic fins of N. brasiliensis, as a percentage of the body disc, is larger than that of other ray species (Macesic and Kajiura, 2010). In addition, because wind-tunnel experiments on a model of a shark showed the importance of pelvic fins in controlling flow over the body (Harris, 1938), we measured the elevation, e (cm), of the pelvic fins as the distance between points d₃ and p (Fig. 2). Using the same experimental approach, Harris (Harris, 1936) also showed the hydrodynamic importance of the orientation and posture of the body; thus we measured (Fig. 2) the tail's pitch angle, α_T (deg), as the angle between the line segments between points d₁, t₁ and t₁, t₂; the camber of the disc, κ (deg), as the angle between points d₁, d₂, d₃; and the amplitude of the disc's yaw, ω (deg), as the angle between the midline and the axis of progression. The two other candidate control variables, the coefficients of lift and drag, C_L and C_D, respectively, are described in the next section.

Lift and drag properties of gliding animals such as rays can be estimated using the hydrodynamic resultant that balances gravity (Vogel, 1994). This steady-state method assumes that: (1) the glider is moving at constant speed along the glide path, (2) the glide path has a constant glide angle, the (3) body is not self-propelling and (4) the body is not changing shape. We only analyzed the section of the glide when the rays had reached top speed on the glide path, U_max (m s^–2), and were not beating their tails. Because of its large size and flattened shape, relative to the small size of the tail, we also assumed that the body disc was the primary lifting body during gliding. This steady-state method provides a first approximation of the lift and drag coefficients.

Our primary analytical goal was to test the predictions that the candidate control variables were significantly correlated with the performance variables associated with horizontal BCF swimming and gliding (Fig. 3). We used least-squares linear regression, with performance variables treated as dependent responses to the candidate control variables, which are treated as independent variables. For each of the four performance variables – U_h and St for horizontal BCF swimming and U_p and θ for gliding – we performed two different tests, univariate-regression and multivariate-regression screening. Even though multivariate regression offers a more appropriate test of the predictions, univariate regressions are included because many researchers continue to use them to build simple mathematical relationships.

In the first test, univariate-regression screening, we separately regressed each performance variable onto a single candidate control variable. With this approach, we run the risk of overestimating the number of significant control variables if the control variables are correlated with each other (Zar, 1999). To understand the extent of this risk in our case, we performed all possible pairwise correlations among the independent variables.

In the second test, multivariate-regression screening, we separately regressed each performance variable onto the entire pool of candidate control variables using a two-step process. Step one involved a stepwise linear regression, mixed direction, with P=0.25 for entry to and withdrawal from the iterative process. Stepwise regression examines various combinations of candidate variables to find the one combination that yields the best model, as measured by the coefficient of determination, r². In step two, the candidates selected by the stepwise method were used as independent factors in a linear multiple regression model. With this approach, we run the risk of overestimating the number of control variables if variables selected by the stepwise approach fail to be significant in the multiple regression. To mitigate this risk, we report the full multivariate models but focus only on significant variables. We also performed all possible pairwise correlations among the independent variables. All statistical tests were run in JMP 8.0 with a significance level of 0.05.

Data are presented as means ± s.e.m.

When two dependent variables are mathematically related to a common independent variable, the parameter, the two so-called `simultaneous' equations can be combined into a single equation by substitution. This yields a new equation in which the independent variable is hidden; it is `parameterized' within the new function. Recognizing and simplifying families of parametric equations in functional systems described by simultaneous univariate equations is one way to construct multidimensional functional models. Another way, noted in the previous section, is to create multivariate regression models. The difference between the approaches, as we use them here, is that our parametric models relate multiple responses to a common independent variable whereas our multivariate regression models relate a single response to multiple independent variables. In this way, the two approaches are complementary.

All values of the ρ_b, except one (Table 1), were above the ρ_w value of 1024 kg m^–3 (Vogel, 1994). The mean value of ρ_b for electric rays was 1169±57.2, which was significantly greater than the ρ_w value of 1024 (one-tailed t-test, t=2.57, P=0.0177, N=9). During horizontal swimming, the Reynolds number (Re) ranged from 1554 to 4916, with a mean of 2848. During gliding, the Re ranged from 3219 to 5052, with a mean of 5052. Re was calculated as the product of the characteristic length (the length of the body disc), the density of seawater and the velocity, divided by the dynamic viscosity of seawater.

In sinking trials, all post mortem rays sank in trajectories that were vertical. No forward or backward gliding was observed in over 30 trials on nine neonates.

Once they initiated swimming, neonate rays swam steadily and horizontally at different depths in the water column. Neonates initiated glides when they were high in the water column. During a glide, they occasionally beat their tails.

This regression model supports the prediction that horizontal BCF swimming speed is correlated with f, as expected by elongate-body theory, but also with disc orientation and body shape, both of which are associated with lifting-body theory (Fig. 3).

This regression model supports the prediction that the St of horizontal BCF swimming is correlated with disc orientation and body shape, both of which are associated with lifting-body theory (Fig. 3).

During horizontal swimming, 10 out of a possible 36 significant (P<0.05) pairwise correlations among the candidate control variables were detected: (1) f and α_T; (2) A and α_T; (3) Ψ_A and Ψ_P; (4) Ψ_A and κ; (5) Ψ_P and α_D; (6) Ψ_P and e; (7) Ψ_P and α_D; (8) Ψ_P and ω; (9) α_D and α_T; and (10) C_D and ω. As determined by univariate regression, the response variables U_h and St were statistically correlated (adjusted r²=0.174, F_1,74=16.83, P=0.0001; U_h regressed onto St with a slope of –0.0416). Information from regressions and correlations about the control of gliding performance are integrated in a modified path diagram (Fig. 6A).

In order to model how each control variable may have simultaneous effects on multiple response variables, we created three linear parametric equations in the U_h–St performance space (Fig. 6B) using the U_h and St univariate regressions (Table 4) for α_D, α_T and Ψ_A, respectively:

We created two additional parametric equations for f using its univariate regression with U_h (Table 4), the equation for St (Eqn 2) and the minimum (0.00183 m) and maximum (0.0108 m) measured values of A, respectively:

Please note that minimum and maximum values of A were used because A was not correlated with U_h or St and both A and f occur in Eqn 2.

During all but two of the 13 glides examined, the neonates would intermittently beat their tails. These undulations, however, were never continuous throughout the trial. Thus the values of f reported for gliding should be interpreted with this intermittent beating in mind. Within-glide variation existed in all of the response and candidate control variables. Speed along the path, U_p (m s^–1), for example, had standard deviations for a trial that varied from a low of 3% to a high of 9% of the mean U_p. For this first-approximation analysis, we used the mean values over the course of a glide, thus oversimplifying the actual gliding behavior so that we could perform a steady-state analysis.

Gliding performance, as measured by U_p, was predicted by a single control variable, α_T, out of 11 possible, from the univariate-regression screening (Table 5; Fig. 7A). By multivariate-regression screening, U_p was predicted by α_T and three additional control variables, Ψ_p, α_D and ω (Table 5). The resulting multivariate control equation is as follows (adjusted r²=0.813, F_4,12=14.08, P=0.0011):

or, when the coefficients are centered by their means and scaled by half their range:

This regression model supports the prediction that gliding speed is correlated with disc orientation, disc yaw and body shape, all of which are associated with lifting-body theory (Fig. 3). We had not predicted that the flexion of the tail, a variable associated with elongated-body theory, would be correlated.

This regression model supports the prediction that glide angle is correlated with disc orientation and body shape, both of which are associated with lifting-body theory (Fig. 3). We had not predicted that f, a variable associated with elongated-body theory, would be correlated.

During gliding, five out of a possible 53 significant (P<0.05) pairwise correlations among the candidate control variables were detected: (1) A and e; (2) A and α_T; (3) C_L and C_D; (4) C_L and κ; and (5) C_D and κ. As determined by univariate regression, the response variables U_p and θ are statistically independent (adjusted r²=0.121, F_1,12=2.66, P=0.1312). Information from regressions and correlations about the control of gliding performance are integrated (Fig. 8).

From gliding data, C_L regressed onto C_D yields the following linear equation (Fig. 9A), through the origin (r²=0.419):

In negatively buoyant neonate electric rays N. brasiliensis, kinematic correlations of the motion, orientation and shape of the body with swimming performance provide evidence that electric rays combine BCF undulatory propulsion and actively modulated body lift to control locomotion. Electric rays modulate horizontal swimming speed, U_h, by varying undulatory frequency, f, the tail's posterior flexion, Ψ_P, and two variables involved with lifting-body mechanics, the body disc's pitch, α_D, and the tail's pitch, α_T (Figs 4, 6). Electric rays modulate speed on the glide path, U_p, with three of the same control variables, α_T (Fig. 7A), α_D and Ψ_P, and a different control variable, ω (Fig. 8). This overlap in the control of all four performance variables (Fig. 10A) suggests that the same underlying mechanics, modulated in different ways, produces both horizontal swimming and controlled gliding (Fig. 7). We now understand how electric rays create two different swimming modes and how they control their performance within each of those modes.

As mentioned above, one mechanism used in both swimming modes is lift generated by the body disc. The disc functions as an underwater wing, or, strictly speaking, a lifting body, in both horizontal BCF swimming and gliding. Because of the importance of a static wing's angle of attack, α_A, in determining its lift-to-drag ratio, we expected α_A to parameterize serially within the glide polar of the wing (Vogel, 1994). Instead, over α_A that ranged from 3 to 28 deg, the body disc of electric rays showed no regular relation between α_A and the coefficients of lift and drag, C_L and C_D, respectively (Fig. 9A). This is not the case in the airborne wings of birds, butterflies, bees, engineered airfoils or even the fins of air-gliding flying fish (Park and Choi, 2010). The difference between these airborne wings and the underwater wings of electric rays is that the body discs of rays change shape, altering camber, κ, of the disc. The disc's κ is an excellent predictor of the lift and drag on the ray's body disc (Fig. 9B) and, hence, creates a new type of glide polar (Fig. 9C), which, in its usual form, characterizes changes in a wing's lift-to-drag ratio as a function of angle of attack (Vogel, 1994).

Low aspect-ratio wings, with short spans relative to their chord, suffer poor performance because of induced drag; flaps on aircraft wings counter this by altering the wing's camber during take-off and landing (Arakeri, 2009). The underwater wing of electric rays likewise reconfigures by altering camber.

When electric rays swim by gliding, the most important control variable is α_T, because it predicts both U_p and the glide angle, θ, in univariate and multivariate models (Fig. 8). The importance of α_T can be seen in observations of the transition from horizontal swimming to gliding: the depression of the tail, increasing α_T, appears to initiate the glide by altering the drag on the body and levering the body disc into position to increase its lift as the animal descends on the glide path (Fig. 11). In this model, the tail, because it has only vertically oriented fins, does not function as a traditional trailing edge control surface like a horizontal elevator on the tail of an airplane. For the same reason, we also suspect that the tail of electric rays does not function like the tail of birds to reduce the parasite drag induced by the body by controlling the bound vorticity and separation (Maybury and Rayner, 2001; Thomas, 1996; Tobalske, 2007). Instead, we propose that by pitching downward, the tail increases the drag acting upon it, and that drag force is used to create leverage to control the attitude of the whole body.

Although we are proposing a new mechanism for controling the amount of lift generated by the ray's body, the generation of lift by the body is not a novel idea. Rosenberger reported observations of steady-wing gliding in the cownose ray, Rhinoptera bonasus, and the smooth butterfly ray, Gymnura micrura (Rosenberger, 2001). Lift may also come from the body, as has been shown in the smoothhound dogfish, Mustelus canis (Harris, 1936; Harris, 1938), the leopard shark, Triakis semifasciata (Wilga and Lauder, 2000) and the northern spearnose poacher, Agonopsis vulsa (Nowroozi et al., 2009). Lift from the body is thought to be a primary thrust mechanism in cownose rays (Heine, 1992).

Like the electric rays, the Atlantic guitarfish, Rhinobatos lentiginosus, is a batoid that swims with BCF undulations, increasing its U_h by increasing its f (Rosenberger, 2001). For BCF swimmers in general, f, paired with increases in tail-tip amplitude, A, at slow speeds (Bainbridge, 1958), has long been thought to be U_h's primary control variable (Wardle, 1975). The situation is not so simple for electric rays. The neonate electric rays also modulate their Ψ_P, α_D and α_T (Fig. 4). In addition, the Strouhal number, St, is a performance variable that co-varies with U_h and is modulated by similar control variables (compare Eqns 7(7b) and 8(8b); Fig. 6A); it appears that the electric ray is selecting its swimming performance based not just on a desired U_h but also on its position in the U_h–St performance space (Fig. 6B).

Using parametric equations derived from univariate regressions, nearly all of the area occupied in the U_h–St performance space of electric rays is accounted for by five control variables: f, A, Ψ_P, α_D and α_T (Fig. 6B). Interestingly, different combinations of f and A alone could account for nearly any U_h–St position that was measured. However, A is correlated with neither U_h nor St, even though it varies by a factor of nearly six, from a minimum of 0.00183 m to a maximum of 0.0108 m. At the same time, f is correlated only with U_h (Fig. 6A).

To traverse the U_h–St space diagonally, either α_D, α_T or Ψ_A suffice over a limited range (Fig. 6B); all three appear to be involved in the transition from the optimum St range, where 46% of the trials and the greatest range of U_h occurred, to a region of low U_h and high, less-than-optimal St. Because α_D and α_T help control lifting-body properties (Fig. 2), it appears that disc-generated lift is most important at low U_h from approximately 0.04 to 0.06 m s^–1, the same range, interestingly, that covers the bulk of path speeds measured during gliding (Fig. 7A).

The importance of lift at low speeds is at first paradoxical, given that lift is proportional to the square of velocity (Eqn 3). However, because most of the slow U_h values have St values above the optimal range, it is possible that lift from the disc is providing thrust either to augment thrust lost by the inefficient interaction of shed vortices, indicated by high St, or to help stabilize the body with the loss of inertial stabilization. Electric rays may swim at low speeds with instability in order to enhance maneuverability, as has been shown in sea lions, Zalophus californianus (Fish et al., 2003).

The electric ray serves as a target for bioinspired design (Krishnamurthy et al., 2009). The self-propelled RayBot has a static body disc and a BCF-style undulatory tail (Krishnamurthy et al., 2010). Even though RayBot can swim and glide, its performance is limited: it can only modulate f and A to control U_h (Long, 2011). To increase its range of locomotor abilities, RayBot would clearly benefit from the addition of two key features possessed by electric rays: (1) a cambering body disc (see Fig. 9) and (2) a pitching tail (see Fig. 11). Both of these features would allow RayBot to modulate its performance as an underwater wing.

Adjustments to lift have been successfully achieved in a cephalic wing used to control vertical movements in a self-propelled fish-like robot (Mohammadshahi et al., 2008). More generally, hydrodynamic flow control is essential in allowing marine mammals to maneuver (Fish et al., 2003; Fish et al., 2008). Unlike marine mammals, however, electric rays do not swim or accelerate rapidly. Instead, in terms of biomimetic designs, electric rays offer advantages for a different set of engineering goals for underwater robots (Bandyopadhyay, 2005; Fish, 2006): reduced energy consumption, a high level of stability and stealth. Electric rays reduce energy consumption by using gliding to move horizontally, in much the same way that engineered underwater gliders do (Rudnick et al., 2004). In terms of stability and stealth, electric rays are slow-moving benthic feeders, punting (Koester and Spirito, 2003) along the bottom in search of prey and remaining in contact with the substrate during feeding (Dean and Motta, 2004). This stable, low-profile behavior is made possible because of the ray's negative buoyancy, which also gives them the ability to hold station passively. Slow movements near a solid surface provide for hydrodynamic stealth, keeping fluid perturbations well within the boundary layer.

Interest in batoid-inspired underwater robots of various types is high. The stingray's undulatory body disc is a target of interest (Clark and Smits, 2006; Low and Willy, 2006), as is the oscillatory wing of manta rays (Punning et al., 2004; Gao et al., 2007; Wang et al., 2009). Robotic fins and the robots they propel provide important tools not only for engineers but also for biologists testing functional hypotheses (Lauder et al., 2007; Long, 2007). Engineers looking for biomimetic principles (sensuFish, 2006) that can be applied to biologically inspired underwater robots may appreciate the mathematical control laws provided here that link the motor output of specific structures with the locomotor performance of the whole swimming ray.

We thank the Vassar Animal Care and Biology Department staff for their support on this project. Jonathan Hirokawa and Sonia Roberts generously devoted their time to digitizing videos. Mason Dean and Laura Macesic provided useful consultation on neonate Narcine husbandry. We thank two anonymous reviewers who helped improve the quality of this manuscript.