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 dorsoventrally 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 lifttodrag 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 bodygenerated lift will be maximized by winglike shapes, we hypothesized that dorsoventrally 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.
Dorsoventrally 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 dorsoventrally 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–pairedfin 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–pairedfin mode (Rosenberger and Westneat, 1999), relatively little is known about the swimming of batoid species using the BCF mode.
For BCFswimming 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. 0810B 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 HC5 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: (1) 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 volumedisplacement 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 onesample ttest (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 middisc body point, m_{2}, from the ventral view (mean U_{h} from ventral and lateral views did not differ statistically). The middisc 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 midpoints (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).
For horizontal BCF swimming, we also calculated the Strouhal number, St, the dimensionless frequency of vortex shedding at the caudal fin, which is calculated as follows (Rohr and Fish, 2004): (2) where A (m) is the mean lateral peaktomidline 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).
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_{2} from the lateral view. For gliding, we also calculated the mean glide angle by trigonometry using U_{h} and U_{p} (Fig. 1B,C).
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) largeamplitude elongatedbody 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 thrustpower 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).
Liftingbody theory (Harris, 1936; Vogel, 1994) inspired seven more candidate control variables (Fig. 3 and Table 3). Crucial to the liftanddrag 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 nonzero θ 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 windtunnel 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.
Liftingbody 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 steadystate 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 selfpropelling 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 steadystate method provides a first approximation of the lift and drag coefficients.
To calculate C_{L}, we used the following equation (Vogel, 1994): (3) 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): (4) where R (N) is the hydrodynamic resultant calculated as the product of M_{b}, defined in Eqn 1, and gravitational acceleration, g (Fig. 1B). Because this is a steadystate 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.
To calculate C_{D}, we used the following equation (Vogel, 1994): (5) where D is the steadystate drag, calculated as follows (Bishop, 2007): (6)
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 leastsquares 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, univariateregression and multivariateregression 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, univariateregression 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, multivariateregression screening, we separately regressed each performance variable onto the entire pool of candidate control variables using a twostep 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 socalled `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 (onetailed ttest, 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.
Horizontal BCF swimming performance, as measured by U_{h}, was predicted in univariateregression screening by four of nine possible control variables (Table 4; Fig. 4): f, Ψ_{p}, α_{D} and α_{T}. By multivariateregression 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): (7a) or, when the coefficients are centered by their means and scaled by half their range:
(7b)This regression model supports the prediction that horizontal BCF swimming speed is correlated with f, as expected by elongatebody theory, but also with disc orientation and body shape, both of which are associated with liftingbody theory (Fig. 3).
Horizontal swimming performance, as measured by St, was predicted in univariateregression screening of two of out of seven possible control variables (Table 4; Fig. 5): Ψ_{A} and α_{D}. By multivariateregression 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): (8a) or, when the coefficients are centered by their means and scaled by half their range:
(8b)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 liftingbody 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:
(9) (10) (11)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:
(12) (13)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. Withinglide 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 firstapproximation analysis, we used the mean values over the course of a glide, thus oversimplifying the actual gliding behavior so that we could perform a steadystate analysis.
Gliding performance, as measured by U_{p}, was predicted by a single control variable, α_{T}, out of 11 possible, from the univariateregression screening (Table 5; Fig. 7A). By multivariateregression 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):
(14a)or, when the coefficients are centered by their means and scaled by half their range:
(14b)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 liftingbody theory (Fig. 3). We had not predicted that the flexion of the tail, a variable associated with elongatedbody theory, would be correlated.
Gliding performance, as measured by θ, was predicted by three out of 11 possible control variables in the univariateregression screening (Fig. 7B): 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 nonsignificant control variable, α_{D}. The resulting multivariate control equation is as follows (adjusted r^{2}=0.702, F_{4,12}=8.07, P=0.0065): (15a) or, when the coefficients are centered by their means and scaled by half their range:
(15b)This regression model supports the prediction that glide angle is correlated with disc orientation and body shape, both of which are associated with liftingbody theory (Fig. 3). We had not predicted that f, a variable associated with elongatedbody 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):
(16)We created a glide polar by plotting the two coefficients and then attempting to parameterize α_{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: (17) (18) 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):
(19)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 liftingbody 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 lifttodrag 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 airgliding 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 lifttodrag ratio as a function of angle of attack (Vogel, 1994).
Low aspectratio 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 takeoff 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 steadywing 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 tailtip 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 covaries with U_{h} and is modulated by similar control variables (compare Eqns 7 and 8; 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, lessthanoptimal St. Because α_{D} and α_{T} help control liftingbody properties (Fig. 2), it appears that discgenerated 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 electricrayinspired underwater robots
The electric ray serves as a target for bioinspired design (Krishnamurthy et al., 2009). The selfpropelled RayBot has a static body disc and a BCFstyle 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 selfpropelled fishlike 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 slowmoving 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, lowprofile 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 batoidinspired 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.
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.
FOOTNOTES

This work was funded by NSF DBI0442269 and IOS0922605.
 © 2011.