SciPost Physics

Submission

Inferring leader-follower behavior from presence data in themarine environment: a case study on reef manta rays
Juan Fernández-Gracia1? , Jorge P. Rodríguez2 , Lauren R. Peel3 , Konstantin Klemm1 , Mark G.
Meekan3 ,and Víctor M. Eguíluz1
1 Instituto de Física Interdisciplinar y Sistemas Complejos IFISC (CSIC-UIB), E07122 Palmade Mallorca, Spain
2 IMEDEA
3 Australian Institute of Marine Science, Indian Ocean Marine Research Centre (IOMRC),
University of Western Australia (M470), 35 Stirling Highway, Crawley, WA 6009, Australia?

juanf@isc.uib-csic.es
April 12, 2022

1

Abstract

224

Social interactions are ubiquitous in groups of animals, including humans. These interactions might be of different nature, for example, competitive, mutualistic, or kinship;
and their global structure is naturally studied with the tools of complex network theory.
Traditionally, it has been challenging to examine social interactions in the marine environment given the difficulties associated with data collection, however, developments inacoustic telemetry technologies now present a novel way to remotely examine such behaviour. Here, we propose a method to extract leader-follower networks from presence
data collected by a single acoustic receiver at a specific location. The method is basedon the Kolmogorov-Smirnov distance between the distribution of lag times between theconsecutive presence of an individual followed by the presence of another and its conjugate distribution. After characterising the method through controlled generated data,
we apply it to detection data collected for acoustically tagged reef manta rays (Mobula
alfredi) by a single acoustic receiver positioned at a frequently visited site. First, weshow that the presence of reef manta rays in the vicinity of the cleaning station displays
several temporal heterogeneities, such as a circadian rhythm, as well as burst-like behavior, where the time between consecutive presences follows a power-law distribution.
Second, we infer the leader-follower network of manta rays and characterize individuals in terms of their position on this network relative to their sex and size. We find,in agreement with biological and ecological insights, that (i) female reef mantas follow
more males than expected; (ii) male reef manta rays follow fewer females than expected,
but with a stronger association to certain individuals; (iii) reef manta rays follow eachother with a weaker association than larger individuals do, while the rest of the reported
interactions between individuals appear to be random.

25

1

26

Social interactions among animals mediate many processes, such as the transmission of information [1, 2] or diseases [3]; collective behavior [4–6] (flocking [7], social learning [8],predator avoidance [9], cooperative foraging [10]); selection of phenotypes [11]; mating [10];or the emergence and maintenance of cooperation [10,12]. The structure of these social interactions as a whole is typically studied with the methods from network theory [13], which have

3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

27
28
29
30

Introduction

1

SciPost Physics

31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

Submission

recently and increasingly been applied to animal groups [1–10, 14] and represent a promising
tool in the field of movement ecology [15].
The first step of studying these interactions involves finding or collecting the data that describe them, before the underlying networks can be extracted in a meaningful and statisticallysignificant way. For humans, the availability of large amount of digital traces (notably, call detail records and social media accounts [16]) and the ease at which they can be accessed, has
facilitated the detailed statistical description of human interactions. Tracking animal social
systems [17] in a similar manner, however, has remained a challenge, given the difficulties
associated with collecting a sufficient amount of data from the network of interest. This isparticularly true in the marine environment, where the large and three-dimensional nature ofanimal movements impose logistical, technological and financial constraints on data collection
capacities.
Traditionally, the collection of presence data in the marine environment relied on techniques that were applied within discrete sampling periods (e.g. photo-identification [18] or
capture-mark-recapture [19, 20]). This temporal limitation restricted the extent to which social interactions between sampled individuals could be considered. Recent advances in thefield of acoustic telemetry, however, have since provided a means to overcome these issues,
allowing the movement and residency patterns of marine fauna to be monitored continuously,
and over long periods of time (up to ten years).
In studies using passive acoustic telemetry, acoustic tags that emit a unique sound signal
are externally attached to [21], or surgically implanted into [22], study animals. These tags
are subsequently detected by acoustic receivers placed at specific locations within the study
area (collectively, an acoustic array), and the timestamp of each detection is logged by thedevices from which data can be downloaded at a later date. Acoustic arrays and tagging programs have been established worldwide [23–25], and the frequency of their use is increasing
as the affordability of the required equipment improves [23, 26]. While the primary motivation of such studies is often to obtain spatial data to inform the development of conservation
measures [21, 23, 27], the potential to use these presence data to examine leader-follower
behavior in marine species is yet to be extensively explored.
Leadership behavior has been reported in many animal groups, including insects (e.g. ants)
and birds (e.g. migrating storks), and has been found to be of paramount importance inachieving coordination among individuals [28–30]. Should a population lack clear leaders ora hierarchical structure [31], individuals may have preferences on whom they should interact
with, subsequently modifying the strength of the social interaction. Additionally, interactions
within a social network may not be reciprocal, generating directionality in the flow of information. This allows for the introduction of focal individuals and their leading and following
connections, such that the focal individuals’ dynamics are coupled with that of its leaders,
while that of its followers are coupled with the focal individual’s.
Within a typical social network, individuals are represented by nodes, and the interactions
recorded between them are represented by connections (or edges). The power of networks torepresent social interactions lies within the ability to characterize edges relative to the type ofinteraction that has occurred. Interactions can be symmetrical, representing exchanges fromone peer to the other and vice versa, or directed, originating from a source towards a target.
In this sense, leadership behavior can be examined with network tools, too, assuming that thesources are leaders and the targets are followers. However, appropriate methods for detecting
these asymmetric relations are currently lacking.
One can find several social network inference methods in the literature that are well suited
for acoustic telemetry data [14, 32–38]. The majority of these rely on the co-occurrences ofindividuals, and thus have the problem of the Gambit of the group (GoG [14]). That is, every
pair of animals observed in the same group is treated as having interacted with one another,

