# Scale-free behavioral cascades and effective leadership in schooling fish

Julia Múgica  Departament de Física, Universitat Politècnica de Catalunya,
Campus Nord B4, 08034 Barcelona, Spain  Jordi Torrents  Departament de Física,
Universitat Politècnica de Catalunya, Campus Nord B4, 08034 Barcelona, Spain
Departament de Física de la Matèria Condensada, Universitat de Barcelona,
Martí i Franquès 1, 08028 Barcelona, Spain  Javier Cristín  Grup de Física
Estadística, Departament de Física, Facultat de Ciències, Universitat Autònoma
de Barcelona, 08193 Barcelona, Spain  Andreu Puy  Departament de Física,
Universitat Politècnica de Catalunya, Campus Nord B4, 08034 Barcelona, Spain
M. Carmen Miguel  Departament de Física de la Matèria Condensada, Universitat
de Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain  Universitat de
Barcelona Institute of Complex Systems (UBICS), Universitat de Barcelona,
Barcelona, Spain  Romualdo Pastor-Satorras  Departament de Física, Universitat
Politècnica de Catalunya, Campus Nord B4, 08034 Barcelona, Spain

###### Abstract

Behavioral contagion and the presence of behavioral cascades are natural
features in groups of animals showing collective motion, such as schooling
fish or grazing herbivores. Here we study empirical behavioral cascades
observed in fish schools defined as avalanches of consecutive large changes in
the heading direction of the trajectory of fish. In terms of a minimum turning
angle introduced to define a large change, avalanches are characterized by
distributions of size and duration showing scale-free signatures, reminiscent
of self-organized critical behavior. We observe that avalanches are generally
triggered by a small number of fish, which act as effective leaders that
induce large rearrangements of the group’s trajectory. This observation
motivates the proposal of a simple model, based in the classical Vicsek model
of collective motion, in which a given individual acts as a leader subject to
random heading reorientations. The model reproduces qualitatively the
empirical avalanche behavior observed in real schools, and hints towards a
connection between effective leadership and avalanche behavior in collective
movement.

##  I Introduction

Collective motion is an ubiquitous phenomenon in nature, observed in a wide
variety of different living systems and on an even wider range of scales, from
mammal herds and fish schools, to bacteria colonies and cellular migrations
Vicsek and Zafeiris (2012); Sumpter (2010); Ramaswamy (2010). The study of
collective motion allows scientists to infer the intricate interaction
mechanisms governing the diversity of behaviors found in natural grouping
species Vicsek and Zafeiris (2012); Rosenthal  _et al._ (2015); Chen  _et al._
(2016); Calovi  _et al._ (2018). Identifying the most relevant traits will
prove essential if we ever want to take advantage of nature wisdom for
engineering applications such as in swarm robotics Brambilla  _et al._ (2013)
or in driver-less cars. Social animals group and travel together to gain
several benefits, from better foraging and more efficient offspring training,
to improved navigational accuracy and reduced risk of predation Krause and
Ruxton (2002). Examples illustrating the emergence of ordered collective
motion in social animal groups can take the spectacular form of wildebeest
herds crossing deserts in Africa, or huge fish schools running away
coordinately from predators Procaccini  _et al._ (2011). From a more mundane
perspective, the seemingly simple movement of a sheep herd crossing a road
also arises as a result of the collective, coordinated motion of individual
sheep Ginelli  _et al._ (2015).

The kind of coordination required to produce such impressive patterns of
collective motion requires an efficient transfer of information among the
group components. In this regard, leadership is sometimes brought about to
rationalize the cooperative movements of animal groups by single individuals
that appear to have a strong influence on the flock behavior Couzin  _et al._
(2005). The effects of leadership have been considered in several contexts,
including crowd behavior Aubé and Shield (2004), hierarchical leadership Shen
(2008), linear response theory in flocking systems Pearce and Giomi (2016);
Kyriakopoulos  _et al._ (2016) or the emergence of complex patterns of
cooperation and conflict Smith  _et al._ (2016). Leadership can arise as a
natural instinct in some animals, which form a permanent hierarchical
structure, but in other cases it can exhibit a switching dynamics that can
even depend on context Nagy  _et al._ (2010); Flack  _et al._ (2012); Nagy
_et al._ (2013); Chen  _et al._ (2016). In this sense, effective leadership
can come from individuals having useful information about their environment,
such as the position of food or predators, not visible to the rest of the
flock Couzin  _et al._ (2011); Ward  _et al._ (2011). Another important aspect
of collective animal motion is the existence of spontaneous individual-level
behavioral variations, which may be transmitted to the group as if those
particular individuals were effective group leaders. As a result of abrupt
changes of dynamic behavior at the level of one or a few individuals, animal
groups can exhibit intermittent collective rearrangements, or can even undergo
state transitions at the macroscopic level. Collective behavioral oscillations
or waves in groups have been reported in golden shiners Rosenthal  _et al._
(2015), which were related to an underlying or hidden communication network
Rosenthal  _et al._ (2015). On the other hand, sheep herds have been shown to
pass from slow group dispersive motion while grazing, to rapid aggregation
induced by sudden, individual changes of speed Ginelli  _et al._ (2015) in the
absence of nearby threatening sources. Most interestingly, these experimental
studies emphasize that animal rearrangements can either spread extensively
within the group or extinguish rapidly, leading to an avalanche-like type of
response with a broad-tailed distribution of avalanche magnitudes Rosenthal
_et al._ (2015); Ginelli  _et al._ (2015). This sort of avalanche behavior is
well known in the physics literature Fisher (1998), where it has been
discussed in magnetic materials Zapperi  _et al._ (1998), superconductors
Altshuler and Johansen (2004), plastic deformation of crystalline materials
Miguel  _et al._ (2001), fracture phenomena Zapperi  _et al._ (1999), or
earthquakes Kawamura  _et al._ (2012).

In this paper we examine the interplay between effective leadership and
behavioral cascades (avalanching behavior) by means of an empirical analysis
of the movement of black neon tetra fish _Hyphessobrycon herbertaxelrodi_ ,
and through the theoretical analysis of a variation of the classical Vicsek
model Vicsek  _et al._ (1995) that includes an explicit leader. In our
empirical analysis, we define avalanches in terms of changes in the fish
heading above a given turning angle threshold, which lead to a sudden
reorientation of the global trajectory of the school. We observe that the
distributions of size and duration of the measured avalanches show scale-free
signatures in analogy with self-organized critical processes Pruessner (2012)
that can be described in terms of a set of characteristics scaling exponents.
We explore the possible presence of leadership by considering the statistics
of avalanche initiators, observing that some fish have an anomalous large
probability of starting an avalanche, acting thus as effective leaders
promoting substantial school rearrangements. In order to check the general
effects of leadership in avalanche behavior, we consider a Vicsek-like model
in which a global leader alternates a directed motion, unaffected by other
individuals, with sudden variations of its direction of motion, in the spirit
of run-and-tumble locomotion Méndez  _et al._ (2014). Akin to our experimental
observations, the model exhibits intermittent scale-free avalanche-like
behavior, not present in the original model. Our results confirm the presence
of scale-free signatures in behavioral cascades in collective motion Rosenthal
_et al._ (2015) and highlight the role of effective leadership interactions in
the emergence of this sort of collective behavior.

##  II Empirical analysis of schooling fish

We have analyzed the avalanche behavior in black neon tetra (_Hyphessobrycon
herbertaxelrodi_), a small freshwater fish (adult mean body size of $2.5$cm)
that have a strong tendency to form compact and highly polarized schools
Gimeno  _et al._ (2016). Experiments, performed by the group of F. S. Beltrán
and V. Quera at the Institute of Neurosciences, University of Barcelona
(Spain), consisted in groups of $40$ individuals, freely swimming in an
experimental rectangular tank of dimensions $100\times 93$cm and $5$cm of
depth. Videos of the fish movement were recorded at $20$ frames per second
with a resolution of $1072\times 1004$ pixels. Three independent recordings,
each of length $T=12000$ frames were performed. The path of individual fish
was digitized using a custom-made tracking software. The paths obtained were
later visually inspected to correct for a few anomalously large jumps, due to
the switching of fish identities in the tracking process. The trajectories
obtained were finally smoothed by applying a Savitzky-Golay filter of window
length $7$ and polynomial order $3$ Press and Teukolsky (1990). The final
result were three digital data series, labeled A, B and C, that we analyzed
numerically.

