Fish larvae tackle the complex fluid-structure interactions of undulatory swimming with simple actuation

Most fish swim with body undulations that result from fluid-structure interactions between the fish’s internal tissues and the surrounding water. As just-hatched larvae can swim effectively without a fully-developed brain, we hypothesise that fish larvae tackle the underlying complex physics with simple actuation patterns. To address this hypothesis, we developed a dedicated experimental-numerical approach to calculate the lateral bending moment distributions, which represent the system’s net actuation. The bending moment varies over time and along the fish’s central axis due to muscle actions, passive tissues, inertia, and fluid dynamics. Our 3D analysis of a large dataset of swimming events of larvae from 3 to 12 days after fertilisation shows that these bending moment patterns are not only relatively simple but also strikingly similar throughout early development, and from fast starts to periodic swimming. This suggests also similar muscle activation patterns, allowing fish larvae to produce swimming movements relatively simply, yet effectively, while restructuring their neuromuscular control system.


25
Swimming is a vital component of the fitness of a fish because fish swim to search for food, 26 hunt prey, escape from predators, migrate and disperse, and manoeuvre through complex 27 environments. Many fish species swim by performing body undulations that result from an 28 interaction between body tissues and the surrounding water [1,2]. Understanding these 29 complex fluid-structure interactions is crucial to gain insight into the mechanics and control 30 of fish swimming [3]. 31 To analyse the fluid-structure interactions during swimming, we need to understand the 32 external fluid mechanics (water), the internal solid mechanics (skin, muscle, skeleton), and 33 their coupling. During swimming, the body of the fish moves through the water, inducing a 34 flow around it [4,5]. The resulting fluid dynamic forces interact with the body tissues via the 35 skin, resulting in a change in deformation. This deformation will change the motion of the 36 surface of the body, which influences the fluid dynamic forces, thus forming a loop of tight 37 coupling between the fluid mechanics and the internal solid mechanics [2,6,7]. This complex 38 two-way fluid-structure interaction creates the typical travelling wave pattern observed in 39 swimming fish [8,9]. 40 The complexity of the physics would suggest that fish need a sophisticated control system 41 to produce swimming motions reliably. Zebrafish (Danio rerio, Hamilton 1822) larvae, the 42 subject of this study, would seem to contradict this: they can swim immediately after hatching, 43 at considerable speed and tail beat frequency [3, 10,11]. Furthermore, the spinal cord can 44 initiate swimming motions even when severed from the brain [12]. This suggests that a 45 relatively simple system can produce reliable undulatory swimming, despite the non-linear 46 governing physics. Throughout the first days of development, the larvae refine their control 47 of swimming [10,13,14] and improve swimming performance [10,11]. These improvements 48 raise the question: do fish larvae produce swimming differently across early development, and 49 at different swimming speeds and accelerations? 50 To answer this question, we need insight into the internal mechanics of the axial muscles 51 and passive tissues that actuate the motion. Muscle activation patterns can be measured 52 directly with electromyography, where electrodes are inserted in the muscles to measure 53 potential differences [15][16][17]. However, this technique may incur considerable changes in 54 swimming behaviour. Especially for fish larvae, it requires them to be paralysed [18] or fixed 55 hence its speed. The resultant power (Fig 2H), defined as the sum of the fluid and kinetic 115 power, is dominated by the fluid power. 116

