Leadership emergence in a data-driven model of zebrafish shoals with speed modulation

Models of collective animal motion can greatly aid in the design and interpretation of behavioural experiments that seek to unravel, isolate, and manipulate the determinants of leader-follower relationships. Here, we develop an initial model of zebrafish social behaviour, which accounts for both speed and angular velocity regulatory interactions among conspecifics. Using this model, we analyse the macroscopic observables of small shoals influenced by an “informed” agent, showing that leaders which actively modulate their speed with respect to their neighbours can entrain and stabilise collective dynamics of the näıve shoal.


Introduction
Social fish exhibit a diverse and complex repertoire of emergent behaviours, which are exemplified by their ability to maintain stable and coherent dynamical states and to react collectively to stimuli from within the group and their environment [1][2][3][4][5][6][7][8]. The richness of the spatio-temporal patterns exhibited by fish schools originates from the continuous interaction between members of the group [9][10][11][12][13][14]. The evolutionary advantages of this coordinated, group-level behaviour are often associated with predator avoidance [15], foraging [11], and migration [16]. In each instance, interaction within the group enables rapid information sharing so that individual exposure to stimuli can trigger a collective response which benefits the entire shoal.
Sometimes, the group behaviour resulting from these individual responses is intuitive, for example a fish turning away from a predator may elicit a similar response from its neighbours such that the entire group is spared. Less intuitive is the emergence of distributed sensing in which shoals navigate a complex environment via a combination of attraction between group members and their individual response to local cues [17]. In both instances, we see how unequally distributed information, initially limited to a few individuals, can propagate throughout the group eliciting an appropriate response. To comprehend the emergence of group-level behaviour, it is thus important to understand and quantify how information propagates via a e-mail: m.dibernardo@bristol.ac.uk interactions between individuals. This has motivated a recent surge in research with the specific aim of inferring the precise mechanisms underpinning interactions in animal groups [18,19], and specifically schooling fish [8,20,21].
The dynamics of informed individuals and their effects on collective behaviour has been the subject of numerous theoretical and experimental studies. Both have found that a small subset of the population, trained with pertinent knowledge, such as the location of a food source, is able to influence the swimming direction of the shoal [22,23]. Effects of controlled external stimuli have also been studied, specifically the use of biomimetic robots (ethorobotics) to elicit responses both at the individual level [24][25][26] and in the collective behaviour of fish shoals [27][28][29][30]. Meanwhile, purely computational studies have explored the dynamic response of biologically inspired mobile agents in the presence of an informed agent or group of agents [31][32][33].
Here, we take a data-driven, bottom-up approach, aimed at deriving from experimental observations, models of the interactions between conspecifics to investigate the effects of speed regulation on leader-follower interactions in zebrafish shoals. Building on recent efforts to infer details of the specific interactions among teleost fish [18][19][20][21], we extend the individual zebrafish model formulated in [34] to study their social behaviour in small shoals. Specifically, the model we describe incorporates speed regulation, which is both fundamental to the locomotory patterns of individuals [35,36] and has been recently proposed as a central mechanism for explaining collective behaviour of similar teleosts [17,20,21,37,38].
To investigate the emergence of leadership, we consider a single informed agent moving on a linear trajectory, which interacts with the 'naïve' shoal by modulating its forward speed. We find optimum values for the strength of the interaction, which enable an informed agent to lead and stabilise naïve group dynamics more effectively than one which moves at a constant speed. We adapt a variety of global observables (polarisation, cohesiveness, and milling) to explore the effect of the informed agent and assess its ability to "control" collective behaviour.