###  II.1 Avalanche analysis

The heading of each fish is defined in terms of its instantaneous velocity
$\vec{v}_{i}(t)$. Given the path of a fish as a function to time
$\vec{r}_{i}(t)$, $t=1,2,\ldots,T$ (time measured in frames), we define the
velocity at time $t$ using a Richardson extrapolation of order $4$, in terms
of the expression Fornberg (1988)

$\vec{v}_{i}(t)=\frac{1}{12}\left[\vec{r}_{i}(t-2)-8\vec{r}_{i}(t-1)+8\vec{r}_{i}(t+1)-\vec{r}_{i}(t+2)\right],$
(1)

for $t=3,4,\ldots,T-2$.

In our recordings, fish tend to move in a coherent and highly polarized
fashion, swimming along a common and slowly changing average velocity. In this
regime, the heading of the fish velocities changes almost continuously.
However, at some instants of time, we can recognize swift rearrangements of
the individuals’ headings, that lead to a change of the average orientation of
the school, accompanied by an increase of the average velocity and a decrease
and a delayed increase of the global order of the school. We interpret these
sudden rearrangements of individual heading as triggeres of _avalanches_ of
activity. In order to quantify them, we fist examine the angular variations in
the heading of individual fish, defined as the _turning angle_
$\varphi_{i}(t)$ formed by the velocity vectors $\vec{v}_{i}(t+1)$ and
$\vec{v}_{i}(t)$, and computed as

$\varphi_{i}(t)=\left\lvert\arctan\left\\{\frac{\lVert\vec{v}_{i}(t)\times\vec{v}_{i}(t+1)\rVert}{\lVert\vec{v}_{i}(t)\cdot\vec{v}_{i}(t+1)\rVert}\right\\}\right\rvert,$
(2)

where $\times$ stand for the vectorial product and $\lVert\cdot\rVert$
represents the vector modulus. For symmetry reasons, the angles, computed in
the interval $[-\pi,\pi]$, are projected onto $[0,\pi]$.

Figure 1: (a) Probability density of turning angles $\varphi$ of individual
fish. (b) Number of fish turning an angle larger than $\varphi_{\mathrm{th}}$
in a sequence of $1000$ frames in data series A, for different values of the
turning threshold. (c) Number of avalanches observed as a function of the
turning threshold. (d) Relative fluctuations of the avalanche size ($s$) and
duration ($t$) distributions as a function of the turning threshold.

In Fig. 1(a) we plot the distribution of turning angles $P(\varphi)$ for the
different data series (hollow symbols) and for the aggregation of all three of
them (red line). Discarding some extremely small values of $\varphi$, which
can be attributed to imprecision of the tracking algorithm when following an
essentially straight segment, the distributions show an extended plateau for
small turning angles, corresponding to stretches of time in which fish barely
change their heading and are thus compatible with movement along a smoothly
winding trajectory. Instead, for values larger than $0.01$ radians, the
distribution starts to decrease sharply. These large turns correspond to the
sparsely distributed large rearrangements of direction observed in the fish.
In order to quantitatively identify avalanches, as is customarily done in
condensed matter physics Laurson and Alava (2006), we define a _turning
threshold_ $\varphi_{\mathrm{th}}$ that distinguishes small turns
$\varphi<\varphi_{\mathrm{th}}$, associated to smooth trajectories, from large
turns $\varphi>\varphi_{\mathrm{th}}$, associated to sudden rearrangements
that trigger an avalanche. In Fig. 1(b) we plot, for a given value of the
threshold, the number of _active_ fish, defined as those performing a turn
larger than $\varphi_{\mathrm{th}}$, as a function of time. Here we can see
the actual presence of turning avalanches, defined as trains of consecutive
frames in which more than one fish is active, delimited by two frames (one at
the beginning and another at the end of the train) with no active fish. These
curves highlight the intermittent and heterogeneous character of avalanches,
which may be rather small or can also reach relatively large sizes.

Our experimental data only allows the identification of a limited number of
avalanches. Indeed, in Fig. 1(c) we plot the total number of recorded
avalanches as a function of the turning threshold. From here we observe that
the range of values of $\varphi_{\mathrm{th}}$ that lead to at least $1000$
avalanches range approximately in the interval $[0.20,1.50]$. To study the
statistics of avalanches, we compute their duration $t$ and size $s$, defined
as the number of consecutive time steps (frames) with at least one active
fish, and the sum of the number of active fish at each time step of an
avalanche, respectively. A first broad statistical characterization is given
by the relative size and duration fluctuations, measured as the standard
deviation $\sigma$ divided by the corresponding average value $\mu$. In Fig.
1(d), we plot these relative fluctuations for both $s$ and $t$, respectively.
From this plot, we observe that relative fluctuations are only larger than $1$
for threshold values within the interval between $0.1$ and $1.2$. We therefore
restrict our analysis to the conservative threshold interval $[0.20,1.20]$.

Figure 2:  (a) Cumulative probability distribution of the avalanche sizes
$P_{c}(s)$ for different values of the turning threshold
$\varphi_{\mathrm{th}}$. (b) Cumulative probability distribution of the
avalanche durations $P_{c}(t)$ for different values of the turning threshold
$\varphi_{\mathrm{th}}$. Here and in (b), filled symbols represent the actual
empirical distributions, while hollow symbols correspond to distributions
obtained by randomizing the turning angles in the trajectory of each fish. (c)
Check of the scaling of the cumulative size distribution with the turning
threshold, as given by Eq. (5). (d) Check of the scaling of the cumulated time
distribution with the turning threshold, as given by Eq. (5). (e) Average size
$\bar{s}_{t}$ of avalanches of fixed duration $t$ as a function of $t$. The
main plot shows the empirical data. The inset presents the results from a
randomization of the turning angles in each fish trajectory. In this case, we
plot the average duration as a function of the theoretical prediction $Nqt/Q$,
see Eq. (13).

We consider the shape of the probability distributions of avalanche sizes,
$P(s)$, and durations, $P(t)$, focusing on the cumulative distributions,

$P_{c}(s)=\sum_{s^{\prime}=s}^{\infty}P(s^{\prime})\qquad\mbox{and}\qquad
P_{c}(t)=\sum_{t^{\prime}=t}^{\infty}P(t^{\prime}).$ (3)

In Fig. 2(a) we plot in full symbols the cumulative size distribution obtained
for different values of $\varphi_{\mathrm{th}}$. From the double logarithmic
scale in the plot, we can see that the size distributions show long tails,
compatible with a power-law behavior of the form $P(s)\sim s^{-\tau_{s}}$ for
small values of $s$. This power-law behavior is due to the correlated nature
of turns in the fish school, a feature that can be corroborated from the
comparison of these results with the avalanche distributions obtained from
trajectories reconstructed by randomizing the sequence of turning angles of
each fish. In the latter case, one obtains a clear exponential decay, as shown
in hollow symbols in Figs. 2(a) and (b); see Appendix A.

Upon closer scrutiny, we can also observe that, for sufficiently large
$\varphi_{\mathrm{th}}$, the initial power-law behavior of the size
distributions is followed by a faster decay for $s$ larger than a
characteristic size $s_{c}$ that appears to be a decreasing function of the
threshold $\varphi_{\mathrm{th}}$. Inspired by the observations in other
avalanche systems Laurson  _et al._ (2009) and in models of self-organized
criticality Pruessner (2012), we can assume that, for different values of the
threshold, the size distributions exhibit a scaling behavior of the form

$P(s)=s^{-\tau_{s}}G_{s}\left(\frac{s}{s_{c}(\varphi_{\mathrm{th}})}\right)$
(4)

where the scaling function $G_{s}(z)$ is constant for small $z\ll 1$ and
decays rapidly to zero for $z\gg 1$. In analogy with avalanches in condensed
matter and critical phenomena Yeomans (1992); Cardy (1996) we make the ansatz
for the behavior of the size cut-off
$s_{c}(\varphi_{\mathrm{th}})\sim\varphi_{\mathrm{th}}^{-\sigma_{s}}$, where
$\sigma_{s}$ is a characteristic exponent. We can estimate the values of the
exponents by noticing that Eq. (4) implies, for the cumulative distribution,
$P_{c}(s)=s^{-\tau_{s}+1}F_{s}\left(s\varphi_{\mathrm{th}}^{\sigma_{s}}\right)$,
where $F_{s}(z)$ is another scaling function. The previous expression can be
rewritten as

