Swarming, Schooling, Milling: Phase diagram of a data-driven fish school model

We determine the basic phase diagram of the fish school model derived from data by Gautrais etal (PLoS Comp. Biol. 8, e1002678 (2012)), exploring its parameter space beyond the parameter values determined experimentally on groups of barred flagtails (Kuhlia mugil}) swimming in a shallow tank. A modified model is studied alongside the original one, in which an additional frontal preference is introduced in the stimulus/response function to account for the angular weighting of interactions. Our study, mostly limited to groups of moderate size (in the order of 100 individuals), focused not only on the transition to schooling induced by increasing the swimming speed, but also on the conditions under which a school can exhibit milling dynamics and the corresponding behavioral transitions. We show the existence of a transition region between milling and schooling, in which the school exhibits multistability and intermittency between schooling and milling for the same combination of individuals parameters. We also show that milling does not occur for arbitrarily large groups, mainly due to a distance dependence interaction of the model and information propagation delays in the school, which cause conflicting reactions for large groups. We finally discuss the biological significance of our findings, especially the dependence of behavioural transitions on social interactions, which were reported by Gautrais etal to be adaptive in the experimental conditions.

individual parameters. We also show that milling does not occur for arbitrarily large groups, mainly due to a distance dependence interaction of the model and information propagation delays in the school, which cause conflicting reactions for large groups. We finally discuss the biological significance of our findings, especially the dependence of behavioural transitions on social interactions, which were reported by Gautrais et al to be adaptive in the experimental conditions. S Online supplementary data available from stacks.iop.org/NJP/16/015026/ mmedia

Introduction
Transitions between different types of collective behaviour play a major role in the adaptiveness of animal groups [1,3] (see [2] 8 ). These transitions are commonly observed when animal groups shift from one collective behaviour to another either spontaneously or in response to a threat. For instance in fish schools, individuals may adopt different spatial patterns when they are travelling, feeding or displaying defensive behaviour [4][5][6][7][8][9]. It is commonly observed that periods of swarming associated with feeding behaviour, where the group remains cohesive without being polarized, are often interspersed with brief periods of schooling during which fish search for food and move from one place to another [4]. The actual mechanisms and behavioural rules that trigger such transitions are still poorly understood.
Several models of collective motion have been introduced in this context. For instance, in the Aoki-Couzin model [10,11], sharp transitions in collective behaviour are observed for small parameter changes. In this zonal model, slight variations in the width of the alignment zone which controls the alignment behaviour of fish to their neighbours can yield drastic changes of the structure and polarization of the school, such as transitions between schooling, milling (collective vortex formation) and swarming. In the physics literature, transitions between different types of collective behaviour are often described as phase transitions. A well studied case is the transition to orientational order observed in the popular Vicsek model [12][13][14]. Similar transitions have been described in self-propelled particle models by varying the value of parameters such as the blind angle [15], individual speed [16] or a 'strategy parameter' [17], which ponders a behavioural compromise between aligning with neighbours and reacting to their direction changes.
All the above models suffer from a lack of experimental validation, having often been developed in a very general context not particularly in relation to experiments/observations. However, a methodology to build models for animal collective motion from the quantitative analysis of trajectories in groups of increasing sizes has been recently proposed [18]. Using video tracking of groups of barred flagtails (Kuhlia mugil) in a shallow tank, the stimulus/response function governing an individual's moving decisions in response to the position and orientation of neighbours was extracted from two-dimensional trajectory data. It was found that a gradual weighting between alignment (dominant at short distances) and attraction (dominant at large distances) best accounted for the data. It was also found that the parameters governing these two interactions depend on the mean speed of the fish, leading to an increase in group polarization with swimming speed, a direct consequence of the predominance of alignment at high speed.
Here we investigate the basic phase diagram of the Gautrais et al (2012) model, exploring its parameter space beyond the parameter values determined from the data. This is done in order to explore the full range of collective behaviour that could be displayed by the model. Our study is limited to groups of moderate size (in the order of 100 individuals) and focuses, beyond the transition to schooling induced by increasing the swimming speed, on the conditions under which a school can exhibit milling dynamics. The non-dimensionalized model reveals that changes in swimming speed are equivalent to changing the two social interaction parameters while maintaining the noise constant. Exploring the parameter plane defined in this way, we find four regions with distinctive collective behaviour. We show in particular that, in this model, as in most others, milling does not occur for arbitrarily large groups.