Data-driven stochastic model of individual zebrafish locomotion
The model used to describe the individual (isolated) swimming dynamics of zebrafish was originally developed in [34]. Similarly to recent efforts on the so-called persistent [8,39] or jump persistent [40] turning walker, the model uses two coupled stochastic differential equations (SDEs) in time t, one for the forward speed U (t), and one for the turning (angular) rate Ω(t) 1 where the coupling function f c is given as follows: Here, the parameter σ 0 defines a maximum volatility, where the function f c decays to zero passing through f c = σ ω when U (t) = μ u . This relationship was estimated from experimental observations which indicated a strong negative correlation between the speed of the fish and the volatility of the turning rate [34]. The smooth random walks constructed after integrating these equations in time have time-averaged speed that approaches μ u and directionally unbiased turning, i.e. null mean turning rate. The mean-reversion (autocorrelation) rates θ u and θ ω for each process along with the variances of the stochastic Wiener processes dW (t) and dZ(t), denoted by parameters σ u and σ ω respectively, are all calibrated directly from experimental trajectory data using a maximum likelihood method (also described in [34]).
To ensure forward only motion of the simulated fish, we limit the resulting integration of (1a) to positive values by setting any negative outputs to zero. Finally, we prescribe a maximum value for the turning speed Ω max = 15 rad s −1 , used to truncate the output of the SDE in (1b) such that turning behaviour remains within experimentally observed limits.
2.2 Multi-agent dynamics: Data-driven modelling of zebrafish shoals

Experimental observations of zebrafish pairs
The experimental data used to validate and calibrate the model developed in this work was obtained at the Dynamical Systems Laboratory (NYU Polytechnic School of Engineering, NY, USA) 2 . Using a high-resolution overhead camera (1280×960 pixels), we recorded videos of pairs of wild-type zebrafish (Danio rerio) swimming together in a shallow circular tank with radius of 45 cm and water depth of 10 cm, such that that swimming is restricted to a quasi two-dimensional plane. Fish were 6-8 months of age, with a mean body length (BL) around 3 cm. In total, 18 pairs of unique individuals were observed for 30 min each including 10 min of habituation to the novel environment which was discarded from subsequent analysis [41]. The protocol is similar to that described in [42].
Video image analysis and multi-target tracking was achieved using software developed in MATLAB (R2011a, MathWorks), sampled at 30 Hz [30]. Raw experimental data consisted of two-dimensional Cartesian positions x(t) = [x, y](t) for each fish. Position data was smoothed using a third-order Savitsky-Golay filter with a moving window of 15 frames (0.5 s) prior to the numerical computation of velocity and acceleration via successive application of central differences. The forward speed u(t) = v(t) and the turning rate ω(t) were computed from the velocity vector; for further details see [34,43].
To infer interactions between fish, we use the same force mapping method described in [20] to compute the reaction force (acceleration) of a focal fish as a function of the relative position of its neighbour. For each data frame, the acceleration a i of a focal fish i is decomposed into two components: a speeding acceleration a s i tangential to the direction of motion, and a radial turning acceleration a t i in the perpendicular direction. The position of a neighbouring fish j is then expressed in terms of its front-back distance d F B and its left-right distance d LR , relative to the focal fish at the origin, with its velocity vector aligned with the positive y-axis.
For this analysis, we discretise the relative positional space into bins defining a square 30 × 30 grid extending up to 4 BL (12 cm) along both axes. Selecting each of the two fish as the focal fish, values of a s i and a t i are accumulated as separate two-dimensional histograms. The resulting histograms for both speeding and turning components of the acceleration are then normalised by the bin population density and averaged over the two focal-neighbour combinations. We exclude potentially strong  wall effects by rejecting frames in which either fish is closer than 2 BL to the boundary. Frames in which the speed of either fish is less than 0.5 BL s −1 are also rejected to reduce excessive fluctuations from tracking errors at low (or stationary) fish speeds.
The force (acceleration) maps shown in Fig. 1a (top two panels) are composite histograms for the complete (18 × 20 min) data set of 18 pairs of zebrafish, normalised appropriately to provide average component accelerations in units of BL s −2 . For both components we find that variations in acceleration occur in a single directional axis, namely the front-back (FB) axis for a s i and the left-right (LR) axis for a t i . With respect to the speeding acceleration, we find that within a small region around the focal fish (< 1 BL) there is a repulsive interaction, characterised by deceleration (a s i < 0) when a neighbour is just ahead of the fish, and positive acceleration (a s i > 0) when immediately behind. At greater distances where d F B > 1 BL, an attractive interaction dominates, with the focal fish accelerating in response to a neighbour far ahead and, conversely, slowing down when a neighbour is further behind. Along with very similar results golden shiners, and mosquitofish [20,21], we conclude that speed modulation is an important factor for zebrafish interaction. The turning acceleration map closely resembles that of the speeding maps, albeit varying primarily as a function of d LR . Within a small radius (<0.5 BL) around the focal, we find that the focal fish accelerates in the opposite direction to its neighbour's position; adopting the convention that negative a t i represents an acceleration to the left and vice versa. At greater distances, we see a strongly attractive region in which the turning "force" on the focal fish is in the direction towards its neighbour.
Since both speeding and turning forces are found to be principally governed by the separation in orthogonal axes (front-back and left-right respectively), we can summarise the interaction by projecting both forces along the relevant axis, shown in Fig. 1(b) (black dashed lines). These force projections clearly show both linear and repulsive regions of the interactions and provide a more convenient means of calibrating the interaction model discussed in the following section.