$\varphi_{\mathrm{th}}^{\sigma_{s}(1-\tau_{s})}P_{c}(s)=F^{\prime}_{s}\left(s\varphi_{\mathrm{th}}^{\sigma_{s}}\right),$
(5)

where $F^{\prime}_{s}(z)=z^{-\tau_{s}+1}F_{s}(z)$. Eq. (5) implies that, when
plotting the rescaled distribution
$\varphi_{\mathrm{th}}^{\sigma_{s}(1-\tau_{s})}P_{c}(s)$ as a function of the
rescaled size $s\varphi_{\mathrm{th}}^{\sigma_{s}}$, with the correct
exponents $\tau_{s}$ and $\sigma_{s}$, plots for different values of
$\varphi_{\mathrm{th}}$ should collapse onto the same universal function
$F^{\prime}_{s}(z)$. Using this idea, one can estimate numerically the
exponents $\tau_{s}$ and $\sigma_{s}$ as those that provide the best collapse
of the data rescaled using Eq. (5) for the different values of
$\varphi_{\mathrm{th}}$, see Appendix B.

Following this approach, using values of $\varphi_{\mathrm{th}}$ in the
interval $[0.7,1.2]$, we estimate the exponents $\tau_{s}\simeq 2.03$ and
$\sigma_{s}\simeq 3.03$. In Fig. 2(c) we show the data collapse for Eq. (5)
obtained for the cumulated size distributions using these values. Different
intervals of the turning threshold provide slightly different values of the
exponents, from which we estimate the average exponents quoted in Table 1. The
same procedure can be applied to the duration distribution, see Fig. 2(b),
where now the cumulative duration distribution $P_{c}(t)$ fulfills Eq. (5)
with the corresponding exponents $\tau_{t}$ and $\sigma_{t}$. In the same
interval of thresholds we find $\tau_{t}\simeq 2.33$ and $\sigma_{t}\simeq
1.59$, see Fig. 2(d), while the average exponents are given in Table 1. We can
check the validity of these results considering that, for small values of $s$
and $t$, the distributions $P(s)\sim s^{-\tau_{s}}$ and $P(t)\sim
t^{-\tau_{t}}$ imply that the average size of avalanches of duration $t$,
$\bar{s}_{t}$, takes the form

$\bar{s}_{t}\sim t^{m},\qquad\qquad\mathrm{with}\quad
m=\frac{\tau_{t}-1}{\tau_{s}-1}.$ (6)

In Fig. 2(e) we represent the empirical average avalanche size $\bar{s}_{t}$
as a function of the duration $t$. For the different values of the turning
threshold considered, we estimate numerically that $\overline{s}\sim
t^{1.31}$. This observation is in good agreement with the expression in Eq.
(6), which, using the values from Table 1, yields $m=1.4(2)$. In Fig. 2(e) we
also show the average avalanche size observed in randomized avalanches, which
shows a linear dependence as expected theoretically, see Appendix A. This last
result highlights the relevant effect of turning angle correlations in real
fish.

| Schooling fish
---|---
| $\tau_{s}$ | $\sigma_{s}$ | $\tau_{t}$ | $\sigma_{t}$
| 2.0(1) | 3.1(3) | 2.4(1) | 1.70(4)
|  | Vicsek model with a perturbed leader
---
$\varphi_{\mathrm{th}}(\eta)=2.5\pi\eta$
| $\tau_{s}$ | $D$ | $\tau_{t}$ | $z$
$\eta=0.2$ | 1.73(5) | 2.01(2) | 4.03(5) | 0.39(5)
$\eta=0.3$ | 1.69(5) | 2.01(2) | 3.49(5) | 0.42(5)
|  | Vicsek model with a perturbed leader
---
$\varphi_{\mathrm{th}}(\eta)=2.8\pi\eta$
| $\tau_{s}$ | $D$ | $\tau_{t}$ | $z$
$\eta=0.2$ | 0.99(5) | 2.06(2) | 0.26(4) | 0.50(5)
$\eta=0.3$ | 1.04(5) | 2.03(2) | 0.57(5) | 0.52(5)
Table 1: Summary of scaling exponents for the avalanche size and duration
distributions obtained from observations of a real fish school and from the
Vicsek model with a perturbed leader.

###  II.2 Effective leadership and avalanche behavior

Figure 3:  Leadership probability $\chi_{i}$ computed for each fish $i$ in the
three data series considered A, B and C (rows, left to right) and for
different values of the turning threshold $\varphi_{\mathrm{th}}=0.25,0.5,1.0$
and $1.2$ (columns, top to bottom). Symbols are color-coded with the number
$N_{\alpha}$ of actual avalanches in which each fish participates. Full lines
represent the average leadership probability in a null model of uncorrelated
avalanches. The shadowed regions represents the $99\%$ confidence interval of
this value.

In order to explain the origin of the avalanche behavior observed in our
empirical data, we consider the possibility that avalanches are triggered by
some initiator or effective _leader_ , which consistently starts the large
turning rearrangements that lead to the formation of an avalanche. While
several definitions of leadership have been proposed within the field of
collective animal motion Strandburg-Peshkin  _et al._ (2018), here we use a
measure explicitly devised to detect the presence of preferential initiators
of avalanches. We consider the originator of an avalanche as the fish that
performs the first large heading turn in the evolution of the avalanche. As
more than one fish can be active in any frame, we consider as initiators all
active fish in the first frame of an avalanche. We define the _leadership
probability_ $\chi_{i}$ of fish $i$ in a given data series as the ratio of the
number of avalanches in which the fish $i$ is active in the first frame,
divided by the total number of avalanches in which fish $i$ participates. The
calculation is restricted to sufficiently large avalanches, of duration larger
than $5$ frames. In Fig. 3 we plot the value of $\chi_{i}$ computed for each
one of the $N=40$ fish in each series, for different values of the turning
threshold $\varphi_{\mathrm{th}}$. As we can see, the leadership probability
shows an important variation among fish. Moreover, for the largest values of
$\varphi_{\mathrm{th}}$ considered, the leadership probability can take values
up to $0.60$, indicating that some fish initiate more than half of the
avalanches in which they participate.

In order to quantify the relevance of the values of $\chi_{i}$ obtained, and
ascertain that they are not the effect of random fluctuations in the activity
of the fish, given our small populations, we compare our empirical estimates
with the results obtained in a null model in which the turns performed by fish
are completely independent, see Appendix C. The continuous lines in Fig. 3
represent the null model average leadership probability, while the shaded
region represents its $99\%$ confidence interval. Our results in Fig. 3
indicate that, with the exception of series C, in all series and for all
values of the turning threshold considered, several fish have an unusually
very large probability to initiate an avalanche, much larger than the value
expected from pure random fluctuations. We can associate them to effective
leaders of the school, which initiate with large probability the avalanches in
which they participate.

More information can be obtained by considering the evolution of the
leadership probability as a function of the turning threshold for each fish in
each time series, see Supplementary Figure SF1. From this plot we can confirm,
first of all, that some fish never initiate an avalanche ($\chi_{i}=0$) for
large values of the turning threshold, while others consistently start much
less avalanches than they should by mere random fluctuations. Some other fish
behave as initiators for some range of values of the turning threshold.
Finally, some fish reliably initiate a large number of avalanches, much more
than they should by pure randomness. These fish can be identified as
consistent effective leaders, which trigger a large majority of the
avalanches, independently of the value of the threshold used to quantitatively
define them.

##  III Modeling avalanches in the presence of leaders

To explore the effects of the presence of leadership in the avalanche behavior
of schooling fish, we consider as the simplest modeling scenario a variation
of the classic Vicsek model Vicsek  _et al._ (1995) in which we introduce an
effective leader.

###  III.1 Model definition