2

SciPost Physics

Submission

102

although this might appear due to coincidence in visitation patterns and not social affiliation.
Early analytical methods also required an appropriate temporal window to be defined, within
which co-occurrences were considered [36–39]. This restriction, and that of the GoG, was resolved by the introduction of the ‘GMMevents’ method by Psorakis and collaborators [32, 33].
These methods allowed researchers to investigate the complex social structure of different animal groups, uncovering their relation to individual characteristics, such as sex, size, personality
or genetic traits [11, 34, 40, 41].
Here, we aim to develop an analytical method to provide evidence of leader-follower interactions from acoustic telemetry data, alongside estimates for the statistical significance of theobserved interactions. We first describe a novel method for extracting leader-follower interactions from established social networks. We then apply this method to real-world, presence
data collected for reef manta rays (Mobula alfredi) at a single location using passive acoustic
telemetry. Reef manta rays present an ideal study species for work of this nature for a numberof reasons: (1) they show repeated and prolonged visitation to key sites on coral reefs known
as cleaning stations, where individual mantas socialize, and cleaner fishes remove parasites
from their bodies [42]; (2) individuals are long-lived and can be accurately identified fromunique pigmentation patterns displayed on their bodies throughout their lifetime [18,43]; (3)reef manta rays exhibit two distinct sexes and the size of an individual can be used as a proxyfor its maturity status, allowing leader-follower behaviors to be considered relative to thesebiological traits. Finally, we discuss the applicability of this novel analysis procedure to existing acoustic telemetry datasets, and the potential to expand these works to investigate theleader-follower behavior of populations across entire acoustic arrays.

103

2

104
105

All the methods described here have been implemented computationally using Python 3 andare available through the page in Ref. [44].

106

2.1

107

We propose a measure for quantifying leadership in multiple individual time series. In particular, our method is relevant for examining the collective movement patterns of entities,
be it active matter, robots, humans or animals, where presence data of individuals has beenrecorded at a specific location.
The data for each individual consists of an ordered set of timestamps at which the individual is located in the vicinity of that particular location. In order to infer if individual i typicallyfollows individual j, we hypothesize that the time delay, or lag time, between the consecutive
detection of first i and then j will be longer than the reverse (i.e. the consecutive times ofdetecting j followed by i). We therefore extract the lag times for j followed by i, {t i j }, andi followed by j, {t ji }, and their corresponding distributions pi j (t) (p ji (t)), such that pi j (t)d t(p ji (t)d t) measures the probability that a randomly selected lag time of i after j ( j after i)lays in the interval (t, t + d t). We then compute the distance between the distributions of lagtimes as the Kolmogorov-Smirnof distance [45,46] (DKS ) between them, where DKS is definedas the maximum distance between the two cumulative distribution functions P1 and P2 as

81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101

108
109
110
111
112
113
114
115
116
117
118
119
120

Materials and methods

Leader-follower interactions from presence data

DKS (P1 , P2 ) = max(|P1 (t) − P2 (t)|).
t

121
122
123

(1)

We use a signed version of DKS to distinguish which of the two cumulative distributions isgreater than the other. Since this quantity is no longer a metric, we now call it the KolmogorovSmirnof arrow AKS .
AKS (P1 , P2 ) = P1 (τ) − P2 (τ),(2)
3

SciPost Physics

Submission

133

where τ = arg max t (|P1 (t) − P2 (t)|). For measuring AKS , we find the time, τ, that maximizes
the distance between the two compared cumulative distributions, and then draw an arrowfrom the second to the first distribution (first and second refers to their order in the argumentsof AKS ). If the arrow goes upwards, AKS = DKS and otherwise AKS = −DKS . For two interevent
sequences 1 and 2, AKS (P1 , P2 ) > 0 implies that times associated with process 1 are shorter
than those associated with process 2. Thus, if individual j is following individual i, we will
expect AKS (Pi j , P ji ) > 0, while AKS (Pi j , P ji ) < 0 reflects the opposite (i following j). The
strength of the interaction will be given by the absolute value of the arrow. The value of thearrow is related to the excess of probability to have smaller values than τ of the cumulative
distribution that lies above.

134

2.2

135