Model
Let us first recall the findings of Gautrais et al [18,19]. Barred flagtails (K. mugil) in a circular 4 m diameter tank with a depth of 1.2 m were video-recorded from above for a few minutes in groups of one to 30 individuals. In this quasi-two-dimensional geometry, very few 'crossing' events (with one fish passing under another) were observed. At least five replicates were performed for each group size. In each experiment, fish synchronized their (mean) swimming speed, and their instantaneous speed (tangential velocity), although fluctuating by approximately 10-20% around its mean, was found uncorrelated to their turning speed (angular velocity). This led to the representation of single fish behaviour as a two-dimensional smooth random walk in which an Ornstein-Uhlenbeck process [20] acts on the fish angular velocity ω. The angular velocity is mathematically described by the stochastic differential equation where v is the (constant) swimming speed, dW refers to a standard Wiener process, ξ = 0.024 m is a characteristic persistence length and σ = 28.9 m −1 s −1/2 controls the noise intensity, where both ξ and σ have been estimated from experimental data on barred flagtail trajectories. The large-scale properties of this model have been extensively studied by Degond and Motsch [21,22] In the above equation, ω * (t) is the stimulus/response function based on the individuals response to the proximity of the tank wall-irrelevant in the following as we shall consider groups evolving in an infinite domain-and on the reaction to the neighbouring fish orientation and distance. It was found that for groups of two fish, these social interactions are well described by a linear superposition of alignment and attraction with weights depending on d i j , the distance between the two fish. It was found further that in groups of more than two fish, many-body interactions could be safely approximated by the normalized sum of pair interactions with those individuals forming the first shell of Voronoi neighbours of the focal fish ( figure 1). Mathematically, one has where V i is the Voronoi neighbourhood of fish i containing N i + 1 individuals, φ i j = φ j − φ i is the angle between the orientational headings of fish i and j, and θ i j is the angle between the heading of fish i and the vector linking fish i to fish j (figure 1). The functional forms involving the sine function, chosen for simplicity, were also validated by the data. They ensure that interactions vanish when the fish are already aligned or in front of each other as can be seen in figure 1(c). Note there is a speed dependence on the alignment, and hence a tendency to produce polarized groups at large speeds. A similar positive correlation between group speed and polarity has also been found in groups of giant danios (Devario aequipinnatus) and golden shiners (Notemigonus crysoleucas) [23,24].
Note also that the strength of the attractive positional interaction, scaled by k P , increases linearly with d i j , which is rather unrealistic: if two fish were very far apart, they would turn at extremely large speeds towards each other. Moreover, attraction is of the same amplitude whether fish j is in front or behind fish i, which may also look unrealistic, as one might expect fish i to be less 'aware' and attracted to fish j if j is behind i. Nevertheless this form of attraction interaction was found to account well for the data, because in the experimental tank fish never go far away from each other.
However, a careful analysis of the data in pairs of two fish revealed that the behavioural reactions of a focal fish depend on the angular position of its neighbour (θ i j ) even if the resulting impact on collective motion in larger groups cannot be detected in a confined environment. So we will consider the possibility of an angular weighting of the interactions in the form of an extra multiplicative factor to the interaction term ω * (t), so that (1) is rewritten as The function i j was designed to ensure a frontal preference and some kind of rear blind angle. Its specific form (which can be seen in figure 1(d)) was chosen for its simplicity and its unit angular average. In the following, we perform numerical simulations of the model, varying its key parameters: the swimming speed v, and k P and k V , the behavioural parameters controlling attraction and alignment. Most of our study is restricted to moderate-size groups of N = 100 fish. However, we also present results as a function of group size when investigating the robustness of milling (section 6). Typically, a set of 100 random initial conditions is simulated over a period of 1000 s, where the first half of each series was discarded to account for transient behaviour.
The polarization of the school is quantified by the global polar order parameter which reaches order 1 values for strongly polarized, schooling groups. Milling behaviour is detected via the global normalized angular momentum which will tend towards unity when the school is rotating in a single vortex. As we shall see, the two scalar, positive, quantities P and M often play symmetric roles in the ordered phases of the groups: a schooling group will show large P and negligible M, and the opposite is true for a milling group.

