How zebrafish turn: analysis of pressure force dynamics and mechanical work

Animals in different environments have evolved a variety of locomotor strategies for specific purposes. In water, many fish navigate and propel themselves forward using body undulations, and most studies in the past have focused on steady swimming of fish to better understand how patterns of body deformation contribute to thrust generation (e.g. Lauder, 2015; Videler, 1993; Webb, 1975). Corresponding measurements of metabolic rate via respirometry in flow tanks have contributed to our understanding of how fish use metabolic substrates and how fish metabolism varies with swimming speed and distance (van Ginneken et al., 2005; Blank et al., 2007; Di Santo et al., 2017; Korsmeyer et al., 2002; Sepulveda et al., 2007). However, many fishes also routinely perform unsteady maneuvers during navigation, feeding, foraging or predator avoidance. The generation of locomotor forces and the production of mechanical work by the body during such maneuvers is largely unknown in part because they are difficult to study with classical techniques such as respirometry (but see Schakmann et al., 2020; Roche et al., 2014). Given that most fish perform unsteady maneuvers repeatedly throughout their lives, it is important to understand mechanical work production by the body and its effect on the surrounding fluid.

Previous analyses of the pattern of force production by the body during steady rectilinear locomotion have used a work loop technique (Johnson et al., 1994; Josephson, 1985) and analysis of the timing of muscle activation relative to body bending kinematics to assess the extent of positive and negative work done along the body axis by segmental red and white myotomal fibers (Ellerby et al., 2001; Jayne and Lauder, 1995; Rome et al., 1993). Studies by Coughlin (2000) and Coughlin and Rome (1996), among other work on this topic, have shown that most net-positive locomotor power in fish is provided by middle and posterior regions of segmental body musculature during steady swimming.

Swimming locomotion emerges as a result of the interaction between the body and fluid (Tytell et al., 2010). Corresponding interaction forces, which originate from muscular forces moving the body, do mechanical work that results in an energy exchange between the body and the fluid. By estimating patterns of pressure change on the fish body surface that result in water movement either away from or toward the fish, we can estimate the exchange of energy along the body and sum these effects to estimate total energetic costs for locomotor maneuvers. By analyzing the magnitude of these effects along the body in conjunction with body kinematics, we can estimate which body regions do positive or negative work on the fluid, and make such estimates in a time-dependent manner throughout the maneuver. This general approach has the advantage of not requiring electrode implantation into body musculature at multiple sites in small fish such as zebrafish that otherwise are an excellent experimental system in which to study unsteady locomotor behaviors. However, to accomplish such pressure-based analyses it is necessary to obtain detailed measurements of both body motion (kinematics) and forces (dynamics) that result from body surface pressure changes, and surface pressures need to be estimated around the swimming fish throughout the maneuver. Fortunately, recent work has provided a method to successfully derive pressure fields from velocity fields in freely moving animals and has validated this approach against canonical cases (Dabiri et al., 2014; Gemmell et al., 2015; Lucas et al., 2017, 2020). This paves the way for determining distributed forces along the body around swimming fish during maneuvers, estimating longitudinal patterns of positive and negative work, and quantifying overall dynamics and energetics during routine turning maneuvers in live animals. With this method, we can in particular identify periods of negative work, in which the body movements are supported by the fluid and in addition we can determine where along the body this occurs.

Previous computational and kinematic studies of unsteady locomotion have shown that during rapid maneuvers, in contrast to steady swimming, the posterior body region in fish is bent in a direction opposite to the direction of movement of that body area (Borazjani et al., 2012; Jayne and Lauder, 1993; Tytell and Lauder, 2008). This strongly suggests that the high degree of fluid loading on the posterior body region during rapid movements is greater than the force generated by the posterior body musculature, which indicates negative work (resisting loading) in this axial location. But experimental estimates of longitudinal work and fluid–structure interactions in unsteady fish maneuvers are not currently available to address this question.

Therefore, in this study, we focused on testing the hypothesis that unsteady maneuvers in zebrafish involve the production of negative work in specific body regions, thus allowing the body to move supported by the fluid in periods of the turn. In order to accomplish this, we first present a modified flow visualization technique based on particle image velocimetry (PIV) that uses infrared light and small algae particles as tracers to provide a detailed two-dimensional visualization of fluid motion around the body of maneuvering zebrafish. We then compute pressure fields around the body and the pattern of axial distributed forces and torques to paint a picture of longitudinal body dynamics during turning. Using these data we then estimate local and global positive and negative work contributions of different body segments, and compare our results with results on steady swimming and measurements of metabolic cost obtained in previous work in swimming zebrafish.

In the following, we describe the workflow (Fig. 1) to estimate the hydrodynamic forces along the body of a freely swimming fish. This process involves the extraction of velocity fields by means of a modified PIV technique (Fig. 2) and the computation of pressures along the body boundaries which are multiplied with corresponding surface area segments.