The Vicsek model Vicsek  _et al._ (1995); Ginelli (2016) is defined in terms
of $N$ self-propelled particles (SPPs), characterized by a position
$\vec{r}_{i}$ and a velocity $\vec{v}_{i}$ of constant modulus $v_{0}$,
evolving in a two dimensional space, and thus being fully characterized by the
heading angle $\theta_{i}$ defined by the velocity vector. Particles interact
among them by trying to align their instantaneous velocity with the average
velocity of the set of nearest neighbors inside a circular region of radius
$R$ centered in the considered SPP. A noise source of strength $\eta$,
representing physical or cognitive difficulties in gathering or processing
local information, allows the formation of an ordered (_flocking_) phase at
low noise intensities, and of disordered (_swarming_) states at high enough
noise values. See Appendix D for further details and simulation conditions.

In the variation we study, see Appendix D, we consider that a given particle,
say particle $1$, plays the role of an effective global leader with a long
range influence over the rest of the SPPs. The velocity of the leader,
$\vec{v}_{1}(t)=\vec{v}_{L}$ is not affected by the behavior of its neighbors,
and therefore its heading could remain constant $\theta_{1}(t)=\theta_{L}$
over time. The other SPPs can, on the other hand, feel the orienting effect of
their local neighborhood as well as that of the leader, independently of their
relative distance. Therefore they take it into account when computing the
average velocity of their neighbors, to which they try to align.

The presence of an unperturbed leader has the effect of suppressing the
disordered phase exhibited by the classic Vicsek model. As we can see in Fig.
SF2, while for the classical model the transition becomes sharper when
increasing the systems size $L$, the leader induces an ordered state for any
value of $\eta$, with an order parameter (see Appendix D) fairly independent
of system size and vanishing only in the limit of maximum noise $\eta=1$.

###  III.2 Avalanche behavior in response of leader perturbations

In this section, we focus our attention on the system-wide perturbations that
are induced by changes in the preferred direction of motion of the leader. To
analyze them, we consider a random reorientation of the leader’s heading by an
angle $\Delta\theta_{L}$, performed in the steady state corresponding to a
given value of the noise intensity $\eta$, and measure the subsequent
rearrangements that this perturbation induces in the heading of the rest of
fish, as given by the turning angle
$\varphi_{i}(t)=\theta_{i}(t+1)-\theta_{i}(t)$ projected on the interval
$[0,\pi]$.

Figure 4:  (a) Probability density of the turning angles $\varphi$ for the
Vicsek model with a leader, for different values of the noise intensity
$\eta$, in a system of size $L=220$. Dashed lines correspond to a non-turning
leader. Full lines represent a leader perturbed periodically every $250$ time
steps. Vertical lines, indicating the departure of the distributions for
perturbed and non-perturbed leaders, are estimated at a value
$\varphi_{c}(\eta)\simeq 2.4\pi\eta$ (b) Number of SPPs turning an angle
larger than $\varphi_{\mathrm{th}}(\eta)=2.5\pi\eta$ in a sequence of $10000$
simulation time steps in the Vicsek model with a periodically perturbed
leader, for different values of $\eta$. (c) Cumulated distribution of sizes
$P_{c}(s)$ of avalanches induced by a periodically perturbed leader in a
system of size $L=220$ with turning threshold $\varphi_{\mathrm{th}}(\eta)$,
for different values of the noise intensity. (d) Cumulated distribution of
durations $P_{c}(t)$ of avalanches induced by a periodically perturbed leader
for different values of the noise intensity.

In Fig. 4(a) we represent the probability density of SPPs turning angles
$P(\varphi)$ in the steady state, for different values of the noise intensity
$\eta$. In this plot we consider the model with a fixed, non-turning leader
(dashed lines), and the case of a periodically perturbed leader (full lines),
in which the leader experiences a random rotation $\Delta\theta_{L}$ of its
heading, uniformly distributed in the interval $[-\pi,\pi]$, every 250 time
steps, a time lapse larger than the maximum avalanche duration recorded in
simulations. As we can see, for fixed $\eta$, the two distributions are almost
identical for small $\varphi$, while they differ drastically regarding the
behavior of the tails beyond a given cut-off turning angle
$\varphi_{c}(\eta)$. A numerical analysis performed for different values of
$L$ allows to estimate this cut-off as $\varphi_{c}(\eta)\simeq 2.4\pi\eta$.
The presence of this turning angle cut-off, not available in empirical data,
permits to distinguish the changes of heading due to the effect of the leader
perturbations, and suggests that the proper definition of avalanches should
consider turning thresholds larger than the cut-off $\varphi_{c}(\eta)$. In
the following, we will fix the value of the threshold to
$\varphi_{\mathrm{th}}(\eta)=2.5\pi\eta$. We notice that for large $\eta\geq
0.5$, the angular distributions with and without perturbations are identical,
compatible with a large noise masking external perturbations and making
avalanches non discernible.

In Fig. 4(b) we plot a sample of the number of SPPs that turn an angle larger
than $\varphi_{\mathrm{th}}(\eta)$ as a function of time. This curve
emphasizes the heterogeneous character of the avalanche sizes in response to
the leader’s changes of direction, akin to what is observed in fish schools:
Sometimes a perturbation is followed by a small number of SPPs reorientations;
but other times, it triggers the reorientation of a large number of particles.
As expected, the strength of the effects of the leader perturbations decreases
with increasing noise, indicating that interesting avalanche behavior will
only occur for moderate levels of noise.

We compute the cumulative probability distributions $P_{c}(s)$ and $P_{c}(t)$,
Eq. (3), of observing an avalanche of size and duration larger than $s$ and
$t$, respectively, plotted in Fig. 4(c) and (d) for a turning threshold
$\varphi_{\mathrm{th}}(\eta)$ and different values of $\eta$. As we can see
from these plots, the values $\eta=0.2$ and $0.3$ lead to size an time
distributions analogous to that observed in rearrangement avalanches in real
fish schools, with a shape that can be approximated by a power-law form for
intermediate values, followed by a crossover to a sharp decrease for large $s$
and $t$ above a characteristic size or time. The behavior for $\eta=0.1$ is
more complex, probably due to the fact that for small noise one expects a
fairly homogeneous response with many SPPs following a leader perturbation. We
thus discard this value in the following analysis.

Figure 5:  (a) Cumulative probability distribution of size $P_{c}(s)$ of
avalanches induced by a perturbed leader in a system with $\eta=0.2$, turning
threshold $\varphi_{\mathrm{th}}(\eta)$ and different values of $L$. (b)
Cumulative probability distribution of durations $P_{c}(t)$ of avalanches
induced by a perturbed leader. (c) Check of the scaling of the cumulated size
distribution as given by Eq. (5). (d) Check of the scaling of the cumulated
time distribution as given by Eq. (5). Statistics are performed over at least
$10^{5}$ different avalanches.

The fact that we work now with a numerical model, allows us to explore the
behavior of the system for different systems sizes $L$ at a fixed turning
threshold, which was not possible in our fixed size empirical data. In Fig.
5(a) and (b) we plot the cumulative size and duration distributions in
avalanches in the Vicsek model with a turning leader for a turning threshold
$\varphi_{\mathrm{th}}(\eta)$, $\eta=0.2$ and different system sizes. As we
can observe, the behavior of the distributions is analogous to that observed
in real fish schools, compatible with a power-law decay but that are now
truncated by a size and time cut-offs that are functions of the system size
$L$. Inspired again by self-organized criticality Pruessner (2012), we can
assume now that the distributions obey a finite-size scaling form

$P(s)=s^{-\tau_{s}}G_{s}\left(\frac{s}{L^{D}}\right),\quad
P(t)=t^{-\tau_{t}}G_{t}\left(\frac{t}{L^{z}}\right).$ (7)

where $D$ and $z$ are new characteristic exponents that define the
characteristic size $s_{c}(L)\sim L^{D}$ and time $t_{c}(L)\sim L^{z}$ as a
function of the system size Pruessner (2012); Cardy (1988). The better
statistics in numerical simulations allow to estimate the characteristic
exponents applying the more precise moments analysis technique De Menech  _et
al._ (1998), see Appendix E. Application of this method leads to the
characteristics exponents reported in Table 1. We check the accuracy of these
values performing a data collapse analogous to that performed for the
avalanches in real fish, which, for the cumulated distributions, takes the
form,

$L^{D(\tau_{s}-1)}P_{c}(s)=F^{\prime}_{s}(sL^{D}),\quad
L^{z(\tau_{t}-1)}P_{c}(t)=F^{\prime}_{t}(tL^{z}),$ (8)