We create synthetic presence data for a pair of individuals (hereafter, Individual A and Individual B) proposing a model of temporal sequences. The first set of times, corresponding
to Individual A, is given by a homogeneous Poisson process of rate λ1 , i.e., Individual A appears randomly at a constant rate. The temporal sequence for Individual B follows from a
non-homogeneous Poisson process with rate λ2 (t), which depends on time in such a way thatit has a constant rate λ∗ plus an excess rate δ in the interval (−λ∗ , ∞) during an amount oftime ∆t immediately after the presence of Individual A is recorded. This is similar to a Hawkes
process [47] and to the model for correspondence by Malmgren et al. [48], but in this casethe excitations are given by an independent source and are not self-excitations. If ∆t = 0,
we just have two independent homogeneous Poisson processes of rates λ1 and λ2 (t) = λ∗ ,wheras for ∆t → ∞, we have two independent homogeneous Poisson processes of rates λ1and λ2 = λ∗ +δ. For an intermediate ∆t, Individual B may attempt to actively be present at, or
avoid, the location for a period of time ∆t, shortly after Individual A was there for a positive,
or negative, value of δ, respectively. If δ = 0, we again have two independent Poisson processes of rates λ1 and λ2 = λ∗ . In the following, we choose λ1 = 1 without loss of generality.
From the time series for Individuals A and B we extract the sets of lag times {t 12 } and {t 21 }and compute the Kolmogorov-Smirnof arrow AKS ({t 12 }, {t 21 }) so that if AKS > 0, Individual
A is following Individual B, and if AKS < 0 Individual B is following Individual A. Several realizations of the process were completed, and the distribution of arrows compiled to examine
whether the KS arrow is able to capture the aforementioned leader-follower behaviours (seeFigures 1 and 2).
For two independent Poisson processes of different rates (δ = 0, λ1 6= λ2 ), the distributionof KS arrows are always symmetric, with low values centered around, but not equal to, zero.
This reflects the fact that the sequences are uncorrelated, and that it is not possible to discern
if one individual is following the other. The absence of values around zero is related to thefact that the KS arrow will only be zero if the compared distributions are exactly equal. Thisequality will be only occur when the sampling is infinite, which was confirmed by drawing
longer and longer time series sequences and observing that the distribution collapsed towards
a single delta function at zero (not shown here). Regarding the correlated case, we checked
how the distribution of arrows changes as we change the parameters δ and ∆t, setting λ∗ = 1(not shown here). If δ = 0, ∆t = 0 or ∆t much bigger than 1/λ1 the distribution is symmetricaround 0 with values close to it, reflecting again that there is no evidence of one individual
following the other, as expected for independent presence sequences. For δ < 0, the distribution is non-zero only for positive values of the arrow, indicating that the first individual follows
the second (although the model is just assuming that the second one is avoiding the first one),while for δ > 0 the distribution is shifted to negative values of the arrow, showing that thesecond individual is following the first.
Lastly, we examined the behavior of the distribution of KS arrows when the two time series

124
125
126
127
128
129
130
131
132

136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172

Application to a presence toy model

4

SciPost Physics

Submission

185

sequences were reshuffled in such a way that both individuals retained the same number of
events, but without correlations, if there were any. For both the correlated and uncorrelatedcases, the result was the same: the distribution of KS arrows for the reshuffled sequences istypical of a pair with no leader-follower interactions, i.e., with low arrow values, centeredaround but not equal to zero (see Fig. 3). Actually if the pair of sequences was coming from
two independent random sequences, the distribution of arrows for the reshufflings is identicalto the distribution of KS arrows for an ensemble of independent pairs of sequences (Fig. 3A); while for the case of correlated sequences both distributions differ (Fig. 3 B). This result
opens the opportunity to discern whether a KS arrow AKS is significant, as one can computethe p-value of the KS arrow found for the real pair of sequences by calculating the amount ofprobability that is above |AKS | and below −|AKS | in the distribution of KS arrows coming fromthe reshufflings. This is the probability that a value higher than |AKS | or lower than −|AKS |comes from a pair of random sequences.

186

2.3

187

190

We applied the described leader-follower analysis methodology to presence data collected foracoustically-tagged reef manta rays (Mobula alfredi) at a cleaning station at a remote coral
reef in Seychelles. The resulting social network and reported leader-followers behavior werethen examined relative to the sex and size of tagged manta rays.

191

2.3.1

192

210

Presence data for reef manta rays were collected using passive acoustic telemetry, whereby asingle acoustic receiver (VR2W; VEMCO) was placed in close proximity to a known cleaning
station at D’Arros Island, Seychelles [21]. Twenty-five acoustic tags (n = 4, V12; n = 21,
V16-5H) were externally deployed on free-swimming reef manta rays using a using a modified Hawaiian sling between April 2013 and March 2014. Tags anchors were made of either
titanium (n = 4) or stainless steel (n = 21), and were positioned towards the posterior dorsal
surface of each individual. Prior to tag deployment, a photograph was taken of the unique andconsistent spot pattern present on the ventral surface of each manta ray [43, 49] to allow for
individuals to be identified in future monitoring surveys. The sex of each individual was determined by the presence (male) or absence (female) of claspers, and size (i.e. wingspan) was
visually estimated to the nearest 0.1 m. Individuals were classified into one of the following
two size groups based on their estimated size: small (1 female, 9 males), and big (9 females,
6 males).
Detection data from the cleaning station receiver were downloaded every six months,
and once a year the battery was replaced and the receiver inspected for damage or clock
drift. Preliminary range testing indicated a detection radius of approximately 150m (mean =(165 ± 33)m [27]) at the cleaning station receiver. Upon import of the detection data into a
Microsoft Access database, false positive detections were removed through filtration for active
tags and any receiver clock-drift time corrections were calculated assuming linear drift [27].

211

2.3.2

212

To construct the leader-follower network, we analyzed the interactions occurring between
pairs of tagged manta rays. Interactions were defined as the presence of one individual and aconsecutive appearance of the other. For statistical purposes, we considered only the pairs withmore than 50 interactions (n = 12 individuals). We then assessed value and direction of the
KS arrows between the distribution of lag times for each pair, as described above. We used a
global reshuffling scheme to compute the p value associated with each KS arrow, and discarded

173
174
175
176
177
178
179
180
181
182
183
184

188
189

193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209

213
214
215
216
217

A case study: reef manta rays (Mobula alfredi) at a cleaning station

Data collection

Leader-follower network

5

SciPost Physics

Submission

226