Multi-agent model
Our interaction model is adapted from the previous work of Gautrais et al. [8] and subsequent modifications by Calovi and others in [33,44]. We derive a minimal model in which each contributing term, and associated parameters, are constrained by features observed directly from the experiment.
Component accelerations of fish i, revealed from the force mapping analysis, are explicitly included in the SDEs governing the time evolution of speed and angular velocity in (1) by incorporating response functions U * i (t) and Ω * i (t), resulting in the following (multi-agent) SDEs: The response functions are described by normalised linear superposition of pairwise interactions between the focal fish i and each of its neighbours, as follows: Here, N i is the number of fish in the first shell of a periodic Voronoi neighbourhood . Relative spatial metrics are described by three variables: the signed angle θ ij between the focal fish, given its current velocity vector and its neighbour's position; the signed relative heading angle φ ij between focal fish i and neighbour j; and the distance d ij between focal fish i and neighbour j. Attractive and repulsive regions of the speeding and turning interactions are introduced in both equations in (4), with a term proportional to the separation distance d ij , offset by a constant value r u or r ω respectively. The resulting forces are therefore linear with distance but will have an opposite sign (repulsion) with respect to relative position of the neighbouring fish when the separation is between ±r u/ω . The directionality of the speeding and turning response are modulated by a factor of cos θ ij and sin θ ij respectively, such that |U * i | is maximised when a neighbour is located directly ahead or behind, whilst |Ω * i | is maximised when a neighbour is located directly to the left or right, relative to the focal fish's orientation v i . The strength of these terms is controlled with gain factors K s for speeding and K p for turning.
An extra term in the turning interaction (4b), dependent on the relative heading angles sin φ ij , provides a turning response acting to align the focal fish with its neighbour, with associated gain parameter K v . Evidence for this explicit alignment interaction, to be presented in a future study, is revealed by exploring how the angular accelerationω of zebrafish varies as a function of the relative pair orientation φ ij . The turning interaction function is further supplemented by an additional prefactor (1 + cos θ ij ) which imparts an angular weighting on the turning interaction, maximised when the neighbour is directly in front with respect to the viewing angle of the focal fish. As described in [33,44], this extra position-dependence breaks the otherwise symmetric force response allowing for stable rotational (milling) and winding collective phases to emerge, producing qualitatively better fits to the group behaviour of small zebrafish shoals.
In contrast to the previous models mentioned in [8,33], we also provide a function f d which restricts the range of both the speed and turning response, defined as follows: Here, the constant d c prescribes a hard, radial cut-off distance to the interaction, where f d decays from a maximum at one, dropping to zero with a rate given by parameter δ. With f d as a multiplicative factor, we thus obtain responses which are roughly proportional to the distance d ij at short distances, and decay rapidly at the cut-off distance d c -preventing unrealistically strong interactions between distant neighbours [8]. Finally, both functions in (4) are normalised by the associated SDE mean-reversion rates θ, where θ u and θ ω are constant for all simulations of the model, discussed in what follows.