see Figs. 5(c) and (d) for $\eta=0.2$; the case $\eta=0.3$ is presented in
Fig. SF3. As we can see from these values, the exponents show a dependence on
the value of the noise $\eta$, although the size exponents appear to be
compatible within error bars. It is important to notice that the presence of a
rotating leader in necessary to obtain scaling avalanche distributions. Even
in the absence of a leader, the heading fluctuations due to noise and
interactions in the standard Vicsek model allow to define avalanches for a
given threshold. These avalanches, however, show a simple, short ranged
exponential distribution, as shown in Fig. SF4.

We have finally checked the effects of changing the turning threshold in the
scaling of the distributions as a function of the system size. In Fig. SF5 and
SF6 we show the results for a turning threshold
$\varphi_{\mathrm{th}}=2.8\pi\eta$, summarized in Table 1. As we can see, the
scaling exponent $\tau_{s}$ and $\tau_{t}$ in our model depend on the value of
the threshold. This fact is in contrast with the behavior of the fish school,
in which the characteristics exponents appear to be independent of the
threshold, and thus allowing for a scaling solution of the form given by Eq.
(4). Interestingly, the exponents $D$ and $z$ appear to be rather detail
independent, taking the approximate values $D\simeq 2$ and $z\simeq 1/2$ for
any value of $\eta$ and $\varphi_{\mathrm{th}}$, which would indicate that
avalanches in this model are compact Pruessner (2012).

##  IV Discussion

Behavioral cascades, taking the form of intermittent rearrangements
(avalanches) in the patterns of movement are an important, albeit sometimes
neglected, feature of collective motion in animals. Here we have shown that
behavioral cascades can be observed in the rearrangement dynamics of swimming
fish schools. Such avalanches, defined in terms of a turning threshold for the
heading of the fish, have distributions of sizes and times exhibiting a
scaling behavior compatible with a power-law tail truncated by a cut-off that
is an increasing function of the turning threshold. A data collapse analysis
allows to determine the exponents characterizing the scaling form. We
conjecture that such avalanche behavior can be due to the presence of
effective leadership in the schools. In order to support this conjecture, we
introduce a measure of leadership, based in the concept of avalanche
initiators, and observe that, indeed, some fish have consistently an unusually
large probability to initiate any avalanche in which they participate. These
predominant initiators can be interpreted as effective leaders, determining
the start of sudden rearrangements of the school headings. Leadership in the
context of avalanche initiation could account for individuals having sudden
behavioral changes or specific information about the environment, such as the
proximity of a wall.

To check whether the presence of leaders is enough to induce avalanche
behavior in collective motion, we have considered a very simple model,
consisting in a variation of the classical Vicsek model with the addition of a
global leader, which influences the movement of the particles, subject to
random changes in its heading. Interestingly, this simple model displays an
intermittent behavior qualitatively similar to that observed in real schools,
with avalanche size and duration distributions displaying a self-similar
scaling form.

Our results provide a new perspective on the avalanche behavior observed in
real collective motion situations Ginelli  _et al._ (2015), which can be
associated to a simple mechanism of leadership, observed in many natural
situations, indicating the possibility of a direct relation between these two
phenomena. Leadership in the present context of a moving school corresponds to
those individuals that first react to any external input, or that first
exhibit a random behavioral change, and preferentially start sudden
rearrangements of the trajectories of other fish in the school. Such
interpretation is validated by the numerical results from our model. It is
also worth emphasizing that, while it does not offer a perfect quantitative
prediction of the characteristic exponents, it nevertheless allows to
reproduce the scaling form of the avalanche distributions within a minimalist
modeling framework.

Different venues of future research stem from the results presented here. From
an empirical perspective, it would be interesting to further study the nature
of the avalanches observed in real schools, and to correlate them with other
physical properties measured in similar systems Herbert-Read  _et al._ (2011);
Katz  _et al._ (2011), as well as with other measures of leadership devised in
other contexts of collective motion Strandburg-Peshkin  _et al._ (2018). From
a numerical point of view, our results present new challenges in the
understanding of the properties of the proposed model. Indeed, a clearly open
question remains to ascertain the ultimate origin of the scaling behavior
observed in avalanches in a system in which no apparent critical transition
exists. Another interesting question regards the effects of leader switching
strategies. We expect the scale-free nature of the observed avalanches to be
preserved, provided that the influence of the leader, sensory wise, remains
rather long-ranged. In this sense, as we have numerically checked (data not
shown), a short-ranged leader, only with local influence over its nearest
neighbors, is not able to induce system-wide orientation rearrangements. On
the other hand, the value of the exponents associated to the size and duration
cutoffs are apparently independent of the noise intensity imposed on the
system. These observations hint towards a possible partial universality, which
is not shared, however, by the power-law decay exponent. Further work in this
direction is clearly needed in order to clarify these issues.

##  Appendix A A null model of fish avalanches

In the absence of any sort of dynamical correlations between the turning
angles of fish, the evolution of avalanches is purely determined by the
independent turning probability $P(\varphi)$ of each fish. As a null model of
avalanche behavior, we consider the case in which each fish independently
turns a angle $\varphi$ at each time step. Consider an avalanche of duration
$t$ and size $s$, starting at time $t^{\prime}=1$. If the avalanche lasts $t$
time steps, it means that at least one fish turned an angle larger than
$\varphi_{\mathrm{th}}$ every frame from $t^{\prime}=1$ to $t^{\prime}=t$, and
that no fish turned an angle larger than $\varphi_{\mathrm{th}}$ at frame
$t^{\prime}=t+1$. Under these conditions, the probability that a fish turns an
angle larger than $\varphi_{\mathrm{th}}$ in any frame is

$q=\int_{\varphi_{\mathrm{th}}}^{\pi}P(\varphi)\;d\varphi,$ (9)

and the probability that at least one fish turns an angle larger than
$\varphi_{\mathrm{th}}$ in a given frame is

$Q=1-(1-q)^{N},$ (10)

where $N$ is the number of fish. Therefore, the normalized probability that an
avalanche lasts for $t$ frames in this null model is

$P_{0}(t)=\frac{(1-Q)Q^{t}}{\sum_{t^{\prime}=1}^{\infty}(1-Q)Q^{t^{\prime}}}=(1-Q)Q^{t-1},\;t=1,2.\ldots,\infty,$
(11)

where we consider that avalanches have a minimal duration of one frame. That
is, in the uncorrelated null model, the avalanche duration distribution has an
exponential form, with average avalanche duration
$\langle{t}\rangle_{0}=\sum_{t=1}^{\infty}tP_{0}(t)=1/(1-Q)$

Consider now a frame in an avalanche of finite duration. In this frame, at
least one fish turned an angle larger than $\varphi_{\mathrm{th}}$, therefore
the the probability of observing the $s_{1}$ large turns in this frame is

$p_{1}(s_{1})=\frac{1}{Q}\binom{N}{s_{1}}q^{s_{1}}(1-q)^{N-s_{1}},\;s_{1}=1,2,\ldots,N,$
(12)

If the avalanche has duration $t$, at each frame a number $s_{1}$ of fish,
distributed with the probability Eq. (12), will turn a large angle. Therefore,
the distribution of sizes in avalanches of duration $t$, $P_{0}(s|t)$ will be
given the convolution of the probability Eq. (12) $t$ times with itself. The
form of this expression is hard to compute. However, we can approximate the
avalanche size distribution as follows: Since Eq. (12) is similar to a
binomial distribution, it is bell-shaped and centered at the average value

$\bar{s}_{1}=\sum_{s_{1}=1}^{N}s_{1}p_{1}(s_{1})=\frac{Nq}{Q}.$ (13)

Therefore, the average size of an avalanche of duration $t$ is

$\bar{s}_{t}=\frac{Nq}{Q}t,$ (14)

linear with $t$. Assuming that the relation between size $s$ and duration $t$
is tight, given the bounded distribution $p_{1}(s_{1})$, we can use relation
$s\simeq\frac{Nq}{Q}t$ and the distribution $P_{0}(t)$ from Eq. (11) to obtain
the probability transformation $P_{0}(t)dt=P_{0}(s)ds$, leading to

$P_{0}(s)\simeq\frac{1-Q}{Nq}e^{-sQ\ln(1/Q)/(Nq)},$ (15)