those for which p > 0.002. This reshuffling scheme allowed for the number of detections foreach individual to remain fixed, but for the timing of the events to be reshuffled from the poolof all detections. This differs from a local reshuffling scheme, which was also trialed and where
only the two sequences involved in the calculation of the arrow are involved in the reshuffling,
however, both schemes lead to similar results. In order to assess the correlation of sex andindividual size with the topology of the network, we compared two basic quantities: (1) the
number of edges present in the network of each type and (2) the average weight of the edges.
Comparisons were drawn across 100,000 reshufflings of the sexes and sizes, respectively, but
keeping the network structure fixed.

227

3

Results

228

3.1

Case study: reef manta rays at a cleaning station

229
230

A total of 41,607 detections of acoustically-tagged reef manta rays were recorded at the cleaning station receiver between April 2013 and January 2016.

231

3.2

232

240

The manta detection data were not distributed evenly in time, but rather in a quite peculiar
fashion. In Fig. 4A we show the distribution of events for mantas with more than 100 records(n = 18 individuals). Detections are most commonly recorded around noon, with an evident
circadian pattern (Fig. 4 B). The timing between consecutive events of single individuals,
hereafter ’interevent times’, are distributed following a heavy-tailed distribution, with a tailconsistent with a power law of exponent 1.38 (Fig. 4 C). The exponent has been fitted using
the ’powerlaw’ package available for python [50]. This implies that the temporal appearancesof the manta rays at the cleaning station are burst-like, with periods of high activity followed
by long times of inactivity and individual absences.

241

3.3

242

The final network constructed for reef manta rays at the cleaning station consisted of 12 individuals; five females (one small, four big) and seven males (3 small and 4 big). A total of
33 edges, representative of calculated KS arrows, were included in the network to generate asingle component (Fig. 5 A). Regarding the sexes, mixed edges are on average stronger than
expected, much more so for males following females, although there are many fewer males
following females than expected, in contrast to females following males, which are overrepresented. For same sex edges, the results reflect a randomized scenario. That is, while thenumber of edges present is as expected, their strength is smaller than expected (see Fig. 5 B).
As for the sizes, the most salient result is that the average weight of the interactions of small
individuals following small individuals is much lower than expected. The other types of edges
are only slightly off the values that are expected. So, small individuals following big individuals do so at a stronger frequency than expected and there is one link less than expected. For
big reef manta rays following small ones, the strength is slightly weaker than expected andthere is only one extra edge. For big individuals following other big conspecifics, the strengthof the interaction is slightly higher than expected, and there are only a couple more edges
than expected (see Fig. 5 C). In Fig. 6 we show the average appearance rate averaged over all
pairs. The appearance rate is computed for one individual as a function of the time since thelast presence event of the other individual, but restricting to pairs for which there is a connection in the inferred network. One can see that the appearance rate for followers is bigger than

218
219
220
221
222
223
224
225

233
234
235
236
237
238
239

243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260

Temporal heterogeneities

Leader-follower network.

6

SciPost Physics

Submission

264

for leaders up to a certain time lag of around 200 minutes. This enhanced appearance rateof the following individuals after the appearance of the individual that they follow confirms
our measurement of leader-follower behaviors and discards that we are not actually detecting
avoidance situations.

265

4

266

Leadership analyses represent a source of behavioral information. For example, the presenceor absence of leadership between adult and juvenile individuals provides insight into the learning strategies of a species, which can be biased towards learning from adults (i.e., juveniles
following them), or towards learning through trial and error (i.e., displaying independent behavior) [8]. Additionally, leadership analyses can be used to reveal long-term hierarchies inanimal groups and aid in characterizing the structure of populations [10].
Here, we introduced a statistical method for the inference of leadership patterns recorded
in presence data collected at a single location. The method relies on randomizations of the
input data, which provide estimates for the p-values associated with the leader-follower interactions found. We have tested this methodology against intuitive models of leading behavior,
and applied it to a real-world scenario involving presence data for reef manta rays at a cleaning station recorded using passive acoustic telemetry. Our method allowed us to constructa directed network from the manta detection data, which represents leader-follower interactions within the tagged population, and examine whether covariates such as sex and size affectthe position of the individuals in the network. Female manta rays were found to follow moremales than expected, while males followed fewer females, but with a stronger association.
With respect to size, smaller mantas followed other small mantas with a weaker interaction
strength, whereas the rest of interactions within the network appeared with a similar pattern
as that expected from the randomized case. Furthermore, we examined whether there is anenhanced probability of the appearance of a follower after the appearance of a leader at thecleaning station appeared at the location, and confirmed that the associations found within
the study are of the a leader-follower nature.
Comparison with other similar works. Ref. [34] examines the associations among reef manta
rays based on sighting data recorded using photo-identification techniques at five locations inIndonesia for five years. This presents a unique opportunity to compare sampling, network
analysis techniques, and findings between two populations of reef manta rays. While thephoto-identification data set of Ref. [34] spans a longer temporal range than that of the present
study, data collection was restricted to daylight hours only. In contrast, the automated datacollection facilitated by the acoustic receivers and tags of this study allowed for presence data
to be collected across the entire diel cycle. This allowed for a more detailed and dynamic
approach to the examination of leader-follower interactions in reef manta rays to be developed,
and made it possible to move beyond traditional associations defined using the “Gambit of theGroup” (i.e. assuming that all individuals observed together are associated). Note that thisdifference is crucial, since in our framework, two individuals, even if observed always together,might not be related by a leader-follower interaction depending on the timing of the presence
events.
Nevertheless, it is interesting to compare results between populations and datasets in termsof the influence of sex and size on network position. Perryman et al. found that male reefmanta rays tend not to associate with other males, and for them avoidance is common, whereas
females may be associated significantly with other females. The highest percentage of preferred dyadic associations was given between different sex pairs, however, there was partial
sexual segregation. Regarding maturity status, Perryman et al. found that juveniles tend to