Four wild-type zebrafish, Danio rerio (F. Hamilton 1822) (∼3–6 months old), with an average total length of 22.4±2.4 mm (mean±s.d.) and an estimated mass of 37.8±10.5 mg were used for the experiments in this study. The mass of each fish was computed with the model m=4.14×10⁻⁶ BL^3.17 from McHenry and Lauder (2006) considering body length as BL=0.79×L (from measurements of the side view) of the total length L, thus excluding the caudal fin. We report the mass of each fish with corresponding 95% confidence intervals in Table S1. Fish were recorded in room darkness performing spontaneous turns at Reynolds numbers of Re=685±199 (maximum 984, minimum 259), based on Re=UL/ν with the speed U defined as center of mass displacement over the turning duration, total length L and kinematic viscosity of v=10⁻⁶ m² s⁻¹. Experiments were carried out at a water temperature of 28°C. All experiments were performed in accordance with Harvard animal care and use guidelines, IACUC protocol number 20-03 to G.V.L. We recorded a total number of 80 turning maneuvers, of which N=18 were selected for the final complete analysis of kinematics and dynamics. We selected the turns based on two criteria: (1) only turns where fish swam in a plane parallel to the water surface were selected to measure accurate flow fields (see ‘PIV’, below); (2) only turns where fish had their pectoral fins held against the body were selected, as we were interested in the contribution of body undulations during a turn. Individual zebrafish exhibit large intra-individual variation in the speed and angle of turns and we recorded a variety of such turns from each individual. While turns ranged from small-angle to more rapid large-angle maneuvers, recording these latter turns was significantly more problematic given our requirement of analyzing only turns that occurred within a single horizontal plane as determined with a second lateral view camera (see below). Large-angle turns frequently involved out-of-plane motion that rendered them unsuitable for analyzing body kinematics and calculating pressure fields. Large-angle turns are particularly interesting as they approach the timing and body amplitude motions seen in zebrafish C-starts, and inclusion of these turns, even with a reduced sample size, allows comparison of these large-angle maneuvers with previous work on C-start escape responses.

To measure the interactions of the surrounding fluid with a turning zebrafish, we employed PIV. Classic PIV approaches to analyze fish locomotion (Drucker and Lauder, 1999; Müller et al., 1997) involve the illumination of particles in a laser sheet. However, the use of a laser comes with two key drawbacks: (1) non-transparent objects in the fluid cast shadows and occlude the movement of particles therein, which often requires the use of two aligned laser sheets (e.g. Lucas et al., 2020); (2) bright laser light sheets can alter fish behavior, especially in the case of smaller fish such as zebrafish, where laser light sheet thickness is of the order of 20–30% of body depth, and our goal was to study spontaneous turning behavior. Behavioral alteration of fish swimming behavior by laser light is insufficiently studied, and we observed zebrafish avoiding laser light sheets in preliminary experiments. In an effort to overcome these drawbacks, we tested and used an alternative technique that is inspired by PIV at the micro-scale (Gemmell et al., 2014). Instead of recording isolated brightly illuminated particles that appear in a laser sheet, we recorded particles in a narrow optical depth of field using infrared illumination only (not visible to zebrafish, and not using a laser light sheet), thus capturing only the movements of particles within a thin slice of fluid. In addition, we used algal cells as particles, and avoided introducing larger silver or plastic particles into the testing arena. Using an array of infrared lights (Fig. 2) also avoided shadows that are often present in laser-imaged PIV studies and provided full-field illumination of the flow around the body of individual fish.

The experimental setup in Fig. 2 shows the arrangement of the high-speed cameras (Photron FASTCAM Mini AX200 for the ventral view, and Photron AX50 for the side view). We used a Nikon Micro Nikkor, 105 mm, f/2.8 macro lens and recorded with a resolution of 1024×1024 pixels, a frame rate of 1000 Hz and a shutter speed of 1/1000 s. The use of a macro lens allowed us to obtain sharp images at a close distance, while the small f-stop decreased the depth of field for a given field of view (a smaller field of view will decrease the depth of field). In our experiments we captured a field of view of 44 mm×44 mm. We placed fish in a clear plastic water tank (13 cm×8.5 cm, and 3 cm water depth), seeded the water with neutrally buoyant algae particles (Tetraselmis, ∼12 μm), and illuminated them by means of infrared lights. The reflected light of the moving particles in the fluid along with the fish body could then be recorded by the cameras. The use of infrared lights allowed us additionally to brightly illuminate the fish body, which greatly simplified the later image segmentation of body and fluid.

In addition to the ventral view, we recorded fish swimming from a side view to select for maneuvers that occurred exclusively in the plane of focus. We analyzed the undulatory mechanism of turning in this study and thus selected only maneuvers with pectoral fins closely held to the body, and in which the fish stayed precisely within the horizontal plane of focus throughout the turn.

We processed the high-speed videos from the ventral view in DaVis 8.3 (LaVision Inc.), by automatically masking the fish body (Fig. 1B) in the fluid with a series of filters. Subsequently, we obtained velocity fields of size 64×64 vectors (Fig. 1C) using cross-correlation. The computation of the masks served two purposes: (1) to obtain more accurate velocity vectors along the fish body boundaries and (2) to use them for the pressure computations, which requires the exclusion of non-fluid points in the analyzed velocity field.