that is, an exponential decay with a characteristic size

$s_{c}=\frac{Nq}{Q\ln(1/Q)}.$ (16)

##  Appendix B Numerical data collapse analysis

We start from a set of avalanche size (or duration) distributions, that we
assume to fulfill the scaling relation, at the level of the cumulated
distributions,

$P_{c}(s)=s^{-\tau_{s}+1}F_{s}\left(s\varphi_{\mathrm{th}}^{\sigma_{s}}\right).$
(17)

In order to find the exponents $\tau_{s}$ and $\sigma_{s}$, we proceed as
follows: We consider general exponents $x_{s}$ and $y_{s}$, from which we can
write the new rescaled expressions

$\varphi_{\mathrm{th}}^{y_{s}(1-x_{s})}P_{c}(s)=F^{\prime}_{s}\left(s\varphi_{\mathrm{th}}^{y_{s}}\right),$
(18)

where $F^{\prime}_{s}(z)=z^{-\tau_{s}+1}F_{s}(z)$. Plotting
$\varphi_{\mathrm{th}}^{y_{s}(1-x_{s})}P_{c}(s)$ as a function of
$s\varphi_{\mathrm{th}}^{y_{s}}$, the curves for different values of
$\varphi_{\mathrm{th}}$ will collapse onto the universal function
$F^{\prime}_{s}(z)$ when $x_{s}=\tau_{s}$ and $y_{s}=\sigma_{s}$. We can
estimate this exponent by considering the difference of the curves for the
different values of $\varphi_{\mathrm{th}}$ and choosing the exponents
$\tau_{s}$ and $\sigma_{s}$ as the values of the exponents $x_{s}$ and $y_{s}$
that minimize this difference. To compute the difference, we locate the
interval of values of $s\varphi_{\mathrm{th}}^{y_{s}}$ common for all
$\varphi_{\mathrm{th}}$. In this interval, we compute a spline of order $k$
for each quantity $\varphi_{\mathrm{th}}^{y_{s}(1-x_{s})}P_{c}(s)$ and
interpolate a fixed number $n$ of equispaced points. The difference is defined
as the sum of the variances of the values of
$\varphi_{\mathrm{th}}^{y_{s}(1-x_{s})}P_{c}(s)$ in each point of the
interpolation, for the different values of $\varphi_{\mathrm{th}}$. In the
results presented here, we consider splines of order $k=2$ and interpolate
$n=10$ points for each $P_{c}(s)$ curve.

##  Appendix C Leadership probability in the null model of avalanche behavior

In the avalanche null model defined above, consider a fish that participates
in a given avalanche. To estimate its leadership probability we have to
compute the probability that it leads the avalanche (i.e. it is active in its
first time step), provided that it participates in it. To compute it, we use
Bayes’ theorem to write

$P(p)P(l|p)=P(l)P(p|l),$ (19)

where $P(p)$ is the probability that the fish participates in a given
avalanche, $P(l|p)$ the probability of leading an avalanche in which it
participates (the probability we are seeking), $P(l)$ the probability of
leading an avalanche, and $P(p|l)$ the probability that a fish participates in
an avalanche provided that it leads it. Obviously, $P(p|l)=1$. To estimate the
rest of probabilities, we need information about the duration $t$ of the
avalanche. Thus, we have $P(p)=1-(1-q)^{t}$, the probability that the fish
turns at least once in the development of the avalanche, and $P(l)=q$, the
probability that the fish is active (performs a large turn) in the first time
step of the avalanche. Therefore, from Eq. (19) we obtain

$P(l|p)=\frac{P(l)P(p|l)}{P(p)}=\frac{q}{1-(1-q)^{t}}.$ (20)

Within this null model, consider a fish that participates in $N_{a}$
avalanches, each of duration $t_{\alpha}$, $\alpha=1,\ldots,N_{a}$. The
probability of leading any of these avalanches is
$p_{\alpha}=q/[1-(1-q)^{t_{\alpha}}]$. Therefore, the probability $P(\ell)$ of
leading $\ell$ of the $N_{a}$ avalanches is given by a Poisson binomial
distribution, representing the probability distribution of a sum of
independent Bernoulli trials that have different success probabilities
$p_{\alpha}$ Hoeffding (1956). The Poisson binomial distribution has a rather
convoluted form, but its mean and variance can be easily expressed as

$\displaystyle\mu$ $\displaystyle=$ $\displaystyle\sum_{\ell}\ell
P(\ell)=\sum_{\alpha}p_{\alpha},$ (21) $\displaystyle\sigma^{2}$
$\displaystyle=$
$\displaystyle\sum_{\ell}\ell^{2}P(\ell)-\left[\sum_{\ell}\ell
P(\ell)\right]^{2}=\sum_{\alpha}p_{\alpha}(1-p_{\alpha}).$ (22)

The average leadership probability of a fish in this null model is thus given
by $\chi_{0}=\langle{\ell}\rangle/N_{a}$, where
$\langle{\ell}\rangle=\sum_{\alpha}p_{\alpha}$ is the average number of
avalanches led by the fish. In Fig. 3 we show the actual values of $\chi_{i}$
computed for each fish. The full line and shaded region represents the null-
model average value $\chi_{i}^{0}$ computed for each fish, taking into account
the number of avalanches in which it participates, and its $99\%$ confidence
interval, respectively.

##  Appendix D Vicsek model with leadership

In the classic Vicsek model Vicsek  _et al._ (1995); Ginelli (2016), $N$ self-
propelled particles (SPPs) move in a two dimensional space. The dynamics is
overdamped and defined in discrete time, with the instantaneous position
$\vec{r}_{i}(t)$, $i=1,\ldots,N$, of each particle being related with its
velocity $\vec{v}_{i}(t)$ by

$\vec{r}_{i}(t+\Delta t)=\vec{r}_{i}(t)+\vec{v}_{i}(t)\Delta t,$ (23)

where $\Delta t$ is an integration time step, arbitrarily fixed to $\Delta
t=1$. Velocities have a constant modulus, $|\vec{v}_{i}(t)|=v_{0}$, and thus
are fully determined by their direction, given by the heading angle
$\theta_{i}(t)$ that the velocity forms with, say the $x$ axis, such that
$\vec{v}_{i}(t)=(v_{0}\cos\theta_{i}(t),v_{0}\sin\theta_{i}(t))$. Heading is
assumed to belong to the interval $[-\pi,\pi]$. In this model each particle
$i$ tends to orient its direction of motion along the average direction
$\vec{V}_{i}(t)$ of the particles located inside a circular area
$\mathcal{V}_{i}$ of radius $R$ centered at its own position and including
itself, i.e.

$\vec{V}_{i}(t)=\frac{1}{n_{i}(t)}\sum_{j\in\mathcal{V}_{i}}\vec{v}_{j}(t),$
(24)

where $n_{i}(t)$ is the number of particles in the neighborhood
$\mathcal{V}_{i}$ at time $t$. This dynamics is implemented in the update rule
for the heading angle

$\theta_{i}(t+\Delta
t)=\Theta\left[\mathbf{V}_{i}(t)\right]+\eta\;\xi_{i}(t),$ (25)

where the function $\Theta[\vec{V}]$ represents the angle of vector $\vec{V}$
and $\xi_{i}(t)$ a random noise, uniformly distributed in the interval
$[-\pi,\pi]$, and $\eta\in[0,1]$ a parameter gauging the strength of the noise
term.

This model exhibits an order-disorder phase transition defined in terms of an
order parameter (polarization) given by

$\phi(\eta)=\frac{1}{v_{0}N}\left\langle\left|\sum_{i=1}^{N}\mathbf{v}_{i}(t)\right|\right\rangle_{t},$
(26)

the brackets representing a temporal average. The transition separates an
ordered phase for noise strength smaller than a critical value $\eta_{c}$,
corresponding to a flocking (schooling) phase, from a disordered phase for
$\eta>\eta_{c}$, corresponding to a swarming, disordered phase.