261
262
263

267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307

Conclusion

7

SciPost Physics

308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340

341
342
343

344
345
346
347
348
349
350
351
352
353
354

Submission

associate in the short term with other juveniles and mature adults. It is challenging to determine whether the differences in findings between this, and the present study, are the result of
natural behavioral differences or are artifacts of differing analysis procedures. Should acoustic
data become available for reef manta rays at any of the five Indonesian field sites used by Perryman et al., future works should aim to replicate the leader-follower analyses outlined here
and to compare the results with the findings of the photo-identification-based study.
Ref. [32] describes a method for extracting social structure information from data similar
to that of the present study, but again, relies on the Gambit of the Group theory. The major drawback of this method (available within the package GMMevents) is that it relies onclustering the observations of assocations into different event windows, as observations occur
in bursts. In many situations, these events windows may be very difficult to define, as theinter-event times of the observations appear power-law distributed, without a clear cut forthe clusters, and it may be that association data is lost in the process. Ref. [14] builds upon
the methodology described in in Ref. [32], and allows for ‘leadership’ behavior to be inferred.
This is achieved by investigating which animal appears first in the association events that theclustering method has identified, however, by doing so, eliminates all of the information about
the dynamics contained in the timings of observations inside the event. Our method does notseparate times into classes (such as being part of an event), and instead, considers the complete temporal range of the data, making the analysis procedure less prone to ambiguities andinformation loss while allowing for leadership behaviors to be examined.
Limitations and future work. The analytical method proposed in this study allows only
for the detection of paired (i.e. dyadic) interactions from presence data collected at a single location, however, cannot currently be used to detect more complex, collective behaviors.
The expansion of these analysis to encompass data collected by multiple acoustic receivers
placed at different locations (i.e. acoustic arrays; e.g [21]) will allow for leadership networks
to be described using a multilayer approach [51],whereby each receiver is considered as adifferent layer. Such an approach may reveal if leadership patterns are coherent among different locations that may be significant for different ecological reasons (e.g. cleaning, feeding),
if followers continue to follow their leader across wider spatial scales, or if interactions are
more complex and followers in one location become leaders in another. For now, however,the present methodology can be applied to existing acoustic telemetry data sets around theglobe, and may reveal previously unidentified leader-follower behaviors and patterns in marine species across various ecosystem types.
Author contributions JFG, JPR, and VME designed the research. LRP provided acoustic
telemetry data. JFG, JPR, KK and VME developed the methods and analysed the results. All
authors contributed to the writing of the manuscript.
Funding information JFG acknowledges funding from the University of the Balearic Islands through its postdoc program and the postdoctoral program Vicenç Mut from the Direcció General de Política Universitària i Recerca from the government of the Balearic Islands.
KK acknowledges funding from MINECO through the Ramón y Cajal program and throughproject SPASIMM, FIS2016-80067-P (AEI/FEDER, EU). JFG, JPR and VME acknowledge funding from the Office of Sponsored Research at King Abdullah University of Science and Technology (KAUST) through project CAASE. Field research in Seychelles was approved by, andconducted with the knowledge of, the Seychelles Bureau of Standards, the Seychelles Ministryof Environment, Energy and Climate Change. Fieldwork and data collection was funded by
the Save Our Seas Foundation (SOSF) and supported by The Manta Trust and SOSF – D’Arros
Research Centre (SOSF-DRC).

8

SciPost Physics

Submission

A

B

C

Figure 1: (Color online) Two independent homogeneous Poisson processes of rate
λ = 1. The sequences were generated up to a maximum time of 100 time units. A
Raw event data. Individual A in blue and Individual B in red. B Cumulative distribution of lag times and Kolmogorov-Smirnof arrow, AKS , obtained from the sequences
shown above. In blue is the cumulative distribution of lag times of Individual A
following Individual B, while in red is its conjugate (lag times of B following A). C
Distribution of Kolmogorov-Smirnof arrows from 105 realizations of pairs generated
as in the sequences above. The distribution shows two peaks around 0, signature of
no leader-follower relation.

Figure 2: (Color online) Correlated time sequences. Individual A (blue) performs ahomogeneous Poisson process of rate λ1 = 1, while Individual B (red) follows thecorrelated non-homogeneous Poisson process described in the text with parametersλ∗ = 1, ∆t = 0.2 and δ = 4, i.e., it performs a Poisson process of rate λ2 = 1
except for 0.2 units of time after an event of Individual A, when it performs a
Poisson process of rate λ2 = 5. The sequences were generated up to a maximum
time of 100 time units. A Raw event data. B Cumulative distribution of lag times
and Kolmogorov-Smirnof arrow for the sequences shown above. C Distribution of
Kolmogorov-Smirnof arrows from 105 realizations of pairs generated as in the sequences above.

9

SciPost Physics

Submission

Figure 3: (Color online) Assessing significance of the leader-follower relation. A Random uncorrelated sequences with λ1 = λ2 = 1 and t max = 103 . In dashed black lines
the distribution of KS arrows for the ensemble of those sequences (104 independent
realizations). In green is the distribution of KS arrows for 104 reshufflings of the pairof sequences that gives rise to the arrow marked in red. B Correlated sequences with∆t = 0.2, δ = 4 and a maximum time of 103 time units. The KS arrow distribution
for these ensemble of sequences is shown in black (104 independent realizations),
while in green is the distribution of arrows coming from 104 reshufflings of the sequence pair that gives rise originally to the value of the arrow signaled at the red
line.