Based on the velocity fields, we determined the corresponding pressure fields (Fig. 1E) in the fluid using Queen 2.0 software (Dabiri et al., 2014) running in MATLAB (Mathworks, Inc.). The pressure fields were obtained from paths of integrated pressure gradients that are estimated based on consecutive velocity fields. The proposed method is suitable for our problem, as it is able to cope with fluid–solid interfaces based on the median averaging that is used across different integration paths. Previous research has validated the pressure fields obtained using this approach for both static and dynamic examples in aquatic locomotion (Dabiri et al., 2014; Gemmell et al., 2015; Lucas et al., 2017, 2020). The influence of the grid size on the spatial convergence of the underlying algorithm is discussed in appendix 1 of Dabiri et al. (2014). Whereas too large a grid spacing decreases the accuracy in the pressure estimations, smaller grid spacing increases the accuracy until a limit is approached that is due to inherent model errors and numerical round-off errors. Our grid spacing was determined by the resolution of the velocity vector field and led to roughly 30 pressure points along the length of the fish body. This ensured sufficient accuracy in the pressure estimations (see appendix 1 of Dabiri et al., 2014: D/16 leads to an error of 5% in instantaneous pressure, where D is body dimension).

Hydrodynamic forces acting on the swimming fish are determined by the pressure at the interface between the fluid and body. Therefore, it was crucial to precisely define this interface (Fig. 1D) in each frame, as well as in the transition between frames. In particular, it is important to smooth the body outline to minimize noisy pressure estimates at the interface. We automatically processed the interfaces according to these criteria in MATLAB: (i) import of masks from DaVis software to MATLAB; (ii) distance transform of the masks to initially smooth and ensure that all solid parts are excluded from the fluid (use bwpropfilt and bwboundaries functions from the image processing toolbox in MATLAB); (iii) resampling of body outline points to 500 points in each frame; and (iv) smoothing of body outlines with Savitzky–Golay filter (3rd order, 41 frame length).

Using the Queen 2.0 software, we calculated the pressure fields from consecutive velocity fields that were spaced 5 ms apart. We chose this larger sampling time (recordings at 1000 Hz), as we encountered noisy pressure estimations due to numerical errors for lower sampling times (see appendix 2 of Dabiri et al., 2014). Furthermore, we included the viscous term in the pressure calculations (pressure gradient is dependent on material acceleration and viscosity related to the local fluid velocity curvature), and used the temporal spline smooth filtering of input velocity fields, and the nearest neighbor spatial smoothing of both pressure gradient and pressure fields.

Once the pressure field was computed, we determined the pressure at the fluid–structure interface (body outline) based on a nearest neighbor approach (Fig. 3).

We divided the fish body into segments to analyze the local interaction dynamics with the fluid and the kinematics of turning. For this purpose, we manually tracked body midline using a custom-written MATLAB script (Fig. 1F). For each video frame, we placed seven modifiable markers along the body, which were automatically connected with a spline fit. The midline was then resampled to 1000 points. This way, we traced midlines with a sampling time of 10 ms, and subsequently linearly interpolated between frames to maintain midline estimates at 1000 Hz.

The estimation of interaction forces F from the pressure p at the fluid–structure interface requires inclusion of a corresponding surface area A, as by definition F=p×A. Furthermore, we made the assumption that pressure is equally distributed along the depth of the body, and thus that forces could be obtained by multiplying the computed pressure at the mid-section of the body with corresponding surface area slices perpendicular to the midline. This approach has been discussed further and validated by Lucas et al. (2017) (see also Lucas et al., 2020).

We used a similar approach to that in McHenry and Lauder (2006) to estimate the surface area at a given point at the fluid–structure interface. The surface area corresponding to a point at the fluid–structure interface was therefore computed as the area of a half-ellipse cylinder (Fig. 4). The half-widths of the ellipses were estimated based on the depth and width profiles of the fish, which were obtained from side and ventral view photographs using a custom-written MATLAB script. Note that the mapping of corresponding body depth and width at a given point on the fluid–structure interface was computed for each frame separately as a deformable body was considered. We validated this approach (Table S1) by computing the total surface area for each of our four individuals and comparing them with the model in McHenry and Lauder (2006).

For each point on the fluid–structure interface, we computed force based on pressure and segment slice area that were estimated as described before. The force direction was determined by the normal direction to the body–fluid interface, because pressure-based forces act normal to the surface. The underlying assumption for this approach is that shear forces are small (Lucas et al., 2017, 2020), and that the hydrodynamic interaction forces largely result from pressure-based effects. This assumption is well justified for swimming at larger Re numbers (Re≈10,000). Our turning experiments appeared at intermediate Re numbers in the range of several hundred, but still support the omission of shear forces based on a computational fluid dynamics simulation from Li et al. (2012). In that study, the authors simulated a C-start of a larval zebrafish (Re=550) and reported both pressure and shear stresses during this maneuver. They showed that shear forces were smaller than pressure forces and localized the main contributions of shear at the head during the preparatory stroke acting as drag forces. They further found shear contributions at the tail, although much smaller and acting as thrust. These findings justify the omission of shear stresses in our computation of hydrodynamic forces and allow us to capture the main dynamics of the considered spontaneous turning maneuvers.

In the following, we refer to the turn side as the lateral side of the body that faces the direction of the turn and call the contralateral side the opposite side.