Model calibration and parameter estimation
In this study, we provide a nominal set of approximate model parameters (see Table B.1 in Appendix B), estimated via direct comparison to the available experimental data. An exhaustive analysis in which we explore and fine tune the current model and its parameters will be presented elsewhere. For the calibration of the individualistic model and associated parameters described in Sect. 2.1, we refer the reader to our previous work [34]. One key difference is that the zebrafish pairs observed in this study are found to have an average speed of around 3 BL s −1 , compared with 4.6 BL s −1 reported for isolated zebrafish. For this reason, we select equilibrium speed parameter μ u = 3 BL s −1 as the present nominal value. The remaining free parameters concern those of the interaction response functions: gain strengths K s (speeding attraction / repulsion), K p (turning attraction / repulsion), and K v (turning alignment); repulsion radii in each axis: r u and r ω ; and distance decay function parameters d c and δ.
The plots in Fig. 1b show axial projections of the speeding and turning acceleration maps through the origin in the front-back (cos θ ij = 1) and left-right (sin θ ij = 1) axis respectively, comparing both experimental and simulated data sets. The response curves through these primary axes demonstrate the range of the repulsive region, where we estimate the radii r u ≈ 1.4 BL and r ω ≈ 0.8 BL from the zero-crossings either side of the y-axis. Beyond these radii we observe approximately linearly increasing accelerations with distance. Data suggests that both acceleration responses subsequently decay to zero beyond ≈ 10 BL, however the decreasing density of data at this range hinders an accurate estimation of the cut-off distance d c . For the purposes of this study, we set d c = δ = 10 BL.
Lastly, we fit the gain parameters K s , K p , and K v , by performing simulations of the model with durations up to that of the complete experimental data set (360 min with sampling frequency of 30 Hz) such that sufficient comparison can be made using the same force mapping technique. The results of such a simulation and comparison to experimental data are shown in Fig. 1, where we select K s = 2.5 s −2 and K p = 3 cm −1 s −2 , to best match the response projections for both a s i and a t i . Importantly, we find that the choice of gain parameters for each of the interaction terms is largely independent and can be tuned separately to match the desired amplitudes whilst obtaining a good degree of fit in the repulsive regions. Estimation of the parameter K v is less straightforward, due to the dependence of the alignment term in (4b) on the relative orientation φ ij , which is not represented in the force maps we have calculated (Fig. 1). In this study, we select an approximate alignment strength after fitting K s and K p as described above, tuning the value of K v such that resulting group level dynamics (e.g. distributions of P , M , and C) correspond well with experimental data. For 2-fish simulations, we find K v ≈ 8 s −2 provided the best fit to collective dynamics described by the global observables, with an increase to K v ≈ 15 required to maintain comparable dynamics for N ≥ 5 fish.

Collective dynamics and the role of speed regulation in leadership emergence
Next, we use the model derived above to investigate the role of speed modulation in the emergence of leadership within small groups of fish. We consider the collective dynamics of a group of N naïve agents which are perturbed by a single "informed" agent (IA), moving along a predetermined trajectory, modulating its forward speed through a feedback mechanism with its neighbours.
The model parameters of the agents belonging to the naïve group are presented in Table. B.1 (Appendix B). The internal dynamics of the informed agent (labelled agent 0) is entirely deterministic, with (constant) zero turning speed such that it translates linearly across the simulation cell. Its motion is driven by the same equations in (3), with speed U 0 (t), but where Ω 0 (t) = 0 and σ u dW (t) = 0, ∀t. A compatible interaction network between all agents is chosen based on a periodic equivalent to the non-metric Voronoi neighbourhood used in related models [8,13,33,45]. A more detailed description and a diagram of the construction of this topology is shown in Fig. C.1. The set of system parameters we explore is restricted to independent variations of the equilibrium speeds of the naïve group μ n u and of the IA μ 0 u , along with the degree to which the IA moderates its speed with respect to its naïve neighbours via K 0 s . Interaction gain parameters K [s,p,v] for the naïve group are fixed to nominal values throughout the study.