Figure 4: (Color online) Temporal heterogeneities in the data. A Raw event data for
acoustically-tagged reef manta rays (Mobula alfredi) with more than 100 detection
events at a receiver located near a cleaning station in Seychelles. B Appearance
probability as a function of the hour in the day. The cyan lines correspond to different
males, while the pink lines correspond to different females. The red and blue lines
are the averages for females and males respectively. The black curve corresponds tothe appearance probability of all individuals pooled. We can see that the appearanceof reef manta rays at the cleaning station is concentrated around noon. C Interevent
times distribution, i.e., the distribution of times between consecutive presence eventsof the same individual. The long tail, compatible with a power law of exponent 1.38
makes clear the burst-like presence of manta rays at the cleaning station. The peak inthe distribution is located at a time of one day, which is a reflection of the circadianpattern (i.e. 24 hr) of the manta rays.

10

SciPost Physics

Submission

Figure 5: (Color online) Leadership network for acoustically-tagged reef manta rays
(Mobula alfredi) detected by a single acoustic receiver placed at a cleaning station
at a remote coral reef in Seychelles relative to the sex (Female, pink; Male, blue)and size (node size) of individuals. A Leader-follower network of manta rays at aconfidence value p = 0.002. B and C Types of edges depending on sex and size. Anedge of type xy stands for an individual of type x following an individual of type y,where x and y can be f (female) or m (male) for the sexes, or s (small) and B (big)for the sizes. The black symbols are the expected values from 104 reshuffling of thesexes/sizes of the individuals. The x-axis shows the number of such edges, while they-axis shows the average strength of the edges. Deviation of the real data in the x
direction indicates a difference of the number of edges of that kind found in the real
network. Deviation in the y direction signals a difference with the expected strengthof the leader-follower relation. The arrows show how the results from the original
network differ from the results of randomizations.Results for the randomization ofsexes (B). Mixed edges are on average stronger than expected, and there are muchless males following females than expected, in contrast to females following males,
which are overrepresented. For same sex edges the results are similar, the number ofedges is approximately the one expected, while the strength is smaller than expected.
Results for the randomization of sizes (C). The most salient result is that the average
weight of the interactions of males following males is much lower than expected.
The other types of edges are only slightly off the values that are expected. So small
individuals following big ones follow a bit stronger than expected and there is one
link less than expected. For big ones following small ones, the strength is slightly
weaker than expected and there is only one extra edge. For big ones following big
ones the strength is slightly bigger than expected and there are only a couple moreedges than expected.

11

SciPost Physics

Submission

Figure 6: (Color online) Appearance rate for the follower right after the appearanceof the leader (blue) or for the leader after the appearance of the follower (red). The
rates show that the follower actually has a higher rate of appearing right after an
event of the leader. The curves show the average for all the pairs present in the
network. The dashed lines show the standard errors.

355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378

References[1] J. Ren, W. Sun, D. Manocha, A. Li and X. Jin, Stable information transfer network facilitates
the emergence of collective behavior of bird flocks, Physical Review E 98(5), 1 (2018),
doi:10.1103/PhysRevE.98.052309.
[2] T. Gernat, V. D. Rao, M. Middendorf, H. Dankowicz, N. Goldenfeld and G. E. Robinson,
Automated monitoring of behavior reveals bursty interaction patterns and rapid spreading
dynamics in honeybee social networks, Proceedings of the National Academy of Sciencesof the United States of America 115(7), 1433 (2018), doi:10.1073/pnas.1713568115.
[3] S. Chen, B. J. White, M. W. Sanderson, D. E. Amrine, A. Ilany and C. Lanzas, Highly dynamic animal contact network and implications on disease transmission, Scientific Reports
4, 1 (2014), doi:10.1038/srep04472.
[4] I. D. Couzin, J. Krause, R. James, G. D. Ruxton and N. R. Franks, Collective memoryand spatial sorting in animal groups, Journal of Theoretical Biology 218(1), 1 (2002),
doi:10.1006/jtbi.2002.3065.
[5] I. D. Couzin and J. Krause, Self-Organization and Collective Behavior in Vertebrates, Advances in the Study of Behavior 32, 1 (2003), doi:10.1016/S0065-3454(03)01001-5.
[6] D. J. Sumpter,
The principles of collective animal behaviour,
PhilosophicalTransactions of the Royal Society B: Biological Sciences 361(1465), 5 (2006),
doi:10.1098/rstb.2005.1733.
[7] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte,
A. Orlandi, G. Parisi, A. Procaccini, M. Viale and V. Zdravkovic, Interaction ruling animalcollective behavior depends on topological rather than metric distance: Evidence from a fieldstudy, Proceedings of the National Academy of Sciences of the United States of America
105(4), 1232 (2008), doi:10.1073/pnas.0711437105.

12

SciPost Physics

379
380
381
382
383
384
385
386
387
388
389
390

Submission[8] B. G. GALEF and K. N. LALAND,
Social Learning in Animals: Empirical Studies and Theoretical Models,
BioScience 55(6), 489 (2005),
doi:10.1641/00063568(2005)055[0489:sliaes]2.0.co;2.
[9] A. J. W. Ward, J. E. Herbert-Read, D. J. T. Sumpter and J. Krause, Fast and accurate
decisions through collective vigilance in fish shoals, Proc. Natl. Acad. Sci. USA 108, 2312
(2011).
[10] A. Sih, S. F. Hanser and K. A. McHugh, Social network theory: New insights and issues for behavioral ecologists, Behavioral Ecology and Sociobiology 63(7), 975 (2009),
doi:10.1007/s00265-009-0725-6.
[11] J. Krause, R. James and D. P. Croft, Personality in the context of social networks, Philosophical Transactions of the Royal Society B: Biological Sciences 365(1560), 4099 (2010),
doi:10.1098/rstb.2010.0216.

