## SUMMARY

Unlike most batoid fishes, electric rays neither oscillate nor undulate their body disc to generate thrust. Instead they use body–caudal–fin (BCF) locomotion. In addition, these negatively buoyant rays perform unpowered glides as they sink in the water column. In combination, BCF swimming and unpowered gliding are opposite ends on a spectrum of swimming, and electric rays provide an appropriate study system for understanding how the performance of each mode is controlled hydrodynamically. We predicted that the dorso-ventrally flattened body disc generates lift during both BCF swimming and gliding. To test this prediction, we examined 10 neonate lesser electric rays, *Narcine brasiliensis*, as they swam and glided. From video, we tracked the motion of the body, disc, pelvic fins and tail. By correlating changes in the motions of those structures with swimming performance, we have kinematic evidence that supports the hypothesis that the body disc is generating lift. Most importantly, both the pitch of the body disc and the tail, along with undulatory frequency, interact to control horizontal swimming speed and Strouhal number during BCF swimming. During gliding, the pitch of the body disc and the tail also interact to control the speed on the glide path and the glide angle.

## INTRODUCTION

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

## MATERIALS AND METHODS

### Animals

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.

### Swimming trials

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.

### Physical properties of *post mortem* neonates

*Post mortem*, we measured the body density, ρ

_{b}(kg m

^{–3}), of each neonate. The mean value of ρ

_{b}for the neonates was used to calculate the buoyant mass,

*M*

_{b}(kg), of each individual: where

*M*(kg) is the mass of the body in air and ρ

_{w}(kg m

^{–3}) is the density of seawater. The ρ

_{b}of each

*post mortem*neonate was measured using the volume-displacement method. Each specimen was weighed and then ρ

_{b}was calculated as the ratio of

*M*to volume. The mean of the neonates' ρ

_{b}was compared with ρ

_{w}, 1024 kg m

^{3}at 20°C (Vogel, 1994), using an independent one-sample

*t*-test (JMP, version 8.0, Cary, NC, USA).

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.

### Kinematic analysis

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.

### Performance variables

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*_{2}, 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*_{2} yields a more conservative estimate of *U*_{h} because it varies less than *m*_{1} or *m*_{3}. The *U*_{h} over all trials was significantly lower for *m*_{2} 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*_{2} (*F*_{2,268}=65.71, *P*<0.0001). During gliding, the lowest variance also occurred at *m*_{2}, as well as *m*_{3} (*F*_{2,39}=7.31, *P*=0.002).

*St*, the dimensionless frequency of vortex shedding at the caudal fin, which is calculated as follows (Rohr and Fish, 2004): where

*A*(m) is the mean lateral peak-to-midline amplitude of the tip of the tail over three tailbeat cycles and

*f*(Hz) is the undulatory frequency measured at the tip of the tail over three tailbeat cycles. Experiments on foils, fish and dolphins show that the optimum

*St*range for propulsion is from 0.20 to 0.35 (Triantafyllou et al., 2002). In cetaceans, the

*St*range of 0.225 to 0.275 correlates with peak swimming efficiency (Rohr and Fish, 2004).

### Candidate control variables

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*_{3}, *m*_{4}, *m*_{5} and *m*_{4}, *m*_{5}, *m*_{6}, 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*_{1} and *d*_{3} 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*_{3} 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*_{1}, *t*_{1} and *t*_{1}, *t*_{2}; the camber of the disc, κ (deg), as the angle between points *d*_{1}, *d*_{2}, *d*_{3}; 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.

### Lifting-body properties during gliding

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.

*C*

_{L}, we used the following equation (Vogel, 1994): where

*S*

_{w}(m

^{2}) is an approximation of the disc's wetted surface area, calculated as twice the product of disc's length and width, both of which were measured in each trial using photographs digitized in ImageJ (version 1.14, NIH, Bethesda, MD, USA) (Table 1). The use of

*S*

_{w}rather than a planform area decreases coefficient estimates by approximately one half;

*S*

_{w}is used for normalization in eel propulsion (Tytell and Lauder, 2004). The variable

*L*(N) is lift, which was calculated as follows (Bishop, 2007): where

*R*(N) is the hydrodynamic resultant calculated as the product of

*M*

_{b}, defined in Eqn 1, and gravitational acceleration,

*(Fig. 1B). Because this is a steady-state analysis, and gliding speed varies, we compared inertial forces with lift forces. In the worst case, when the largest rays are paired with the largest accelerations (estimated from change in velocity), inertia is 11% of lift.*

**g***C*

_{D}, we used the following equation (Vogel, 1994): where

*D*is the steady-state drag, calculated as follows (Bishop, 2007):

### Statistical analyses

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*^{2}. 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.

### Parametric equations

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.

## RESULTS

### Physical properties of the neonates

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.

### Horizontal swimming

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.

*U*

_{h}, was predicted in univariate-regression screening by four of nine possible control variables (Table 4; Fig. 4):

*f*, Ψ

_{p}, α

_{D}and α

_{T}. By multivariate-regression screening,

*U*

_{h}was predicted by two significant control variables,

*f*and α

_{D}, and one marginally significant (0.05<

*P*<0.10) control variable, α

_{T}. Thus the screening produces the following multivariate equation (adjusted

*r*

^{2}=0.631,

*F*

_{3,73}=44.08,

*P*<0.0001): or, when the coefficients are centered by their means and scaled by half their range:

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

*St*, was predicted in univariate-regression screening of two of out of seven possible control variables (Table 4; Fig. 5): Ψ

_{A}and α

_{D}. By multivariate-regression screening,

*St*was predicted by two control variables, one of which was different: α

_{D}and α

_{T}(Fig. 5). The resulting multivariate control equation is as follows (adjusted

*r*

^{2}=0.171,

*F*

_{2,72}=44.08,

*P*<0.0001): or, when the coefficients are centered by their means and scaled by half their range:

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*

^{2}=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.

### Gliding

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*^{2}=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.

*f*, α

_{T}and ω. Because of the small sample size associated with gliding, we also report the marginally significant (

*P*=0.055) control variable α

_{D}. By multivariate screening, θ was predicted by three significant control variables,

*f*, α

_{T}and

*e*, and one non-significant control variable, α

_{D}. The resulting multivariate control equation is as follows (adjusted

*r*

^{2}=0.702,

*F*

_{4,12}=8.07,

*P*=0.0065): or, when the coefficients are centered by their means and scaled by half their range:

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*^{2}=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).

### The body disc as an underwater wing

From gliding data, *C*_{L} regressed onto *C*_{D} yields the following linear equation (Fig. 9A), through the origin (*r*^{2}=0.419):

_{A}along the regression line (Vogel, 1994). α

_{A}does not, however, plot linearly along the regression line (Fig. 9A). Alternatively, from the path diagram (Fig. 8), we know that κ is significantly correlated with both

*C*

_{L}and

*C*

_{D}, and these linear relationships (Fig. 9B) produce the following regressions: with

*r*

^{2}values of 0.499 and 0.479, respectively. By solving for κ, Eqns 13 and 14 can be combined to produce a glide polar equation in which κ is linearly parameterized over the range of κ measured in these experiments (Fig. 9C, maroon line):

## DISCUSSION

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.

### The ray's body disc is a reconfiguring wing

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.

### The ray's tail initiates and controls gliding performance

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

### Electric rays control swimming speed and performance in an unusual way

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

### Design of electric-ray-inspired underwater robots

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 (*sensu*Fish, 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.

## FOOTNOTES

This work was funded by NSF DBI-0442269 and IOS-0922605.

## Acknowledgements

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.