Emergence of in-line swimming patterns in zebrafish pairs

Graphical Abstract Mathematical models promise new insights into the mechanisms underlying the emergence of collective behaviour in fish. Here, we establish a mathematical model to examine collective behaviour of zebrafish, a popular animal species in preclinical research. The model accounts for social and hydrodynamic interactions between individuals, along with the burst-and-coast swimming style of zebrafish. Each fish is described as a system of coupled stochastic differential equations, which govern the time evolution of their speed and turn rate. Model parameters are calibrated using experimental observations of zebrafish pairs swimming in a shallow water tank. The model successfully captures the main features of the collective response of the animals, by predicting their preference to swim in-line, with one fish leading and the other trailing. During in-line swimming, the animals share the same orientation and keep a distance from each other, owing to hydrodynamic repulsion. Hydrodynamic interaction is also responsible for an increase in the speed of the pair swimming in-line. By linearizing the equations of motion, we demonstrate local stability of in-line swimming to small perturbations for a wide range of model parameters. Mathematically backed results presented herein support the application of dynamical systems theory to unveil the inner workings of fish collective behaviour.


Introduction
The majority of fish species live part of their lives in a group; this affords several advantages, which include enhanced abilities to escape from predators, search for food and find correct migratory routes (Larsson, 2012). Common to life in a group is swimming together in a school, where fish synchronize their motion in a crystallized formation of tightly-packed individuals (Miller & Gerlai, 2011). The mechanisms underlying the formation of these schooling patterns and their potential hydrodynamic implications have long been the subject of intense debate, which is yet to be resolved (Ashraf et al., 2017;Partridge & Pitcher, 1979;Weihs, 1973).
While experimental research is the cornerstone against which new hypotheses must be tested, mathematical modelling offers a promising approach for detailing the inner workings of collective behaviour. Since the seminal work of Aoki (1982), a number of modelling approaches have been proposed to study collective behaviour (Vicsek & Zafeiris, 2012). Reaching beyond agent-based models implementing behavioural rules in discrete time, several authors have proposed modelling approaches based on stochastic differential equations, such as Butail et al. (2016), Calovi et al. (2014Calovi et al. ( , 2015 and Zienkiewicz et al. (2015a,b).
Particularly promising is the recent model by Filella et al. (2018), which considers, for the first time, the presence of hydrodynamic interactions between swimming fish. These interactions are modelled by associating each fish with a vortex dipole, which encapsulates both self-propulsion and the velocity field generated by a swimming fish. The velocity field generated by one fish induces advection and rotation on other animals in the vicinity, and vice versa, thereby providing a mechanism for bidirectional hydrodynamic interactions between fish. While paving the way toward the mathematical integration of social and hydrodynamic interactions, their study was not grounded in experimental observations and all of the claims were based on numerical simulations.
Here, we seek to address these limitations by tailoring the model of Filella et al. (2018) to real data on pairs of zebrafish, a species of choice in the study of collective behaviour in laboratory settings (Kalueff et al., 2013;Orger & de Polavieja, 2017). To achieve this goal, we include realistic turn rate dynamics that capture the burst-and-coast swimming style of zebrafish (Kalueff et al., 2013). We calibrate the model on real data from Zienkiewicz et al. (2015a) and demonstrate its predictive power on multiple measures of collective behaviour. Going beyond numerical simulations, we offer analytical insight into the model by elucidating the stability of in-line swimming. This swimming pattern has been experimentally and numerically studied in the technical literature on collective behaviour with multisensory cues, but its stability has never been analytically investigated (De Bie et al., 2020;Kato et al., 2004;Laan et al., 2017;Perna et al., 2014;Sadaat et al., 2021).
The rest of the paper is organized as follows: § 2 introduces the modelling framework; § 3 offers analytical insight into the stability of schooling patterns; § 4 discusses calibration of model parameters on experimental results; § 5 presents and discusses the results; and § 6 summarizes the main findings of the work and identifies avenues for future research.

Modelling Collective Behaviour of a Zebrafish Pair
We model zebrafish as self-propelled bodies swimming in an unbounded two-dimensional plane. At time t, fish f ( f ∈ {1, 2}) is identified by the position of its centroid r f (t) and swimming directionv f (t) with respect to a global reference frame (Figure 1a). Working with a Cartesian coordinate system with unit vectorsî andĵ for x and y, respectively, we write r Figure 1a). Toward the analytical treatment of the problem, we hypothesize that each fish has a constant selfpropulsion speed v f . When swimming in a background flow, the velocity of the animal in the global reference frame is the superposition of self-propulsion and advection. In particular, swimming in a pair causes the speed of fish f to dynamically change with respect to the fixed coordinate system owing to the flow created by the other individual in the pair. Specifically, the velocity of the focal fish f is given by (Filella et al., 2018) where a superimposed dot means a time derivative and U f (r f (t)) is the advection velocity experienced by fish f at time t owing to the presence of the other fish in the pair. As a first approximation to the fluid flow generated by the other fish in the pair, we employ the potential flow framework developed by Filella et al. (2018), Gazzola et al. (2016), Kanso & Tsang (2015) and Tchieu et al. (2012), wherein vorticity shedding owing to swimming is neglected and the fluid is modelled as incompressible, inviscid and irrotational. From an application perspective, the most restrictive of these assumptions are those of irrotationality and no vorticity shedding. Irrotationality will be violated in any sort of shear flow, which includes flow near solid boundaries, thereby challenging the use of the model to examine wall-hugging behaviours that could be associated with stress (Kalueff et al., 2013). Furthermore, neglecting vorticity shedding precludes the possibility to examine hydrodynamic interactions when the fish are in close proximity, interacting with each other's wake (Weihs, 1973).
Following the finite-dipole paradigm, originally proposed by Tchieu et al. (2012), fish f in the pair is represented as a dipole in the far field, wherein the vortices comprising the dipole are positioned orthogonal to the swimming direction at a distance r 0, f . As demonstrated in Figure 1 of Tchieu et al. (2012), this representation emulates the far-field streamlines around a self-propelled body in an inviscid fluid. At position r f (t), the other fish, denoted asf , creates the advective field (Filella et al., 2018) where (t) = |rˇf (t) − r f (t)| is the distance between the two fish. In (2), the unit vectorsêˇf (t) andêˇf (t) define a local polar coordinate system centred at fishf , in whichêˇf (t) points toward the other fish in the pair and forms a right-handed coordinate system withêˇf (t) andk =î ×ĵ. The viewing angleˇf ,f (t) is the angle between the swimming direction of fish j andêˇf (t) (Figure 1a). The distance r 0, f should be of the order of the amplitude of the tail beating and its relationship with the speed of the animal provides the circulation of each vortex in the dipole, f = 2πr 0, f v 0, f . For the analysis, it is more convenient to express the velocity field (2) in the global Cartesian coordinate system. From Figure 1a, we deducê where Δx(t) = x 1 (t) − x 2 (t) and Δy(t) = y 1 (t) − y 2 (t), such that (t) = Δx(t) 2 + Δy(t) 2 . Similarly, we can express the sines and cosines of the viewing angle in (2) in terms of individual orientations as where the plus (minus) sign is for f equal to 1 (2). By replacing (3) and (4) into (2), we establish To complete the model of the fish pair, we must prescribe governing equations for the orientation f (t), f ∈ {1, 2}. Different from the work of Filella et al. (2018), we employ a dynamic model that captures inherent delays in zebrafish to respond to visual or hydrodynamic input. Based on earlier models of social interaction, such as those by Calovi et al. (2014Calovi et al. ( , 2015, we propose Here, f is the mean reversion rate that defines the characteristic time scale for fish to respond to external stimuli, thereby condensing neurological processes in decision-making and physical delays in reaction to external stimuli into a single damping coefficient per unit moment of inertia (in a kinematic model like in the work of Filella et al. (2018), f → ∞); ★ (t) is the imposed turn rate from social interactions; f (t) is the imposed turn rate owing to hydrodynamic interactions; and N f (t) is added noise in the form of a zero-mean stationary stochastic process that captures non-deterministic forcing from the surroundings (called 'free will' by Filella et al., 2018). All these contributions are linearly superimposed in the present model. The modelling choice of including reversion on the turn rate dynamics while assuming constant speed is a first approximation, motivated by earlier work that demonstrated the need of a rich dynamics for the turn rate to faithfully reproduce key elements of the behavioural phenotype (Gautrais et al., 2009;Mwaffo et al., 2015).
The response function ★ f (t) accounts for social interactions, which primarily arise from visual cues and result in the tendency of the animals to align their bodies during swimming (schooling) and swim in the vicinity of others (shoaling), see, for example, the paper by Miller & Gerlai (2011). We write the response function for fish f , as originally proposed by Calovi et al. (2014Calovi et al. ( , 2015, as where k v,f and k p,f are gain parameters for alignment and attraction, respectively (Figure 1a), which summarize multisensory cues that underlie social response. These gains are the means by which one can simulate social behaviour, in the form of an animal orienting its swimming direction to align with and approach others. The prefactor (1 + cos f ,f (t)) accounts for the cone of vision of zebrafish (Pita et al., 2015), which preferentially biases interactions toward animals that are swimming in front rather than those behind. As discussed by Calovi et al. (2015), this prefactor 'introduces a strong asymmetry between the force exerted by f onf and the one exerted byf on f , and hence breaks the (Newtonian) action-reaction principle which is most familiar in the context of purely physical force, such as gravitation'.
The response function is such that the contribution to the imposed turn rate from the alignment rule tends to orient the swimming direction of f to align withf . Thus, it is zero when the animals are swimming in parallel ( f ,f (t) = 0 or π) and it is maximized in magnitude when they are swimming orthogonal to each other ( f ,f (t) = ±π/2 and f ,f (t) = 0). We choose to scale the attraction gain by the velocity, which is a constant, to acknowledge that faster fish will adjust their swimming direction more rapidly in response to external stimuli. However, the contribution from the attraction rule is such that f will tend to correct its swimming direction to approachf : it is zero when f ,f (t) = 0 and maximized in magnitude when f ,f (t) = ±π/2. In this model, alignment is taken to be independent of the relative distance between the animals, while attraction increases the further the animals are from one another. These approximations are expected to be valid when fish swim within approximately ten body lengths; otherwise, it may be appropriate to include a decay function that reduces the strength of the interaction for large distances (Zienkiewicz et al., 2015b).
Within the finite-dipole framework, the turn rate imposed by hydrodynamic interactions f (t) arises from differences in the velocity experienced by each vortex comprising the dipole of fish f . Thus, in the far-field limit, the gradient of streamwise velocity in the direction orthogonal tov f induces the hydrodynamic turning is the unit vector orthogonal tov f (t) forming a local right-handed coordinate system and ∇ is the gradient operator in the global coordinate system. With respect to the paper by Filella et al. (2018), we note the difference in sign between (8) and (5) therein, which arises from our choice to work with transverse versus aligned dipoles, that is, T-dipoles instead of A-dipoles, as defined by Kanso & Tsang, 2014. One can verify that the sign is indeed negative by considering a fish swimming in the ; in this case, the turn rate in (8) is − such that the fish will turn clockwise, as one would expect also based on the work of Tchieu et al. (2012). For completeness, we report results also for the case of A-dipoles in the supplementary material. By using the expression of the velocity U f (r f (t)) in (5), along with the Cartesian representation ofv ⊥ f (t) andv f (t), we establish where the minus (plus) sign is for f equal to 1 (2). Comparing (9) with (7), we register a contrasting dependence on the distance between the two fish, whereby social interaction linearly increases with (t), while the hydrodynamic turning decreases with (t) 3 . Added noise in (6b) captures the uncertainty in the behaviour of the animal, which may not necessarily balance between social and hydrodynamic interactions. Zebrafish locomotion is characterized by a unique burst-and-coast swimming style in which 'fish move forward (burst) in a single motion and glide (coast) to a slow speed, or stop from which they burst forward again' (Kalueff et al., 2013). This swimming style contributes to a rich repertoire of locomotory patterns, which include sudden U-and C-turns that alternate with instances of straight swimming (Kalueff et al., 2013). As a result, uncertainty in turn rate dynamics cannot be captured by Gaussian white noise (that is, increments of a Wiener process), as traditionally advocated for other swimming animals (Gautrais et al., 2009). This issue has been addressed by Mwaffo et al. (2015) by modelling added noise as the sum of two terms: (i) a Wiener process f W (t), where W (t) is a standard Wiener process and f is a positive constant (note that for a standard Wiener process, dW (t) ∼ N (0, dt), where N ( , 2 ) is a normal distribution with mean and variance 2 , so that multiplying by f scales the variance of the normal distribution by 2 being a Bernoulli distribution with parameter q). The Wiener process represents the baseline uncertainty with animal locomotion. The compound Poisson process allows the capture of jumps in the turn rate dynamics, which are associated with C-and U-turns. A sample video of a zebrafish pair swimming is included as a supplementary movie. Overall, added noise is modelled through three independent parameters: f , setting a baseline activity for the animal; f , quantifying the frequency of the jumps; and f , measuring the strength of the jumps.
In light of the choice of the form of added noise, (6b) can be viewed as a jump reverting mean diffusion process (Tankov, 2003), in which the turn rate relaxes toward a time-varying mean, given by the sum of the turn rates imposed by social and hydrodynamic interactions in (7) and (8), respectively. The original model by Filella et al. (2018) has neither the feature of reversion nor the presence of jumps. Hence, at every time, the turn rate of each fish matches the sum of the turn rates imposed by social and hydrodynamic interactions, subject to added Gaussian white noise (note that in the model by Filella et al. (2018) where f has units rad s −1 rather than rad s −3/2 as in our study).

Analytical Insight into the Model
Predicting the behaviour of a zebrafish pair requires the time integration of (1) and (6). While the general mathematical treatment of this system may not be feasible, important insight can be gathered by examining its linearized dynamics about salient configurations that have previously been documented (De Bie et al., 2020;Kato et al., 2004;Laan et al., 2017;Perna et al., 2014;Sadaat et al., 2021). To further simplify the problem, we assume the two fish to be identical, such that r 0, Following standard practice in the literature on synchronization of coupled dynamical systems (Pikovsky et al., 2003), we decompose the dynamics of the fish pair into average and error dynamics. The former captures the evolution of the centre of mass of the pair and of the average orientation, whereas the latter describes changes in relative distance and orientation. Hence, we write where the plus (minus) sign is for f equal to 1 (2), a superimposed bar indicates the average dynamics and the symbol Δ is consistently used to label relative dynamics.

E7-7
By using this coordinate transformation, we consider the eight-dimensional state variable z(t) = [x(t),ȳ(t),¯(t),¯(t), Δx(t), Δy(t), Δ (t), Δ (t)] T , where T indicates matrix transposition. The time evolution of the state is given by where u(t) ∈ R 2 encapsulates added noise to the turn rate dynamics of the animals and the nonlinear function F : R 8 × R 2 → R 8 captures the first-order dynamics of the turn rate in (6b) along with the coupling between the state variables of the two fish through the advection velocity (2), social interactions (7) and hydrodynamic interactions (8).
The class of nominal solutions about which we perform a linearization is part of the schooling patterns identified by Filella et al. (2018), wherein fish swim along the same, constant direction, that is, Δ (t) = 0 and¯(t) = c for some c . Fish will follow a straight line (Δ (t) = 0 and¯(t) = 0) and rigidly translate one with respect to each other, as shown in Figure 1b. In this in-line swimming pattern (sometimes called 'tandem'), one fish follows the other, such that Δx(t) = d cos c and Δy(t) = d sin c , where d is the constant distance between the fish. Without loss of generality, we set c = 0, such that the x-axis of the Cartesian coordinate system is aligned with the schooling direction.
It is easy to verify that for any choice of d, in-line swimming constitutes a nominal solution of (11) in the absence of added noise. Specifically, both the imposed turn rates owing to social and hydrodynamic interactions in (7) and (8) are identically zero, such that the orientation of each fish remains constant , 2}, such that each fish will be subject to the same advective velocity, thereby preserving the relative configuration. The centre of mass of the pair will move at a constant speed given by v 0 (1+r 2 0 /d 2 ), which yields the increase in speed arising from hydrodynamic interaction observed through simulations on large groups (Filella et al. 2018).
The homogeneous, linearized dynamics of the system is informative of the local stability of the inline configuration. By taking the gradient of the nonlinear vector field in (11) and evaluating it at the nominal solution, we obtain the following state matrix describing the evolution of perturbations with respect to the nominal, in-line solution: Here and henceforth, a superscript identifies perturbations with respect to the nominal solution.
Equation (12) demonstrates several interesting features of the model. First, perturbations on the position of the centre of mass do not enter the dynamics of any of the other state variables. Second, the relative distance along the schooling direction remains constant, independent of any other state variable. Third, the perturbation on the relative distance Δy (t) depends only on the relative rotation Δ (t). Fourth, the response of the fish pair is asymmetric, owing to the presence of social interactions which are sensitive to the relative orientation of the animals (leading versus trailing fish); asymmetry is evident in the turn rate dynamics by combining (12d) and (12h) to obtain governing equations for 1 (t) and 2 (t).
To further clarify the structure of the system, we perform a coordinate transformation that accounts for the specific nature of in-line swimming. Specifically, we define¯★ (t) =¯(t) − Δy (t)/d to filter out rigid-body motion in the definition of the average angle; in fact, in-line swimming would not be affected if both fish moved along the y-axis while rotating their bodies to preserve straight swimming. Consistently, we define¯★ (t) = ¯★ (t). By using these variables in (12) and ordering the components of the state vector as z (t) = [x (t),ȳ (t), Δx (t), Δy (t),¯★ (t),¯★ (t), Δ (t), Δ (t)] T , we can write the state matrix of the system in the following block-triangular form: The stability of in-line swimming to external perturbations is ascertained by examining the spectrum of the four-by-four block on the bottom diagonal, A 22 . In-line swimming is asymptotically stable if, and only if, this sub-matrix is Hurwitz. If A 22 is Hurwitz, the orientations converge to a common value (Δ (t), Δ (t) → 0), such that the two fish swim along a straight line (¯★ (t),¯★ (t) → 0). The other blocks have no role on stability, whereby the location of the centre of mass and the relative distance Δx (t) are inconsequential with respect to in-line swimming, and Δy (t) converges to a constant value that depends on the common orientation attained by the fish pair.
While a precise evaluation of the stability of in-line swimming requires the calculation of the spectrum of A 22 , important insight can be gathered from the application of the Routh-Hurwitz criterion. Briefly, application of the criterion requires the construction of an ancillary table of coefficients, which is computed from the coefficients of the characteristic polynomial (Meinsma, 1995). By counting the number of sign changes in the first columns of this table, it is possible to exactly count the number of eigenvalues with positive real part. Going through simple algebra, we obtain the following expression for the first column of the Routh table R = {R 4 , R 3 , . . . , R 0 }: Flow E7-9 In general, it is difficult to tell the sign of the terms in (14), but several interesting limit cases can be explored. For example, for a sufficiently large d, we discover that there are no sign changes (all terms being positive), which demonstrates that in-line swimming is asymptotically stable for sufficiently large separation distances. However, for small values of d, there is one sign change (R 1 → −∞), which indicates the presence of one eigenvalue with a positive real part.
We can also examine the case in which fish interact only hydrodynamically, such that both k p and k v approach zero. In this case, we register at least one sign change, whereby R 0 = −12 2 r 2 0 v 2 0 (d 2 + 2r 2 0 )/ d 6 < 0. This suggests that visual cues are needed to ensure the stability of in-line swimming, which supports the simulation results in Figure 7b of Tchieu et al. (2012) and echo analytical predictions on large-head microswimmers (Kanso & Tsang, 2015) and suspensions of phoretic particles (Kanso & Michelin, 2019).
Another interesting case is the limit of large values of , which is reminiscent of the set-up by Filella et al. (2018); in this case, R 0 / 2 = −6(d 2 + 2r 2 0 )r 2 0 v 0 + (d 2 + 5r 2 0 )d 4 k p + 6d 3 r 2 0 v 0 k v > 0 is the necessary and sufficient condition for the matrix to be Hurwitz, because all the other terms are positive. For small values of d, the leading term is −12r 4 0 v 0 , which confirms the instability of in-line swimming, unless k v and k p are sufficiently large to dominate R 0 . Predictably, large values of k v and k p will favour stability for any choice of d. Should hydrodynamic interactions be dismissed by setting r 0 = 0, we would always obtain asymptotic stability.
Another swimming pattern that is often considered in the literature is side-by-side swimming (De Bie et al., 2020; Kato et al., 2004;Laan et al., 2017;Perna et al., 2014), where the fish swim parallel to each other with Δx(t) = −d sin c and Δy(t) = d cos c . For the present model, this swimming pattern does not constitute a nominal solution: hydrodynamic interactions will have an equivalent effect as alignment and attraction, thereby causing the fish to approach one another.
The analysis presented thus far relies on the classical finite-dipole paradigm of Tchieu et al. (2012) to study fish schooling, which was later referred to as the transverse (T) dipole to distinguish from the aligned (A) dipole (Kanso & Tsang, 2014). As noted by Kanso & Tsang (2014), the two models only differ in 'how an individual dipole aligns itself with the local flow gradient'. The T-dipole aligns itself based on the difference in the flow velocity in the direction transverse to the dipole, while the A-dipole aligns itself according to the difference in the flow velocity along the dipole direction, similar to the active swimmers in microfluidic devices considered by Brotto et al. (2013). For completeness, we report analytical insight into the model also for the case of A-dipoles in the supplementary material.

Experimental Data and Model Calibration
We calibrated the model using experimental data from a previous study (Zienkiewicz et al., 2015a); the same dataset was also used by Butail et al. (2016) and is included in the supplementary material. The dataset contains 17 different pairs of adult zebrafish with an average body length (BL) of 30 mm swimming in a circular arena of 0.9 m diameter and water depth of 0.1 m. The shallow depth encouraged the fish to swim in a quasi-two dimensional plane. Fish motion was filmed at 30 frames s −1 for 20 min, and a multi-target tracking system (Ladu et al., 2014) was used to track the centroids of each individual. The first 10 min were treated as a habituation phase, in which the animals explored the novel tank and acclimatized to the environment.
The last 10 min of the videos were used for the analysis. Each pair of experimental trajectory data consists of positions r f (t) and orientations f (t) in the instantaneous direction of motion of the two animals sampled at a time step Δt = 1/30 s. We filtered the data in two steps. Time windows where the two fish were closer than half a BL from each other were removed to avoid diving instances, which would violate the assumption of planar swimming. To mitigate the effect on unmodelled hydrodynamic interactions with the walls, we also removed instances when either fish swam closer than two BLs to the wall. Overall, these pre-processing steps curtailed approximately half of the dataset.
This same dataset was also used to validate the predictions of the model, by evaluating salient metrics of collective behaviour that emerge within the pair. Specifically, we considered the average polarization, which measures the relative orientation of the pair: Pol = (1/T) ∫ T 0 ( v 1 (t) +v 2 (t) /2) dt, where T is the duration of the observation (Calovi et al., 2014). A polarization value of one indicates that the swimming directions of both fish are always aligned during he observation window, whereas a value of zero identifies that the animals are swimming in opposite directions. While polarization helps to quantify schooling, it does not distinguish between in-line and side-by-side swimming, which we capture through a pursuit index defined as PI = (1/T) ∫ T 0 ((v 1 (t) +v 2 (t)) · (r 1 (t) − r 2 (t))/2 (t)) dt. A value of 1 (−1) indicates in-line swimming, with the first (second) fish swimming in front of the second (first) fish. A zero value corresponds to the pair swimming side-by-side. The average inter-individual distance between the pair of fish, Dist = (1/T) ∫ T 0 | (t)| dt, is used to quantify their shoaling tendency, that is, their preference to stay close to one another. Finally, we measured the average speed of the centre of mass Speed = (1/T) ∫ T 0 ( r 1 (t) + r 2 (t) /2) dt to further detail the effect of interactions on locomotion. For each fish, we identified eight model parameters (characteristic length, r 0 ; speed, v 0 ; baseline activity, f ; reversion rate, f ; jump frequency, f ; jump intensity, f ; and attraction and alignment gains, k p,f and k v,f , respectively) from the equations of motion (1) and (6). We performed the identification in two sequential steps. First, we calibrated characteristic length, r 0 , and speed, v 0 , by minimizing the sumsquare-error cost function (SSECF) calculated from the experimentally measured speed and the speed in (1), defined as SSECF = k∈Data f ( r f (t k ) − v fvf (t k ) + U f (r f (t k )) ) 2 , where 'Data f ' identifies the subset of indices that remain after the removal presented above for fish f . The optimization was solved by constraining each parameter within a physically plausible range. Specifically, based on prior data on fish swimming in isolation and the typical body length of a zebrafish (Ladu et al., 2014), we selected v f ∈ [0, 0.30] m s −1 and r f ∈ [0, 30] mm. Calibrated parameters are presented in Figure 2a,b.
Second, we identified the remaining six parameters k v,f ] T using a maximum likelihood approach on the turn rate dynamics in (6b). The maximum likelihood estimation was subject to further data preprocessing, whereby we removed instances when fish were swimming approximately in-line (AI > 0.75); at these instances, the prefactor in the social interaction (7) would cause excessively small values of the response function, which would challenge robust identification.
Hence, we estimated the model parameters by minimizing the sum of negative logarithms of the likelihood function L( f (t)| f (t − Δt), Θ f ) over experimental time-series using the FMINCON function in MATLAB (Butail et al., 2016). The likelihood function was constructed as the sum of two normal probability distributions associated with the two forms of added noise, where and H(•, , 2 ) is the probability mass function of N ( , 2 ).

Figure 2. Calibrated model parameters: (a) characteristic length; (b) speed; (c) baseline activity; (d) mean reversion rate; (e) jump frequency; (f) jump intensity; (g) gain parameter for attraction; and (h) gain parameter for alignment. The coloured area in each violin plot is the probability density and coloured circles are individual calibrations. Thick grey bars indicate first and third quartiles; thin grey bars identify minimum and maximum values; and white circles are the median. Calibrated parameters for each fish are reported in the supplementary material.
Based on the social model of Butail et al. (2016), we constrained the parameter ranges for the search as follows. First, we set f ∈ [1, 2] s −1 , f ∈ [1, 5] rad s −3/2 and f ∈ [1, 5] rad s −1 . Then, we explored the range in jump frequency of f ∈ [0.9, 1.1] 0, f , where 0, f was an estimate of the frequency of extreme events in the time-series, defined as the ratio between the number of time steps above three standard deviations from the mean of the turn rate and the total number of time steps. With respect to interactions parameters, we explored the broad range k p,f ∈ [0, 100] rad m −1 s −1 and k v,f ∈ [0, 100] rad m −1 to capture inter-individual variability in social behaviour.

Results and Discussion
By applying the calibration procedure explained above, we obtained the complete set of calibrated parameters, as shown in Figure 2. The values of r 0 were one order of magnitude less than the animal body length (r 0 = 3.1 ± 3.2 mm), which corresponded to the typical tail beat amplitude of a cruising zebrafish . The self-propulsion speed ranged from 2BL to 6BL per second (v 0 = 0.094 ± 0.027 m s −1 ), which suggested variability among animals with respect to their locomotion. Interestingly, these two variables were highly correlated, such that larger distances between vortices comprising the dipole will be accompanied by larger speeds, thus implying larger values of the vortices' circulation, see Table 1. This trend was analogous to the relationship between tail beat amplitude and the circulation of vortices shed during tail beating in the real experiments by Mwaffo et al. (2017).
The baseline activity ( = 1.66 ± 0.61 rad s −3/2 ), mean reversion rate ( = 1.77 ± 0.32 s −1 ), jump frequency ( = 0.64 ± 0.12 s −1 ), jump intensity ( = 2.63 ± 0.77 s −1 ), and gain parameters for attraction (k p = 3.08 ± 3.43 rad m −1 s −1 ) and alignment (k v = 18.12 ± 13.12 rad m −1 ) were comparable with those of Butail et al. (2016), which did not include a model of hydrodynamic coupling. These parameter values were not all independent of each other (Table 1). We determined a positive relationship between E7-12 M. Porfiri et al. Table 1. Coefficient of determination (R) and associated p-value in parentheses between model variables. Significant correlations at a level of 0.05 are indicated in bold (should one correct for multiple comparisons using a conservative Bonferroni correction, some significant correlations would become trends). baseline activity and jump intensity, which suggested that animals exhibiting higher activity levels also displayed stronger bursts. We also recorded a positive correlation between r 0 and both the noise intensities, and , which suggested that larger animals will tend to be more active and display more pronounced body movements during C-and U-turns. Finally, we discovered that the alignment gain was positively correlated with the baseline activity, which suggested that more active animals might exhibit a stronger social response, perhaps pointing at a deeper association between locomotory activity and social response that might be mediated by stress (Kalueff et al., 2013). Overall, the extent to which hydrodynamic interactions improve the quality of the fit can be ascertained using the Akaike information criterion (Akaike, 1998). Specifically, for the present model in (6b), this quantity is equal to 52.10 ± 33.33. Should we compute the same for the model by Butail et al. (2016), which excluded hydrodynamic interactions, we would obtain a strikingly similar result of 46.74 ± 32.54. Hence, accounting for hydrodynamic interactions does not lead to a tangible improvement on the fit of the turn rate, although it allows the capture of speed-based interactions that would be otherwise neglected. Evidence in favour of the existence of speed-based interactions has been documented in a number of experimental studies on freshwater fish (see, for example, Herbert-Read et al., 2011;Katz et al., 2011), but their connection to hydrodynamics has remained elusive. The supplementary material contain calibrated model parameters for A-dipoles, which are highly comparable to those presented herein.
To demonstrate the model predictive power, we performed a series of in silico experiments, which employed the calibrated parameters from Figure 2; a sample video is included as a supplementary movie. This analysis aimed to ascertain whether the model predicts equivalent statistical outcomes to real experiments, thereby allowing for partial replacement of animal trials. We ran the model 17 times to mirror experimental observations, using identical parameters for each fish in the pair (drawn from normal distributions approximating those in Figure 2) and random initial conditions (orientations drawn uniformly in [0, 2π] and distance within 4r 0 , similar to Tchieu et al. (2012)). For each metric of collective behaviour, we performed t-test comparisons between chance and either real or in silico experiments. Chance was computed by randomly shuffling fish across the 17 experimental observations 10 000 times and taking the mean.
As shown in Figure 3, the model was successful in anticipating a strong schooling tendency of the pair (t (1,16) = 20.67, p < 0.001 (experiments); t (1,16) = 6.09, p < 0.001 (in silico)). During schooling, animals favour in-line swimming, which was accurately predicted by the model (t (1,16) = 14.74, p < 0.001 (experiments); t (1,16) = 10.15, p < 0.001 (in silico)). Although to a lesser quantitative extent, the model was also able to anticipate shoaling tendency of the pair (t (1,16) = −31.20, p < 0.001 (experiments); t (1,16) = −3.94, p = 0.001 (in silico)) and the increase in speed arising from animal interactions (t (1,16) = 3.31, p = 0.004 (experiments); t(1, 16) = 2.12; p = 0.05 (in silico)). k v (rad m -1 ) k p (rad m -1 s -1 ) k p (rad m -1 s -1 ) k p (rad m -1 s -1 ) Experimental subjects tend to swim within 2 to 3 BLs, which correspond to approximately one-third of the theoretical predictions, thereby leading to an underestimation of short-range interactions. The source of the discrepancy lies with the absence of social interactions in the speed dynamics, which could, however, be included in the model by following Zienkiewicz et al. (2015b). Interestingly, using A-dipoles does not challenge the predictive power of the model, as presented in the supplementary material. Overall, Figure 3 points to a strong tendency of live fish to school in an in-line pattern, which is accurately captured through in silico experiments. Further insight can be garnered through the analysis of the eigenvalues of A in-line,22 in (13), which was derived from linearization of the proposed model. Specifically, we used the calibrated values of r 0 , v 0 and from Figure 2 to analyse the stability of in-line swimming with respect to k p , k v and d in Figure 4. Numerical computation of the eigenvalues confirmed the theoretical proposition from the Routh-Hurwitz criterion that small values of d hinder stability. For sufficiently large values of d, there was a wide range of combinations of the gains that exhibited stability. Using calibrated values of k p and k v , we always attained stability for any choice of d above approximately one BL, in agreement with observations on inter-individual distance in Figure 3. In contrast to the kinematic model of Filella et al. (2018) that suggests a beneficial role of large values of k p and k v , we discovered a richer scenario where large values of k p may destabilize the swimming pattern. A sample video of the effect of perturbations on in-line swimming is included as a supplementary movie.
For completeness, in Figure 4, we also present the results obtained by setting r 0 = 0, that is, dismissing hydrodynamic interactions without modifying any other model parameter. Hydrodynamic interaction effectively leads to a repulsion zone between animals swimming in-line that would otherwise not be present. Repulsion has been widely documented in the literature on mathematical models of animal groups, but generally attributed to a social response (Vicsek & Zafeiris, 2012), with the exception of  who also documented hydrodynamically-based repulsion. For a given value of d, we also registered a wider range of gain combinations that guaranteed stability; however, the extent of the difference became negligible for distances above one body length. Comparable insight can be gathered for A-dipoles, as presented in the supplementary material.

Conclusions
Here, we have presented and analysed a mathematical model of collective behaviour of zebrafish pairs. The model combines the jump persistent turning walker of Mwaffo et al. (2015) with the interaction model by Filella et al. (2018) toward a realistic description of zebrafish locomotion and of social and hydrodynamic interactions. The model consists of a system of coupled nonlinear stochastic differential equations for the speed and turn rate of each fish. Model parameters were calibrated using the real data on zebrafish pairs swimming in shallow water of Zienkiewicz et al. (2015a).
The calibrated model was successful in predicting all of the main features of the collective response of live animals. Not only did the animals exhibit strong shoaling and schooling tendencies, by swimming in close physical proximity and aligning their bodies, but they also opted for an in-line pattern where one fish would lead and the other would follow. A contribution of this study was to demonstrate local stability of this specific collective pattern through a detailed analysis of the linearized dynamics. Above a critical value of the inter-individual distance, in-line swimming becomes stable for a wide parameter range. This is in stark contrast with side-by-side swimming, which does not emerge as a viable solution of the model.
Taken in toto, these results question some of the assumptions that are often used in modelling fish behaviour (Vicsek & Zafeiris, 2012), which explains the mechanisms underlying schooling (Partridge & Pitcher, 1979;Weihs, 1973) and draws inference of the leader-follower relationship from experimental observations of their orientation (Krause et al., 2000). First, our theoretical results indicate that repulsion emerges from hydrodynamic interactions, rather than being a social rule pursued by individuals. Second, the long debate regarding the determinants of schooling patterns, with some fish occupying frontal positions and other rear ones, may reduce to a simple stability problem for the group. Neither reduced swimming costs nor foraging benefits are needed for fish pairs to opt for inline swimming. Third, the inference of leadership traits from the frequency of instances in which an animal initiates a turn that is followed by others could be skewed by the stability of in-line swimming. Specifically, should the animals differ in their noise parameters, we would detect a false leader-follower association.
The main limitation of the model is that it neglects the reverse von Kármán wake generated by self-propelled swimmers, which results in a poor representation of the flow field immediately behind the dipole. The flows around the fore-body and to the sides of the dipole are, however, representative of those produced by a swimming fish. This is likely relevant in the study of in-line swimming, by adding further dynamics into the collective response (Liao et al., 2003) and potentially mitigating the short-range decay of hydrodynamic interactions observed in the present model. While the vortex dipole model for fish swimming has been used in a variety of analytical models of swimming, there exists a need for explicit experimental validation.
There are several additional directions for the potential improvement of the model, which include: (i) the development of a dynamic model for the speed in the form of a stochastic different equation, similar to that of Zienkiewicz et al. (2015a,b) to overcome the approximation of constant speed and include additional mechanisms of social interaction; (ii) the extension to hydrodynamic interactions with walls in the tank and background water currents toward an improved understanding of the transition from in-line to side-by-side swimming as a function of the flow speed (De Bie et al., 2020); and (iii) the extension of individual differences in the stability analysis that could help explain the modulatory role of fish body size on schooling and shoaling (Karakaya et al., 2020;Reebs, 2001).