393[12] H. Ohtsuki, C. Hauert, E. Lieberman and M. A. Nowak, A simple rule for the evolution of cooperation on graphs and social networks, Nature 441(7092), 502 (2006),
doi:10.1038/nature04605.

394[13] M. Newman, Networks: An Introduction, OUP Oxford, ISBN 9780199206650 (2010).

395[14] D. M. Jacoby and R. Freeman, Emerging Network-Based Tools in Movement Ecology, Trends
in Ecology and Evolution 31(4), 301 (2016), doi:10.1016/j.tree.2016.01.011.

391
392

396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420

[15] J. Krause, D. Lusseau and R. James, Animal social networks: An introduction, Behavioral
Ecology and Sociobiology 63(7), 967 (2009), doi:10.1007/s00265-009-0747-0.
[16] V. D. Blondel, A. Decuyper and G. Krings, A survey of results on mobile phone datasets
analysis, EPJ Data Science 4(1), 1 (2015), doi:10.1140/epjds/s13688-015-0046-0,
1502.03406.
[17] J. Krause, S. Krause, R. Arlinghaus, I. Psorakis, S. Roberts and C. Rutz, Reality
mining of animal social systems., Trends in ecology & evolution 28(9), 541 (2013),
doi:10.1016/j.tree.2013.06.002.
[18] E. S. Germanov and A. D. Marshall, Running the Gauntlet: Regional movement patternsof Manta alfredi through a complex of parks and fisheries, PLoS ONE 9(10), 1 (2014),
doi:10.1371/journal.pone.0110071.
[19] B. Albanese, P. L. Angermeier and C. Gowan, Designing Mark–Recapture Studies to Reduce
Effects of Distance Weighting on Movement Distance Distributions of Stream Fishes, Transactions of the American Fisheries Society 132(5), 925 (2003), doi:10.1577/t03-019.
[20] L. W. Cameron, W. Roche, P. Green, J. D. Houghton and P. J. Mensink, Transatlantic
movement in porbeagle sharks, Lamna nasus, Fisheries Research 207(June), 25 (2018),
doi:10.1016/j.fishres.2018.05.014.
[21] L. Peel, G. Stevens, R. Daly, C. Keating Daly, J. Lea, C. Clarke, S. Collin and M. Meekan,
Movement and residency patterns of reef manta rays mobula alfredi in the amirante islands,
seychelles, Mar Ecol Prog Ser 621, 169 (2019), doi:10.3354/meps12995.
[22] R. B. McAuley, B. D. Bruce, I. S. Keay, S. Mountford, T. Pinnell and F. G. Whoriskey,
Broad-scale coastal movements of white sharks off Western Australia described by passive acoustic telemetry data, Marine and Freshwater Research 68(8), 1518 (2017),
doi:10.1071/MF16222.
13

SciPost Physics

421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463

Submission[23] J. C. Carrier and C. A. Simpfendorfer, Shark Research: Emerging Technologies and Applications for the Field and Laboratory, CRC Marine Biology Series. CRC Press, Taylor &Francis Group, ISBN 9781138032927 (2018).
[24] A. Steckenreuter, X. Hoenner, C. Huveneers, C. Simpfendorfer, M. J. Buscot, K. Tattersall,
R. Babcock, M. Heupel, M. Meekan, J. Van Den Broek, P. McDowall, V. Peddemors et al.,Optimising the design of large-scale acoustic telemetry curtains, Marine and Freshwater
Research 68(8), 1403 (2017), doi:10.1071/MF16126.
[25] G. T. Crossin, M. R. Heupel, C. M. Holbrook, N. E. Hussey, S. K. Lowerre-Barbieri, V. M.
Nguyen, G. D. Raby and S. J. Cooke, Acoustic telemetry and fisheries management, Ecological Applications 27(4), 1031 (2017), doi:10.1002/eap.1533.
[26] J. D. Stewart, F. R. Jaine, A. J. Armstrong, A. O. Armstrong, M. B. Bennett, K. B. Burgess,L. I. Couturier, D. A. Croll, M. R. Cronin, M. H. Deakos, C. L. Dudgeon, D. Fernando
et al., Research priorities to support effective Manta and Devil Ray conservation, Frontiersin Marine Science 5(SEP), 1 (2018), doi:10.3389/fmars.2018.00314.
[27] J. S. Lea, N. E. Humphries, R. G. von Brandis, C. R. Clarke and D. W. Sims, Acoustic
telemetry and network analysis reveal the space use of multiple reef predators and enhancemarine protected area design, Proceedings of the Royal Society B: Biological Sciences283(1834) (2016), doi:10.1098/rspb.2016.0717.
[28] A. J. King, D. D. Johnson and M. Van Vugt, The Origins and Evolution of Leadership,
Current Biology 19(19), R911 (2009), doi:10.1016/j.cub.2009.07.027.
[29] A. Flack, M. Nagy, W. Fiedler, I. D. Couzin and M. Wikelski, From local collective behavior to global migratory patterns in white storks, Science 360(6391), 911 (2018),
doi:10.1126/science.aap7781.
[30] J. Garland, A. M. Berdahl, J. Sun and E. M. Bollt, Anatomy of leadership in collective
behaviour, Chaos 28(7) (2018), doi:10.1063/1.5024395, 1802.01194.
[31] B. Corominas-Murtra, J. Goñi, R. V. Solé and C. Rodríguez-Caso, On the origins of hierarchy in complex networks, Proceedings of the National Academy of Sciences of the UnitedStates of America 110(33), 13316 (2013), doi:10.1073/pnas.1300832110, 1303.2503.
[32] I. Psorakis, S. J. Roberts, I. Rezek and B. C. Sheldon, Inferring social network structure in
ecological systems from spatiotemporal data streams, Journal of the Royal Society Interface
9(76), 3055 (2012), doi:10.1098/rsif.2012.0223.
[33] I. Psorakis, B. Voelkl, C. J. Garroway, R. Radersma, L. M. Aplin, R. A. Crates, A. Culina,
D. R. Farine, J. A. Firth, C. A. Hinde, L. R. Kidd, N. D. Milligan et al., Inferring social
structure from temporal data, Behavioral Ecology and Sociobiology 69(5), 857 (2015),
doi:10.1007/s00265-015-1906-0.
[34] R. J. Perryman, S. K. Venables, R. F. Tapilatu, A. D. Marshall, C. Brown and D. W. Franks,Social preferences and network structure in a population of reef manta rays, Behavioral
Ecology and Sociobiology 73(8) (2019), doi:10.1007/s00265-019-2720-x.
[35] H. W. Lauw, E. P. Lim, H. Pang and T. T. Tan, Social network discovery by mining
spatio-temporal events, Computational and Mathematical Organization Theory 11(2),97 (2005), doi:10.1007/s10588-005-3939-9.
[36] H. Whitehead, Analyzing Animal Societies: Quantitative Methods for Vertebrate Social
Analysis, University of Chicago Press, ISBN 9780226895246 (2008).
14