Swimming effort and vigour 117
We reconstructed 3D kinematics from 113 video sequences of fast-start responses followed by 118 swimming, calculated flow fields throughout the sequence with CFD, and fitted distributions 119 of internal forces and moments. These swimming sequences hardly contain periodic 120 swimming. To analyse the data despite its aperiodicity, we subdivided it into half-beats based 121 on zero-crossings of the bending moment in the mid-point along the centreline (Fig 1F,G). For 122 each of these 285 half-beats, we calculated the period length, mean speed, mean acceleration 123 of the centre of mass to the next half-beat, and peak (95 th percentile) bending moment. 124 To reduce the number of parameters for the analysis, we identified combinations of 125 parameters with high explanatory capacity. To control swimming, the fish has two main 126 parameters to change the bending moment (see section below): its amplitude and the duration 127 of each half-beat. We define the swimming effort as = peak half −1 -higher bending moments 128 and shorter periods increase . We fitted a generalised linear model ( We expect the net propulsive force to scale with the mass of the fish, its acceleration, and 134 its squared speed (from the dynamic pressure). Based on this, we define swimming vigour as 135 = ( 2 + ), where m is body mass, v is swimming speed, and a is mean acceleration (i.e. 136 change of speed to the next half-beat per unit of time). The coefficient is calculated with an 137 optimisation algorithm that minimised the sum of squared errors of a linear fit of vigour to 138 effort with total least squares. The fitted value of 517.7 m -1 results in a clear trend of vigour as 139 a function of effort (Fig 3B; generalised linear model fit with gamma distribution and log link 140 function: P < 0.0001), collapsing the broad clouds of speed and acceleration (Fig 3C,D). 141

Bending moment distributions are similar across swimming vigour and development 142
To assess how bending moment patterns differ across vigour and development (indicated by 143 body length [27]), we compared bending moment patterns normalised by their amplitude. We 144 normalised the bending moment distribution of each half-beat by dividing by the peak 7 a curve through the effort landscape, a function of peak bending moment and half-beat 177 duration (Fig 5). In general, high peak bending moments are only produced for tail beats of 178 short duration. As the duration decreases (i.e. frequency increases), the bending moment 179 amplitude decreases. Higher efforts generally lead to higher speeds (Fig 5B), unless the larva 180 is accelerating strongly. Strong accelerations are mostly found with slow-swimming larvae 181 using short half-beat durations and high peak bending moments (Fig 5A). For high-effort tail 182 beats, the larvae are generally either swimming fast, or accelerating: high-effort, low-speed tail 183 beats show high accelerations, while high-effort, low acceleration tail beats show high speeds. 184 Swimming vigour tends to increase with increasing effort (Fig 5C, see also Fig 3B). 185

186
In this study, we analysed the actuation used during swimming of developing zebrafish larvae 187 by calculating bending moment patterns with inverse dynamics. We found that larvae use 188 similar bending moment patterns across development. They adjust their swimming vigour, a 189 combination of speed and acceleration, by changing the peak bending moment and tail-beat 190 duration. At higher speeds and accelerations, the larvae produce the required fluid-dynamic 191 forces by increasing bending moment amplitude and/or decreasing tail-beat duration. 192 Previous inverse-dynamics approaches for the internal mechanics of swimming used 193 simplified models for both the fluid and the structure. The structure was modelled with linear 194 bending theory, assuming small deformations of the centreline [22,24]. The effects of a large-195 amplitude correction to these was expected to be small for adult fish that swim periodically 196 In our analysis of swimming sequences across development, we do not assume periodicity. 210 Periodic motion is a common assumption in the analysis of fish swimming [11,32]. For 211 zebrafish, cyclic swimming occurs most often after a fast start, and rarely spontaneously 212 [10,33], and generally only for a few tail beats. To analyse aperiodic motion, we subdivided 213 each swimming sequence in tail beats based on zero-crossings of the bending moment in the 214 middle of the body. During aperiodic swimming, the mean speed varies between successive 215 tail beats. For this reason, we define a parameter "swimming vigour", that combines the effects 216 of acceleration and swimming speed as = ( 2 + ) . This approach for analysing 217 aperiodic swimming could be of general use in swimming research. The subdivision in half-218 beats can also be done with quantities other than the bending moment, for example body 219 curvature. This enables similar analyses of aperiodic swimming from pure kinematics without 220 inverse dynamics. 221 We define the swimming effort exerted by the fish based on the amplitude of the bending 222 moment and the duration of the tail-beat, as = peak half −1 . This quantity correlates with 223 resultant power, indicating that it is indeed an indicator of amount of effort exerted by the 224 larva ( Fig 3A). The speed and acceleration fall on broad clouds as a function of effort (Fig  225   3C,D), since the required power depends on their combination, rather than their individual 226 values. The effort-speed landscape ( Fig 3C) shows a two-pronged distribution, one branch 227 showing high effort but low speed, and the other, broader branch showing increased effort 228 with speed. This distribution is mainly explained by the acceleration, showing high values in 229 the lower branch-fish only accelerate strongly from low speeds and use high effort to do so. 230 This is reflected in the effort-acceleration landscape (Fig 3D), low (including negative) 231 acceleration are found at high speeds, and vice versa. 232 When speed and acceleration are combined into the swimming vigour, these clouds 233 collapse more closely to a trend line (Fig 3B). Variation in this trend may be partly caused by 234 turning behaviour and contributions of the pectoral fins. We can estimate the relative 235 contribution of the speed and acceleration to the swimming vigour, giving an indication of 236 their relative cost. If we assume force production to maintain speed and to accelerate are 237 equally costly, we can estimate a drag coefficient from the coefficient in the vigour equation as 238 ,estimated = ( 1 2 0.061, which is considerably lower than the value of 0.26 calculated from a previous CFD study 240 on larval zebrafish [34]. This means our equal-cost assumption does not hold: the contribution 241 of the speed term is relatively low compared to the acceleration term. Since swimming vigour 242 correlates with swimming effort, this indicates that acceleration is more costly to achieve 243 compared to maintaining swimming speed-the larvae need to invest more effort to produce 244 force to accelerate than to swim steadily. 245 Most of the resultant power produced by the fish is used to increase the energy in the fluid, 246 rather than the kinetic energy of the body (Fig 2H). The energy spent on the water is likely lost 247 bending moment was found to be higher than the overall bending moment, while the wave 258 speed was found to be lower. However, the overall dynamics look reasonably similar. If we 259 assume similar distributions of passive tissues inside the fish across the considered 260 developmental stages [27], similar total bending moment patterns will require a similar muscle 261 contribution. Furthermore, the difference in amplitude between similar bending moment 262 patterns must originate from the muscle moment, since it is the only net source of power in 263 the system -the work done by the fluid and passive tissues indirectly comes from the muscles. 264 We found that the bending moments follow a similar pattern across development and 265 swimming vigour (i.e. speed and acceleration). The only significant coefficient in the linear 266 models is the phase of the centre of volume of the bending moment patterns as a function of 267 swimming vigour (Fig 4F), but the effect is limited. More vigorously swimming fish generate 268 the peak bending moment slightly earlier in the half-beat. The mean pattern looks qualitatively 269 similar to earlier calculations done for adult fish [24]. It is a single-peaked distribution, with 270 the maximum around the bulk of the muscle (Fig 2E, Fig 4A), and a fast-travelling wave 271 character. Muscle electromyograms (EMG) done on paralysed zebrafish also looked similar to 272 adult activation patterns [18]. This suggests that this simple pattern of bending moments is 273 common to fish across species and developmental stage. Even though fish larvae swim in the 274 intermediate regime [3], and adult fish often swim in the inertial regime [39], the differences 275 in fluid dynamics do not seem to require fundamentally different bending moment patterns. 276 Since the bending moments look similar along the body and over the phase for each half-277 beat, the two parameters left for the larvae to adjust for each half-beat are its duration and the 278 amplitude of the bending moment. Young larvae use a relatively narrow range of amplitudes 279 ( Fig 4H) and durations (Fig 4G), which broadens as the fish develop. Older larvae are able to 280 generate higher peak bending moments, likely correlated to development of their muscle 281 system [40]. Furthermore, older larvae use a broader range of tail beat durations than young 282 larvae, suggesting that older larvae have more freedom to control their swimming vigour.

296
An in-depth mathematical treatment of the methods is given in the Supplemental Information. 297 We made high-speed video recordings of fast starts of three separate batches of 50 zebrafish 299 larvae from 3-12 days post fertilisation (dpf). The camera setup was identical to the setup 300 described in Voesenek et al. [26], with three synchronised high-speed video cameras, recording 301 free-swimming larvae at 2000 frames per second. To reconstruct the swimming kinematics 302 from the recorded high-speed video, we used in-house developed automated 3D tracking 303 software [26] in MATLAB (R2013a, The Mathworks). For every time point in a multi-camera 304 video sequence, the software calculates the best fit for the larva's 3D position, orientation and 305 body curvature to the video frames. These parameters are then used to calculate the position 306 of the larva's central axis and the motion of its outer surface (Fig 1A,B). 307

Subdividing motion 308
We calculated phase-averaged quantities for an individual swimming sequence to look in 309 detail at the generated bending moments and powers. We determined whether a (subset of a) 310 sequence is periodic with a similar approach to Van Leeuwen et al. [11]. For every possible 311 subset of a swimming sequence, we calculated the sum of absolute difference with a time-312 shifted version of the curvature, similar to an autocorrelation. We then calculated extrema in 313 this function -if extrema are detected, their maximum value determines the "periodicity" of 314 the sequence. We then selected the longest possible subsequence that has a periodicity value 315 higher than a threshold of 35 -this is a swimming sequence for a 3 dpf fish. We divided this 316 sequence in half-phases based on peaks in the body angle [11,26]. The curvature, bending 317 moment, fluid power, kinetic power, and resultant power were then phase-averaged based on 318 these subdivisions. 319 Most of the swimming of larval zebrafish is aperiodic, but there is an alternating pattern in 320 the bending moments. For this reason, we analysed swimming per half-beat, based on the 321 bending moment. We found the zero crossings of the bending moment at 0.5ℓ. Since some of 322 these points are related to noise, we evaluated every possible permutation of zero crossings 323 per sequence on several criteria with a custom MATLAB (R2018b, The Mathworks) program. 324 We eliminated zero-crossings with neighbouring sections with an amplitude of less than 5% 325 of the peak half-beat amplitude in the sequence, as they are most probably noise. We required 326 more than three zero crossings to have at least two half-beats to be able to calculate a mean 327 acceleration. Extreme values in each half-beat should alternate direction to eliminate noisy 328 zero crossings: the larva beats its tail left and right, so therefore bending moment must 329 alternate. Finally, we eliminated half-beats with a duration shorter than 2.5 ms (equivalent to 330 200 Hz tail-beat frequency) -the maximum tail-beat frequency observed for zebrafish larvae 331 is 95 Hz [11]. From all permutations that met the criteria, we selected the permutation with the 332 smallest standard deviation in half-period length across the sequence. This left the longest 333 possible, least noisy sequence of half-beats for every swimming bout. 334 Out of 113 swimming sequences, we selected 398 half-beats with this procedure. For each 335 of these half-beats, we calculated the duration, mean speed, and peak bending moment. We 336 determined the mean acceleration by calculating the difference in mean speed between the 337 following and current half-beat. Since we could not calculate mean acceleration for the last 338 half-beat in each sequence, 285 half-beats remained for which we computed all quantities. 339

Calculating fluid-force distributions 340
To calculate fluid-force distributions, we used the immersed boundary method Navier-Stokes 341 solver IBAMR [46]. We converted the tracked video data into a three-dimensional point cloud 342 model in the fluid solver. We exclusively used swimming sequences where the larvae start 343 from rest in quiescent water, so we do not need to consider history in the wake. The solver 344 time step was much smaller than the time step between video frames, so we interpolated the 345 reconstructed kinematics with a quintic spline [47]. Using this interpolated state, we updated 346 the location of the point cloud representing the surface of the fish. This resulted in a three-347 dimensional velocity and pressure field at every point in time ( Fig 1C). To verify the accuracy 348 of the method, we compared reconstructed bending moments from IBAMR to an 349 experimentally validated CFD solver [31,34,48], showing only small differences, see the 350 Supplemental Information. 351 We extracted force distributions by interpolating the pressure and velocity gradient tensor 352 components to the centre of each face of an offset triangulated representation of the fish 353 surface. We then integrated these values into contributions to the pressure force and the shear 354 force at every face of the surface (Fig 1D). By further integration, we calculated the force at 355 every point along the centreline in a coordinate system attached to the larva's head ( Fig 1E). 356

Calculating bending moments 357
To calculate bending moments, we represented the fish by its central axis only. Effects of 358 muscles, spine, and other tissues were combined for every transversal slice along this axis. two-dimensional space. We derived the equations of motion for this beam (see Supplemental 361 Information) in an accelerating and rotating coordinate system attached to the fish's head [49]. These data were also used to calculate the resultant power from the fluid-dynamic forces 383 and the changes in kinetic energy. Note that we cannot calculate a meaningful value for the 384 internal power, because we do not separate passive and active effects. The net power when 385 combining these effects does not correspond to the actual power consumption of the fish, as 386 the passive and active components might partly compensate each other. 387

Calculating muscle cross-sectional area 388
We performed micro-computed-tomography (μCT) images of a 3 dpf zebrafish larva at the scans with a resolution of 650×650×650 nm per voxel. From these data, a centreline was 392 extracted by finding the centre of area of each slice, segmented with simple thresholding. 393 Finally, the muscle area was manually digitised in 51 planes, for which the image data was 394 interpolated in a plane perpendicular to the centreline. 395            In this study, we calculated internal forces and moments for swimming zebrafish larvae. The   a suitable coordinate transformation. In summary, we model the fish as a beam undergoing 21 large bending deformations in two dimensions. 22 We describe the deformation of the beam with the displacement of each infinitesimal beam 23 element with respect to the reference configuration (Fig. 1A). It is defined as a function ξ(s, t) = where (x 0 (s), y 0 (s)) is the reference configuration of the beam. We define the reference config- This displacement results in a local deformation angle θ(s, t) for each beam element (Fig. 1A).

29
It can be calculated from the displacements with: To derive the equations of motion, we consider a free body diagram of an infinitesimal 32 beam element of length ds (Fig. 1B). The force balance in xand y-direction for a beam element 33 ρ fish A(s) ds ∂ 2 ξ ∂t 2 (s, t) = N (s, t) cos θ(s, t) − N (s + ds, t) cos θ(s + ds, t) − Q(s, t) sin θ(s, t) + Q(s + ds, t) sin θ(s + ds, t) where ρ fish is the beam density, A is the beam cross-sectional area, N is the internal normal , and dropping the explicit f (s, t) notation, yields the 46 equations of motion: 50 We reconstructed the motion of the fish from video in three-dimensional space, but described 51 the equations of motion in a two-dimensional plane. However, in the video-tracking method, 52 we assumed that the fish deforms in a single plane. Hence, we can create a coordinate system 53 aligned to this plane and obtain the equations in two dimensions only. We define this head 54 reference frame as fixed to the snout of the fish and rotating along with the stiff head region in 55 the deformation plane. Any pointx in world coordinates (denoted with a circumflex) can be 56 expressed in the fish coordinate system at time t as:

Fictitious forces
where R(t) is the time-dependent rotation matrix expressing the orientation of the snout and 58x snout is the position of the snout.

59
When we transform the motion to the non-inertial fish reference frame, additional equation 60 terms accounting for the effect of the translation and rotation of the frame must be considered.

61
These additional acceleration terms for any point x in the rotating reference frame are given 62 by 2 : whereâ head is the acceleration of the origin,ω is the angular velocity of the rotating coordinate 64 system, v is the velocity of the considered point, andα is the angular acceleration of the rotat-2 Calculating fluid forces from kinematics 73 The equations of motion include an external force distribution, produced by the water on 74 the skin. Since this is exceedingly difficult to measure directly and non-invasively, we mod-75 elled the fluid dynamics numerically. We used two independent computational fluid dynam-76 ics (CFD) methods to calculate fluid-dynamic forces. We used an experimentally validated 77 method to validate the second method to assess its accuracy when calculating internal forces 78 and moments. 79 We performed computational fluid dynamics using a Navier-Stokes solver based on over- The surface of the fish was described as a cloud of Lagrangian points, moving over the 97 Eulerian fluid solution mesh. The motion of these points was prescribed based on quintic 98 spline interpolation 7 of the tracked kinematics, with a custom-developed add-on to IBAMR.

99
The density of the point cloud was chosen such that the mean distance between the points is 100 0.75× the smallest mesh level. This ensured that each cell inside the fish body will have at least 101 one point in it, and not much more.

102
skin of the fish with a custom Python 3 program. In this program, we interpolated 8 the pres-104 sure and velocity gradients to a triangulated surface slightly offset from the fish skin. These 105 were then used to calculate the local stress on each triangular face as: where p is the pressure, τ is the shear stress tensor, and n is the outward facing normal of the 108 face. Under the assumption of a Newtonian fluid, the shear stress tensor is defined as: where µ is the dynamic viscosity, and u, v, w are the velocity components in respectively x-, y-, 110 and z-direction. These surface stress distributions were then grouped into segments along the 111 fish to calculate the local net fluid force in the moving reference frame.

112
3 Reconstructing internal forces and moments with inverse dynamics 113 We reconstructed bending moments from the motion of the fish and its simulated external fluid 114 force distribution, an approach commonly called inverse dynamics. This section describes the 115 optimisation procedure we used to reconstruct the internal forces and moments.

116
In our inverse dynamics approach, we cannot separate the effects of the active and passive 117 tissues inside the fish: the internal forces and moments that we compute include the effects of 118 both. Considering this, the moment equation becomes: where ∂M * ∂s = ∂M ∂s − m muscle . This combined moment is what we reconstruct from the motion 120 with our optimisation procedure. From here onwards, we drop the asterisk notation and refer 121 to the combined active and passive internal bending moment as the 'bending moment'.   137 We calculated the resultant power on the fish from two source: the power exerted on the fluid, 138 and the changes in kinetic energy. Both quantities are computed in the inertial reference frame. 139 We calculated the power per unit length that the fish exerts on the water as:

Calculating resultant power
wherev is the velocity of the centreline in world coordinates. We negated the power since we 141 are considering the power that the fish exerts on the water, rather than the inverse.

142
The kinetic energy per unit length at any point in time was calculated as: The kinetic power per unit length is the time derivative of the kinetic energy: Note that we do not calculate powers based on the internal forces and moments. We do  The normal forces and moments can be calculated from the displacements using constitutive 158 equations. To generate the reference data, we assumed a Hookean material, resulting in the 159 following equations: where Y is the Young's modulus. The strain ε can be computed from the displacements as 162 follows: These relations complete the set of equations required to calculate the accelerations of each 164 point on the beam with equations 7 and 8. 165 We used the backward Euler method to integrate the beam accelerations to velocities, and 167 velocities to displacements. To calculate the velocities and displacements in x-direction at the 168 time step i, we used: and analogous expressions for η.
The scaled equations for temporal integration of ξ then become: and analogously for the other three equations related to integration of ξ and η.
The scaled moment balance becomes: Finally, the constitutive and kinematic relations become: 185 We integrated the scaled equations of motion with multivariate Newton-Raphson, an iterative 186 solution method for non-linear partial differential equations. Based on the solution at the pre- points in the beam. To calculate the value for the next iteration, we solved:

Solution method
where the subscripts i and i + 1 denote the current and next iteration, J F is the Jacobian of the 192 residuals, F is the vector of residuals, and ϕ is the vector of variables. 193 We calculated the Jacobian numerically by perturbing each variable with a fixed-step size 194 and then calculating one-sided finite-differences. We solved the system with a direct solver 9 .

209
We simulated a beam with the same geometry and loading, but in a time-dependent man- f y,fluid (s, t) = A y cos(2πf t), and the muscle moment by: The resulting motion is shown in Fig. 2B, demonstrating a swimming-like motion. We selected 222 a single cycle from the motion (Fig. 2C), after the motion had become reliably periodic after 223 11 cycles. The total internal bending moment, which is what we reconstruct, consisted of an 224 elastic contribution (Fig. 2D), and a contribution from the 'muscle' moment (Fig. 2E). 225 We reconstructed the total internal forces and moments of the reference 'swimmer' with from IBAMR, we also calculated the internal forces and moments to compare to the reference 245 from the validated solver.