...

Collective Sensing and Information Transfer in Animal Groups

by user

on
Category: Documents
3

views

Report

Comments

Transcript

Collective Sensing and Information Transfer in Animal Groups
Collective Sensing and Information
Transfer in Animal Groups
Sara Brinton Rosenthal
A Dissertation
Presented to the Faculty
of Princeton University
in Candidacy for the Degree
of Doctor of Philosophy
Recommended for Acceptance
by the Department of
Physics
Adviser: Iain D. Couzin
March 2015
c Copyright by Sara Brinton Rosenthal, 2015.
All rights reserved.
Abstract
This dissertation addresses several topics in collective motion in animal groups.
Coordination among social animals often requires rapid and efficient transfer of information among individuals, which may depend crucially on the underlying structure
of the interaction network used for communication. In two experimental systems and
in one simulation study, we study the nature of interactions and interaction networks,
and how these interactions scale up to global order in the system. The resulting collective properties allow animal groups access to “collective computation” such that
they can process and respond to stimuli in ways in which a single individual cannot.
First, we study collective evasion maneuvers, manifested through rapid cascades of
behavioral change (a ubiquitous behavior among taxa), in schooling fish (Notemigonus
crysoleucas). We automatically track the positions, and calculate visual fields of all
individuals in schools of approximately 150 fish, and determine the functional mapping between socially generated sensory input and motor response during collective
evasion. We find that individuals employ simple, robust measures to assess behavioral
changes in neighbors, and that the resulting networks by which behavior propagates
throughout groups are complex; being weighted and directed. By studying these
interaction networks, we reveal the (complex, fractional) nature of social contagion,
and establish that individuals with relatively few, but strongly-connected, neighbors
are both most socially influential, and most susceptible to influence.
Next, we study the relationship between emergent periodic synchronization in ant
colonies, and interactions between ants in different behavioral states. We investigate
the factors driving fluctuations in the overall level of synchronization observed in the
colony, and find that flexible behavioral responses to interactions can explain these
fluctuations, in both real data and in a simulated ant colony.
Finally, we use a model of collective movement, motivated by experimental data,
where agents are coupled by speed to their environment to demonstrate that simple
iii
behavioral rules can allow the group to maximize performance in a dynamical search
task. Additionally, within the constraints of the model, the behaviors that optimize
performance place the population near a state transition. The model reveals that
individuals readily locate and track dynamic resources by splitting and fusing to
form groups that match the length scale of these resources. This occurs even when
individuals have no way of evaluating the sizes of resources or determining the sizes of
groups to which they belong. Our model demonstrates that fission-fusion dynamics
can allow social animals to balance the exploration-exploitation tradeoff, to best take
advantage of dynamic resources.
iv
Acknowledgements
I would first like to express my deepest thanks to Iain Couzin, my advisor, who is
continually inspiring in his passion for science. I always came out of meetings with
Iain with many more ideas than I came in with, and with much more excitement for
that matter. He has instilled in me the importance of finding the story hidden in the
science.
I would also like to extend thanks to the members of my thesis committee- William
Bialek and Frank Calaprice, as well as to Josh Shaevitz for reading my thesis. I am
honored to have the chance to work with such brilliant minds. Thanks also to Thomas
Gregor, who first introduced me to the world of collective behavior. Thanks to Herschel Rabitz, who got me started on my first research project, as an undergraduate.
Thanks to Naomi Leonard, Will Happer, and Ned Wingreen for their valuable input and advice. Thanks to the physics department for the flexibility which allowed
me to find and pursue my true interests, even though they turned out to be in a
non-traditional area of research.
None of the projects discussed here would have been possible without many collaborators. Thanks to Colin Twomey and Andrew Hartnett, who were my two deskneighbors as well as collaborators, and who were never too busy to give thoughtful
answers to my many questions. Thanks to Andrew Hein and George Hagstrom; the
collective sensing project has truly been a rewarding collaboration. Thanks to Catherine Offord, who has shared her vast knowledge of the ant world with me. And thanks
to everyone else in the Couzin Lab, for many valuable discussions, and for making
daily life fun and intellectually engaging.
My family has been enormously supportive of my grad school career. Thanks to
my mom, Abby Rosenthal, for all of the encouragement (and for help with taxes every
year!). Blake, my sister, thank you for always being there for me, through the good
v
times and bad. Thank you to my aunts, Sara Rosenthal and Julie Prazich, who have
always been so proud and supportive of our academic achievements.
Of course I could not have made it this far without the love and support of my
future husband, Thomas Weber, who can always find a way to make me laugh. Thank
you for everything Tom!
Five plus years at Princeton has really flown by, and it would not have been nearly
as enjoyable without all the amazing friendships made here. Thanks to my friends and
house-mates over the years- Bridget Berglind, Jessica Shang, Pip Coen, Ben Safdi,
Loren Alegria, Cynthia Gerlein-Safdi, Diego Pacheco, and Ann Duan. Thanks also
to Gretchen Wendel, with whom a great friendship transcends time and distance.
And finally, I would like to thank my dad, Carey Rosenthal, for teaching me to
love math and science, for teaching me to think about a simple problem before a
complicated one, and for writing math problems as notes in my lunch box as a kid.
vi
To my dad, Carey Rosenthal.
vii
Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
iii
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
v
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
xi
1 Introduction
1
2 Revealing the hidden networks of communication in schooling fish
6
2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.2
Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.3
2.2.1
Functional mapping between sensory input and behavioral change 15
2.2.2
Reconstruction and analysis of the quantitative interaction network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
Technical details . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
2.3.A Experimental methods . . . . . . . . . . . . . . . . . . . . . .
20
2.3.B Stereotypy of startle response in golden shiners . . . . . . . .
21
2.3.C Description of tracking software . . . . . . . . . . . . . . . . .
30
2.3.D Description of vision software . . . . . . . . . . . . . . . . . .
32
2.3.E Features used in model selection . . . . . . . . . . . . . . . . .
34
2.3.F Multi-model inference . . . . . . . . . . . . . . . . . . . . . .
36
2.3.G L1-penalized logistic regression for feature selection . . . . . .
42
3 Predicting complex behavioral contagion
viii
46
3.1
Predicting behavioral cascades and the relationship between spatial
position and social influence in groups . . . . . . . . . . . . . . . . .
46
3.2
The nature of the contagious process . . . . . . . . . . . . . . . . . .
51
3.3
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
3.4
Technical details . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
3.4.A Group spatial measures . . . . . . . . . . . . . . . . . . . . . .
55
3.4.B Spatial structure of the network . . . . . . . . . . . . . . . . .
57
3.4.C Exterior visual field . . . . . . . . . . . . . . . . . . . . . . . .
61
3.4.D Fitting local clustering coefficient model of cascade size . . . .
62
3.4.E Susceptibility . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
3.4.F Simple, numeric, and fractional model simulations . . . . . . .
70
4 Collective Oscillations in Leptothorax acervorum
86
4.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
86
4.2
Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
89
4.3
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
95
4.3.1
Role of interactions-model fitting . . . . . . . . . . . . . . . .
95
4.3.2
Responsiveness over time . . . . . . . . . . . . . . . . . . . . . 100
4.3.3
Simulating collective oscillations . . . . . . . . . . . . . . . . . 103
4.4
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
4.5
Technical details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
4.5.A Modifications to the MCA model . . . . . . . . . . . . . . . . 107
4.5.B Best-fit value of g . . . . . . . . . . . . . . . . . . . . . . . . . 110
4.5.C Colonies without brood . . . . . . . . . . . . . . . . . . . . . . 111
5 Exploration and exploitation in dynamic environments
112
5.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
5.2
Model development . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
ix
5.3
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
5.3.1
Social interactions help individuals exploit resource peaks . . . 123
5.3.2
Individuals in large groups fail to explore the environment . . 125
5.3.3
Individuals maximize performance by balancing exploration
and exploitation . . . . . . . . . . . . . . . . . . . . . . . . . . 128
5.3.4
Populations of optimal individuals exhibit fission-fusion dynamics128
5.3.5
Optimal individual behavior places populations near a state
transition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
5.4
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
5.5
Technical details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
5.5.A Tracking dynamic signals . . . . . . . . . . . . . . . . . . . . . 136
5.5.B Alternative peak shapes . . . . . . . . . . . . . . . . . . . . . 138
5.5.C Sensitivity of optimum value . . . . . . . . . . . . . . . . . . . 140
5.5.D Resource depletion . . . . . . . . . . . . . . . . . . . . . . . . 147
5.5.E Evolutionary Algorithm . . . . . . . . . . . . . . . . . . . . . 149
5.5.F Theoretical observations behind phase transitions in model and
balance of exploration with exploitation . . . . . . . . . . . . 153
6 Conclusions and outlook
164
Bibliography
167
x
List of Figures
2.1
Collective evasion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.2
Comparison of triggered and spontaneous behavioral cascades . . . .
14
2.3
Features used in model selection . . . . . . . . . . . . . . . . . . . . .
16
2.4
The sensory basis of social contagion . . . . . . . . . . . . . . . . . .
18
2.5
Interaction networks are complex, with weighted, directed pathways of
communication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
S2.1 Spontaneous startle characteristics . . . . . . . . . . . . . . . . . . .
24
S2.2 Spatial distributions of initiators, first-responders, and non-responders
26
S2.3 Experimental setup for triggered alarms . . . . . . . . . . . . . . . .
29
S2.4 Comparison of raw and log-transformed features . . . . . . . . . . . .
35
S2.5 Feature transforms . . . . . . . . . . . . . . . . . . . . . . . . . . . .
36
S2.6 Results from multi-model inference . . . . . . . . . . . . . . . . . . .
38
S2.7 Feature pair conditional probabilities . . . . . . . . . . . . . . . . . .
40
S2.8 Sensory network edge weight distribution . . . . . . . . . . . . . . . .
41
S2.9 Result from L1 Regularization . . . . . . . . . . . . . . . . . . . . . .
43
S2.10Fidelity in spatial position . . . . . . . . . . . . . . . . . . . . . . . .
45
3.1
Local weighted clustering coefficient predicts the magnitude of behavioral cascades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
S3.1 Spatial measures . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
S3.2 Network statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58
xi
S3.3 Spatial structure of influence and local density . . . . . . . . . . . . .
60
S3.4 Spatial correlation of clustering coefficient . . . . . . . . . . . . . . .
61
S3.5 Exterior visual field . . . . . . . . . . . . . . . . . . . . . . . . . . . .
62
S3.6 Variance explained by clustering coefficient . . . . . . . . . . . . . . .
64
S3.7 Alternative network measures . . . . . . . . . . . . . . . . . . . . . .
68
S3.8 Behavioral cascade cross validation . . . . . . . . . . . . . . . . . . .
70
S3.9 Log likelihood for numeric and fractional models . . . . . . . . . . . .
73
S3.10Maximal and marginal likelihoods for contagion models . . . . . . . .
75
S3.11Likelihoods as functions of response threshold . . . . . . . . . . . . .
77
S3.12Best fit parameters for contagion models . . . . . . . . . . . . . . . .
80
S3.13Likelihoods as functions of response threshold, normal distribution . .
81
S3.14Contagion model results, normal distribution . . . . . . . . . . . . . .
83
S3.15Sensitivity of simulated cascade size to width of response distribution
85
4.1
Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . .
90
4.2
Ants transition profile . . . . . . . . . . . . . . . . . . . . . . . . . .
91
4.3
Activity waves in ant colony . . . . . . . . . . . . . . . . . . . . . . .
93
4.4
Activity waves in ant colony, magnified . . . . . . . . . . . . . . . . .
94
4.5
Coherence of activity waves in 3 colonies . . . . . . . . . . . . . . . .
95
4.6
Frequency of interaction types . . . . . . . . . . . . . . . . . . . . . .
96
4.7
Relative transition probabilities . . . . . . . . . . . . . . . . . . . . .
99
4.8
Colony inflow rate
4.9
Coherence and interaction coefficients . . . . . . . . . . . . . . . . . . 101
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.10 Relative contribution to variation in coherence . . . . . . . . . . . . . 102
4.11 Spatial fidelity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.12 Activity waves in a simulated colony . . . . . . . . . . . . . . . . . . 104
4.13 Coherence and coupling coefficients, simulations . . . . . . . . . . . . 106
S4.1 Speeds of interacting ants . . . . . . . . . . . . . . . . . . . . . . . . 108
xii
S4.2 Best-fit gain parameter . . . . . . . . . . . . . . . . . . . . . . . . . . 110
S4.3 Homing and coherence . . . . . . . . . . . . . . . . . . . . . . . . . . 111
5.1
Critical angle for capture by peak . . . . . . . . . . . . . . . . . . . . 121
5.2
Performance and resource tracking by social agents . . . . . . . . . . 124
5.3
Group size and performance in a two-peak environment . . . . . . . . 127
5.4
Peak size matching and group fission-fusion . . . . . . . . . . . . . . 130
5.5
Optimal behavioral rules occur near a critical point . . . . . . . . . . 133
5.6
Default parameters used in simulations . . . . . . . . . . . . . . . . . 136
S5.1 Noisy resource tracking . . . . . . . . . . . . . . . . . . . . . . . . . . 137
S5.2 Exploration/exploitation with Gaussian resource peaks . . . . . . . . 139
S5.3 Sensitivity to peak mobility . . . . . . . . . . . . . . . . . . . . . . . 141
S5.4 Sensitivity to peak size . . . . . . . . . . . . . . . . . . . . . . . . . . 142
S5.5 Sensitivity to number of peaks . . . . . . . . . . . . . . . . . . . . . . 143
S5.6 Sensitivity to number of agents . . . . . . . . . . . . . . . . . . . . . 144
S5.7 Sensitivity to ψ1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
S5.8 Explore/exploit with depleting resources . . . . . . . . . . . . . . . . 149
S5.9 Evolving the optimal ψ0 . . . . . . . . . . . . . . . . . . . . . . . . . 152
S5.10Potential energy and group radius . . . . . . . . . . . . . . . . . . . . 156
S5.11Hysteresis in state transitions . . . . . . . . . . . . . . . . . . . . . . 157
S5.12Agent entering region where it feels peak influence . . . . . . . . . . . 159
S5.13Solutions for critical angle, ∆c . . . . . . . . . . . . . . . . . . . . . . 163
xiii
Chapter 1
Introduction
In physics, we strive to reduce complex phenomena to constituent elements, to explain the world in which we live. Such efforts by great minds over the past century have brought about revolution upon revolution in our understanding of how
the world around us works. From the theory of relativity to the discovery of the
Higgs boson, physics has triumphed in connecting theoretical ideas with experimental evidence. Recently, the frontier of biological systems has been tackled by
physical theorists and experimentalists alike, where much success has been found,
but where the complexity of living organisms has often stymied the desire for clean
and elegant solutions. The problem of complexity is even more apparent when considering not just single organisms, but systems of interacting animals.
Be that
as it may, there have been clues indicating that the complex global phenomena
and efficient communication we see in coordinated animal groups can be generated from simple, local, individual-level interactions, without any centralized control (to name just a few examples; human societies [Castellano et al., 2009], flocking starlings [Ballerini et al., 2008], schooling fish [Berdahl et al., 2013], social insects
[Mersch et al., 2013], and bacterial colonies [Miller and Bassler, 2001]). These global
behaviors can be highly synchronized [Mirollo and Strogatz, 1990, Strogatz, 2003,
1
Boi et al., 1999], and appear to be robust to error while remaining highly sensitive to
relevant stimuli [Couzin, 2007], even in the absence of leadership. Reviews of experimental and theoretical findings in collective motion of mobile animal groups may be
found in [Couzin and Krause, 2003] and [Vicsek and Zafeiris, 2012].
Physicist have long been captivated by the concept of emergence- of local interactions leading to global order- beginning with the ubiquitous Ising model, first
discussed in [Lenz, 1920] and [Ising, 1925], a mathematical model of a system of
magnetic dipoles on a lattice each interacting with a set of near neighbors. In twodimensions this model identifies a critical temperature at which there is an abrupt
transition between order, where all dipoles have identical spins, and disorder, where
there is no spatial correlation between spins [Onsager, 1944]. Additionally, these
Ising models, and the related principle of maximum entropy [Jaynes, 1957], which
help in fitting Ising model parameters to data, have found varied uses in applications to biological systems, such as in analyzing collective properties of networks of
neurons [Hopfield, 1982, Tkacik et al., 2006] (see [Mora and Bialek, 2011] for a review). A recurrent theme in recent work (reviewed in [Mora and Bialek, 2011], and
[Chialvo, 2010]) is that emergent systems appear to be poised at a critical point,
balancing between states of oversensitivity and undersensitivity.
Unlike in networks of neurons, where advances in technology have allowed for
direct recording and measurement of physical and functional connections, in animal
collectives it is impossible to measure directly the social connection between any
two individuals. However, having a good estimate of the communication network
among organisms is critical for understanding how collective properties emerge from
local interactions. A number of interaction metrics have been used as proxies in
place of direct measurements, such as assuming that individuals interact with all
neighbors within a certain metric (Euclidean) distance, or topological (N nearest
neighbors) distance. One study finds evidence that topological distance interactions
2
rather than metric distance interactions better explain the observed behavior of starling flocks [Ballerini et al., 2008]. However, new advances in imaging and tracking
technology have allowed for more precise measurements of individual motion, and
have thus made possible a more biologically motivated estimate of interaction connections between individuals (for example, in a recent study the authors find evidence that information transfer among individuals in fish schools is visually mediated
[Strandburg-Peshkin et al., 2013]). In Chapter 2, experimental data is used to estimate one such interaction network during collective evasion in freely schooling fish,
by mapping sensory input to behavioral output. Then in Chapter 3, this network
is used to understand how waves of alarm spread through groups. We predict the
magnitudes of these behavioral cascades, and identify the most socially influential
individuals. We also reveal the complex and fractional nature of social contagion,
where the most socially influential individuals are those which have few, but strongly
connected, neighbors.
A different type of interacting organism is considered in Chapter 4; we study the
relationship between emergent activity synchronization in ant colonies (Leptothorax
acervorum), and coupling between interacting ants. Similar to many types of neurons, individual ants do not display periodic activity cycles when isolated, but when
they are grouped in a colony, activity synchronization emerges. In continuing the
analogy with neural systems, where active neurons are able excite those inactive neurons they are connected to, active ants have the ability to excite those stationary ants
they interact with, causing them to become active. Rhythmic pulses in the brain are
known to be beneficial for a number of reasons, and while it is not known exactly what
function these synchronized activity waves play in ant colonies, it has been hypothesized that they allow for coordinated foraging bouts, and simulations have shown
that coordinated activity can result in greater task efficiency [Hatcher et al., 1992].
In our data we find that behavioral responses to interactions are flexible, and that this
3
response flexibility plays a role in moderating the level of synchronization in colony
activity.
In contrast to projects using experimental data, many modeling approaches start
with a set of behavioral and interaction rules, and evaluate the resulting collective
properties. In studies of mobile animal groups, simple models of collective motion
have shown that somewhat realistic global properties of collective behavior can be
generated from local interactions, where these models are composed of a population
of agents each of which has a certain autonomous self-propulsion. One of the first
such ‘self-propelled particle’ models is found in [Vicsek et al., 1995], where the authors find that agents (moving with constant velocity) move in a disordered state if
the level of noise is high, but that velocities align and motion is ordered once the level
of noise is low enough (they draw the analogy between their model and a model of
ferromagnetic interaction, where alignment of spins is represented by agents aligning
their directions of motion, and temperature is represented by noise). Since this paper,
there have been many variations on this model, with varying degrees of biological realism (e.g. [Couzin et al., 2002, D’Orsogna et al., 2006], and see [Czirók et al., 1999]
for a review). These models are useful in simulating collective behavior using only a
minimal set of rules. In Chapter 5 I present a simulation study using a self-propelled
particle model, motivated by experimental data from two recent studies of fish schooling [Katz et al., 2011, Berdahl et al., 2013], in which agents exert repulsive forces on
one another when they are closer than a certain distance, and exert attractive forces
outside this radius of repulsion [Katz et al., 2011]. This attractive force decreases
in magnitude at longer distances. Critically, this model does not include explicit
alignment with neighbors, as Katz et. al. found no evidence that fish match their
orientation to their neighbors. In addition to the social rule, agents are also coupled by speed to their environment. This coupling is motivated by the results of
[Berdahl et al., 2013], where the authors found that schooling fish respond strongly
4
to the level of light they sense at their position. We use this model of collective
motion to explore the ramifications of passive group size regulation when agents are
coupled to their environment, which contains sparse, dynamic resources. We find that
with just these social rules and the simple environmental coupling, agents match their
group size to the scale of the environmental resources, and groups dynamically fuse
and break apart in order to balance exploration with exploitation of the environment.
5
Chapter 2
Revealing the hidden networks of
communication in schooling fish
Chapters 2 and 3 present two sections of a project on which I was a first author,
along with Colin R. Twomey, and secondary authors Andrew T. Hartnett, Hai Shan
Wu, and corresponding author Iain D.Couzin. This work is currently in press in Proceedings of the National Academy of Sciences, USA [Rosenthal et al., 2015]. These
chapters contain lightly edited versions of the main text and supplementary information found in the paper. In this chapter I describe the experimental technique we
used to map sensory input to behavioral output, and thus to estimate the pairwise
interactions used by fish in propagating waves of alarm. In Chapter 3, I explain how
we use the interaction network we estimated in Chapter 2 to subsequently predict
the magnitude of these behavioral cascades, and to identify socially influential and
susceptible individuals in the group. I also describe a simulation model we built to
further investigate our hypotheses about the nature of social contagion.
SBR, CRT and IDC designed the concept for the study. CRT developed the
machine vision software, HSW developed the tracking software, ATH designed and
conducted the mechanically triggered experiments. SBR and CRT designed and per6
formed numerical and statistical analyses, SBR and CRT designed and executed simulation study. SBR, CRT and IDC wrote the paper.
2.1
Introduction
The social transmission of behavioral change is central to collective animal behavior. In many mobile animal groups, such as schooling fish and flocking birds,
social contagion is often fast, resulting in dramatic waves of response, as evident, for example, when individuals are under threat of attack from predators
[Procaccini et al., 2011]. Despite the ubiquity, and importance, of behavioral contagion, and the fact that survival depends on how individual interactions scale to
collective properties [Handegard et al., 2012], we still know very little about the
sensory basis and mechanism of such coordinated collective response.
In the early 20th century Edmund Selous proposed that waves of turning of large
flocks of birds result from a direct transference of thoughts among individuals: “They
must think collectively, all at the same time, or at least in streaks or patches a
square yard or so of an idea, a flash out of so many brain” [Selous, 1931]. By the mid
1950s, however, attention had turned from telepathy to synchrony arising from the
rapid transmission of local behavioral response to neighbors, with some of the first
experimental studies of cascading behavioral change undertaken by Dimitrii Radakov
(summarized in [Radakov, 1973]). Projecting cine-film of large fish schools responding to an aversive stimulus (including moving objects within, or above, the water)
Radakov hand-traced the paths of each fish, frame-by-frame, revealing that the speed
of the “wave of agitation” propagated much faster (12-15 times faster) than the
maximum swim speed of individuals [Radakov, 1973]. Using similar methodology
Foster and Treherne [Treherne and Foster, 1981] studied rapid waves of escape response in marine skaters, and they termed this “the Trafalgar effect” in reference to
7
the speed of communication, via signaling flags, among ships in the English Navy’s
fleet at the battle of Trafalgar in 1805 that allowed information regarding the enemy to propagate much faster than the ships themselves (and past the horizon).
Since these studies, similar behavioral cascades have been found in many other
organisms [Couzin and Krause, 2003, Gerlotto et al., 2006, Kastberger et al., 2008,
Krause and Ruxton, 2002].
While describing general macroscopic properties, such as the speed, or direction,
of behavioral waves, is relatively straightforward, revealing the mechanism through
which such information propagates among individuals has proven much more difficult.
In many situations, such as when a predator attacks a group [Procaccini et al., 2011],
or when artificial stimuli are employed (e.g.
objects moving within, above
[Radakov, 1973] or dropped into the water [Marras and Domenici, 2013]), it is not
possible to differentiate between the propagation of behavior via social contagion and
that resulting from direct response to the stimulus, or some combination of both. For
example, the sound of an object dropped into the water [Marras and Domenici, 2013]
creates a near-instantaneous acoustic cue available to all individuals (the speed of
sound in water is approximately 1.5 km/s). This is further exacerbated by the
fact that response latency associated with direct behavioral response increases with
distance to the stimulus [Domenici and Batty, 1997]; thus the null expectation for
asocial response by members of a group to a stimulus would be a fast wave of
response (appearing to travel via contagion) from the stimulus outwards.
In previous studies, therefore, it has not been possible to isolate the social component of rapid collective evasion.
While simulations can qualitatively
reproduce phenomena reminiscent of such waves [Vabø and Nøttestad, 1997,
Wood and Ackland, 2007], the underlying assumptions made may be incorrect.
For example a predominant paradigm has been to consider individuals
as “self-propelled particle”, which (inspired by collective processes in physical
8
systems) interact with neighbors through social forces [Cucker and Smale, 2007,
Mogilner et al., 2003, Couzin et al., 2002, Vicsek et al., 1995]. In such models it is
usually assumed that they do so with neighbors within a fixed distance (a “metric
range” [Couzin et al., 2002, Vicsek et al., 1995]) or with a fixed number of nearneighbors regardless of their distance (“topological range” [Camperi et al., 2012]),
assumptions that, although mathematically convenient to researchers, do not
necessarily represent what is convenient, or appropriate, for neural sensing and
decision-making. Furthermore, it has been shown that metric and topological representations reflect poorly the sensory information employed during social response
in schooling golden shiner fish [Strandburg-Peshkin et al., 2013], the focal species of
the present work.
A major challenge in the study of collective animal behavior is that the pathways of
communication are not directly observable. In the study of isolated organisms, it has
long been realized that mapping the physical and functional connectivity of neural networks is essential to developing a quantitative, and predictive, science of how individual behavior is generated. By contrast, in the study of mobile animal groups the analogous issue of determining the structure of the sensory networks by which interactions,
and the resulting group behavior, are mediated remains to be explored. This is despite the fact that the structure and heterogeneity of networks is known to profoundly
impact contagious processes, from spreading neural electrical activity [Chialvo, 2010],
innovations [Young, 2011] or disease [Moore and Newman, 2000, Volz et al., 2011], to
propagating internet [Motter and Lai, 2002] or power grid [Kinney et al., 2005] failure. In all such scenarios predicting the magnitude of contagion, and identifying
influential nodes (either in terms of their capacity to instigate or inhibit widespread
contagion) is crucial.
Here we focus on studying rapid waves of behavioral change in the context of collective evasion, employing strongly schooling fish (golden shiners) as an experimental
9
study system. To uncover the process by which this behavioral change spreads, we
exploit the fact that shiners, like other fish, exhibit fast-start behavior when they perceive an aversive stimulus (e.g. via the visual, acoustic or mechanosensory system)
[Eaton, 1984], and occasionally do so in the absence of any external stimulus. Studying fast-start evasion resulting from spontaneous startle events, instead of presenting
a stimulus visible to multiple individuals, offers us the opportunity to unambiguously
identify the initiator of escape waves and to avoid confounding social and asocial
factors.
Since fast-start is mediated by a reflex circuit involving a pair of giant reticulospinal neurons, the Mauthner cells [Eaton et al., 2001], it may be expected that
individual fish would be unable to establish the causal factor for escape in others, i.e.
whether it resulted from a real threat or not. We first test this hypothesis by comparing evasion response resulting from spontaneous startles to those resulting from an
experimentally controlled alarming mechanosensory stimulus. We find no difference
in either individual response, or the propagation of waves of evasion, suggesting that
when responding to fast-start behavior shiners do not differentiate between threatinduced and spontaneous startles. This is consistent with previous experiments on
birds [Lima, 1995], and also with theoretical predictions which suggest that the risk
of predation makes it simply too costly for vulnerable organisms, like golden shiners,
to wait to determine if the escape motion of others is associated with real danger
[Giraldeau et al., 2002].
2.2
Results and Discussion
We studied free-swimming groups of juvenile golden shiners (Notemigonus crysoleucas), a widespread species of freshwater fish that swim close to the water surface
[Whittier et al., 2000] and whose schooling behavior is well characterized by vision
10
(Strandburg-Peshkin et al., 2013); in this species the lateral line is thought to contribute minimally to schooling [Burgess and Shaw, 1981, Hanke and Lauder, 2006].
Like many juvenile fish [Domenici and Blake, 1997] shiners experience very high mortality due to predation (abundance is limited by predation rather than food availability [Johannes et al., 1989]). Under our experimental conditions, that minimize
external visual and acoustic cues (Section 2.3.A), we find that shiners occasionally
exhibit spontaneous fast-starts in the absence of any apparent cue. Spontaneous
fast-starts may arise from spontaneous discharge, or so-called “synaptic noise”, in
the reflex escape circuit [Hatta and Korn, 1999], or from the reflex system being triggered somehow in only a single individual, such as via local conspecific motion, or
perhaps ripples on the surface of the water. Such startles are rare, occurring on
average approximately once every 3.3 hours, per individual (Section 2.3.B) and are
easily identifiable due to a sudden acceleration well outside the norm associated with
general schooling behavior (Fig. 2.1A, Fig. S2.1; Section 2.3.B).
We found no evidence that golden shiners differentiate between fast-starts arising
from an induced alarming stimulus (sudden mechanosensory stimulation, Fig. S2.3)
and those that arise spontaneously (Fig. 2.2). Following initiation, evasion behavior
can propagate quickly (approximately twice the maximum swim speed of any individual, Fig. 2.1D, Fig. S2.1B; Section 2.3.B) within groups, away from the point of
initiation (and thus rapidly away from any potential threat).
Waves of evasion can spread extensively through groups (Fig. 2.1A), or may
rapidly die out, resulting in a broad distribution of cascade magnitudes (number of individuals that respond), a property shared with other spreading processes, e.g. neural
activity [Beggs and Plenz, 2003], and human communication [Leskovec et al., 2009],
and with propagation speeds that decelerate over time (Fig. 2.1D, Fig. 2.2A). Consequently, predicting whether a startle event will result in a widespread behavioral
cascade is challenging. Furthermore, we find that variability in the behavioral re11
t = 0 ms
t = 167 ms
t = 833 ms
t = 333 ms
startle onset
A
t = 0 ms
B
C
2
number of events
10
1
10
0
10 0
10
1
10
cascade size
2
10
t = 833 ms
D
60
40
20
0
0
0.5
1
time after initiation (s)
speed of propagation (cm/s)
t = 667 ms
cascade spatial extent (cm)
t = 500 ms
200
150
100
50
0
0
0.5
1
time after inititiation (s)
Figure 2.1: Collective evasion. (A) Six sequential frames 167 ms apart of an
escape response cascade in schooling fish, school size of 154 fish; color indicates the
time of startle initiation with lines showing the trajectory of startled individuals.
(B) Distribution of observed cascade sizes on a log-log scale. (C) Spatial extent of
behavioral cascades: the average distance between initiator position and position of
individuals that responded time t after initiation. (D) Average speed of cascade wave
propagation, calculated as the change in the spatial extent of the cascade per second,
averaged over all cascade events. Data points in (C) and (D) show the mean, with
error bars for ± the standard error.
sponse of the initiator does not influence the nature of contagion (the influence of
the initiator cannot, for example, be attributed to its maximum speed, turning rate,
or turning acceleration, see Section 2.3.B), emphasizing the need to consider explicitly the sensory information available to potential responders, who may observe the
initiator differently due to their different relative positions throughout the group.
12
150
100
50
0
0.2
0.4
0.6
0.8
1
time3(s)
C
80
1
successful3trials
control3trials
70
standardized3speed
speed3of3propagation3(cm/s)
triggered3startles
spontaneous3startles
200
speed3of3fastest3fish3(cm/s)
B
A
60
50
40
30
20
10
triggered3startles
spontaneous3startles
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0
0.1
0.2
0.3
0.4
0.5
0.6
time3after3triggering3(s)
0.7
0.1
−0.2
0
0.2
0.4
0.6
time3relative3to3t03(s)
Figure 2.2: Comparison of triggered and spontaneous behavioral cascades (A) Average speed of cascade wave propagation, calculated as the change in the spatial extent
of the cascade per second, averaged over all triggered cascade events. Compare to
blue curve, for speed of spontaneous behavioral cascade. (B) Comparison of response
speeds of successful trials (red) to control trials (black). Speed profiles of fastest
fish in school following alarm triggering are averaged over all events, and plotted for
both successful and control trials. No significant increase in speed is observed for
control trials, and speed in successful trials was significantly higher than in control
trials (Wilcoxon rank-sum test, p < 0.00001). (C) Average speed profile for triggered
startle response (centered on zero so that the max speed is one, and occurs at t = 0.
Curve is created by averaging over all initiators. Compare to blue curve, for spontaneous startles. Data points in (A), (B), and (C) show the mean, with error bars for
± the standard error.
2.2.1
Functional mapping between sensory input and behavioral change
In order to reveal the hidden pathways of interaction, we must first establish a
functional mapping between sensory input and the relevant motor output (escape response). To obtain this mapping, we considered only the first responder
to the initiator, for which we have the clearest causal relationship. As expected
[Strandburg-Peshkin et al., 2013], visual contact (Fig. 2.4A) is a strong predictor of
response, with fish being unlikely to respond if the initiator is occluded from view,
regardless of their proximity (likelihood ratio test, p < 0.0001; Section 2.3.D).
13
To determine which properties of visible neighbors are used in decision-making
by the first responder, we considered both absolute and comparative features of the
neighborhood [Vlaev et al., 2011]. Absolute features include Euclidean (metric) distance to visible neighbors, and their angular position and area subtended on the retina
(Fig. 2.4B, for a complete list see Table 2.3). These properties depend only on the
relationship between the observer and the initiator, but are invariant to other visible
neighbors. Comparative features, on the other hand, do contain implicit information regarding multiple visual neighbors, and include properties such as topological
distance, ranked angular area, and relative metric distance (Fig. 2.4C, Table 2.3).
Log-transformed predictors were also included in the analysis, as many sensory systems operate on a logarithmic scale in order to facilitate scale invariance and to
increase dynamical range, and thus response sensitivity, reflected by the well-known
Weber-Fechner law [Dehaene, 2003].
Two approaches were employed to find the combination of features which best
predicts the identities of the first responders in the data. First we used multi-model
inference based on exhaustive (i.e. not step-wise) evaluation of all feature subsets,
where we compute a logistic regression (generalized linear model, or GLM, with a
logistic link function) for every possible combination of features and compare their
performance (Fig. 2.4D) [Johnson and Omland, 2004]. Secondly, we used an L1 regularized logistic regression [Park and Hastie, 2007] (see Section 2.3.G and SI Fig. S2.9)
of all features.
We selected the two most predictive features of behavioral response for inclusion
in the model, where we exclude other features because of low feature relevance, or
because they contain similar information as features with higher relevance (see Section 2.3.F for details on feature selection). The two model selection methods found
consistent results, with small differences due to the non-exhaustive search used by
the L1 regularization method. The two most predictive features are: (1) the loga14
Absolute features
Metric distance (,
log())
Euclidean distance from focal fish to neighboring
fish
()/
Speed ()
Angular position ()
Relative angular position of neighboring fish from
direction of heading of focal fish
/
Angular speed ()
Angular area (, log()) Angular area subtended on the retina of the focal
fish by neighboring fish
Loom (, log())
Heading ()
/
Difference between heading of focal fish and
heading of neighboring fish
/
Heading change ()
Comparative features
Topological distance
(, log())
Ranked Euclidean distance from focal fish to
neighbor fish
Relative metric distance
(, log())
Euclidean distance to from focal fish to nearest
neighbor divided by Euclidean distance to neighbor
fish.
Ranked angular area
(, log())
Rank of neighbor sorted by size of subtended
angular area on retina of focal fish
Relative angular area
(, log())
Angular area of neighbor fish on retina of focal fish
divided by maximum angular area observed by focal
fish
Table 2.3: Description of all features used in model selection.
rithm of the metric distance to visible neighbors (Fig. 2.4E), and (2) the ranking of
the initiator, with respect to others, in terms of the angular area subtended on the
retina (Fig. 2.4F) (Fig. S2.6, Section 2.3.F). This can be formulated as a probability
of response, P (i|j), in fish i given that fish j has startled.
P (i|j) = (1 + exp(−β1 − β2 LMD − β3 AR))−1
(2.1)
Where β1 , β2 , and β3 are the model coefficients fit to the data (0.30235, −1.4213,
and −0.12608, respectively), LMD is the log of the metric distance between fish i and
15
fish j, and AR is the ranked angular area of fish j subtended on the eye of fish i. The
logarithm of the metric distance is a scale-dependent absolute feature; the initiator of
a startle behavior is much more likely to trigger evasion in an observer if it is visible,
and in close proximity (Fig. 2.4E). Ranked angular area is a scale-independent and
comparative feature, ranking an individual’s neighbors by their visual size (the area
they subtend on the retina of the focal fish), such that the neighbor with the largest
angular area observed by a given individual will have an AR = 1, the neighbor with
the second largest angular area will have an AR = 2, and so on (Fig. 2.4F). Ranking
in this way encodes implicit information about other neighbors, and thus facilitates
inhibition. For example if an initiator has a relatively low rank - that is, there are
individuals who subtend higher areas on the retina who do not respond - then this
inhibits the focal individual from responding. This does not imply counting, and
is consistent with the fast computation of preattentive features of relative stimulus
conspicuity, or salience [Koch and Ullman, 1987].
16
Figure 2.4: The sensory basis of social contagion. (A) Neighbors (blue) in
direct visual contact with the initiating individual (red), just prior to the initiator’s startle. School of 151 fish. Direct lines of sight from the eyes of neighbors
to the body of the initiator are shown in white, with overlapping fields of view
drawn additively. (B) Illustration of the angular position relative to heading (θ)
and angular area (AA) estimation for a neighboring fish. (C) Metric distance, topological distance, and ranked angular area for two neighboring fish. (D) Relative
feature importance [Johnson and Omland, 2004], for the top five predictive features
plus three others previously considered important for information transfer in animal
groups [Dill et al., 1997]. (E) Probability of an individual’s response as a function of
the top predictor, log(metric distance), holding ranked angular area constant. (F)
Probability of an individuals response as a function of the second predictor, ranked
angular area, holding log(metric distance) constant. In (E) and (F), the solid blue
line is the fit of the model with the top two feature variables to the first responder
data, while the shaded blue area represents 95% confidence intervals. Histograms
above are distributions of first-responders in the data, while histograms below are
distributions of non-responders.
17
2.2.2
Reconstruction and analysis of the quantitative interaction network
This mapping between sensory input and behavioral response (parameterized only
by the first responders to each event) allows us to subsequently construct a network
predictive of how behavior will spread through the group, where the weight of an edge
between individual i and j, wij = P (si |sj ), is the probability of a behavioral response
by individual i if individual j exhibits behavioral change.
Since neither visibility nor the ranked angular area are necessarily symmetric
properties (wij need not equal wji ), social connectivity is both directed and weighted
(Fig. 2.5). A visualization of this network just prior to an initiation event is shown
in Fig. 2.5, where a magnified view of a subset of the group highlights the weighted,
directed interconnectivity. Due to individuals’ spatial embedding, the role of distance, and the spatial regulation exhibited with respect to near-neighbors, networks
exhibit a relatively lattice-like structure with few strong, but many weak, long-range
connections (Fig. S2.8). In Chapter 2 we use this reconstructed interaction network
to predict the magnitudes of behavioral cascades.
18
Figure 2.5: Interaction networks are complex, with weighted, directed pathways of communication. (A) An interaction network, where link weights are
determined as a function of the log metric distance and ranked angular area. (B)
Directedness of the network. While many links are reciprocal (green), there also exist
connections which are much stronger in one direction than the other (yellow/blue).
2.3
2.3.A
Technical details
Experimental methods
Groups of 150 ± 4 juvenile golden shiners (Notemigonus crysoleucas) were allowed to
school freely in a 2.1 × 1.2 m experimental tank, in 4.5 – 5 cm deep water. Measures
were taken to ensure the experimental arena was acoustically and visually isolated
from external stimuli. Two layers of sound insulation were placed under the tank
to provide acoustic isolation, and the tank was enclosed in a tent made of featureless white sheets, within an environmental chamber, within which no persons were
present during filming. We chose to work with this species because they display
highly coordinated schooling behavior [Katz et al., 2011, Tunstrøm et al., 2013]. Juvenile fish approximately 5 cm in length were purchased from I. F. Anderson Farms
(http://www.andersonminnows.com), and kept in a climate controlled laboratory
19
space for two months before they were used in experiments. There were approximately 1,000 fish in total, housed in seven separate 20 gallon tanks, at a density of
approximately 150 fish per tank. The tank water was conditioned, de-chlorinated,
oxygenated, and filtered continuously. Fifty percent of tank water was exchanged
twice per week. Nitrates, nitrites, pH, saline and ammonia levels were tested weekly.
The room temperature was controlled at 16 ◦ C, with 12 hours of light and 12 hours of
dark, using dawn-dusk simulating lights. The fish were filmed overhead at 30 frames
per second, using a Sony EX1, at a resolution of 1920 × 1080 pixels. All fish were fed
flake food three times per day, and were filmed in the experimental tank 2–4 hours
after feeding. The 138 spontaneous startle events were collected from 5 different free
swimming groups of 150 (± 4) fish, with each group composed of the members of a
randomly selected tank. Each group was recorded for 53 minutes. The experimental
procedure is identical to that of [Katz et al., 2011], where more details can be found.
All experimental procedures were approved by the Princeton University Institutional
Animal Care and Use Committee.
2.3.B
Stereotypy of startle response in golden shiners
Like many fish, golden shiners exhibit a startle escape response to the presence of
certain stimuli [Domenici and Blake, 1997]. There are several types of startle response
and “fast-start” motion (e.g. C and S starts) exhibited by fish, but for this work we
focused on a simple characterization of response given by the sudden acceleration of
an individual, deviating well outside the norm of normal swimming behavior. Faststarts are known to be induced through vision or sound alone, or by multi-modal
combinations of unexpected stimuli [Eaton, 1984]. We identified 138 instances of
spontaneous startle events, in which a single individual within the group performs an
initiating fast-start in the absence of any known cue. Such spontaneous startles were
rare, occuring, per fish, approximately every 3.3 hours. This rate will likely depend on
20
how secure the fish feel in the environmental arena. We tracked a window of 60 frames
(2 seconds) around the startle event itself, where the tracking window was determined
by the length of time between initiator and last responder. These fish are capable of
very rapid accelerations, and sustain speeds of about 50 cm/s (Fig.S2.1B) for a short
time after the onset of a startle. We observe in our data, in both spontaneous and
triggered alarms, that the cascade propagates outwards at about 100 cm/s (Fig. 2.5C,
Fig. 2.2A). Given that our schools, to a rough approximation, measure 100 cm along
the long axis when occupying the polarized state, it would only take about 1 second
for the alarm to reach every individual, were the alarm sustained. Since the false
alarms we study never propagate to the entire school, 1.6 seconds is long enough to
observe the entire cascade.
In 71 spontaneous startle events, one or more individuals are classified as responding to the initiator. We most often observed cascade sizes of 1 or 2 individuals, with
the alarm only rarely spreading to a large fraction of the school (Fig. 2.5B). We classified startles and responses using a speed and turning rate threshold of two standard
deviations above the baseline swimming speed and turning rate (this baseline varies
from event to event). Initiators and first responders have speeds and turning rates
higher than the school baseline (Wilcoxon rank sum test, initiators p < 0.00001; first
responders p < 0.00001). The later responders tend to be slower, however, making it
necessary to compare to the school’s average swimming speed to identify those fish
which are responding.
In many cases individuals have been found to be unable to distinguish, from the
behavior of conspecifics, the reason for sudden behavioral change (i.e. whether a
threat is present or absent [Lima, 1995]). This is especially true for organisms like
fish whose startle behavior is relatively stereotyped [Domenici and Blake, 1997]. Variability in the magnitude of the initial startle is not predictive of resulting cascade size,
supporting the hypothesis that neighboring fish process an observed startle event as a
21
binary behavior, and cannot distinguish between real alarms and false alarms, which
we test explicitly below. Qualities such as the maximum speed (likelihood ratio test,
χ2 (df = 1, N = 138) = 2.04, p = 0.15), maximum acceleration (likelihood ratio test,
χ2 (df = 1, N = 138) = 2.89, p = 0.09), maximum turning rate (likelihood ratio test,
χ2 (df = 1, N = 138) = 1.61, p = 0.21), and maximum turning acceleration (likelihood ratio test, χ2 (df = 1, N = 138) = 2.67, p = 0.11) are not good predictors
of the influence of the initiator, where the magnitude quantities were fitted to the
cascade sizes using Poisson regression because the dependent variable (cascade size)
is a count variable. Response latency of first responder increases with distance from
initiator, consistent with [Domenici and Batty, 1997], Fig. S2.1A. The difference between the swimming speeds of startlers, responders, and non-responders can be seen
in Fig. S2.1B. The change in speed over time characteristic of a startle is shown in
Fig. S2.1C, where we plot the mean normalized speed and standard error of all observed initial startle events, relative to time of maximum speed. The startles are
characterized by a rapid acceleration followed by a slow deceleration.
B
C
0.07
0.2
0.06
frequency
0.25
0.15
0.1
0.05
0.05
0.04
0.03
0.02
0
20
10
0
initiators
first−responders
non−responders
0 10 20
count
0
5
10
15
20
25
distance/Ccm)/from/s0/to/s1/at/t0
0.01
0
0
20
40
60
80
normalized/speed/±/SE
0.08
count
latency/Cs)/from/s0/to/s1
A
100
1
0.8
0.6
0.4
0.2
maximum/sustained/speed/Ccm/s)
0
−0.2
0
0.2
0.4
time/relative/to/t0/Cs)
0.6
Figure S2.1: (A) The latency between the initial startle and first response grows with
distance between the initiator and first responder. (B) The distribution of speeds of
initiators (blue), first responders (red), and non-responders (yellow). Data are shown
in bars, while the kernel-smoothed distribution is overlaid in dotted lines. Initiator
and first responder speeds are significantly higher than non-responder speeds. (C)
Average shape of speed-profile for startle response (centered on zero and standardized
so that max speed is one, and occurs at t = 0. Curve is created by averaging over all
initiators, with error bars showing the standard error).
22
The spatial extent of the behavioral cascade was calculated by measuring, for each
frame, the largest distance between locations of all responses that occurred in that
frame and the location of the cascade initiation. We averaged over all behavioral
cascades in the data, and plotted the mean ± the standard error (Fig. 2.5C). The
speed of propagation was calculated by measuring the difference in spatial extent per
frame, for each behavioral cascade in the data. Again, we averaged over all cascade
events and plotted the mean ± the standard error (Fig. 2.5D).
Golden shiners have been found to swim in both ordered and disordered states.
Ordered states are characterized by high polarization (all fish swimming in a common direction), or high rotation (where fish are locally polarized, but are rotating
around a common center), while disordered swimming is characterized by both low
rotation and low polarization [Tunstrøm et al., 2013]. Groups of 150 fish only very
rarely enter the disordered regime, predominantly exhibiting polarized or rotating
motion [Tunstrøm et al., 2013]. We find that initiators occur nonrandomly in the
group structure. For example, we find that initiators tend to be found closer to the
group boundary than non-responders (Wilcoxon rank sum test, p < 0.0001 all group
states), and in polarized states, initiators are more likely to be found closer to the
group front than non-responding fish (Wilcoxon rank-sum test, p = .008, polarized
groups). However, locations of initiators don’t tend to happen in more or less dense
regions than non-responding fish (initiators: Wilcoxon rank sum test, p = .09, all
group states). See Fig. S2.2 for details. We found evidence that cascades propagate
slightly further in polarized groups than in rotating groups, a weak, but statistically
significant increasing trend (R = 0.23, p = 0.02).
23
A
B
0.25
30
all5fish,5N5=520,501
s1,5N5=571
s0,5N5=5138
25
s0,5N5=5138
0.2
0.15
15
10
0.05
5
0
5
10
15
20
25
distance5from5boundary5Ncm=
s0,5N5=534
0
s1,5N5=522
0.025
20
0.1
0
0.03
all5fish,5N5=520,501
s1,5N5=571
frequency
0.3
frequency
C
35
frequency
0.35
all5fish,5N5=55038
0.02
0.015
0.01
0.005
0
0.02
0.04
0.06
0.08
0.1
local5density5Nfish/cm2=
0
−50
0
50
distance5front-back5Ncm=
polarized5schools
Figure S2.2: Distributions of initiator, first responder, and non-responder positions
in relation to (A) distance from group boundary, (B) local density, and (C) distance
front-back (in polarized schools).
Response to induced fast-start behavior
Considering the inherently stereotyped nature of fast-start behavior, and the lack
of sensitivity of responders to the magnitude of the initiating startle (see above), we
expect that spontaneous fast-starts are responded to in the same way as induced faststarts (i.e. those that are triggered in direct response to an aversive stimulus). To test
directly this hypothesis, and to validate our approach, it is necessary to compare social
contagion resulting from spontaneous and triggered evasion maneuvers. However,
to avoid confounding asocial and social factors in the response, we must employ a
stimulus that triggers fast-start in focal individuals but is not perceived by potential
responders.
Initially we attempted to employ a laser mounted on a series of servo motors
positioned 2 meters above the tank. Custom-designed software allowed us to aim
quickly at any location in the tank and then to turn on and off the laser, providing
a rapid and relatively localized ‘blink’ of light. While this approach was successful
at inducing startles, and could be easily controlled, we found that it was not a truly
local stimulus since light was scattered off both the water surface and the base of the
24
tank. Furthermore we could not exclude the possibility that the necessary motion of
the laser mounting may have been perceptible to the fish below.
Following this attempt we focused on developing a device that could provide local tactile (mechanosensory) stimulation to induce fast-start evasion behavior. We
designed the monofilament perturbative stimulus illustrated in Fig. S2.3 that allows
us to automatically raise a thin filament (made out of polyethylene) across a channel of dimensions 122 cm × 30 cm within our experimental arena, using a HiTec
HS-5965MG high-speed programmable digital servo to control motion. Following
computer-controlled initiation by an experimenter, the filament was raised from the
base of the tank to a height of 3 cm within the water body, before lowering to the base
again, with the full procedure taking less than 200ms. To remove potential acoustic
stimuli caused by this motion, the design allowed the servo motors controlling the
motion to be placed outside the tank (Fig. S2.3A) and teflon was applied to the areas
where the filament was in contact with the tank structure to minimize friction, and
thus minimize sound production. A hydrophone was placed in the tank allowing us
to validate quiet operation and generally effective sound insulation.
The striking of fish by the monofilament induced fast-starts, whereas releasing the
mechanism in close proximity to fish without physical contact resulted in no response
(Fig. 2.2B), thus demonstrating that perturbations were effectively localized. We
triggered startle responses in schools of 150 fish, recording the motion of individuals
at 150 frames per second, using a Prosilica GX-1050 (with a 9mm F1.4 c-mount lens,
and resolution of 1024 × 768). The fish were allowed to habituate for one hour prior
to triggering startles. Following the habituation time, startles were triggered at each
of the three monofilament locations (see Fig. S2.3B), over the next three hours, with
a minimum of five minutes between triggering events. In total, we collected 67 triggered startle events. Shiners maintain a density of 0.02 fish/cm2 in the experimental
channel, matching what we observe in free swimming schools (Fig. S2.2).
25
A
SupportyBeam
6topyview2
Servo
6sideycut-away2
TankyWall
ChannelyWall
Supportyw/
TeflonyGuides
Elastic
B
TeflonyDonut
TankyFloor
122cm
18cm
30cm
43cm
43cm
61cm
213cm
Figure S2.3: Experimental setup for triggered startles (A) The monofilament
is raised by a servo, pictured from above in top view, which is fixed outside the
light tent. Here we used a HiTec HS-5965MG high-speed programable digital servo,
with a 10cm arm. The monofilament then runs through the white support structure,
through the a white teflon return donut, through the channel wall, and then through
the same structures in reverse order on the other side—shown in side cut-away. The
white teflon donut is attached to another filament and an elastic. When the servo is
triggered the wire raises, stretching the elastics, which in turn pull the monofilament
back to the channel floor when the servo relaxes. (B) All channel experiments were
conducted in the above experimental tank. The fish were constrained a region of the
213cm × 122cm tank consisting of two 43cm × 43cm gravel-bottomed shelters and a
122cm × 30cm channel connecting the shelters. The channel had clear internal walls
with clear monofilaments placed along the tank at three different positions (dashed
blue lines). These monofilaments could be quickly raised and lowered to strike the
underside of passing fish.
26
The trajectories of all fish were obtained via automated tracking (see Appendix B,
below). The timescale of individual response to triggered fast-start evasion (Fig. 2.2C,
for response to spontaneous and induced fast-startles) and the speed of resulting
waves of evasion (Fig. 2.2A, for response to spontaneous and induced fast-starts)
were indistinguishable, providing evidence that our fish do not differentiate between
spontaneous fast-starts and those triggered by exposure to an aversive stimulus. We
performed linear regressions on the speed of propagation for both spontaneous and
triggered startle events, compared the regression coefficients (slope and intercept)
with t-tests, and found that the slopes and intercepts for both types of startles do
not differ, with p-values of p = 0.97, p = 0.49, respectively.
Since the triggered startles were conducted at a frame rate of 150 fps, compared
to a frame rate of only 30fps for the spontaneous startle data, we smoothed the
triggered startles using a 5-frame smoothing window, and plotted both speed profiles
in Fig. 2.2C. Smoothing was necessary in order to be able to compare the maximum
fish swimming speeds between datasets.
2.3.C
Description of tracking software
Automated tracking of large groups of fish is a challenging problem due to their
non-rigid body structure, as well as partial and total occlusions of individuals when
viewed from above, generated by fish swimming over and under one another. Occlusions occur in our experiments approximately twice a minute per individual for
free-swimming groups of 150 fish, which is enough to be problematic for accurate
tracking. The SchoolTracker computer vision software used in this project (developed
by HSW) was created to reliably solve this problem. SchoolTracker is comprised of
three modules: detection, tracking, and track linking. The detection module takes
video recorded from a high-resolution camera mounted above the experimental arena
(providing an overhead view of the school) and identifies possible fish body positions
27
and orientations. This step is made difficult by the presence of occlusions, which is
why most image segmentation methods fail. Instead, the detection module looks for
local features such as corners and line segments, and systematically estimates first
fish head positions from the corners, then body position and orientations from the
optimal assignment between line segments and estimated head positions. It then determines the statistical threshold values based on the distributions of head intensity
or fish image intensity difference and removes the detection outliers based on this
threshold.
Given a set of detections for each frame, the tracking module identifies which
detections from one frame to the next belong to the same individual, forming a track
for that individual over time. The multiple object tracking algorithm we use is able
to track the complex motion of each fish, and resolve occlusions (overlapping bodies)
by matching against a library of learned robust body images. The algorithm first
estimates the motion states of each individual via two types of adaptive alpha-beta
filters (a form of simplified Kalman filter [Blackman and House, 1999], but modified
for tracking fish in this instance), and then associates the predicted future states of
the individuals with detections in the next frame by minimizing the sum of association
costs, which in our implementation was defined by the difference of head positions,
body angle and corresponding learned body images. This is formulated as a linear
assignment problem, and is solved by the data association algorithm that is similar to
what we proposed previously in [Wu et al., 2011]. It scales very well when handling
a large number of fish as it was designed to deal with detection error, temporal
disappearance due to short-term occlusion or missed detections.
Due to some (rare) long-term occlusions (e.g. when two fish completely overlap
for many frames), the tracks constructed by the tracking module may be broken into
shorter segments. Since the number of fish in the experimental arena is constant
throughout the experiment, the track linking module is used to piece together track
28
segments via combinatorial optimization, to construct coherent tracks that span the
length of the full video. Upon inspection, manual relinking is performed to correct
errors. By itself, the automated system recovers over 95% of tracks without misidentifications, and over 97% of tracks without missing frames, minimizing the amount of
post-tracking correction necessary (accuracy evaluated on 10 randomly selected, 120
frame video segments of 150 free-swimming golden shiners recorded at 30 fps).
The SchoolTracker software provides a simple graphical user interface to streamline the process of parameter tuning and track verification/editing, so that the final
results accurately reflect the positions and orientations of each individual fish in the
school over the whole length of the video.
2.3.D
Description of vision software
Visual, acoustic, touch, and lateral line sensory modalities may all contribute to faststart behaviors. Due to the importance of vision to schooling in Golden shiner fish
[Strandburg-Peshkin et al., 2013] and the relatively small contribution of the lateral
line [Burgess and Shaw, 1981], we (CRT) developed a method for approximating the
planar field of view of each individual in a group when individuals are restricted
to an effectively planar environment as in this experiment (restricted depth). After tracking, the positions and orientations of each individual were used as the first
step in estimating the posture of each individual’s body. Posture estimation is trivial when fish are non-overlapping: there is a null-cline in the image Laplacian that
traces the midline of the fish’s body. Points along this null cline are sparsely sampled using a 2d-tree data structure, and linked using a greedy energy minimizing line
following method, where the energy of a line is minimized by reducing the distance
between subsequent points, and minimizing the angle of curvature between adjacent
line segments. Finally, a cubic basis spline is fit to the selected points to approximate the midline curve, while the flanks of the individual are estimated by searching
29
for boundary pixels (after thresholding to remove the background) perpendicular to
the tangent of the curve, and fitting a linear function to their distance relative to
the curve (similar to [Fontaine et al., 2008]; here we use the Thiel-Sen estimator for
robustness). This gives an approximation of the body of each individual summarized
by the coefficients of the basis spline and the slopes and intercepts of the left and
right flanks.
Posture estimation is substantially complicated by the occurrence of partial occlusions. This is why skeletonizing after thresholding out the background fails to work
well for midline estimation, as clumps of individuals yield highly irregular skeletons.
The method above is robust to this, but may still fail when individuals cross while
their bodies are curved. To counter this, information about body posture estimation
on previous frames is used to predict the posture of an individual in the future using
a simple Kalman filter.
After the body posture for every individual in the school has been characterized,
an estimation of the location of each eye, for each individual, is made using the current
heading of the fish, the position of the fish’s head, and the estimated flanks. From
the midpoint position of the head, the transition points from body to background on
either side of the individual (perpendicular to the heading) are determined. Estimates
are improved using a Kalman filter to smooth over brief occlusions and possibly noisy
detections. Line segments encoding the shape of every body are stored in a simple
spatial data structure (in this case a 2d hash table) for fast look-up based on spatial
coordinates. For each eye of each individual, rays are cast (at a resolution of π × 10−3
radians per ray) from that position and the point of first intersection with other
geometry in the frame is recorded (distance and identity). The total number of rays
from an eye intersecting each individual in the scene is recorded as the angular area
subtended by that individual on the focal individual’s eye. The circular mean of
the rays intersecting the other body is recorded as that individual’s angular position
30
relative to the focal individual. An example showing line of sight and the angular
area occupied by one focal individual in the field of view of each member of a group
is shown in 2.4A.
2.3.E
Features used in model selection
Preliminary analysis allowed us to eliminate the Voronoi neighborhood as a predictive
feature in our regression models. The Voronoi tessellation of the group was determined by computing the dual of the Delauney triangulation of individuals’ positions
[Watson, 1981]. The Voronoi neighborhood alone is a good predictor of first responders, as it contains similar information to other features such as metric and topological
distance. However, when it is included in models with these other feature variables,
the Voronoi neighborhood loses significance (logistic regression, likelihood ratio test,
p = 0.15). This is reasonable, as the Voronoi neighborhood is a binary variable, and
doesn’t contain as much information as metric distance or angular area.
In our main model selection procedure, we included 12 feature variables, motivated
by the literature [Dill et al., 1997, Partridge et al., 1980], and listed in Table 2.3. The
features each belong to one of two categories: absolute or comparative. Absolute features are those that contain no implicit information about what other neighbors are
or are not doing. Examples of absolute features are metric distance, relative position, relative heading, and subtended angular area. Comparative features do contain
implicit information about other neighbors. Some examples of comparative features
are ordinal variables, such as topological distance and ranked angular area, and relative features, like relative metric distance (MDmin /MD(i)) and relative angular area
(A(i)/Amax ).
Although there is some collinearity in the predictor variables, as is to be expected
since some predictors such as metric and topological distance are good proxies for
each other, the correlations are not strong enough to merit any predictor’s exclusion
31
from the model selection. We inspect the variance inflation factors (VIF) between
pairs of predictor variables.
V IF =
1
1 − R2
(S2.2)
A good rule of thumb [Chatterjee and Hadi, 2013] is that VIFs less than about 10
are acceptable. Metric distance and topological distance are the most collinear, with
a VIF = 2.55, so we can justify including all pairs of predictors without being overly
redundant.
L
S
dTH
TH
dH
H
RelMD
RelA
AR
MD
TD
AA
0
un−transformed feature
log−transformed feature
50
100
150
individual feature χ2
200
250
Figure S2.4: Comparison of log-transformed features to un-transformed features.
Some features have improved individual performance when using the log-transformed
version. We use both log-transformed and raw (un-transformed) features in our model
selection procedures. We include the log-transformed version of a feature if it has significant predictive power on its own, and if using the log-transform results in an
improvement in individual feature performance.
Log-transformed versions of features were also included as sensory perception is
often on a log-scale [Dehaene, 2003]. We included log-transformed versions of features
for which taking the log improved the maximum likelihood of that individual feature
predicting the data (Fig. S2.4). Along with the log-transforms, two more transformations were needed, in order to change relative angular position (θ), and loom (L)
into forms that made sense for a logistic regression model, such that a monotonic re32
lationship could be expected between the independent and dependent variables. We
transformed θ into cos(θ − φ), where φ was the phase shift resulting in the optimal
log-likelihood for the feature (Fig. S2.5A). Similarly, we selected |L| (absolute value
of loom) instead of L, because |L| had a higher log-likelihood (Fig. S2.5B). Neither
of these features ended up performing well in the model selection procedures.
A
B
−298
log-likelihood
cos(θ-θshift)
−298.5
l ue
va
e
lut
so loom
ab
−299
−299.5
−300
ed
rm
sfo
n
tra loom
un
−300.5
−301
−301.5
0
0.5
1
1.5
2
2.5
3
3.5
θshift
0
5
10
15
χ2
20
25
30
35
Figure S2.5: Transforms of feature variables, fitted to first responder data. (A)
Cosine-transformed and shifted relative angular position. The best shifted value for
angular position is approximately π/2. (B) Comparison of un-transformed loom to
absolute value of loom. The absolute value of loom out-performs the un-transformed
loom.
2.3.F
Multi-model inference
For each feature subset, a generalized linear model (GLM) with a logistic link function is fit to predict the probability that a witness of the startle event would be likely
to be the first responder. We use a multi-model inference technique to find the subset
of the 12 features and their log-transforms used in our analysis that best fit our data.
The important features are not simply the features used in the “best” model, as this is
likely subject to overfitting given that we are testing every possible feature combination, and thus could vary depending on the exact dataset used. Instead, the relative
importance of each feature, ωf , is estimated following [Johnson and Omland, 2004]
33
as
ωf = Z −1
X
δif exp(BICmin − BICi )
(S2.3)
i
ki
BICi = l(θˆi ) − log(N ) + O(1)
2
(S2.4)
≈ log P (s1 |s0 , Mi )
(S2.5)
Where BIC is Schwartz’s Bayesian Information Criterion, Z is a normalization constant, l(θˆi ) is the maximum likelihood estimation for the parameters of model i, ki is
the number of features in model i, and δif is an indicator function that is 1 if model i
contains feature f , and zero otherwise. The BIC of each model i is compared against
the minimum BIC score obtained by any model, BICmin , such that features found in
models close to BICmin are given more weight than features in models significantly
greater than BICmin . P (s1 |s0 , Mi ) represents the probability of a fish being a first
responder (s1 ), given a startle by a neighbor fish (s0 ), using model Mi .
B
logMD
AR
logAR
logTD
logAA
dTH
H
logRelMD
RelMD
dH
RelA
AA
L
MD
logL
RelA
TD
TH
S
topT15Tmodels
allTfeaturesTusedTinTmodelTselection
A
0
0.2
0.4
0.6
0.8
relativeTfeatureTimportance
AR+logMD
logMD+logAR
logAA+logTD
AR+logTD+logMD
AR+dTH+logMD
AR+relMD+logMD
AR+logMD+logRelMD
AR+logTD
AR+H+logMD
logTD+logAR
logTD+logMD+logAR
logTD+LMD
AR+dH+LMD
AR+LMD+logAR
MD+AR+logMD
0
1
0.1
0.2
0.3
0.4
relativeTevidenceTweight
Figure S2.6: Results from multi-model inference. (A) Relative importance of every
feature used in model selection. (B) Relative evidence weight for the top 15 models.
Results from model selection on the 12 untransformed candidate features and 7 logtransformed candidate features, when ranked by feature weight, place log Euclidean
34
distance, ranked angular area, log ranked angular area, and log topological distance
as the top four candidates for inclusion in any model (in that order, Fig. S2.6). However, one important consideration is the extent to which different features exclude
each other. Although two features may have similar relative importance scores, and
thus seem to both be necessary for a good model, the presence of both in the same
model may be redundant with one essentially excluding the other. This is the case
for log metric and log topological distances (Fig. S2.7); one is sufficient for inclusion
in a model, and of the two log metric distance has a higher relevance score (the relative feature importance for metric distance is four times greater than for topological
distance, Fig. S2.6). Similarly, log ranked angular area and ranked angular area mutually exclude each other. The remaining features have relevance scores that are even
lower, more than ten times lower than the top feature, so we exclude them from the
final model. We calculate this joint conditional importance as follows:
P
Cf g = Pi
δif δig exp(BICmin − BICi )
f
i δi exp(BICmin − BICi )
(S2.6)
Taking this into account, as well as the relative ranking of best models (Fig. S2.6B),
we find that a model that includes ranked angular area and log metric distance is
supported best by the data.
35
logMD
AR
logAR
logTD
logAA
dTH
H
logRelMD
RelMD
dH
logRelA
AA
L
MD
RelA
TD
TH
S
TH
A
L
M
D
R
el
A
TD
A
lo dH
gR
el
A
H
M
R D
el
M
D
el
lo
gR
lo
lo
gM
D
AR
gA
R
lo
gT
D
lo
gA
A
dT
H
S
Figure S2.7: Feature pair conditional probabilities. This plot represents the likelihood, weighted by the relative evidence weight of the model, that the column variable
appears in a model given that the row variable is in the model. Some variables exclude
each other, meaning that they are not frequently found in the same highly weighted
models, even though they may show up separately in highly weighted models. This
would be expected for variables that contain redundant information. Topological distance and metric distance exclude each other here, meaning that they contain similar
information, and only one should be included in the final model.
We use this pairwise mapping of sensory input to behavioral response to build
the interaction networks. While the majority of strong connections in these networks
occur within a small local neighborhood, there also exist many long range weaker
connections, connecting more distant parts of the network (Fig. S2.8).
36
p(wij > τ)
A
B
1
τ = 1E-3
τ = 1E-2
τ = 1E-1
0.5
0
log wij
0
−5
−10
−15
1
50
100
150
TD + Ν(0, ε)
Figure S2.8: Sensory network edge weight (wij ) as a function of topological distance
(TD) between i and j, showing a lattice-like emphasis on spatially near-neighbors, but
with many weak, long-range (topologically far) connections. (A) shows a histogram of
the (normalized) frequency that an edge connecting i to j is greater than a threshold
τ (see legend), as a function of the topological distance between i and j. High
probability edges (black) are highly concentrated at low topological distances, while
a large number of low probability edges (blue) can still be found at much larger
distances. In (B), the (log-transformed) edge weights are plotted as a function of
topological distance (plus a small amount of Gaussian noise to prevent points from
otherwise obscuring each other at each discrete interval) to show the full distribution
of points summarized in (A). Data in both panels are from 10 networks subsampled
uniformly at random from the set of all networks (one per trial).
37
2.3.G
L1-penalized logistic regression for feature selection
An alternative to the multi-model inference above is to fit a logistic regression model
with all variables present, but with an L1 norm imposed on the regression problem
to both prevent overfitting and favor sparse models. We employ a logistic function
to model the probability that an individual will be the first to respond to an initial
startle, based on the ensemble of features quantified about the initiator.
To enforce that only non-redundant features are used, θ (the vector of coefficients
weighting the variables present in the logistic function modeling P (s1 |s0 ), where s0
represents the initiator of the startle event, and s1 is the first responder) is chosen
such that the maximum likelihood estimate of the model is balanced against the L1
norm of θ, with the strength of the penalty governed by λ. By varying λ, the tradeoff
between model complexity (the number of features with non-zero coefficients) and
the predictive power of the model can be explored, with higher values of λ yielding
less complex (fewer non-zero coefficients) models of P (s1 |s0 ). The goal is to find
max
θ
N
X
(i)
(i)
log P (s1 |s0 ) − λ||θ||1
(S2.7)
i=1
where i indexes all the candidate initial startle/first responder pairs in the data, and
|| · ||1 is the L1 norm.
Choosing a model based on λmin (the value of λ at which the binomial deviance of
the model is minimized) gives a set of features with nonzero coefficients that represent
the best balance between model performance and the number of features used by the
model. A suggested heuristic to be even more conservative about overfitting is to use
the smallest λ1SE < λmin that is within 1 standard error of λmin . Results from this L1
regularization are shown in Fig. S2.9.
38
A
B
logMD
logAR
logTD
L
AA
RelMD
H
logL
logRelMD
logRelA
logAA
S
dTH
TH
dH
RelA
AR
MD
TD
0
logMD
logAR
RelMD
logTD
AA
L
logL
logRelMD
logRelA
logAA
S
dTH
TH
dH
H
RelA
AR
MD
TD
0.2
0.4
0.6
featurecweightc
atcminimumcdeviance
0.8
0
0.2
0.4
0.6
0.8
featurecweightc
atcminimumcdeviancec+c1/2cSEc
Figure S2.9: Results from L1 regularization. (A) Feature weights for the value of
λ with minimum binomial deviance (λmin ), and (B) feature weights for λ > λmin at
1/2 standard error greater than the binomial deviance for λmin . The top two selected
features are in agreement with the results of the multi-model inference method.
Interaction terms
To ensure that we are not missing an important synergistic effect from feature pairs,
we separately consider models which include interaction terms from the top 4 performing features. Only the top individually performing features were included because using all pairwise interaction terms, let alone all interaction terms, from our
full feature set would be computationally prohibitive. For pairwise interactions, the
number of terms grows with the square of the number of features, while for all possible interactions the number grows exponentially with the number of features. In
an L1 regression, none of the interaction terms were selected as important, leading
us to conclude that interaction effects from pairs of features are not necessary in the
model.
Positioning of individuals within group
On a 10 minute video for 20 verified correct tracks, we asked what is the typical
duration of individual spatial fidelity within the group, with respect to the group
39
boundary. We find that the residency time of any particular position in the group
is relatively short (about 25 seconds) given the decay in the autocorrelation function
for distance from the boundary (Fig. S2.10A,B), and all 20 individuals of the random
subset selected spend time on both the boundary and in the interior of the group,
thus individual motion in golden shiner schools is relatively fluid.
A
B
N = 20
t = 25 s
t = 24 s
Figure S2.10: Lin-log plot of distance from boundary (A) autocorrelation function
(ACF) for 20 individuals in a group of 150, with random reshuffled ACF baseline
(dashed line). Red vertical line marks the decay time (minium lag in autocorrelation
function needed to reach a correlation indistinguishable from noise). In (B), a density
plot shows the decay time for 20 individuals, with the median decay time shown in
red.
40
Chapter 3
Predicting complex behavioral
contagion
3.1
Predicting behavioral cascades and the relationship between spatial position and social influence in groups
Several measures have become prominent predictors of influence in networks, including an individuals degree (number of connections) [Albert and Barabási, 2002] and
betweenness-centrality (the number of shortest paths that pass through a focal individual) [Barthelemy, 2004]. In the study of contagious disease, those individuals who
have a large number of social contacts (a high degree), yet whose contacts do not
form a tight clique (i.e. have a low clustering coefficient [Volz et al., 2011]) have the
capacity to spread an infection extensively (Newman, 2003). An important difference
exists between disease and social contagion, however. Whereas disease transmission
can follow contact with a single infected individual (simple contagion), in many so-
41
cial processes behavioral change depends on reinforcement through multiple contacts
(complex contagion [Centola et al., 2007]).
A central focus is whether we can predict - immediately preceding the event whether a change of behavior by any single individual will propagate extensively,
thus profoundly influencing group behavior, or will result in little, or no, collective response. As discussed previously, variability in the motion characteristics of
the initiator do not explain influence (Section 2.3.B). By contrast, however, the local structural properties of the complex interaction network are highly informative.
Unlike predictions arising from epidemiological contagion, we find that it is those
individuals with a high local weighted, directed clustering coefficient (a measure of
path redundancy / the interconnectedness of a nodes neighborhood [Fagiolo, 2007])
who have the most influence in groups (generalized linear model fit, likelihood ratio test, p < 0.0001), acting as the most effective spreaders of behavioral change
(Fig. 3.1A). As distinguished in [Centola and Macy, 2007], in simple contagion processes such as most epidemiological contagion, multiple connections between neighbors, or a high clustering coefficient, creates redundancies in the network which reduces propagation, but in complex contagion processes , these multiple ties provide
a reinforcement effect, helping the response to propagate further. In contrast to
analysis of social contagion on online social networks, such as Twitter and Facebook
[González-Bailón et al., 2011, Pei et al., 2014], individuals’ proximity to the core of
the network is not predictive of social influence (generalized linear model fit with
weighted k-core, likelihood ratio test; p = 0.18, Table S2). In addition, individuals
with a high betweenness centrality, are less likely to be influential (Section 3.4.D).
These centrality measures are not good predictors of social influence in our study,
because individuals with high centrality tend to be located near the group center,
where they are more inhibited from responding because they have a large number of
non-responding neighbors.
42
Since interaction networks are directed, there is a subtle but relevant distinction
between influence and susceptibility to influence. From an initiator’s perspective, a
higher local weighted clustering coefficient represents a more susceptible neighborhood. That is, once the alarm reaches at least one neighbor, more individuals in the
neighborhood will observe multiple alarms, and thus will be more likely to respond.
From a potential responder’s perspective, it will be more likely to respond (is more
susceptible) if it is strongly connected to the initiator (short path length), and if it
has neighbors which are strongly connected to each other (high clustering coefficient).
The mechanism by which shortest paths function is simple to understand: shortest
paths represent most probable paths. In fact, in our data we find that, even when controlling for Euclidean distance, the shortest path length retains significant power for
predicting responders (logistic regression, (df = 2, N = 5502), p < 0.00001). We also
find that a strong interconnection among an individual’s neighbors (high clustering
coefficient) increases the probability of response (beyond the first responder, Section
3.4.E, Fig. S3.8). This is due to the fact that multiple pathways allow for reinforcement of observations, increasing the likelihood of behavioral change. We find that it
is specifically the structural properties of the interaction network, more so than other
correlated variables (such as the number of visible neighbors, local density, weighted
out-degree, proximity to the edge of the group), that account for individual social
influence and susceptibility within groups (Fig. S3.6, Section 3.4.D).
When considering responders, we find that if we consider only the first responder, the interconnectedness of its neighbors is not important (logistic regression,
df = 2, N = 5502, p = 0.5). However when considering subsequent responders the
interconnectedness of their neighborhoods becomes important in explaining susceptibility (logistic regression; df = 2, N = 5502, p < 0.00001) to influence. This is
because in the latter case, multiple pathways allow for reinforcement of observations,
increasing the likelihood of behavioral change, whereas in the case of the first respon43
der the clustering of its neighborhood does not offer such a possibility since it does
not experience, and thus cannot benefit from a neighborhood structure that benefits,
multiple exposures. For all candidate first responders, it is the strength of the connection between initiator and candidates that dictates the likelihood of response. Thus,
with the exception of the first responder, the same network properties that determine
an individual’s influence also determine its susceptibility to influence.
44
B
1
10
simulation:BmagnitudeBof
behavioralBcascade
experiment:BmagnitudeBof
behavioralBcascade
A
0
10
40
0
10
−5
10
0
local weighted clustering coefficient
immediately preceding event initiation
−5
10
local weighted clustering coefficient
immediately preceding event initiation
C
D
1
social influence
external visual field
1
0.5
standardizedBunits
standardizedBunits
fractional
numeric
simple
1
10
school
front
0
−0.5
0.5
0
−0.5
0
school
back
0.5
normalized
distance
1
school
front
0
5
10
15
20
25
distanceBfromBgroupBboundaryB8cm7
E
social influence
0.0001
1e−05
1e−06
1e−07
1e−08
Figure 3.1: Local weighted clustering coefficient predicts the magnitude of behavioral cascades. (A) Relationship between weighted directed clustering coefficient of
the initiator at time of cascade initiation, and magnitude of behavioral cascade, where the
black dots show the mean value of all cascade sizes falling in each bin, and the error bars
show the standard error. Data are divided into 8 logarithmically spaced bins. The dotted
line and blue shaded region represent the best generalized linear model fit to the data with
a log link function, and the 95% confidence interval. Histogram below is the distribution
of initiator weighted directed clustering coefficients. (B) Relationship between simulated
cascade size and local weighted clustering coefficient, for three contagion models (fractional,
solid line; numeric, dash-dot line; simple, dashed line). Solid and dotted lines are regression
fits to the simulated data, while shaded regions are 95% confidence intervales. Here we use
the best fit parameters from each model, and simulate behavioral cascades on all nodes in
20 different fish configurations. See section 3.4.F for model details. (C) Social influence
(red triangles, approximated by weighted directed clustering coefficient), and visual field
that extends outside the group (blue squares), plotted against normalized distance from
school back to school front, in polarized fish groups; Section 3.4.A. (D) Same as (C), plotted against distance from group boundary, again calculated for polarized groups. Units in
(C) and (D) standardized such that mean is zero and standard deviation is 1. (E) Example
of spatial structure of social influence in a school, school size of 149 fish.
45
3.2
The nature of the contagious process
Considerations of social influence, and susceptibility to social influence (as described
above), both suggest that an individual’s activation probability increases with reinforcement from others, but is inhibited by having a large number of visible neighbors
who have not responded, and that redundant paths amplify information. To investigate further the nature of the contagious process, and to explore further why
we see a strong positive relationship between clustering coefficient and influence, we
simulate three types of contagion processes on the fish sensory interaction networks:
(1) a simple contagion model, where response probability depends only on the link
weight connecting the focal node to the active (startled) node; (2) a numeric threshold model, where the response probability of the focal node depends on the absolute
number of observed active nodes; and (3) a fractional threshold model, where the
response probability depends on the fraction of observed active neighbors. Details of
the simulations are included in the supplemental materials (Section 3.4.F).
All three of these models have been used in different contexts. Simple contagion models have been used to study the propagation of disease on networks
[Keeling, 1999], numeric threshold models have been used in the study of selforganized criticality [Bak et al., 1987, Sachtjen et al., 2000], and bootstrap percolation [Adler, 1991, Solomon et al., 2000], and the fractional threshold model was used
to examine the origin of large, rare cascades in random networks [Watts, 2002], Both
the numeric threshold and fractional threshold contagion are models of complex
contagion, meaning that propagation becomes more likely when multiple responses
are observed, in contrast to simple contagion, where propagation only depends on
the link connecting two nodes. Only in the fractional threshold model, however,
is inhibition from inactive neighbors taken into account, in that the likelihood of
propagation is related to the fraction of neighbors who have responded.
46
We simulate transmission events on each of our experimentally derived networks,
which start from the same initiator individuals (nodes) that occur in our data (optimizing over 5000 sets of parameter values, with 50 replicate simulations for each
set of parameter values on each network). We compare this simulated transmission
with our experimental data in two ways: (1) we ask how well the simulation predicts the magnitude of the behavioral cascade, and (2) how well it predicts the actual
identities of responders. We compare the performance (maximum likelihood) of each
candidate contagion model, along with the performance on randomized versions of
the network. The randomization is achieved by randomly permuting the identities
of the nodes in the network, so that the spatial embedding is removed (nodes which
are spatially close are no longer strongly connected in the randomized version). We
use these models on the randomized identity networks as a baseline with which to
compare model performances using the real responder identities (Section 3.4.F).
Both approaches find that the transmission of behavior in our fish schools is best
described by a fractional contagious process (Fig. S3.10, Fig. S3.11, Section 3.4.F).
The fractional threshold model also has the strongest relationship with local clustering
coefficient, in which we find that the most influential nodes are those with high
clustering coefficients (holding out-degree and local density constant), the alternative
processes of contagion (simple and numeric) cannot account for the experimental
relationship between network properties and social influence (Fig. 3.1B).
The fractional nature of contagion also accounts for the relationship between spatial position and influence (Fig. 3.1C–E). Individuals near the front and side edges
of the group tend to have fewer visible neighbors than those near the group center
(those with high centrality). For a fractional contagion process, the activity of a small
number of neighbors near such periphery individuals can therefore have a relatively
large effect. The barrier to startling is lowered still by the fact that neighbors of these
peripheral individuals tend to have more interconnections than do neighbors of more
47
central nodes (Fig. 3.1D), leading to higher levels of reinforcement via observation of
multiple fast-start evasion maneuvers.
3.3
Conclusions
We demonstrate here that it is possible to reconstruct quantitative interaction networks based on the sensory information employed by organisms in mobile groups
when making movement decisions. Unlike abstracted representations of connectivity,
such as when it is assumed that individuals in close proximity are “connected” in a
social network, here the strength of connections (their weight represented by color
in Fig. 2.5A,B) actually represent the probability of responding to others, should
they change behavioral state. Individual fish employ simple, likely fast to compute
[Koch and Ullman, 1987], measures to assess behavioral change, which would give
the social behavior robustness in the face of factors such as changing light conditions, or turbidity, when compared to other possible visual measures. Fish appear
not to regulate their response based on the exact nature of fast-starts (that is they
tend to perceive others as having either responded, or not responded, in a binary
manner). This may be because it is too time consuming to make a more complex
assessment, and/or because a fast-start is inherently ambiguous, typically being a
reflex response. The consequence of this is that, as has been predicted theoretically
[Giraldeau et al., 2002], false alarms may be an inevitable corollary of adaptive collective response in a dangerous and uncertain world. This is consistent with previous observations of animal groups which found false alarms to be common, often accounting
for a large proportion of overall alarms [Cresswell et al., 2000, Blumstein et al., 2004,
Kahlert, 2006, Beauchamp, 2010, Beauchamp, 2012] .
Analysis of the local structural properties of the interaction networks within
schools proves to be highly informative, allowing us to understand, and thus to pre48
dict, how behavioral change is transmitted through groups, as well as to explain a
subtle, but important, distinction between an individual’s social influence, and its
susceptibility to social influence. We reveal the vital role of strong interconnection
among an individual’s neighbors (high local clustering coefficient). This local property of the network increases the probability of response due to the fact that multiple
pathways allow for reinforcement of observations, increasing the likelihood of behavioral change. Such properties would not be observed for simple contagion, and our
experimental data, and simulations, provide strong evidence for behavioral transitions
spreading as complex fractional contagion (response not to a fixed threshold number,
or “quorum”, but rather to the fraction of neighbors who have responded, regardless
of the number of neighbors).
We find that individuals who find themselves near the leading and side edges of the
group tend to be the most socially influential, and most susceptible to social influence.
In addition, individuals in such positions have best access to visual cues outside
the group, adding to their informational importance. This may counteract costs
associated with their proximity to the group boundary, such as increased exposure to
the risks of predation, and may help explain why some predators preferentially attack
highly coordinated groups from their rear [Handegard et al., 2012].
Thus, despite the inherently probabilistic nature of individual behavior, revealing
the structure of sensory networks of interaction allows us to identify which individuals
in a group are the most influential, and to predict cascades of behavioral change
at their moment of initiation. This demonstrates that a precise and quantitative
understanding of collective behavior in large, self-organized animal groups can be
achieved if the hidden pathways of communication are effectively revealed.
49
3.4
3.4.A
Technical details
Group spatial measures
Local density
We draw a circle of radius R around the focal fish, and count the number of neighboring fish falling in that circle. We divide this number by the area of the circle that
falls within the group boundary, to get rid of unwanted edge effects.
PN
ρi =
j=1 (dij
< R)
Ci ∩ Agroup
,
(S3.1)
Where dij is the distance from fish i to fish j, R is the local density radius, Ci is the
circle of radius R around fish i, and Agroup is the area of the group. To make sure
the boundary density is not being affected by having a tiny number of fish fall within
the local radius, R, we define this local radius as the distance to the N th nearest
neighbor, thereby ensuring that at least N neighbors fall within R. An example of
local density is shown in a sample fish group in Fig. S3.1A.
Distance from group boundary
We detected the group boundary using an alpha-shape method. A boundary edge
is drawn between any two fish, Pi and Pj in the group where a disk of radius 1/a
can be placed such that it contains all points in the group and Pi and Pj are on
its boundary. We use the alpha shapes implementation available on the matlab file
exchange platform [Lundgren, 2010], with a radius of 25 cm. Once we have identified
the boundary individuals, we define a fish’s minimum distance from the boundary as
the minimum distance to any point on a line connecting two consecutive boundary
fish. Fish located on the boundary have a minimum distance to the boundary of
50
zero. An example of distance from group boundary is shown in a sample fish group
in Fig. S3.1B.
Distance front-back
We calculate the heading direction of the group as the mean heading of all fish. We
then measure the distance front-back as the shortest distance from each fish’s position
to the line drawn through the group center of mass along the group heading. A fish
has a distance front-back of 0 if there are no other fish further behind the group center
of mass. Likewise, a fish has a distance front-back of 1 if no other fish are further
ahead of the group center of mass.
di = (r¯i − r̄COM ) · |H̄|
(S3.2)
Where H̄ is the mean group heading, ri is the fish position, and rCOM is the group
center of mass position. An example of distance front-back is shown in a sample fish
group in Fig. S3.1C.
51
Figure S3.1: Illustration of spatial measures for three sample groups. (A) local
density, (B) distance from group boundary, and (C) distance from the group front.
52
Polarized/Rotational states
We calculate polarized and rotational states following Tunström et al. [Tunstrøm et al., 2013].
The polarization of the group is a value between 0 and 1, where 0 is completely
unpolarized, and 1 is completely polarized. The polarization is the average of all
individual headings (ui ). A group is categorized as polarized if it has a value of
P > .65 and R < .35, and it is categorized as rotating if P < .35 and R > .65.
Pgroup
v
u N
X
1u
= t
u2 + u2yi
N i=1 xi
(S3.3)
The amount of rotation around the group center of mass position is calculated as
follows:
Rgroup
1
=
N
N
X
ui × ri (S3.4)
i=1
Where a group is completely rotational if R = 1.
3.4.B
Spatial structure of the network
Although many global network measures that exist in the literature have been useful
in characterization of these networks, such as efficiency measures and centrality measures, we find that such global metrics are not as applicable to complex contagion
processes as they are to simple contagion processes. These global measures are also
useful when information often reaches the entire network, unlike in our data where information spread is most often localized to a small portion of the network, although it
can spread extensively. Instead we focus our investigation on local network measures.
We include distributions of some of these local measures in Fig. S3.2.
53
Figure S3.2: Descriptive network statistics. (A) Frequency of in-degree; inset: indegree and out-degree are correlated, but not perfectly so. (B) Frequency of number
of neighbors (based on a binary threshold of AA > 0.02 radians). (C) Cumulative
distribution function of the difference between incoming and outgoing edge weights.
Although many edge weights have similar in and out values, there are enough differences that symmetrizing the network and treating edges as undirected is not justified
in general.
A simple measure of path redundancy is the local clustering coefficient, which
compares the number of closed triangles (where all three neighbors are connected)
to the number of total possible triangles (all pairs of neighbors with the focal node).
We use the weighted directed clustering coefficient from [Fagiolo, 2007] as a simple
extension of the clustering coefficient concept to weighted and directed graphs, which
is defined as the following;
P P
1 j k (wij + wji )(wij + wki )(wjk + wkj )
Ci =
tot
↔
2
[dtot
i (di − 1) − 2di ]
X
X
out
din
=
a
,
d
=
aij
ji
i
i
j6=i
(S3.5)
(S3.6)
j6=i
in
out
dtot
i = di + di
X
d↔
=
aji aij
i
(S3.7)
(S3.8)
j6=i
Where wij is an entry in the weighted adjacency matrix, and aij is an entry in the
binary adjacency matrix. We found that some characteristic network properties are
spatially structured. Most notably, the local clustering coefficient tends to be found
54
in regions of local density, near the boundary. In polarized groups, fish with high
clustering coefficients tend to be found near the front of the group (Fig. S3.3, Fig. 3.1,
main text).
Figure S3.3: Spatial structure of local density and social influence. (A)Local density
is not strongly correlated with distance from boundary, but (B) density does tend to
be higher near the group front. (C) social influence (approximated by local clustering
coefficient) varies positively with local density.
Additionally, we find that there exist regions of local clustering in the fish schools
(Fig. S3.4). We calculate the correlation function for fluctuations in clustering coefficient as a function of distance. The correlation length, which is the distance beyond
which there tends to be no similarity between pairs of clustering coefficients, averaged
over all nodes in all networks in our dataset is 32 cm.
P
ψ(r) =
i,j
δ(r < dij < r + )(CCi − |CC|) · (CCj − |CC|)
P
i,j δ(r < dij < r + )
(S3.9)
Where |CC| is the average value of local clustering coefficient over all nodes in the
network, so that we are measuring fluctuations from the mean value, δ is an indicator
function that selects pairs of nodes i and j such that the distance between them is
between r and r + . The denominator is just the normalization by number of node
pairs which are this distance r apart.
55
0.25
real(network
randomized(network
correlation(length
0.2
correlation(function
0.15
0.1
0.05
0
−0.05
−0.1
0
10
20
30
40
50
60
distance(between(nodes(i(and(j((cm)
70
80
Figure S3.4: Correlation function between clustering coefficients of pairs of nodes, as
a function of the distance between them. Nodes which are close together tend to have
more similar values of local clustering coefficients than nodes which are far apart.
3.4.C
Exterior visual field
Access to visual information outside the group is important for obtaining direct visual
information regarding predators. We evaluate this exterior visual field of each fish
by measuring the fraction of its field of view which extends to the tank boundaries
without intersecting with a neighboring fish. As expected, we find that the fraction
of the field of view extending outside the group is much higher for fish located on the
group boundary, as well as those near the group front and group rear (Fig. 3.1C,D
main text, Fig. S3.5).
56
Figure S3.5: Visualization of external field of view; rays (white lines) unobstructed
by other fish are drawn from eye positions to corresponding positions on the tank
wall. Fish on the group boundary can see much more outside the group than those
in the group interior.
3.4.D
Fitting local clustering coefficient model of cascade
size
We determine how well the local clustering coefficient of the initiator describes the
resulting cascade size by fitting a generalized linear model with log link function
(Poisson regression because our dependent variable is a count variable), to the data
using the local weighted clustering coefficient of the initiator as a predictor variable.
We find that this model is highly significant (likelihood ratio test, χ2 = 78.1, (df =
1, N = 138), p < 0.00001). To test whether this effect is truly due to the network
structure, or results from other correlated quantities, such as local density, out-degree,
and distance from group boundary, we estimate the relative contribution of each
feature in two ways (Fig. S3.6). First we evaluate the relative importance of each
feature using L1 regularization. We find that, at the point of minimum deviance
and minimum deviance plus 1/2 standard error, the local clustering coefficient has
57
the highest relative weight. Second, we evaluate each variable’s unique contribution
to the explained variance of the cascade size distribution [Meyers et al., 2006]. The
local clustering coefficient has the highest unique contribution to the 69% of variance
explained by the total model. The sum of the values in Fig. S3.6B is smaller than
69% due to correlations among the variables. The strong performance of clustering
coefficient in these models suggests that the interconnectedness of a node’s neighbors
is essential to understanding the way in which alarm cascades flow through a group.
B
A
numberyofyneighbors
numberyofyneighbors
distanceyfromyboundary
distanceyfromyboundary
localydensity
localydensity
maxyout-link
maxyout-link
weightedyout-degree
weightedyout-degree
clusteringycoefficient
clusteringycoefficient
0
0.1
0.2
0.3
featureyweight
minimumydeviancey+y1/2ySE
0.4
0
0.02
0.04
0.06
0.08
0.1
uniqueycontributionytoyexplainedyvariancey
Figure S3.6: Relative contribution of each variable to the model’s fit to the cascade size distribution. (A) relative feature weight from L1 regularization. (B) Each
variable’s unique contribution to total variance explained by the model. The local
clustering coefficient has the highest contribution to the model in both cases.
Along with fitting models to the cascade size data, we also plot the mean and
standard error of cascade size as a function of clustering coefficient. We assign the
data to logarithmically spaced bins, of sizes such that all but the very smallest contain
four or more data points, and calculate the mean and standard error of the cascade
sizes in each bin. We find that most of the mean values fall within the 95% confidence
interval for the generalized linear model fit (blue shaded region, Fig. 3.1A, main text).
Alternative versions of weighted, directed local clustering coefficient
The generalization from the binary undirected clustering coefficient to weighted directed clustering coefficient (W DCC) is not entirely obvious. We use the most widely
58
accepted definition [Fagiolo, 2007] for the majority of the analysis. However, one
drawback to this version is that it does not distinguish between in-connections and
out-connections, but rather sums the in-strength and out-strength. For the purposes
of considering influential and susceptible individuals this may not be the best version
to use. Only the out-connections of a focal node, and not the in-connections will have
bearing on that node’s influence. Likewise, only the in-connections will be relevant
to a node’s susceptibility. For this reason we propose two slightly modified versions
of the weighted clustering coefficient, defined as follows.
Ciout =
1 XX
tot
↔
(wji )(wki )(wjk + wkj )/[dtot
i (di − 1) − 2di ]
2 j k
(S3.10)
Where Ciout is equivalent to (10), except for removing the in-connections from the
focal node.
Ciin =
1 XX
tot
↔
(wij )(wik )(wjk + wkj )/[dtot
i (di − 1) − 2di ]
2 j k
(S3.11)
Where Ciin is also equivalent to (10), now removing the out-connections from the
focal node and only keeping the in-connections. Using the modified C out results
in a slightly better model when fit to the cascade size data (likelihood ratio test,
χ2 = 94.7, (df = 1, N = 138), p ≤ 0.00001, compared to χ2 = 78.1, (df = 1, N =
138), p ≤ 0.00001 for the unmodified W DCC), as well as a slight improvement when
we look at responder susceptibility (see the following section), (likelihood ratio test,
χ2 = 57.8, (df = 1, N = 5502), p ≤ 0.00001, compared to χ2 = 57.7, (df = 1, N =
5502), p ≤ 0.00001). The difference between in-connections and out-connections tends
to be relatively small compared to the weight of the connections themselves, which is
why the effect of using these clustering coefficient variations is not tremendous.
59
Alternative network measures
We compare the local clustering coefficient to network measures that have been used in
different contexts to identify influential nodes [Barthelemy, 2004, Kitsak et al., 2010,
Pei et al., 2014].
A summary of the results is included in Table S3.7.
In un-
weighted networks, the k-core is the set of nodes which have at least k-connections
[Dorogovtsev et al., 2006]. A node’s k-core number is the largest k-core to which it
belongs. As a definitive weighted version of the k-core hasn’t yet been established,
we approximate a weighted k-core by averaging k-core numbers of each node in unweighted networks created by binarizing the weighted network at different threshold
values. We also normalize the k-core numbers such that they fall between 0 and 1,
to account for networks that have been binarized using higher thresholds having a
smaller number of connections, and thus a smaller number of cores. The unweighted
k-core number of initiator nodes is negatively correlated with resulting cascade size.
This is likely due to the fact that central nodes tend to have higher k-core numbers,
and these central nodes experience more reaction inhibition, as a result of having
more neighbors. The model using weighted k-core number of initiator nodes is not
significant, (GLM fit, likelihood ratio test, p = 0.179).
We also consider the weighted betweenness centrality of the initiator nodes as
predictors of cascade size, given by:
g(i) =
X σst (i)
σst
s6=i6=t
(S3.12)
Where σst is the number of shortest paths passing from node s to node t, and σst (i)
is the subset of those paths which also pass through node i. Nodes with high betweenness centrality are less likely to be influential. This gives more support to the
hypothesis that cascades in our system spread through fractional contagion. Nodes
60
with high betweenness centrality are located in the center part of the group, where
they are more inhibited from responding due to having a larger number of neighbors.
The edge weight variation is calculated as the variance in edge connections per
node, as a measure of how homogeneous or heterogeneous a certain node’s connections
are. This measure is also a negative predictor of cascade size, but the significance is
most likely explained by spatial effects (central nodes tend to have more edge-weight
variation).
��
Network feature
Coefficient
p value
LocalWweightedWclustering
coefficient
0.37
78.1
<0.0001
K-core
-0.11
9.05
0.0026
Weighted K-core
0.05
1.81
0.179
Weighted betweenness
-0.35
52.4
<0.0001
EdgeWweightWvariation
0.14
15.4
0.0008
Table S3.7: Comparison of model performance using alternative network measures.
3.4.E
Susceptibility
We find that the same network properties that determine an individual’s influence also
determine its susceptibility, where we define susceptibility as an individual’s likelihood
of responding to an observed startle event. This concept of susceptibility results from
looking at responses from an individual’s perspective, while influence results from
looking at responses from an initiator’s perspective. An individual will be more likely
to respond (is more susceptible) if it is strongly connected to the initiator (short
path length), and if it has neighbors which are strongly connected to each other
(high clustering coefficient). The mechanism by which shortest paths are important
is simple to understand: shortest paths represent most probable paths. In fact, in our
61
data we find that, even when controlling for Euclidean distance, the shortest path
length retains significant power for predicting responders (logistic regression, (df =
2, N = 5502), p < 0.00001). We also find that a node is more susceptible if it is in a
highly clustered region, that is, if it is in a region with strong interconnections between
neighbors. A node that is a member of a highly clustered region can take advantage
of observing multiple responders. In a simple contagion process, these multiple ties
would be redundant, but in complex contagion, they are helpful in instigating a
response. Here we define susceptibility as the likelihood of a fish responding given
that it observes the initiator.
We measure the likelihood of late responders (all responders after the first responder) as a function of clustering coefficient (controlling for in-degree, local density,
distance from boundary, and edge-weight connecting observing fish to startled fish),
and find that nodes with a higher W DCC are more likely to respond to an observed
startle event than those with a lower W DCC (logistic regression, (df = 5, N = 5502),
p < 0.00001). We consider only subsequent responders (after the first) because these
are the fish that have the chance to benefit from being in a highly clustered region by
observing multiple responders. When we repeat the same analysis for just the first
responder, the clustering coefficient has no significant effect in the model (logistic
regression, (df = 5, N = 5502), p = 0.50). Since the first responder does not benefit
from multiple exposures, we only expect to see the complex propagation effect in
individuals responding after the first responder. Both perspectives support the main
point, that redundant paths amplify information, by asking on the one hand, “given
the local neighborhood of an initiator, how many fish are likely to respond?”, and on
the other hand, “given the local neighborhood of any potential responder, how likely
is a response from that fish?”.
We use k-fold cross-validation to prevent overfitting when evaluating influence and
susceptibility effects. This was done by partitioning the data into k = 10 subsets,
62
training a model on k − 1 subsets, and testing it on the remaining subset left out
of the training sample. This procedure is repeated k times in total, withholding a
different test subset each time. We plot the values predicted from the training model
on the test data and compare with observed data (Fig. S3.8). The range of predicted
values falls within the 95% confidence interval for the best-fit model using all the
data. This gives us confidence that we are not overfitting to the current dataset, and
that we will observe the effect in a new set of data.
Figure S3.8: Cross validating the model. The data was repeatedly divided into 10
subsets, trained on 9 of them, and tested on the remaining one (plotted in black diamonds). 95% confidence interval from model fits are shown in the blue shaded areas.
(A) Cross-validation for the relationship between cascade size and local clustering
coefficient. (B) Cross-validation for relationship between susceptibility to influence
and local clustering coefficient. Histograms show the frequencies of positive responses
(above) and no responses (below) as functions of clustering coefficient.
3.4.F
Simple, numeric, and fractional model simulations
We simulated complex and simple propagation on our estimated fish interaction networks, using both numeric thresholds and fractional thresholds. Activation using
numeric thresholds requires an absolute number of neighbor activations, while activa63
tion using fractional thresholds requires that a fraction of neighbors become activated.
Because our data only contain cascades initiated by a single fish, we only consider
cascade simulations seeded by one node. We activate a seed node i at time t0 . In
the following time-step each neighbor, j, of i has a chance of observing i as active,
with probability wij , which is the strength of the link between the two nodes. Node
j itself becomes active when the number or fraction of observed active nodes exceeds
a threshold, for numeric and fractional contagions, respectively. To compensate for
the increased barrier to activation with an increasing threshold, we include an amplification parameter C, which represents the number of chances each node has to
observe a neighbor’s activation in each time-step t (alternatively the number of coinflips weighted by wij , where at least one positive result is needed to connect nodes
i and j). The simulation continues until no more activations occur. Seeding the
cascade with only a single node presents a problem for homogeneous-threshold based
contagions, however, because any threshold greater than 1 requires that more than
one fish begin in the active state. Consequently, we assign thresholds randomly on
nodes in the graph, selected uniformly (although our results do not depend on this
form of distribution, as we will show in the following section) from integers [1, φnmax ].
For example, for an average threshold of 3, there will be equal numbers of nodes with
thresholds 1, 2, 3, 4, or 5. Some nodes with threshold 1 are necessary to allow the
cascade to propagate beyond the seed node. Similarly for fractional thresholds, we
distribute thresholds uniformly from (0, φfmax ] on the graph.
The simulation proceeds as follows, continuing over t time-steps, until no more
activations occur. Only one activation (Ai changes state from 0 to 1) is allowed per
node, per simulation, and any node that becomes active during the course of the
64
simulation remains active for the remainder.
Ai (t) = 1 if



 Ai (t − 1) = 1
(S3.13)

PN


j=1 δi,j > φi
δi,j = (wij Aj (t − 1) > min(Uj ))
(S3.14)
Where Uj is a vector of Cf or Cn uniform random numbers, N is the total number of
fish in the school, and φi is a threshold defined by the type of contagion used in the
model.
φi,simple = 1,
φi,numeric = ni ,
φi,fraction = ni
X
Mij
(S3.15)
j
Where ni is the particular activation threshold of node i, Mij is the binary adjacency
P
matrix, and j Mij is the number of neighbors of node i.
We estimate the prior distributions of parameters as follows. Clearly, φsimple = 1,
because simple contagions only require one connection to propagate. The average
node in our graph has 30 visible neighbors, so we allow φnmax to be drawn from the
uniform distribution [1,30]. φfmax is drawn from the uniform distribution (0,1], so
that on average the number of neighbors required for activation at the highest value
of φfmax equals the number required for the highest value of φnmax . The coin-flip
parameter, C, is intended to compensate for impedance to propagation due to higher
thresholds. Preliminary analysis reveals that across all values of φnmax and φfmax , the
optimal value for Cn is twice the value of φnmax , and Cf is twice the value of φfmax × 30.
So we set Cn to be 2 × φnmax , and Cf to 2 × φnmax × 30 for the remainder of our
analysis (Fig. S3.9).
65
Figure S3.9: Log likelihood for (A) numeric and (B) fractional models are maximized
at a ratio t/C of 0.5, and 0.5/30, respectively.
To evaluate the performance of the three contagion models–simple, numeric, or
fractional–on our derived interaction networks, we seed the propagation simulation
with the initiator nodes (individuals) observed from the cascade data, and evaluate
the results using two scoring metrics. First, we find the model that best predicts the
identities of the responder nodes. We calculate the maximum likelihood over a random
selection of parameter values, drawn from the prior distributions, averaged over 10
parameter sets of 500 parameters each. For each set of parameters for each node
we averaged 50 iterations of the simulation to remove any effect from the random
distribution of thresholds on the graph. We compare our three candidate model
performances on our fish network to the performances on randomized versions of the
network. The randomization is achieved by swapping the identities of responder nodes
with the same number of randomly selected nodes in the group. The randomized
models provide a baseline for comparing the simple, numeric, and fractional contagion
models, and all three candidate models at a minimum clearly better predict responder
identity than the baseline (Fig. S3.10).
66
We also calculate the marginal likelihood, which is the integration of the likelihoods over the prior distribution of φ, and find qualitatively the same results. The
marginal likelihood is calculated because a reasonable model should not depend on
fine-tuning of parameters (Fig. S3.10B) [Mann et al., 2013].
B
−2360
marginaltlikelihood
−2300
−2380
−2400
−2420
−2410
−2420
−2430
al
nd
io
n
ic
ra
ct
et
pl
si
m
fra
er
−2400
−2440
fra
m
pl
e
nu
m
si
nd
al
tra
nd
ra
io
n
fra
ct
pl
et
m
si
om
−2440
−2390
om
ltr
an
do
m
si
m
pl
e
nu
m
er
ic
fra
ct
io
na
l
−2360
−2380
na
−2340
−2370
ct
io
−2320
om
maximumtlog-likelihood
A
Figure S3.10: (A) Results from the simple, complex numeric, and complex fractional
contagion model fits. The models are evaluated by how well they predict the identities of the responded individuals, seeded by the initiators from the data. We show
the mean values of the maximum likelihoods, averaged over 10 iterations of 500 sample parameter values. Simple random and fractional random are the results from
randomized versions of the data, and represent the baseline maximum log-likelihood.
Both complex contagion models (numeric and fractional) fit our data better than the
simple contagion model. (B) Marginal likelihoods for candidate models (alternative
to maximum likelihood (A)). Error bars indicate standard error (note that some are
smaller than the marker size). Marginal likelihoods were calculated by averaging the
log-likelihoods from 10 trials using 500 parameters in each trial, fit to the responder
identity data.
Next, we consider an alternative scoring metric, where a model is evaluated by how
accurately it predicts the resulting cascade length. The best fit model for predicting
the identities of the responders is also the best fit model for predicting the size of
the cascade (Fig. S3.11). The best fit threshold for the fractional complex model
67
is approximately φf = 0.25, with agreement between both scoring metrics. Best fit
values for both scoring metrics are included in Table S3.12. While the numeric and
fractional models perform similarly using the first scoring metric, the numeric model
is much worse at predicting cascade size, with a negative correlation persisting for
most of the threshold range. In the fractional model, however, we see a strong positive
correlation between predicted cascade size and experimentally observed cascade size.
This is likely due to the difference in the two contagion types near the group boundary.
The numeric threshold will predict smaller cascades from the boundary than from the
bulk, due to boundary nodes having fewer neighbors. In fractional contagion, cascades
from the boundary will be larger relative to the bulk, for the same reason (boundary
nodes have fewer neighbors).
68
A
B
model4predicting4cascade4size
marginal4likelihood4±4SE
model4predicting4identities
marginal4likelihood4±4SE
−2340
−2350
−2360
−2370
−2380
−2390
−2400
−2410
−2420
0
0.2
0.4
0.6
0.8
1
complex4fractional4φmax
C
model4predicting4cascade4size
marginal4likelihood4±4SE
model4predicting4identities
marginal4likelihood4±4SE
−2360
−2370
−2380
−2390
−2400
−2410
−2420
0
10
20
30
complex4numeric4φmax
−696
−698
−700
−702
−704
−706
−708
−710
−712
0
40
0.2
0.4
0.6
0.8
1
complex4fractional4φmax
D
−2350
−2430
−694
−696
−698
−700
−702
−704
−706
positive4correlation4with
response4variable
negative4correlation4with
response4variable
−708
−710
−712
−714
0
10
20
30
complex4numeric4φmax
40
Figure S3.11: Optimal threshold values (φ), organized by row for complex fractional
(A, B) and complex numeric (C, D) contagion models, and by column for identity
(A, C) and cascade size (B, D) model scoring metrics.
We purposefully chose to use basic models to allow us general insights into complex
and simple propagation events evolving on our networks, and to determine if there
exists support for one contagion mechanism over another. Further work is needed to
fully understand the interplay between the structure of the network and the type of
contagion mechanism by simulating cascades on specific types of networks.
In order to determine which of the considered models best supports the finding
that the final cascade size increases with the local clustering coefficient of the initiator, we simulate the fractional, numeric and simple threshold models initiated from
randomly selected nodes throughout the network, using the best-fit parameters for
69
Complex fractional model
Parameter
Mean
SE
φfmax ; identities
0.2314
0.0172
φfmax ; cascade size
0.2438
0.0270
Complex numeric model
Parameter
Mean
SE
φnmax ; identities
8.900
0.3797
φnmax ; cascade size
24.57
0.5102
Table S3.12: Optimal parameter values for both evaluation metrics and both contagion types. Optimal parameter values do not coincide exactly with marginal likelihood parameter values. The optimal values were obtained by finding the location
of the maximum log-likelihood in each of the 10 simulation runs, while the marginal
likelihood values in Fig. S3.11 are obtained by averaging every log-likelihood value
from each value of φmax on the x-axis. The marginal likelihood parameters are best
on average, while the optimal parameter values are best overall. The difference between them is accounted for by there being more spread in the data near the optimal
parameter values.
φ and C (Fig. 3.1B, main text). When controlling for out-degree and local density
(these are both correlated with clustering coefficient, and can affect the final cascade
size), we find that simulated cascade size correlates strongly with local clustering coefficient only in the fractional threshold model. In the numeric model, cascade size
is slightly correlated positively with clustering coefficient, but not as strongly as in
the fractional model. In the simple model, cascade size is essentially independent of
clustering coefficient, when controlling for out-degree and local density. These results
are robust to the choice of response threshold distribution (Fig. S3.15A,B), as we
describe in the next section.
70
Alternative response threshold distributions
We consider alternative threshold distributions in this section, since we do not know
the form of variation arising from sensory noise, or from individual variability in
sensitivity to visual cues, for example. We explore the effects on simulated cascade
size when we relax the constraint of a uniform distribution of response thresholds
in the group, to see if the effects we observe are consistent with a wider range of
threshold distributions. We now set the response thresholds of each individual by
drawing from a Gaussian distribution with mean φn,f and standard deviation σn,f ,
while restricting the distribution of thresholds to be greater than zero.
We repeat the analysis described above, this time using a Gaussian distribution
for response thresholds, to find values for φn,f and σn,f which best predict responder
identities and experimental cascade size. In Fig. S3.13, we plot model performance
as heatmaps with the variation in response threshold on the x-axis, and the mean
response threshold on the y-axis (this is a 2-dimensional version of Fig. S3.11). We
find that, for fractional contagion, models that best predict both cascade size and
responder identities have low values of mean and standard deviation (Fig. S3.13A,B).
For numeric contagion, models that best predict responder identities also have low
values of mean and standard deviation (Fig. S3.13C), however numeric contagion is
a poor predictor of cascade size, as for much of parameter space simulated cascade
size is negatively correlated with experimental cascade size (Fig. S3.13D), for spatial
reasons we described in the preceding section (see Fig. S3.11D). We also find that
the maximum likelihood is greatest for fractional contagion models irrespective of the
type of distribution (Fig. S3.14A and Fig. S3.10 for normal and uniform distributions,
respectively).
71
φf
σf
σf
D
LLH: Cascade size
φn
LLH: Cascade size
LLH: Responder ID
φf
C
Numeric contagion
B
LLH: Responder ID
Fractional contagion
A
φn
σn
σn
Figure S3.13: Gaussian threshold distribution: best fit parameters. Heatmaps show
the average log-likelihood value at each value of φ and σ (5000 total parameter values). This figure is a 2-dimensional version of Figure S3.13. (A) Fractional contagion,
model predicting responder identities. (B) Fractional contagion, model predicting
cascade size. (C) Numeric contagion, model predicting responder identities. (D) Numeric contagion, model predicting cascade size. Simulated cascade size is negatively
correlated with experimental cascade size, for numeric contagion, which explains why
there is no clear optimal parameter region here.
Next, we chose parameters from the optimal region of phase space, for the three
types of contagion (φf = 0.05, σf = 0.06; φn = 1.5 σn = 1; φs = 1, σs = 0), and
simulated cascades initiated by every node, in 20 randomly selected networks. We
confirm our results from Fig. 3.1B (main text); that cascade size varies most strongly
with clustering coefficient when using fractional contagion, as compared with numeric
or simple contagion (Fig. S3.14B). Note that the mean value of the cascade size for
numeric contagion here is higher than in Figure 4B, main text. This is because the
72
parameters that best fit the data for the Gaussian distribution result in a slightly
higher average cascade size, for numeric contagion.
A
B
−2300
1
10
simulation::magnitude:of:
behavioral:cascade
maximum:log−likelihood
−2280
−2320
−2340
−2360
−2380
−2400
−2420
−2440
numeric
fractional
simple
0
nu e
m :ra
fra eric ndo
:
m
ct
io ran
na do
l:r
an m
do
sim m
nu ple
m
fra eric
ct
io
na
l
10
−6
10
−4
10
sim
pl
local:weighted:clustering:coefficient
immediately:preceding:event:initiation
Figure S3.14: Gaussian threshold distribution results (A) Maximum log-likelihood
for models predicting responder identities in experimental cascades. Mean values of
the maximum likelihoods are shown, with standard errors averaged over 10 iterations
of 500 sample parameter values (note that error bars are smaller than marker size).
Simple random, numeric random, and fractional random are the results from randomized versions of the data, and represent the baseline maximum log-likelihood. Similar
to Figure S3.14, fractional contagion is the best model. (B) Relationship between
simulated cascade size and local weighted clustering coefficient, for all three contagion
models, using gaussian distributions. Solid and dotted lines are regression fits to the
simulated data, while shaded regions are 95% confidence intervals.
Finally, we performed a sensitivity analysis to the width of the Gaussian distribution, by selecting a mean threshold value and varying the standard deviation of
the distribution from which individual thresholds are chosen. We selected mean values for the distributions to be small enough to allow responses to propagate, even
for small values of σ, as well as resulting in consistent values between fractional and
numeric contagion. Since the average fish in the group has 30 visible neighbors, having a mean numeric threshold (φn ) of 3 corresponds to a mean fractional threshold
(φf ) of 0.1. From 10 randomly selected networks in our data we simulated a cascade
73
initiated from every node in the network, and calculated the resulting cascade size
as the average from five trials, where each node’s threshold is drawn independently
in each trial. For both fractional and numeric contagion as σ increases the average
cascade size increases, because larger values of σ allow for more individuals with lower
thresholds, facilitating propagation. We find that simulated cascade size scales much
more strongly with fractional contagion than with numeric contagion, invariant to the
width of the distribution from which response thresholds are drawn (Fig. S3.15A,B).
A
B
2
2
10
simulated cascade size
numeric contagion
simulated cascade size
fractional contagion
10
1
10
0
10
φf = 0.1, σf = 0.03
φf = 0.1, σf = 0.06
1
10
0
10
φn = 3, σf = 1
φn = 3, σf = 2
φf = 0.1, σf = 0.1
φn = 3, σf = 3
φf = 0.1, σf = 0.16
−6
10
φn = 3, σf = 5
−1
−1
10
10
−4
−6
10
10
clustering coefficient
−4
10
clustering coefficient
Figure S3.15: Simulated cascade size plotted against clustering coefficient of initiator,
for (A) fractional contagion, (B) numeric contagion. Four different levels of standard
deviation in threshold response distributions are tested, while the mean threshold
response is held constant. Solid lines represent best fit of cascade size to clustering
coefficient, after controlling for node out-degree and local density, while shaded regions
represent the 95% confidence interval around the line of best fit. Cascade size varies
more strongly with clustering coefficient when using fractional contagion (A) than
when using numeric contagion (B).
74
Chapter 4
Collective Oscillations in
Leptothorax acervorum
In this chapter we switch from studying collective motion in fish schools to collective
oscillations in colonies of ants. Although the two organisms are behaviorally disparate,
our motivation in studying them remains the same. As in chapters 2 through 3, we
seek to understand a complex emergent behavior by studying local interactions.
This is an ongoing project which is in collaboration with Catherine Offord. Catherine collected the data, and we worked on the analysis together. I designed and implemented the model-fitting procedure for establishing the effects of different interaction
types, as well as the simulation experiments.
4.1
Introduction
The phenomenon of synchronization is ubiquitous in nature and a fundamental trait
of many self-organized biological systems, found in networks of neurons, populations
of cells, and populations of organisms [Pikovsky et al., 2002, Strogatz, 2003]. Ant
colonies in the genus Leptothorax are known to display synchronized activity cycles,
which have drawn attention in recent years as an intriguing example of emergence
75
in self-organized animal groups [Bonabeau et al., 1997, Sumpter, 2010].
In these
cycles, individual ants synchronize their behavior, resulting in periodic oscillations in
collective activity inside the nest which peak approximately every twenty minutes. A
natural analogy can be drawn between neural systems and these ant colonies, as individual elements in each system do not display periodic behavior, but when placed in
a group synchronization emerges, resulting from local interactions between elements
[Couzin, 2009, Cole, 1991a, Rabinovich et al., 2006, Buzsáki and Draguhn, 2004].
Rhythmic oscillations in the brain have been shown to be beneficial for a number of
reasons [Buzsáki and Draguhn, 2004], including being a mechanism for periodically
elevating the system’s sensitivity to stimuli, providing discrete windows of opportunity for response without the need to exist in a permanently heightened state
of sensitivity. While it is not currently known why ant colonies synchronize (and
has in fact been posited that the synchronization is merely an inevitable result of
interactions [Cole, 1991b]), there could exist a similar advantage for ant colonies
that display periodic windows of heightened sensitivity. Specifically, synchronization could provide a mechanism for efficient foraging so that all foragers exit in a
coordinated fashion, or for brood care to ensure that brood items are not neglected
[Bonabeau et al., 1998].
Individual-based models of colony behavior have demonstrated that periodic
activity is attainable with only simple rules governing motion and the activating or
inhibiting influence of tactile interactions between ants [Cole and Cheshire, 1996,
Miramontes et al., 1993, Solé et al., 1993b]. A number of further factors affecting
ant colonies’ propensity to oscillate have been investigated theoretically and experimentally, including the total density of ants in the nest [Cole and Cheshire, 1996,
Solé et al., 1993a], the quantity and type of brood items present [Cole and Hoeg, 1996,
Cole and Cheshire, 1996], starvation [Franks et al., 1990] and colony age [Cole, 1992].
76
However, even in constant conditions, ant colonies are not consistent oscillators.
While at times they exhibit highly synchronized behavior, they also often exhibit
chaotic or randomly fluctuating collective activity, even in the absence of perturbation [Franks et al., 1990]. Despite frequent observation of this variation in coherence
of oscillations in activity levels, only a few studies have directly addressed the question
of how and why it occurs. In one such study, [Boi et al., 1999] found that preventing the return of foragers to the nest resulted in an increased coherence of periodic
oscillations in collective activity, revealing that returning forager ants disrupt colony
synchronization.
While experimental studies have mainly focused on colony-level phenomenon
(mainly out of necessity due to technical limitations), the majority of theoretical
studies have examined primarily the combination of conditions necessary to achieve
synchronized and periodic activity cycles. The most widely used type of model treats
individual ants as ‘mobile cellular automata’ (MCA) moving randomly on a lattice,
using fixed behavioral rules, which integrate interactions with neighbors to govern
transitions between active and inactive states. In these spatially explicit models, for
a nest of a given ant density, synchronization can be achieved at the collective level
for certain regions of parameter space. Outside this range of parameters, simulated
colonies exhibit less coherent behavior.
For example, [Cole and Cheshire, 1996]
noted that when active ants were less influential, simulated colonies could not
achieve periodic cycles and instead exhibited irregular fluctuations in collective
activity. However, these theoretical hypotheses about behavioral rules cannot be
tested experimentally without individual level data. Such individual data has been
previously unavailable, due to the impracticability of tracking individual ants across
the relevant time scales.
In this study, we directly investigate individual behavior and the role of interactions in maintaining activity cycles in colonies of Leptothorax acervorum. Using a
77
system of unique barcode tags, we track individuals in four separate colonies inside
their nests for a duration of approximately 36 hours, with a temporal resolution of 2
frames per second. Our method provides unprecedented volumes of behavioral data
to test impact of interactions between ants in maintaining activity cycles. Furthermore, we demonstrate the novel result that ants have flexible behavioral responses
to interactions, and that when they are strongly responsive, the activity waves inside
the nest are more coherent. Finally, in simulation experiments using a mobile cellular automaton model, we show that coherence of simulated activity cycles can be
controlled by the responsiveness of individuals, in a way that qualitatively matches
our data.
4.2
Methods
Colonies of Leptothorax acervorum were collected in June 2014 in the Alps near Aigle,
Switzerland. Four colonies of 100 ants, with between 1 and 4 queens per colony,
were created from these original colonies for tracking. Ants were tagged with unique
barcodes according to the protocol described in [Mersch et al., 2013]. Each square
tag measured 0.7mm in length (plus a 0.1mm white border) and was applied to the
thorax of all ants with superglue.
Colonies were housed in nests consisting of a foam wall pressed between two
microscope slides to create a chamber measuring 42mm x 63mm x 2mm. A 2mm
wide entrance situated in the middle of one of the longer nest walls allowed access
to a larger foraging area approximately 40cm x 20cm supplied with excess water and
food [Bhatkar and Whitcomb, 1970]. Fluon-coated Perspex walls surrounding the
microscope slides were used to prevent ants from walking over the top of the nest
from the foraging area (Fig. 4.1A). Colonies were allowed to acclimatize to the new
nest for between four and five hours before filming began.
78
For tracking, nests and foraging arenas were placed in environmentally sealed
boxes (temperature, 20◦ C, 67% relative absolute humidity) using the same set-up as
[Mersch et al., 2013] to obtain real-time tracking of individual ants at 2 frames per
second. Colonies were tracked for approximately 36 hours between 16:30 and 09:00
two days later, and the lights were turned on and off at 07:00 am and 19:00 every
day.
Time-discretized trajectories were extracted from tracked data using custom-built
software written in C++ (as described in [Mersch et al., 2013]; see Fig. 4.1B for one
frame of tracking output). Interactions between ants were inferred using orientation
and proximity measures acquired from the trajectory data. Only pairwise interactions were recorded. Ants within one average ant length (2mm) of each other, and
oriented with a total angle of less than 70 degrees between them were considered to
be interacting. A pair of interacting ants is shown in Fig. 4.1C. Subsequent analyses
of the data were performed in MATLAB.
79
Figure 4.1: Experimental setup. (A) Entire experimental arena is shown in light
blue, with water and food sources located outside the nest. Filming was conducted in
the nest area (white box). (B) Output from tracking software, in infrared, showing
colony inside the nest area. Individuals are identified from their barcodes and assigned
unique number ids. (C) A pair of ants interacting. Ants are approximately 2mm long,
and the barcodes affixed to their backs are 0.7mm on each side.
We find that ants exhibit sharp speed transitions between active and inactive
states (Fig. 4.2). We chose a speed threshold of 0.1538mm per second to distinguish
between active and inactive ants, and will conduct a sensitivity analysis on this parameter to insure that our results do not depend on it. The first seven hours of
tracking data were omitted to avoid measuring effects of acclimatization (although
[Boi et al., 1999] shows that colonies of Leptothorax acervorum begin exhibiting oscillations in activity even when only acclimatized for one hour).
80
A
B
Deactivations
0.9
colony920
colony933
colony939
colony940
0.8
0.8
0.7
speed9(mm/s)
speed9(mm/s)
0.7
0.6
0.5
0.4
0.3
0.6
0.5
0.4
0.3
0.2
0.2
0.1
0.1
0
−10
−5
0
Activations
0.9
5
0
10
time9relative9to9transition9point9(s)
−10
−5
0
5
10
time9relative9to9transition9point9(s)
Figure 4.2: Speed profiles of ants during transition period from active to inactive (A)
and from inactive to active (B). Curves are shown for each of the four colonies in the
dataset.
In order to characterize the periodicity of oscillations occurring in our ant colonies,
we define a measure for the coherence of these oscillations, as the power of the largest
frequency component in the Fourier-transformed time-series, divided by the sum of
all frequency components, as in [Cole and Cheshire, 1996]. We find that this measure
of coherence reliably selects regions of high periodicity (see Fig. 4.3, and Fig. 4.4
for a magnified view of oscillations). Since this analysis is focused on short term
oscillations (on an hourly scale), rather than long timescale fluctuations (on a daily
scale), we subsample our data; we measure the coherence of oscillations for samples
of 10,000 frames (i.e. 5,000 seconds), measured at every 1,000 frames (this introduces
some smoothing, as much of every sample overlaps with other samples). The average transition rate is 0.0036 transitions per frame, so on average there are 36 state
transitions per individual in every 10,000 frames.
81
fraction of
colony active
A
1
0.8
0.6
0.4
0.2
0
00
coherence of oscillations
B
03
06
09
12
03
06
09
12
15
18
21
00
03
06
15
18
21
00
03
06
0.5
0.45
0.4
0.35
0.3
0.25
0.2
00
time (hours)
fraction colony active
Figure 4.3: Coherence of activity waves in ant colony 20. (A) Fraction of colony
active shown in blue, with the time of day shown on the x-axis. The time of lights
being turned on and off is indicated by the vertical red lines. (B) Coherence of
oscillations, with time of day on the x-axis.
1
0.8
0.6
0.4
0.2
0
21
00
03
time(hours)
Figure 4.4: Activity waves in ant colony 20, magnified view. Here we show the
collective oscillations from hour 20:00 to 05:00 on the second night of filming.
82
4.3
Results
All four colonies showed evidence of periodic behavior in activity cycles (Fig. 4.5),
with periods ranging from 10 to 20 minutes long (although we will focus on one colony
for the majority of our analysis here – colony 20, which displays the most coherent
oscillations). As in [Franks et al., 1990], we find no evidence for a refractory period,
i.e. a period of compulsory inactivity following a bout of activity. Refractory periods
have been posited as a mechanism for synchronization in some models of activity
cycles [Goss and Deneubourg, 1988], but is not necessary for synchronization with
MCA models as used in this work.
0.5
colony 33
colony 39
colony 40
coherence of oscillations
0.45
0.4
0.35
0.3
0.25
0.2
00
03
06
09
12
15
18
21
00
03
06
09
12
time (hours)
Figure 4.5: Coherence of activity waves in colonies 33, 39 and 40.
4.3.1
Role of interactions-model fitting
Since we classify ants into one of two possible states (active or inactive), it follows that
there are four possible types of interactions; (1) inactive-inactive (I-I), (2) inactiveactive (I-A), (3) active-inactive (A-I), and (4) active-active (A-A). Our notation here is
such that the first element represents the activity state of the focal ant, and the second
83
element represents the activity state of the interacting ant. We observe examples of all
four types of interactions in our data, although they occur with different frequencies
(Fig. 4.6). The I-I interactions are the least frequent, by virtue of the fact that inactive
ants are slow-moving and are unlikely to encounter other inactive ants. Although we
treat I-A and A-I interactions separately in our model-fitting, we only show one bar
for both types in Fig. 4.6, to avoid double-counting. We measure the effects of each of
the four types of interactions on the state transition likelihood. We call these effects
the coupling coefficients, and we denote them by CII , CIA , CAI , CAA .
18000
16000
frequency
14000
12000
10000
8000
6000
4000
2000
0
tiv
e-A
Ac
ve
cti
e
tiv
-
ve
ti
Ac
c
Ina
Ina
na
e-I
ve
cti
v
cti
Figure 4.6: Frequency of interaction types, in colony 20. Bars show number of activeactive, active-inactive, and inactive-inactive interactions
We expect that different types of interactions may have different effects on transition probabilities. To characterize the difference in responsiveness to different types
of interactions, we estimate the coupling coefficients by fitting a mobile cellular automaton model to our data. A similar method was used in [Cole and Cheshire, 1996],
although the amount of data in that experiment limited the possible analysis. This
model has similar qualities to an integrate and fire model of interacting neurons (as
in [Gerstner and Kistler, 2002]), as it combines effects of multiple stimuli (where here
the stimuli are interactions with other ants), and includes a parameter to control
84
memory of stimuli. We define the internal state of ant i, Si (t) (which reflects the
probability of changing state at time t) as a function of interactions experienced by
ant i, as follows:
Si (t) =
X
Cb tanh(g(δib + Si (t − 1)))
(4.1)
b
Where Si is the internal state of ant i, and higher values of Si indicate a higher
likelihood of transitioning between states. The sum runs over interaction types b,
where b can be any of the four interaction types described above, and Cb are the
coupling coefficients for each interaction type. δib is an indicator function that is 1 if
ant i has an interaction of type b, at time t, and 0 otherwise. g is a gain parameter,
and it controls the slope of the hyperbolic tangent function, and is a proxy for the
length of time an interaction is remembered. Higher values of g lead to a slower falloff
of Si , which means the interactions are remembered for longer, while lower values of
g lead to a faster decay in Si , which means the interactions are forgotten faster. Note
that only values of g less than or equal to 1 make biological sense, because g > 1 leads
to Si (t) that increases with time instead of decays with time. Using a tanh function is
convenient because it constrains Si to be between 0 and 1, so that many interactions
in rapid succession have diminishing effects on increasing Si (t). This model is slightly
modified from the one used in [Solé et al., 1993b] and [Miramontes et al., 1993], as
we here pulled the coupling coefficients, and sum, outside of the tangent function,
to allow for easier evaluation of the relative effects of different interaction types (see
Section 4.5.A for details). We understand that the exact behavioral rules used by ants
may have a different functional form from the one used here, but this model is useful
for comparing relative effects of interaction types on state transition probabilities.
85
To estimate the coupling coefficients, and to establish which types of interactions
have the largest effect in predicting when activity transitions occur, we consider separate state functions for each interaction type, by dividing Si into components. That
is, Si,II is the internal state of ant i due to only inactive-inactive interactions (it is
increased by an amount proportional to CII when ant i experiences an I-I interaction, and decays back to zero otherwise) and so on for the other interaction types.
Then, we fit a logistic regression model to the activity transition times, using each
Si,b as a predictor variable; so if Si,b is consistently higher when ant i transitions than
when ant i does not transition, interactions type b are predictive of when a transition
will occur. Thus, the resulting weights on the predictor variables correspond to the
relative importance of each interaction type on instigating a transition, and we use
them as estimates for the coupling coefficients in Equation 4.1 (note that these are
only estimates of the relative importance of interaction types, and are not exactly
substitute-able for the coefficients in Equation 4.1). Also note that we only consider
relevant types of interactions for each type of transition (we only use interaction types
where the focal ant is in the correct state for the transition type being analyzed). That
is, we use Si,II and Si,IA as predictor variables for activating transitions, and Si,AI
and SI,AA as predictor variables for deactivating transitions.
Si,II (t) = tanh(g(δi,II + Si,II (t − 1)))
(4.2)
Si,IA (t) = tanh(g(δi,IA + Si,IA (t − 1)))
(4.3)
Si,AI (t) = tanh(g(δi,AI + Si,AI (t − 1)))
(4.4)
Si,AA (t) = tanh(g(δi,AA + Si,AA (t − 1)))
(4.5)
Then, we use linear combinations of these separated components to fit the response
variable, R+/− , a binary variable that is equal to 1 if an ant is activated/deactivated
86
at time t, and 0 otherwise. Here data from all ants are combined in the response and
predictor variables, so the coefficients we fit represent the colony average. We hope
to extend the analysis to fit individual behavioral rules.
+/−
Where C0
R+ = logit−1 (C0+ + CII SII + CIA SIA )
(4.6)
R− = logit−1 (C0− + CAI SAI + CAA SAA )
(4.7)
are the baseline transition probabilities given no interactions. We
selected the gain parameter, g, by optimizing over 1000 randomly selected values;
choosing the value which resulted in the highest log-likelihood. The optimal value for
g is 1 (see Section 4.5.B).
Of the four interaction types, we find that the inactive-active (I-A) interactions
are overwhelmingly the most predictive of when activating transitions occur. That
is, ants are most responsive to those interactions they have with active ants, when
they themselves are inactive. For deactivations, although the effect is smaller, both
active-inactive (A-I) and active-active (A-A) interactions have inhibitory effects. In
other words, an active ant that interacts with another ant, in either state, is more
likely to remain active than an ant which doesn’t interact at all. We illustrate these
effects in Fig. 4.6, in which we plot the relative transition probability, where we
compare the transition likelihood resulting from each type of interaction, with the
baseline transition probability of non-interacting ants. More specifically, the relative
transition probability in Fig. 4.6 is the predicted response given the best-fit coupling
coefficients, and an internal state Sb = x, divided by the model intercept, with Sb = 0.
For example, to calculate the relative transition probability for I-A interactions, for
+
all SIA = x and for a given constant value of SII , PIA,relative
= (C0+ + CII ∗ SII + CIA ∗
x)/(C0+ + CII ∗ SII ). An inactive ant which has interacted with other active ants such
87
that it has a value of SIA = x will be y times more likely to transition in any given
frame than an ant which has not experienced any interactions with active ants.
B
25
25
Activations
relative transtion probability
Activations
relative transition probability
A
20
15
10
5
0
0
0.2
0.4
0.6
0.8
20
15
10
5
0
1
0
0.2
0.4
SII
D
1.2
Deactivations
relative transition probability
Deactivations
relative transition probability
C
1
0.8
0.6
0.4
0.2
0
0
0.2
0.4
0.6
0.8
1
0.6
0.8
1
SIA
0.6
0.8
1
SAI
1.2
1
0.8
0.6
0.4
0.2
0
0
0.2
0.4
SAA
Figure 4.7: Relative transition probabilities, with activation probabilities in (A)
and(B), and deactivations in (C) and (D). The relative transition probabilities compare probability of changing state at S = x, compared to S = 0.
4.3.2
Responsiveness over time
One possible explanation for why the ant colonies are at times very synchronized,
and at others not synchronized at all in their activity, is that there exists a flexibility
in behavioral response to interactions (e.g. the coupling coefficients we fit in the
preceding section are not constant, but rather change over time). This variability
could be attributed to preoccupation with other tasks, such as feeding, or caring for
the brood. If this is the case, then coherence of oscillations should decrease when
there is a larger rate of influx into the nest (foragers returning, carrying food, for
88
example). We do in fact find a strong negative relationship between coherence of
oscillations and rate of influx to the nest (see Fig. 4.8). This finding is in agreement
with [Boi et al., 1999], where they find more pronounced oscillations when foragers
are prevented from returning to the nest.
0.8
R9=9−0.44967,9P9=92.6794e−12
0.7
coherence
0.6
0.5
0.4
0.3
0.2
0.1
2
3
4
5
inflow9rate9(ants/frame)
6
−3
x910
Figure 4.8: Colony inflow rate (colony 20). Coherence of oscillations decreases as
more ants flow into the nest
We can measure how the behavioral response to interactions changes over time,
and compare this responsiveness to how coherent activity waves are, at different
points throughout the activity time-series. We assess the variation in responsiveness
by measuring the coupling coefficients from fitting the MCA model, as above, but
instead of using the entire time-series, we fit the model to smaller windows in the
data; at every 1,000 frames, we use samples of 10,000 frames (measured at the same
values as for oscillation coherence, as described in Section 4.2). Thus we can measure
how the coupling coefficients vary over time.
We find that there is considerable variability in CIA over time (the responsiveness
to I-A interactions,) and that CIA strongly correlates with the coherence of activity
waves (Fig. 4.9). Although there is some variation in values of the other coupling
coefficients (CII , CAI , CAA ), all values are quite close to zero, indicating that they are
small effects, and none correlate strongly with coherence (Fig. 4.9A,C,D).
89
A
B
R = −0.032775, P = 0.62954
R = 0.5941, P = 2.766e−22
0.5
coherence
coherence
0.5
0.4
0.3
0.2
0.1
0
−0.1
0.4
0.3
0.2
0.1
0
0.1
C
0.2
CII
0.3
0
−0.1
0.4
0.2
CIA
0.3
0.4
0.6
R = 0.17357, P = 0.01007
R = 0.088343, P = 0.19277
0.5
coherence
0.5
coherence
0.1
D
0.6
0.4
0.3
0.2
0.1
0
−0.1
0
0.4
0.3
0.2
0.1
0
0.1
0.2
CAI
0.3
0
−0.1
0.4
0
0.1
0.2
CAA
0.3
0.4
Figure 4.9: Coherence and interaction coefficients. Increases in responsiveness to I-A
interactions correspond to increases in coherence of activity waves.
While we observe that CIA is strongly correlated with coherence in activity waves,
we want to check it against alternative explanations for this variation in coherence,
such as density or spatial effects. We assemble a list of possible effects, including all
four coupling coefficients, the rate of interactions, the average group activity level,
the average distance to near neighbors, and the average distance to the colony center.
We examine the relative importance of these features in predicting coherence in two
ways, similar to the analysis of effects in predicting behavioral cascade magnitude, in
Chapter 3.
First, we perform an L1 penalized regularization on the set of features
[Park and Hastie, 2007].
This method is useful because it penalizes more com90
plicated models, and drives unimportant features to zero. Standard practice is to
measure predictor strength at the point of minimum deviance, plus one standard
error. At this location, all features except CIA are driven to zero, leaving CIA as a
positive predictor of oscillation coherence (Fig. 4.10A). The strength of the predictor
can be thought of in the following way. If a predictor variable has strength x in a
model, a unit increase in that predictor will result in an increase in x in the response
variable. (Note that we use a standardized coherence here, so that the maximum
value is 1, and standardized predictor variables, so that the resulting predictor
strengths are comparable to one another).
Secondly, we analyze the unique contribution to the variance in response variable
that is explained by the full model (where the full model is the linear combination
of all proposed features with coefficients fit to the response variable, using a linear
link function). The full model explains 55% of the variance in the response variable
here. Each bar in Fig. 4.10B displays a feature’s unique contribution to this explained
variance. The unique contribution is calculated by subtracting the variance explained
by a model which leaves out one feature, from the variance explained by the full
model. The feature whose absence causes the biggest decrease in variance is the one
that explains the largest portion of it. Note that the sum of the bars in Fig. 4.10B
is less than 0.55, because some portion of the explained variance is not unique to a
single predictor variable. Using this method, we again find that CIA is the largest
contributor to the variance explained by the model.
91
B
minimum2deviance2+21SE
L1 regularization
A
cAA
cAA
cAI
cAI
cIA
cIA
cII
cII
int2rate
int2rate
activity
activity
dist2center
dist2center
density
density
0
0.05
0.1
0
predictor2strength
0.1
0.2
0.3
unique2contribution2to2explained2variance
Figure 4.10: Relative contribution to variation in coherence. (A) relative feature
weight from L1 regularization. CIA is the strongest effect. (B) Each variable’s unique
contribution to the total variance explained by the model. Again, CIA is the strongest
effect.
In future work we will compare the performance of the MCA model to other
candidate models, such as the integrate and fire model, in which memory is allowed
to be infinite, in fitting the data. We will also consider the feasibility of taking an
individual-based, instead of a colony-averaged approach. That is, we will evaluate
the responsiveness of individual ants to the four different types of interactions, to
investigate whether the collective oscillations are driven by a subset of very responsive
ants, or whether all individuals in the colony participate equally.
4.3.3
Simulating collective oscillations
Having observed the relationship between responsiveness to interactions and oscillation coherence in the data, we explore it further by conducting simulation experiments on mobile cellular automaton models, similar to those in [Solé et al., 1993b,
Miramontes et al., 1993]. In the MCA simulations, agents change state based on ‘interactions’ (mediated by proximity) with neighboring agents, updating their internal
state according to Equation 4.1. An agent becomes active when interactions with
neighbors result in its internal potential crossing a minimum activation threshold.
92
Likewise, the agent becomes inactive when its internal potential decreases below the
activation threshold. Unlike in the model-fitting described above, in the simulations,
there exists the need for a movement rule, governing how agents behave when in
the active or inactive states. In previous models, agents moved randomly on a lattice when active (to a randomly selected, unoccupied, space on the lattice one step
away from its current location), and remained stationary when inactive. While this
random movement was sufficient to produce global oscillations in simulations, this
description of individual movement does not accurately reflect the behavior of ants
within the nest, which exhibit a tendency to occupy ‘home ranges’ or ‘spatial fidelity
63 mm
zones’ [Sendova-Franks and Franks, 1995] (Fig. 4.11).
42 mm
Figure 4.11: Spatial fidelity of 9 randomly selected ants from colony 20. The heatmaps
represent the fraction of each ant’s active time they spent at each position.
To incorporate this homing tendency, we implement a more realistic movement
rule here, where agents are attracted to the central brood with a homing parameter,
93
pH , such that ant i has probability pH of attempting to move in the direction of the
colony center at every time-step. In the event that the desired space is occupied by
another ant, the focal ant resumes random movement. If there are no unoccupied
spaces one step away, the agent remains stationary for that time step. Changing pH
controls the density of the colony.
0.9
0.8
fraction simulated
colony active
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
200
400
600
800
1000
timesteps
Figure 4.12: Activity waves in a simulated colony, using optimal interaction coupling
coefficients.
We perform simulation experiments on groups of 100 agents, with random initial
positions on a 50 unit square lattice. We explore the parameter space of the model,
first finding the combination of parameters that results in the largest oscillation coherence, by comparing 1000 different sets of coupling coefficients. In Fig. 4.12 we plot
the activity cycles in a simulated colony using these optimal parameters. Note that
the ants begin the simulation dispersed randomly on the lattice, and the coherence
of oscillations in activity increases as they become more clustered around the center.
Next, starting from this optimal set, we vary each of these parameters in turn, and observe how sensitive the resulting oscillation coherence is to each coefficient (Fig. 4.13).
94
The optimal values for CII and CAI fall close to zero, such that higher values of these
parameters lead to less coherent oscillations. The optimal value for CAA falls at a
small, but nonzero number (about 0.15), with a steep drop off in coherence beyond
0.25. Only in CIA do we observe a strong increase in coherence with the value of the
coefficient, up to about 0.7, at which point it levels off. This is consistent with what
we observe in the data- that coherence of oscillations in activity level is controlled by
the strength of responsiveness to inactive-active interactions. In the data we don’t
observe the entire range of coupling coefficients, so we can’t measure the relationship
between coherence and CIA for high values of CIA , to see if we would eventually
observe the same leveling off.
Here we have selected the other parameters in the simulations such that they allow
for periodic activity waves. In the future we would like to more carefully match these
parameters to our observations of the real ant colonies, as well as to analyze how
sensitive our results are to parameter values. In general, a more thorough exploration
of this model’s parameter space is necessary to fully understand the driving factors
of collective oscillations.
95
A
B
0.25
coherence (simulation)
coherence (simulation)
0.25
0.2
0.15
0.1
0.05
0
0
0.2
0.4
0.6
0.8
C
0.05
0
0.2
0.4
0.6
0.8
1
0.8
1
CIA
0.25
coherence (simulation)
coherence (simulation)
0.1
D
0.25
0.2
0.15
0.1
0.05
0
0.15
0
1
CII
0.2
0
0.2
0.4
0.6
0.8
1
0.2
0.15
0.1
0.05
0
0
0.2
0.4
0.6
CAA
CAI
Figure 4.13: Coherence and coupling coefficients, simulations. Increasing CIA (responsiveness to inactive-active interactions) increases coherence, which is consistent
with what we measure in the data.
4.4
Conclusions
This study employed a high-resolution dataset of individual behavior to address variation in coherence in ant colony activity cycles. We use individual level data to
estimate the behavioral responses to interactions of different types (Fig. 4.7). We
further demonstrate the novel result that the behavioral response to interactions,
specifically those between inactive and active ants, varies across time and correlates
strongly with the coherence of oscillations in collective activity (Fig. 4.9). We show
that this variation in responsiveness explains the variation in coherence better than
other possible factors, including the rate of interactions, the average group activity
level, average distance to near neighbors, and average distance to colony center. In
96
addition, we demonstrate that coherence of oscillations in a simulated colony can be
driven by this same coupling coefficient, using a spatially explicit MCA model.
While it has been suggested that activity cycles are merely an epiphenomenon of
communication and the capacity for positive feedback in ant colonies [Cole, 1991b],
a potential benefit to synchronization could arise if colonies are more efficient when
ants are active simultaneously [Robinson, 1992]. Another possibility is that by synchronizing periods of activity, worker ants within the nest are more evenly spread,
resulting in greater task efficiency [Hatcher et al., 1992, Delgado and Sole, 2000]. It
is difficult to design experiments particularly with colonies in natural habitats to
assess fitness advantages, but investigating the effects of potential regulators of coherence may provide more insight into the adaptive function of activity cycles in ant
colonies, if such a function exists.
In this study, we have shown using large volumes of individual data from real
colonies how interactions drive synchronized and unsynchronized activity inside the
nest. We have provided evidence for behavioral flexibility in individual ants in response to interactions, and demonstrated the role of this flexibility in moderating
the coherence of oscillations in collective activity. Our study is the first to examine
variation in individual ant responsiveness to interactions, and our results underline
the importance of understanding subtleties in individual behavior when modeling the
emergence of synchronization in biological systems.
4.5
4.5.A
Technical details
Modifications to the MCA model
The original MCA model, used in [Miramontes et al., 1993, Solé et al., 1993a] has the
following form:
97
Si (t) = tanh(g(
X
Cb Sb (t − 1) + JSi (t − 1)))
(S4.8)
b
We have made a few modifications to this model for use in our simulation experiments (the form of the model we use can be found in Equation 4.1). First, instead
of summing over the internal states of neighboring ants, we sum over their external
states (either active or inactive), replacing the internal state with an indicator function which is 1 if an interaction of type b happened, and 0 otherwise. This is a more
plausible scenario, biologically, as it is likely that communication is on a coarse scale
(see Fig. S4.1 where we show that the speed of interacting ant has negligible effect
on the speed of the focal ant).
A
B
0.4
0.3
0.7
0.6
Inactive−Active
∆ focalaspeed
Active−Active
∆ focalaspeed
0.2
0.1
0
−0.1
−0.2
0.5
0.4
0.3
0.2
−0.3
0.1
−0.4
−0.5
0
2
4
6
neighboraspeed (mm/s)
8
0
0
2
4
6
neighboraspeed (mm/s)
8
Figure S4.1: Speeds of interacting ants have negligible effect on the change in speed
of the focal ant, in colony 20. Here we show the average change in speed, 20 seconds
after an interaction occurred, for (A) active-active interactions and (B) inactiveactive interactions. The remaining two interaction types are not shown here because
the speed of inactive ants is negligible.
Secondly, we have removed the self-interaction parameter, setting J = 1, as we
found it added unnecessary complexity for our model-fitting purposes. Third, we
extracted the coupling coefficients from the argument of the tanh function, in order
98
to make comparison easier with the results from the logistic regression analysis. In
the remainder of this section, we will show that by writing the model in this form, it
functions in the same way as when the coupling coefficients remain inside the tanh
argument.
First, consider the model at a time, t, when no interactions are currently happening. The model, Equation 4.1, is then reduced to:
Si (t) = (CII + CIA + CAI + CAA ) tanh(gSi (t − 1))
(S4.9)
Here δi,b = 0, because no interactions are present. Thus when there are no interactions, Si decreases according to the value of g. The sum of the coupling coefficients
just controls the magnitude of Si , but does not affect the rate of decrease.
Secondly, consider a case where there is a single interaction, at time t. We select
an I-I interaction, as an example, but can be easily generalized to any interaction
type. In this case the model becomes:
Si (t) = CII tanh(g(1 + Si (t − 1)) + (CIA + CAI + CAA ) tanh(gSi (t − 1))
(S4.10)
We can rewrite this as:
Si (t) =CII (tanh(g(1 + Si (t − 1)) − tanh(gSi (t − 1)))+
(CII + CIA + CAI + CAA ) tanh(gSi (t − 1))
(S4.11)
(S4.12)
The last term in the equation is the same as in the case of no interactions- this
causes the decrease in potential that happens at every step. The first term, multiplied
99
by CII represents the amount added to Si when there is an interaction of type I-I. So
the larger CII is, the more responsive are ants to interactions of type I-I.
Third, we consider the case where an ant experiences multiple simultaneous interactions. We generalize this fully, to any number of interactions of any type, although
in reality simultaneous interactions are rare.
Si (t) =CII (tanh(g(nII + Si (t − 1)) − tanh(gSi (t − 1)))+
(S4.13)
CIA (tanh(g(nIA + Si (t − 1)) − tanh(gSi (t − 1)))+
(S4.14)
CAI (tanh(g(nAI + Si (t − 1)) − tanh(gSi (t − 1)))+
(S4.15)
CAA (tanh(g(nAA + Si (t − 1)) − tanh(gSi (t − 1)))+
(S4.16)
(CII + CIA + CAI + CAA ) tanh(gSi (t − 1))
(S4.17)
Where nII , nIA , nAI , nAA are the respective number of interactions of each type
occuring at time t.
From this algebra you can see that the main effects of the model are preserved
when the coupling coefficients are extracted from the argument of the tanh function.
That is, there is a constant decrease in Si at each time step, regardless of the number
of interactions occurring. When interactions do occur, they simply add to Si (t), an
amount proportional to the value of the coupling coefficient, and dependent on the
value of Si (t − 1).
4.5.B
Best-fit value of g
In our model-fitting procedure, before fitting the interaction coefficients, we first found
the value of the gain parameter (g) which resulted in the best model performance.
Out of 1000 randomly selected values of g, we find that the model performance is
maximized at the highest possible value, g = 1, (Fig. S4.2). This suggests that ants
100
remember prior interactions, to some extent. In the future we plan to compare the
MCA model performance against alternative models, particularly those which allow
for infinite memory in the system, so we can more accurately establish the duration
for which memory of an interaction persists.
B
A
−8000
−8190
−8010
−8192
−8194
−8030
Deactivations
Log−likelihood
Activations
Log−likelihood
−8020
−8040
−8050
−8060
−8070
−8198
−8200
−8202
−8204
−8206
−8080
−8090
−8196
−8208
0.2
0.4
0.6
0.8
−8210
1
g
0.2
0.4
0.6
0.8
1
g
Figure S4.2: Best-fit gain parameter in colony 20 for (A) activating transitions, and
(B) deactivating transitions. For both transition types the best-fit gain parameter is
1.
4.5.C
Colonies without brood
In addition to the four normal colonies used in analysis in the main text of this
chapter, we have data from two colonies from which the brood was removed. These
brood-less colonies display similar patterns of activity waves, although they are less
pronounced than in the colonies with brood present. Another characteristic of ants
in these non-brood colonies is that they tend to display less spatial fidelity (they tend
to be less clustered around the nest center than ants in colonies with brood). We
have observed in simulations that having a higher homing parameter leads to more
coherent colony oscillations. This is consistent with what we find in the data, with
brood colonies having both a higher coherence and a higher ‘homing’ than non-brood
101
colonies (Fig.S4.3, where here we define homing as the variance in distance from
colony center).
variance in position
0.9
0.8
0.7
0.6
0.5
0.4
Colony 20NB
no brood
0.3
Colony 20
brood
0.2
0.1
0.284785
0.37811
<coherence>
Figure S4.3: Homing and coherence in two colonies- one with brood and one with
no brood. Here we use the variance in distance from colony center as a proxy for
homing. We find that ants in the non-brood colony tend to have higher variance in
position than ants in the brood colony, and that the non-brood colony has a lower
average coherence than the brood colony.
102
Chapter 5
Exploration and exploitation in
dynamic environments
In this chapter, we investigate collective behavior in a different ecological context,
where simulated individuals in a group use collective computation to optimally benefit
from environmental resources (such as tracking moving food patches, or sheltered
regions). They must balance a desire to stay in a group, to take advantage of the many
benefits associated with a large group of conspecifics, with a desire to find and exploit
the best resources in the environment, and thus avoid overcrowding. In this chapter
we study this balance of exploration with exploitation of environmental resources in
a social, experimentally motivated, self-propelled particle model of collective motion,
in which agents are coupled to their environment.
The work in this chapter is the result of a collaboration with Andrew M. Hein
and George Hagstrom, along with co-authors Andrew Berdahl, Colin J. Torney, and
corresponding author Iain D. Couzin. This chapter contains a lightly edited version
of the main manuscript and supplemental information currently in prep. AMH, SBR,
GH, CJT, and IDC developed the concept for the study. AMH, SBR, and GH devel-
103
oped and performed simulations and conducted mathematical and statistical analyses.
AMH, SBR, GH and IDC wrote the manuscript.
5.1
Introduction
Many of the benefits of grouping arise from enhanced access to information. Fish in
shoals perceive predator attacks that occur far away in the shoal [Handegard et al., 2012],
foraging bats localize prey by eavesdropping on the echolocation calls of their group
mates [Cvikel et al., 2015], and human problem solvers discover favorable regions of
solution space by observing and mimicking other members of their social networks
[Mason and Watts, 2012].
These and similar observations have led to the view
that animal groups can act as a distributed sensors, affording each group member
access to information at a larger spatial scale than it could perceive were it alone
[Berdahl et al., 2013, Torney et al., 2009].
Because groups that comprise more individuals often occupy more space, there
is a strong link between group size and the enhancement in performance that distributed sensing affords. However, from the perspective of an individual animal,
there are limits to the benefits of being in a large group. Individual performance
in a variety of tasks improves as group size increases, up to a point, and then begins to decline for larger groups (e.g. [Cvikel et al., 2015, Kao and Couzin, 2014]).
For example, if groups are too large relative to the spatial extent of environmental
gradients, the gradient tracking ability of individuals decreases [Berdahl et al., 2013,
Torney et al., 2009]. Strong resource competition and other effects of crowding compound these costs in large groups. From the perspective of individuals, then, there
ought to be a particular group size that maximizes performance [Sibly, 1983]. In the
scenario of resource gathering, for example, individuals must continuously evaluate
their current state, to decide if it would be better to continue exploration of the envi104
ronment to possibly find a more attractive resource, or to commit to exploitation of
the resource in their current location. In addition, social animals must make many
judgment calls on whether to leave the benefits and relative safety of a group in order
to have better access to resources.
How can an individual ensure that it is in a group of optimal size? The ubiquitous
fission-fusion dynamics displayed by groups of social animals [Couzin and Krause, 2003]
are de facto evidence that actions of group members can cause large and sometimes
frequent changes in group size. Yet, in animal groups such as pelagic fish schools,
bird flocks, and insect swarms, individuals likely have access only to local information
about group mates and it is unclear how an individual could evaluate the size of the
group to which it belongs, let alone choose to leave or join a group to achieve an
optimal group size. Moreover, classical evolutionary models of group size regulation
produce predictions that seem spurious; for example, that groups should increase in
size until the fitness of social and solitary individuals is equal, and that all groups in
a given environment should be the same size [Sibly, 1983]. The discrepancy between
predictions of evolutionary models and data on group size has become known as the
“group size paradox” [Beauchamp and Fernández-Juricic, 2005].
Here we present a new approach to studying group size regulation, that explicitly considers the rules individuals use to make movement decisions. We develop a
simple model that captures two fundamental constraints on the grouping behavior of
social animals: the limited ability of individual animals to perceive distance regions of
their environment (i.e., limited perception), and the spatial nature of the movement
decisions that ultimately determine group membership. In the model, individuals
respond to the locations of their near neighbors and also to local measurements of a
dynamic environmental cue, which is in this case the level of a resource that individuals seek to maximize. While our work is motivated by detailed studies of schooling
fish [Berdahl et al., 2013, Katz et al., 2011], the structure of the model is generic and
105
could be applied to social foraging more generally, or to autonomous vehicles solving
an environmental monitoring problem (e.g. [Bachmayer and Leonard, 2002]).
We demonstrate that simple decision rules based only on local information allow individuals to balance exploration and exploitation [Sutton and Barto, 1998] to
maximize a reward. Individuals accomplish this balance by dynamically breaking
away from large groups when their instantaneous level of reward is low, and fusing to form larger supergroups when their level of reward is high. Fission-fusion
dynamics, which are widely exhibited by animal groups [Couzin and Krause, 2003,
Couzin and Laidre, 2009], arise naturally from two simple behavioral rules. By contrast to past studies of group size regulation, we show that the decision rules that
maximize individual performance confer the flexibility to form groups of many different sizes depending on the environmental context, rather than leading to a single
equilibrium group size [Slobodchikoff, 2013]. Our results suggest that individual performance in dynamic decision-making problems can be greatly improved through
social interactions even though individuals may lack the ability to actively regulate
group size.
5.2
Model development
We model a population of animals interacting with one another and with a sparse,
non-consumable resource (we discuss the case of a consumable resource in Section
5.5.D; results do not differ quantitatively from those presented below). To capture
the spatial nature of movement decisions, we model individuals, each of which moves
according to an environmental response rule and a social response rule. Both rules
are derived from experimental data on the movement decisions of freely swimming
golden shiners (Notemigonus crysoleucas) [Katz et al., 2011, Berdahl et al., 2013].
Thus, an important distinction between the model developed here and many self106
propelled particle models used to study collective behavior (e.g. [Vicsek et al., 1995,
Couzin et al., 2002]) is that we begin with individual behavioral rules based on experimental data rather than choosing interaction rules that produce realistic-looking
collective behavior.
To model the response of individuals to the environment, we develop an environmental response rule based on responses of golden shiners to environmental cues.
Berdahl et al. measured movements of individual fish and individuals in schools in a
dynamic light environment; they measured how each fish responded both to neighbors
and to the properties of the light field in its immediate vicinity [Berdahl et al., 2013].
Each fish responds strongly to absolute level of light at its position by swimming
faster when the light level is high and slower when the light level is low. Fish respond
only weakly to the spatial gradient in light levels and appear to adjust their headings
based on social cues rather than features of the light field. However, at the level of an
entire school, differences in speed across the school causes the school to turn toward
darker regions of the environment. The school exhibits emergent taxis behavior and
descend light gradients even though each individual fish responds primarily to the
absolute light level rather than the local light gradient. To capture these observed
behaviors, we define the ith individual’s environmental response as a function of the
level of an environmental cue at its current position:
m
vi
dva,i
= (ψ0 − ψ1 S(xi ) − η|vi |2 )
,
dt
|vi |
(5.1)
Where m is the individual mass (taken to be 1 hereafter), xi and vi are the position
and velocity of the ith individual respectively, and η is a damping term that limits
individuals to a finite speed. ψ0 is a constant that determines the default agent speed
in the absence of social or environmental interactions (and functions as a proxy for
107
how explorative each individual is), ψ1 is a ‘gain’ term that determines the response
of the ith individual to the resource value S(xi ) at its position. This resource is
dynamic, and moves in a directed drift motion (defined below, in 5.5). Individuals
simply change their speeds in response to the resource; they do not change their
headings. This type of simple environmental response provides a minimal mechanism
for response to environmental cues that can be studied as a benchmark for comparison
with other more complicated methods, such as gradient sensing. In what follows, we
refer to “cue” and “resource” interchangeably as we model the case where the cue
is the resource itself (see for example, [Hein and McKinley, 2012, Torney et al., 2009]
for cases where the cue is not a resource).
In addition to responding to the environment, individuals respond to a set of near
neighbors according to a social response rule, based on observed pairwise interactions
among freely swimming golden shiners [Katz et al., 2011]. Katz et al. used a force
matching method to infer the effective “social forces” exerted on a focal fish by a
neighboring fish. The movement behavior of the focal fish is represented by a force
exerted by its neighbor that can accelerate, decelerate, and turn the focal fish. This
approach allows one to compute the effect of each individual on the movement decisions of the other individual as a function of the relative positions of the two. Katz
et al. found that a near neighbor exerts a repulsive force on the focal individual,
causing the focal individual to accelerate if the neighbor is behind it or decelerate if
the neighbor is in front. As the distance between the focal individual and neighbor
increases, the force exerted on the focal individual gradually changes from repulsive to
attractive, with maximum attraction occurring at a distance of two-four body lengths.
For longer distances, the force exerted on the focal individual is still attractive but
decreases in magnitude. Critically, there was no evidence that individuals explicitly
match the body orientation of their neighbors, as has often been assumed in past
models (e.g., [Couzin et al., 2002]). Another study [Herbert-Read et al., 2011] iden108
tified similar interaction rules in mosquitofish (Gambusia holbrooki) and also found
little evidence that individuals match body orientations of neighbors.
To capture the features of observed social interactions among fish, we model the
response of the individual i to its neighbors by the gradient of a social potential
(where the social potential is a Morse function).
#
"
X
dvs,i
m
= −∇
Cr e−|xi −xj |/lr − Ca e−|xi −xj |/la ,
dt
j∈N
(5.2)
i
where ∇ is the two-dimensional spatial gradient, the term in brackets is a social
potential, Ca , Cr , la and lr are constants, and the set Ni is a set of the k nearest
neighbors of the ith individual, where a neighbor is an individual within a distance
lmax of the focal individual. Setting a finite number of near neighbors and finite
detection distance ensures that individuals can only respond to a limited number of
their neighbors in crowded regions of space, and provides a simplification of sensorybased models of social interactions [Strandburg-Peshkin et al., 2013].
The social interaction rules specified by 5.2 are consistent with observed social
interactions of a wide variety of social animals: at short distances individuals repel one
another such that an individual will speed up if it is directly in front of a neighbor and
slow down if it is directly behind a neighbor. At intermediate distances, individuals
tend to adjust their headings to move towards one another (although there is no
explicit alignment in the model). At long distances, individuals cannot detect one
another [Katz et al., 2011].
Combining the social and environmental response rules yields two equations that
govern each individual’s movement:
dxi
= vi
dt
109
(5.3)
m
dva,i dvs,i
dvi
=
+
.
dt
dt
dt
(5.4)
The full parameter space in a similar model is explored in [D’Orsogna et al., 2006],
with ψ1 = 0 (that is, the agents do not interact with the environment); here we
focus on a parameter regime that gives short-range inter-individual repulsion and
intermediate-range attraction.
To study the behavior of populations of individuals, we simulate a discretized
version of the system described by Equations 5.3 and 5.4. In particular, we choose
a time step, τ , within which the acceleration due to social influences (Equation 5.2)
and resource value S are assumed to be constant. Positions, speeds, and accelerations
of all individuals at time t + τ are then given by the solutions to Equations 5.3 and
5.4 at time t + τ with the values of S(xi ) and |xi − xj | determined at time t. A small
navigational noise vector of magnitude ω and uniform heading (0 to 2π) is added
to the velocity of each agent at each time step. Taking the limit as τ goes to zero
means that individuals are constantly acquiring information from the environment
and instantaneously altering their actions in response. A deterministic version of this
yields the fully continuous system described by Equations 5.3 and 5.4, which can be
analyzed alongside discrete simulations (Section 5.5.F).
To determine how populations of animals interact with a resource environment,
we studied populations of N individuals in a two-dimensional environment containing one or more resource peaks. We computed each individual’s performance as
the mean value of the resource it experienced over a foraging bout. By evaluating the performance of individuals as a function of the parameters in Equations 5.1
and 5.2, we could determine how the behaviors of individuals affect their performance. We could also study the sizes of groups that form, splitting and fusing of
110
these groups, and the relationship between the locations of groups and the locations
of resources without invoking a priori assumptions about the relationship between
group size and performance as is commonly done in studies of group size regulation
(e.g. [Beauchamp and Fernández-Juricic, 2005, Sibly, 1983]. For reference, we compare the performance of social individuals to the performance of asocial individuals
(i.e. Ca and Cr = 0 in Equation 5.2).
The social interaction rule allows us to build an interaction network for the entire
population. Two individuals are connected if at least one of them influences the other
through Equation 5.2. We define a “group” as a set of individuals that belong to the
same connected component in this network.
Resource environment
We model individuals in a large two-dimensional environment. We assume the boundary of environment is periodic such that individuals, inter-individual potentials, and
the environmental signal (i.e. the resource) are all projected onto a torus. The resource is distributed over the environment as a set of one or more resource peaks, each
of which decays exponentially in value with increasing distance to the peak center.
The value of the resource in a single peak at a location Xi , is given by
S(x, xs ) = λ0 exp−1|x−xs |/λ1
(5.5)
Where λ0 is a constant that determines the resource value at the peak center,
λ1 is a decay length parameter, and xs is the location of the centroid of the peak
of interest. The total resource value an individual experiences is the sum over all
peaks in the environment, however, we are interested in cases where resources are
sparse (i.e. the environment varies spatially in quality) and individuals will typically
111
only experience an appreciable resource value from a single peak. Each peak moves
stochastically, according to Brownian motion, but with a directed drift. In the cases
of multiple peaks, all peaks drift in parallel motion, to avoid peaks becoming near
enough to cause interference between groups. The decay length, λ1 is much smaller
than the length of the environment such that it is reasonable to calculate the effective
size of a peak as the total peak mass, 2πλ0 λ21 .
There exists a critical angle for an agent to be captured by a resource peak of
given size and given occupancy (number of agents currently occupying the peak).
This behavior is reminiscent an object in orbit around a planetary mass, where in
this case the mass is increased if the resource peak has a higher occupancy. We derive
this critical angle (∆c ) in Section 5.5.F, and include here a vector field for illustration,
with peak occupancy N = 8 (Fig. 5.1), with the red lines representing the trajectories
of particles approaching the peak on the edge of the capture region.
Figure 5.1: Critical angle for capture by peak. The green region corresponds to
zero velocity, and the light yellow interior region corresponds to the region where the
potential is stronger than centrifugal forces, which for pi/2 < ∆ < π/2 represents the
trapping region. The two red circles correspond to the unstable equilibrium points,
and the red trajectories represent the boundaries of the capture region.
112
Simulations
We performed two types of numerical experiments to study the search performance
of individuals. In experiment 1, we allowed individuals to begin the simulation in a
single social group, centered at the location of a signal peak. This environment only
contains one signal peak, so individuals that leave the region containing the peak will
experience a low signal value. We then recorded the signal value experienced by each
individual through time, as the peak moved stochastically through the environment.
The goal of this experiment was to determine whether individuals that act according
to the rules described above are capable of tracking a dynamic signal through space. In
experiment 2, we again started individuals in a single group, centered at the location
of a signal peak. In this experiment, however, the environment contained a second
signal peak, located far from the first. Peaks move according to the same stochastic,
directed drift rules, and move in parallel, to ensure that the distance between them
remains large. In this second experiment, individuals that leave the first peak can
benefit if they are able to locate and track the second peak. We discuss more complex
environments in Section 5.5, but qualitative results are unchanged.
In both experiments, we fixed some parameters at constant values such that the
individuals experience social attraction and repulsion that is consistent with the rules
observed in real animal groups (e.g. short range repulsion and long range attraction,
see [Katz et al., 2011]). See [D’Orsogna et al., 2006] for a full exploration of parameter space using a similar model. Unless otherwise noted, all simulations used the
parameter values or ranges listed in Table 5.6. Each simulation was run for 10,000
time steps. We measured the performance of each individual, in each simulation as
the mean value of the resource (< Si >), experienced over the second half of the simulation, where we exclude the first half to avoid transient effects of initial conditions.
Individuals that spend all their time far from peaks have a value of < Si > near zero,
113
whereas individuals that spend their time near a peak will have a value of < Si >
near λ0 , the maximum value of the resource peak.
5.3
5.3.1
Results
Social interactions help individuals exploit resource
peaks
Individuals that responded both to neighbors and to the environment dramatically
out-performed asocial individuals (Fig. 5.2A). Asocials slow down when the mean
value of the resource is high, which allows them to perform well over a short period
of time (Fig. 5.2B). However, as the resource peak moves, asocials cannot track it
and the number of individuals near the peak center quickly declines (Fig. 5.2B). In
contrast to asocials, individuals that respond to the environment and social cues
effectively track the peak as it moves through the environment (Fig. 5.2C,D) and
aggregate more quickly near the peak even in a static environment. Critically, social
individuals are able to remain near a moving resource peak without any knowledge
of how the peak is moving, how large the peak is, or where they are in relation to
the peak’s center. Mean resource value of asocials falls rapidly as the drift speed of
moving peaks increases (Fig. 5.2B). Social agents, on the other hand, experience a
much slower decline in performance over a broad range of peak speeds.
114
mean=signal=value
A
C
1.4
1.2
1
0.8
mean=signal=value
t===116
t===212
t===308
t===404
t===500=
t===20
t===116
t===212
t===308
t===404
t===500=
0.4
0.2
0
B
t===20
0.6
social,
sensing
social,
asocial,
sensing non-sensing
3
D
social
asocial
2.5
2
1.5
1
0.5
0
0
0.05
0.1
0.15
peak=drift=speed
0.2
Figure 5.2: Performance and resource tracking by social agents. (A) The mean resource value (mean taken over 1000 replicate simulations, of 300 individuals and the
second 2500 steps of the total 5000 steps) of individual searchers. Error bars are
the standard error over 1000 simulations. (B) Mean resource value of social (orange)
and asocial (blue) agents as a function of peak drift speed (see Methods, 5 peaks
in environment). Note rapid decrease in performance of asocial agents as peak drift
speed increases. Each point represents results from a single simulation. (C) Sequence
of asocial individuals (blue points) interacting with exponentially decaying resource
peak (resource density in grayscale). View zoomed in to region surrounding peak.
Each individual’s trajectory is represented by the tail attached to each point; the
length of tail is proportional to speed. Individuals decelerate near peak but accelerate as peak moves. (D) Social individuals (orange points) interacting with resource
peak (peak as described in C). Peak centroid moves according to 2D Brownian motion
with drift. Direction of peak motion indicated by grey arrow. Tracking occurs through
a dynamic process. As peak moves, individuals at peak’s trailing edge experience a
resource value that becomes weaker through time. As resource value becomes weaker,
these individuals accelerate as dictated by Equation 5.1, but turn toward neighbors
on the peak and end up traveling in an arc toward the leading edge of the moving
peak. This allows the group to track the peak as it moves through the environment.
115
5.3.2
Individuals in large groups fail to explore the environment
One disadvantage of being in a large group is that, the larger the group, the smaller
the fraction of group members that can be near the peak center at any point in
time. As in real animal groups, crowding increases resource competition. This leads
to an immediate question: by remaining near their neighbors, do social individuals
over-exploit nearby resources?
To answer this question, we studied populations in environments with two identical
resource peaks. Individuals begin the simulation near one of the peaks (referred to as
the starting peak) and both peaks move according to the same stochastic rules, defined
above. We performed simulations over a range of values of the ψ0 parameter, which
determines individuals’ default speeds in the absence of environmental cues. When ψ0
is small, the magnitude of the first term on the r.h.s. of Equation 5.1 will always be
small and changes in velocity (and speed) will be dominated by social influences (the
second term on the r.h.s. of Equation 5.1). When ψ0 is large, individuals accelerate
rapidly when resource value is low and are therefore more likely to break away from
groups. Below, we describe results of simulations with two peaks, which is the simplest
case where individuals face an explore-exploit tradeoff. In Section 5.5, we show that
the same qualitative results hold in more complex environments.
Fig. 5.3 shows that populations exhibit several behavioral regimes as a function
of ψ0 . For small values of ψ0 , individuals form a small group near the starting
peak (Fig. 5.3A, blue points), and experience a relatively low mean resource value
(Fig. 5.3B, blue points). As ψ0 increases, the number of individuals that track the
starting peak increases rapidly (Fig. 5.3A) leading to an increase in mean resource
value of the individuals near the peak center (Fig. 5.3B), however, as ψ0 reaches
roughly 1, the performance of individuals in the group near the peak begins to decline (Fig. 5.3B, ψ0 roughly 1-2), although the number of individuals that track the
116
peak continues to increase. The reason for this decrease in performance is intuitive:
as the group near the starting peak becomes larger, more individuals are forced to
remain around the periphery of the peak, where they experience low resource values.
Individuals fail to aggregate near the second peak in this regime (Fig. 5.3A red points,
group size on second peak is near zero for ψ0 < 2.2).
117
A
mean.resource-.groups.on.peaks
B
group.size.on.peaks
300
group.size.on.peak.A
group.size.on.peak.B
250
200
150
100
50
0
0
1
2
3
4
5
6
3.5
signal.of.group.on.peak.A
signal.of.group.on.peak.B
3
2.5
2
1.5
1
0.5
0
0
1
2
3
ψ0
4
5
6
ψ0
C
mean.resource,.all.agents
1.5
explore-exploit
explore
exploit
1
0.5
0
3
4
5
6
0
1
2
3
ψ0
4
5
6
Figure 5.3: Group size and performance in a two-peak environment (exploration and
exploitation). (A) The number of individuals on each peak (the starting peak, blue;
the second peak, red) as a function of the exploration parameter, ψ0 . Below a ψ0 value
of roughly 2.2, individuals do not form a group on the second peak. (B) Mean resource
value of individuals on the starting peak (blue) and second peak (red). (C) Resource
value averaged over all individuals in the population (individuals in groups nearest
each peak + all other individuals in the environment). Note maximum value occurs
in the regime where individuals aggregate on both the starting and second peaks
(ψ0 ∼ 2.5 − 3). Semitransparent points are results of 2000 individual simulation runs.
Runs were divided into 50 evenly spaced bins. Bolded points and error bars show
mean of each bin ± 2 standard errors.
118
5.3.3
Individuals maximize performance by balancing exploration and exploitation
As ψ0 reaches approximately 2, individuals begin to form a group on the second
resource peak (Fig. 5.3A, red points). Mean performance of individuals in the group
nearest the second peak rapidly increases (Fig. 5.3B, red points). When performance
is averaged over the entire population, there is a clear optimal value at ψ0 ≈ 2.7
(Fig. 5.3C). In this regime, the sizes of the groups nearest each of the two peaks are
almost exactly the same (Fig. 5.3A). For larger values of ψ0 , the average performance
over all individuals begins to decline (Fig. 5.3C), due to over-expoloration. The
phenomena shown in Figure Fig. 5.3 suggest that individual behavior falls into three
qualitative regimes. At low ψ0 , individuals exploit nearby resources, but fail to find
other peaks in the environment (Fig. S5.5 shows similar results in more complex
environments with more peaks). At high ψ0 , individuals encounter both peaks but fail
to exploit them effectively (Fig. 5.3C). For intermediate values of ψ0 , individuals track
the starting peak but also find and track the second peak (Fig. 5.3A). In this regime,
individuals balance exploration of the environment, which primarily occurs in small
groups, and exploitation of resource peaks, which is more effective in larger groups;
mean individual performance is maximized in this regime. Evolutionary simulations
confirm that ψ0 values in this optimal range quickly dominate a population that
is initially heterogeneous (Section 5.5.E, Fig. S5.9). Individuals with high and low
values of ψ0 are quickly removed from the population through selection.
5.3.4
Populations of optimal individuals exhibit fissionfusion dynamics
In the parameter regime that maximizes individual performance, individuals around
the periphery of the starting peak break away and move independently or form small
119
groups. These individuals begin to fuse to form another group around the second
peak. Fission occurs because individuals that are experiencing low resource values
have high default speeds and the first term on the r.h.s. of Equation 5.1 dominates
social interactions. These individuals break away from the group on the peak. Fission
occurs near the second resource peak because individuals slow down near the peak
and are more strongly influenced by the social interactions represented by the second
term in Equation 5.1.
When the sizes of resource peaks are not equal, individuals form groups whose
sizes correspond to the length scale of the environment (that is, the size of the resource
peaks, see Model development). Fig. 5.4A,B shows the sizes of groups on each of the
two peaks in an environment where the sum of the peak size is held constant. Beyond
a critical size (about 100 in the case of the parameters used in Fig. 5.4A), the number
of individuals on the starting peak increases nearly linearly as a function of the size of
the peak (Fig. 5.4A, red points). At the same time, the number of individuals on the
second peak decreases as the size of that peak decreases (Fig. 5.4B, blue points). The
simple decision rules captured in our model allow individuals to form groups of many
different sizes, but critically, these sizes match the length scales of the environmental
features that individuals interact with. Rather than being rigid, these decision rules
confer enormous flexibility at the level of the population. This flexibility is evident
in Fig. 5.4C,D,E, which shows the sizes of the five largest groups (where the size is
indicated by the width of the bands) over time. In the optimal ψ0 regime, the size of
the largest group shrinks, while the second-largest group grows until the group sizes
are approximately equal (Fig. 5.4C). When ψ0 is below the optimal range, individuals
remain in a single large group (Fig. 5.4D). When ψ0 is above the optimal range,
individuals form many small groups of similar size (Fig. 5.4E). Collective fission and
fusion is even more evident when resources are depleted because the distribution of
resources is constantly changing (Section 5.5.D).
120
C
A
200
150
100
50
0
0
200
400
600
mass6of6peak6A
800
1000
D
E
explore
B
explore-exploit
250
exploit
size6of6group6on6peak
size6of6group6on6peak6A
size6of6group6on6peak6B
relative6sizes6of656largest6groups
300
0
2500 5000 7500 10000
time
Figure 5.4: Peak size matching and group fission-fusion. (A) Peak size matching when
peaks differ in size. Total mass (i.e., resource value integrated over the domain) of two
peaks sum to 1000; the larger the starting peak (x-axis), the smaller the second peak.
Group size is mean size of group nearest the starting peak (blue) and second peak
(red; mean taken over the last 2500 steps of each simulation). Bolded points (and
error bars) represent mean (± 2 standard errors of the mean) of 1000 simulations.
Individual simulation runs shown by semitransparent points. Each point represents
a simulation in which the environment contains two peaks that differ in size, 300
individuals, over 5000 time steps, with ψ0 = 2.7. (B) Illustration of group sizes on
peaks of varying sizees. (C-E) Relative sizes (size indicated by height of bands) of
the five largest groups in populations with (C) ψ0 = 2.7 (in the optimal range), (D)
ψ0 = 2.0 (below the optimal range), and ψ0 = 5.7 (above optimal range). In panels
C-E, at each point in time, thickness bands denote sizes the five largest groups in the
population. Groups ranked from largest (top) to smallest (bottom) in each panel.
Increase in band thickness denotes fusion. Decrease in band thickness denotes fission.
In optimal regime (panel B), note fission of the largest group and simultaneous fusion
into the second-largest group (top band and second from top).
121
5.3.5
Optimal individual behavior places populations near a
state transition
The values of model parameters that maximize individual performance are surprisingly robust; the same model parameters yield maximal or near-maximal performance
over a wide range of environmental conditions (Section 5.5.C), and enable populations
to form groups that match the sizes of resource peaks (Fig. 5.4A,B). In addition, it
is qualitatively evident that the optimal parameter regime places populations near a
state transition (see Fig. 5.5A, which shows the mean inter-individual distance (measured as mean distance between all pairs of individuals in the population) over a range
of values of ψ0 and the ratio of social repulsion to attraction strength, C = Cr /Ca ).
Changing C and ψ0 moves the group state across a state transition (from a grouped,
coherent state, to a gaseous, dispersed state), in a featureless environment (no resources). In fact, the system dynamics exhibit qualities which are reminiscent of a
first order phase transition, with a prominent hysteresis effect evident when changing between states (Fig. S5.11). For a given value of C, there is a small range of
the exploration parameter, ψ0 , over which the mean inter-individual distance changes
dramatically. This same region of parameter space yields the highest mean individual
performance (Fig. 5.5).
To see why this transitional parameter regime is associated with the highest performance, consider a simplified version of Equation 5.1 in which S is constant, and
therefore, dva , i/dt = (ψ − η|vi |2 )vi /|vi |, where ψ = ψ0 − ψ1 ∗ S. In this formulation,
dva,i /dt no longer varies as a function of an individual’s position. When ψ/η is small,
default speed is low and individuals form solid-like groups with well-defined spacing between individuals [D’Orsogna et al., 2006]. For large ψ/η, individuals travel
p
ψ/η and exhibit gas-like behavior; social interactions are
at speeds near |v| =
weak relative to autonomous acceleration, and individuals explore the environment
relatively independently. In the critical regime, the behavior of individuals changes
122
rapidly, from solid-like to gas-like, for small changes in ψ (Fig. 5.5A). Now assume,
as we do above, that resource value varies over the environment (i.e., S(x) is high in
some regions and low in others). By decelerating in response to high resource level,
individuals move from the region of ψ that leads to gas-like behavior to the region in
which individuals form larger coherent groups (i.e., individuals move along the critical
region of curves in Fig. 5.5C). This property allows individuals to explore effectively
but also to fuse and exploit favorable regions of the environment.
123
5
C
4
3
2
1
B
1
1.5
2
2.5
3
3.5
4
4.5
5
5.5
150
100
50
5
mean resource value
4.5
4
C
3.5
3
2.5
2
1.5
mean inter-individual distance
1
C
mean inter-individual distance
A
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
1
1.5
2
2.5
3
3.5
4
4.5
5
5.5
200
C =1
C = 1.4
C = 1.9
C = 2.3
C = 2.8
C = 3.2
C = 3.7
C = 4.1
C = 4.6
C =5
150
100
50
0
0
1
2
3
4
5
6
Figure 5.5: Optimal behavioral rules occur near a critical point. (A) Mean interindividual distances (color scale) in a uniform environment (S = 0) as functions of the
behavioral parameters C = Cr /Ca and ψ0 . Note thin region of parameter space where
inter-individual distance changes rapidly. (B) Mean resource value of population as
functions of parameters C and ψ0 . Note that the region of highest performance
(dark red) matches the region of rapid change in inter-individual distance in (A). (C)
Optimal values of from panel (B) superimposed on curves from cross-sections of panel
(A) showing that optimal values of ψ0 occur in the region of parameter space where
small changes in ψ0 elicit large changes in inter-individual distance.
124
5.4
Discussion
Our results show that group size can be passively regulated by individual behavioral responses to social and environmental cues. Group fission and fusion arise
from these responses in a dynamic environment. The simple behavioral rules captured in our model, which are consistent with observed decision-making behavior
([Berdahl et al., 2013, Katz et al., 2011], allow animals to track dynamic resource
peaks, but also to explore the environment sufficiently to discover new resources.
The performance of groups in search and tracking tasks has been shown to depend
strongly on the relationship between group size and the length scale of environmental
gradients [Berdahl et al., 2013, Torney et al., 2009]. Here, the rules individuals use
to make movement decisions result in group sizes that change passively to match the
length scale of resources in the environment.
Modern techniques for simultaneously observing animals in environmental gradients in situ (e.g. prey distributions, eddies, oxyclines) suggest that these gradients strongly influence how animals are distributed in space [Bertrand et al., 2014,
Bertrand et al., 2008, Kai et al., 2009]. For example, predatory sea birds, fish, and
macrozooplankton in the Northern Humboldt Current System off Peru all aggregate
within small, dynamic oceanographic structures, which become hotspots of feeding
activity [Bertrand et al., 2014]. While there is evidence that collective behavior is
likely involved in these aggregations [Bertrand et al., 2008], we know very little about
precise behavioral mechanisms that could enable animal groups to locate and track
such ephemeral structures. Our results suggest that local, social interactions combined with a simple environmental response may be all that is required. It will be
interesting to merge our approach with continuum models of group fission-fusion (e.g.
[Bonabeau et al., 1999, Ma et al., 2011, Niwa, 2004]) that seek to predict the distribution of group sizes in a population, but do not explicitly consider responses of
individuals to structure in the environment.
125
One of the insights of our model is that individual performance is optimized when
the population undergoes a state transition at the border of favorable regions of the
environment. The fact that populations of agents that exhibit this behavior readily
evolve and remain stable through evolutionary time (see Section 5.5.E) indicates
that this population-level property can arise through selection on the behaviors of
individuals.
The tradeoff between exploration and exploitation is at the heart of decisionmaking problems [Cohen et al., 2007, Sutton and Barto, 1998]. To survive and reproduce, animals must employ behavioral rules that allow them to balance this tradeoff.
The model studied here demonstrates that social interactions allow animals to effectively manage this tradeoff to both find and track dynamic resources without memory
or a cognitive basis for associating measured signal values with spatial information
about the resource. Our results suggest that social individuals can balance exploration and exploitation using behavioral rules that are much simpler, and cheaper
(in terms of neural computation, for example) than they would require if they were
to ignore other individuals. We suggest that fission-fusion dynamics are an outcome
of this balance. Our study scenario bears a strong resemblance to dynamic coverage
problems in distributed control theory [Lekien and Leonard, 2010] and dynamic optimization problems in complex problem spaces [Passino, 2002]. The algorithms we
discuss may therefore be useful starting points for optimization schemes or control
algorithms for autonomous vehicles.
126
Parameter name
Default value
Number of agents (N)
300
Repulsion coefficient (Cr)
1.1
Attraction coefficient (Ca)
1
Repulsion decay length (lr)
.133
Attraction decay length (la)
1
Agent mass (m)
1
Resource response gain (1 )
2
Default speed (0 )
3
Friction ()
1
Resource height (0 )
10
Resource decay length (1 )
2
Agent density ()
0.0025
Number of resource peaks
2
Topological interaction radius
25
Metric interaction radius
20
Figure 5.6: Default parameters used in simulations, unless otherwise specified.
5.5
5.5.A
Technical details
Tracking dynamic signals
Tracking occurs through a dynamic process. As the peak moves, individuals at the
peak’s trailing edge experience a signal that becomes weaker through time. As the
signal becomes weaker, these individuals accelerate as dictated by Equation 5.1- they
are coupled to the environment through acceleration. However, they also turn toward
their neighbors still on the peak, and end up traveling in an arc toward the leading
edge of the moving peak (Fig. 5.2C, trails show path individuals followed to reach
current positions). Because of this dynamic, individuals track the peak as it moves
through the environment. In addition, social agents maintain a significant advantage
over asocial agents, even as the peak’s motion becomes fast (Fig. 5.2B) and stochastic
(Fig. S5.1).
127
3
social agents
asocial agents
mean signal value
2.5
2
1.5
1
0.5
0
0
0.2
0.4
0.6
0.8
1
peak noise
Figure S5.1: Tracking with increasing peak noise (the noise here is the variance
of symmetric 2-D Gaussian distribution from which noise in the peak position is
generated at each time step). Social and asocial agents perform comparably well
(given infinite time, all asocial agents would eventually find the stationary peaks)
for stationary peaks (300 agents in a 5-peak environment). Mean resource value of
asocials decreases rapidly with increasing peak noise. Decrease in performance for
social agents is much less rapid.
Critically, individual searchers are able to remain near a moving resource patch
without any knowledge of where the patch is moving, how large the patch is, or
where they are relative to the patch’s center. The only information individuals have
about the resource is the instantaneous resource value at their own positions. Possessing only local environmental information is certainly a reality for many organisms.
For example, some sensory modalities that animals use to perceive environmental
conditions, such as the fish lateral line system, are capable of measuring conditions
only in the animal’s immediate vicinity [Coombs, 1994]. Moreover, for animals such
as schooling marine fishes, crowding [Partridge, 1982] severely restricts individuals’
abilities to perceive the environment at a distance. Our model predicts that even
when individuals can only detect the scalar value of a resource at their own exact
128
positions (i.e. local environmental detection only), groups of such individuals can still
successfully track a moving resource.
5.5.B
Alternative peak shapes
We check our results using an alternative shape for the resource patches, to ensure
that the results we find are not an artifact of having an explicitly exponential peak.
We consider a Gaussian resource patch, where we modify Equation 5.5 as follows.
Sg (x, xs ) = λ0 exp−|x−xs |
2 /λ
1
(S5.6)
In Fig. S5.2 we plot the mean performances of agents in the 2-peak Gaussian
environment, again varying ψ0 from 0 to 6 (similarly to Fig. 5.3. As expected, the
behavior changes with the choice of peak shape parameters, but we find there exist
parameter values for this alternate shape for which the results from the exponential peak are replicated. We use a larger λ1 for the Gaussian resource, because the
Gaussian falls off faster than the exponential. Agents in environments with Gaussian
resources display the same three regimes- exploiting at low ψ0 , balancing exploitation
with exploration at mid-range ψ0 , and exploring at high ψ0 . Thus we gain even more
confidence that our model is robust to a variety of environmental conditions.
129
Gaussian.resource.peaks
B
200
group.size.on.peak.1
group.size.on.peak.2
group.size.on.resource.peaks
180
160
140
120
100
80
60
40
20
0
0
1
2
3
mean.resource,.all.agents
ψ0
4
5
6
4
ψ0
C
ue.on.peak.1
ue.on.peak.2
3
mean.resource.-.groups.on.peaks
A
5
6
5
signal.value.on.peak.1
signal.value.on.peak.2
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0
0
1
2
4
5
6
3
ψ0
4
5
6
3.5
3
2.5
2
1.5
1
0.5
0
0
1
2
3
ψ0
Figure S5.2: This figure is equivalent to 5.3 in the main text, with 2 Gaussian peaks
(where they are 2 exponential peaks in the main text. The signal obtained by agents
from these Gaussian peaks. (A) Size of group on central peak (peak 1), blue; size of
group on secondary peak (peak 2), red. (B) Mean signal value of group on peak 1,
blue; mean signal value of group on peak 2, red. (C) Mean signal value of all agents
(including those not on either peak). Behavior for Gaussian peaks is essentially
identical to behavior for exponential peaks. Semitransparent points are results of 200
individual simulation runs, divided into 30 evenly spaced bins. Bolded points and
error bars show mean of each bin, ± 2 standard errors.
130
5.5.C
Sensitivity of optimum value
In order to establish that our results are not due to fine-tuning of some of the parameters in our model, we perform a sensitivity analysis to see how the location of
the optimal ψ0 shifts as we vary a selection of parameters. First we look at the sensitivity to characteristics of the environment. In Fig. S5.3, we vary peak speed and
peak noise (in a 2-peak environment), and find that the overall mean signal value of
the group decreases as the peak speed and peak noise increase, as one would expect
because the resource patches become harder to track when they move faster and less
predictably. As peak speed increases, the location of the optimal ψ0 remains essentially unchanged, as long as the peak speed is still low enough for agents to find and
track the secondary peak. If the peak moves faster than a certain level (peak speed
> 0.02), the peak’s motion overtakes the default speed of low ψ0 agents, so that these
individuals fail to keep up with either peaks motion, and are unable to track it. At
very high speeds, only agents with values of ψ0 near the optimum gain any benefit
from their environment. Increasing peak noise, however, does result in a downwards
shifting of the optimal ψ0 , although the effect not drastic. This means that the more
stochastic peak movement is, the more tracking is degraded in agents with high levels
of exploration (high ψ0 ). Once the peak’s stochasticity crosses a certain threshold
(peak noise > 0.25), the agents fail to track the secondary peak, although they are
still able to track the central peak, because their initial positions allow them to locate
it immediately. This results in the disappearance of the optimal “explore-exploit”
regime we observe with less peak stochasticity.
131
B
A
peak7drift7=70
peak7drift7=70.00
peak7drift7=70.01
peak7drift7=70.02
peak7drift7=70.03
peak7drift7=70.04
peak7drift7=70.05
peak7drift7=70.06
peak7drift7=70.07
peak7drift7=70.08
2.5
2
1.5
mean7resource7value,7all7agents
mean7resource7value,7all7agents
3
1
0.5
0
0
1
2
ψ0
3
4
3
peak7noise7=70
peak7noise7=70.04
peak7noise7=70.08
peak7noise7=70.13
peak7noise7=70.17
peak7noise7=70.22
peak7noise7=70.26
peak7noise7=70.31
peak7noise7=70.35
peak7noise7=70.4
2.5
2
1.5
1
0.5
0
0
5
1
2
ψ0
3
4
5
Figure S5.3: Sensitivity to peak mobility. (A) mean signal value plotted against ψ0 ,
for 10 different values of peak drift speed. The location of the optimal ψ0 is not
strongly dependent on how fast the resource patches drift. (B) Mean signal value
plotted against ψ0 , for 10 different values of peak noise. The optimal value for ψ0 here
drifts lower with increasing peak noise, although the effect is not very strong. For
both cases, there is a certain value for peak drift speed and peak noise above which
agents can’t track the secondary resource patch. Once the peak moves too fast, the
location of the peak changes faster than the time required for a large enough group of
agents (to enable tracking) can find it. Once beyond this threshold, the added benefit
from having a second resource patch in the environment disappears.
Next, we consider the effect of changing peak size on the optimal value of ψ0 ,
again in a 2-peak environment. In Fig. S5.4 we plot the mean signal value normalized
by peak size (2πλ0 λ21 , so that we plot mean signal per unit resource, as the total
resource volume is proportional to peak size. We find that as peak size increases,
the location of the optimal ψ0 shifts upward, so that if there is more resource to go
around, it pays off to be more exploratory. This effect is probably an artifact of the
peak shape we’ve chosen here. Because the peaks are exponential, they are sharply
pointed in the center, and agents are slowed down to zero speed before they reach
the center, due to the way agent speed is coupled to the environment . Agents with
higher values of ψ0 will reach more central values of the peak before they are slowed
to a halt, which results in better performance for these agents with higher ψ0 than
132
comparable agents with lower ψ0 , even though both types would be able to track the
peak equally well.
5
0.8
4.5
0.7
mean signal/peak size
peak size
4
3.5
3
2.5
2
1.5
1
0.6
0.5
0.4
0.3
0.2
0.1
0
1
2
3
4
ψ0
5
6
7
8
9
Figure S5.4: Sensitivity to peak size. Mean signal value normalized by peak size
(≈ λ21 from 5.5, main text) is indicated by color, where warmer colors correspond to
higher performance. We normalized by peak size here to account for the increase in
total resource in the environment when peak size is higher (volume of peak increases
with λ21 because the peaks are exponential). The location of the optimal ψ0 increases
with peak size, so that it pays to be more exploratory when there is more resource
available. We ran simulations for 10 different peak sizes for each of 200 different
values of ψ0 . Each simulation consisted of 300 agents in a 2-peak environment for
10,000 timesteps
The final characteristic of the environment that we explore sensitivity to is the
number of peaks. The optimal value for ψ0 is essentially agnostic to number of peaks
in the environment (Fig. S5.5), at least in the tested range of 2 through 10 peaks.
There are diminishing returns to adding more peaks to the environment, because once
the number of peaks surpasses a certain level, there is more than enough resource for
all agents, and thus there is no overcrowding. In the case of an infinite number of
agents, there would be no diminishing returns to adding more peaks; agents would
always benefit from a larger number of peaks.
133
2.5
n peaks = 2
n peaks = 3
n peaks = 4
n peaks = 5
n peaks = 6
n peaks = 7
n peaks = 8
n peaks = 9
n peaks = 10
<signal value>
2
1.5
1
0.5
0
−1
0
1
2
3
ψ0
4
5
6
7
Figure S5.5: Sensitivity to number of peaks. We vary the number of peaks in the
environment from 2 to 10, where each peak moves in parallel motion to avoid overlap,
and compare the average group performance. Returns diminish as more peaks are
added, as once there is a certain amount of resource available in the environment, there
is more than enough for all agents, and crowding ceases to be an issue. The optimal
ψ0 is once again essentially independent of the number of peaks in the environment.
Each simulation consisted of 300 agents in a 2-peak environment for 10,000 timesteps
In addition to varying environmental characteristics, we also explore sensitivity
of the optimal ψ0 to the number of agents in the simulation (where we return to
the 2-peak scenario for this analysis). As the number of agents increases, the mean
performance of the group decreases monotonically (Fig. S5.6A). This is because the
resource level is held constant, so as the number of agents goes up, the amount of
available resource per agent decreases. We see this more clearly if we look at the
total signal value, instead of the mean signal value (Fig. S5.6B). Small groups do not
fully use the available resources, so adding more agents results in increasing levels of
total signal value. However, at larger group sizes, the total signal value saturates,
and adding more agents does not result in more resource consumption. Essentially, a
finite number of agents can fit on the resource patches. Apart from the smallest two
group sizes (N = 50, N = 100), the location of the optimal ψ0 is invariant to group
134
size. In the two smallest group sizes, the agents are unable to track either peak with
low values of ψ0 , because the radius of the group they form is too small compared to
the radius of the resource peak, and the tracking mechanism breaks down.
A
B
3
400
N = 50
N = 100
N = 150
N = 200
N = 250
N = 300
N = 350
N = 400
N = 450
N = 500
2
1.5
350
total resource value
<resource value>
2.5
1
0.5
0
−1
300
250
200
150
100
50
0
1
2
3
ψ0
4
5
6
0
−1
7
0
1
2
3
ψ0
4
5
6
7
Figure S5.6: Sensitivity to number of agents. (A) Mean signal value of all agents
plotted as a function of ψ0 . Average performance of group decreases as the group
gets bigger, because the total available resource remains the same, and has to be
spread amongst a larger pool of individuals. (B) Total signal value (as opposed to
average signal value in panel (A)), plotted against ψ0 . Total signal value increases as
more agents are added, for small group sizes, because these small groups do not fully
exploit the available resources. For large group sizes, however, the total signal value
saturates. Location of optimal ψ0 is essentially agnostic to number of agents. We ran
simulations for 10 different number of agents for each of 200 different values of ψ0 , in
a 2-peak environment for 10,000 timesteps
Finally, we consider sensitivity of the optimal ψ0 value to the value of ψ1 , which
determines the coupling strength of agents to their environment (see Equation 5.1).
We again simulate groups of 300 agents in a 2-peak environment, and measure their
average signal value, varying both ψ0 and ψ1 . If ψ1 is too low (Fig. S5.7), agents
don’t couple strongly enough to the environment to amass a large enough group on
the secondary peak to be able to track it (they are able to track the starting peak
because the simulations were initialized such that the group was placed on top of
it). Once ψ1 is high enough to enable tracking of the secondary peak, increasing ψ1
135
further has the effect of shifting the optimal value of ψ0 slightly higher. This effect
diminishes with high ψ1 , with the optimal ψ0 converging to approximately 3.5, for
the highest value of ψ1 tested (Fig. S5.7).
B
A
1
0.8
0.6
1.2
3.5
1
3
ψ1
1.2
<resource value>
4
ψ1 = 1
ψ1 = 1.3
ψ1 = 1.7
ψ1 = 2
ψ1 = 2.3
ψ1 = 2.7
ψ1 = 3
ψ1 = 3.3
ψ1 = 3.7
ψ1 = 4
<signal value>
1.4
2.5
2
0.4
1.5
0.2
0
−1
0
1
2
3
ψ0
4
5
6
1
7
0.8
0.6
0.4
0.2
0
1
2
3
ψ0
4
5
6
Figure S5.7: Sensitivity to ψ1 . (A) Average signal value of the group plotted against
ψ0 , for 10 values of ψ1 , ranging between 1 and 4. ψ1 dictates the strength of coupling
between the agents and the environment. Higher values of ψ1 result in greater sensitivity to a small environmental signal, and so lead to increased agent performance
and better tracking capabilities. (B) Heat map showing group performance (warmer
colors indicate better performance), as a function of ψ1 and ψ0 . The optimal value for
ψ0 is relatively insensitive to ψ1 . We ran simulations for 10 different number of agents
for each of 200 different values of ψ0 , in a 2-peak environment for 10,000 timesteps
5.5.D
Resource depletion
One natural modification to our model is to add a resource depletion mechanism to
the resource patches. Although the resource patches in our model need not solely
be thought of as food (they could also be shelter, as in [Berdahl et al., 2013]), if
one does consider scenarios in which the resources are gathered (such as a source
of energy/nutrition), the resource will naturally be depleted over time. The rate of
consumption can be reasonably approximated as a function of the size of the group
which is consuming it. Here we explore the effects consumable resource patches
136
have on our results, and ask whether the agents in our simulations retain the ability
to balance exploration with exploitation when they possess a particular trait (i.e.
a particular value of ψ0 ), and if so, how that value is different from the case of
non-consumable resources. We perform this experiment using resource patches with
heights that decrease at a rate (with time constant γ) proportional to the number
of agents currently on each peak (NA ). Once the resource has been depleted below
a threshold level (λmin ), it is destroyed and a new peak is generated at a random
location in the environment, with peak height = λ0 . This depletion for a peak, A, is
organized as follows.
λ0A,0 (t) =



 λ0 , if λ0A,0 (t − 1) < λmin
(S5.7)


 λ0A,0 (t − 1) expNA /γ , otherwise
In this experiment we initiate two peaks at random locations in the environment,
no longer having a central peak as in previous experiments, thus recovering the symmetry between both peaks. We select parameters for γ and λmin such that agents
retain the ability to locate and track peaks, and that peaks are destroyed before becoming too small to be detectable. We find that agents are able to find new resource
peaks relatively quickly , compared to the length of simulation run. After running
simulations for 200 different values of ψ0 , we find that the optimal value is almost
exactly the same in the resource depletion case (Fig. S5.8) as in the non-depleting resource case (Fig. 5.3C). That is, the same level of exploration is required to optimally
gain resources from an environment with consumable resources as from an environment where the resource is infinite. Note that the overall performance of the group
is much lower when there is resource depletion, because the total available amount of
resource is lower, on average, than when the resource is non-consumable.
137
0.6
mean resource value
0.5
0.4
0.3
0.2
0.1
0
0
1
2
3
ψ0
4
5
6
7
Figure S5.8: Explore/exploit with depleting resources. Mean signal for all agents
plotted against ψ0 , in environment with depleting resource patches. Height of each
peak (λ0 in Eq. 5.5, main text) decreases at each timestep, at a rate proportional
to the number of agents on each peak. Once this peak height decreases below a
threshold, it disappears and a new peak is generated at a random location in the
environment. There are two non-central depleting resources in this experiment, so
the asymmetry between peaks 1 and 2 disappears, as there is no longer a constant
central peak. The location of the optimal value for ψ0 is essentially unchanged, even
in this resource depletion experiment, although the mean performance of the group
is much lower, as there is less total available resource in the environment due to the
depletion.
5.5.E
Evolutionary Algorithm
In order to test whether optimal behavior can evolve in a heterogeneous population, we implement a simple evolutionary algorithm, similar to that used in
[Guttal and Couzin, 2010]. In the first generation, N (N = 300) agents are initiated,
in a 2-peak environment, with heterogeneous values of ψ0 , drawn from a uniform
random distribution ranging from 0 to 6. The number of agents remains constant
across generations, and the generations are non-overlapping. Each generation con138
sists of a simulation run for 10,000 timesteps, as in prior numerical experiments
described. At the end of each generation, N selections are made from the population,
with an agent’s probability of being selected for reproduction proportional to its
relative performance across the simulation, compared to the other agents. Individuals
which perform well are more likely to be selected to reproduce, and are likely to
produce more offspring, than those individuals which perform poorly. The selection
probability of agent i, P (ai ), is defined as follows.
P (ai ) =< S(ai ) >t /
N
X
< S(aj ) >t
(S5.8)
j=1
Where S(ai ) is the instantaneous signal value of agent i, and angular brackets with
subscript t denote time-averaging over generation number t. If an agent is selected
for reproduction, a child is produced in the next generation with ψ0,child = ψ0,parent ,
with a small mutation. The mutated ψ0 is equal to the original ψ0 , added to a
normally distributed random number, ω, with width, σ, that controls how fine-tuned
the resulting population can become.
ψ0 (a0i ) = ψ0 (ai ) + ω/σ
(S5.9)
Where the prime indicates a child. In Fig. S5.9, we plot the evolution of the
exploration parameter, ψ0 , as a function of the generation number. We plot this as
a heatmap, with ψ0 on the y-axis, and where warmer colors indicate a more popular
choice of parameter. White areas represent sections with no data; that is, where
no agents in a given genereation were selected to have that particular parameter
value. Distributions of ψ0 from three selected generations are shown as insets in
Fig. S5.9. Agents with values of ψ0 < 2, and ψ0 > 3.5 are quickly extinguished
139
from the population, occurring in only a few generations. As the generation number
increases, the population quickly collapses, so that the large majority of agents possess
values of ψ0 = 2.9 in the last generation run (generation 755). This corresponds almost
exactly to the location of the optimal ψ0 in Fig. 5.3C, main text. As the evolutionary
simulation is allowed to progress, we do notice some small fluctuation in the value of
the most popular ψ0 . This happens because of the mutation step- in each generation
some children are produced with values differing from their parents, and because
the performance near the optimal value is not that much worse than performance
exactly at the optimum. This means that some agents with near optimal values of
ψ0 can persist for a few generations. Ideally we would like to run this evolutionary
algorithm for a much larger number of generations (on the order of thousands instead
of hundreds), but our current implementation is limited in the number of generations
it can simulate in a reasonable amount of time. Since each generation relies on the
results of the agents in the previous generation, options for parallelizing the code
are limited. However, we feel confident that the results would not change very much,
even in the case of thousands of generations, due to the speed with which the exploreexploit range (from Fig. 5.3, main text) values for ψ0 are selected. Even by the 20th
generation, the large majority of agents have values of ψ0 between 2.5 and 3.5.
140
Figure S5.9: Evolving the optimal ψ0 . Heat map showing distribution of agents with
ψ0 falling within the bin on the y-axis (data divided into 25 bins), for each generation,
on the x-axis. Warmer colors represent more popular values of ψ0 , and white regions
are areas with no data (no agents have that particular value of ψ0 ). Insets show the
distribution of agent values of ψ0 , for generation 1, generation 378, and generation
755. The population quickly collapses so that most agents have near optimal values
of ψ0 .
5.5.F
Theoretical observations behind phase transitions in
model and balance of exploration with exploitation
In the model of collective motion described in this chapter, agents exert forces on
neighbors via the gradient of a social potential Equation 5.1. In the absence of
interactions with environmental resources, with only social forces, agents exist in two
phases, depending on the interaction coefficients and the amount of self-propulsion
present. If the self-propulsion of agents is low enough (in Equation 5.2, agents will
group together in a cohesive phase. If the self-propulsion becomes too high, however,
the agents break away and enter a gaseous, dispersed phase. Critically, when agents
modulate their speed in response to environmental resources, the phase they occupy
depends on their default speed as well as the value of resource they experience (Ψ =
ψ0 − ψ1 S(x)). Agents in the dispersed gaseous phase are able to rapidly explore
their entire environment, unlike agents in the cohesive phase, while agents in cohesive
groups are able track resource peaks that are both drifting and diffusing. A simple
141
mathematical model is derived using the original self-propelled particle model that
shows that the agents maximize their ability to take advantage of their environment
when their parameter values are such that a phase transition occurs at the boundary
of the favorable parts of the environment. The existence of separate phases encode
collective decisions of the group, emerging only from the interaction laws of the agents
and local measurements of the environment.
In the cohesive phase, agents settle into configurations which minimizes this social
potential. We find numerical solutions to this minimized social potential per agent, as
a function of group size, as well as the characteristic group radius for this minimized
potential.
First, assume that there are a group of N agents, in the cohesive state, which
occupy a circular region with radius l, and thus having a density, ρ of
ρ=
N
πl2
(S5.10)
Given this density, we define another radius, lk as the radius within which we
expect to find K agents. lk is related to l in the following way
r
lk = l
K
N
(S5.11)
Then we can calculate the potential of an agent, i, which interacts with all neighbors within lk , by integrating the social potential function 5.1 over the area of the
circle.
142
Z
Φi = ρ
(−Ca e−|xi −x|/la + Cr e−|xi −x|/lr )dA
(S5.12)
|xi −x|<lk
This integral evaluates to:
−2Ca 2
Φi =
[(la +
l2
r
√K l
K
2Cr
lla )e N la ] + 2 [(lr2 +
N
l
r
√K l
K
llr )e N lr ] − Ca la2 + Cr lr2
N
(S5.13)
We numerically compute the radius l, which minimizes the potential energy per
agent, Φi , selecting parameters consistent with what were used in the simulation
experiments. In Fig. S5.10A we plot this minimum potential energy per particle as
a function of total number of agents (N ). The potential energy steadily decreases
until N = K (where K = 25 as in our simulations), at which point the potential
energy becomes invariant to group size. Similarly, in Fig. S5.10B, we plot the total
radius of the group of N agents, as a function of N , which minimizes the potential.
When N < K, the group radius is invariant to N , while once N > K, the group
√
density becomes constant, and the group radius increases with N , as N . This is
characteristic of topological interaction rules, as it leads to a fixed constant density.
With metric interactions, the density increases as the number of agents. Agents which
stay in the cohesive phase do not explore their environment effectively, as they only
occupy a small fraction of the available area. Though these cohesive groups may at
times drift slowly through the environment, the total explored area is much less than
when agents are in the dispersed gaseous phase. However, being in a cohesive group
is essential for being able to track moving resource peaks, as explained in the tracking
section, above.
143
B
0
11
2
10
4
9
group radius
potential energy per agent
A
6
8
10
12
7
6
5
14
16
8
4
0
50
100
150
200
number of agents
3
0
50
100
150
number of agents
200
Figure S5.10: (A) Potential energy per particle as a function of number of agents in
group. Cutoff at N = K indicates constant density.
(B) Group radius as a function
√
of number of agents. Radius is proportional to N once N > K
In fact, the system dynamics show a strong hysteresis effect. Agents that start
in a state of high energy (high Ψ, where Ψ = ψ0 here because the environment is
featureless), and thus start in the disperse state condense into a cohesive group at
a low value of Ψ as Ψ is decreased gradually. Meanwhile, agents that start with
low energy, and thus start in the cohesive state, maintain their cohesion through
relatively large values of Ψ, as Ψ is increased gradually. Here the order parameter
is the average potential energy per particle, and it is a proxy for the local density of
agents. Low potential energy indicates that agents are in the cohesive state, while
high potential energy means they are in the gaseous state. Agents in the gaseous
state need a nucleation site (such as would be supplied by a resource peak) in order
to be able to change states.
144
<potential energy>
0
−5
−10
−15
0
0.5
1
1.5
2
ψ
2.5
3
3.5
4
Figure S5.11: Hysteresis in state transitions as energy (Ψ) is raised and lowered. Low
potential energy indicates a cohesive state, while high potential energy indicates a
dispersed gaseous state.
Crystallization
When Ψ < 0, such as when agents are located in the interior portions of a resource
peak, the default kinetic energy of the agents approaches, and the self-propulsion
force is always frictional. The agents arrange themselves in a crystal structure which
minimizes their potential. In simulations with environmental resource peaks, agents
in the crystallized state form directly over the resource peaks, and are surrounded
by a cloud of particles in the cohesive state. As they have no motion, agents in the
crystal state do not explore their environment.
Critical angle for capture by resource peak
In this section, we explore further a simplified model of collective exploration and
exploitation, and calculate the critical angle for which an agent will be captured by
145
a resource peak currently occupied by N agents. In order to do this we first have to
make some assumptions:
1. The environmental response function, Ψ(r) is composed of a default background
speed (u, where u2 = ψ0 /η from 5.1) experienced by the agents in the absence
of a resource. This default speed is slowed when a resource is encountered;
Ψ(r) = ψ0 − Ae−r/λ1 , (where A = ψ1 λ0 /η from 5.1 and 5.5).
2. Particles can be in one of two phases; either captured by the peak as part
of a cohesive group, or off the resource peak in a dispersed phase. Particles
in the dispersed phase are distributed uniformly in space and direction. The
magnitude of their velocity is u when r > rM , where rM is the maximum
interaction radius of the potential of a resource peak. Particles in the cohesive
group produce an attractive potential on the peak, Φpeak (r).
3. Particles in the dispersed phase only interact with particles in the cohesive
group.
4. Friction is strong enough that we can assume particles always travel at
their speed resulting from interaction with the resources in the environment;
√
u2 − Ae−r/λ1 for r < rM and u for r > rM . That is, the potential from the
peak occupancy can only change the direction of the particle, but not its speed.
5. A particle will be captured by the peak if it has a trajectory that reaches
the radius of zero velocity. That is, it reaches the point where Ψ(r) = 0;
r0 = λ1 log(A/u2 ). This radius marks the location of the phase transition
between cohesive phase and dispersed gaseous phase.
6. The cohesive phase has a maximum density ρ, and the peak has a maximum
size given by r0 .
146
Imagine a particle reaches the radius of influence of a resource peak, rM , where it
begins to feel the peak’s force (see Fig. S5.12). If the angle of the particle’s trajectory
is sufficiently directed towards the resource peak, the particle will be captured. What
is this angle, ∆c where the particle will be captured?
Δ
r0
rM
Figure S5.12: Agent entering region where it feels peak influence. There is a group
of agents on the resource peak centered on the origin. These agents are contained
within a circle of radius r0 , corresponding to the zero velocity region. An agent enters
the region of peak influence, inside the radius rM , from the left. The angle of velocity
of the agent relative to the peak is ∆.
To solve this, we first write down the equations of motion in polar coordinates
√
for a particle traveling with speed v(r) = u2 − Ae−r/λ1 , in direction γ, and which
experiences the potential from other particles currently occupying the resource peak
of −∇Φpeak (r). This allows us to answer the following question: given initial radius
r = rM , initial angle relative to the origin θ0 , and initial direction (heading) γ0 , does
the particle reach the zero velocity radius, r0 ? The equations of motion are as follows:
dr
= v(r) cos(θ − γ)
dt
dθ
v(r)
=
sin(γ − θ)
dt
r
Φ0peak (r)
dγ
=
sin(γ − θ)
dt
v(r)
147
(S5.14)
(S5.15)
(S5.16)
Next we remove an angle by defining the angle ∆ = −π − γ + θ. This is the angle
between the heading vector and the vector directed from the position of the agent to
the origin. Thus we can rewrite our equations of motion in terms of r and ∆ alone,
resulting in the following:
dr
= v(r) cos(∆)
dt
Φ0peak (r) v(r)
d∆
= −(
−
) sin(∆)
dt
v(r)
r
(S5.17)
(S5.18)
We can discover a few things by examining this system of equations.
1. If ∆ = 0, then
implies that
dr
dt
d∆
dt
= 0, so the particle’s heading will remain 0. Similarly, ∆ = 0
< 0, so the particle’s radius will decrease until it reaches the
radius of zero velocity (r0 ).
2. If − π2 < ∆ < π2 , then
3. If both (
Φ0peak (r)
r
dr
dt
< 0, and the particle will move closer to the peak.
− v(r)
) and − π2 < ∆ < π2 , then the value of ∆ will become closer
r
to 0. That is, if both conditions are met, and ∆ > 0, then
and if ∆ < 0, then
d∆
dt
d∆
dt
will be negative,
will be positive.
We make one last assumption, which has been true in most practical cases, that
allows us to make progress with the analysis. We assume that there is exists a point
within such that the potential felt by the peak is balanced by the centripetal force of
the approaching particle, so that dΦpeak (r)/dr − v(r)2 /r has one sign change on the
interval (r0 , ∞). Thus we have the following inequalities:
Φ0peak (r)
v(r)
<
r
r
for r < r∗,
Φ0peak (r)
v(r)
>
r
r
148
for r > r∗
(S5.19)
This allows us to divide the (r, ∆) plane into four regions:
1. − π2 < ∆ < π2 , and r < r∗
2. ∆ < − π2 or
π
2
< ∆ and r > r∗
3. − π2 < ∆ < π2 , and r > r∗
4. ∆ < − π2 or
π
2
< ∆ and r < r∗
In region 1, the particles more towards the peak and the potential force is stronger
than the centrifugal force. In region 2, the agents move away from the peak, and the
potential force is weaker than the centrifugal force. In region 3, agents move towards
the peak but the potential force is weaker than the centrifugal force. In region 4, the
agents move away from the peak but the potential is stronger than the centrifugal
force. We conclude that any trajectory that enters region 1 will eventually reach the
zero velocity radius, and any trajectory that enters region 2 will escape to ∞.
Consider again the hypothetical agent in region 3 and at radius r = rM . The
agent either enters region 1, region 2, or the boundary points between the two regions
(which are unstable equilibria). The points (r∗, ± π2 ) are unstable equilibrium points,
each corresponding to a periodic orbit around the environmental peak. There are two
values of ∆ such that the solution of the initial value problem with initial conditions
(rM , ∆) reach these equilibria. We call these angles ±∆c , and any particle with initial
conditions (rM , ∆), with |∆| < ∆c will be captured, and with |∆| > ∆c will escape.
Instead of considering the time dependent differential equation, we will search for
an equation that describes the shape of a trajectory. That is, we will use the original
equation in combination with the chain rule to write a differential equation for ∆.
This method will only work in regions where ∆ is a single valued function of r, and
thus we must restrict our integration to the region r > r∗.
First we simplify the equations of motion by taking their quotient:
149
dr
=
d∆
cot(∆)
Φ0peak (r)
v(r)2
−
(S5.20)
1
r
This differential equation can be integrated, and we do so with the assumption
that a trajectory begins at (rm , ∆i ) and ends at (r∗, π/2). The resulting expression
is
Φ0peak (r) 1
− )dr =
(
v(r)2
r
rM
Z
r∗
Z
π/2
cot(∆)d∆
(S5.21)
∆c
In order to simplify the resulting expressions, we define the function F (r) to be
the anti-derivative of Φ0peak (r)/v(r)2 . For our choice of potential function and environmental signal, this integral can be evaluated exactly in terms of hypergeometric
functions:
Z
F (r) =
r
Ca N re−r/la
Ca N λ1 r( λ1 − l1 )
λ1 r/λ1 u2
λ1
a
1
=
e
)
,
2
−
,e
F
(1,
1
−
2 1
la (u2 − Ae−r/λ1 )
A(la − σ)
la
la
A
(S5.22)
The integral equation then simplifies to:
F (r∗) − F (rM ) − log(
r∗
sin(∆c
) = −log(
)
rM
sin(π/2)
We can solve this equation for the critical angle ∆c .
150
(S5.23)
∆c = sin−1 (
r∗ F (rM )−F (r∗)
e
)
rM
(S5.24)
This critical angle ∆c is a function of the parameters defining the agent behavior,
such as the background velocity, u, and the parameters defining the cohesive group,
such as the peak occupancy, N . When u is low, trajectories are more influenced
by the potential of the particles on the resource peak, so they are more likely to be
captured. As u gets very small, ∆c approaches π/2. As N increases, the potential
from the peak gets stronger, and ∆c gets larger for all values of u. However, when u is
too large, ∆c goes to 0, and no particles can be captured by the peak. In Fig. S5.13,
we include numerical solutions for ∆c , for different values of background speed u and
peak occupancy N .
1.6
N=0
N=1
N=5
N = 10
N = 15
N = 20
Crit ical Angle ∆ c
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
u
Figure S5.13: Solutions for critical angle, ∆c , as a function of the background speed,
u, for 6 different values of peak occupancy, N .
151
Chapter 6
Conclusions and outlook
The topics addressed in this work are the result of an underlying fascination with
all types of emergent systems. I was first attracted to the subject, as many are,
by personal observations in nature, of the impressive acrobatics of starling flocks,
and of the captivating collective movement of fish. These animals seemingly moved
with a single mind, making decisions and reacting to perturbations in mesmerizing
waves of motion, in the absence of any obvious leadership. In fact, as we investigate
here, a great amount may be learned about collective motion if we know more about
the interactions between individuals. Notably, in many animal systems, interactions
among individuals appear to produce remarkably robust collective patterns, and allow
for collective computation and efficient spread of information.
In Chapters 2 and 3 I demonstrated that by finding a functional mapping from sensory input to behavioral output, and reconstructing the resulting interaction network
in freely swimming golden shiner fish (Notemigonus crysoleucas), we gain insight into
how behavioral alarm cascades flow through groups. Notably, the interconnections
among an individual’s neighbors prove to be critical to whether the alarm will spread
extensively through the group, or will die out quickly. Simulation studies confirm
that in our experimentally derived networks, cascades that propagate through com152
plex fractional contagion best explain our data, and deepen our understanding of the
relationship between network properties, spatial structure, and social influence.
Ongoing work by colleagues in the Couzin lab group is focused on more controlled experiments, where these groups of fish are pushed into a more sensitive behavioral regime. By controlling environmental stimuli we may be able to explore
how alarms propagate through groups when individual fish are oversensitive to false
alarms. Golden shiners produce an alarm pheromone to which they react strongly
when it is present in high enough concentration. When this pheromone is added to
water, behavioral false alarms are observed to occur both more frequently and to
spread to a larger fraction of the group.
In Chapter 4 I described an ongoing project studying the collective oscillations
in activity observed in colonies of Leptothorax acervorum where preliminary results
have revealed that individual ants influence the activity states of those they interact
with. Additionally, we observed that these behavioral responses to interactions are
flexible, and vary over the course of the experiment, and that this response flexibility
can explain the variation in coherence we observe in the activity waves. Simulations
using a mobile cellular automaton model confirm that the phenomena we observe may
be driven by flexibility in responsiveness to interactions. Further work is needed to
complete this project, however, to determine the ramifications that flexible behavioral
responses have for collective oscillations. Furthermore, a more individual based approach may be possible with this data. That is, we may be able to determine whether
the flexibility we observe is consistent across the entire colony, or if it is driven by
certain types of ants, such as the brood-workers or the foragers. By coupling these
individual differences to the spatial fidelity zones of specific ants, we may be able to
gain a better understanding of how the activity waves work, as well as their functional
significance in regulating and amplifying aspects of colony life.
153
Chapter 5 explored how with a simple model of collective motion, where agents
are coupled minimally to their environment, groups can optimize the tradeoff between
exploration and exploitation of their environment. Additionally, groups dynamically
fuse and break apart, to result in group sizes which match the length scale of the environmental resources at hand. Furthermore, this optimal level of exploration places
groups near a state transition, allowing them to switch spontaneously between a cohesive group state when they are in a good part of the environment, and a dispersed
gaseous state when a resource is overcrowded, or they are in a bad part of the environment.
In each of the projects discussed in this work, we strive to explain global collective phenomena by understanding the local interactions among the elements of the
collective. Through careful experimental design and statistical analysis, and by using numerical simulations to translate our ideas into a mathematical form, we have
caught a glimpse into the inner workings of collective phenomena, offering hope that
we may one day find universal principles. As Feynman remarked, “if we watch long
enough, we may eventually catch on to a few of the rules” [Feynman et al., 1963],
governing the remarkably diverse emergent systems we observe in nature.
154
Bibliography
[Adler, 1991] Adler, J. (1991). Bootstrap percolation. Physica A: Statistical Mechanics and its Applications, 171(3):453–470.
[Albert and Barabási, 2002] Albert, R. and Barabási, A.-L. (2002). Statistical mechanics of complex networks. Reviews of modern physics, 74(1):47.
[Bachmayer and Leonard, 2002] Bachmayer, R. and Leonard, N. E. (2002). Vehicle
networks for gradient descent in a sampled environment. In Decision and Control,
2002, Proceedings of the 41st IEEE Conference on, volume 1, pages 112–117. IEEE.
[Bak et al., 1987] Bak, P., Tang, C., and Wiesenfeld, K. (1987). Self-organized criticality: An explanation of the 1/f noise. Physical review letters, 59(4):381.
[Ballerini et al., 2008] Ballerini, M., Cabibbo, N., Candelier, R., Cavagna, A., Cisbani, E., Giardina, I., Lecomte, V., Orlandi, A., Parisi, G., Procaccini, A., et al.
(2008). Interaction ruling animal collective behavior depends on topological rather
than metric distance: Evidence from a field study. Proceedings of the national
academy of sciences, 105(4):1232–1237.
[Barthelemy, 2004] Barthelemy, M. (2004). Betweenness centrality in large complex
networks. The European Physical Journal B-Condensed Matter and Complex Systems, 38(2):163–168.
[Beauchamp, 2010] Beauchamp, G. (2010). Determinants of false alarms in staging
flocks of semipalmated sandpipers. Behavioral Ecology, 21(3):584–587.
[Beauchamp, 2012] Beauchamp, G. (2012). Flock size and density influence speed of
escape waves in semipalmated sandpipers. Animal Behaviour, 83(4):1125–1129.
[Beauchamp and Fernández-Juricic, 2005] Beauchamp, G. and Fernández-Juricic, E.
(2005). The group-size paradox: effects of learning and patch departure rules.
Behavioral Ecology, 16(2):352–357.
[Beggs and Plenz, 2003] Beggs, J. M. and Plenz, D. (2003). Neuronal avalanches in
neocortical circuits. The Journal of neuroscience, 23(35):11167–11177.
[Berdahl et al., 2013] Berdahl, A., Torney, C. J., Ioannou, C. C., Faria, J. J., and
Couzin, I. D. (2013). Emergent sensing of complex environments by mobile animal
groups. Science, 339(6119):574–576.
155
[Bertrand et al., 2008] Bertrand, A., Gerlotto, F., Bertrand, S., Gutiérrez, M., Alza,
L., Chipollini, A., Dı́az, E., Espinoza, P., Ledesma, J., Quesquén, R., et al.
(2008). Schooling behaviour and environmental forcing in relation to anchoveta
distribution: an analysis across multiple spatial scales. Progress in Oceanography,
79(2):264–277.
[Bertrand et al., 2014] Bertrand, A., Grados, D., Colas, F., Bertrand, S., Capet, X.,
Chaigneau, A., Vargas, G., Mousseigne, A., and Fablet, R. (2014). Broad impacts
of fine-scale dynamics on seascape structure from zooplankton to seabirds. Nature
communications, 5.
[Bhatkar and Whitcomb, 1970] Bhatkar, A. and Whitcomb, W. (1970). Artificial diet
for rearing various species of ants. Florida Entomologist, pages 229–232.
[Blackman and House, 1999] Blackman, S. and House, A. (1999). Design and analysis
of modern tracking systems. Boston, MA: Artech House.
[Blumstein et al., 2004] Blumstein, D. T., Verneyre, L., and Daniel, J. C. (2004). Reliability and the adaptive utility of discrimination among alarm callers. Proceedings
of the Royal Society of London. Series B: Biological Sciences, 271(1550):1851–1857.
[Boi et al., 1999] Boi, S., Couzin, I., Del Buono, N., Franks, N., and Britton, N.
(1999). Coupled oscillators and activity waves in ant colonies. Proceedings of the
Royal Society of London. Series B: Biological Sciences, 266(1417):371–378.
[Bonabeau et al., 1999] Bonabeau, E., Dagorn, L., and FREon, P. (1999). Scaling in
animal group-size distributions. Proceedings of the National Academy of Sciences,
96(8):4472–4477.
[Bonabeau et al., 1998] Bonabeau, E., Theraulaz, G., and Deneubourg, J.-L. (1998).
The synchronization of recruitment-based activities in ants. BioSystems, 45(3):195–
211.
[Bonabeau et al., 1997] Bonabeau, E., Theraulaz, G., Deneubourg, J.-L., Aron, S.,
and Camazine, S. (1997). Self-organization in social insects. Trends in Ecology &
Evolution, 12(5):188–193.
[Burgess and Shaw, 1981] Burgess, J. W. and Shaw, E. (1981). Effects of acousticolateralis denervation in a facultative schooling fish: A nearest-neighbor matrix
analysis. Behavioral and neural biology, 33(4):488–497.
[Buzsáki and Draguhn, 2004] Buzsáki, G. and Draguhn, A. (2004). Neuronal oscillations in cortical networks. science, 304(5679):1926–1929.
[Camperi et al., 2012] Camperi, M., Cavagna, A., Giardina, I., Parisi, G., and Silvestri, E. (2012). Spatially balanced topological interaction grants optimal cohesion
in flocking models. Interface focus, 2(6):715–725.
156
[Castellano et al., 2009] Castellano, C., Fortunato, S., and Loreto, V. (2009). Statistical physics of social dynamics. Reviews of modern physics, 81(2):591.
[Centola et al., 2007] Centola, D., Eguı́luz, V. M., and Macy, M. W. (2007). Cascade dynamics of complex propagation. Physica A: Statistical Mechanics and its
Applications, 374(1):449–456.
[Centola and Macy, 2007] Centola, D. and Macy, M. (2007). Complex contagions and
the weakness of long ties. American Journal of Sociology, 113(3):702–734.
[Chatterjee and Hadi, 2013] Chatterjee, S. and Hadi, A. S. (2013). Regression analysis by example. John Wiley & Sons.
[Chialvo, 2010] Chialvo, D. R. (2010). Emergent complex neural dynamics. Nature
physics, 6(10):744–750.
[Cohen et al., 2007] Cohen, J. D., McClure, S. M., and Angela, J. Y. (2007). Should i
stay or should i go? how the human brain manages the trade-off between exploitation and exploration. Philosophical Transactions of the Royal Society B: Biological
Sciences, 362(1481):933–942.
[Cole, 1991a] Cole, B. J. (1991a). Is animal behaviour chaotic? evidence from the
activity of ants. Proceedings of the Royal Society of London. Series B: Biological
Sciences, 244(1311):253–259.
[Cole, 1991b] Cole, B. J. (1991b). Short-term activity cycles in ants: generation of
periodicity by worker interaction. American Naturalist, pages 244–259.
[Cole, 1992] Cole, B. J. (1992). Short-term activity cycles in ants: age-related changes
in tempo and colony synchrony. Behavioral ecology and sociobiology, 31(3):181–187.
[Cole and Cheshire, 1996] Cole, B. J. and Cheshire, D. (1996). Mobile cellular automata models of ant behavior: movement activity of leptothorax allardycei. American Naturalist, pages 1–15.
[Cole and Hoeg, 1996] Cole, B. J. and Hoeg, L. (1996). The influence of brood type
on activity cycles in leptothorax allardycei (hymenoptera: Faruicidae). Journal of
insect Behavior, 9(4):539–547.
[Coombs, 1994] Coombs, S. (1994). Nearfield detection of dipole sources by the goldfish (carassius auratus) and the mottled sculpin (cottus bairdi). The Journal of
experimental biology, 190(1):109–129.
[Couzin, 2007] Couzin, I. (2007). Collective minds. Nature, 445(7129):715–715.
[Couzin, 2009] Couzin, I. D. (2009). Collective cognition in animal groups. Trends
in cognitive sciences, 13(1):36–43.
[Couzin and Krause, 2003] Couzin, I. D. and Krause, J. (2003). Self-organization and
collective behavior in vertebrates. Advances in the Study of Behavior, 32:1–75.
157
[Couzin et al., 2002] Couzin, I. D., Krause, J., James, R., Ruxton, G. D., and Franks,
N. R. (2002). Collective memory and spatial sorting in animal groups. Journal of
theoretical biology, 218(1):1–11.
[Couzin and Laidre, 2009] Couzin, I. D. and Laidre, M. E. (2009). Fission–fusion
populations. Current biology, 19(15):R633–R635.
[Cresswell et al., 2000] Cresswell, W., Hilton, G. M., and Ruxton, G. D. (2000). Evidence for a rule governing the avoidance of superfluous escape flights. Proceedings
of the Royal Society of London. Series B: Biological Sciences, 267(1444):733–737.
[Cucker and Smale, 2007] Cucker, F. and Smale, S. (2007). Emergent behavior in
flocks. Automatic Control, IEEE Transactions on, 52(5):852–862.
[Cvikel et al., 2015] Cvikel, N., Berg, K. E., Levin, E., Hurme, E., Borissov, I., Boonman, A., Amichai, E., and Yovel, Y. (2015). Bats aggregate to improve prey search
but might be impaired when their density becomes too high. Current Biology.
[Czirók et al., 1999] Czirók, A., Vicsek, M., and Vicsek, T. (1999). Collective motion of organisms in three dimensions. Physica A: Statistical Mechanics and its
Applications, 264(1):299–304.
[Dehaene, 2003] Dehaene, S. (2003). The neural basis of the weber–fechner law: a
logarithmic mental number line. Trends in cognitive sciences, 7(4):145–147.
[Delgado and Sole, 2000] Delgado, J. and Sole, R. V. (2000). Self-synchronization
and task fulfilment in ant colonies. Journal of theoretical biology, 205(3):433–441.
[Dill et al., 1997] Dill, L. M., Holling, C., and Palmer, L. H. (1997). Predicting the
three-dimensional structure of animal aggregations from functional considerations:
The role of information. New York: Cambridge University Press.
[Domenici and Batty, 1997] Domenici, P. and Batty, R. (1997). Escape behaviour
of solitary herring (clupea harengus) and comparisons with schooling individuals.
Marine Biology, 128(1):29–38.
[Domenici and Blake, 1997] Domenici, P. and Blake, R. (1997). The kinematics
and performance of fish fast-start swimming. Journal of Experimental Biology,
200(8):1165–1178.
[Dorogovtsev et al., 2006] Dorogovtsev, S. N., Goltsev, A. V., and Mendes, J.
F. F. (2006). K-core organization of complex networks. Physical review letters,
96(4):040601.
[D’Orsogna et al., 2006] D’Orsogna, M. R., Chuang, Y.-L., Bertozzi, A. L., and
Chayes, L. S. (2006). Self-propelled particles with soft-core interactions: patterns,
stability, and collapse. Physical review letters, 96(10):104302.
158
[Eaton et al., 2001] Eaton, R., Lee, R., and Foreman, M. (2001). The mauthner cell
and other identified neurons of the brainstem escape network of fish. Progress in
neurobiology, 63(4):467–485.
[Eaton, 1984] Eaton, R. C. (1984). Neural mechanisms of startle behavior. Springer
Science & Business Media.
[Fagiolo, 2007] Fagiolo, G. (2007). Clustering in complex directed networks. Physical
Review E, 76(2):026107.
[Feynman et al., 1963] Feynman, R. P., Leighton, R. B., and Sands, M. (1963). The
Feynman Lectures on Physics, Desktop Edition Volume I, volume 1. Basic Books.
[Fontaine et al., 2008] Fontaine, E., Lentink, D., Kranenbarg, S., Müller, U. K., van
Leeuwen, J. L., Barr, A. H., and Burdick, J. W. (2008). Automated visual tracking
for studying the ontogeny of zebrafish swimming. Journal of Experimental Biology,
211(8):1305–1316.
[Franks et al., 1990] Franks, N. R., Bryant, S., Griffiths, R., and Hemerik, L. (1990).
Synchronization of the behaviour within nests of the antleptothorax acervorum
(fabricius)i. discovering the phenomenon and its relation to the level of starvation.
Bulletin of Mathematical Biology, 52(5):597–612.
[Gerlotto et al., 2006] Gerlotto, F., Bertrand, S., Bez, N., and Gutierrez, M. (2006).
Waves of agitation inside anchovy schools observed with multibeam sonar: a way
to transmit information in response to predation. ICES Journal of Marine Science:
Journal du Conseil, 63(8):1405–1417.
[Gerstner and Kistler, 2002] Gerstner, W. and Kistler, W. M. (2002). Spiking neuron
models: Single neurons, populations, plasticity. Cambridge university press.
[Giraldeau et al., 2002] Giraldeau, L.-A., Valone, T. J., and Templeton, J. J. (2002).
Potential disadvantages of using socially acquired information. Philosophical
Transactions of the Royal Society of London. Series B: Biological Sciences,
357(1427):1559–1566.
[González-Bailón et al., 2011] González-Bailón, S., Borge-Holthoefer, J., Rivero, A.,
and Moreno, Y. (2011). The dynamics of protest recruitment through an online
network. Scientific reports, 1.
[Goss and Deneubourg, 1988] Goss, S. and Deneubourg, J.-L. (1988). Autocatalysis
as a source of synchronised rhythmical activity in social insects. Insectes Sociaux,
35(3):310–315.
[Guttal and Couzin, 2010] Guttal, V. and Couzin, I. D. (2010). Social interactions,
information use, and the evolution of collective migration. Proceedings of the National Academy of Sciences, 107(37):16172–16177.
159
[Handegard et al., 2012] Handegard, N. O., Boswell, K. M., Ioannou, C. C., Leblanc,
S. P., Tjøstheim, D. B., and Couzin, I. D. (2012). The dynamics of coordinated
group hunting and collective information transfer among schooling prey. Current
biology, 22(13):1213–1217.
[Hanke and Lauder, 2006] Hanke, W. and Lauder, G. (2006). Fish schooling: 3d
kinematics and hydrodynamics. In Integrative and Comparative Biology, volume 46,
pages E54–E54.
[Hatcher et al., 1992] Hatcher, M., Tofts, C., and Franks, N. (1992). Mutual exclusion
as a mechanism for information exchange within ant nests. Naturwissenschaften,
79(1):32–34.
[Hatta and Korn, 1999] Hatta, K. and Korn, H. (1999). Tonic inhibition alternates
in paired neurons that set direction of fish escape reaction. Proceedings of the
National Academy of Sciences, 96(21):12090–12095.
[Hein and McKinley, 2012] Hein, A. M. and McKinley, S. A. (2012). Sensing and
decision-making in random search. Proceedings of the National Academy of Sciences, 109(30):12070–12074.
[Herbert-Read et al., 2011] Herbert-Read, J. E., Perna, A., Mann, R. P., Schaerf,
T. M., Sumpter, D. J., and Ward, A. J. (2011). Inferring the rules of interaction
of shoaling fish. Proceedings of the National Academy of Sciences, 108(46):18726–
18731.
[Hopfield, 1982] Hopfield, J. J. (1982). Neural networks and physical systems with
emergent collective computational abilities. Proceedings of the national academy
of sciences, 79(8):2554–2558.
[Ising, 1925] Ising, E. (1925). Beitrag zur theorie des ferromagnetismus. Zeitschrift
für Physik A Hadrons and Nuclei, 31(1):253–258.
[Jaynes, 1957] Jaynes, E. T. (1957). Information theory and statistical mechanics.
Physical review, 106(4):620.
[Johannes et al., 1989] Johannes, M. R., McQueen, D. J., Stewart, T. J., and Post,
J. R. (1989). Golden shiner (notemigonus crysoleucas) population abundance correlations with food and predators. Canadian Journal of Fisheries and Aquatic
Sciences, 46(5):810–817.
[Johnson and Omland, 2004] Johnson, J. B. and Omland, K. S. (2004). Model selection in ecology and evolution. Trends in ecology & evolution, 19(2):101–108.
[Kahlert, 2006] Kahlert, J. (2006). Factors affecting escape behaviour in moulting
greylag geese anser anser. Journal of Ornithology, 147(4):569–577.
160
[Kai et al., 2009] Kai, E. T., Rossi, V., Sudre, J., Weimerskirch, H., Lopez, C.,
Hernandez-Garcia, E., Marsac, F., and Garçon, V. (2009). Top marine predators track lagrangian coherent structures. Proceedings of the National Academy of
Sciences, 106(20):8245–8250.
[Kao and Couzin, 2014] Kao, A. B. and Couzin, I. D. (2014). Decision accuracy in
complex environments is often maximized by small group sizes. Proceedings of the
Royal Society B: Biological Sciences, 281(1784):20133305.
[Kastberger et al., 2008] Kastberger, G., Schmelzer, E., and Kranner, I. (2008). Social waves in giant honeybees repel hornets. PLoS One, 3(9):e3141.
[Katz et al., 2011] Katz, Y., Tunstrøm, K., Ioannou, C. C., Huepe, C., and Couzin,
I. D. (2011). Inferring the structure and dynamics of interactions in schooling fish.
Proceedings of the National Academy of Sciences, 108(46):18720–18725.
[Keeling, 1999] Keeling, M. J. (1999). The effects of local spatial structure on epidemiological invasions. Proceedings of the Royal Society of London. Series B: Biological Sciences, 266(1421):859–867.
[Kinney et al., 2005] Kinney, R., Crucitti, P., Albert, R., and Latora, V. (2005). Modeling cascading failures in the north american power grid. The European Physical
Journal B-Condensed Matter and Complex Systems, 46(1):101–107.
[Kitsak et al., 2010] Kitsak, M., Gallos, L. K., Havlin, S., Liljeros, F., Muchnik, L.,
Stanley, H. E., and Makse, H. A. (2010). Identification of influential spreaders in
complex networks. Nature Physics, 6(11):888–893.
[Koch and Ullman, 1987] Koch, C. and Ullman, S. (1987). Shifts in selective visual
attention: towards the underlying neural circuitry. In Matters of intelligence, pages
115–141. Springer.
[Krause and Ruxton, 2002] Krause, J. and Ruxton, G. D. (2002). Living in groups.
OUP Oxford.
[Lekien and Leonard, 2010] Lekien, F. and Leonard, N. E. (2010). Nonuniform coverage and cartograms. In Decision and Control (CDC), 2010 49th IEEE Conference
on, pages 5518–5523. IEEE.
[Lenz, 1920] Lenz, W. (1920). Beitrag zum verständnis der magnetischen erscheinungen in festen körpern. Z. Phys., 21:613–615.
[Leskovec et al., 2009] Leskovec, J., Backstrom, L., and Kleinberg, J. (2009). Memetracking and the dynamics of the news cycle. In Proceedings of the 15th ACM
SIGKDD international conference on Knowledge discovery and data mining, pages
497–506. ACM.
[Lima, 1995] Lima, S. L. (1995). Collective detection of predatory attack by social
foragers: fraught with ambiguity? Animal Behaviour, 50(4):1097–1108.
161
[Lundgren, 2010] Lundgren, J. (2010). Alpha shape of 2d/3d point set. Matlab.
http://www.mathworks.com/matlabcentral/fileexchange/28851-alpha-shapes.
[Ma et al., 2011] Ma, Q., Johansson, A., and Sumpter, D. J. (2011). A first principles derivation of animal group size distributions. Journal of theoretical biology,
283(1):35–43.
[Mann et al., 2013] Mann, R. P., Perna, A., Strömbom, D., Garnett, R., HerbertRead, J. E., Sumpter, D. J., and Ward, A. J. (2013). Multi-scale inference of
interaction rules in animal groups using bayesian model selection. PLoS computational biology, 9(3):e1002961.
[Marras and Domenici, 2013] Marras, S. and Domenici, P. (2013). Schooling fish
under attack are not all equal: Some lead, others follow. PloS one, 8(6):e65784.
[Mason and Watts, 2012] Mason, W. and Watts, D. J. (2012). Collaborative learning
in networks. Proceedings of the National Academy of Sciences, 109(3):764–769.
[Mersch et al., 2013] Mersch, D. P., Crespi, A., and Keller, L. (2013). Tracking individuals shows spatial fidelity is a key regulator of ant social organization. Science,
340(6136):1090–1093.
[Meyers et al., 2006] Meyers, L. S., Gamst, G., and Guarino, A. (2006). Applied
multivariate research: Design and interpretation. Sage.
[Miller and Bassler, 2001] Miller, M. B. and Bassler, B. L. (2001). Quorum sensing
in bacteria. Annual Reviews in Microbiology, 55(1):165–199.
[Miramontes et al., 1993] Miramontes, O., Solé, R. V., and Goodwin, B. C. (1993).
Collective behaviour of random-activated mobile cellular automata. Physica D:
Nonlinear Phenomena, 63(1):145–160.
[Mirollo and Strogatz, 1990] Mirollo, R. E. and Strogatz, S. H. (1990). Synchronization of pulse-coupled biological oscillators. SIAM Journal on Applied Mathematics,
50(6):1645–1662.
[Mogilner et al., 2003] Mogilner, A., Edelstein-Keshet, L., Bent, L., and Spiros, A.
(2003). Mutual interactions, potentials, and individual distance in a social aggregation. Journal of mathematical biology, 47(4):353–389.
[Moore and Newman, 2000] Moore, C. and Newman, M. E. (2000). Epidemics and
percolation in small-world networks. Physical Review E, 61(5):5678.
[Mora and Bialek, 2011] Mora, T. and Bialek, W. (2011). Are biological systems
poised at criticality? Journal of Statistical Physics, 144(2):268–302.
[Motter and Lai, 2002] Motter, A. E. and Lai, Y.-C. (2002). Cascade-based attacks
on complex networks. Physical Review E, 66(6):065102.
162
[Niwa, 2004] Niwa, H.-S. (2004). Space-irrelevant scaling law for fish school sizes.
Journal of theoretical biology, 228(3):347–357.
[Onsager, 1944] Onsager, L. (1944). Crystal statistics. i. a two-dimensional model
with an order-disorder transition. Physical Review, 65(3-4):117.
[Park and Hastie, 2007] Park, M. Y. and Hastie, T. (2007). L1-regularization path
algorithm for generalized linear models. Journal of the Royal Statistical Society:
Series B (Statistical Methodology), 69(4):659–677.
[Partridge, 1982] Partridge, B. L. (1982). The structure and function of fish schools.
Scientific american, 246(6):114–123.
[Partridge et al., 1980] Partridge, B. L., Pitcher, T., Cullen, J. M., and Wilson, J.
(1980). The three-dimensional structure of fish schools. Behavioral Ecology and
Sociobiology, 6(4):277–288.
[Passino, 2002] Passino, K. M. (2002). Biomimicry of bacterial foraging for distributed optimization and control. Control Systems, IEEE, 22(3):52–67.
[Pei et al., 2014] Pei, S., Muchnik, L., Andrade Jr, J. S., Zheng, Z., and Makse, H. A.
(2014). Searching for superspreaders of information in real-world social media.
Scientific reports, 4.
[Pikovsky et al., 2002] Pikovsky, A., Rosenblum, M., Kurths, J., and Hilborn, R. C.
(2002). Synchronization: A universal concept in nonlinear sciences, volume 2.
Cambridge University Press Cambridge.
[Procaccini et al., 2011] Procaccini, A., Orlandi, A., Cavagna, A., Giardina, I., Zoratto, F., Santucci, D., Chiarotti, F., Hemelrijk, C. K., Alleva, E., Parisi, G., et al.
(2011). Propagating waves in starling, sturnus vulgaris, flocks under predation.
Animal behaviour, 82(4):759–765.
[Rabinovich et al., 2006] Rabinovich, M. I., Varona, P., Selverston, A. I., and Abarbanel, H. D. (2006). Dynamical principles in neuroscience. Reviews of modern
physics, 78(4):1213.
[Radakov, 1973] Radakov, D. V. (1973). Schooling in the ecology of fish.
[Robinson, 1992] Robinson, G. E. (1992). Regulation of division of labor in insect
societies. Annual review of entomology, 37(1):637–665.
[Rosenthal et al., 2015] Rosenthal, S. B., Twomey, C. R., Hartnett, A. T., Wu, H. S.,
and Couzin, I. D. (2015). Revealing the hidden networks of interaction in mobile
animal groups allows prediction of complex behavioral contagion. Proceedings of
the National Academy of Sciences, in press.
[Sachtjen et al., 2000] Sachtjen, M., Carreras, B., and Lynch, V. (2000). Disturbances
in a power transmission system. Physical Review E, 61(5):4877.
163
[Selous, 1931] Selous, E. (1931). Thought-Transference-or What?-in Birds. Constable
& Company.
[Sendova-Franks and Franks, 1995] Sendova-Franks, A. and Franks, N. (1995). Spatial relationships within nests of the ant leptothorax unifasciatus and their implications for the division of labour. Animal Behaviour, 50(1):121–136.
[Sibly, 1983] Sibly, R. (1983). Optimal group size is unstable. Animal Behaviour,
31(3):947–948.
[Slobodchikoff, 2013] Slobodchikoff, C. N. (2013). The ecology of social behavior.
Elsevier.
[Solé et al., 1993a] Solé, R., Miramontes, O., and Goodwin, B. (1993a). Emergent
behavior in insect societies: Global oscillations, chaos and computation. In Interdisciplinary approaches to nonlinear complex systems, pages 77–88. Springer.
[Solé et al., 1993b] Solé, R. V., Miramontes, O., and Goodwin, B. C. (1993b). Oscillations and chaos in ant societies. Journal of Theoretical Biology, 161(3):343–357.
[Solomon et al., 2000] Solomon, S., Weisbuch, G., de Arcangelis, L., Jan, N., and
Stauffer, D. (2000). Social percolation models. Physica A: Statistical Mechanics
and its Applications, 277(1):239–247.
[Strandburg-Peshkin et al., 2013] Strandburg-Peshkin, A., Twomey, C. R., Bode,
N. W., Kao, A. B., Katz, Y., Ioannou, C. C., Rosenthal, S. B., Torney, C. J.,
Wu, H. S., Levin, S. A., et al. (2013). Visual sensory networks and effective information transfer in animal groups. Current Biology, 23(17):R709–R711.
[Strogatz, 2003] Strogatz, S. (2003). Sync: The emerging science of spontaneous
order. Hyperion.
[Sumpter, 2010] Sumpter, D. J. (2010). Collective animal behavior. Princeton University Press.
[Sutton and Barto, 1998] Sutton, R. S. and Barto, A. G. (1998). Introduction to
reinforcement learning. MIT Press.
[Tkacik et al., 2006] Tkacik, G., Schneidman, E., Berry, I., Michael, J., and Bialek,
W. (2006). Ising models for networks of real neurons. arXiv preprint q-bio/0611072.
[Torney et al., 2009] Torney, C., Neufeld, Z., and Couzin, I. D. (2009). Contextdependent interaction leads to emergent search behavior in social aggregates. Proceedings of the National Academy of Sciences, 106(52):22055–22060.
[Treherne and Foster, 1981] Treherne, J. and Foster, W. (1981). Group transmission
of predator avoidance behaviour in a marine insect: the trafalgar effect. Animal
Behaviour, 29(3):911–917.
164
[Tunstrøm et al., 2013] Tunstrøm, K., Katz, Y., Ioannou, C. C., Huepe, C., Lutz,
M. J., and Couzin, I. D. (2013). Collective states, multistability and transitional
behavior in schooling fish. PLoS computational biology, 9(2):e1002915.
[Vabø and Nøttestad, 1997] Vabø, R. and Nøttestad, L. (1997). An individual based
model of fish school reactions: predicting antipredator behaviour as observed in
nature. Fisheries oceanography, 6(3):155–171.
[Vicsek et al., 1995] Vicsek, T., Czirók, A., Ben-Jacob, E., Cohen, I., and Shochet, O.
(1995). Novel type of phase transition in a system of self-driven particles. Physical
review letters, 75(6):1226.
[Vicsek and Zafeiris, 2012] Vicsek, T. and Zafeiris, A. (2012). Collective motion.
Physics Reports, 517(3):71–140.
[Vlaev et al., 2011] Vlaev, I., Chater, N., Stewart, N., and Brown, G. D. (2011). Does
the brain calculate value? Trends in cognitive sciences, 15(11):546–554.
[Volz et al., 2011] Volz, E. M., Miller, J. C., Galvani, A., and Meyers, L. A. (2011).
Effects of heterogeneous and clustered contact patterns on infectious disease dynamics. PLoS computational biology, 7(6):e1002042.
[Watson, 1981] Watson, D. F. (1981). Computing the n-dimensional delaunay tessellation with application to voronoi polytopes. The computer journal, 24(2):167–172.
[Watts, 2002] Watts, D. J. (2002). A simple model of global cascades on random
networks. Proceedings of the National Academy of Sciences, 99(9):5766–5771.
[Whittier et al., 2000] Whittier, T. R., Halliwell, D. B., and Daniels, R. A. (2000).
Distributions of lake fishes in the northeast-ii. the minnows (cyprinidae). Northeastern Naturalist, 7(2):131–156.
[Wood and Ackland, 2007] Wood, A. J. and Ackland, G. J. (2007). Evolving the
selfish herd: emergence of distinct aggregating strategies in an individual-based
model. Proceedings of the Royal Society B: Biological Sciences, 274(1618):1637–
1642.
[Wu et al., 2011] Wu, H. S., Zhao, Q., Zou, D., and Chen, Y. Q. (2011). Automated
3d trajectory measuring of large numbers of moving particles. Optics express,
19(8):7646–7663.
[Young, 2011] Young, H. P. (2011). The dynamics of social innovation. Proceedings
of the National Academy of Sciences, 108(Supplement 4):21285–21291.
165
Fly UP