Speed-induced transition to schooling
We first report the basic transition from swarming to schooling observed when increasing the swimming speed of individuals. Already observed in [18], it is a natural consequence of the predominance of the alignment interaction at high speeds. Here, we confirm the existence of the transition in the absence of walls confining the group.  (c) Average distance to the neighbours for both sets of simulations. We performed ten replicates for each different speed, and all simulations had 10 7 iterations, where we discarded the first half of the data, and used only one datapoint for every ten iterations.
We find that the transition to schooling occurs at a speed range of zero to four body lengths per second, a reasonable one considering the barred flagtail swimming ability ( figure 2(a)). The transition is smoother with angular preference (red curves) and consequently requires larger speed values to reach similar results as the original model. Even though no pronounced milling is found in either cases of figure 2(b), a more pronounced result is achieved with the angular preference, which seems to directly influence the school distance to the centroid in figure 2(c).

Non-dimensionalization of the model
In order to get more insight into how speed affects the relative importance of noise, directional and positional behavioural reactions, we rewrite the model in a non-dimensionalized form. After a straightforward manipulation we are left with only three independent parameters: The time step dt = dt/(ξ 2 σ 2 ) has no influence on the simulations provided it is small enough. For all simulations we used dt = 0.1α for α < 0.1 or dt = 0.01 for the remaining cases. The full set of non-dimensionalized equations reads where dW refers to a standard Wiener process. The termω * can be rewritten as Note that in our equivalent system, time is now counted in units of a typical timet = v(ξ σ ) −2 , defined with the correlation length ξ and noise amplitude σ . The distance is expressed in units of a typical lengthx = v(ξ σ ) −2 , which is simply the displacement of a fish during timet with speed/velocity v. The angular speed and thus the evolution of the fish orientation are determined by the relative weight of equivalent alignment (β) and positional (γ ) intensity over the equivalent angular inertia term α.
The expressions of the non-dimensionalized coefficients show that increasing speed, while keeping other parameters constant, leads to an increase in the equivalent positional and alignment interactions. At the same time it reduces the angular inertia while keeping the equivalent noise amplitude/intensity constant. The speed induced transition can thus be seen as a competition between noise and social interactions.
Hereafter we perform an extensive study of the non-dimensionalized version of the model, with particular attention paid to the emergence of significant milling.