Quantifying collective behaviour and responses
In order to quantify the effects of the IA, we consider global observables both with respect to the centre-of-mass (CoM) frame of the naïve group and with respect to the IA position [28,31,33].
In the CoM frame, we focus on two observables (order parameters) varying in the range [0-1]: (i) Polarisation P which describes the relative alignment of the group, maximised when every member is oriented in the same direction, and (ii) Cohesion  C which is a measure of the spread of positions around the CoM, as exponentially decreasing function relative to a fixed scale length r coh = 3 BL. Following the description in [33], we also compute associated susceptibilities χ P and χ C of the naïve group polarisation and cohesion -a normalised measure of the respective variances, which we use to describe the relative stability of the group dynamics.
Using the same definitions for the order parameters above, we define a similar set of observables in the IA frame: (i) Polarisation P 0 describes the alignment of the group average orientation with that of the IA, and (ii) Cohesion C 0 from the IA position at x 0 . As described in [28], we also compute a 'stretching' order parameter S 0 , maximised when the velocities of all naïve group members points in the direction of the IA's position.
A full description and mathematical derivation of each of these observables can be found in Appendix A.

Numerical experiments
For this study we consider a group of N = 5 naïve fish, interacting with a single IA (Fig. 2). Average values of the global observables are compared across trial simulations in order to assess the perturbing effects of the IA. Simulations are performed with a single set of initial conditions per trial in a periodic cell of side length L = 40 BL (120 cm) for 1000 s, discarding the first 100 s to remove initial transient dynamics.
The multi-agent model is integrated using a Euler-Maruyama method with a timestep Δt = 0.1 s. Naïve population initial positions, x n , are randomly generated within a circle around the origin, such that their mean separation is of the order 1 BL, with randomly assigned heading directions and initial speeds all equal to μ u . The IA's initial position x 0 is set at the origin for every trial, with heading (fixed) along the positive x-axis with initial speed equal to μ 0 u . We also compute the distributions of the (CoM frame) observables over the same time interval as a function of equilibrium speed μ u , in the absence of the IA. This data, averaged over 10 simulated replicates, is used to provide baseline values to compare the effects of the IA with the unperturbed dynamics of the naïve group (Fig. 3).
Data resulting from numerical experiments associated with the IA perturbations consist of a matrix for each global observable we describe in Sect. 3.1, where each element is the mean value computed over the last 900 s for the combination of parameters [μ u , μ 0 u , K 0 s ]. Each of the resulting matrices are subsequently processed using an algorithm described in [46], providing multi-dimensional smoothing based on a penalised least-squares method via the discrete cosine transform. . We focus on cross sections of the parameter space in which we keep the value of the naïve group equilibrium speed μ u fixed at a value close to the observed average swimming speed (3 BL s −1 ), and at twice this value (6 BL s −1 ). We report the temporal average values of the order parameters and spatial observables, as a function of the ratio of equilibrium speeds μ 0 u /μ u and compare the resulting response curves for different values of the IA (speeding) interaction strength K 0 s . We discuss separately the special case in which the IA moves with fixed speed (no modulation: K 0 s = 0), followed by a description of the effects of increasing K 0 s such that the IA modulates its speed as a function of that of its neighbours.
Videos of the simulated collective dynamics for both the fixed speed IA (μ u = 2μ u , K 0 s = 0) and speed modulating IA (μ u = 2μ u , K 0 s = 4) scenarios are available in the on-line supplementary material.