In the variation of the Vicsek model we consider, one of the SPPs, say
particle $1$ plays the role of a leader which influences the orientation of
the rest of SPPs in the system, independently of their relative distance.
Therefore, in the heading update rule Eq. (25), the average velocity of the
neighbors is replaced by the average $\vec{V}_{i}^{L}(t)$ computed in the set
$\mathcal{V}_{i}^{L}=\mathcal{V}_{i}\;\cup\;\\{1\\}$, including the global
leader and all the particles in the local neighborhood of $i$. The heading and
velocity of this leader is constant in time, $\theta_{1}(t)=\theta_{L}$, and
it represents a privileged direction it wants to follow.

Simulations of the model are performed in square boxes of different size $L$
with periodic boundary conditions. The we fix the density of particles
$\rho=N/L^{2}=1$, the radius of interaction $R=1$, and the constant speed of
the SPPs $v_{0}=0.03$.

##  Appendix E Moments analysis technique

The finite-size scaling (FSS) method Cardy (1988) assumes that the dependence
on system size $L$ of the avalanche size and time distributions is of the form

$\displaystyle P(s,L)$ $\displaystyle=$ $\displaystyle
s^{-\tau_{s}}\mathcal{F}_{s}\left(\frac{s}{s_{c}(L)}\right),$ (27)
$\displaystyle P(t,L)$ $\displaystyle=$ $\displaystyle
t^{-\tau_{t}}\mathcal{F}_{t}\left(\frac{t}{t_{c}(L)}\right),$ (28)

where $\mathcal{F}_{x}(z)$ are scaling functions that are approximately
constant for $z<1$, and decay very fast to zero for $z>1$. The quantities
$s_{c}(L)$ and $t_{c}(L)$ are the cut-off characteristics size and time, which
are assume to depend on system size as $s_{c}(L)\sim L^{D}$ and $t_{c}(L)\sim
L^{z}$, thus defining the standard critical exponents $\tau_{s}$, $\tau_{t}$,
$D$ (the fractal dimension) and $z$ (the dynamic critical exponent) Pruessner
(2012).

Assuming the scaling form given by Eqs. (27) and (28), we can compute
numerically the associated critical exponents applying the moment analysis
technique De Menech  _et al._ (1998). One starts by defining the $q$-th moment
of the avalanche size distribution on a box of size $L$ as

$\displaystyle\langle{s^{q}}\rangle_{L}$ $\displaystyle=$
$\displaystyle\sum_{s}s^{q}\;P(s,L)\simeq\int
ds\;s^{-\tau_{s}+q}\mathcal{F}_{s}\left(\frac{s}{L^{D}}\right)$ (29)
$\displaystyle=$ $\displaystyle L^{D(q+1-\tau_{s})}\int
dx\;y^{-\tau_{s}+q}\mathcal{F}_{s}(x)\sim L^{\sigma_{s}(q)},$

where we have introduce the FSS form in Eq. (27), and taken the continuous
approximation for the $s$ and $t$ variables. The exponents
$\sigma_{s}(q)\equiv D(q+1-\tau_{s})$ can be estimated as the slope of the
numerical evaluation of $\langle{s^{q}}\rangle_{L}$ as a function of $L$ in a
double logarithmic plot. Then, for sufficiently large values of $q$, we can
perform a linear fit of the exponent $\sigma_{s}(q)$ to the form

$\sigma_{s}(q)=Aq+B,$ (30)

with $A=D$ and $B=D(1-\tau_{s})$, from where $D$ and $\tau_{s}$ can be
directly estimated. Along the same lines, the exponents associated to the
avalanche time distribution can be evaluated considering the $q$-th moment of
the time distribution, $\langle{t^{q}}\rangle_{L}\sim L^{\sigma_{t}(q)}$, with
$\sigma_{t}(q)\equiv z(q+1-\tau_{t})$.

##  Appendix F Supplementary Figures

Supplementary Figure SF1: Plot of the leadership probability $\chi_{i}$ as a
function of the turning threshold $\varphi_{\mathrm{th}}$, for the different
fish in each series A, B, C (top to bottom). Full lines represent the average
leadership probability in a null model of uncorrelated avalanches. The
shadowed regions represents the $99\%$ confidence interval of this value.
Notice that in some plots certain values of $\varphi_{\mathrm{th}}$ are
missing. This is due to the fact that the corresponding fish do not
participate in any avalanche, and therefore its leadership probability is not
defined. Supplementary Figure SF2: Average order parameter
$\langle{\phi(\eta)}\rangle$ as a function of noise intensity in the classic
standard Vicsek model (SVM) and the Vicsek model with a non-rotating global
leader for different system sizes. Supplementary Figure SF3:  (a) Cumulative
probability distribution of size $P_{c}(s)$ of avalanches induced by a
rotating leader in a system with $\eta=0.3$, turning threshold
$\varphi_{\mathrm{th}}(\eta)$ and different values of $L$. (b) Cumulative
probability distribution of durations $P_{c}(t)$ of avalanches induced by a
rotating leader in a system with $\eta=0.2$, turning threshold
$\varphi_{\mathrm{th}}(\eta)$ and different values of $L$. In both cases,
statistics is performed over at least $10^{5}$ different avalanches. (c) Check
of the scaling of the cumulated size distribution with turning threshold
$\varphi_{\mathrm{th}}(\eta)$, as given by Eq. (8) in the main paper. (d)
Check of the scaling of the cumulated time distribution with turning threshold
$\varphi_{\mathrm{th}}(\eta)$, as given by Eq. (8) in the main paper.
Supplementary Figure SF4:  (a) Cumulative probability distribution of size
$P_{c}(s)$ of avalanches in the standard Vicsek model with $\eta=0.2$, turning
threshold $\varphi_{\mathrm{th}}(\eta)$ and different values of $L$. (b)
Cumulative probability distribution of duration $P_{c}(t)$ in the standard
Vicsek model with $\eta=0.2$, turning threshold $\varphi_{\mathrm{th}}(\eta)$
and different values of $L$. (c) Cumulative probability distribution of size
$P_{c}(s)$ of avalanches in the standard Vicsek model with $\eta=0.3$, turning
threshold $\varphi_{\mathrm{th}}(\eta)$ and different values of $L$. (d)
Cumulative probability distribution of duration $P_{c}(t)$ in the standard
Vicsek model with $\eta=0.3$, turning threshold $\varphi_{\mathrm{th}}(\eta)$
and different values of $L$. Statistics are performed over at least $10^{5}$
different avalanches. Supplementary Figure SF5:  (a) Cumulative probability
distribution of size $P_{c}(s)$ of avalanches induced by a rotating leader in
a system with $\eta=0.2$, turning threshold $\varphi_{\mathrm{th}}=2.8\pi\eta$
and different values of $L$. (b) Cumulative probability distribution of
durations $P_{c}(t)$ of avalanches induced by a rotating leader in a system
with $\eta=0.2$, turning threshold $\varphi_{\mathrm{th}}(\eta)$ and different
values of $L$. In both cases, statistics is performed over at least $10^{5}$
different avalanches. (c) Check of the scaling of the cumulated size
distribution with turning threshold $\varphi_{\mathrm{th}}(\eta)$, as given by
Eq. (8) in the main paper. (d) Check of the scaling of the cumulated time
distribution with turning threshold $\varphi_{\mathrm{th}}(\eta)$, as given by
Eq. (8) in the main paper. Supplementary Figure SF6:  (a) Cumulative
probability distribution of size $P_{c}(s)$ of avalanches induced by a
rotating leader in a system with $\eta=0.3$, turning threshold
$\varphi_{\mathrm{th}}=2.8\pi\eta$ and different values of $L$. (b) Cumulative
probability distribution of durations $P_{c}(t)$ of avalanches induced by a
rotating leader in a system with $\eta=0.3$, turning threshold
$\varphi_{\mathrm{th}}(\eta)$ and different values of $L$. In both cases,
statistics is performed over at least $10^{5}$ different avalanches. (c) Check
of the scaling of the cumulated size distribution with turning threshold
$\varphi_{\mathrm{th}}(\eta)$, as given by Eq. (8) in the main paper. (d)
Check of the scaling of the cumulated time distribution with turning threshold
$\varphi_{\mathrm{th}}(\eta)$, as given by Eq. (8) in the main paper.

###### Acknowledgements.

We acknowledge financial support from the Spanish
MCIN/AEI/10.13039/501100011033, under Projects No. PID2019-106290GB-C21 and
No. PID2019-106290GB-C22. We thank F. S. Beltrán and V. Quera, from the
Institute of Neurosciences, University of Barcelona (Spain), for providing the
raw data of fish recordings.