Phase diagram in the space of behavioural parameters
We explored the parameter plane formed by coefficients β and γ in the domain Two cases are presented, without and with the angular weighting function (left and right panels respectively). In both cases, a transition to schooling is observed as β is increased  (1), while the right column graphs represent simulations with angular preference as seen in (4). The statistics here represented are: (a) and (b) polarization parameter; (c) and (d) milling parameter; (e) and (f) average neighbour distance. Each datum point is the result of an average over 100 different simulations, where the first half of each simulation was discarded to avoid transitional behaviour, resulting in a total of 2×10 6 data points for every parameter configuration.
(figures 3(a) and (b)). Such a transition is virtually independent of γ considering the original model ( figure 3(a)) while presenting a functional form for simulations using the angular preference (figure 3).
A strong qualitative difference is nevertheless observed between the left and the right panels in figures 3(c) and (d). With the angular dependence (right panels), outright schooling is preceded by a large region where milling is observed (3(d)), whereas no significant milling occurs without angular weighting (3(c)). This may not seem too surprising given that favours the formation of (local) columns in which fish follow those in front of them.
Lastly, we observe a peculiar property in figures 3(e) and (f), which show the average distance to neighbours. With the angular preference ( figure 3(f)), the region of low β exhibits a weak ridge bordering the swarming region (γ and β ≈ 0) where the average distance is expected to diverge. This ridge extends to high γ values for which one would infer small-size schools and shorter neighbour distances. One should be aware that the colour codes presented on these last two figures are presented on a logarithmic scale. This was introduced since a divergence occurs in the swarming region, which would mask secondary ridges, such as the one found in figure 3(f).
The explanation, despite being initially counterintuitive, is rather simple. Despite the γ -induced tendency to stay together, a low β means that fish have a lower tendency to disturb their alignment to match that of their neighbours, meaning that fish will rely mostly on the positional interaction. This implies that if a neighbour is directly ahead or behind the focal fish, there is no need to change the direction of motion. Also, both cases (with or without angular preference) have zero interaction for anti-parallel neighbours. This means that when γ dominates the interactions, the stable fish configuration is a line. Having noise in (1) or (4) leads to fish organizing themselves in a quasi-unidimensional configuration, and the fish in front, not having additional neighbours on their sides, will eventually turn due to noise. After a slight turn, this fish will now have new Voronoi neighbours, with fish previously located behind it now located slightly to its side, and make a 180 • turn to recover the optimal configuration. This mechanism leads fish to organize themselves in two columns going in opposite directions connected at the end of the school by this 'U turn' (see snapshot in figure 4(c), panel III).
To visualize the above results in a synthetic way, we took advantage of the fact that the main behavioural regions apparent in figures 3(b), (d) and (f) are distinguished by mutually exclusive quantities: when schooling is strong, milling is weak and the distance to neighbours (D) is small; when milling is strong, schooling is weak and D is small; when D is larger, both schooling and milling are weak. This mutual exclusion of schooling and milling is clear in the time series shown in figure 4(a). We thus constructed the following composite order parameter: which takes complex values such that the phase codes for the behaviour and the modulus for the intensity of this behaviour, where the average distance was normalized according to the largest value observed. As a result, the information of the right panels of figure 3 is represented synthetically in figure 4(b). Three distinct regions are clearly apparent as each is represented in a different colour (red for polarization, blue for milling and green for average neighbour distance).  regions II and III are similar, having only changed the width/length ratio. Nevertheless, region III has very weak rotational behaviour at odds with full blown milling. Note finally the intermittence between schooling and milling in the transition region between regions I and II.
Finally, we come back to the effect of the swimming speed. The above analysis in the (β, γ ) plane was repeated at different values of α corresponding to speed values of 0.4, 0.8 and 1.2 m s −1 , building a three-dimensional phase diagram. In figure 5 we show that these transition lines superimpose onto a master curve delimiting the schooling and milling zones. Thus, no qualitative changes were observed for these values of α, indicating the robustness of the phase diagram in figure 4 for speeds in the experimental range. We studied this phase boundary marking the transition, via intermittent regimes, between schooling and milling more quantitatively. Defining the transition points as those where both the polarization and milling parameters spend more than 40% of the time above the threshold value of 0.8, we find that the where A and B take values independent of α in the ranges studied ( figure 5).
Although we do not have a quantitative explanation for the above functional form, these fits indicate the possibility of a rather simple theory accounting for the full phase diagram of our model. Meanwhile, a qualitative explanation for the dependence on γ is given by the need to counteract the global polarization tendency given by β, while maintaining a local one. Also, figure 3(a) indicates minimum β × γ values to overcome the noise. In the supplementary material (available from stacks.iop.org/NJP/16/015026/mmedia) there are two videos (movies 5 and 6) of the two dimensional histogram evolution for the schooling and milling parameters as we change β or γ while maintaining the other one constant (γ = 17.28 and β = 13.30), meaning we have respectively vertical and horizontal transition cross-sections in figure 5.