Informed agent with constant speed (K
The results shown in Fig. 4 report naïve group responses in terms of its CoM polarisation (P ), cohesion (C), and associated susceptibilities (χ P and χ C ). In order to highlight the resulting dynamical perturbations, we display on the vertical axis, the difference between averages of the perturbed system and baseline values (P , C , etc.) computed for the unperturbed five fish school for corresponding values of μ u . Note that when the IA moves at a constant speed, the response characteristics of the different observables shown in each panel are indicated by the darkest trace, corresponding to K 0 s = 0. Our first observation is that when the equilibrium speed of the IA, μ 0 u , is less than that of the naïve group (shaded regions where μ 0 u /μ u < 1), the presence of the IA

3352
The European Physical Journal Special Topics has the effect of disrupting the dynamics of the naïve group, whereby we find reduced polarisation and cohesion compared with the baseline values indicated. Within this region, we also find the corresponding susceptibility measures χ p and χ C to be above their baseline values, indicating an increase in the temporal variance of both order parameters.
As the equilibrium speed of the IA approaches that of the naïve group, the response of each of the observables becomes more pronounced -we consistently observe values increasing (or decreasing) towards the respective baseline. Each response curve then reaches a maximum (or minimum) at a critical value of the ratio μ 0 u /μ u > 1, which is a function of the naïve equilibrium speed μ u . In terms of the group polarisation, we find that baseline values of P and χ P are obtained almost exactly when μ 0 u = μ u , for both examples of fixed naïve group equilibrium speeds shown. This suggests that when the IA is translating linearly at a constant speed, provided the speed is well matched with the naïve group, we recover baseline average polarisation without any decrease in the stability of the group dynamical state.
When the equilibrium speed of the IA is increased beyond μ u , the polarisation continues to rise, exceeding baseline values, accompanied by a proportional decrease in susceptibility χ P . Maximum variation ΔP and Δχ P from the unperturbed baselines are reached simultaneously at approximately 1.4μ u and 1.2μ u for fixed naïve group equilibrium speeds of 3 BL s −1 and 6 BL s −1 respectively, after which average values once again fall below baselines at approximately 2.0μ u and 1.5μ u . Past this value of the IA speed, the dynamics effectively becomes invariant to further increases in μ 0 u . Importantly, we find that the position of this peak in the response curve is almost identical across each of the observables, in terms of the IA equilibrium speed required to produce the maximum response.
A similar analysis of the naïve group cohesion reveals a more subtle group response to the motion of the informed agent, which is found to be enhanced as a function of the equilibrium speed of the naïve population. In the lower speed group (Fig. 4a), we find that an informed agent moving with constant speed always disrupts cohesion with respect to the baseline values (C < C , χ C > χ C ). However, for the faster group (Fig. 4b), the baseline cohesion susceptibility χ C is attained when μ 0 u ≈ μ u , achieving maximal cohesion fractionally above the unperturbed value when μ 0 u = 1.2μ u , somewhat less than the IA equilibrium speed required for maximal polarisation. Further analysis of our data indicates that when perturbed by a fixed speed IA, baseline cohesion can only be exceeded when μ u ≥ 3.6 BL s −1 .
We now consider the response of the naïve group to the relative position and orientation of the IA itself. The IA frame observables (P 0 , C 0 , and S 0 ) characterise the perturbations of the IA specifically in terms of its ability to entrain the naïve group. We define optimal entrainment when P 0 , C 0 , and S 0 are maximised, i.e. when the naïve group: (i) is highly aligned with the IA orientation, (ii) has a compact distribution around the IA position, and (iii) on average has a group velocity directed towards the IA.
The response curves for the IA frame observables are shown in Fig. 5, where again we focus initially on the K 0 s = 0 (darkest) trace. The features of all curves are almost identical, exhibiting peaks in the same position in terms of the IA equilibrium speed, similar to the CoM frame observables. In the IA frame however, we find that as the IA speed is increased from zero, the values of P 0 , C 0 , and S 0 increase monotonically towards their peak values, attained when the (constant) speed of the IA is just faster than the naïve group equilibrium μ u . Above 2.0μ u and 1.5μ u , respectively for the slower and faster naïve group, values decrease and ultimately saturate with further increases in IA speed. Although both P 0 and S 0 saturate to a value just above those found when μ 0 u < μ u , C 0 is further reduced. Most likely, this is due to the periodic translation of the IA, which, moving at speeds too high for the naive group to maintain stable cohesion with its position, is still able to influence the orientation of the group as it passes on each transit.

Feedback modulation of informed agent speed (K
We proceed by considering the effect of an IA which is able to modulate its speed via its interaction with the naïve group. For each of the observables presented in Figs. 4 and 5, the effects of increasing the strength of IA speed modulation are represented by overlaying the set of traces corresponding to discrete increases in K 0 s . As K 0 s is increased, the IA becomes more reactive, now able to accelerate from its equilibrium speed according to (4a), depending on the position of its neighbours.
The family of shifting curves for each observable indicates a clear trend, characterised by a broadening of the peak response in the horizontal axis, corresponding to an increase in IA equilibrium speed relative to the naïve group. In the (shaded) regions where μ 0 u < μ u , we find that the value of K 0 s has a negligible effect on the perturbed dynamics which cannot be actively controlled by varying the parameters of the IA. However, when μ 0 u > μ u , we note that the response peak varies as a function of K 0 s , resulting in enhanced naïve group dynamics (increased P, C) with increased stability (reduced variances χ P and χ C ). Specifically, by comparing the response curves obtained when the IA moves at constant speed, and those where K 0 s > 0, we observe the following features: (i) an approximately linear increase of the equilibrium speed of the IA μ 0 u is required to produce a maximal dynamical response (of all observables), as a function of K 0 s , and (ii) the variation of the peak value of the observables increases as a function of K 0 s (decreasing for χ P and χ C ) until above approximately K 0 s > 4 after which the peak enhancement of speed modulation decreases. Plots indicating the peak observable values and corresponding equilibrium speed ratio are provided in Figs. C.2 and C.3 (Appendix C).
The first observation suggests that as we increase the strength of the interaction, and thus the maximum allowable acceleration of the IA, we also need to increase its equilibrium speed μ 0 u in order to maximise the dynamical response of the naïve group. An intuitive explanation can be given by thinking of the coupling between the informed agent and the naive group in terms of a "spring" connecting the IA to the naive group. In this analogy, the parameter K 0 s represents the stiffness of the connecting spring. As this parameter increases, so does the restoring force pulling the IA towards the group and vice versa. Therefore an increase in the equilibrium speed of the IA is needed to counterbalance such a force and make the IA able to pull away from the naive group so as to influence its collective motion. The second observation indicates that the IA with feedback speed modulation is able to better enhance the dynamics of the naïve group, with increased polarisation, cohesion and entrainment, compared with the IA moving at constant speed. The enhancement provided by increasing the IA's ability to modulate its speed is limited however, where we find an optimal value K 0 s ≈ 4, independent of the naive group equilibrium speed. The effect of introducing feedback control of the IA speed is therefore to magnify the responses when compared to those obtained in the presence of an IA moving at constant speed, provided that we balance increasing values of K 0 s with an appropriate value of the IA equilibrium speed. We shall refer to this region as the "controllable" region where varying the IA speed and interaction strength K 0 s allows to vary at will, within a certain range, the value of each of the observables. This controllability however is limited within a range of K 0 s < 4, above which, increasing the interaction strength eventually reduces the enhancement or entrainment.

Conclusions
After deriving an experimentally consistent model of zebrafish shoals which include novel speed regulatory interactions, we investigated the effects of an informed agent on a naïve group. We explored both the case where the informed agent moves at constant speed and where it modulates its speed as a feedback function of the position of its neighbours.
Our findings suggest that the presence of an IA can increase global polarisation, reduce susceptibility, as well as elicit entrainment along its trajectory, particularly when the IA is able to modulate its speed as a function of its neighbours' positions. Specifically, we find that the peak response occurs when the parameters of the IA enable it to regulate its speed, moving just faster than the naïve group. This allows the IA to maintain high connectivity with the group, and take up a position at the front of the group with respect to their average velocity.

Appendix A -Quantifying collective behaviour and responses
Here, a full description is given of the global system observables which are used to describe the (global) macro-scale, or collective dynamics of both experimental and simulated fish schools. In order to quantify the effects of the informed agent (IA) we consider the global observables both with respect to the centre of mass (CoM) frame of the naïve group, and with respect to the IA position.

Dynamical order parameters
In general, we define the centre of mass (CoM) X(t) at time t of the group of N naïve agents as the mean vector of agent positions x i (t). Computation of this value is complicated by the unbounded, toroidal geometry of the periodic environment. However, a suitable algorithm which efficiently yields the correct quantity is described in [47] and implemented here. Positions relative to the CoM given are thus given by r i (t) = x i (t) − X(t). In this CoM reference frame, we then define two order parameters which describe the dynamic collective state of the shoal: (i) Polarisation P (t) describes the degree to which the orientations of the agents are aligned, maximised when agents all share a common heading direction. (ii) Cohesion C(t) provides a measure of the spread of agents from the CoM, as an exponentially function with respect to a fixed scale-length r coh = 3 BL. The expression for each order parameter are shown below where we henceforth omit explicit time dependences for clarity: where v(t) = (1/N) v i (t) is the mean group velocity of the naïve agents. By substituting CoM position X(t) with the position of the (single) informed agent x 0 we obtain relative positions r i0 , and calculate the cohesion with respect to the IA position C 0 (IA reference frame), see [31]. We also define the polarisation of the naïve group with respect to the orientation of the IA P 0 by considering the alignment between the naïve group (mean) velocity vector and that of the IA v 0 . From (A.1), we obtain where v n is the mean velocity vector of the naïve group. An additional "stretching" order parameter S 0 measured in this reference frame, derived previously in [28],

3356
The European Physical Journal Special Topics considers the alignment of naïve group velocities with the vectors −r i0 towards the IA: where S 0 is maximised when the all naïve group agents are converging directly on the IA's position, reducing as their combined velocities becomes increasingly misaligned (stretching away from its position).

Susceptibilities
Following the recent work of Calovi et al. [33], our analysis also includes a measure of the linear response (fluctuations) of measured quantities (e.g. P , M , and C) to small perturbations, quantified by the susceptibility χ, calculated for the polarisation P response as follows: where angle brackets refer to the time average of the corresponding quantity. We refer the reader to [33] in which the measurement of these fluctuations as a function of the model parameters is used to identify transitions between the various dynamical collective states. In the present study where most model parameters remain fixed for the naïve population, we utilise susceptibilities as a way of quantifying the relative stability of the population dynamics as we vary the behaviour of the IA and its influence on the naïve group. Correspondingly, a reduction in the quantity χ measured for one or more of the order parameters reflects an increased stability of the dynamical state. Appendix C -Additional Figures   Fig. C.1. Computing a Voronoi interaction neighbourhood in a periodic plane.

Appendix B -Tables
The connectivity A = {a ij } ∈ R N×N between agents in this model are derived from the first-shell of the Voronoi tessellation of 9N positions, replicated from N positions in the L × L primary image (pink square) across the eight additional virtual images forming the periodic system. The network topology in two-dimensions is first computed via the Delaunay triangulation of all 9N points (blue edges), of which we extract all edges to or from one of the N vertices in the primary image. From this subgraph we remove edges representing self interactions across periodic images and, where vertex pairs are duplicated (primary-virtual or primary-primary), we choose the edge connecting to the closest of the two vertices, resulting in the final interaction network, or adjacency matrix: A (red edges). Finally, the inter-agent distance matrix D = {d ij } ∈ R N×N is calculated from the appropriate distances from the set of 9N vertices, given A.

3358
The European Physical Journal Special Topics