## References

  * Vicsek and Zafeiris (2012) T. Vicsek and A. Zafeiris, Phys. Rep. 517, 71 (2012).
  * Sumpter (2010) D. J. Sumpter,  _Collective Animal Behavior_ (Princeton University Press, New Jersey, 2010).
  * Ramaswamy (2010) S. Ramaswamy, Annual Review of Condensed Matter Physics 1, 323 (2010).
  * Rosenthal  _et al._ (2015) S. B. Rosenthal, C. R. Twomey, A. T. Hartnett, H. S. Wu, and I. D. Couzin, Proc. Natl. Acad. Sci. 112, 4690 (2015).
  * Chen  _et al._ (2016) D. Chen, T. Vicsek, X. Liu, T. Zhou, and H. Zhang, Europhys. Lett. 114, 60008 (2016).
  * Calovi  _et al._ (2018) D. S. Calovi, A. Litchinko, V. Lecheval, U. Lopez, A. Pérez Escudero, H. Chaté, C. Sire, and G. Theraulaz, PLOS Computational Biology 14, 1 (2018).
  * Brambilla  _et al._ (2013) M. Brambilla, E. Ferrante, M. Birattari, and M. Dorigo, Swarm Intelligence 7, 1 (2013).
  * Krause and Ruxton (2002) J. Krause and G. Ruxton,  _Living in Groups_ (Oxford University Press, Oxford, 2002).
  * Procaccini  _et al._ (2011) A. Procaccini, A. Orlandi, A. Cavagna, I. Giardina, F. Zoratto, D. Santucci, F. Chiarotti, C. K. Hemelrijk, E. Alleva, G. Parisi, and C. Carere, Animal Behaviour 82, 759 (2011).
  * Ginelli  _et al._ (2015) F. Ginelli, F. Peruani, M.-H. Pillot, H. Chaté, G. Theraulaz, and R. Bon, Proceedings of the National Academy of Sciences 112, 12729 (2015).
  * Couzin  _et al._ (2005) I. D. Couzin, J. Krause, N. R. Franks, and S. A. Levin, Nature 433, 513 (2005).
  * Aubé and Shield (2004) F. Aubé and R. Shield, in  _Lect. Notes Comput. Sci._, Vol. 3305 (Springer Berlin / Heidelberg, 2004) pp. 601–611.
  * Shen (2008) J. J. Shen, SIAM Journal on Applied Mathematics 68, 694 (2008), https://doi.org/10.1137/060673254 .
  * Pearce and Giomi (2016) D. J. G. Pearce and L. Giomi, Phys. Rev. E 94, 022612 (2016).
  * Kyriakopoulos  _et al._ (2016) N. Kyriakopoulos, F. Ginelli, and J. Toner, New Journal of Physics 18, 73039 (2016).
  * Smith  _et al._ (2016) J. E. Smith, S. Gavrilets, M. B. Mulder, P. L. Hooper, C. E. Mouden, D. Nettle, C. Hauert, K. Hill, S. Perry, A. E. Pusey, M. van Vugt, and E. A. Smith, Trends in Ecology & Evolution 31, 54 (2016).
  * Nagy  _et al._ (2010) M. Nagy, Z. Akos, D. Biro, T. Vicsek, Z. Ákos, D. Biro, and T. Vicsek, Nature 464, 890 (2010), arXiv:1010.5394 .
  * Flack  _et al._ (2012) A. Flack, B. Pettit, R. Freeman, T. Guilford, and D. Biro, Animal Behaviour 83, 703 (2012).
  * Nagy  _et al._ (2013) M. Nagy, G. Vásárhelyi, B. Pettit, I. Roberts-Mariani, T. Vicsek, and D. Biro, Proceedings of the National Academy of Sciences of the United States of America 110, 13049 (2013).
  * Couzin  _et al._ (2011) I. D. Couzin, C. C. Ioannou, G. Demirel, T. Gross, C. J. Torney, a. Hartnett, L. Conradt, S. a. Levin, and N. E. Leonard, Science 334, 1578 (2011).
  * Ward  _et al._ (2011) A. J. W. Ward, J. E. Herbert-Read, D. J. T. Sumpter, and J. Krause, Proceedings of the National Academy of Sciences 108, 2312 (2011).
  * Fisher (1998) D. S. Fisher, Physics Reports 301, 113 (1998), arXiv:9711179 [cond-mat] .
  * Zapperi  _et al._ (1998) S. Zapperi, P. Cizeau, G. Durin, and H. E. Stanley, Phys. Rev. B 58, 6353 (1998).
  * Altshuler and Johansen (2004) E. Altshuler and T. H. Johansen, Rev. Mod. Phys. 76, 471 (2004).
  * Miguel  _et al._ (2001) M. C. Miguel, A. Vespignani, S. Zapperi, J. Weiss, and J.-R. Grasso, Nature 410 (2001), 10.1038/35070524.
  * Zapperi  _et al._ (1999) S. Zapperi, P. Ray, H. E. Stanley, and A. Vespignani, Phys. Rev. E 59, 5049 (1999).
  * Kawamura  _et al._ (2012) H. Kawamura, T. Hatano, N. Kato, S. Biswas, and B. K. Chakrabarti, Rev. Mod. Phys. 84, 839 (2012).
  * Vicsek  _et al._ (1995) T. Vicsek, A. Czirok, E. Ben-Jacob, I. Cohen, and O. Shochet, Phys. Rev. Lett. 75, 1226 (1995).
  * Pruessner (2012) G. Pruessner,  _Self-Organised Criticality: Theory, Models and Characterisation_, Self-organised Criticality: Theory, Models, and Characterisation (Cambridge University Press, Cambridge, UK, 2012).
  * Méndez  _et al._ (2014) V. Méndez, D. Campos, and F. Bartumeus,  _Stochastic Foundations in Movement Ecology_ (Springer Verlag, Berlin, Heidelberg, 2014).
  * Gimeno  _et al._ (2016) E. Gimeno, V. Quera, F. S. Beltran, and R. Dolado, Journal of Comparative Psychology 130, 358 (2016).
  * Press and Teukolsky (1990) W. H. Press and S. A. Teukolsky, Computers in Physics 4, 669 (1990).
  * Fornberg (1988) B. Fornberg, Mathematics of Computation 51, 699 (1988).
  * Laurson and Alava (2006) L. Laurson and M. J. Alava, Phys. Rev. E 74 (2006).
  * Laurson  _et al._ (2009) L. Laurson, X. Illa, and M. J. Alava, J. Stat. Mech.: Theo. Exp. P01019 (2009).
  * Yeomans (1992) J. M. Yeomans,  _Statistical mechanics of phase transitions_ (Oxford University Press, Oxford, 1992).
  * Cardy (1996) J. Cardy,  _Scaling and Renormalization in Statistical Physics_, Cambridge Lecture Notes in Physics (Cambridge University Press, 1996).
  * Strandburg-Peshkin  _et al._ (2018) A. Strandburg-Peshkin, D. Papageorgiou, M. C. Crofoot, and D. R. Farine, Philosophical Transactions of the Royal Society B: Biological Sciences 373, 20170006 (2018).
  * Ginelli (2016) F. Ginelli, Eur. Phys. J. Spec. Top. 225, 2099 (2016).
  * Cardy (1988) J. L. Cardy, ed.,  _Finite Size Scaling_ , Current Physics-Sources and Comments, Vol. 2 (North Holland, Amsterdam, 1988).
  * De Menech  _et al._ (1998) M. De Menech, A. L. Stella, and C. Tebaldi, Phys. Rev. E 58, R2677 (1998).
  * Herbert-Read  _et al._ (2011) J. E. Herbert-Read, A. Perna, R. P. Mann, T. M. Schaerf, D. J. T. Sumpter, and A. J. W. Ward, Proceedings of the National Academy of Sciences 108, 18726 (2011).
  * Katz  _et al._ (2011) Y. Katz, K. Tunstrøm, C. C. Ioannou, C. Huepe, and I. D. Couzin, Proceedings of the National Academy of Sciences 108, 18720 (2011), https://www.pnas.org/content/108/46/18720.full.pdf .
  * Hoeffding (1956) W. Hoeffding, The Annals of Mathematical Statistics 27, 713 (1956).