Group-size-induced transition
In this section, we investigate the influence of group size on the robustness of both milling and schooling behaviour in our model using the angular weighting function. Apart from the general theoretical context mentioned above, this is also of direct relevance for fish, as it has been reported recently that in golden shiners increasing group size significantly increases the amount of time spent in a milling state [24]. We performed a series of simulations at fixed γ ≈ 4.5 with an angular inertia term α ≈ 0.014 (equivalent to a speed v = 0.8 m s −1 ). Figure 6 shows the milling and schooling order parameters for simulations of groups from 10 to 4000 fish. These results indicate that there exists an inferior limit of approximately 60 fish in order to achieve significant milling and that this behaviour progressively disappears at large sizes. Although these results must be taken with care given the unbounded character of the variation of positional interaction with distance to neighbours, they show that for up to 100 fish, increasing group size induces a transition from schooling to milling for β 8. For larger values of the alignment parameter, despite having higher overall values for the milling parameter, the transition does not occur. This indicates that in the model studied here, like in many of those proposed before [25], milling dynamics does not emerge in the arbitrarily large infinite-size limit.
Although such a statement has obvious interest for physicists, we believe it bears some importance regarding animal group behaviour. Also, it is worth noting that the variation in the milling parameter seen for β 12 in figure 6(a) happens as the duration of the milling and schooling behaviour increases with N . Such an increase is so intense that for simulations with more than 300 fish, usually only one type of behaviour can be seen for every initial condition, giving a very large standard deviation and a highly fluctuating average. Furthermore, some simulations for N > 1000 displayed more than one vortex at the same time, patterns which the milling parameter cannot account for. One can see in figure 6(b) that the schooling behaviour is affected by the size as well. This again is due to the effect of large schools. This stretch on school extensions enables information propagation delays, which in turn, cause reductions on the global polarization parameter. It is easy to see how a three dimensional configuration could achieve lower extensions with large fish quantities, which would in turn minimize these effects.

Conclusions
Understanding how complex motion patterns in fish schools arise from local interactions among individuals is a key question in the study of collective behaviour [26,27]. In a previous work, Gautrais et al have determined the stimulus/response function that governs an individual's moving decisions in barred flagtails (K. mugil) [18]. It has been shown that two kinds of interactions controlling the attraction and the alignment of fish are involved and that they are weighted continuously depending on the position and orientation of the neighbouring fish. It has also been found that the magnitude of these interactions changes as a function of the swimming speed of fish and the group size. The consequence being that groups of fish adopt different shapes and motions: group polarization increases with swimming speed while it decreases as group size increases.
Here we have shown that the relative weights of the attraction and alignment interactions play a key role in the emergent collective states at the school level. Depending on the magnitude of the attraction and the alignment of fish to their neighbours, different collective states can be reached by the school. The exploration of the parameter space of the Gautrais et al model reveals the existence of two dynamically stable collective states: a swarming state in which individuals aggregate without cohesion, with a low level of polarization, and a schooling state in which individuals are aligned with each other and with a high level of polarization. The transition between the two states is induced by an increase of the swimming speed. Furthermore, the addition in the model of a frontal preference to account for the angular weighting of interactions leads to two other collective states: a milling state in which individuals constantly rotate around an empty core thus creating a torus and a winding state, in which the group self-organizes into a linear crawling structure. This last group structure is reminiscent of some moving patterns observed in the Atlantic herring (Clupea harengus) [28].
Of particular importance is the transition region between milling and schooling. In this region, the school exhibits multistability and it regularly shifts from schooling to milling for the same combination of individual parameters. This particular property was recently reported in experiments performed in groups of golden shiners [24]. Our results show that the transition region can be described by a simple functional form describing the respective weights of the alignment and attraction parameters and is independent of the fish swimming speed in the experimental range. The modulation of the strength of the alignment and attraction may depend on the behavioural and physiological state of fish [29]. In particular, various environmental factors such as a perceived threat may change the way fish respond to their neighbours and hence lead to dramatic changes in collective motion patterns at the school level [30,31].
Finally, we show that collective motion may dramatically change as group size increases. The absence of a milling state in the largest groups is a natural consequence of the spatial constraints exerted upon individuals' movements as the number of fish exceeds some critical value. These constraints could be much less stringent in three dimensions, as testified by the observation of cylindrical vertical milling structures involving several thousands of fish in bigeye trevally (Caranx sexfasciatus). With an additional dimension, the same number of fish could result in a mill of much smaller diameter, not only avoiding the problems of distance dependence, but also minimizing information propagation delays and the emergence of competing behaviour.
A more in depth physics study of the transitions presented here is necessary. Unfortunately, the present model is restricted to schools of biologically relevant sizes (N ≈ 100 fish) preventing such analysis for the current social interactions. In this manner, additional changes to the model are required, among them are: (i) eliminating the linear distance dependence on the interactions, by either implementing a saturation or a decay to this interaction; (ii) changing the boundary conditions for periodic ones to avoid evaporating schools once the distance dependence is removed; and (iii) a three dimensional study of the model to check the impact of an additional level of freedom on the school behaviour.