We present our results in the different sections with respect to 20 segments along the anteroposterior body axis. Segment 1 represents the anterior-most segment at the head, and segment 20 the posterior-most segment including the tail. We show the influence of the number of segments regarding the power computations in our figshare dataset (https://doi.org/10.6084/m9.figshare.c.4986179).

Turning is part of the burst-and-coast behavior in zebrafish. In this context, turns happen between gliding periods, in which the body displays a straight long-axis posture with close to zero body bending (straight line) from head to tail. Consequently, we can define the start of a turn as the time instance at which movement away from the straight-line posture of the body is first detectable in the high-speed videos. We defined the end of a turn as the time instance where turning movement stopped and the body bending remained still from frame to frame.

Following the definition of the start and end of a turn, we defined and measured the turning angle between the starting straight-line position and the ending straight-line position of the fish in all turns.

Bending angle was defined as the relative joint angle between segments along the anteroposterior body axis. Corresponding angular speed was obtained by differentiation after smoothing with a moving average filter.

According to the definition of power in Eqn 1, positive power refers to the rate of mechanical work done by the body on the fluid. Negative power refers to the rate of mechanical work done by the fluid on the body. The same definitions apply for mechanical work.

We defined the total work as the sum of absolute positive and absolute negative work. Net work was obtained via the integration of mechanical power over time.

We use the interquartile range (IQR) as a measure of variability around the median value of a distribution. The IQR is defined as the distance between the 25% and 75% quantile.

The D. rerio in our experiments were recorded swimming in a burst-and-coast style. This behavior is characterized by repeated sequences of full-body undulatory propulsive strokes with changes in heading (turning) followed by straight-line gliding periods. Our focus was on the turning portion of this behavior, which we could readily isolate by analyzing the period of heading change starting from a straight-line body posture back to a final straight-line position.

A change of heading direction requires an inherently asymmetric behavior that we could detect in the kinematics of the body and in body dynamics (see below). Fig. 5 highlights these asymmetries appearing in the midline kinematics in the recorded turns. They become especially apparent in the representation with all the midlines aligned at the head (Fig. 5A–C, bottom row). Two phases can be distinguished. Phase 1: bending of the mid-body to one side, and the posterior part to the other side; overall, the body forms an S-shaped midline. Phase 2: extension back to the straight position, but with all body parts bent to the same side; overall, the body forms a J-shaped midline.

We further report the overall characteristics of turning angle and duration in Fig. 5 and found that our selected experiments (for the kinematic and dynamic analysis) are a good representation of the overall turning behavior for all the sequences we recorded. For our analysis, we did not distinguish between left or right turns and considered only the magnitude of the turning angle. We found that 90% of all the recorded experiments (N=80) had a turning angle below 70 deg, and that the corresponding distribution peaked at 21.4 deg, indicating the preferred turning angle of zebrafish in our study. Only a weak positive linear correlation (P=0.006, R²=0.09) between turning angle and duration was found.

We noticed different dynamical features in our analysis of frequently occurring turns below turning angles of 70 deg compared with turning angles above 70 deg (Table S2), and therefore show one representative example of a low- and high-magnitude turn in Fig. 6. Complete snapshot sequences (every 15–20 ms) of both velocity and pressure fields for all investigated turns are provided in figshare (https://doi.org/10.6084/m9.figshare.c.4986179) (see also Fig. 6).

The patterns of fluid flow are shown in Fig. 6A,C. We consistently observed a bow wake at the head of the fish in all the experiments, which is indicated by the velocity vectors pointing away from the anterior body.

A second characteristic feature that appeared in all turns was the strong flow of water towards the body on the turn-side during phase 1 (bending of mid-body) of the turn (Fig. 6A2,3,C2,3). Interestingly, this flow pattern resulted in different pressures at the mid-body section, as we describe in ‘Pressure distribution’, below.

The velocity fields further revealed the formation of two vortices behind the turning fish. The first of these was shed into the fluid at the end of phase 1 when the tail changed movement direction from the turn side to the opposite side (Fig. 6A3,C3). The second vortex appeared at the end of the turn with the tail flick towards the straight-line position (Fig. 6A5,C5).

Zebrafish turns at all magnitudes of angular excursion showed similar overall patterns of body surface pressure change. On the turn side, a strong suction zone was created in the middle of the body and positive pressure zones emerged at the anterior and posterior ends of the body (Fig. 6B2,D2). The strong negative pressure in the middle on the turn side is a consequence of the mid-body bending to the opposite side, therefore creating a region for water to flow in toward the body. As the turn progressed, the water continued to flow towards the mid-body on the turn-side; however, as the bending movement slowed down, water began to collect against the body. Consequently, the pressure in the middle of the body on the turn side changed from negative to positive.

Furthermore, we found positive pressure regions at the head of the fish in phase 2 of the turning maneuvers. This confirmed the presence of the bow wake that we had detected in the velocity field before. Another pressure characteristic at the anterior turn-side part of the body involved the transition from a positive pressure at the beginning to a negative pressure as the turn progressed. This was a result of the anterior body first pushing the water to the side as it accelerated and transitioning to be pulled by the water as it decelerated. In phase 2 of the turn, we observed a distinct pressure pattern at the tail. This phase was characterized by a tail flick, where the posterior part of the body was bent perpendicular to the rest of the body and flicked back to a straight line. In the process, a positive pressure region was created at the posterior end of the body, which at the same time generated forces pointing to the final heading direction (Fig. 6B4,5,D4,5).

In addition to the flow and pressure fields, we analyzed segmental forces along the body and the time evolution during different turning maneuvers (Fig. 7). We found that the force dynamics were qualitatively different as turn angle increased. For small-angle turns, we found a sequential timing of force maxima (Fig. 7, left). In other words, the forces along the body showed a traveling wave pattern on both the turn side and the opposite side of the body. However, the large-angle turns showed a double-burst pattern where force maxima were reached in bursts around the same time, once in the first half and once in the second half during the turning maneuver (Fig. 7, right).

Whereas the force patterns showed varying characteristics with turn magnitude, the bending patterns (bending angle over time) appeared as traveling waves for all turns (Fig. 7). Fig. 7 also shows that the body bending pattern changed with the amplitude of the bending angle. To understand the influence of bending on turning angle in different parts along the anteroposterior body axis, we looked at corresponding correlations of maximum bending angle and angular speed (Fig. 8A,B). We found positive correlations between turning angle and maximum bending angle in the anterior (P<0.001, R²=0.69) and middle section (P<0.001, R²=0.64) of the body. No correlation was found between maximum bending angle in the posterior body and turning angle (P=0.66, R²=0.01). Our analysis also revealed positive correlations between turning angle and angular speed of bending in the anterior (P<0.001, R²=0.82) and in the middle part (P=0.004, R²=0.42) of the body. Again, no corresponding correlation was found in the posterior body (P=0.83, R²=0.003).

To provide insight on the local stresses that are applied on the body from the fluid during a turning maneuver, we analyzed the maximum hydrodynamic forces. To obtain a finer resolution along the body, this analysis was carried out for 20 body segments. As shown in Fig. 8C, we measured maximum force peaks in the middle of the body (at 50% total length) and in the posterior body (at 90% total length). The posterior force peak occurred in a more variable way (IQR=0.30, s.d.=0.17) compared with the mid-body force peak (IQR=0.09, s.d.=0.06), which appeared consistently across all the selected turning maneuvers.

Based on the force computations and the estimated midline kinematics, we quantified the mechanical power that was exchanged between the fish and fluid in order to obtain a measure for the mechanical cost to perform a turning maneuver.

We first describe detailed spatial and temporal power distributions for two representative turns, as shown in Fig. 9, which illustrates stacked power curves of the different segments along the anteroposterior fish body axis. Power was analyzed with respect to two aspects: (1) positive power, which describes the rate at which the body does mechanical work on the fluid and increases energy of the fluid in the process; and (2) negative power, which describes the mechanical work rate of the fluid on the body, thus increasing the energy of the body. Negative power gives us an idea of how and when the body could potentially reduce energetic effort as it moves supported by the fluid. Fig. 9 shows that the majority of the mechanical power flows from the body to the fluid, and that negative power (fluid→body; Fig. 9B,D) is an order of magnitude smaller than positive power (Fig. 9A,C). Moreover, we observed three consistent characteristic peaks from the stacked power curves as shown in Fig. 9. The first peak appeared as positive power (Fig. 9A,C) in phase 1 as part of the S-shaped body bending. We found differences in the peak magnitude and the corresponding area below the power curves that varied with turn magnitude. This indicates that larger amounts of work are done by the body for turns with larger turn angles. For all types of turning maneuvers, we found that the first positive power peak contained contributions from all body segments along the anteroposterior axis, indicated by the full spectrum of colors from green to red in Fig. 9. The second power peak during turning appeared as negative power with main contributions from the anterior and middle part of the body. This peak in power indicates that both of these body regions move in supported by the hydrodynamic forces from the fluid at the end of phase 1 in the transition from the S-shaped to the J-shaped body posture. The next peak again occurred in the form of positive power and largely incorporated the posterior part of the body. This peak is related to the push of the tail in the direction of the final heading direction in phase 2 and was smaller than the first positive power peak, indicating that a smaller amount of mechanical work was done in this phase. In the final moments of the turning maneuvers, we found small peaks in negative power for some of the turns as shown in Fig. 9. However, these peaks were much less consistent and not observed in all turns.

We determined the cost of a turning maneuver from the mechanical work computed by integration of the power at each segment over time. Our measurements revealed negative work done by the water across all the analyzed recordings. As a proportion of the total work, we found that negative work ranged from 10.8% to 20.1%, with a median of 13.3% as shown in Fig. 10B. In addition, a negative linear correlation between maximum COM speed and the fraction of negative work was found (P<0.001, R²=0.66; Fig. 10A), indicating that slower turns had larger portions of negative work. Analyzing the contributions of the different segments along the anteroposterior body axis allowed us to locate the regions where negative work was done during the turning maneuvers. For this purpose, we computed the work at each segment for each turn by normalizing by the total work done in the respective turn in order to compare different turns. Then, we visualized the corresponding distribution across all turns with box plots (Fig. 10C); negative work showed the highest values in the anterior and mid-body, confirming the observations related to Fig. 9. We also analyzed the contributions of different body parts with regard to net mechanical work (Fig. 10D). Here, we found positive linear correlations between turning angle and net mechanical work in the anterior part (P<0.001, R²=0.84) and in the middle part of the body (P<0.001, R²=0.67). Interestingly, we could not find a corresponding significant linear correlation for the posterior part (P=0.11, R²=0.15), which indicates that the mechanical work done by the tail does not depend on the turning angle. Finally, in order to compare the cost of a zebrafish turn with both steady swimming behavior and other fish species, we computed the mass-specific TC for each recorded maneuver. We present the results of these computations in Table S3. The TC varied over a range from 5.13 to 71.35 mJ kg⁻¹ s⁻¹. We found a negative linear correlation between turning duration and TC (P<0.001, R²=0.49), and a positive correlation between turning angle and TC (P<0.005, R²=0.39), suggesting that costs increase for faster turns and larger turning angles.

Our results in this study were obtained with a modified PIV technique using algae as particles (following Gemmell et al., 2014) and infrared light, an approach that confers several benefits: we avoided intense laser illumination of fish and consequent behavioral modification, we were able to use small near-neutrally buoyant natural algal particles, and we avoided the creation of shadows from directional laser illumination. Furthermore, we used two high-speed cameras to ensure that we only analyzed sequences in which fish were precisely located in the camera plane of focus. We noted that larger angle turns in particular frequently displayed significant out-of-plane motion, which can adversely affect pressure calculations, and we thus restricted our analysis only to zebrafish turns that occurred in a horizontal plane as verified by a second lateral view camera, even though this reduced our sample size for larger amplitude turns. Analysis of pressure profiles around the zebrafish body during turns showed that even slight (∼20 deg) misalignment of the body angle with respect to the horizontal plane of camera focus substantially changed the estimated pressure profiles due to out-of-plane particle motion and calculated velocity profiles that do not reflect full body motion. We thus recommend that future studies include some method of ensuring that fish are centered in the focal plane without reducing water depth, which itself introduces surface and boundary effects. We found that the quality of the captured particle movements using algae in the fluid around turning zebrafish was comparable to or even exceeded traditional PIV approaches which use a laser, as evidenced by the high-resolution velocity vector fields shown in Fig. 6. Furthermore, our velocity fields are in general agreement with previous analyses on zebrafish exhibiting burst-and-coast behavior. For example, Müller et al. (2000) found the same overall vortex patterns across different scales (larvae and adults) that we identified in our analysis for turning zebrafish.

We hypothesized, based on the kinematics of the burst-and-coast swimming mode and the pattern of body curvature during the active burst phase, that certain body regions in turning zebrafish could perform negative work on the water. Specifically, we hypothesized that negative mechanical work occurs in posterior regions of the body during early phases of the turn near the time of maximal body curvature. Using our workflow (Fig. 1) to quantify the interaction between fish body motion and surface pressure patterns, we were able to quantify patterns of force and torque production along the body (Fig. 6) and then calculate mechanical power and work done along the body (Figs 9 and 10C). Our analysis supports the overall hypothesis that periods of negative work exist throughout the turn, but indicates that it is the anterior and middle body segments that experience the greatest overall proportion of negative work. Examination of the temporal pattern of work production in conjunction with patterns of body bending (Fig. 9), confirmed that these body regions consistently (cf. power curves of all analyzed turns in figshare: https://doi.org/10.6084/m9.figshare.c.4986179) experience times when the fluid is performing work on the body. The corresponding negative power peak consistently occurred subsequent to the first positive power peak. For example, from 0.04 s until 0.1 s in the turn shown in Fig. 9D, both the anterior and middle body regions were experiencing net negative work and the fluid was causing body deformation and displacement of these body segments. Although Fig. 9B,D indicates that the very last posterior segment (tail tip) is also contributing to negative work, this observation was not consistent over all the analyzed turns.

In our analysis, we focused on the interaction between a deforming fish body and the surrounding fluid environment, and quantified this interaction by means of mechanical power and work. This analysis differs from considerations of mechanical work produced by muscular forces that induce body deformations (Müller and van Leeuwen, 2006; van Leeuwen, 1995). Whereas our analysis provides a quantification of the external forces acting on a swimming body, mechanical work related to muscles involves the internal forces acting inside the body. Nonetheless, external and internal forces are interlinked as they emerge through the interaction of the body and fluid. This allowed us to relate metabolic rate to our mechanical work estimates based on external forces, as we show below. To get a more detailed view on the muscular dynamics during the turning maneuver would require measurement of muscle activity along the body. Given the small size of zebrafish, it is not possible to directly record electrical activity in body musculature from multiple locations along both sides of the body during free-swimming maneuvers. We can, however, formulate a hypothesis for a possible pattern of mechanical work by muscular forces during a turn.

Previous computational and kinematic studies of unsteady fish locomotion have shown that during rapid maneuvers, the posterior body region in fish is bent in a direction opposite to the direction of movement of that body area (Borazjani et al., 2012; Jayne and Lauder, 1993; Tytell and Lauder, 2008). This strongly suggests that the forces in this body region that arise from fluid loading are greater than those that arise from muscular activity in posterior body myotomes, which have been shown to be active during this phase of the maneuver (Jayne and Lauder, 1993). Muscle activity in posterior myotomes on the side of the body that is undergoing lengthening is reflective of net negative work at this time and in this particular body region. This contrasts with the prevailing view established by a large number of studies of fish steadily swimming (Coughlin, 2000; Coughlin and Rome, 1996; Ellerby et al., 2001; Jayne and Lauder, 1995; Johnson et al., 1994; Rome et al., 1993): measurements of myotomal muscle activity, in conjunction with body bending kinematics or sonomicrometry (to determine muscle strain), clearly show that locomotion is powered by positive work in the posterior body region. Thus, unsteady burst-and-coast maneuvers in zebrafish show a rather different pattern of energy exchange between internal muscular dynamics and external fluid forces than has been shown previously for steady swimming.

Computational fluid dynamic analyses can provide an informative comparison to the experimental results reported here. Borazjani et al. (2012) used computational fluid dynamics to study the escape response of bluegill sunfish. Although that study was on the rapid C-start behavior, it nonetheless provides an interesting comparison and contrast to our results on zebrafish maneuvers. Fig. 11 of Borazjani et al. (2012) illustrates surface pressure patterns on the body throughout the maneuver, and their computations show similar patterns of body surface pressure to those estimated from our experimental approach. Body surface pressure calculated for bluegill shortly after the start of the initial body bend (fig. 11B in Borazjani et al., 2012) showed a large negative pressure region on the inner curved tail surface, very similar to the pattern that we observed in zebrafish (Fig. 6B2,D2). As the maneuver proceeds and maximal body curvature is achieved, a strongly positive pressure is visible on the outer curved body surface, with strong negative pressure on the inner mid-body region (compare fig. 11C of Borazjani et al., 2012, with Fig. 6D3).

Li et al. (2012, 2014) studied flow patterns produced by maneuvering zebrafish larvae, and although the anatomy and muscular development of larval zebrafish differ substantially from those of mature individuals (McHenry and Lauder, 2006), instructive comparisons are still possible between their computational and our experimental results. Pressure fields computed for a larval zebrafish rapid maneuver similar to a C-start escape response (fig. 19 of Li et al., 2012: compare 50 ms and 60 ms times) show negative pressure on the inner curved tail surface, and positive pressure on the outer curved surface that propagates toward the tail during this part of the maneuver. Experimental data from our zebrafish maneuvers were generally similar, with positive pressure forming on the outer curved body surface as the body bends and the tail begins its contralateral stroke, but we did not observe the anterior movement of this high-pressure region and it remained localized at the posterior body (Fig. 6D3,D4). These differences could be due to differences in body flexibility and kinematics between larval and more anatomically mature zebrafish.

Our initial collection of 80 turning maneuvers showed that 90% of the turning angles were below 70 deg, which we classified as small-angle turns. Within a selection of the collected recordings, we evaluated and compared detailed kinematics and dynamics of angle turns above 70 deg against these more frequently occurring smaller turning maneuvers. The results can be summarized as follows: we found overall similar flow and pressure patterns for all turns, with large-angle turns differing mostly in displaying higher pressure magnitude. Additionally, we found differences in the force dynamics along the body, which show a traveling wave pattern for small-angle turns and a burst pattern for large-angle turns. This indicates that the underlying motor patterns could be different for the two cases, a suggestion that could be evaluated by recording patterns of muscle electrical activity in zebrafish during turning. The smaller angle turns resemble half-steady swimming bouts, whereas the larger angle turns show more similarities to C-start escape behavior.

In our analysis, we only found correlations of maximum bending and bending speed with respect to turning angle in the anterior and middle part of the body. The corresponding characteristics of the caudal region did not correlate with turning angle. This was especially noticeable in a few cases of similar turning angles of about 19 deg, which showed large variation in maximum bending and bending speed in the posterior body region. This variability is linked to turning duration, where faster maneuvers often showed larger maximum bending angle and speed. From these results we hypothesize that the anterior and middle parts of the body control the turning angle of these maneuvers.

Our analysis of body surface pressure and the subsequent processing steps used to generate work and power estimates allowed us to estimate the ‘cost of turning’ (TC). In order to put these costs in the context of energy consumed during swimming maneuvers in fish, we compared our results with metabolic rates measured in previous studies.

Plaut and Gordon (1994) measured oxygen consumption in swimming zebrafish at different controlled speeds and a water temperature of 28.0±0.2°C. In many of their experiments, they found no significant difference between the resting and active metabolic rate, which indicates that they estimated the net metabolic rate related to swimming as close to zero. In the experiments with significantly increased active metabolic rate compared with resting metabolic rate, net metabolic rate was found to be 100±42 mJ kg⁻¹ s⁻¹ (mean±s.d., N=4) at a swimming speed of 1 BL s⁻¹ (or 0.034 m s⁻¹), computed with an oxy-calorific equivalent of 18.89 J ml⁻¹ O₂ (Elliott and Davison, 1975; van Ginneken et al., 2005).

Our estimated TC across all experiments (30.53±19.52 mJ kg⁻¹ s⁻¹, mean±s.d., N=18) comprise 25.2% and 30.5% of the net metabolic rate from Lucas and Priede (1992) and Plaut and Gordon (1994), respectively. Previous work has reported the conversion of chemical free energy into mechanical work in muscles at 12–23% in isolated fish fast muscle fibers (Moon et al., 1991) and in isolated muscles of other animals (Smith et al., 2005), from 14% in mice up to 41% in toads. Our mechanical work estimates from body pressure profiles are in general accordance with previous metabolic rate measurements and the low mechanical efficiency of skeletal muscles.

We consider a turning maneuver to be efficient when a change in direction is achieved with a low energy effort, acknowledging that estimates of organismal efficiency are complicated and difficult to compare among species and behaviors. To our knowledge, there has been only one recent study (Schakmann et al., 2020) that has quantified the cost of turning in another fish species, namely in the goldring surgeonfish (Ctenochaetus strigosus). Although it is mostly a pectoral-fin swimmer, Schakmann et al. (2020) provide a series of experiments in which repeated turning by means of body bending (but with additional pectoral fin movements) is induced in a bidirectional oscillating flow. By taking the metabolic rate that they found in unidirectional flows and looking at the difference with the rate at bidirectional flows (table 3 of Schakmann et al., 2020), we obtain an estimate for the metabolic rate of turning for the different flow conditions: 75.9 mJ kg⁻¹ s⁻¹ at an oscillating flow frequency of 0.1 Hz and average speed of 1.5 BL s⁻¹; 385 mJ kg⁻¹ s⁻¹  at 0.1 Hz and 3.0 BL s⁻¹; and 559 mJ kg⁻¹ s⁻¹ at 0.3 Hz and 3.0 BL s⁻¹. Based on these measurements, the turning costs at the lowest oscillating flow frequency (0.1 Hz) and flow speed (1.5 BL s⁻¹) are lower in goldring surgeonfish than in our zebrafish turning experiments. However, they are larger for more unsteady and faster flows. Given that the zebrafish turns were performed in still water, we hypothesize that zebrafish might have a higher cost of turning. A possible explanation could be linked to a beneficial use of pectoral fins to support body bending during turns. Nonetheless, this remains to be tested as these experiments were performed under different conditions. Future experiments applying our PIV and pressure computation approach in an unsteady flow as proposed in Schakmann et al. (2020) could further clarify the underlying mechanistic dynamics of turning in these conditions.

In addition to the comparison with data from other fishes, our analysis in this study revealed periods of negative work, which lowers the net mechanical work. This can be considered an energy-efficient feature of the turning maneuver, because as part of the second power peak during a turning maneuver, the anterior and mid-body movements are supported by the fluid. The resulting dynamics are a direct consequence of the first power peak, where fluid is set in motion. We found in our analysis that negative work represents between 10% and 20% of the total mechanical work, which raises the question why positive and negative work are distributed this way. There are many ways to spend the same amount of mechanical work, e.g. by increasing both positive and negative work accordingly, the same net mechanical work could be achieved. It remains unclear whether or why the ratio of these two quantities is constrained. We propose the following hypotheses that could explain our observations: (1) the ratio of positive and negative work is not related to energy efficiency but to the outcome of the natural dynamics of turning in a fluid at the observed Re numbers; (2) the observed balance between positive and negative work results from zebrafish aiming to maximize the forward velocity (and coast distance) at the end of the turn (similar to Gazzola et al., 2012); (3) the ratio of positive to negative work represents an optimum in order to perform an energy-efficient turn; and (4) the proportions of positive and negative work contribute to energy optimality within a burst-and-coast sequence.

The computation of the forces and the mechanical power/work as presented in this study relies on a few assumptions that consequently lead to a number of limitations. The forces computed at the interface between body and fluid are based on the obtained pressure field around the moving fish. Here, we assumed that the corresponding forces are dominated by pressure forces that act normal to the body surface, and that shear forces related to viscous effects are small in comparison. This assumption is supported by computational fluid dynamic simulations of a larval zebrafish C-start (Re=550) carried out by Li et al. (2012). The range of Re numbers in our study was similar, but reaches limits for slow and small-angle turns where smaller Re numbers occur and shear forces can be more significant. Li et al. (2012) show that the (much smaller) shear stresses appear most strongly at the head, which indicates that our approach most likely underestimates drag forces at the anterior part during turning maneuvers.

In order to extrapolate from pressure to force we multiplied pressure with corresponding body surface slice area. We approximated the varying pressure fields along the depth of the fish body by a constant pressure field measured at the mid-section of the body. Lucas et al. (2017, 2020) showed in their analysis of a tail-shaped flapping foil and fish swimming that the resulting estimations of force and torque based on pressure at the mid-section of the foil were in good agreement with experimental measurements from force and torque transducers, and from computational analyses of fish swimming. This and our comparison with previous measurements of metabolic rate in zebrafish supports our general approach.

Finally, the approach presented in this study considers exclusively locomotor patterns resulting from body undulations. Experiments with extended pectoral fins were excluded from the PIV and pressure analysis as mapping from the measured pressures around a fin in addition to quantifying pressures around the bending body would require a detailed analysis of the 3D kinematics of the fins (see Danos and Lauder, 2007, for a discussion of zebrafish fin dynamics). Computational fluid dynamic approaches may be more suited to estimating body pressure patterns in such multi-fin complex 3D locomotor behaviors (Liu et al., 2017).

Many thanks to Dave Matthews for help with the experimental setup and members of the Lauder Lab for their assistance with this project.