SciPost Physics

464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505

Submission[37] K. P. Oh and A. V. Badyaev, Structure of social networks in a passerine bird: Consequences
for sexual selection and the evolution of mating strategies, American Naturalist 176(3)(2010), doi:10.1086/655216.
[38] S. Gero, D. Engelhaupt, L. Rendell and H. Whitehead, Who cares? Between-group variation in alloparental caregiving in sperm whales, Behavioral Ecology 20(4), 838 (2009),
doi:10.1093/beheco/arp068.
[39] G. Krings, M. Karsai, S. Bernhardsson, V. D. Blondel and J. Saramäki, Effects of time
window size and placement on the structure of an aggregated communication network, EPJ
Data Science 1(1), 1 (2012), doi:10.1140/epjds4.
[40] L. M. Aplin, D. R. Farine, J. Morand-Ferron, E. F. Cole, A. Cockburn and B. C. Sheldon,
Individual personalities predict social behaviour in wild networks of great tits (Parus major),
Ecology Letters 16(11), 1365 (2013), doi:10.1111/ele.12181.
[41] L. J. Brent, S. R. Heilbronner, J. E. Horvath, J. Gonzalez-Martinez, A. Ruiz-Lambides,
A. G. Robinson, J. H. Pate Skene and M. L. Platt, Genetic origins of social networks inrhesus macaques, Scientific Reports 3, 1 (2013), doi:10.1038/srep01042.
[42] G. M. Stevens, J. P. Hawkins and C. M. Roberts, Courtship and mating behaviour of manta
rays Mobula alfredi and M. birostris in the Maldives, Journal of Fish Biology 93(2), 344(2018), doi:10.1111/jfb.13768.
[43] A. D. Marshall, C. L. Dudgeon and M. B. Bennett, Size and structure of a photographically
identified population of manta rays Manta alfredi in southern Mozambique, Marine Biology
158(5), 1111 (2011), doi:10.1007/s00227-011-1634-6.
[44] J. Fernández-Gracia, Leadership package,leadership_package (2020).

https://github.com/juanfernandezgracia/[45] A. N. Kolmogorov and A. N. Shiryaev, Selected Works II: Probability Theory and Mathematical Statistics, Mathematics and its Applications. Springer Netherlands, ISBN
9789401122603 (2012).
[46] M. Statistics, On the Comparison of Two Empirical Distribution Functions Author ( s ):Lajos Takacs Source : The Annals of Mathematical Statistics , Vol . 42 , No . 4 ( Aug ., 1971), pp . 1157-1166 Published by : Institute of Mathematical Statistics Stable URL : https:
42(4), 1157 (2019).
[47] A. G. Hawkes, Spectra of Some Self-Exciting and Mutually Exciting Point Processes,
Biometrika 58(1), 83 (1971), doi:10.2307/2334319.
[48] R. D. Malmgren, D. B. Stouffer, A. S. Campanharo and L. A. Amaral,
On
universality in human correspondence activity,
Science 325(5948), 1696 (2009),
doi:10.1126/science.1174562, 0909.4724.
[49] A. D. Marshall, L. J. V. Compagno and M. B. Bennett, <strong>Redescription of thegenus <em>Manta </em>with resurrection of <em>Manta alfredi </em>(Krefft, 1868)(Chondrichthyes; Myliobatoidei; Mobulidae)</strong>, Zootaxa 2301(1), 1 (2009),
doi:10.11646/zootaxa.2301.1.1.
[50] J. Alstott, E. Bullmore and D. Plenz,
Analysis of Heavy-Tailed Distributions,doi:10.1371/journal.pone.0085777.
15

powerlaw:
A Python Package for
PLoS ONE 9(1), e85777 (2014),

SciPost Physics

506
507
508

Submission[51] J. Mourier, E. J. I. Lédée and D. M. P. Jacoby, A multilayer perspective for inferring spatial and social functioning in animal movement networks, bioRxiv p. 749085 (2019),
doi:10.1101/749085.

16