Modeling Collective Motion in Animal Groups: from Mathematical Analysis to Field Data by Ryan J. Lukeman B.A., St. Francis Xavier University, 2003 M.Sc., Dalhousie University, 2005 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Mathematics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2009 c© Ryan J. Lukeman 2009 Abstract Animals moving together cohesively is a commonly observed phenomenon in biology, with bird flocks and fish schools as familiar examples. Mathematical models have been developed in order to understand the mechanisms that lead to such coordinated motion. The Lagrangian framework of modeling, wherein individuals within the group are modeled as point particles with po- sition and velocity, permits construction of inter-individual interactions via ‘social forces’ of attraction, repulsion and alignment. Although such models have been studied extensively via numerical simulation, analytical conclu- sions have been difficult to obtain, owing to the large size of the associated system of differential equations. In this thesis, I contribute to the modeling of collective motion in two ways. First, I develop a simplified model of mo- tion and, by focusing on simple, regular solutions, am able to connect group properties to individual characteristics in a concrete manner via derivations of existence and stability conditions for a number of solution types. I show that existence of particular solutions depends on the attraction-repulsion function, while stability depends on the derivative of this function. Second, to establish validity and motivate construction of specific models for collective motion, actual data is required. I describe work gathering and analyzing dynamic data on group motion of surf scoters, a type of diving duck. This data represents, to our knowledge, the largest animal group size (by almost an order of magnitude) for which the trajectory of each group member is reconstructed. By constructing spatial distributions of neighbour density and mean deviation, I show that frontal neighbour preference and angular deviation are important features in such groups. I show that the observed spatial distribution of neighbors can be obtained in a model incor- porating a topological frontal interaction, and I find an optimal parameter set to match simulated data to empirical data. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii Statement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xviii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Biological Background . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Mathematical Modeling . . . . . . . . . . . . . . . . . . . . . 2 1.3 Survey of the Literature : Models . . . . . . . . . . . . . . . . 2 1.3.1 Eulerian Models . . . . . . . . . . . . . . . . . . . . . 2 1.3.2 Lagrangian Models . . . . . . . . . . . . . . . . . . . . 3 1.3.3 Unifying Lagrangian Approaches . . . . . . . . . . . . 7 1.4 The Need for Empirical Data . . . . . . . . . . . . . . . . . . 7 1.5 Survey of the Literature: Empirical Studies . . . . . . . . . . 8 1.5.1 Small Groups, No Tracking . . . . . . . . . . . . . . . 8 1.5.2 Small Groups with Tracking . . . . . . . . . . . . . . . 9 1.5.3 Large Groups, No Tracking . . . . . . . . . . . . . . . 10 1.5.4 What Empirical Data is Needed . . . . . . . . . . . . 11 1.6 Contributions of this Thesis . . . . . . . . . . . . . . . . . . . 11 1.6.1 Model Analysis . . . . . . . . . . . . . . . . . . . . . . 12 1.6.2 Empirical Data . . . . . . . . . . . . . . . . . . . . . . 13 iii Table of Contents Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2 Minimal Mechanisms for School Formation in Self-Propelled Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.2 The Model of Interacting Self-Propelled Particles . . . . . . . 24 2.2.1 Assumptions and Equations . . . . . . . . . . . . . . . 24 2.2.2 Classification of Schools . . . . . . . . . . . . . . . . . 27 2.3 Schools in One Dimensional Space . . . . . . . . . . . . . . . 29 2.3.1 Schools Formed by Following One Immediate Neigh- bour . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.2 Schools Formed by Interactions with Two Nearest Neigh- bours. . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.3.3 Examples of Some Schooling Forces . . . . . . . . . . 34 2.3.4 The Stability of a School Solution . . . . . . . . . . . 36 2.3.5 Numerical Simulations in One Dimension . . . . . . . 41 2.4 The Soldier Formation in 2D Space . . . . . . . . . . . . . . . 43 2.4.1 Existence Conditions for the Soldier Formation . . . . 45 2.4.2 Stability Analysis . . . . . . . . . . . . . . . . . . . . 46 2.4.3 Numerical Simulation of the Soldier Formation . . . . 50 2.5 Schools in the Form of 2D Arrays . . . . . . . . . . . . . . . . 53 2.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3 A Conceptual Model for Milling Formations in Biological Aggregates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.2 The Model of Interacting Self-Propelled Particles . . . . . . . 66 3.2.1 Relationship to Previous Models . . . . . . . . . . . . 67 3.2.2 The Milling Formation . . . . . . . . . . . . . . . . . . 68 3.3 Existence Conditions . . . . . . . . . . . . . . . . . . . . . . . 69 3.4 Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . 72 3.4.1 Eigenvalues . . . . . . . . . . . . . . . . . . . . . . . . 76 3.5 Numerical Investigation . . . . . . . . . . . . . . . . . . . . . 77 3.5.1 Numerical Simulations . . . . . . . . . . . . . . . . . . 78 3.5.2 Moving Mill Formations . . . . . . . . . . . . . . . . . 84 3.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 iv Table of Contents Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4 A Field Study of Collective Behaviour in Surf Scoters: Em- pirical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.1.1 Difficulties in Obtaining Data . . . . . . . . . . . . . . 97 4.1.2 Previous Work . . . . . . . . . . . . . . . . . . . . . . 98 4.1.3 Tracking . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.1.4 Challenges Avoided in This Study . . . . . . . . . . . 99 4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.2.1 Location and Materials . . . . . . . . . . . . . . . . . 99 4.2.2 Surf Scoter Behaviour in Field Study . . . . . . . . . . 100 4.2.3 Experimental Technique . . . . . . . . . . . . . . . . . 100 4.2.4 Calibration and Testing . . . . . . . . . . . . . . . . . 100 4.2.5 Postprocessing Images . . . . . . . . . . . . . . . . . . 107 4.2.6 Extracting Positions . . . . . . . . . . . . . . . . . . . 107 4.2.7 Tracking . . . . . . . . . . . . . . . . . . . . . . . . . 108 4.2.8 Edge Effects . . . . . . . . . . . . . . . . . . . . . . . 110 4.2.9 Processed Events . . . . . . . . . . . . . . . . . . . . . 113 4.2.10 Body Alignment Versus Velocity . . . . . . . . . . . . 114 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 4.3.1 Basic Statistics . . . . . . . . . . . . . . . . . . . . . . 116 4.3.2 Nearest-Neighbour Distance Distributions . . . . . . . 117 4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5 Ducks in a Row: Inferring Interaction Mechanisms from Field Data of Surf Scoters . . . . . . . . . . . . . . . . . . . . . 124 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.1.1 Specific Goals of this Work . . . . . . . . . . . . . . . 125 5.1.2 Using Spatial Distributions . . . . . . . . . . . . . . . 125 5.2 Empirical Methods . . . . . . . . . . . . . . . . . . . . . . . . 125 5.3 Empirical Results . . . . . . . . . . . . . . . . . . . . . . . . 126 5.3.1 Relative Location of Neighbors . . . . . . . . . . . . . 126 5.3.2 Relative Deviation of Neighbors . . . . . . . . . . . . 129 5.4 Building a Model . . . . . . . . . . . . . . . . . . . . . . . . . 131 5.4.1 Model Framework . . . . . . . . . . . . . . . . . . . . 133 v Table of Contents 5.5 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 5.5.1 Fixed Parameters . . . . . . . . . . . . . . . . . . . . 141 5.5.2 Free Parameters . . . . . . . . . . . . . . . . . . . . . 143 5.5.3 Parameter Effects on Radial Distributions . . . . . . . 143 5.6 An Optimal Parameter Set . . . . . . . . . . . . . . . . . . . 146 5.6.1 Goodness-of-fit Measure . . . . . . . . . . . . . . . . . 146 5.6.2 Optimization Process . . . . . . . . . . . . . . . . . . 146 5.6.3 Optimization Results . . . . . . . . . . . . . . . . . . 147 5.6.4 Detailed Parameter Exploration Near Optimal Set . . 147 5.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 6.2 Analysis of Perfect School Solutions . . . . . . . . . . . . . . 157 6.2.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . 157 6.2.2 Perfect School Analysis in Context . . . . . . . . . . . 158 6.3 Empirical Data and Modeling . . . . . . . . . . . . . . . . . . 159 6.3.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . 159 6.3.2 Empirical Data and Modeling in Context . . . . . . . 161 6.4 Summary: Main Contributions . . . . . . . . . . . . . . . . . 162 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 Appendices A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 A.1 The Stability Analysis in 1D . . . . . . . . . . . . . . . . . . 166 A.2 Stability Analysis for the Soldier Formation in 2D . . . . . . 168 B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 B.1 Transformation to Relative Coordinates . . . . . . . . . . . . 171 B.2 Derivation of the Linearized Perturbed System . . . . . . . . 174 B.3 Eigenvalue Equation . . . . . . . . . . . . . . . . . . . . . . . 176 B.4 Derivation of Inequality (3.37) . . . . . . . . . . . . . . . . . 177 vi Table of Contents C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 C.1 Data Sequences: Details . . . . . . . . . . . . . . . . . . . . . 179 C.2 Data Analysis: Details on Currents . . . . . . . . . . . . . . . 179 vii List of Tables 4.1 Basic statistics: data sequences . . . . . . . . . . . . . . . . . 117 4.2 Average speed and nearest neighbor distances for the 14 sep- arate sequences of data that were analyzed. Units are body- lengths (BL). Speed and NND was averaged over all individ- uals in all frames. The last column is the standard deviation of NND. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.1 Model : summary of interaction forces . . . . . . . . . . . . . 142 5.2 Summary of parameters: Using Model V together with an op- timization routine, an optimal parameter set is found that best matches simulated neighbor distributions to observed neigh- bor distributions. Parameters are dimensionless unless units are given. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 C.1 Data sequences: details . . . . . . . . . . . . . . . . . . . . . . 180 C.2 Current statistics: by sequence . . . . . . . . . . . . . . . . . . 181 viii List of Figures 2.1 A schematic diagram of model vectors in 2 dimensions. . . . . 26 2.2 A schematic diagram of particles in 1 dimension, showing our subscripting convention. Grey arrows indicate distance- dependent interaction forces. . . . . . . . . . . . . . . . . . . . 29 2.3 An example of a distant-dependent force in the form of a Hill function with a decay past xf , to depict a finite sensing range. 35 2.4 Four cases of neighbour interaction: [a] neighbour j travels faster than i in the same direction, (vji > 0), [b] neighbour j travels slower than i in the same direction, (vji < 0), [c] neigh- bour j diverging from i, (vji > 0), [d] neighbour j converging to i, (vji < 0). . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.5 [a] An example of interaction forces where the net force has positive slope at d. Here, g+(x) = 2[3x 2/(x2 + 42) − 1], and g−(x) = 3x/(x + 4) − 1. [b] An example of interaction forces where the net force has negative slope at d. Here, g+(x) = (1/2)[3x3/(x3 + 23)− 1], and g−(x) = 3x/(x+ 8)− 1. . . . . . 37 2.6 1D simulations showing position in a moving frame moving at velocity v as indicated. Interaction forces are as in [a] Fig. 2.5.a, and [b] Fig. 2.5.b. Self propulsion terms are [a.i] al = 0.5, [a.ii] al = 0.7, [b.i] al = 0.9, [b.ii] al = 0.8, while a = 0.1 in all cases. In [a], higher school speed resulted in an increase in NND from 3.26 in [a.i] to 3.49 in [a.ii]. The opposite occurred in [b], where higher school speed caused a decrease in d from 3.73 in [b.i] to 3.175 in [b.ii]. . . . . . . . . . . . . . . . . . . . 41 2.7 An example of an approximate soldier formation in Atlantic bluefin tuna, courtesy of Dr. M. Lutcavage, Large Pelagics Research Center, UNH. . . . . . . . . . . . . . . . . . . . . . . 43 ix List of Figures 2.8 A schematic diagram of individuals in soldier formation with frontal nearest-neighbour detection. Grey arrows indicate the direction of schooling force, while black arrows indicate direc- tion of motion. . . . . . . . . . . . . . . . . . . . . . . . . . . 44 2.9 Numerical simulations in 2D for [a] 8 particles with ~al = [0.1, 0.1], ~a = [0,−0.2] and [b] 6 particles with ~al = [0.1, 0.1], ~a = [0, 0]. Trajectories are plotted through time. Note in [a] that the solution line connecting individuals is oriented in the direction of ~al − ~a = [0.1, 0.3], while the school velocity at steady state is ~v = [0.2, 0.2]. In [b], both the solution line and school velocity are in the same direction, with ~al = [0.1, 0.1], ~v = [0.2, 0.2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 2.10 Numerical simulations in 2D, for [a] 10 particles and [b] 25 par- ticles with ~al = ~a = [0.1, 0.1], and g+(x) = 0.5(exp(−x/50)− 2 exp(−x/1)). Trajectories are plotted through time, and grey arrows indicate interactions between nearest neighbours. Note that individuals have constant spacing between nearest neigh- bours, but are not restricted to a particular relative angle with nearest neighbours. . . . . . . . . . . . . . . . . . . . . . . . . 52 2.11 Numerical simulations in 2D with g+(x) = 1−x/2, which cor- responds to (rather non-physical) short-range attraction and long-range repulsion. Trajectories are plotted in time for [a] 4 particles evolved to t = 40, and [b] 6 particles evolved to t = 300. Individuals group in pairs (though due to overlap- ping, some pairs appear as a single particle), and these pairs follow unpredictable trajectories. . . . . . . . . . . . . . . . . 52 2.12 Numerical simulation of two nearest-neighbour interactions. Trajectories are plotted through time. Note that the particles, initially perturbed, evolve to a regular tiling of the plane with individual distance d = 0.707. . . . . . . . . . . . . . . . . . . 55 3.1 An example of milling in Atlantic bluefin tuna, courtesy of Dr. M. Lutcavage, Large Pelagics Research Center, UNH. . . . . . 65 3.2 A schematic diagram of the milling formation. Black arrows indicate direction of motion (tangential to the circle), while grey arrows indicate direction of schooling force. . . . . . . . . 69 3.3 A schematic diagram showing the decomposition of the inter- action force into tangential and centripetal components. . . . 70 x List of Figures 3.4 A schematic diagram of particles in the mill formation show- ing angles introduced in the text, and relative position and velocity vectors. . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.5 Equation (3.17) represented graphically for γ = 0.5, n = 5, and g(x) = A exp(−x/a) − B exp(−x/b). Horizontal axis: inter-individual distance d, vertical axis: distance-dependent force magnitude. A = 1.5, a = 10, B = 3. In panel [a], b = 1.5, and milling is possible whereas in panel [b], b = 2 and no milling solution exists. . . . . . . . . . . . . . . . . . . . . 72 3.6 A plot of the left-hand and right-hand sides of (3.17) for a range of n values, from n = 4 to n = 12. Note that no intersection occurs for n = 4, but two occur for n = 5. . . . . 73 3.7 Numerical solution to (3.34) over a range of values of g′(d) for four particles (n = 4) and γ = 0.5. The real parts of eigenvalues are plotted for i = 0, . . . , 3. The asterisks indicate branches of Re(λ) with multiplicity equal to 2. Note that the region over which all eigenvalues have non-positive real parts is approximately [-0.26,0.25]. . . . . . . . . . . . . . . . . . . . 79 3.8 As in Figure 3.7, but for five particles (n = 5). The real part of eigenvalues are plotted for i = 0, . . . , 4. The region over which all eigenvalues have non-positive real parts is now approximately [-0.11,0.19]. . . . . . . . . . . . . . . . . . . . 80 3.9 As in Figure 3.7, but for six particles (n = 6). The real part of eigenvalues are plotted for i = 0, . . . , 5. The region over which all eigenvalues have non-positive real parts is now ap- proximately [-0.08,0.17]. . . . . . . . . . . . . . . . . . . . . . 81 3.10 A plot of the condition for mill solutions to exist [a1] or not [b1] given by Equation (3.17), and corresponding simulated particle tracks [a2] and [b2], with γ = 0.5. Note that in the case of an intersection, a stable mill forms at d ≈ 2 (whereas d ≈ 0.5 is unstable). However, when no intersection exists, the mill stops rotating and stationary particles are then arranged with d such that g(d) = 0 (final direction is outward due to small repulsive forces felt just before the particles stop). Arrows indicate direction facing at the end of the simulation. The interaction function used is as in Figure 3.5 with A = 0.5, a = 5, B = 1, and [a] b = 0.5, [b], b = 0.8. . . . . . . . . . . . 83 xi List of Figures 3.11 [a1] An interaction function with g′(d) lying outside the sta- bility region for both intersections (implying no stable mill formation). Simulations using this function are shown in [a2]. Note that particles initially form a mill-like solution which is destroyed as time evolves. [b1] An interaction function with only one intersection, whose slope is too large (g′(x) > s for all x). Simulations in [b2] show that the radius of the mill increases in time. . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.12 The irregular periodic mill formation. [a] Individual distance d versus time for one particle in a system of 5 particles. For the 4 times denoted in [a], snapshots of particles are shown in [b], showing the variation in d between particles 1 and 2. The interaction function used was g(x) = 1 − 0.11x, so that g′(d) = −0.11 lies on the stability boundary. . . . . . . . . . . 86 3.13 A plot of trajectories in time for six particles with~a = (0.1, 0.1)T . Note that the entire mill moves in the direction of ~a. . . . . . 88 3.14 A plot of trajectories in time for five particles with [a] ~a = (0.223, 0.446)T , and [b] ~a = (0.224, 0.448)T . Note that for a small parameter change, the system behaviour is fundamen- tally different. In these simulations, g(x) is as in Figure 3.5, with A = 1.5, a = 10, B = 3, and b = 1.6. The corresponding existence condition is shown in the inset to [b]. . . . . . . . . . 89 4.1 A schematic diagram of the experimental setup. θ is the cam- era axis angle, while φ is the angle-of-view of the camera lens. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.2 A schematic diagram of the setup with angles used in vertical transformation. φ̂ represents the angle corresponding to real distance y (in pixels). L is the real vertical extent of the image (in pixels). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 4.3 The horizontal calibration matrix, associating each (x, y) im- age coordinate with the horizontal real pixel value. . . . . . . 104 4.4 The calibration grid used to test the image transformation. The pixel location of the upper-right corner of each grid vertex was marked and transformed in MATLAB (see Fig. 4.5). . . . 105 4.5 42 Grid points marked in Fig. 4.4 plotted as ‘x’ marks, with reconstructed positions as ‘o’ marks. . . . . . . . . . . . . . . 106 xii List of Figures 4.6 An example of one image from a time series of images collected in the field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4.7 Processing an image to obtain positions. . . . . . . . . . . . . 109 4.8 An example frame with reconstructed velocities (grey), show- ing partially incorrect and missing velocities. In reality, indi- viduals are highly polarized in this frame. . . . . . . . . . . . 111 4.9 The example frame of Fig. 4.8 showing velocities (grey) ob- tained by re-tracking the event. . . . . . . . . . . . . . . . . . 111 4.10 Example trajectories for 4 groups. Starting positions are plot- ted in green, final positions in red. . . . . . . . . . . . . . . . . 112 4.11 In the first frame, the individual at the origin has nearest neighbours filling each quadrant, and so is not an edge indi- vidual. In the second frame, the individual has no nearest neighbours (of the first 8) in the third quadrant, and thus is considered an edge individual. . . . . . . . . . . . . . . . . . . 113 4.12 A plot of a group with edge individuals (as determined by the algorithm outlined in the text) in red, with all other individ- uals in black. . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 4.13 Current vector x-value cx vs. waste tracer drift speed, from raw images, with linear least-squares fit. . . . . . . . . . . . . 116 4.14 Nearest-neighbour distributions for the first 8 nearest neigh- bours. Successive means are also plotted (◦) above distributions.119 4.15 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 4.16 Nearest neighbour distributions of the first neighbour, with fitted probability density function q(d) (dashed) overlaid. . . . 120 5.1 A typical analyzed frame: input image (a) is processed, trans- formed to real space, and individual positions are reconstructed (b). Positions in successive frames are linked using particle- tracking software, giving velocities (b). . . . . . . . . . . . . . 127 5.2 Trajectories for 4 sequences, showing surf scoters approaching the dock (y = 0) in a highly polarized fashion. Starting po- sitions are plotted in green, final positions in red. Units are given in pixels, where 1 BL = 46 pixels. . . . . . . . . . . . . . 128 xiii List of Figures 5.3 Neighbour distributions from empirical data, normalized to have maximal value 1. Neighbor positions are calculated rela- tive to the heading of a focal individual (in the direction [0 1] in the plot, indicated by the white graphic). A frontal bias is seen in neighbor positioning. . . . . . . . . . . . . . . . . . . . 130 5.4 Average density of neighbours in a circular shell at the pre- ferred distance, as a function of angle. A distinct region of higher density is centered at 90◦. . . . . . . . . . . . . . . . . 131 5.5 Spatial distribution of heading deviation, from empirical data. Radial bands are plotted at 1 BL, 2 BL, and 3 BL for reference. The region where deviation is 0 corresponds to the repulsion zone where no individuals are found. . . . . . . . . . . . . . . 132 5.6 A schematic diagram partitioning local space around a focal individual into sectors of location preference, high tendency to align, and high tendency to deviate rapidly (move away from) the focal individual. Schematic was assembled according to trends shown in the data in Figs. 5.3 and 5.5. . . . . . . . . . 133 5.7 Attraction-repulsion function g(x) is negative in repulsion re- gion R, zero in alignment region AL, and positive in attraction region ATT. Magnitudes of attraction and repulsion are con- trolled by weighting parameters ωatt and ωrep. . . . . . . . . . 135 5.8 Model I: (a) A schematic diagram representing the interaction zones, here being simply repulsion. (b) Spatial distributions of neighbors, relative to a focal individual oriented in the di- rection [0 1] (as indicated by the white graphic at the origin), where density has been normalized to have maximum value equal to 1. White dashed lines are superimposed at radial dis- tances 1, 2 and 3 BL. (c) Angular density distributions at the preferred distance are plotted, following the convention of Fig. 5.4. I compare distributions for data (green) and the model (blue). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 5.9 As in Fig. 5.8, but for Model II. . . . . . . . . . . . . . . . . . 138 5.10 As in Fig. 5.8, but for Model III. . . . . . . . . . . . . . . . . 139 5.11 As in Fig. 5.8, but for Model IV. . . . . . . . . . . . . . . . . 140 5.12 Attraction-repulsion function gf(x) is negative in repulsion re- gion R, and positive and constant beyond R (up to a limit of ratt). The magnitude of gf(x) is controlled by the weighting parameter ωfront. . . . . . . . . . . . . . . . . . . . . . . . . . 141 xiv List of Figures 5.13 Average radial neighbor density, from data. . . . . . . . . . . . 144 5.14 Effect of parameter variation on radial neighbor distribution in basic model simulations. Unless varied as indicated, param- eters are ωrep = 5, ωatt = 1, ωal = 1, ωξ = 0.15. 50 Simulations of 100 individuals were performed to t = 100, and average den- sities over all simulations (after quasi-equilibrium was reached, t = 50 to t = 100) were calculated. . . . . . . . . . . . . . . . 145 5.15 A comparison of spatial neighbor density for [a] the best-fit simulations of the most appropriate model, Model V , and [b] data (repeating Fig. 5.3). Essential features of the data are observed in simulated data, including the spatial extent of neighbors, and the frontal bias. . . . . . . . . . . . . . . . . . 148 5.16 A comparison of average angular neighbor density at the pre- ferred distance for data (green), and simulation of Model V with optimal parameters (blue). . . . . . . . . . . . . . . . . . 149 5.17 Error as each parameter is varied about the optimal value. Each parameter is scaled in the plot so that optimal val- ues coincide (blue dotted line). Parameter values explored were as follows: ωrep = 5, 7.5, 10, 12.5, 15, ωatt = 0.5, 1, 1.5, 2, ωal = 0.25, 0.5, 0.75, 1, ωfront = 0.05, 0.1, 0.15, 0.2, and ωξ = 0.275, 0.3, 0.325, 0.35, 0.375. In each case, the optimal value is a minimum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 5.18 An example of a follow-the-leader formation observed in surf scoters approaching a dock to forage on mussels. Here, inter- actions to the front are dominant. . . . . . . . . . . . . . . . . 152 B.1 A schematic diagram of particles at the threshold of mill break- ing. The angle between particles is π/2, beyond which particle i no longer senses particle i + 1. ~vmi indicates the velocity of particle i in the absence of autonomous self-propulsion. . . . 177 xv Acknowledgements First and foremost, I wish to express my sincere gratitude for the efforts of my co-supervisors, Leah Edelstein-Keshet, and Yue-Xian Li. Their guidance and input were essential for the completion of this thesis. Perhaps more im- portantly, they were instrumental in creating opportunities for me to develop my academic career, and they supported this development with much time and encouragement, and meticulous care. I would like to thank the Math Biology faculty, who have given me valu- able feedback throughout my years at UBC, and who stand as a perfect example of a friendly and supportive academic community. I would like to thank the IAM faculty, for their excellent teaching, and for creating a welcoming, exciting academic environment. I also wish to thank the Mathematics Department, and in particular the staff, for all the help along the way. To my fellow IAM graduate students, thanks for the friendly environment you helped create, and for all the entertainment that served as a counterbal- ance to work. Finally, to my family, and especially my wife, thank you for your encour- agement, support, and love at every step of the way. To Nova Scotia we go! xvi Dedication For my beautiful wife, Sionnach. xvii Statement of Co-Authorship Aspects of this thesis resulted from close collaboration with my supervisors, Yue-Xian Li and Leah Edelstein-Keshet. To facilitate details of this collab- oration, I divide the work into two subdivisions: first, the analytical work resulting in Chapters 2 and 3, and the field data collection and analysis resulting in Chapters 4 and 5. In the first subdivision, the research program was identified and designed (in preliminary form) by my advisors. I carried out the analysis (shared with advisors in Chapter 2, independently in Chapter 3), and coded, performed, and interpreted numerical simulations. I shared manuscript preparation for the publication associated with Chapter 2, and was the primary author for the publication associated with Chapter 3. In the second subdivision, I was responsible for the identification and design of the research (in consultation with my advisors). I carried out the research, including experimental design, data gathering, data processing and analysis, and model implementation and analysis. I prepared the manuscripts (to be published) associated with Chapters 4 and 5. xviii Chapter 1 Introduction 1.1 Biological Background Across many species, animals often aggregate into cohesive, ordered groups. Particularly fascinating are the aggregations in which coordination and pat- tern are evident at the group level, in contrast with the many individuals comprising the group. Familiar examples range from schools of fish, herds of ungulates, flocks of birds, and swarms of insects. Clearly, biological aggre- gates form at many different scales. The unifying property of such aggrega- tions is the emergence of order at the group level through interactions among individuals within the group. Through this mechanism, complex group struc- ture can arise from the many relatively simple interactions that simultane- ously take place among group members. Understanding the mechanisms of interaction between individuals that permit cohesive group patterns forms the central focus of this thesis. Aggregative behaviour in animals is generally viewed as a behaviour that optimizes individual fitness in the face of selective pressures acting on the in- dividuals. The particular motivations and pressures vary according to species and situation, but most commonly, membership within a group confers ben- efits through reduction of predatory risk or facilitation in foraging. However, individuals within groups also experience costs associated with the aggre- gate, primarily through increased competition for resources. The size of the group is determined through the competing influences of selective cost and advantage conferred by the aggregate. Considerable work has been done to address these questions of why animals aggregate but this is not the focus of this thesis. Here instead, I address the ‘how’ of aggregations, to focus on what particular mechanisms of inter-individual interaction lead to what group-level properties. In this way, we seek to understand the emergent or- der that is observed in large groups, in terms of individual properties and behaviours. 1 1.2. Mathematical Modeling 1.2 Mathematical Modeling Mathematical modeling is the primary tool for exploring the connection be- tween individual properties and group properties. Using models, researchers are able to test sets of interaction mechanisms and visualize the resulting group behaviour, generally not predictable directly from the rules alone. Al- though particular implementations may differ dramatically, interactions are generally modeled as a combination of any or all of the following three in- fluences: repulsion away from individuals very close, attraction to individ- uals far away, and alignment with nearby neighbours. Lagrangian models [3, 13, 24, 25, 30, 36, 47, 56, 60], based on tracking the positions and veloci- ties of individuals contrast with Eulerian models [1, 18, 19, 20, 22, 29, 31, 34, 39, 53, 54] which describes the local flux of individuals via population den- sity [54]. Lagrangian models are thus comprised of high-dimensional coupled ordinary differential equations (ODEs), whereas Eulerian models are formu- lated in terms of advection-diffusion partial differential equations (PDEs). Eulerian models are most suited to very large populations (insects, bacte- ria, etc.), while Lagrangian models are best suited for smaller groups with distinguishable individuals. Eulerian models are derived from stochastic versions of Lagrangian mod- els (under suitable approximations). Eulerian models are more easily an- alyzed in many cases using tools of partial differential equations, whereas Lagrangian models provide more information about individuals, at a cost of analytical difficulty due to their high dimension [21]. 1.3 Survey of the Literature : Models 1.3.1 Eulerian Models Interactions in Eulerian models can be either local (e.g., [17, 31]) or non- local (e.g., [18, 19, 34, 53, 54]). Local interactions represent responses to only immediate neighbours; in the continuum, interactions at some spatial location depend only on the density at that single point in space [17]. In constrast, non-local interactions represent responses to neighbours within some region (of non-zero radius), which leads to integro-differential equations. 2 1.3. Survey of the Literature : Models A typical 1-D non-local Eulerian model is given by ∂ρ ∂t = ∂ ∂x ( D ∂ρ ∂x − v(ρ)ρ ) , where x is the spatial coordinate, ρ is density of organisms, D the diffusion coefficient, and v(ρ) the velocity [54]. To incorporate non-local effects, v is typically modeled as a spatial convolution. The details of interaction choices for a given model are given by the formulation of this convolution, in par- ticular the type of weighting function used. Further modifications include a density-dependent diffusion coefficient D(ρ) [1] instead of constant D [29]. Eulerian models have been used to analyze a variety of solutions, including traveling band solutions [17], vortex-type solutions [53] stationary clump so- lutions [54], and numerically generate an array of more complex solutions such as zig-zag movement and ‘breathers’ [18]. However, the abstraction of individuals to a local density obscures important information regarding individual properties within groups, and analysis is almost entirely limited to one-dimensional cases. In this thesis, we focus entirely on the Lagrangian modeling perspective, which maintains the integrity of individuals within the group. 1.3.2 Lagrangian Models Lagrangian models have the benefit of being very intuitive: individual in- teractions of attraction, repulsion, and alignment are encoded explicitly as forces acting between individuals, together with some autonomous forces (i.e, drag, self-propulsion). Many Lagrangian models were developed to simulate behaviour of particular animal groups. Sakai [47] and Suzuki and Sakai [49] constructed a simulation model based on equations of motion, with forces of thrust, drag, attraction/repulsion, and alignment. By varying the relative strength of noise to interaction forces, these studies showed the existence of different solution types : amoeba-like groups, doughnut (milling) shaped groups, and rectilinear motion (polarized groups). Although examination of the model was limited, this work established a model framework for collective motion in animal groups. For example, self-propulsion and drag forces used in the model of Chapter 2 and 3 are based on thrust and drag from [47] and [49]. Aoki [2, 3, 4] simulated the movement of a fish school by modeling speed and direction of individuals as stochastic variables, with interactions with one 3 1.3. Survey of the Literature : Models neighbour, chosen probabilistically (as a function of distance). Aoki’s model was zonal, in that it partitioned local space around an individual into radial regions of repulsion, alignment, attraction, and searching response. In [2], individual properties were shown to mediate positioning within the group. Numerous zonal models based on the same principle of interactions restricted to distinct spatial regions followed, including the model I introduce in Chap- ter 5 for flocks of ducks swimming collectively. Huth and Wissel [25, 26] compared the decision-based algorithm of [3] versus an averaging of inter- actions with some number of neighbours, showing the averaging mechanism to lead to more robust schools, and realistic polarity, whereas decision-based schools were less organized. Huth and Wissel’s model was zonal (following [3]), and included a ‘blind angle’ behind an individual where no interaction takes place. From different motivations (computer animation), Reynolds [46] presented a model of flocking for ‘boids’, (artificial birds) based on attraction, repul- sion and alignment with neighbours. Though more complicated forces were included (e.g., banking) for realism, the model showed how realistic flocking behaviour can result from purely local interactions, without external control or ‘leaders’. Niwa [36, 37, 38] developed a stochastic differential equation model for fish schooling, composed of a locomotory force (derived from swimming en- ergetics), a ‘grouping force’ of attraction/repulsion, and alignment. Group properties were investigated as the level of noise relative to interactions was varied. Niwa showed that the group exhibits transitions between behaviours as noise levels change, from amoeba-like, disorganized groups, to polarized groups moving rectilinearly (similar to behaviors described in [3]). Niwa also showed that by tuning randomness to a particular level, noise can facilitate schooling by reducing the time of onset of schooling. Warburton and Lazarus [60] investigated the relationship between the shape of interaction force (attraction/repulsion) on group properties, show- ing that although stable groups form across most functions with short-range repulsion and long-range attraction, the shape of these response curves affects nearest-neighbour distance, group shape, and level of cohesion. Although limited to generic interaction curves and conclusions via simulation, [60] es- tablished the importance of interaction function shape on group properties, a link I study analytically in Chapters 2 and 3. Mogilner et. al [35] considered specific types of interaction functions (i.e., exponential, inverse power) and analytically derived conditions on function 4 1.3. Survey of the Literature : Models parameters to guarantee stable, cohesive well-spaced groups. In contrast to this thesis, in [35] individuals are assumed to be coupled all-to-all (i.e., global sensing of neighbors) and particular functional forms were used to obtain analytical conclusions. Gueron et. al. [24] developed a zonal model for animal groups based on distance-dependent interactions. Between zones of repulsion and attraction, the model included a neutral zone where neighbours imposed no forcing. By investigating the level of cohesion permitted by varying widths of the neutral zone, the authors argued that neutral zones of certain width permit cohesion, but reduce the amount of acceleration and deceleration based on neighbour position, which would confer an energetic benefit. Vabo and Nottestad [55] presented a cellular automata model to investigate group dynamics in the presence of a predator, and was able to generate behaviours similar to those seen in nature, including ‘vacuole’, ‘fusion’, and ‘fountain’ formations. Viscido et. al. [58] investigated the dependence of group dynamics on number of interacting neighbours. They found that when the number of in- fluencing neighbours is similar to total population size, static, milling-type solutions occur, whereas when population size is much larger than number of influencing neighbours, mobile, polarized groups form. In [59], a variety of metrics were proposed to investigate the parameter-dependence of group properties, which suggested that a neutral zone, relative magnitude of align- ment, and number of influential neighbours are crucial to determining the group structure. In Chapter 5, I investigate relative magnitudes of alignment, attraction, and repulsion that best reproduce features of observed groups. Couzin et. al. [13], using a fixed-speed zonal model (repulsion, align- ment, attraction), investigated group behaviour as zone widths are changed. By changing width of the alignment zone, group structure changed from disorganized, to milling, to polarized. Interestingly, the authors showed that these transitions exhibited hysteresis: the previous history of group structure influences collective behaviour as individual interactions change. Stated an- other way, the current geometry of the group, in conjunction with individual properties, determines group response to parameter changes. In Chapters 2 and 3 of this thesis, I show co-existence of two different solutions over certain parameter regimes; whether a group is found in one (polarized) or the other (milling) depends on the geometry of the group. Couzin et. al [12] used a similar model as in [13], but included some proportion of ‘informed individ- uals’, who had some desired direction of travel together with an alignment force. They showed that as group size increases, the proportion of individuals 5 1.3. Survey of the Literature : Models needed to guide the entire group decreases. Lagrangian Models in Physics. Though the Lagrangian models above are typically biologically inspired, there have also been significant contribu- tions from the physics literature. Vicsek et. al. [15, 56] used a fixed-speed model of particles interacting via heading matching of neighbours (without attraction/repulsion interactions) showed transitions in group behaviour. In particular, solutions where no net transport occurs (quasi-stationary) switch to finite-transport solutions through symmetry-breaking of the rotational symmetry. These changes in group structure were invoked by tuning particle density and level of noise in the system. Levine and Rappel [30] showed vortex-type solutions (mills) that self-organized without any external in- fluence or chemotaxis, and then derived an advection-diffusion continuum analogue to the model via coarse-grain averaging. D’Orsogna et. al. [16] categorized different solution types via approximate analogy with canonical dissipative systems. This permitted a statistical mechanical analysis of the model, and parameters of the Morse-type interaction function were used to derive conditions for H-stable (well-spaced) versus catastrophic increasing density with increasing population) solutions. Chuang et. al. [11] derived a 2D continuum analogue to [16], and performed a linear stability analysis on stationary solutions. Lagrangian Models in Engineering. From the engineering perspective, work has been done in developing decentralized control rules for systems of agents subject to various movement constraints, with application to, e.g., cooperative robotics [9], and unmanned aerial vehicles [28]. Tanner et. al [50, 51] developed control laws on interacting agents that give rise to tight formations and collision avoidance. Sepulchre et. al. [48] proposed a design methodology to stabilize parallel and circular motion of particles moving at unit speed in the case of each agent interacting with all other agents. Marshall et. al. [33] studied cyclic pursuit from a control perspective, analyzing the particular solution of fixed-speed agents following one another in a circular manner. Existence and (linear) stability conditions were derived based on control inputs. The solution analyzed in [33] is similar to the milling solution I study in Chapter 3, although individuals are subject to different constraints such as fixed velocity, which affects the anaylsis significantly. 6 1.4. The Need for Empirical Data 1.3.3 Unifying Lagrangian Approaches Although all the models described above share the common property that individuals interact with a combination of repulsion, attraction, and align- ment forcing, clearly there is a wide array of specific implementations. The variety of models stems not only from species-specific choices, but also in ad-hoc choices made by modelers where no data exists to guide such choices. In an attempt to provide a unifying framework to such modeling pursuits, Parrish et. al. [40] have organized a number of previous Lagrangian mod- els according to a number of model features. Below, model properties are categorized in a similar way to [40], with some modification. • Population size: the number of individuals simulated by the model. • Speed: models vary according to whether speed is constant (e.g. [13, 25, 26, 55]), varying according to equations of motion (e.g., [16, 30, 36, 37, 38]), or drawn from a probability distribution (e.g., [3]). • Interaction (spatial dependence): Interactions can be zonal, (e.g., [3, 12, 13, 24, 25, 26]) or based solely on linear distance (e.g., [16, 30]. Zonal models provide a framework for imparting a hierarchy of decisions in the interaction, (e.g., aligning and attracting only if no individuals are within the repulsion zone). • Neighbor scaling: given a set of influencing neighbours, relative weight- ing of influence imparted by each neighbour can be constant or distance- weighted. • Rule size: models vary according to how many neighbours are ‘sensed’ by a given individual, ranging from 1 (nearest-neighbour interaction) to every individual (all-to-all coupling). 1.4 The Need for Empirical Data Clearly, given the variety of choices above, there are many possible models. In order to choose between implementations, (in)validate models, and choose realistic parameters, empirical data of group motion is required. However, gathering data on moving groups of animals is a typically difficult under- taking, whereas a mathematical model with arbitrary choices of interaction 7 1.5. Survey of the Literature: Empirical Studies can be implemented with ease. Obtaining empirical data in the field is com- plicated by the typically three-dimensional nature of animal aggregations (resulting in occlusion of interior individuals), the difficulties of calibrating measurement equipment at varying locations, high speeds of movement, and the transient behaviour (both in space and time) of animals in the field in general. Overcoming these difficulties by moving to a controlled labora- tory setting places restrictions on the spatial extent and size of groups, and (depending on the species and experimental setup) can create difficulty in obtaining natural behaviours within an artificial environment. Furthermore, for many species (e.g., flocking birds), a laboratory setting is simply not feasible. 1.5 Survey of the Literature: Empirical Studies 1.5.1 Small Groups, No Tracking Despite the difficulties mentioned above, significant efforts have been made to record positions and movements of individuals within groups. Early work focused on recovering individual positions in relatively small groups, but did not link individuals across frames to construct trajectories (and thus had no dynamic component). Cullen et. al. [14] used stereo photography to reconstruct 3D positions of 10 fish over 11 photographs, showing that pilchards tend to prefer neighbours at diagonal positions, both in bearing and elevation. Major and Dill [32], also using stereo photography, reconstructed 3D positions of flocks of dunlin and starlings, in flocks of 25-76 individuals. They found that dunlin have a tighter, more cohesive structure than starlings, and that dunlin have a propensity for nearest neighbours behind and below a focal individual, while few nearest neighbours occupy relative positions in front and below. Budgey [8], motivated by birdstroke tolerance in aircraft, reconstructed positions of starlings, rock doves, lapwings, mixed gulls, and Canada geese within flocks of 21-61 birds using the stereo method in the field. Nearest-neighbour distances were calculated and a linear relationship between nearest-neighbour distance and wingspan was postulated. Aoki and Inagaki [5] gathered three-dimensional positions of anchovy groups at night via a stereo pair of cameras dragged by a drifting vessel, and showed that anchovies have a tendency at night to swim at similar depths as 8 1.5. Survey of the Literature: Empirical Studies neighbours, and also to have neighbours in front of or behind. Furthermore, the nearest-neighbour distance at night was increased as compared to day- time measurements. Ikawa et. al. [27] modified the orthogonal method by adding a third camera to reconstruct positions of 5-10 mosquitoes in a swarm, and calculated nearest-neighbour distance and some speeds for individuals in the swarm. 1.5.2 Small Groups with Tracking By linking positions of an individual through successive frames in time, a trajectory can be reconstructed, giving dynamic information on interactions. Partridge [41] studied fish (20-30 saithe) schooling in an annular tank, cap- tured by a camera mounted on a rotating gantry, with positions calculated by the shadow method, where overhead positions and shadows cast by a light source are matched to calculate 3-D position. Individual fish were uniquely ‘freeze-tagged’ and so individual trajectories could be tracked in successive frames. The data was used to compare two volume-estimation methods (the total volume method versus the individual ‘free-space envelope’ method). Using the same experimental setup as [41], Partridge et. al. [44] compared an obligate schooler (herring), a strongly facultative schooler (saithe), and a weakly facultative schooler (cod), showing differences in spacing, angular preference of neighbours, and school shape. Also, saithe were shown to decrease NND as either school speed or school size increased. Partridge [43] further analyzed saithe schools via correlations between heading and speed of individuals and their nearest neighbours as a function of lag, to investigate response times of a focal fish to the behavior of its neighbour. This analysis showed significant speed correlation between a focal fish and its first nearest neighbour, maximal at small response times (less than 0.3 sec), but low mean heading correlation due to high variability among individuals. Partridge [42] studied small groups (2-6) of minnows in a tank using the shadow method to reconstruct positions, and calculated ve- locities by tracking individuals through successive frames. Partridge found that two-fish schools behave differently compared to 3-6 fish schools, in terms of neighbour correlations, spacing, and relative positions of neighbours, sug- gesting that results from two individuals cannot be generalized to schools. A natural extension of this observation is to question whether or not small groups (10-20 individuals) accurately represent behavior found in large ag- gregations, the answer to which lies in data sets of large groups. 9 1.5. Survey of the Literature: Empirical Studies Pomeroy and Heppner [45] used orthogonal photography (where images are captured simultaneously by orthogonal cameras) to capture the move- ments of flocks of 12-16 rock doves during a turn, finding that rock doves reposition themselves relative to the flock throughout a turn such that the arc ascribed by each individual varies less than if the group turned as a rigid body. More recently, Grunbaum et. al. [23] used pairs of video cameras record- ing movements of giant danios in a tank, and then tracked individuals using purpose-built tracking software, to generate long-time three-dimensional tra- jectories of groups of 4-8 fish. By generating probability density functions (PDF) of nearest-neighbour distance, and calculating empirical orthogonal functions of the PDF, advective fluxes were calculated for an advection- diffusion equation for nearest-neighbour position. Using the same experi- mental setup, Viscido et. al. [57] compared empirical observations of schools with simulation data, concluding that in order to reproduce realistic features of schools, the magnitude of alignment forcing in models should be almost two orders of magnitude smaller than attraction/repulsion forcing. Tien et. al. [52] tracked the movements of 40 fish (creek chubs and blac- knose dace) enclosed in a shallow natural pool in a creek via an overhead camera. The shallow water essentially restricted movement to two dimen- sions. By analyzing movements of fish relative to nearest neighbours, it was shown (to some degree of correlation) that these fish respond to neighbours in a series of radial zones: a nearby region of mutual repulsion, then a region of neutrality, followed by a region of attraction, where the attraction response is made only by the focal fish (i.e, not mutual attraction). In the presence of a simulated predator, the radial zones decreased in size. 1.5.3 Large Groups, No Tracking A significant breakthrough in empirical data collection was made by Ballerini et. al. [6] who captured three-dimensional positions of starling flocks of up to 3000 individuals at 10 frames per second, using a trifocal stereo photography method. Their (static) analysis of these flocks found a strong anisotropy in the relative position of nearest neighbours. By quantifying this anisotropy, the number of neighbours with which an individual interacts (on average) was determined. Importantly, it was found that the number of interacting neighbours was not dependent on density, suggesting that individual interac- tion is ruled by topological (i.e., a specific number of neighbours) interaction, 10 1.6. Contributions of this Thesis as opposed to metric (i.e., all neighbours within a fixed sensing region) in- teraction. In a related work based on the same data set, Cavagna et. al. [10] introduce statistical measures of density as a function of sensing radius of a focal individual. It was shown here that starling flocks have higher densities at the edge, and individuals have weak ‘shells’ of neighbours - in between random and crystalline. As I will show in Chapter 5, similar ‘shells’ of neigh- bors are observed in scoter flocks I study. In another related work, Ballerini et. al. [7] reported group statistics (morphology, orientation, turning dynam- ics, density, nearest-neighbour distance) of the starling flocks. Although this data set represented a significant step forward, incomplete reconstructions of positions (only 80% of individuals were reconstructed in any given frame) prevented any tracking of individuals in time. 1.5.4 What Empirical Data is Needed Although there have been numerous empirical studies involving a few in- dividuals tracked in time (i.e., having reconstructed trajectories), or many individuals without tracking (i.e., data limited to static images), there is a void of empirical data for tracking many individuals simultaneously. The main obstacle in obtaining such data is the need for reconstructing all indi- viduals in each time frame in order to track individuals effectively. The 80% position reconstruction rate of [6] precludes trajectory reconstruction across usable time scales, as the probability of an n-frame trajectory is 0.8n, which implies a trajectory reconstruction success rate of 10% over 10 frames (cor- responding to one second, in [6]). Yet, tracking of only few individuals raises the issue of whether or not the number of individuals is sufficient to trigger the schooling behaviour in groups. Furthermore, in small groups, edge ef- fects (where individuals at the exterior of a group have different statistical properties as compared to interior individuals) dominate analysis, and any spatial description of interactions is limited by the small spatial extent of small groups (especially so in three dimensions). 1.6 Contributions of this Thesis In this thesis, collective motion in animal groups is investigated in two ways: by developing and applying analytical tools to explore a Lagrangian model, where previously, exploration was limited to numerical simulation; 11 1.6. Contributions of this Thesis and through gathering and analyzing empirical data against which models can be tested. 1.6.1 Model Analysis In the first part of this thesis, we study a general Lagrangian model of n (i = 1, . . . , n) self-propelled particles, of the form ~̇xi = ~vi, (1.1) ~̇vi = ~fa + ~fi, (1.2) where ~̇x ≡ d~x/dt denotes a derivative with respect to time, ~fa is a force that is generated autonomously by each individual without the influence of other individuals and ~fi is a force that is generated due to the interaction with oth- ers. Because ~fi represents interindividual interactions, Eqs. (1.1)-(1.2) are coupled, often with relatively complex connection structures. Furthermore, terms in Eq. (2.2) are typically nonlinear. Thus, model equations cannot be solved explicitly, and steady-state existence and stability analyses are typi- cally not possible. For these reasons, Lagrangian models have been studied primarily by simulation; i.e., quantifying changes in model output resulting from parameter variation, and categorizing different types of behavior exhib- ited by the model. Although such efforts are useful, they lack the clarity afforded by analytical observations, where the relationship between model parameters and group-level properties is clear-cut. In this thesis, we consider a special class of solutions, termed perfect schools : geometrically simple organizations of individuals with fixed spac- ing. For these solutions, the interindividual coupling and equilibrium states can be written explicitly, leading to existence conditions. Then, the simpli- fied coupling in perfect school solutions is exploited to obtain a matrix system of Eqs. (2.1)-(2.2) where the coefficient matrix is structured. This structure permits a linear stability analysis on the system for particular perfect school solutions. These analytical results give clear expressions for how group struc- ture and stability depend on model inputs (i.e., individual properties). This complete characterization of model outputs from inputs is missing in simu- lation studies, which can only draw inferences from limited parameter space exploration. In particular, I derive existence conditions dependent on self- propulsion, drag, and interaction function strength. I then derive conditions on stability of perfect school solutions in terms of the interaction function, and its derivative. 12 1.6. Contributions of this Thesis The complexity of Lagrangian models lies in the many interactions among individual units. Though interactions may be relatively simple, a complex web of interactions creates feedbacks in the system that make emergent group properties difficult to predict. Herein lies the strength of an analytical treat- ment: analysis, as in this thesis, eliminates the unpredictability of group structure that emerges from model inputs because the analytical results link emergent group properties to individual properties via straightforward math- ematical relationships. Although the types of solutions considered in this work are overly sim- ple when compared to real-world animal groups, the analysis nevertheless is useful for understanding qualitative aspects of such animal groups. For example, general properties of attraction and repulsion between individuals can be deduced based on observed spacing, velocity, and group geometry. In Chapter 2, one-dimensional solutions, and ‘soldier formations’ in two dimen- sions (wherein individuals are organized in a linear formation) are studied. In Chapter 3, milling formations (equidistant individuals following one another around a closed circular path) are studied. I show how different solutions can co-exist for sets of parameters, depending on the arrangement of individuals (i.e., the model state). 1.6.2 Empirical Data In the second part of this thesis, to address the void of dynamic data for relatively large groups of animals, empirical data is gathered on groups of a few hundred surf scoters (diving ducks) swimming collectively on the surface of the water in English Bay, near Vancouver, BC, Canada. This data is significant, as it represents a 10-fold increase in group size wherein individuals are still reliably tracked. Thus, for the first time (to our knowledge), complete trajectories (and therefore individual velocities) are known for all members of a relatively large group, allowing for interindividual responses to be measured in time. Because models of collective motion are naturally dynamic (e.g., time-dependent ODEs), trajectories provide a necessary and valuable data set against which models can be formulated and validated. Data processing tools are created that allow for complete reconstruction of positions of individual ducks in each frame, and trajectories of individuals are reconstructed using existing particle-tracking software. Methods are de- veloped to correct observations for camera angle and currents in the water, and to correct tracking errors caused by deficiencies in the existing software. 13 1.6. Contributions of this Thesis A complete description of empirical methods as well as basic group statistics (mean velocity, mean spacing, etc.) are reported in Chapter 4. In Chapter 5, two-dimensional spatial distributions are used to characterize individual interactions, and a model that captures some aspects of the empirical data is presented and studied via numerical simulation. As I show in this chapter, to account for observations, we infer that ducks have a hierarchy of rules, have an angular bias in their attraction to neighbors, and exhibit different avoidance manoeuvres depending on relative angle between individuals. A concluding chapter follows in Chapter 6. 14 Bibliography [1] W. Alt, Degenerate diffusion equations with drift functionals modeling aggregation, Nonlinear Anal. 9 (1985), 811–836. [2] I. Aoki, An analysis of the schooling behavior of fish: internal organi- zation and communication process, Bull. Ocean. Res. Inst. 12 (1980), 1–65. [3] , A simulation study on the schooling mechanism in fish, Bull. Jpn. Soc. Sci. Fish. 48 (1982), 1081–1088. [4] , Internal dynamics of fish schools in relation to inter-fish dis- tance, Bull. Jpn. Soc. Sci. Fish. 50 (1984), 751–758. [5] I. Aoki and T. Inagaki, Photographic observations on the behaviour of japanese anchovy engraulis japonica at night in the sea, Mar. Ecol. Prog. Ser. 43 (1988), 213–221. [6] M. Ballerini, N Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Gi- ardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic, Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study, Proc. Nat. Acad. Sci. 105 (2007). [7] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giar- dina, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic, An empirical study of large, naturally occurring starling flocks: a bench- mark in collective animal behaviour, arXiv:0802.1667 (2008). [8] R. Budgey, Three dimensional bird flock structture and its implications for birdstroke tolerance in aircraft, International Bird Strike Committee IBSC 24/WP 29 (1998). 15 Chapter 1. Bibliography [9] Y. Cao, A. Fukunaga, A. Kahng, and F. Meng, Cooperative mobile robotics: antecedents and directions, Proc. IEEE/RSJ Int. Conf. Intelli- gent Robots and Systems (1995), 226–234. [10] A. Cavagna, A. Cimarelli, I. Giardina, A. Orlandi, G. Parisi, A. Procac- cini, R. Santagati, and F. Stefanini, New statistical tools for analyzing the structure of animal groups, Mathematical Biosciences 214 (2008), 32–37. [11] Y. L. Chuang, M. R. D’Orsogna, D. Marthaler, A. L. Bertozzi, and L. S. Chayes, State transitions and the continuum limit for a 2d interacting, self-propelled particle system, Physica D 232 (2007), 33–47. [12] I. D. Couzin, J. Krause, N. R. Franks, and S. A. Levin, Effective lead- ership and decision-making in animal groups on the move, Nature 433 (2005), 513–516. [13] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and N. R. Franks, Collective memory and spatial sorting in animal groups, J. Theor. Biol. 218 (2002), 1–11. [14] J. Cullen, E. Shaw, and H. Baldwin, Methods for measuring the three- dimensional structure of fish schools, Can. J. Zool. 71 (1993), 1494–1499. [15] A. Czirók, M. Vicsek, and T. Vicsek, Collective motion of organisms in three dimensions, Physica A 264 (1999), 299–304. [16] M. R. D’Orsogna, Y. L. Chuang, A. L. Bertozzi, and L. S. Chayes, Self- propelled particles with soft-core interactions: Patterns, stability, and collapse, Physical Review Letters 96 (2006), 104302. [17] L. Edelstein-Keshet, J. Watmough, and D. Grunbaum, Do travelling band solutions describe cohesive swarms? An investigation for migratory locusts, J. Math Biol 36(6) (1998), 515–549. [18] R. Eftimie, G. de Vries, and M. A. Lewis, Complex spatial group pat- terns result from different animal communication mechanisms, Proc. Nat. Acad. Sci. 104 (17) (2007), 6974–6979. [19] R. Eftimie, G. de Vries, M. A. Lewis, and F. Lutscher, Modeling group formation and activity patterns in self-organizing collectives of individ- uals, Bull. Math. Bio. 69 (2007), 1537–1565. 16 Chapter 1. Bibliography [20] G. Flierl, D. Grunbaum, S. Levin, and D. Olson, From individual to aggregations: the interplay between behaviour and physics, J theor Biol 196 (1999), 397–454. [21] D. Grunbaum, Translating stochastic density-dependent individual be- havior with sensory constraints to an eulerian model of animal swarm- ing, Journal of Mathematical Biology 33 (1994). [22] D. Grunbaum and A. Okubo, Modelling social animal aggregation., Frontiers in Mathematical Biology (S. Levin, ed.), Springer, N.Y., 1994, pp. 296–325. [23] D. Grunbaum, S. Viscido, and J. Parrish, Extracting interactive control algorithms from group dynamics of schooling fish, Cooperative Control LNCIS 309 (2004), 103–117. [24] S. Gueron, S. A. Levin, and D. I. Rubenstein, The dynamics of herds: From individual to aggregations, J theor Biol 182 (1996), 85–98. [25] A. Huth and C. Wissel, The simulation of movement of fish schools, J. Theor. Biol. 156 (1992), 365–385. [26] , The simulation of fish schools in comparison with experimental data, Ecological Modelling 75/76 (1994), 135–145. [27] T. Ikawa, H. Okabe, T. Mori, K. Urabe, and T. Ikejoshi, Order and flexibility in the motion of fish schools, Journal of Insect Behavior 7(2) (1994), 237–247. [28] E. W. Justh and P. S. Krishnaprasad, Equilibria and steering lawes for planar formations, Systems and Control Letters 52(1) (2004), 25–38. [29] K. Kawasaki, Diffusion and the formation of spatial distributions, Math. Sci. 16(183) (1978), 47–52. [30] H. Levine, W.J. Rappel, and I. Cohen, Self-organization in systems of self-propelled particles, Physical Review E 63 (2001), 017101. [31] F. Lutscher, Modeling alignment and movement of animals and cells, J. Math. Bio. 45 (2002), 234–260. 17 Chapter 1. Bibliography [32] P. F. Major and L. M. Dill, The three-dimensional structure of airborne bird flocks, Behavioural Ecology and Sociobiology 4 (1978), 111–122. [33] J. A. Marshall, M. E. Broucke, and B. A. Francis, Formations of vehicles in cyclic pursuit, Automatic Control, IEEE Transactions on 49 (11) (2004), 1963–1974. [34] A. Mogilner and L. Edelstein-Keshet, A non-local model for a swarm, J Math Biol 38 (1999), 534–570. [35] A. Mogilner, L. Edelstein-Keshet, L. Bent, and A. Spiros, Mutual in- teractions, potentials, and individual distance in a social aggregation, J Math Biol 47 (2003), 353–389. [36] H.-S. Niwa, Self-organizing dynamic model of fish schooling, J theor Biol 171 (1994), 123–136. [37] , Newtonian dynamical approach to fish schooling, J theor Biol 181 (1996), 47–63. [38] ,Migration of fish schools in heterothermal environments, J theor Biol 193 (1998), 215–231. [39] A. Okubo, Dynamical aspects of animal grouping:swarms, schools, flocks, and herds, Adv. Biophys. 22 (1986), 1–94. [40] J. Parrish, S. Viscido, and D. Grunbaum, Self-organized fish schools: an examination of emergent properties, Biol. Bull. 202 (2002), 296–305. [41] B. L. Partridge, Fish school density and volume, Marine Biology 54 (1979), 383–394. [42] , The effect of school size on the structure and dynamics of min- now schools, Animal behaviour 28 (1980), 68–77. [43] , Internal dynamics and the interrelations of fish in schools, J Comp Physiol 144 (1981), 313–325. [44] B. L. Partridge, T. Pitcher, J. M. Cullen, and J. Wilson, The three- dimensional structure of fish schools, Behav Ecol Sociobiol 6 (1980), 277–288. 18 Chapter 1. Bibliography [45] H. Pomeroy and F. Heppner, Structure of turning in airborne rock dove flocks, Auk 109(2) (1992), 256–267. [46] C. W. Reynolds, Flocks, herds, and schools: a distributed behavioural model, SIGGRAPH: Computer Graphics 21(4) (1987), 25–34. [47] S. Sakai, A model for group structure and its behavior, Biophysics Japan 13 (1973), 82–90. [48] R. Sepulchre, D. Paley, and N. Leonard, Stabilization of planar collec- tive motion: all-to-all communication, Automatic Control, IEEE Trans- actions on 52(5) (2007), 811–824. [49] R. Suzuki and S. Sakai, Movement of a group of animals, Biophysics Japan 13 (1973), 281–282. [50] H. Tanner, A. Jadbabaie, and G. Pappas, Stable flocking of mibile agents, part ii: dynamic topology, 42nd IEEE Conference on Decision and Control 2 (2003), 2016–2021. [51] , Stable flocking of mobile agents, part i: fixed topology, 42nd IEEE Conference on Decision and Control 2 (2003), 2010–2015. [52] J. Tien, S. Levin, and D. Rubenstein, Dynamics of fish shoals: identify- ing key decision rules, Evolutionary Ecology Research 6 (2004), 555–565. [53] C. Topaz and A. Bertozzi, Swarming patterns in a two-dimensional kine- matic model for biological groups, SIAM J. Appl. Math. 65 (1) (2004), 152–174. [54] C. Topaz, A. Bertozzi, and M. Lewis, A nonlocal continuum model for biological aggregation, Bull. Math. Bio. 68 (7) (2006), 1601–1623. [55] R. Vabo and L. Nottestad, An individual based model of fish school reac- tions: predicting antipredator behaviour as observed in nature, Fisheries Oceanography 6(3) (1997), 155–171. [56] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Phys. Rev. Lett. 75(6) (1995), 1226–1229. 19 Chapter 1. Bibliography [57] S. Viscido, J. Parrish, and D. Grunbaum, Individual behaviour and emer- gent properties of fish schools: a comparison of observation and theory, Mar. Ecol. Prog. Ser. 273 (2004), 239–249. [58] S. V. Viscido, J. K. Parrish, and D. Grunbaum, The effect of population size and number of influential neighbors on the emergent properties of fish schools, Ecological Modelling 183 (2004), 347–363. [59] , Factors influencing the structure and maintenance of fish schools, Ecological Modelling 206 (2007), 153–165. [60] K. Warburton and J. Lazarus, Tendency-distance models of social cohe- sion in animal groups, J. theor. Biol. 150 (1991), 473–488. 20 Chapter 2 Minimal Mechanisms for School Formation in Self-Propelled Particles 1 2.1 Introduction A fundamental problem in biology is how behaviour and interactions at the level of the individual impact the emergent behaviour at the level of the group as a whole. In this paper, we examine the connection between individ- ual behaviour and group properties within the context of social organisms that form a school or a flock. We ask how the parameters associated with each individual, together with the effective forces that represent their mutual interactions influence the existence of a school, its geometry, the speed and direction of movement, as well as the stability of the school structure. Our perspective in this paper is based on the Lagrangian viewpoint, that is, we follow individual particles, rather than densities of organisms. We formulate equations of motion based on the Newtonian approach, i.e. de- scribing changes in the velocities and positions of the particles under forces of propulsion and interaction. Such an approach leads to a generalized system of ordinary differential equations that is given below. ~̇xi = ~vi, (2.1) ~̇vi = ~fa + ~fi, (2.2) where ~̇x ≡ d~x/dt denotes a derivative with respect to time. Here, we roughly classified all forces into two general categories: ~fa that is generated au- tonomously by each individual without the influence of other individuals 1A version of this chapter has been published as ‘Y.-X. Li, R. Lukeman, and L. Edelstein-Keshet, 2008, Minimal mechanisms for school formation in self-propelled parti- cles, Physica D, 237(5), p 699-720 21 2.1. Introduction and ~fi that is generated due to the interaction with others. Depending on the specific choice of these forces, the model can have very different proper- ties. For example, in [1, 5, 12, 13, 14], ~fa = (α− β|~vi|2)~vi was studied. This model was originally obtained by minimizing the specific energy cost of a swimming fish [26]. In [8], ~fa = −β~vi which is simply a drag force that helps to stabilize the motion of a particle. In both these models, the movement of an individual is purposeless in the absence of the interaction with other par- ticles. It either stays motionless or keeps moving at a constant speed along its initial direction of motion. In the present study however, there is some preferred direction, and each individual knows where to go and whether it is at a leading position in the group. The interaction force ~fi that is directly responsible for many important characteristics of the emergent aggregate pattern plays a key role in defining a distinct model. Pairwise interaction is the most commonly studied force, although in [8], a localized mean field interaction was used to describe the influence of nearest neighbours on the movement of a particle. Interactions with infinite range and all-to-all coupling have been considered in most stud- ies cited above as well as the study of non-Newtonian particles moving in viscous medium [11]. An important question considered in these studies is how specific features of the aggregate change as the number of particles is in- creased. Whether the nearest neighbour distance (NND) collapses to zero in the limit of an infinitely large group of particles has been shown to influence strongly the possible school patterns that can occur in the system [1, 5, 11]. In the present study, our focus is very different from these studies. We are not concerned with the effects of an infinite size of the group but rather focus on groups with finite number of particles and with short range interactions. In most cases, we only consider nearest neighbour interactions. Our main goals are to find analytic insights, based on simplifications and assumptions that render the problem amenable to analysis. These insights allow us to pinpoint the minimal mechanism that determines each specific feature of a school pattern. This approach is also different from a number of other studies in which numerical insights were the main focus [2, 3, 4, 6, 7, 17, 20, 25]. In order to arrive at analytic results, we consider classes of special solu- tions that we refer to as perfect schools. These are configurations of particles in which the spacing is constant and identical, and where individuals in the group move at a constant speed. We further subdivide perfect schools into those in which individuals share a common heading, and those in which they do not. In this paper, we are primarily concerned with the former. This is 22 2.1. Introduction the Lagrangian analogue to motion as a traveling wave with uniform interior density. It is well known that aggregates occur if particles mutually attract at long distances and repel at short distances. When the interaction is all-to- all, the NND can approach zero as the size of the group increases depending on whether the pairwise interaction satisfies certain conditions [1, 5, 11]. A large number of possible transient and static patterns can occur in particles that interact with each other in such a way, and it is beyond the scope of this paper (and likely impossible) to predict all possible patterns. In this paper, we focus mainly on regular arrays, both in 1D and 2D, that are reached after the transient state is over. Our approach is to exploit the simplicity of these regular schools in an attempt to address the following questions analytically: (i) Under what conditions on the school forces do such perfect schools occur (existence)? (ii) Given their existence, what further conditions guarantee that these patterns are not destroyed by random perturbations (stability)? While real schools and flocks are much more plastic and irregular, focusing on these patterns leads to clear-cut results that can be derived analytically. The relative simplicity of our model and the precise definition of these special solutions allow us to reach several new conclusions that are difficult to obtain by simulations alone. We will show that by looking for perfect school solutions, we arrive at conditions linking self-propulsion forces with interaction forces evaluated at the NND. Stability depends on the slope of the interaction force as a function of NND. The structure of the interaction matrix that determines how many neighbours interact with each individual is crucial for achieving explicit stability analysis. At the same time, in order to extend our results to cases that are too challenging to analyze fully, or to more complex cases, we complement our study with numerical simulations. Our main goal in this paper is to understand which factors lead to which of the properties of a group. As will be shown, the questions posed above can be answered conclusively only when analytical solutions of the schools and explicit expressions for determining their stability can be achieved. Many realistic school structures observed in nature are irregular under the influ- ence of noise from various sources. The range of structures that we explore, however, belong to a class of simple and regular patterns that allow us to achieve analytical results and predictions. Analogues of these structures are often observed in nature, and based on our theory presented in this paper, are the most likely realized patterns given simple rules of nearest neighbour interaction. 23 2.2. The Model of Interacting Self-Propelled Particles The questions posed above can be answered conclusively in the case of the special solutions that form the class of interest. This would not be true if we were to consider all possible (dynamic and steady state) group configurations that satisfy our model equations. For this reason, we do not here consider a variety of more complicated and more realistic group structures that occurs in nature. Nevertheless, the range of structures we do explore are consistent with groups of autonomous robots, or formations of vehicles on a highway or in military environments, topics that have emerged as recent areas of application of such theory. Understanding the possible behaviour of such artificial autonomous self-propelled agents forms an additional motivation for this paper. In Section 2, we introduce the basic assumptions of the model we use in this study and the differential equation system that describes the motion of the particles. In Section 3, we study schools in one dimensional space analytically and numerically. The simplicity of these school solutions allows us to develop a basic theoretical framework. We also introduce the concept of stressed and stress-free schools. In Section 4, a soldier formation solu- tion in two dimensions is studied analytically using a similar approach as in Section 2. Numerical simulations in 2D are then used to verify analytical results on soldier formations. In Section 5, we discuss some numerical results that demonstrate the formation of other types of 2D arrays. Results are summarized and discussed in Section 6. 2.2 The Model of Interacting Self-Propelled Particles 2.2.1 Assumptions and Equations The following basic assumptions apply to the schools that we study in this paper. (A1) All particles are identical and obey the same set of rules. (A2) Each particle is capable of sensing if it is located in the interior or at the edge (front/rear) of the school. (A3) The force experienced by one particle from a given neighbour is com- posed of a distance-dependent term and a velocity-dependent term. 24 2.2. The Model of Interacting Self-Propelled Particles Only pairwise interactions are considered. (A4) The forces exerted by different neighbours are additive. (A5) The magnitude of the forces between a pair of particles depends on the distance between the pair and their relative velocities. (A6) The distance-dependent force between a pair is a vector parallel to the vector connecting the pair while the velocity-dependent force is parallel to the difference between the velocities of the pair. Consider a group of n self-propelled particles, each identified by a sub- script i = 1, 2, . . . , n, with position and velocity denoted by ~xi and ~vi, re- spectively. We assume particles are polarized, that is, each has a front and a back. For simplicity, we assume that the body alignment of the ith particle is identical to its direction of motion, i.e. the direction of ~vi. The movement of the ith particle is governed by the classical Newton’s Law of Motion which yields ~̇xi = ~vi, (2.3) ~̇vi = ~ai − γ~vi + ~fi, (2.4) where the mass of each individual is scaled so that mi = 1 for all particles. ~ai is an autonomous self-propulsion force generated by the i th particle which may depend on environmental influences and on the location of the particle in the school. This is the term that makes this model different from some previously studied models [1, 5, 11, 12]. The term γ~vi is a damping force with a constant drag coefficient γ (> 0) which assures that the velocity is bounded. Eq. (2.4) implies that if a particle stops propelling, its velocity decays to zero at the rate γ. ~fi is the interaction or schooling force that moves the ith particle relative to its neighbour(s). If one particle described by (2.3)-(2.4) is moving in D dimensional space (D = 1, 2, 3), the number of first-order equations for each particle is 2D. For a group of n such particles, there are 2Dn first-order equations. In the following analysis, ~ai is typically a constant vector which may be different for a leader, follower or individual in the tail (i.e., the last particle) of the group. Based on the model assumptions (A3-A6), the schooling force is given by ~fi ≡ ~fxi + ~f vi = ∑ j g±(|~xji|) ~xji|~xji| + ∑ j h±(|~vji|) ~vji|~vji| , (2.5) 25 2.2. The Model of Interacting Self-Propelled Particles where ~xji = ~xj − ~xi; ~vji = ~vj − ~vi. ~fxi and ~f vi are the position and velocity dependent forces, respectively. Index j sums over all the particles that influence the movement of particle i. g+ and h+ (g− and h−) are chosen when the influence comes from the front (behind), i.e. if ~xji · ~vi > 0 (< 0). g+ 6= g− and h+ 6= h− imply that the forces from particles in front and those from particles behind are different or asymmetric. g±(x) are defined for x ≥ 0, i.e., are functions of the absolute distance between two particles. To produce a nonzero spacing distance, g±(x) should typically be positive for large x and negative for small x, indicative of short-range repulsion and long-range attraction. This feature is universal in almost all models of aggregate formation. (See review of attraction-repulsion models in [11].) In most results that we studied in this paper, g±(x) possesses such a feature, unless specified otherwise. Furthermore, the forces that each particle can generate should be bounded. Based on these definitions, g+(·) > 0 (< 0) means that particle j which is in front of particle i exerts an attractive (repulsive) force on the latter; while h+(·) > 0 (< 0) implies that the velocity-dependent force makes the velocity of particle i converge to (diverge from) that of particle j (see Fig. 2.1). Note that the position-dependent and velocity-dependent components of the schooling force can point to different directions and possess different magnitudes. All constants and functions that appear in (2.3)-(2.5) are model inputs that are assumed to be given. Figure 2.1: A schematic diagram of model vectors in 2 dimensions. Different choices of the functions g± and h± can generate large numbers of school formations that are qualitatively different. Further classification of the subtypes of these models is possible based on specific selection of these 26 2.2. The Model of Interacting Self-Propelled Particles functions. We do not intend to take on such a task in this paper. However, examples that we present in the following sections clearly demonstrate the crucial importance of the form of these functions in determining the school patterns that emerge. Using the classical Newton’s Law of Motion to describe a system of self- propelled particles is an approach adopted in numerous previous studies [8, 15, 16] although, as we shall discuss later, differences in the dynamics that govern the autonomous self-propulsion and in the range of pairwise in- teractions can lead to significant diffferences in the school patterns that can emerge. 2.2.2 Classification of Schools Having given the definition of perfect schools in Section 1, we note that such schools can be further divided into two types: those in which each individual has an identically constant direction, denoted as Type I schools, and those in which the individuals change direction in time, denoted as Type II schools. Type I schools are the main focus of this paper. Type II schools are typically more difficult to analyze. A good example is a mill formation, in which particles rotate around an invariant circle (see, e.g., [10]). Here, headings are changing continuously as particles rotate, and each individual has a different heading at any given time. Other Type II schools include groups collectively performing a turning manoeuvre. Note that in both Types I and II perfect schools, the group moves collectively as a rigid body under transformation (e.g., translation, rotation, etc). For a perfect school, we can define the center of mass and the school velocity by ~X(t) = 1 n n∑ i=1 ~xi(t), ~V (t) = 1 n n∑ i=1 ~vi. (2.6) The combined autonomous self-propelling force and schooling force are: ~A = 1 n n∑ i=1 ~ai, ~F = 1 n n∑ i=1 ~fi. (2.7) Using these definitions, we can write down the equations of motion of the school as a whole. ~̇X = ~V , (2.8) ~̇V = ~A + ~F − γ~V . (2.9) 27 2.2. The Model of Interacting Self-Propelled Particles In the special case ~vi = ~V , where ~V is a constant vector and identical for all i (i.e., a Type I school), the equations of motion for each individual must be identical to those of the whole school. We distinguish between two types of individuals: those denoted as leaders, and their followers. In general, we define a leader particle as an individual that does not see any other particles in its frontal visual field, where the sign of ~vi · ~xji distinguishes the frontal and rear visual fields. All others are denoted as followers. We associate distinct self-propelling forces ~al and ~a with the leaders and the followers respectively (where superscripts are labels, not exponents). In one-dimensional schools, we consider the special case of interactions both to the front and rear of an individual. In this case, we additionally distinguish between a tail particle (i.e., an individual that does not see any particles in its rear visual field) and other interior followers, and associate a distinct self-propelling force ~at with the tail particle. In our model of school formation, we consider as inputs the leader and follower autonomous self-propulsions ~al and ~a, and the interaction functions g± and h±. Model outputs are the individual distance d between particles, the school velocity ~v and (where applicable) the autonomous self-propulsion of the tail ~at. The perfect schools that we study in this paper are idealizations of some observed social aggregates in nature. Fish schools in the form of surface sheets or linear soldier formations are often encountered. Single-file trail following occurs in ants, caterpillars, as well as birds; linear and ∨-shaped flocks are observed in formations of many bird species. The analytical methods we adopt here are very similar to those used in the study of synchronized clustered states in neuronal networks [9]. We substitute a specific school solution of interest into the system of nonlin- ear differential equations and determine the conditions for the existence of such a solution. Then, we introduce a small perturbation in the equations near this school solution and linearize the system. We find the eigenvalues of the linearized system, and we determine conditions that make their real parts negative. We then verify these existence and stability conditions using numerical simulations of the system, and explore patterns that arise from arbitrary initial conditions. 28 2.3. Schools in One Dimensional Space 2.3 Schools in One Dimensional Space Aggregates of animals moving on a single trail are examples of schools in one dimensional space (1D). Vehicles moving on a single-lane highway can form similar aggregates. Consider a group of n particles moving in 1D. The equations of motion are now reduced to the scalar equations: ẋi = vi, (2.10) v̇i = ai − γvi + fi, (2.11) where i = 1, 2, . . . , n. For convenience, we associate the label 1 with the leading particle, and sequentially label the remaining n−1 particles (see Fig. 2.2). Then a1 ≡ al and ai ≡ a for i = 2, . . . , n− 1, and an ≡ at. Figure 2.2: A schematic diagram of particles in 1 dimension, showing our subscripting convention. Grey arrows indicate distance-dependent interac- tion forces. Once a perfect school in 1D is formed, only two quantities are required to completely characterize the school: its velocity, v, and nearest-neighbour distance, d. These will be determined by the choice of specific interaction forces and parameter values. According to Eq. (2.5), the schooling force is given by fi = ∑ j g±(|xji|) xji|xji| + ∑ j h±(|vji|) vji|vji| , (2.12) where xji = xj − xi is the distance between neighbours and vji = vj − vi is the relative velocity of particles i and j, and j sums over all particles that influence the movement of particle i. We change into a coordinate system that is moving at the school speed v. Let pi and qi be the position and velocity of the i th particle in this moving coordinate. Using the relation pi = xi − vt and qi = vi − v, we obtain the dynamic equations in the moving frame: ṗi = qi, (2.13) q̇i = ai − γ(qi + v) + fi. (2.14) 29 2.3. Schools in One Dimensional Space In the moving frame, a perfect school solution corresponds to a stationary distribution of the particles in space. Thus, every perfect school must be a steady state solution of (2.13)–(2.14); i.e., ṗi = 0, q̇i = 0, for which qi = 0, (2.15) ai − γv + fi = 0, (2.16) for all i. For such a steady state, fi = ∑ j g±(|j−i|d) sgn(j−i)+ ∑ j h±(0)0 = ∑ j g±(|j−i|d) sgn(j−i), (2.17) where j sums over all particles that influence particle i and where sgn(j − i) = j − i|j − i| . Eq. (2.17) implies that for a steady-state solution in which interaction con- nections are fixed, all the schooling forces are functions of d only. 2.3.1 Schools Formed by Following One Immediate Neighbour The simplest case of a perfect school in one dimension occurs when each particle follows only one nearest neighbour at constant distance d in front of it. Then the leader feels no schooling force, i.e., f1 = 0, so a l − γv = 0 while all followers including the last particle experience the same schooling force fi = g+(d). Because the tail particle interacts in the same way as interior particles, we let a = at, so an ≡ a. Evaluating Eq. (2.16) at steady state for the leader gives v = al/γ. Substituting this expression for v into Eq. (2.16) and evaluating at steady state for a follower gives g+(d) = a l−a. The value of d is determined by the number of intersection points between the function g+(x) and the horizontal line at the value a l − a. There could be zero, one, two, or more intersection points depending on the functional form of g+(x), implying that the number of school formations that exist could be zero, one or more. As we shall illustrate later, not all solutions are stable. The stability depends on the slope of g+(x) at each intersection point. Note that if all particles have identical autonomous self-propulsion, i.e. ai = a, then fi = g+(d) = 0 when the school is formed, implying that 30 2.3. Schools in One Dimensional Space all schooling forces vanish at the steady state perfect school. Only when the particles deviate from their ‘preferred’ positions in the school would the schooling forces become nonzero. Then the effect of these forces brings the particles back to their correct positions (if stable) or amplifies the deviation (if unstable). A school in which schooling forces vanish at steady state will be denoted a stress-free school. However, if the leader has a distinct self- propulsion, i.e. al 6= a, then fi = al − a 6= 0 which means that the schooling forces do not vanish when the school is formed. We shall call such a school a stressed school. 2.3.2 Schools Formed by Interactions with Two Nearest Neighbours. Consider a 1D school with particles labeled as in Fig. 2.2, such that a1 = a l, ai = a for i = 2, . . . , n − 1 and an = at. Now, we assume that each particle interacts with its nearest neighbours ahead and behind it, giving the following interaction forces: f1 = −g−(|x2 − x1|) + h−(|v2 − v1|) sgn(v2 − v1), (2.18) fi = g+(|xi−1 − xi|)− g−(|xi+1 − xi|) + h+(|vi−1 − vi|) sgn(vi−1 − vi) + h−(|vi+1 − vi|) sgn(vi+1 − vi), i = 2, · · · , n− 1, (2.19) fn = g+(|xn−1 − xn|) + h+(|vn−1 − vn|) sgn(vn−1 − vn). (2.20) When a school is formed, vi = v and xi−1 − xi = d (and thus xi+1 − xi = −d) for all i, which implies that f1 = −g−(d), (2.21) fi = g+(d)− g−(d), (i = 2, · · · , n− 1), (2.22) fn = g+(d). (2.23) Next, we substitute these forces into the steady-state equation (2.14) which gives γv = al − g−(d), (2.24) γv = a + g+(d)− g−(d), (2.25) γv = at + g+(d), (2.26) 31 2.3. Schools in One Dimensional Space where an = a t denotes the autonomous self-propulsion of the tail particle. Solving these equations for v, al, and at we obtain v = 1 γ [al + at − a], (2.27) al = a+ g+(d), (2.28) at = a− g−(d), (2.29) which leads to the following conclusions: 1. The speed of the school depends explicitly only on the self-propulsion forces al, at and a and γ, among which al and a are parameters with fixed values. Eq. (2.28) determines the value of d as well as the number of possible school solutions based on the given form of g+(d). Once the value of d is determined, Eq. (2.29) determines the value of at, which is the self-propulsion force the tailing particle must adopt in order to keep up with the school. 2. Individuals at the front or rear of the group have fewer interactions. Thus, to keep a fixed nearest-neighbour distance, some compensation for missing forces is needed. Otherwise, it would be impossible to maintain a perfect school. We note that in many biological examples, which deviate significantly from our perfect school assumptions, there is crowding at the edges of a herd or swarm. Examples of this type occur in herds of wildebeest shown in [19] and in locust hopper bands [24]. Both examples show a gradual increase of density of individuals towards a front edge, and sharp boundary at that edge. In the case of interactions with both nearest neighbours, a stress-free steady-state solution implies that g+(d) = g−(d) = 0. From Eqs. (2.27)– (2.29), the stress-free condition implies that at = al = a, and v = a/γ at steady state; that is, each individual must have the same autonomous self-propulsion, and the school velocity is dependent only on the common autonomous self-propulsion and the friction coefficient. Alternatively, for a stressed solution with nearest-neighbour distance d and v > 0, (2.27)–(2.29) imply that al > g−(d). (2.30) The school velocity is given by v = [al − g−(d)]/γ, and thus (2.30) assures that v > 0. Further, Eq. (2.30) means that the leading particle must generate 32 2.3. Schools in One Dimensional Space a sufficiently large autonomous self-propelling force to overcome attraction (or effective drag) due to the particle behind it, while the difference between the two forces determines the school speed. Eq. (2.28) means that the differ- ence between al and a together with the functional form of g+(x) determine the number, and values of allowable group spacing d. For fixed al and a, the autonomous self-propulsion of the tail particle must be adjusted to satisfy Eq. (2.29). We note that a stressed school can occur even if al and a have different signs; i.e., if the self-propulsion of the leading particle and interior particles are in opposite directions. In this case, the schooling forces over- power the intrinsic motion of the interior particles. In this sense, interior individuals are all followers, with the tail particle a special type of follower whose autonomous self-propulsion is distinct from that of the other follow- ers. We note that Eqs. (2.28)–(2.29) place a further restriction on g±(x): that al − a ≤ maxx(g+(x)), and that a− at ≤ maxx(g−(x)) to guarantee the existence of at least one perfect school solution. As one example, suppose that a = 0 for all interior particles i.e., indi- viduals base responses on interactions alone, instead of intrinsic autonomous forces. In such a school, to keep a nearest-neighbour distance d, at = −g−(d) for the tail particle. The spacing distance d is determined from (2.28), i.e., al = g+(d). Then for a monotonic increasing g+(x), increasing a l results in a larger nearest-neighbour distance d. However, because steady-state velocity is v = [al − g−(d)]/γ, increasing al does not necessarily increase v. Instead, Eq. (2.27) implies that v = 1 γ (g+(d)− g−(d)) . (2.31) Eq. (2.31) illustrates the relationship between school velocity and spac- ing. Specifically, if the slope of g+(x) − g−(x) is positive, the distance d between neighbours in the school increases as school velocity v increases. On the other hand, if the slope of g+(x) − g−(x) is negative, d decreases as v increases. Examples of these two different cases are shown in Fig. 2.5. This result is immediately applicable to experimental observations. For ex- ample, if the NND d increases as the school moves with a higher speed, there is reason to believe that the case in Fig. 2.5.a might be occurring. It is almost impossible to determine the functional form of g±(x) for re- alistic animals, but by observing the relationship between the changes in school speed and the NND, we can comment on the nature of the slope of g+(x) − g−(x). Note that a = 0 is not required for this result; in the case 33 2.3. Schools in One Dimensional Space a 6= 0, v = a/γ + (1/γ) (g+(d)− g−(d)) , which is identical to Eq. (2.31), except for a shift in velocity of a/γ. 2.3.3 Examples of Some Schooling Forces Distance-dependent forces. The theory developed so far applies to any functions g±(x), h±(x). However, requiring realistic schooling forces restricts our choice of functions. For distance-dependent forces, we take as a specific example Hill functions, i.e., saturating monotonic increasing functions of x: g(x) = c [ (1 + km/dm0 )x m xm + km − 1 ] , (m ≥ 1), x ≤ xf . (2.32) m > 0 and d0 > 0 represent the steepness and x-intercept of g(x), respec- tively. Note that g(0) = −c and c[1+km/dm0 ] are, respectively, the minimum and the saturation values of g(x). Eq. (2.32) defines a four-parameter family of functions from which g+(x) and g−(x) might be chosen. A decay beyond some finite distance can represent a limited sensing range (0 ≤ x ≤ xf ) of the individuals. Fig. 2.3 shows a typical distance-dependent schooling force based on a Hill function. Other typical forces used include exponential forces of the form g(x) = A exp(−x/a)−R exp(−x/r), (2.33) or inverse power forces g(x) = A xa − R xr , where A,R, a, and r are parameters (see [11]). Velocity-dependent forces. In one dimension, the velocity-dependent forces h±(v) have a simple interpretation. In this case, all individuals have the same absolute direction (though they may differ by a sign, depending on whether movement is along the axis in the positive or negative direc- tion). Then the velocity-dependent force given by Eq. (2.12) summed over influencing neighbours, ĥ±(|vji|) vji|vji| , can be represented as a single function of vji, since in 1D the direction term vji/|vji| is simply sgn(vji), i.e., ±1. We thus define h±(v) = ĥ±(|v|) v|v| , 34 2.3. Schools in One Dimensional Space Figure 2.3: An example of a distant-dependent force in the form of a Hill function with a decay past xf , to depict a finite sensing range. where h(v) is defined for all v ∈ R. There are a number of cases to consider when determining reasonable forms of h(v). From Fig. 2.4, the influence of a Figure 2.4: Four cases of neighbour interaction: [a] neighbour j travels faster than i in the same direction, (vji > 0), [b] neighbour j travels slower than i in the same direction, (vji < 0), [c] neighbour j diverging from i, (vji > 0), [d] neighbour j converging to i, (vji < 0). neighbour j on the velocity vi of individual i should cause acceleration when vji > 0, but deceleration when vji < 0. We assume that there is no influence 35 2.3. Schools in One Dimensional Space when vi = vj , i.e., h(0) = 0. We will later show that under the assumption h(0) = 0, h(v) has no effect on existence of a perfect school solution. Thus, a reasonable form of h(v) is a monotonically increasing function with h(0) = 0. This type of assumption is incorporated in models for swarming in 2D and 3D, discussed in [12, 15, 16, 17, 20] as “arrayal forces”, since these forces act to array individuals in parallel. As an example, the function we use in this paper is h(v) = C tanh(v/ρ), where C > 0 and ρ are parameters governing the amplitude and steepness of h(v), respectively. If ρ < 0, h(v) is monotonically decreasing; we include this generality, but we will later show that only the monotonically increasing case (ρ > 0 in this example) enhances stability of the perfect school solution. As in the distance-dependent functions, particular choices of parameters lead to different shapes of functions for h+(v) and h−(v). In [8], h(v) = αv (α > 0 is a constant) in the case when only nearest neighbour interactions are considered. The number of reasonable choices of the functions g±(x) and h± that yield qualitatively identical results as the expressions we adopted in this study is infinite. Any information that leads extra restrictions of these functions would be valuable for making our choices more realistic. 2.3.4 The Stability of a School Solution We next explore the stability of a 1D perfect school. To do so, we consider a perfect type I school, moving with speed v and nearest-neighbour distance d: vsi = v and x s i = x 0 1 + vt − (i − 1)d (where x01 is the position of the leading particle at t = 0). This corresponds to the steady-state solution psi = x 0 1 − (i− 1)d and qsi = 0 in the moving frame coordinates. Using the moving-frame analogues of Eqs. (2.18)–(2.20) and (2.12), and the expression for v from Eq. (2.27) in Eqs. (2.13)-(2.14), we obtain the moving-frame equations of motion for nearest-neighbour interactions ahead and behind individuals: dpi dt = qi, (2.34) dq1 dt = al − g−(|p2 − p1|) + h−(|q2 − q1|) sgn(q2 − q1) 36 2.3. Schools in One Dimensional Space Figure 2.5: [a] An example of interaction forces where the net force has positive slope at d. Here, g+(x) = 2[3x 2/(x2 +42)− 1], and g−(x) = 3x/(x+ 4)− 1. [b] An example of interaction forces where the net force has negative slope at d. Here, g+(x) = (1/2)[3x 3/(x3+23)−1], and g−(x) = 3x/(x+8)−1. −γq1 − (al + at − a), (2.35) dqi dt = a+ g+(|pi−1 − pi|)− g−(|pi+1 − pi|) +h+(|qi−1 − qi|) sgn(qi−1 − qi) +h−(|qi+1 − qi|) sgn(qi+1 − qi)− γqi − (al + at − a), (2.36) dqn dt = at + g+(|pn−1 − pn|) + h+(|qn−1 − qn|) sgn(qn−1 − qn) −γqn − (al + at − a). (2.37) 37 2.3. Schools in One Dimensional Space To investigate the linear stability of the school solution, we introduce a small perturbation to the equilibrium values of the velocity and the position of each particle in the moving frame: pi(t) = p s i + δi(t), (2.38) qi(t) = q s i + ωi(t), (2.39) where ṗsi = 0, and q̇ s i = 0 for all i. Substituting Eqs. (2.38) and (2.39) into Eqs. (2.34)–(2.37), and expanding nonlinear terms in a Taylor series up to first order, we obtain dδi dt = ωi, (2.40) dω1 dt = −g−(d)− g′−(d)[δ1 − δ2]− h′−(0)[ω1 − ω2] −γω1 + (a− at), (2.41) dωi dt = g+(d)− g−(d) + g′+(d)[δi−1 − δi]− g′−(d)[δi − δi+1] +h′+(0)[ωi−1 − ωi]− h′−(0)[ωi − ωi+1] −γωi + (2a− al − at), (2.42) dωn dt = g+(d) + g ′ +(d)[δn−1 − δn] + h′+(0)[ωn−1 −ωn]− γωn + (a− al), (2.43) where i = 1, · · · , n in Eq. (2.40) and i = 2, · · · , n − 1 in Eq. (2.42) and where we have used h±(0) = 0. Using Eqs. (2.28)–(2.29), we can reduce these equations to the following simplified form governing perturbations: dδi dt = ωi, (2.44) dω1 dt = −g′−(d)[δ1 − δ2]− h′−(0)[ω1 − ω2]− γω1, (2.45) dωi dt = g′+(d)[δi−1 − δi]− g′−(d)[δi − δi+1] +h′+(0)[ωi−1 − ωi]− h′−(0)[ωi − ωi+1]− γωi, (2.46) dωn dt = g′+(d)[δn−1 − δn] + h′+(0)[ωn−1 − ωn]− γωn, (2.47) 38 2.3. Schools in One Dimensional Space where i = 1, · · · , n in (2.44) and i = 2, · · · , n− 1 in Eq. (2.46). The eigenval- ues of the coefficient matrix associated with this linear system give stability information. It is not straightforward to obtain explicit expressions for the eigenvalues in the general case of Eqs. (2.44)–(2.47). However, when we simplify the system further for the case of following only one immediate neighbour in front, we obtain the system dδi dt = ωi, (2.48) dω1 dt = −γω1, (2.49) dωi dt = g′+(d)[δi−1 − δi] + h′+(0)[ωi−1 − ωi]− γωi, (2.50) dωn dt = g′+(d)[δn−1 − δn] + h′+(0)[ωn−1 − ωn]− γωn. (2.51) When written in matrix form, the coefficient matrix corresponding to Eqs. (2.48)–(2.51) is a 2 × 2 block matrix with n × n blocks. To obtain the eigenvalues of this matrix we first simplify the determinant equation using the following result from [18]: given the block matrix [ A B C D ] , where A,B,C, and D are m×m matrices, if A and B commute, then det [ A B C D ] = det [DA−CB] . (2.52) This result allows us to reduce the dimension of the matrix in the determinant formula for the eigenvalues. Assuming only forward (or, similarly, backward) sensing allows a further reduction to an upper (lower) triangular matrix, for which the determinant can be obtained from the product of the diagonal entries. The details of this calculation are given in Appendix A, where it is shown that the associated 2n eigenvalues are λ1 = 0, (2.53) λ2 = −γ, (2.54) λ± = −(h′+(0) + γ)/2± √ [(h′+(0) + γ)/2]2 − g′+(d), (2.55) 39 2.3. Schools in One Dimensional Space where eigenvalues λ+ and λ− each have multiplicity n−1. The zero eigenvalue in Eq. (2.53) characterizes the translational invariance of the solution, while λ2 in Eq. (2.54) is always negative for γ > 0. However, in the absence of damping (γ = 0), the multiplicity of the zero eigenvalue is 2, indicating that the school solution is neutrally unstable. From Eq. (2.55), we obtain the stability conditions (guaranteeing λ± < 0) h′+(0) + γ > 0, (2.56) g′+(d) > 0. (2.57) Therefore, the contribution of the speed-dependent force is stabilizing if h′+(0) > 0, and destabilizing if h ′ +(0) < 0. However, even if the speed- dependence is destabilizing (i.e. if h′+(0) < 0), the solution can still be stable if the damping, γ, is large enough to compensate in Eq. (2.56). This suggests that a non-zero damping term is indispensable for ensuring stability of a perfect school solution. We have already shown that conditions for existence of the perfect school are independent of velocity-dependent interaction forces; condition (2.56) also illustrates that h+(v) is not essential for stability of the school solution. If h+(v) = 0, i.e., there is no velocity-dependent force at all, then λ± = −γ/2 ± √ (γ/2)2 − g′+(d). Then, assuming Eq. (2.57) holds, the school solution remains stable if γ > 0. On the other hand, the distance-dependent force is crucial, both for exis- tence (as previously shown) and for stability of a perfect school solution. If g+(x) = 0, λ± = 0, −(h′+(0) + γ). Recalling that the multiplicity of λ+ is n − 1, the school solution becomes neutrally unstable in the absence of the force g+(x), as the multiplicity of the zero eigenvalue is then n. Also note that the system can exhibit an oscillatory approach to steady state in the presence of complex conjugate eigenvalues; i.e., when [(h′+(0) + γ)/2] 2 < g′+(d). The results on school formation in 1D reveal that damping and the pres- ence of a distance-dependent schooling force are key factors for the existence of a stable school pattern. Velocity-dependent schooling forces can make the school solution more or less stable depending on the sign of h′+(v), as shown in Eq. (2.56). We have shown that the school speed is determined by the autonomous self-propelling forces of particles at different locations within the school and is independent of the schooling forces. We also showed that school formation occurs even when each particle follows only its neighbour immediately ahead. This ‘follower’ school presents one minimal model of per- fect school formation in 1D, and its simplicity has permitted an analytical 40 2.3. Schools in One Dimensional Space description of the existence and stability of such schools. 2.3.5 Numerical Simulations in One Dimension Figure 2.6: 1D simulations showing position in a moving frame moving at velocity v as indicated. Interaction forces are as in [a] Fig. 2.5.a, and [b] Fig. 2.5.b. Self propulsion terms are [a.i] al = 0.5, [a.ii] al = 0.7, [b.i] al = 0.9, [b.ii] al = 0.8, while a = 0.1 in all cases. In [a], higher school speed resulted in an increase in NND from 3.26 in [a.i] to 3.49 in [a.ii]. The opposite occurred in [b], where higher school speed caused a decrease in d from 3.73 in [b.i] to 3.175 in [b.ii]. We simulate the one-dimensional model (2.10)–(2.11) in MATLAB using the solver ODE45, which integrates a system of ODEs with an embedded 41 2.3. Schools in One Dimensional Space Runge-Kutta method. The algorithm is based on the Dormand-Prince 4(5) formulae, which uses a 4th and 5th order pair of methods. We show the case of distance-dependent interactions only; i.e., h± = 0. Data is plotted in the moving frame, relative to the average velocity at t = tend. Particles are initialized randomly in the interval [0,11]. The model inputs g±(x), a l and a are chosen to satisfy Eq. (2.28), which in turn determines d. When d is known, at follows from Eq. (2.29), giving the self-propulsion that the last particle must produce to keep up with the school. The absolute magnitudes of model inputs were chosen sufficiently small to avoid large changes in model variables at each time step. Relative magnitudes of model inputs were cho- sen so that motion was not dominated by any of the forces of autonomous propulsion, interaction, or drag. We test two distinct examples of interaction forces to illustrate our results. In case 1 (Fig. 2.6.a) we use forces as in Fig. 2.5.a, whereas in case 2 (Fig. 2.6.b), the forces are as in Fig. 2.5.b. Increasing al in both cases generates a larger distance d, but, as noted in our analytic results, we expect velocity to increase in case 1 and decrease in case 2 as d is increased. These results are borne out by simulations. Comparing results presented here to those obtained in other models in which all-to-all interactions were considered [1, 5, 8, 11], there are several important differences. First, the NND is determined by the self-propelling forces al, a and the interaction force function g+(x) (see Eqs.(2.28-2.29). It is independent of the number of particles in the group. This is because only nearest-neighbour interactions are considered here. Increasing the number of particles does not increase the forces each particle experiences. Second, the “edge effect” that was typically observed in the other models is eliminated by making the edge particles behave differently based on their location in the group. Perfect schools cannot occur if the edge particles are not capable of sensing their locations in the group and behaving accordingly. Third, the interaction functions that we use here are monotonic functions that typically violate the so called H-stability conditions studied in [1, 5]. H-stability is characterized by the conditions under which the NND approaches a finite non-zero value as the number of particles approaches infinity in an all-to-all coupled system. The conditions do not apply to the results presented here as only nearest-neighbour coupling is considered. Fourth, in a perturbed perfect school solution, the rate at which the system evolves toward the stable school pattern is governed by the real parts of the eigenvalues of the linearized system given by Eqs.(2.48-2.51). These eigenvalues are determined by the 42 2.4. The Soldier Formation in 2D Space Figure 2.7: An example of an approximate soldier formation in Atlantic bluefin tuna, courtesy of Dr. M. Lutcavage, Large Pelagics Research Center, UNH. parameter γ and the slopes of the interactions forces g′±(d) and h ′ ±(0). These eigenvalues are explicitly given by Eqs.(2.53-2.55) when each particle follows only the immediate neighbour to its front. The rate at which a perfect school solution is approached from an arbitrary initial condition is determined by both the linear and nonlinear effects of interaction forces. 2.4 The Soldier Formation in 2D Space Among the numerous types of schooling behaviour found in higher dimen- sions in nature is the soldier formation. Such schools are characterized by a roughly linear organization of individuals with a common heading in some (nonaxial) direction. A school of this type formed by bluefin tuna is shown in Fig. 2.7. In the following treatment, we consider an idealization of soldier formations found in nature, such that all individuals have the same heading and velocity, and whose positions lie on a straight line (see Fig. 2.8). The simple linear structure of such a school facilitates our analysis. We apply the frontal nearest-neighbour interaction rule (i.e., g− = 0) to a soldier formation of n particles, labeled so that 1 is the leader, and the remaining n−1 are followers. Under the chosen interaction regime, ~x1 is the 43 2.4. The Soldier Formation in 2D Space Figure 2.8: A schematic diagram of individuals in soldier formation with frontal nearest-neighbour detection. Grey arrows indicate the direction of schooling force, while black arrows indicate direction of motion. nearest neighbour of ~x2, ~x2 is the nearest neighbour of ~x3 and so on. Note that the coupling of particles is not predetermined, but instead is determined by the group geometry. We neglect velocity-dependent interactions (i.e., h± = 0). The corresponding equations of motion for the system are d~x1 dt = ~v1, (2.58) d~v1 dt = ~al − γ~v1, (2.59) and d~xi dt = ~vi, (2.60) d~vi dt = ~a− γ~vi + ~fi, (2.61) for i = 2, . . . , n. The schooling force ~fi depends on the relative position of the ith individual with respect to its nearest neighbour, (i− 1), and thus has arguments ~fi = ~fi(~xi−1 − ~xi). Because we assume that individuals sense only what is in front of them, we have ~f1 = ~0 for the leader. 44 2.4. The Soldier Formation in 2D Space 2.4.1 Existence Conditions for the Soldier Formation In this section, Eqs. (2.58)–(2.61) are used to obtain existence conditions for the soldier formation. In the soldier formation d~vi/dt = 0, and so, we consider the formation as a ‘moving school’ steady state for the system (henceforth referred to simply as steady state). We note that this steady state differs from a full steady state, wherein also d~xi/dt = 0, i.e., where the aggregate is stationary, rather than moving. The equation of motion (2.59) for the leader at steady state implies ~v = ~al γ , (2.62) where ~v is the steady-state velocity of the leader. Note that because every individual moves at the same velocity in a steady-state school solution, ~v is also the velocity of all individuals in that school. Substituting this into the steady-state equation of motion for the ith individual gives the existence condition for a soldier formation: ~f si = ~a l − ~a, i = 2, . . . , n. (2.63) We express ~f si as a unit direction vector (oriented along the axis of the school at steady-state) and a magnitude as in (2.5), giving ~f si = g+(d)~u, where ~u = ~xi−1,i/|~xi−1,i|, and d = |~xi−1,i| is the distance between individuals at steady state. We consider two cases: first, if ~al 6= ~a, then by (2.63), g+(d)~u = ~al − ~a, i.e., the direction of axis of the school, ~u, must be parallel to ~al − ~a. Similar to the case of a 1D school, the magnitude of ~al − ~a determines the value(s) of g+(d) that guarantee the existence of such a school formation, i.e., g+(d) = |~al − ~a|. (2.64) For a monotonic function g+, this uniquely determines the inter-individual distance d along the axis of the school at steady state. This is an example of a stressed school in which interaction forces do not vanish at steady state. Once the parameters ~al and ~a are fixed and g+(x) is specified, ~u, d, ~v, and the bearing angle are all determined. On the other hand, if ~al = ~a, then at steady state, ~fi = ~a l − ~a = ~0, i.e., the school is stress-free and there is no restriction on the bearing angle of the axis of the school. 45 2.4. The Soldier Formation in 2D Space 2.4.2 Stability Analysis In the soldier formation, the cohesiveness of the group can be examined in terms of the stability of the system to small perturbations about the steady state; i.e., the local stability. The stability analysis of the soldier formation in 2D follows closely the analysis of 1D schools. The geometry of a soldier formation allows some simplification in analysis by a judicious choice of orientation of axes. Specifically, we align our coor- dinate system so that the x-axis lies along the line of individuals. Then the y−components of the particle positions are zero at the steady-state soldier formation solution. Note, however, that unlike 1D motion, the individuals here move (in general) in some direction other than the x-axis. Thus we write ~xs1 = (0 , 0), ~x s 2 = (d , 0), . . . , ~x s i = ((i− 1)d , 0) , . . . , ~xsn = ((n− 1)d , 0). Henceforth, we display only calculations on the ith equation, but this implicitly represents all i = 1, . . . , n. We transform to a frame of reference moving with the school velocity at steady state (Eq. (2.62)). We let ~xi = ~pi + ~vt, ~vi = ~qi + ~v. We substitute these expressions into (2.58)–(2.61) to obtain d(~p1 + ~vt) dt = ~q1 + ~v, (2.65) d(~q1 + ~v) dt = ~al − γ (~q1 + ~v) , (2.66) d(~pi + ~vt) dt = ~qi + ~v, i = 2, . . . , n (2.67) d(~qi + ~v) dt = ~a− γ (~qi + ~v) + ~fi, i = 2, . . . , n (2.68) Simplifying these equations, noting that γ~v = ~al, we obtain the moving frame equations of motion: d~p1 dt = ~q1, (2.69) 46 2.4. The Soldier Formation in 2D Space d~q1 dt = −γ~q1, (2.70) d~pi dt = ~qi, i = 2, . . . , n, (2.71) d~qi dt = ~a− ~al + ~fi − γ~qi, i = 2, . . . , n (2.72) The arguments of ~fi are now ~fi = ~fi(~pi−1 − ~pi), though we continue to omit these arguments unless needed for clarity. To summarize, ~pi represents the position of individual i in the moving frame, and ~qi represents the velocity of individual i in the moving frame. The soldier formation is a steady-state solution to this system, and thus at steady-state, ~qi = 0 and ~a−~al + ~fi = 0, which gives the existence condition in (2.63), and ~psi = ê1(i − 1)d, where ê1 is the unit vector (1 , 0). Next, we perform a stability analysis on our moving-coordinate 2D variables ~pi and ~qi. We introduce small perturbations in position ~δi and velocity ~ωi about the steady state; i.e., ~pi(t) = ~p s i + ~δi(t) = ê1(i− 1)d+ ~δi(t), and ~qi(t) = ~q s i + ~ωi(t) = ~ωi(t), where ê1 is a unit vector in the x direction. Substituting these perturbations into Eqs. (2.69) – (2.72) and noting that time derivatives of steady-state variables are 0, we get d~δ1 dt = ~ω1, (2.73) d~ω1 dt = −γω1, (2.74) d~δi dt = ~ωi, i = 2, . . . , n, (2.75) d~ωi dt = (~a− ~al) + ~fi − γωi, i = 2, . . . , n. (2.76) We proceed by expanding the nonlinear term, ~fi in a Taylor series. We first rewrite ~fi as (fix, fiy) = ~fi(~zi) = g+(|~zi|) ~zi|~zi| = g+(|~zi|) ( zix |zi| , ziy |zi| ) , 47 2.4. The Soldier Formation in 2D Space where 〈zix, ziy〉 ≡ ~zi = ~pi−1 − ~pi = −ê1d+ ~δi−1 − ~δi, (2.77) and |~zi| = √ z2ix + z 2 iy. That is, ~fi is oriented in the direction of the vector connecting two neigh- bouring individuals, with a magnitude dependent on the distance between the two individuals. The Taylor expansion of ~fi is ~fi = ~fi(−ê1d) + ∂ ~fi ∂~zi |~zi=~zsi · [~zi − ~zsi ] + . . . , (2.78) where the expression in braces is a Jacobian matrix, i.e., the partial deriva- tives of fix and fiy with respect to zix and ziy coordinates. The Jacobian matrix in Eq. (2.78) evaluated at steady state is given by Df ≡ ∂ ~fi ∂~zi |~zi=~zsi = [ g′+(d) 0 0 g+(d)/d ] . Details of this calculation are contained in Appendix B. Df is independent of the equation index i, and is thus identical for all individuals i = 2, . . . , n. Substituting Df into the Taylor expansion of ~fi in Eq. (2.78), and noting from Eqs. (2.77) and (A.8) (in Appendix) that ~zi − ~zsi = ~δi−1 − ~δi, (2.79) we now have ~fi = ~fi(−ê1d) +Df · [~δi−1 − ~δi] + . . . . (2.80) Furthermore, because at steady state we have ~fi(−ê1d) = ~al − ~a, it follows that for the small perturbation, ~fi ≃ ~al − ~a+Df · [~δi−1 − ~δi]. We now use this linear approximation for ~fi in Eq. (2.76), resulting in d~ωi dt = ~a− ~al + (~al − ~a+Df · [~δi−1 − ~δi])− γ~ωi. 48 2.4. The Soldier Formation in 2D Space Canceling out terms and combining with Eqs. (2.73)–(2.75) gives the follow- ing system: d~δ1 dt = ~ω1, (2.81) d~ω1 dt = −γω1, (2.82) d~δi dt = ~ωi, i = 2, . . . , n (2.83) d~ωi dt = Df · [~δi−1 − ~δi]− γ~ωi, i = 2, . . . , n. (2.84) The matrix expression of this system is a 2 × 2 block matrix with 2n × 2n blocks. The eigenvalue equation is formulated in determinant form, and the commutativity of the blocks of the associated matrix allows for a simplifica- tion of the determinant using Eq. (2.52). The eigenvalue calculation is thus reduced to calculating the determinant of an n × n block-triangular matrix with 2× 2 blocks, which is given by the product of the determinants of the diagonal blocks. The matrix form of Eqs. (2.81)–(2.84) and eigenvalue cal- culations are given in Appendix B. The coefficient matrix associated with Eqs. (2.81)–(2.84) has eigenvalues λ = 0, with multiplicity 2, (2.85) λ = −γ, with multiplicity 2, (2.86) λ = −γ 2 ± √ γ2 − 4g′+(d) 2 , each with multiplicity n− 1 (2.87) λ = −γ 2 ± √ γ2 − 4g+(d)/d 2 , each with multiplicity n− 1. (2.88) The first eigenvalues, λ = 0 with multiplicity 2, indicates a translational invariance in two-dimensional space. The next distinct eigenvalue, λ = −γ is always negative when damping is present, γ > 0. This further indicates the importance of a drag term for stability of solutions. The remaining eigenvalues have negative real parts provided both g′+(d) > 0, (2.89) g+(d) > 0. (2.90) 49 2.4. The Soldier Formation in 2D Space Figure 2.9: Numerical simulations in 2D for [a] 8 particles with ~al = [0.1, 0.1], ~a = [0,−0.2] and [b] 6 particles with ~al = [0.1, 0.1], ~a = [0, 0]. Trajectories are plotted through time. Note in [a] that the solution line connecting individuals is oriented in the direction of ~al − ~a = [0.1, 0.3], while the school velocity at steady state is ~v = [0.2, 0.2]. In [b], both the solution line and school velocity are in the same direction, with ~al = [0.1, 0.1], ~v = [0.2, 0.2]. Interpreted in terms of stability of the school solution, Eqs. (2.89)–(2.90) implies that the schooling force acting to maintain the school must be at- tractive at the NND d. Additionally, this attractive force should increase if individuals with spacing d try to move further apart. 2.4.3 Numerical Simulation of the Soldier Formation We perform simulations in two dimensions using single-neighbour interac- tions, with neighbour detection only in the frontal plane (g− = 0). We include no velocity-dependent interactions (h± = 0). We choose an exponen- tial interaction function following Eq. (2.33), as g+(x) = 0.5(exp(−x/50) − 2 exp(−x/1)). Eqs. (2.58)–(2.61) are evolved using ODE45 as in the 1D simulations. Initial positions are assigned randomly in [0, 4]× [0, 4], and we set γ = 0.5. We first test two distinct cases of ~al and ~a values. In case 1 (Fig. 2.9.a) ~al 6= ~a 6= ~0, while in case 2 (Fig. 2.9.b), ~al 6= ~a = ~0. As noted in our analytical results, we expect the soldier formation to move in the direction of ~al, while the solution line connecting individuals is oriented in the direction of ~al −~a. In case 1, this corresponds to a school directed non-axially to the solution line, while in case 2, the solution line and school direction coincide, leading 50 2.4. The Soldier Formation in 2D Space to a trail-following group. The NND in these solutions is determined by Eq. (2.64). For the cases in Fig. 2.9a and 2.9.b, d = 1.7949 and 1.055 respectively. The difference in d between the two cases is due to the difference in the value of |~al − ~a| but not the difference in the number of particles. Because only nearest-neighbour interactions are considered here, this distance remains unchanged as the number of particles is increased. These analytical results are confirmed by simulations shown in Fig. 2.9. Although our stability results are only local, the simulations suggest that the soldier formation is a globally stable pattern for the specific choice of parameter values in Fig. 2.9a. In a follow-up study [10], we found that milling formation could coexist with the soldier formation when the force function g+(x) is bi-phasic. In that case, different initial conditions determine which pattern will eventually emerge. Unlike the soldier formation, the existence of a milling solution and the NND both depend on the number of particles for a fixed choice of the force function. In Section 4.1, we mentioned the case when al = a, which implies no re- striction on the bearing angle of individuals. Eq. (2.64) implies that the zero of g+(x) gives the individual distance d at steady state for such groups. Fig. 2.10 shows examples of such configurations. Note that individual distance at steady state remains constant at d = 0.7, though the relative angles between neighbours is determined by initial conditions. In previous simulations, we used interaction functions that allowed stable schools to arise. However, if we choose g+(x) so that g ′ +(x) is negative for all x > 0, then we fail to satisfy the stability condition in Eq. (2.89). We explore the motion of particles under such conditions in simulations (Fig. 2.11). In these simulations, particles tend to organize in moving pairs which remain together due to close-range attraction. The path traced out by these pairs appears to be very sensitive to initial conditions, and does not follow a predictable pattern. Further simulations with such ‘stability-breaking’ inter- action forces reveal the presence of numerically stable mill solutions in which particles rotate around an invariant circle. These solutions require a closed curve of interaction connections between particles in order to evolve. The existence and stability of such mill solutions form the topic of Chapter 3. 51 2.4. The Soldier Formation in 2D Space Figure 2.10: Numerical simulations in 2D, for [a] 10 particles and [b] 25 par- ticles with ~al = ~a = [0.1, 0.1], and g+(x) = 0.5(exp(−x/50) − 2 exp(−x/1)). Trajectories are plotted through time, and grey arrows indicate interactions between nearest neighbours. Note that individuals have constant spacing between nearest neighbours, but are not restricted to a particular relative angle with nearest neighbours. Figure 2.11: Numerical simulations in 2D with g+(x) = 1 − x/2, which corresponds to (rather non-physical) short-range attraction and long-range repulsion. Trajectories are plotted in time for [a] 4 particles evolved to t = 40, and [b] 6 particles evolved to t = 300. Individuals group in pairs (though due to overlapping, some pairs appear as a single particle), and these pairs follow unpredictable trajectories. 52 2.5. Schools in the Form of 2D Arrays 2.5 Schools in the Form of 2D Arrays We have thus far been primarily concerned with 1D and quasi-1D forma- tions of particles. However, extending these ideas to schools in the form of regular 2D arrays (i.e., a perfect school in 2D) is both a logical step in the mathematical analysis, and is more relevant to many observed biological aggregates. This extension to 2D brings a number of challenges to our ana- lytical approach. First of all, similar to the case of 1D schools in which the tail particle should be aware of its special position in the group and behave differently to maintain the school, particles located at the edges of a 2D array should also behave differently than particles in the interior. Confounding the issue further, the number of edge particles change as the specific shape of the array changes. Thus, it is very difficult to study 2D arrays in a general form since there exist many possible 2D shapes for a group composed of a large number of particles. For groups of small numbers of particles, some analytical insights can be made, though these are shape-dependent and not generalizable. However, numerical simulations can provide insight into the mechanisms of forming schools of 2D arrays that tile the plane in a regular fashion. If we restrict our discussion to the models in which each particle only follows one closest neighbour, we must apply the restriction ~al = ~a, as otherwise a soldier formation would emerge (assuming the corresponding stability conditions are satisfied). For a group of particles that follow only one nearest neighbour in its frontal visual field with ~al = ~a 6= ~0, many stress-free patterns can occur depending on the initial conditions. The NND remains fixed in all these patterns since it is determined by the value of d that makes g+(d) = 0 (note that g+(x) has only a single zero in our simulations). Two examples are shown in Fig. 2.10. For the g+(x) function chosen in these simulations, d = 0.707 in both panels. Since interactions only occur between nearest neighbours, this distance does not change as the number of particles increases. The bearing angle between neighbours is sensitive to initial conditions, such that small changes in initial conditions can lead to large changes in the shape of the group that evolves. To generate a perfect school in 2D we could initialize particles in a regular array, and this array would persist in time so long as it is not subject to noise. However, in both irregular and regular arrays, perturbations can alter the configuration after such groups have been formed. This is due to the fact that the magnitude of the interacting force is only 53 2.5. Schools in the Form of 2D Arrays distance-dependent, and thus there exists a circular arc of positions around a preceding particle which its follower can occupy while remaining equidistant from its neighbour. Then upon undergoing a perturbation, a particle is free to return to any point on this arc, and not necessarily to its original point in the undisturbed array. We also note that velocity-dependent forces act to align the group, but do not contribute to maintaining a fixed bearing angle. Then the question arises as to what minimal mechanisms can generate a stable 2D perfect school solution. We have found that by increasing both the number of neighbours sensed by each particle and the sensing region, we can generate schools that are stable in numerical simulations. Consider a group of particles, each interacting with two nearest neigh- bours. We no longer restrict interactions to the front of a particle, and in- stead assume individuals can sense neighbours in any direction. By removing the restriction of sensing in the frontal region, we eliminate the edge issues discussed earlier. We use a single force for all directions of interaction, given by g(x) = 0.5(exp(−x/50) − 2 exp(−x/1)), as in earlier examples. Particles are initialized in a hexagonal tiling, where each hexagon contains an individ- ual in the middle, such that the entire school can be broken into groups of particles forming equilateral triangles. This tiling is then perturbed, and this perturbed tiling is taken as the initial condition. Such a system was found in numerical simulations to evolve to a regular tiling of fixed interindividual distance d that moves as a rigid body (see Fig. 2.12). Stability analysis for such a solution is quite difficult using analytical techniques introduced here for two reasons: it is difficult to label individuals to generate structure in stability matrices, and due to the geometry, individuals can change which particle they interact with as a result of a small perturbation. This latter fact actually serves to stabilize the solution by maintaining proper distance between individuals even when a direct interaction may not exist, as a per- turbation closer to a neighbour would cause interaction with this neighbour and consequent return to the initial relative position. Simulations of 2D schools that yield a large variety of different school patterns with different shapes and features have been carried out in the presence of small noise, using various interaction schemes. Some groups form perfect schools, while some exhibit variation at the boundary of the school, both spatially and temporally. A detailed exploration of all of these numerically observed patterns lies outside the focus of this paper. 54 2.5. Schools in the Form of 2D Arrays Figure 2.12: Numerical simulation of two nearest-neighbour interactions. Trajectories are plotted through time. Note that the particles, initially perturbed, evolve to a regular tiling of the plane with individual distance d = 0.707. 55 2.6. Discussion 2.6 Discussion The formation of cohesive schools is a phenomenon that occurs in a number of aggregating organisms. Among the most striking and common examples are flocks of birds and schools of fish. These groups provide motivation for a host of mathematical models which seek to lend insight into the behavioral algorithms used to generate school patterns. In contrast to the Eulerian modeling approach (see, e.g., [21, 22, 23]) which gives mean field proper- ties and is useful in modeling the global dynamics of large populations, we use a Lagrangian approach. By modeling individual particles, we gain in- formation on spacing and heading that is lost in the Eulerian approach at the expense of greatly increasing the dimension of the associated system of ODEs. However, in the case of some particular solutions, these large systems of nonlinear ODEs permit analytical insights. In our study, we are concerned with relatively basic interactions between particles, and geometrically simple solutions that such interactions permit. Our modeling approach takes as inputs the inherent behaviour of individ- uals, as well as a functional description of how they interact with one another. With these inputs, the model can then address [a], whether or not the school can exist, [b], if it can exist, the spacing and school velocity that prevails at steady state, [c], the geometry of the school; that is, how individuals align with other individuals, and how they move in relation to this alignment, and [d], whether or not a school solution is linearly stable under small perturba- tions. The simplicity of our approach allows us to directly connect features of existence, stability and school geometry to basic features of interaction functions and model parameters. Our results do not qualitatively depend on the specific choice of the schooling functions g±(x) and h±(x), provided these functions satisfy the criteria established in this paper. In 1D, we studied schools formed by interactions to the front, and then, interactions both to the front and back. We obtained existence conditions for such schools to exist, and analytically obtained school velocity and indi- vidual distance d that prevail at the perfect school solution. We were able to relate how the slope of the “net schooling force function” g+− g− near d de- termines the relationship between velocity and spacing in a school, and this relationship was shown in simulations. We then obtained stability conditions in the case of interactions to the front, with the result that both h′+(0) + γ and g′+(d) be positive to guarantee stability. Using essentially the same analytical techniques as in 1D, we explored 56 2.6. Discussion soldier formations of particles organized in a line, moving in some direc- tion non-axial to the line. We obtained existence conditions which provided analytical descriptions of spacing, heading direction, school speed and the direction of the line of individuals. We performed a stability analysis on such schools, with the result that both g+(d) and g ′ +(d) must be positive to guarantee stability. Following our analysis, we performed numerical simu- lations of the model in 2D, showing that the soldier formation appears to be globally stable with respect to random initial conditions when both the existence conditions and local stability conditions are satisfied, and ~al 6= ~a. In the case of ~al = ~a, we find no restriction on relative angle between neigh- bours, and individuals arrange themselves with constant distance, but with relative angles dependent on initial conditions. We then explored situations in which local stability conditions were not met, which leads to individuals pairing off and each pair following an unpredictable trajectory. Lastly, we presented a minimal mechanism to generate perfect schools in 2D such that particles form a regular array that moves at a constant speed. This solution was shown to be stable in numerical simulations using a typical interaction function. We have shown analytically that one neighbour is sufficient to generate stable perfect schools in 1D, while two neighbours were shown numerically to generate stable perfect schools. This leads to a conjecture that interactions with three neighbours should suffice to generate stable perfect schools in 3D. Exploring such questions analytically and computationally forms the basis of our ongoing work. Although in this paper, we have studied only simple arrangements of par- ticles, our approach clearly demonstrates how the properties of each aspect of the model influences the observed features of the school. Although we have emphasized simplicity in formulating minimal mechanisms for school formation, our results are nevertheless general, in that we have avoided as- suming particular functional forms to obtain results. Also, the simplicity of our approach underscores the considerable variation in the class of solutions, and the sensitivity of these solutions to how particles interact. The study of the formation of social aggregates in self-propelled particles has become increasingly active. A number of recent studies have revealed a large number of interesting patterns in groups of particles coupled through pairwise interactions [1, 5, 8, 11]. In these studies, all-to-all interactions are considered and the limiting behavior of the group as the number of parti- cles approaches infinity has been shown to be crucial in determining some 57 2.6. Discussion important features of the aggregate patterns that can occur. The patterns that emerge in these systems are often not regular and may even be meta- stable or transient states [1]. Results presented in this work represent a rare distinct attempt at focusing on minimal mechanisms that are sufficient for generating regular school patterns that are dynamically stable. We showed that nearest-neighbour interactions are sufficient for generating regular school patterns with some crucial characteristic features that are independent of the number of particles in the group. A detailed analysis of the similarity and difference between the aggregate patterns in all-to-all coupled and nearest- neighbour coupled self-propelling particles could lead to experimental criteria that allow us to check which one is closer to reality. Acknowledgements This work was funded by NSERC discovery grants to Y.X. Li and L. Edelstein-Keshet. R. Lukeman was funded by an NSERC graduate scholarship, and by a UBC University Graduate Fellowship. 58 Bibliography [1] Y. L. Chuang, M. R. D’Orsogna, D. Marthaler, A. L. Bertozzi, and L. S. Chayes, State transitions and the continuum limit for a 2d interacting, self-propelled particle system, Physica D 232 (2007), 33–47. [2] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and N. R. Franks, Collective memory and spatial sorting in animal groups, J. Theor. Biol. 218 (2002), 1–11. [3] A. Czirók, M. Vicsek, and T. Vicsek, Collective motion of organisms in three dimensions, Physica A 264 (1999), 299–304. [4] A. Czirók and T. Vicsek, Collective behaviour of interacting self- propelled particles, Physica A 281 (2000), 17–29. [5] M. R. D’Orsogna, Y. L. Chuang, A. L. Bertozzi, and L. S. Chayes, Self- propelled particles with soft-core interactions: Patterns, stability, and collapse, Physical Review Letters 96 (2006), 104302. [6] A. Huth and C. Wissel, The simulation of movement of fish schools, J. Theor. Biol. 156 (1992), 365–385. [7] , The simulation of fish schools in comparison with experimental data, Ecological Modelling 75/76 (1994), 135–145. [8] H. Levine, W.J. Rappel, and I. Cohen, Self-organization in systems of self-propelled particles, Physical Review E 63 (2001), 017101. [9] Y.-X. Li, Clustering in neural networks with heterogeneous and asym- metrical coupling strengths, Phys. D 80(3) (2003), 210–234. [10] R. Lukeman, Y.-X. Li, and L. Edelstein-Keshet, A mathematical model for milling formations in biological aggregates, Bull. Math. Biol. 71(2) (2009), 352–382. 59 Chapter 2. Bibliography [11] A. Mogilner, L. Edelstein-Keshet, L. Bent, and A. Spiros, Mutual in- teractions, potentials, and individual distance in a social aggregation, J Math Biol 47 (2003), 353–389. [12] H.-S. Niwa, Self-organizing dynamic model of fish schooling, J theor Biol 171 (1994), 123–136. [13] , Newtonian dynamical approach to fish schooling, J theor Biol 181 (1996), 47–63. [14] ,Migration of fish schools in heterothermal environments, J theor Biol 193 (1998), 215–231. [15] A. Okubo, Diffusion and ecological problems: Mathematical models, Springer Verlag, New York, 1980. [16] A. Okubo, D. Grunbaum, and L. Edelstein-Keshet, The dynamics of an- imal grouping, Diffusion and Ecological Problems: Modern Perspectives (A Okubo and S Levin, eds.), Springer, N.Y., 2001. [17] S. Sakai, A model for group structure and its behavior, Biophysics Japan 13 (1973), 82–90. [18] J. R. Silvester, Determinants of block matrices, Maths Gazette 84 (2000), 460–467. [19] A. R. E. Sinclair, The african buffalo : a study of resource limitation of populations, University of Chicago Press, Chicago, 1977. [20] R. Suzuki and S. Sakai, Movement of a group of animals, Biophysics Japan 13 (1973), 281–282. [21] J. Toner and Y. Tu, Long-range order in a two-dimensional xy model: how birds fly together, Phys. Rev. Lett. 75 (23) (1995), 4326–4329. [22] C. Topaz and A. Bertozzi, Swarming patterns in a two-dimensional kine- matic model for biological groups, SIAM J. Appl. Math. 65 (1) (2004), 152–174. [23] C. Topaz, A. Bertozzi, and M. Lewis, A nonlocal continuum model for biological aggregation, Bull. Math. Bio. 68 (7) (2006), 1601–1623. 60 Chapter 2. Bibliography [24] B. P. Uvarov, Locusts and grasshoppers, Imperial Bureau of Entomology, London, 1928. [25] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Phys. Rev. Lett. 75(6) (1995), 1226–1229. [26] D. Weihs, Energetic advantages of burst swimming of fish, J. theor. Biol. 48 (1974), 215–229. 61 Chapter 3 A Conceptual Model for Milling Formations in Biological Aggregates 2 3.1 Introduction In biological aggregations, there are many examples of group-level patterns that emerge from interactions among individuals. Such phenomena have motivated a number of modeling approaches aimed at connecting the in- teractions of individuals to the emergent collective behaviours. Lagrangian models, based on tracking the positions and velocities of individuals contrast with the Eulerian models, formulated in terms of reaction-transport equa- tions. Examples of the Lagrangian approach include [4, 8, 12, 14, 18, 19, 20, 28, 31, 33]. In the Lagrangian framework, each individual is described by several ordinary differential equations (for components of velocity and po- sition in 1-3 spatial dimensions). This means that the complexity and size of the model increases proportionately to the size of the group, rendering such models notoriously difficult to understand analytically. Due to this in- herent challenge, many previous Lagrangian models are largely dominated by simulation studies. Interesting patterns are found and classified, but an open question remains as to which aspects of a given model leads to which emergent feature of the group. Even when exhaustive search of parameter space is undertaken, linking classes of interaction functions to the types of patterns they form is a major challenge. Motivated by the desire to obtain analytic results that characterize the link between individual interactions and the patterns formed by groups, we here explore a minimal model that has the advantage of analytic tractability. 2A version of this chapter has been published as ‘R. Lukeman, Y.-X. Li, and L. Edelstein-Keshet, 2009, A conceptual model for milling in biological aggregates, Bull. Math. Bio., 71(2), p 352-382 62 3.1. Introduction The goals of our paper are to further investigate a framework that is simple enough that we can unequivocally assess which properties of interaction func- tions, parameters, and rules lead to existence and stability of group patterns. Based on the challenges described above, such a task can only be done at this point on a restricted set of simple cases where we can exploit regularity of the pattern, symmetry properties, and limited coupling to achieve analytical power. Fortunately, simulations can be used to complement this analysis, and to expand insights to cases that cannot be solved in closed form. Emergent patterns in nature can help to inform and validate models for group behaviour. One particular pattern that emerges in many of these models, and occasionally in nature, is the milling formation (also known as a ‘vortex formation’), in which individuals move in roughly concentric trajec- tories. In some previous models, milling mechanisms have been engineered in specific confined domains (due to interactions with boundaries), or by in- cluding “rotor chemotaxis” terms, representing a tendency to maintain some angle to the gradient of an attractive force-field. (We compare some previ- ous examples in our Discussion.) Since many competing models can produce similar patterns, the emergence of a mill (or for that matter, any other pat- tern) can not be used to claim “correctness” of a given model. However, it is of interest to assess which attributes of a given model produce the resulting pattern. In this paper we explore what features of our minimal model would be consistent with existence and stability of mill formations. In order to gain insight into the existence and stability properties of milling formations, we take the following approach: we consider a spatially simple solution that still captures the essential behaviour of milling, given by equally-spaced particles moving around a circle of fixed radius, at fixed angular velocity. These particles interact only via distance-dependent forces of attraction and repulsion, neglecting velocity-dependent forces common in similarly motivated models. We exploit the regularity of this idealized milling formation to address the following questions analytically: (i) Under what conditions on the school forces do such idealized formations occur (ex- istence)? (ii) Given their existence, what further conditions guarantee that these patterns are not destroyed by random perturbations (stability)? Under this approach, we do not seek to realistically model the observed behaviours in nature (which are undoubtedly more complex), but instead seek to under- stand the connection between individual interaction and group stability via an abstraction of the milling phenomenon. Focusing on a simpler pattern while using a simple model leads to clear-cut results that can be derived 63 3.1. Introduction analytically. This lends insight into the more complicated group structures found in nature. Our approach is somewhat similar to [17] who consider similar formations from a multi-agent control system perspective, and derive stable steering controls for particles moving with fixed speed. In a companion paper [15], a model for school formation was introduced based on the Lagrangian viewpoint, that is, following individual particles, rather than densities of organisms. Equations of motion were formulated based on the Newtonian approach, i.e., describing changes in the velocities and positions of the particles under forces of propulsion and interaction. This self-propelled particle model was used to study perfect schools, which are con- figurations of particles in which the spacing is constant and identical, and where individuals in the group move at constant speed. One such perfect school in two dimensions, the soldier formation, was studied in [15] using linear stability analysis. Analytic stability bounds were found in terms of the slope of the interaction force as a function of interindividual distance. Numerical simulations were used to validate analytic conditions obtained, and investigate more complicated perfect schools. In this paper, we inves- tigate milling formations using the same theoretical framework as described above. Compared with linear aggregates [15], the dynamics of the milling for- mation are significantly more complicated. The dependence of stability on the slope of the interaction function evaluated at interindividual distance is given in terms of a complex-coefficient quartic polynomial. We investigate solutions to this equation numerically, classifying distinct behaviours in each stability regime. The structure of the interaction matrix that determines how individuals interact is crucial to the stability analysis; consequently, block circulant matrices and their properties enter heavily in our analysis. In Section 2, we introduce the basic assumptions of the model and the differential equation system that describes the motion of the particles. We develop existence conditions for the mill formation in Section 3, and stabil- ity conditions in Section 4. In Section 5, existence and stability conditions are investigated numerically, through numerical solution to the eigenvalue equation, and through simulations. A discussion follows in Section 6. 64 3.1. Introduction Figure 3.1: An example of milling in Atlantic bluefin tuna, courtesy of Dr. M. Lutcavage, Large Pelagics Research Center, UNH. Milling Formations in Nature The milling formation is most often observed in schooling fish. This phe- nomenon was described by A. E. Parr in 1927 [24]. Numerous species of fish, including jack, barracuda, and tuna have been observed to mill in nature [4]. Photographs of jack milling in [25, 26] are presented as examples of emer- gent properties of fish schools. Milling tuna have been observed in the Gulf of Maine (see Figure 3.1), while barracuda mills abound in underwater pho- tography. Basking sharks have also been observed in mills during courting behaviour [10, 35]. Another example of milling is in army ants. In a thorough treatment, Schneirla [29] describes a rare field observation of Eciton praedator exhibiting milling behavior following a rain-induced separation of a group of ants from the main columnar raiding swarm. Although the sensory mechanism and environments differ from fishes to ants, in both cases individuals turn towards a dominant centripetal stimulation [29]. 65 3.2. The Model of Interacting Self-Propelled Particles The above instances of milling in nature are self-organized in the sense that the global pattern emerges solely from interactions among individuals, using only local information [2]. Importantly, the occurrence of such mills is relatively rare in many instances, as they emerge (for example in ants) in unusual or extreme cases. (Schneirla [29] reports that ants in such mills follow one another to death by exhaustion.) Fortunately, given a wide variety of initial conditions, army ants form trunk trail foraging patterns, rather than mills; but when a closed-loop topology occurs by chance, their mill is a stable and long-lasting movement pattern. These observations raise an interesting question concerning the mecha- nisms that ensure the occurrence of self-organized milling patterns. The re- alistic mechanisms remain elusive since we know little about the actual rules of how animals interact. However, a number of models have been proposed in which milling patterns that closely resemble the observed mills have been shown to occur [8, 14]. Although experimental verifications are required to test whether these models are realistic, they provide insights into a number of mechanisms that are capable of producing milling. 3.2 The Model of Interacting Self-Propelled Particles The following basic assumptions apply to the idealized schools that we study in this paper. (A1) All particles are identical and obey the same set of rules. (A2) Particles are polarized, that is, each has a front and a back, and a particle senses other particles only in front of it. (This is a simplification of a classical assumption by [12].) (A3) The “effective force” experienced by one particle from a given neighbour depends only on distance. (A4) A given particle interacts with its nearest neighbour(s). For analytic tractability we consider interactions with at most one neighbor. (A5) The distance-dependent force between a pair of neighbours is a vector parallel to the line segment connecting the pair. 66 3.2. The Model of Interacting Self-Propelled Particles Consider a group of n self-propelled particles (i = 1, 2, . . . , n), with posi- tion ~xi and velocity ~vi. For simplicity, we assume that the body alignment of the ith particle is identical to its direction of motion, i.e. the direction of ~vi. The movement of the i th particle is governed by the classical Newton’s law of motion with equations ~̇xi = ~vi, (3.1) ~̇vi = ~ai + ~fi − γ~vi, (3.2) where ẋ ≡ dx/dt denotes a derivative with respect to time and the mass of each individual is scaled so that mi = 1 for all particles. ~ai is an autonomous self-propulsion force generated by the ith particle that may depend on envi- ronmental influences and on the location of the particle in the school. ~fi is the interaction or schooling force that moves the ith particle relative to its neighbour. The term γ~vi is a damping force with a constant drag coefficient (γ > 0) which assures that the velocity is bounded. Equation (3.2) implies that if a particle stops propelling, its velocity decays to zero at the rate γ. Particles are modeled in 2-dimensional space, so for a group of n particles, there are 4n first-order equations. In the following analysis, ~ai is a constant vector, identical for all particles. Based on the model assumptions (A3-A5), the schooling force is given by ~fi = g(|~xj − ~xi|)(~xj − ~xi)|~xj − ~xi| , (3.3) where the index j refers to the nearest neighbour detected by particle i and g(x) is a function of the absolute distance between particles (i.e., defined for x > 0) that gives the magnitude of the interaction force. To produce a nonzero spacing, g(x) should typically be positive for large x and negative for small x, indicative of short-range repulsion and long-range attraction. This feature is universal in almost all models of aggregate formation. (See review of attraction-repulsion models in [18].) 3.2.1 Relationship to Previous Models The model presented by Equations (3.1)–(3.2) is well-known. It is discussed in [22, 23] and attributed to schooling studied by [28] and [31]. 67 3.2. The Model of Interacting Self-Propelled Particles We briefly compare to other models for self-propelled particles. For com- parison, we write (3.1)–(3.2) as ~̇xi = ~vi, (3.4) ~̇vi = ~faut + ~fint, (3.5) where ~faut is the force that is generated autonomously by each individual and ~fint that is generated due to the interaction with others. In [3, 8, 19, 20, 21], ~faut = (α − β|~vi|2)~vi was studied. This model was originally obtained by minimizing the specific energy cost of a swimming fish [34]. In [14], ~faut = α ~vi |~vi| −β~vi (in the case of no velocity averaging). In both cases, particles tend to a preferred characteristic speed in the absence of interaction (in the first case, v = √ α/β, while in the second, v = α/β). In this paper, faut = −β~vi, and so in the absence of interaction, the speed of particles decays to zero. For the force of interaction, in [3, 8, 14], fint is chosen as an all-to-all coupled distance-dependent Morse-type interaction. In our analysis, we leave the form of fint general, but choose a nearest-neighbour coupling over all-to-all coupling for analytical treatment. We note that the idealized mill formations considered here also occur when using the model equations of [3, 8, 14] under our nearest-neighbour interaction regime. 3.2.2 The Milling Formation In a perfect mill, n particles move continuously around a circle of radius r0 with constant angular velocity ω0. The particles are equally spaced at distance d, labeled sequentially, with particle i sensing particle i + 1, and particle n sensing particle 1. The sector angle defined by adjacent particles in the mill is θ = 2π/n, where n is the number of particles (see Figure 3.2 for labeling conventions). To simplify the presentation of our analysis, we set ~ai = 0 for all i; that is, particles experience no self-propelling force in the absence of interactions. Later, we discuss the case for nonzero ~ai. These assumptions lead to d~xi dt = ~vi, (3.6) d~vi dt = ~fi(~xi+1 − ~xi)− γ~vi, (3.7) 68 3.3. Existence Conditions Figure 3.2: A schematic diagram of the milling formation. Black arrows indi- cate direction of motion (tangential to the circle), while grey arrows indicate direction of schooling force. for i = 1, . . . , n− 1, and d~xn dt = ~vn, (3.8) d~vn dt = ~fn(~x1 − ~xn)− γ~vn. (3.9) Equations (3.8)–(3.9) link the nth and 1st particles into a (periodic) ring formation. 3.3 Existence Conditions The milling solution illustrated in Figure 3.2 represents a steady-state equidis- tant distribution of n particles on a ring of radius r0 when observed in the frame that rotates at constant angular velocity ω0. Solving (3.6)–(3.9) for such a steady state in the rotating frame yields r0, ω0, and interindividual distance d that uniquely determines such a milling solution. Here, we use a simple force-balance argument to obtain the same values. First, decompose the interaction force ~f into a component tangential to the circle, ~ft, and a centripetal component, denoted ~fc (see Figure 3.3). In a steady-state mill 69 3.3. Existence Conditions Figure 3.3: A schematic diagram showing the decomposition of the interac- tion force into tangential and centripetal components. solution, there is no linear acceleration, so tangential forces are in balance. Thus, for each particle, ~ft − γ~v = 0. (3.10) The centripetal force that maintains the motion of a particle of mass m = 1 around a circle of radius r0 at angular speed ω0 is ~fc = ~ac = ω 2 0r0~uc, where ~uc is a radially-directed unit vector. Because |~v| = r0ω0, we have ~fc = |~v|2 r0 ~uc. (3.11) Denote by φ and ψ the angles subtended to tangent and radius, respectively, as in Figure 3.4. Then θ = 2π/n, φ = θ/2 = π/n, and ψ = π/2−π/n. Using the law of cosines and simple trigonometry, we find that d = 2r0 sin ( π n ) , (3.12) and express force magnitudes as |ft| = g(d) cos ( π n ) , |fc| = g(d) sin ( π n ) . 70 3.3. Existence Conditions Figure 3.4: A schematic diagram of particles in the mill formation showing angles introduced in the text, and relative position and velocity vectors. Using (3.10) and (3.11), g(d) cos ( π n ) = γ|~v|, (3.13) g(d) sin ( π n ) = |~v|2 r0 , (3.14) and solving for |~v| and r0 gives |~v| = g(d) cos ( π n ) γ , r0 = g(d) cos2 ( π n ) γ2 sin ( π n ) . (3.15) Then angular frequency, ω0 is related to speed and mill radius, r0 by ω0 = |~v|/r0 = γ tan ( π n ) . (3.16) For given g(x), n and γ, (3.15) and (3.16) give the equilibrium radius and angular velocity of the milling formation. For an n-particle system, these quantities completely characterize the milling solution. Combining (3.12) and (3.15), we obtain the existence condition g(d) = sd, (3.17) 71 3.4. Stability Analysis [a] [b] Figure 3.5: Equation (3.17) represented graphically for γ = 0.5, n = 5, and g(x) = A exp(−x/a) − B exp(−x/b). Horizontal axis: inter-individual distance d, vertical axis: distance-dependent force magnitude. A = 1.5, a = 10, B = 3. In panel [a], b = 1.5, and milling is possible whereas in panel [b], b = 2 and no milling solution exists. where s = γ2/2 cos2 ( π n ) . A steady-state mill formation occurs if and only if there exists a value of d satisfying (3.17) for a given function g(x). Inter- sections of g(d) and the straight line sd are mill formations. In Figure 3.5, intersections indicate potential steady-state values of d for which the mill formation can exist in [a], while in [b] absence of intersections indicate that a steady-state mill formation is not possible. Significantly, the slope s of the right-hand side of (3.17) affects whether or not an intersection exists. For typical g(x), a shallower slope increases the likelihood of an intersection, so decreasing the damping coefficient γ or increasing the number of individuals n increases the likelihood of existence of the milling state. The dependence on n implies the possibility of a mill formation being destroyed when one or more particles leaves the group. Figure 3.6 shows a plot of Equation (3.17) for various values of n. In this example, for fewer than 5 particles, a mill solution will not exist. 3.4 Stability Analysis In this section, we investigate the local stability of the milling solution via a linear stability analysis. Our goal is to linearize (3.6)–(3.9) at the milling so- lution and study the eigenvalues that determine the stability of the solution. 72 3.4. Stability Analysis Figure 3.6: A plot of the left-hand and right-hand sides of (3.17) for a range of n values, from n = 4 to n = 12. Note that no intersection occurs for n = 4, but two occur for n = 5. Here we exploit the simplifications and assumptions of the model, without which stability matrices become unwieldy. The cyclic nature of particle in- teractions is reflected by a block-circulant structure in the stability matrix, which has the form A1 A2 · · · An An A1 · · · An−1 ... A2 · · · An A1 , where the Ai arem×m block matrices. Because an explicit formulation of the determinant is known for an mn×mn block-circulant matrix, we can exploit this structure in calculating eigenvalues. However, here we encounter a chal- lenge: because ~xi+1−~xi terms appear in the stability matrix, and this vector has a different direction for each particle, successive rows of blocks are not identical modulo a shift. To overcome this issue, we transform coordinates. Let ~̂xi = R(φi) (~xi+1 − ~xi) , (3.18) ~̂vi = R(φi) (~vi+1 − ~vi) , (3.19) 73 3.4. Stability Analysis where R(φi) = [ cos(φi) sin(φi) − sin(φi) cos(φi) ] , (3.20) is a rotation matrix with angle φi = 2πi n + ω0t+ π n + π 2 . (3.21) This transformation is best understood as a composition of a number of transformations. First, position and velocity differences are taken to simplify the system of equations. Second, these differences are rotated by the angle ω0t, transforming the system to a rotating frame. In this rotating frame, position and velocities appear to be fixed at the milling solution. Third, an index-dependent rotation of 2πi/n is applied to each particle so that steady- state quantities are independent of index (required for the circulant structure in the stability matrix). Last, a constant rotation of π/n + π/2 is applied so that the steady-state position and velocity vectors lie on the coordinate axes. We note that applying these three rotations does not affect eigenvalues because the linearized coefficient matrix and its analogue without rotation are similar matrices. The transformed system of equations governing rotated position and ve- locity differences is d dt ~̂xi = ω0k× ~̂xi + ~̂vi, (3.22) d dt ~̂vi = ω0k× ~̂vi +R (−2π n ) ~̂xi+1 |~̂xi+1| g ( |~̂xi+1| ) − ~̂xi|~̂xi| g ( |~̂xi| ) − γ~̂vi, (3.23) and steady-state values are given by ~̂x s i = ~̂x s = [ d 0 ] , ~̂v s i = ~̂v s = [ 0 ω0d ] , (3.24) for i = 1, 2, . . . , n with n + 1 identified with 1 (see Appendix A for details). The extra rotation R ( −2π n ) in (3.23) is required to write this equation in terms of the transformed coordinates. Note that both ~̂x s i and ~̂v s i are indepen- dent of the index i; i.e., in this coordinate frame, all particles have the same 74 3.4. Stability Analysis steady-state coordinates ~̂x s and ~̂v s . Next we introduce perturbations to the steady state milling solution in the transformed coordinate system; i.e., ~̂xi = ~̂x s + ~δi(t), (3.25) ~̂vi = ~̂v s + ~ξi(t). (3.26) We substitute (3.25)–(3.26) into (3.22)–(3.23) and linearize interaction terms in (3.23). Then, taking together equations for all i = 1, . . . , n, we can write the full perturbed linear system in matrix notation. To facilitate this, we introduce a number of 2× 2 matrices. We let Ω = [ 0 ω0 −ω0 0 ] , D = ∂ ~f ∂~̂xi |(d,0) = [ g′(d) 0 0 g(d)/d ] , RD = R (−2π n ) [ g′(d) 0 0 g(d)/d ] . (3.27) The full 4n× 4n linear system is then d dt ~δ1 ~δ2 ... ~δn ~ξ1 ~ξ2 ... ~ξn = Ω | I Ω | I . . . | . . . Ω | I −D RD | Ω− γI −D RD | . . . . . . | . . . RD −D | Ω− γI ~δ1 ~δ2 ... ~δn ~ξ1 ~ξ2 ... ~ξn , (3.28) where I is the 2×2 identity matrix. The coefficient matrix is composed of four n× n block matrices, three of which are block-diagonal, and one of which is block-circulant. Each component is itself a 2×2 matrix. If we consider cases with interactions with two or more nearest neighbors, the structure of such matrices becomes intractable using our techniques. For the full derivation of (3.28), see Appendix B. We note that ω0 is given explicitly in terms of n and γ in (3.16), and we can rearrange (3.17) to obtain g(d) d = s, (3.29) 75 3.4. Stability Analysis i.e., expressing g(d)/d explicitly in terms of n and γ. Thus, we reduce the parameters in (3.28) to the drag coefficient γ, the number of particles n, and the slope of the interaction function at steady state, g′(d). 3.4.1 Eigenvalues We denote the coefficient matrix in (3.28) as C, and solve the eigenvalues of C via det (C− λI4n×4n) = 0, where I4n×4n is the 4n× 4n identity matrix. Now we use the following result from [30]: given the block matrix [ A1,1 A1,2 A2,1 A2,2 ] , whereA1,1,A1,2,A2,1, andA2,2 arem×mmatrices, ifA1,1 andA1,2 commute, then det [ A1,1 A1,2 A2,1 A2,2 ] = det [A2,2A1,1 −A2,1A1,2] . (3.30) We partition C − λI4n×4n into four n × n block matrices and apply the above theorem, noting that the corresponding upper left and upper right submatrices commute, to obtain the reduced determinant equation: det (C− λI4n×4n) = det Λ +D −RD 0 · · · 0 0 Λ +D −RD 0 · · · 0 . . . −RD 0 · · · 0 Λ +D = 0, (3.31) where 0 is the 2× 2 zero matrix, and where we have defined Λ ≡ (Ω− λI)(Ω− (γ + λ)I) = [ −λ ω0 −ω0 −λ ] [ −λ− γ ω0 −ω0 −λ− γ ] = [ λ(λ+ γ)− ω20 −γω0 − 2λω0 γω0 + 2λω0 λ(λ+ γ)− ω20 ] . (3.32) Next, we use a diagonalization method for block-circulant matrices from [7]. The diagonalization matrix is given by Fourier matrices; for details see [17] 76 3.5. Numerical Investigation or [7]. The resulting n 2× 2 diagonal blocks are given by D1 = Λ +D+RD, D2 = Λ +D+̟RD, D3 = Λ +D+̟ 2RD, ... Dn = Λ +D+̟ n−1RD, where ̟ = exp (2πj/n) , where j = √−1. Note that the powers of ̟ give the n roots of unity. Because the determinant of a block-diagonal matrix is the product of the determinants of the individual blocks, solutions to (3.31) are solutions to detD1 detD2 · · ·detDn = 0, i.e., the collection of solutions to detDi = 0, i = 1, . . . n. (3.33) Expanding the (i + 1)th equation leads to a 4th-order complex-coefficient polynomial p(λ) = 0 (3.34) to be solved for eigenvalues λ (see Appendix C). Explicit solutions to this equation are not easily obtained. However, given parameter values for γ, n, and g′(d), roots can be found numerically. We have thus far avoided the n = 3 case. Because all three particles are equidistant in a mill formation, a particle can switch nearest neighbours under a small perturbation, breaking the connection topology. Thus the three-particle mill formation is inherently unstable to perturbations. 3.5 Numerical Investigation We are interested in restrictions on the interaction function that guarantee stable mill solutions, i.e., such that linear stability theory points to eigen- values with negative real parts. Note that based on (3.33), the eigenvalue problem involves the matrices Λ, R and DR. For fixed n and γ, these ma- trices depend only on g(d)/d and g′(d), and because at the milling steady state, g(d)/d = s, eigenvalues depend only on g′(d). We thus numerically 77 3.5. Numerical Investigation investigate stability with respect to g′(d) by treating this quantity as a bifur- cation parameter, with the remaining parameters, n and γ, fixed. In Figures 3.7 - 3.9, we plot the real part of λ, Re(λ), for each of the four solutions to Equation (3.34), for each of i = 0, . . . , n − 1, for a total of 4n eigenvalues. We note that for all values of n that were investigated, three eigenvalues with Re(λ) = 0 existed: a zero eigenvalue, associated with the rotational symmetry, and a pair of purely imaginary conjugates, associated with the translational invariance of the center of the mill within a rotating frame. For Figure 3.7, n = 4. There is a range of values of g′(d) in [-0.26,0.25] for which the real part of all eigenvalues are nonpositive, consistent with stability. In Figures (3.8) and (3.9), n = 5 and n = 6, respectively. In contrast to the previous case, the stable range of g′(d) shrinks successively, from [−0.11, 0.19] for n = 5, to [−0.08, 0.17] for n = 6. These results suggest that the mill solution is unstable if the slope of the interaction function at d is either too large or too small (i.e., too steep). Stable mills exist when g′(d) < 0, evidently due to the circular geometry, as linear formations studied in [15] required that g′(d) > 0. In each of these three cases, the upper limit of the stability region coin- cides with the value of g(d)/d at steady state for each n, that is, g′(d) < s, (3.35) for stable mills. (Specifically, for γ = 0.5, and n = 4, 5, 6, we get g(d)/d = 0.25, 0.191, and 0.166 respectively). In each of these cases, the stability boundary is determined by a branch of eigenvalues that crosses from nega- tive to positive real parts at the boundary value of g′(d). It can be shown analytically that a branch of Re(λ) crosses over the g′(d) axis at g(d)/d by examining coefficients of p(λ) when g′(d) = g(d)/d. However, it remains to be proven that this crossover value always represents the upper-bound of the stability region. As n increases, we have found numerically that the lower bound on the stability region approaches 0. 3.5.1 Numerical Simulations A self-organized pattern in a dynamical system refers to a structurally stable solution that occurs as a result of interactions between individuals in a group [8, 14, 15] and does not occur as a result of forcing constraints asserted by other groups (e.g. predators) and/or the environment [9, 11]. Under the same 78 3.5. Numerical Investigation Figure 3.7: Numerical solution to (3.34) over a range of values of g′(d) for four particles (n = 4) and γ = 0.5. The real parts of eigenvalues are plotted for i = 0, . . . , 3. The asterisks indicate branches of Re(λ) with multiplicity equal to 2. Note that the region over which all eigenvalues have non-positive real parts is approximately [-0.26,0.25]. 79 3.5. Numerical Investigation Figure 3.8: As in Figure 3.7, but for five particles (n = 5). The real part of eigenvalues are plotted for i = 0, . . . , 4. The region over which all eigenvalues have non-positive real parts is now approximately [-0.11,0.19]. 80 3.5. Numerical Investigation Figure 3.9: As in Figure 3.7, but for six particles (n = 6). The real part of eigenvalues are plotted for i = 0, . . . , 5. The region over which all eigenvalues have non-positive real parts is now approximately [-0.08,0.17]. 81 3.5. Numerical Investigation interaction rules and model parameters, a number of stable self-organized patterns can occur and coexist. The emergence of one particular pattern is closely related to the initial state or configuration of the group. Large disturbances can also induce a switch from one pattern to another coexisting one. Mills that we studied in this model are obtained in simulations only when rather specific initial conditions were used, whereas most random initial conditions lead to other patterns. Therefore, they are self-organized patterns that occur under stringent initial configurations. This is consistent with the fact that milling is not the most commonly observed pattern in animal aggregates (as is the case in the example of milling army ants). The particular restriction on initial conditions required for the emergence of a milling solution in (3.1)–(3.2) is that the nearest-neighbour connection topology of the group is a simple closed curve. To ensure this geometry, we chose initial positions to be equidistant individuals around a circle of fixed radius r, and then perturbed these (x, y) coordinates by values chosen randomly in [−r/4, r/4]. Initial velocities are tangential to the circle of ra- dius r, with magnitude ω0r. We point out that once this closed connection topology (with arbitrary geometry) is established in our simulations, the cir- cular geometry that emerges is due only to local interactions among particles, without reference to any global information, and thus is self-organized [2]. The model equations (3.1)–(3.2) are evolved using ODE45 in MATLAB. Figure 3.10 compares simulations of n = 6 particles in the cases for which Equation (3.17) has no solution, versus at least one solutions. In Figure 3.10.a1, the left-hand and right-hand side of (3.17) are plotted, and inter- sections indicate individual distances d such that a mill formation can exist. However, our stability analysis indicates that the first intersection corre- sponds to an unstable solution because g′(d) > s. At the second intersection, the slope g′(d) is within the range of stability, and as confirmed by Figure 3.10.a2, particles form a mill with spacing d ≈ 2, and r ≈ 1.9. For the interaction function in Figure 3.10.b, no solution exists for (3.17). Corresponding simulations (Figure 3.10.b.ii) show particles initially forming a mill-like pattern; however, the radius of the group decreases until the indi- vidual distance d is such that g(d) = 0. At this stage, no further interactions occur, and the particles stop moving. We do not technically consider such stationary circular patterns to be “mills”. In the case that there is no zero of the interaction function g, the group collapses to a point. We now investigate the behaviour of the system at the boundary of stabil- ity, and inside the region of instability. Earlier, we showed numerically that 82 3.5. Numerical Investigation Figure 3.10: A plot of the condition for mill solutions to exist [a1] or not [b1] given by Equation (3.17), and corresponding simulated particle tracks [a2] and [b2], with γ = 0.5. Note that in the case of an intersection, a stable mill forms at d ≈ 2 (whereas d ≈ 0.5 is unstable). However, when no intersection exists, the mill stops rotating and stationary particles are then arranged with d such that g(d) = 0 (final direction is outward due to small repulsive forces felt just before the particles stop). Arrows indicate direction facing at the end of the simulation. The interaction function used is as in Figure 3.5 with A = 0.5, a = 5, B = 1, and [a] b = 0.5, [b], b = 0.8. 83 3.5. Numerical Investigation the system is stable over a region of g′(d) which includes g′(d) = 0. Thus, choosing an interaction function whose slope at d is either too positive, or too negative, will result in instability in the mill solution. Just how this instability manifests itself can be seen in simulations. Figure 3.11.a1 shows an interaction function with two intersections indicating possible individual distances. However, at the first intersection, g′(d) is too positive and Equa- tion (3.35) is not satisfied (as indicated by the straight line), while at the second, g′(d) ≈ −0.10 is too negative. Both values are outside the stability region [-0.07,0.17] noted in Figure 3.9 for n = 6. The corresponding sim- ulation in Figure 3.11.a2 shows particles initially forming a circular group, though oscillations in individual distance grow until the connection topology is altered and the mill is broken. If g(d) is linear with slope greater than s, we obtain a single intersection as in Figure 3.11.b1. In this case, simula- tions (Figure 3.11.b2) show a circular group forming, though the radius of the group grows exponentially in time, again confirming instability. Thus there are a number of ways in which the stability of the mill formation is lost, including radial increase in time, evolution to a stationary group, and breaking of the connection topology. The system has yet a different behaviour with g′(d) chosen to lie on the stability boundary. In this case, oscillations in individual distance d develop (as in the connection-breaking instability), but the oscillation magnitude ap- proaches a fixed value and the system evolves to a new steady-state solution, the irregular periodic mill. Despite the periodic oscillation of distance be- tween particles, the sum of the n distances between neighbouring particles is fixed when the irregular periodic mill is established. A plot of d in time showing this periodic behaviour is shown for one particle in an n = 5 system in Figure 3.12.a, while snapshots of the particles in time are shown in Figure 3.12.b. 3.5.2 Moving Mill Formations In our treatment of mill formations up to this point, we have assumed that ~a = 0; i.e., individuals are compelled to move only based on interactions with neighbours, and have no intrinsic desire to move in any direction. In this section, we explore the effect of assigning a non-zero autonomous self- propulsion to each individual. In Figure 3.13, six particles initialized as usual, each with ~a = (0.1, 0.1)T , form a mill that moves in the direction of ~a. Thus, autonomous self-propulsion 84 3.5. Numerical Investigation Figure 3.11: [a1] An interaction function with g′(d) lying outside the stability region for both intersections (implying no stable mill formation). Simulations using this function are shown in [a2]. Note that particles initially form a mill- like solution which is destroyed as time evolves. [b1] An interaction function with only one intersection, whose slope is too large (g′(x) > s for all x). Simulations in [b2] show that the radius of the mill increases in time. 85 3.5. Numerical Investigation Figure 3.12: The irregular periodic mill formation. [a] Individual distance d versus time for one particle in a system of 5 particles. For the 4 times denoted in [a], snapshots of particles are shown in [b], showing the variation in d between particles 1 and 2. The interaction function used was g(x) = 1−0.11x, so that g′(d) = −0.11 lies on the stability boundary. 86 3.6. Discussion in each individual acts to propel the entire mill along the direction of this in- dividual propulsion. However, because the velocity of particles is now jointly composed of autonomous self-propulsion and interaction forces, the mill for- mation can be broken if ~a dominates the motion of individuals. Mill breaking occurs if some particle i no longer senses particle i+ 1; i.e., ~vsi · ( ~xsi+1 − ~xsi ) ≤ 0. (3.36) Substituting known steady-state quantities in (3.36) and using some calculus, the mill-breaking condition reduces to |~a| ≥ g(d) cos2 ( π n ) , (3.37) (see Appendix D for derivation). In Figure 3.14.a, we set ~a = (0.223, 0.446)T , so that |~a| = 0.498 < g(d) cos2 ( π n ) = 0.5, and a stable moving mill is formed. However, in Figure 3.14.b, a small increment in ~a to (0.224, 0.448)T gives |~a| = 0.501 > 0.5, resulting in a fundamental shift in behaviour: the mill formation breaks, and particles form a polarized group with a shift to d ≈ 2 (see the inset of Figure 3.14.b). This new solution is consistent with those studied in [15]. In this example, the increase in the parameter ~a initiated a change in the connection topology of the particles, which in turn instigated an evolution to a different solution. Indeed, for any given set of parameters, there are a number of types of solutions that can occur; it is the arrangement of particles and headings that determines the solution to which the system converges. 3.6 Discussion Milling formations exist in a number of biological aggregates, and provide fascinating examples of self-organized spatial patterns evolving from local interactions among individuals. Such patterns found occasionally in nature are irregular owing to heterogeneity in the population and environment, and likely more complex interactions. To reveal the cause-and-effect relationship between interaction properties and the patterns that emerge, we focused on idealized mill formations of particles rotating in a circular path of fixed radius. Although much simpler than the natural counterpart, this pattern permits mathematical analysis, which enabled us to make clear predictions 87 3.6. Discussion Figure 3.13: A plot of trajectories in time for six particles with ~a = (0.1, 0.1)T . Note that the entire mill moves in the direction of ~a. of how model parameters influence the existence and stability of these pat- terns. Interactions in this model are distance-dependent, and so, likely, do not completely characterize the interactions in natural mill formations. Yet, a detailed understanding of distance-dependent interactions should provide building blocks for understanding models incorporating more complex inter- actions. Using a Lagrangian model based on Newton’s equations of motion, we described mill formations in terms of group radius and angular velocity, and derived the following necessary condition for existence of a mill formation with spacing distance d, where γ is the drag coefficient, and n is the number of individuals: g(d) = sd, where s = γ2/2 cos2 (π/n). Increasing n, and/or decreasing γ increases the likelihood that the existence condition is satisfied. A linear stability analysis on the milling formation resulted in block- circulant matrices that were block-diagonalized by Fourier matrices. This analysis could be done only in a limited, and particularly simple set of cases. 88 3.6. Discussion Figure 3.14: A plot of trajectories in time for five particles with [a] ~a = (0.223, 0.446)T , and [b] ~a = (0.224, 0.448)T . Note that for a small parameter change, the system behaviour is fundamentally different. In these simula- tions, g(x) is as in Figure 3.5, with A = 1.5, a = 10, B = 3, and b = 1.6. The corresponding existence condition is shown in the inset to [b]. 89 3.6. Discussion Notably, the assumption of interactions with a single nearest neighbor was essential to limit the bandwidth of these matrices. The resulting determi- nant equation for eigenvalues yields a system of n complex-coefficient quartic polynomials, whose combined 4n solutions gives the spectrum of the stabil- ity matrix. Through numerical solution of these equations, stability bounds were found for particular cases of n and γ. The geometry of the mill formation differs from soldier formations stud- ied in [15] only through the nth particle being connected to the first particle. Yet, remarkably, this connection adds considerable complexity to the issue of stability of the system. The most interesting stability dependence was seen in the slope of the interaction function at steady-state, g′(d). A series of ranges were computed numerically for fixed γ and various values of n. We found numerically that as n increases, the region of stability decreases, asymptotic to [0, γ2/2]. Interestingly, we found that increasing γ (or decreasing n) in- creases the stability region, but decreases the likelihood of a mill formation existing. Unlike linear schools studied in [15], stable solutions can exist when g′(d) < 0. This is nonintuitive, as g′(d) < 0 means that a given particle i has a weaker attraction as it moves away from its neighbour to the front, i + 1. This would seemingly inhibit restoring the steady-state milling so- lution. However, the circulant connection of particles must compensate for this weaker attraction by inducing the particle to the back, i− 1, to be more strongly attracted to particle i, and as this stronger “pull” propagates suc- cessively back around the mill, the attraction of particle i to i+1 is increased, imparting a stabilizing effect. Numerical simulations validated our calculation of stability regions. Nu- merous routes to instability can occur, including evolution to a motionless steady-state, unbounded increase in radius, and breaking the circulant con- nection topology. Simulations of the system at the boundary of stable and unstable regions of g′(d) indicated the existence a new type of solution, the irregular periodic mill, featuring a stable limit cycle in interindividual dis- tances, resulting in a rotating, irregular, dynamic polygon of individuals with side lengths oscillating in time. Inclusion of an autonomous self-propulsion term, ~a, for each particle was shown to move the entire mill in the direction of ~a, until |~a| ≥ g(d) cos2 ( π n ) , beyond which milling formations evolved to polarized groups similar to those studied in [15]. The model used in this paper is similar to other established self-propelled 90 3.6. Discussion particle models as noted earlier. Slightly different interaction forces were used by some studies [3, 8, 14, 19, 20, 21]. In some of these particles tend to a preferred characteristic speed in the absence of interaction. As our re- sults are analytical and not just numeric, we could consider a more general interaction function fint than adopted in other papers. However, this comes at the expense of restricting attention to (single) nearest-neighbour coupling over all-to-all coupling. Furthermore, our particles sense others only in the forward direction, an assumption which can be thought of as a generalization of the ‘blind zone’ assumptions of fish models (e.g., [4, 12]). Our choice of a 180◦ blind zone could be reduced while still maintaining our mill stability results, so long as the blind zone width is sufficient to prevent a perturbed individual at the back being sensed, breaking the connection topology es- sential for maintaining the mill. Both the all-to-all coupling and our choice of coupling is likely unrealistic for most naturally occurring mills, yet both frameworks lend insight into the phenomenon. As outlined in the introduction, many other competing models can give rise to mills. Milling formations emerge in these models through a number of mechanisms, and an interesting comparison can be made between our mills and those produced by others types of models. In [1, 6], an individual-based model for collective motion of bacteria is used to generate ‘vortices’ via a rotor chemotaxis term in the equations of motion. Including such terms in individual-based ODE’s explicitly enforces rotation: particles moving tan- gentially to a unimodal attractant source distribution or central force field would exhibit rotational motion. Milling can be generated via an exter- nal environmental gradient, as is the case in the individual-based model of swarming Daphnia in [16]. Milling can also emerge in models where particles have a non-zero equilibrium speed in a closed domain; such is the case in the coupled lattice map model of [11] and the spring-dashpot model used in [9]. Meanwhile, several models have exhibited milling solutions using only interactions among individuals; this can occur in all-to-all coupled distance- and velocity-dependent interactions [14], or distance-dependent interactions only [3, 14]. In [4], interactions between individuals are restricted to localized ‘zones’ of attraction, repulsion and alignment, and yet mill formations emerge in certain parameter regimes. From a continuum perspective, in [5, 13, 32] population density models have vortex-type formations as solutions. The milling formations that emerge in the individual-based models dis- cussed above are often qualitatively similar to such formations observed in nature. In [8], conditions are derived that determine the H-stability of milling 91 3.6. Discussion solutions (that is, whether or not the group will collapse to a point as more and more individuals are added). Also, milling formations have been stable to small amounts of noise in numerous studies. However, the analytical de- termination of conditions under which such formations exist and are stable to perturbations of its members remains unanswered due to the spatial com- plexity of these formations. Passing to the continuum limit often facilitates analysis by a reduction in order of the system, but this occurs at the cost of individual properties in favor of mean field properties. In [27], a linear stability analysis was performed on a vortex-type solution of a continuum model, but was inconclusive. It is not clear if any functional benefit is associated with milling forma- tions in nature, or if these are epiphenomena that result from aggregation in special circumstances. In any case, our results on the relationship between model parameters and observed formations may prove useful in furthering the qualitative understanding of these formations. Furthermore, our results are useful for the design of artificial schools (i.e., multi-agent systems) where such rotational coordinated motion is desired, as our modeling framework and stability results provide a simple and robust method of programming such agents. Future work will include obtaining analytical stability bounds on g′(d), which have been presented numerically here. An interesting extension of this work would be to investigate the effect of a velocity-dependent interaction in addition to the distance-dependent interactions studied here. Such an ex- tension would result in a more realistic model for natural groups. As one considers more complicated formations of particles, the issue of labeling par- ticles and determining their interaction topology may present a considerable challenge. 92 Bibliography [1] E. Ben-Jacob, I. Cohen, A. Czirók, T. Vicsek, and D. Gutnick, Chemo- modulation of cellular movement, collective formation of vortices by swarming bacteria, and colonial development, Physica A 238 (1997), 181–197. [2] S. Camazine, J.-L. Deneubourg, N. R. Franks, J. Sneyd, G. Theraulaz, and E. Bonabeau, Self-organization in biological systems, Princeton Uni- versity Press, Princeton, 2001. [3] Y. L. Chuang, M. R. D’Orsogna, D. Marthaler, A. L. Bertozzi, and L. S. Chayes, State transitions and the continuum limit for a 2d interacting, self-propelled particle system, Physica D 232 (2007), 33–47. [4] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and N. R. Franks, Collective memory and spatial sorting in animal groups, J. Theor. Biol. 218 (2002), 1–11. [5] Z. Csahók and A. Czirók, Hydrodynamics of bacterial motion, Physica A 243 (2008), 304–318. [6] A. Czirók, E. Ben-Jacob, I. Cohen, and T. Vicsek, Formation of complex bacterial colonies via self-generated vortices, Phys. Rev. E 54 2 (1996), 1792–1801. [7] P. J. Davis, Circulant matrices, Wiley, New York, 1979. [8] M. R. D’Orsogna, Y. L. Chuang, A. L. Bertozzi, and L. S. Chayes, Self- propelled particles with soft-core interactions: Patterns, stability, and collapse, Physical Review Letters 96 (2006), 104302. [9] D. Grossman, I. S. Aranson, and E. Ben-Jacob, Emergence of agent swarm migration and vortex formation through inelastic collisions, New Journal of Physics 023036 (2008). 93 Chapter 3. Bibliography [10] C. J. Harvey-Clark, W. T. Stobo, E. Helle, and M Mattson, Putative mating behavior in basking sharks off the Nova Scotia coast, Copeia No. 3 (1999), 780–782. [11] J. Hemmingsson, Modellization of self-propelling particles with a coupled map lattice model, J. Phys. A 28 (1995), 4245–4250. [12] A. Huth and C. Wissel, The simulation of movement of fish schools, J. Theor. Biol. 156 (1992), 365–385. [13] V. L. Kulinskii, V. I. Ratushnaya, A. V. Zvelindovsky, and D. Bedeaux, Hydrodynamic model for a system of self-propelling particles with con- servative kinematic constraints, Europhys. Lett. 71 (2) (2005), 207–213. [14] H. Levine, W.J. Rappel, and I. Cohen, Self-organization in systems of self-propelled particles, Physical Review E 63 (2001), 017101. [15] Y.-X. Li, R. Lukeman, and L. Edelstein-Keshet, Minimal mechanisms for school formation in self-propelled particles, Phys. D. 237(5) (2008), 699–720. [16] R. Mach and F. Schweitzer, Modeling vortex swarming in daphnia, Bull. Math. Bio. 69(2) (2007), 539–562. [17] J. A. Marshall, M. E. Broucke, and B. A. Francis, Formations of vehicles in cyclic pursuit, Automatic Control, IEEE Transactions on 49 (11) (2004), 1963–1974. [18] A. Mogilner, L. Edelstein-Keshet, L. Bent, and A. Spiros, Mutual in- teractions, potentials, and individual distance in a social aggregation, J Math Biol 47 (2003), 353–389. [19] H.-S. Niwa, Self-organizing dynamic model of fish schooling, J theor Biol 171 (1994), 123–136. [20] , Newtonian dynamical approach to fish schooling, J theor Biol 181 (1996), 47–63. [21] ,Migration of fish schools in heterothermal environments, J theor Biol 193 (1998), 215–231. 94 Chapter 3. Bibliography [22] A. Okubo, Diffusion and ecological problems: Mathematical models, Springer Verlag, New York, 1980. [23] A. Okubo, D. Grunbaum, and L. Edelstein-Keshet, The dynamics of an- imal grouping, Diffusion and Ecological Problems: Modern Perspectives (A Okubo and S Levin, eds.), Springer, N.Y., 2001. [24] A. E. Parr, A contribution to the theoretical analysis of the schooling behaviour of fishes, Occasional Papers of the Bingham Oceanographic Collection 1 (1927), 1–32. [25] J. Parrish and L. Edelstein-Keshet, Complexity, pattern, and evolution- ary trade-offs in animal aggregation, Science 284 (1999), 99–101. [26] J. Parrish, S. Viscido, and D. Grunbaum, Self-organized fish schools: an examination of emergent properties, Biol. Bull. 202 (2002), 296–305. [27] V. I. Ratushnaya, D. Bedeaux, V. L. Kulinskii, and A. V. Zvelindovsky, Stability properties of the collective stationary motion of self-propelling particles with conservative kinematic constraints, J. Phys. A 40 (2007), 2573–2581. [28] S. Sakai, A model for group structure and its behavior, Biophysics Japan 13 (1973), 82–90. [29] T. C. Schneirla, A unique case of circular milling in ants, considered in relation to trail following and the general problem of orientation, Amer- ican Museum Novitates 1253 (1944), 1–25. [30] J. R. Silvester, Determinants of block matrices, Maths Gazette 84 (2000), 460–467. [31] R. Suzuki and S. Sakai, Movement of a group of animals, Biophysics Japan 13 (1973), 281–282. [32] C. Topaz and A. Bertozzi, Swarming patterns in a two-dimensional kine- matic model for biological groups, SIAM J. Appl. Math. 65 (1) (2004), 152–174. [33] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Phys. Rev. Lett. 75(6) (1995), 1226–1229. 95 Chapter 3. Bibliography [34] D. Weihs, Energetic advantages of burst swimming of fish, J. theor. Biol. 48 (1974), 215–229. [35] S. G. Wilson, Basking sharks schooling in the southern Gulf of Maine, Fish. Oceanogr. 13 (2004). 96 Chapter 4 A Field Study of Collective Behaviour in Surf Scoters: Empirical Methods 3 4.1 Introduction Understanding mechanisms governing collective motion in animal groups re- quires that empirical data be collected. Empirical data provides information against which model hypotheses can be tested. Unfortunately, data on mov- ing animal groups is quite scarce due to a host of difficulties associated with gathering such data in a reliable manner. My goal in this part of the research is to obtain real field observations of animal groups, and infer what are the rules that individuals follow, and what are effective interaction forces that govern their associations. To obtain this data, I develop an experimental method to gather position and velocity data on flocks of surf scoters swim- ming collectively. Data is gathered by photographing flocks from an overhead position in time series. The methods described in this chapter yield a data set that offers dynamic trajectory data of individuals within groups almost an order of magnitude larger than in previous work. 4.1.1 Difficulties in Obtaining Data Obtaining empirical data in the field is complicated by the typically three- dimensional nature of animal aggregations (resulting in occlusion of interior individuals), the difficulties of calibrating measurement equipment at vary- ing locations, high speeds of movement, and the transient behaviour (both in space and time) of animals in the field in general. Overcoming these diffi- culties by moving to a controlled laboratory setting places restrictions on the 3A version of this chapter will be submitted for publication as ‘A field study of collective motion in surf scoters’, R. Lukeman, Y.-X. Li and L. Edelstein-Keshet. 97 4.1. Introduction spatial extent and size of groups, and (depending on the species and experi- mental setup) can create difficulty in obtaining natural behaviours within an artificial environment. Furthermore, for many species (e.g., flocking birds), a laboratory setting is not feasible. 4.1.2 Previous Work Nevertheless, significant efforts have been made to record positions and move- ments of individuals within groups. To reconstruct three dimensional (3D) positions of individuals, researchers have used three primary methods: stereo photography [1, 2, 3, 5, 8], the shadow method, [9, 10, 11, 12] and orthogo- nal photography [6, 7, 13]. Stereo photography uses two parallel cameras, mounted on a beam with some fixed distance of separation. Parallax mea- surements from photographs then allow 3D positions to be reconstructed. The shadow method tracks both the position of individuals, and also the shadow of individuals resulting from shining a light on the group. By follow- ing both the size and position of the shadow in reference to the position of the camera and animal, 3D positions can be calculated. The orthogonal method is similar to stereo photography, except cameras are aligned orthogonally as opposed to parallel, and 3D positions are reconstructed via the geometry of the two cameras and positions in images. 4.1.3 Tracking In order to obtain dynamic information on animal group motion, individuals must be tracked in time; that is, an individual in one frame must be linked to the same individual in the next frame, and so on, to reconstruct the trajectory for each group member. The ability to track individuals within a group becomes more challenging as group size increases. Thus, there have been a number of empirical studies on small groups (2-30 individuals) with tracking [6, 9, 10, 12, 13, 15], and a number of studies of larger groups without any tracking in time ([8], 25-76 individuals, [3],21-61 individuals, [2], more than 1000 individuals). In contrast, the data set described here consists of trajectory data for all individuals within groups of a few hundred individuals. 98 4.2. Methods 4.1.4 Challenges Avoided in This Study This study avoids a number of the challenges inherent in empirical studies of animal groups. First, individuals moved on the surface of the water, which is essentially two-dimensional. There was therefore no need for multi-camera methods for reconstructing positions. It should be noted that unlike some other studies in two dimensions, individuals were not artificially restricted as such; here it was an artefact of the environment in which they naturally move. Second, individuals swimming at the surface moved at relatively slow speeds as compared to, for example, flying birds or swimming fish. Low speeds of movement reduce temporal resolution required for tracking, allow- ing individuals to be successfully tracked using frame rates typically avail- able on consumer cameras. Third, the well-defined spacing of the group, and the high contrast provided by individuals against the water surface, permit mostly automatic detection of individuals in images. Last, the surf scoter groups in this study were located in an area that provided both a convenient overhead location to photograph the group without disturbing their natural behaviour, and also features that facilitated accurate image calibration for a relatively large field area. 4.2 Methods 4.2.1 Location and Materials Data was gathered March 1-12, 2008 in waters adjacent to the Vancouver, BC waterfront, in Burrard Inlet. Overwintering surf scoters were observed foraging near dock pilings around Canada Place, a commercial facility. An elevated promenade around Canada Place provided overhead locations to photograph surf scoter groups swimming on the water surface. A Nikon D70s DSLR camera was used, attached to a Manfrotto 190XPROB tripod. Images were taken in continuous mode, fired at 3 frames per second to generate a time series of surf scoter movements. A Nikon AF-S Nikkor 18-70mm ED lens was used, fixed at the maximal focal length (70mm). Although the Nikon D70s has a maximum resolution of 2000 x 3008 pixels, images were taken in reduced size (1000 x 1504 pixels) to ensure accurate frame rates in continuous mode, and to reduce computation times in data analysis. Aperture was fixed at f4.5, shutter speed varied to optimize exposure (between 1/8000 sec. and 1/250 sec.), and auto-focus was used. 99 4.2. Methods 4.2.2 Surf Scoter Behaviour in Field Study Surf scoters are known to collectively forage primarily for mollusks (e.g., mussels, clams) [14] via synchronous underwater dives. In the field location of this study, overwintering scoters gathered in groups of approximately 200-400 individuals in open water. Upon initiating a dive sequence, the group swam on the water surface from open water towards the dock area surrounding Canada Place, presumably to forage on mussels on or around the dock pilings. After this period of movement to the dock area, individuals then dove in a highly synchronous manner to forage. After a period of foraging underwater, individuals re-emerged at the water surface, again in a highly synchronous manner, and returned back to open water. This process repeated continually throughout the observation period. 4.2.3 Experimental Technique Because the goal of data collection was to track positions of individuals through time, it was necessary to fix the camera elevation, angle, and image frame during data collection events. To this end, the camera was positioned so that the bottom of the viewfinder image aligned with the outer edge of the dock (see Fig. 4.1). Because the boardwalk elevation and dock width were fixed in the region of data collection, the camera elevation and angle were fixed, and easily measured. Furthermore, this set-up allowed the cam- era apparatus to be moved along the boardwalk for a given sequence (but not during), so as to increase the likelihood of encountering individuals in the image frame. Once a group of surf scoters began moving in for a foraging event, the camera was positioned as described above, and upon encountering individuals in the image frame, the camera was fired continuously at 3 fps until individuals had completely left the frame, either through a dive, or by swimming outside the image region. A similar technique was used to capture dive-return movements, although these events were captured much less fre- quently due to the difficulty of estimating the location of re-emergence from a dive, as compared to a group moving in on the surface. 4.2.4 Calibration and Testing Ideally, aerial images are taken directly overhead of the object; in this study, the camera axis was not vertical, but instead approached the water surface at 100 4.2. Methods water dock camera 17.5m 9.0m θ φφ 0 y L Figure 4.1: A schematic diagram of the experimental setup. θ is the camera axis angle, while φ is the angle-of-view of the camera lens. an angle. In this type of photography (oblique aerial photography), positions within the camera frame must be transformed to real-world positions using known measurements of the apparatus setup and camera specifications. Ad- ditionally, commercial cameras (such as the one used in this study) are not designed for precise data collection, and thus, manufacturer’s specifications must be tested for accuracy [4]. Fig. 4.1 gives a schematic diagram of the setup. The camera was at a height of 17.5 m, and the dock was 9 m in width. Because the bottom of the image frame is aligned with the outer dock edge, these two distance measurements allow us to easily calculate the camera axis angle: θ = arctan ( 9.0 17.5 ) = 0.475 rad = 27.2o. The angle-of-view, φ can be calculated from the lens focal length f and the 101 4.2. Methods vertical dimension of the image sensor, dv as 2 arctan ( dv 2f ) . The Nikon D70s has an image sensor of 23.7 x 15.6 mm, and the lens used had a focal length of 70 mm, and thus, we calculate 2 arctan ( 15.6mm 140mm ) = 0.222 rad = 12.7o. However, because the true focal length of a consumer lens can vary from the specified length [4], the angle-of-view was also measured for the camera and lens used in this study by photographing known lengths from a variety of camera axis angles. These measurements gave a true angle-of-view of 0.235 rad = 13.5o. In our image post-processing, the bottom 28 pixels (of 1000) are removed to eliminate the edge of the dock that may appear in the image due to imperfect alignment in the field. Thus, the image size is now 1504 x 972 pixels, and the angle-of-view is φ = 0.972 · 0.235 = 0.229 rad (or 13.1o). To convert positions from an image frame in our setup to real-world po- sitions, two transformations must be made: a vertical transformation due to the camera axis angle, and a horizontal transformation due to projective perspective distortion: lines that are parallel in reality will converge in an image taken at some non-zero camera axis angle. Denote as φ̂ the angle cor- responding to real distance y, according to Fig. 4.2. In the image frame, as φ̂ ranges through [0, φ], pixels range (linearly) through [0, 972]. Using simple trigonometry, real vertical position y in pixels can be recovered by y = 972 φ ( tan ( θ + pφ 972 ) − tan(θ) ) , (4.1) where p is the vertical pixel in the image; i.e, p = 972φ̂/φ. Substituting p = 972 into the above formula gives L, the real vertical extent of the image (in pixels) as L = 1421. Images taken in the field in this study differ from directly overhead images only by rotation in one axis. Borrowing from flight dynamics, the pitch of the camera is changed, while the roll and yaw is fixed at 0 (controlled by the tripod head). As a result, the horizontal perspective distortion in the images is symmetric about the horizontal centre of the image, where there 102 4.2. Methods θ φ y0 φ L Figure 4.2: A schematic diagram of the setup with angles used in vertical transformation. φ̂ represents the angle corresponding to real distance y (in pixels). L is the real vertical extent of the image (in pixels). 103 4.2. Methods x y Horizontal transformation matrix 500 1000 1500 200 400 600 800 0 500 1000 1500 Figure 4.3: The horizontal calibration matrix, associating each (x, y) image coordinate with the horizontal real pixel value. is no perspective distortion. This distortion then increases linearly to the edges. To obtain the horizontal transformation, the camera used in the study photographed a grid at camera axis angle θ, to reproduce the field conditions. In order to deduce the transformation, we need only know the ratio of the real length of the top of the image frame to the bottom of the image frame, which gives the maximal perspective distortion. This ratio was calculated for a number of trials using a ruler to measure real length at each edge of the frame, giving a ratio of 1.197. Using this value, a calibration matrix is constructed giving the horizontal map from image pixels to real pixels, shown in Fig. 4.3. Combining the vertical transformation in (4.1) with the horizontal trans- formation shown in Fig. 4.3, positions in an image are transformed into real positions in space, in units of pixels. In order to test that these spatial trans- formations accurately reconstruct real positions, we photograph a grid with a camera with camera axis angle θ (Fig. 4.4). Grid locations, equidistant 104 4.2. Methods Figure 4.4: The calibration grid used to test the image transformation. The pixel location of the upper-right corner of each grid vertex was marked and transformed in MATLAB (see Fig. 4.5). in reality, will have perspective distortion in the image. These locations are marked on the image, and transformed accordingly. As shown in Fig. 4.5, the transformed positions accurately reconstruct the original grid. We note that there is also a very small amount of pincushion distortion which we assume is negligible. The other camera specification to be tested was the frame rate while shooting in continuous mode. The camera specifications give a rate of 3 fps. This was tested by taking 60 consecutive images of a stopwatch in continuous mode, at a shutter speed of 1/250 sec., using the same image resolution as in the field study. The camera proved to be remarkably accurate in this respect, as time between capturing images was 1/3 sec. (built-in frame rate) + 1/250 sec. (shutter speed) with error less than 1%. Due to the relatively low resolution of the images taken in this study, images could be written to storage as quickly as they were taken, thus the camera storage buffer was never compromised and an accurate frame rate was maintained. 105 4.2. Methods 0 500 1000 15000 200 400 600 800 1000 1200 1400 x y Transformed grid points Figure 4.5: 42 Grid points marked in Fig. 4.4 plotted as ‘x’ marks, with reconstructed positions as ‘o’ marks. 106 4.2. Methods Figure 4.6: An example of one image from a time series of images collected in the field. 4.2.5 Postprocessing Images Data was stored as a series of JPEG images (see Fig. 4.6 for an exam- ple). These images were then processed in MATLAB to isolate individual ducks and estimate the center of mass of each individual. The result is a set of (x, y) coordinate positions characterizing each individual in a frame, for every frame. These positions were then linked frame-by-frame, to create trajectories. We describe this process below. 4.2.6 Extracting Positions Image processing was performed in MATLAB. JPEG images are stored as an m x n x 3 matrix; i.e., 3 color layers (red, green, blue) ofm x nmatrices, where m and n are the horizontal and vertical pixels (1504 and 972 respectively, in this study). Each color value ranges from 0 to 255, indicating the amount of a particular color at a given pixel location. Consistent color differences in images gathered in this study were exploited as a first step in isolating individuals (in this case, ducks tended to have red pixel values greater than blue and green). Images were then thresholded to black and white (Fig. 107 4.2. Methods 4.7.a). Next, morphological operations were performed on the thresholded image, which identified and reinforced objects within the image (Fig. 4.7.b). Then, the thresholded image was overlayed with the original (Fig. 4.7.c), so that the automatic object detection could be manually compared with the original image to detect any errors. In the example in Fig. 4.7, one individual was missed by the algorithm, while a few others are represented by two objects. At this stage, an image editing user interface (purpose-built in MATLAB for this study) was used to fix any errors in the automatic object detection algorithm. Operations included joining objects, marking missed individuals, separating objects that represent two individuals, and modifying objects that poorly represent the center of an individual. There was typically minimal manual editing due to the well-spaced nature of the groups in this study. Following the editing step, the center of mass of each object was calculated, and plotted (Fig. 4.7.d; centers of mass plotted in red). Last, the center-of-mass positions were transformed to real pixel positions according to the perspective transformations described in Section 2.4, and saved. 4.2.7 Tracking Following the image processing of images, the data is reduced to a set of positions for each frame. These data were exported to ImageJ, an image editing software package. The ImageJ ‘Particle Tracker’ plugin was used to associate a position in one frame to the corresponding position in the next frame. Repeating this process through successive frames generated a trajectory through time for an individual, and results were imported back into MATLAB, where speed and heading were calculated from the trajectory data. However, the tracking step presented a number of technical difficulties. Particle tracking algorithms function best when the maximum distance traveled between frames is less than half the minimum distance between individuals. The camera used in this study limited our frame rate to 3 fps, which was not fast enough to ensure this bound on distance traveled between frames. Also, particle tracking codes are generally written for motion without directional bias; i.e., assuming a random-walk type of motion. However, motion of individuals in this study was generally highly polarized, resulting in directed motion. As a result of these technical difficulties, trajectories (and thus velocities) were partially incorrect, and incomplete (see Fig. 4.8). A typical mistake by the algorithm was to associate a position in one frame with 108 4.2. Methods (a) Image after thresholding. (b) Image after morphological oper- ations. (c) Image after overlay. (d) Image after manual editing and center-of-mass marking. Figure 4.7: Processing an image to obtain positions. 109 4.2. Methods a position in the next frame in the opposite direction of travel (i.e., behind the individual). Although it is intuitively clear that this is exceedingly unlikely for individuals in this study, the tracking algorithm cannot incorporate this bias. To overcome this difficulty, an algorithm was developed to correct tracking errors, outlined below. 1 Trajectories are calculated with particle-tracking software, filtered to retain realistic sections (i.e., changes in position and heading less than some threshold), while unrealistic sections are discarded. 2 In a given frame, all individuals associated with a good trajectory are assigned a velocity based on the trajectory. 3 For individuals not associated with a good trajectory, velocity is esti- mated by calculating the mean velocity of neighbours in a local region. 4 The entire sequence is re-tracked, using the estimated velocities as pre- dictions for the following step, and building the trajectory based on the closest individual to the predicted location. After applying this algorithm, trajectory reconstruction was virtually complete (less than 1% of individuals remained without association to a real- istic trajectory). Fig. 4.9 shows an example frame with velocities that have been obtained after the re-tracking step - accuracy is considerably improved as compared to Fig. 4.8. It is useful to note that for the polarized behaviour exhibited by scoter flocks in this study, only a fraction of velocities need to be recovered during the initial tracking step. These first-pass velocities pro- vide a basis for estimating the expected position in the next frame. Then, velocities (and successively, trajectories) can be calculated exactly by find- ing the nearest actual position to the expected position in the next frame. Trajectories calculated after re-tracking were not only more accurate, but also were considerably longer on average, having fewer interruptions due to tracking deficiency. Examples of trajectory results are shown in 4.10, showing convergence, turning and splitting of groups. 4.2.8 Edge Effects When performing statistical analysis of the data collected in this study, it is desirable in certain cases to exclude contributions from edge individuals due to skewing of statistical measures. For instance, relative location of 110 4.2. Methods 0 500 1000 15000 200 400 600 800 1000 1200 1400 Figure 4.8: An example frame with reconstructed velocities (grey), showing partially incorrect and missing velocities. In reality, individuals are highly polarized in this frame. 0 500 1000 15000 200 400 600 800 1000 1200 1400 Figure 4.9: The example frame of Fig. 4.8 showing velocities (grey) obtained by re-tracking the event. 111 4.2. Methods 0 500 1000 15000 200 400 600 800 1000 1200 1400 (a) group converging 0 500 1000 15000 200 400 600 800 1000 1200 1400 (b) group turning 0 500 1000 15000 200 400 600 800 1000 1200 1400 (c) section of the group turns 0 500 1000 15000 200 400 600 800 1000 1200 1400 (d) group splitting Figure 4.10: Example trajectories for 4 groups. Starting positions are plotted in green, final positions in red. 112 4.2. Methods Figure 4.11: In the first frame, the individual at the origin has nearest neigh- bours filling each quadrant, and so is not an edge individual. In the second frame, the individual has no nearest neighbours (of the first 8) in the third quadrant, and thus is considered an edge individual. nearest neighbours all tend to lie on one side of an edge individual. In order to identify and exclude edge individuals from such statistical measures, the following criteria were used. The location of the first 8 nearest neighbours for each individual was calculated. Then, any individuals who do not have at least 1 nearest neighbour in each of the 4 quadrants relative to the individual (see Fig. 4.11) are identified as edge individuals. Also, if the distance to the nearest neighbour in each quadrant is beyond a threshold distance, the (isolated) individual is also identified as an edge individual. Fig. 4.12 shows a plot indicating edge individuals identified within a group. Based on group shape, some sequences were excluded altogether, if a disproportionately large number of individuals were found on the edge (e.g., quasi-1D groups). Further, because edge individuals are not constrained in the same way as interior individuals, these individuals may exhibit behaviors other than those seen in interior individuals. A separate comparative analysis of edge individuals may prove useful in elucidating any differences. 4.2.9 Processed Events We categorize data into sequences, referring to a series of images taken in suc- cession of a particular group movement (approach, return, etc.). A total of 14 sequences were processed, with frames per sequence ranging from 22 to 137. 828 frames were analyzed in all, and a total of 75269 positions were analyzed over all frames (42599 after eliminating edge individuals). Sequence-specific 113 4.2. Methods 0 500 1000 15000 500 1000 Figure 4.12: A plot of a group with edge individuals (as determined by the algorithm outlined in the text) in red, with all other individuals in black. details are listed in Appendix C.1. 4.2.10 Body Alignment Versus Velocity It is normally reasonable to assume that the velocity of an individual is in the same direction as the vector describing the body alignment of the individual. However, when there are environmental influences such as currents, discrep- ancies can arise between alignment and actual velocity. Raw data in this study has a convenient qualitative indicator of the presence of currents: dur- ing collective motion, individual scoters are observed to occasionally excrete waste. After excretion, the waste can be tracked in successive time frames - the direction and speed of this drift provides quantitative information on currents. In this study, currents predominantly were directed parallel to the dock, so we neglect any current components perpendicular to the dock. To quantify currents, we first measure the left/right drift of waste in a test frame from each event (in body lengths per second), giving a relative measure of the strength of the current for each event. Then, a frame from each event is analyzed for body alignment - the alignment of each individual is marked manually, and the average body alignment of all individuals in a frame is 114 4.2. Methods compared to the average velocity for the chosen frame - any discrepancy is assumed to be due to currents. The current-induced discrepancy is a function of the direction of travel - clearly there is little alignment-velocity directional discrepancy if an indi- vidual is moving in the same direction as the current, as opposed to moving perpendicular to the current. Therefore, the current is quantified as a vec- tor which sums with the normalized alignment vector to give the normalized velocity vector. To this end, let b = (bx, by) be the vector of body align- ment of an individual, let v = (vx, vy) be the (normalized) velocity vector (which contains the effects of any current), and let c = (cx, cy) be the current vector that accounts for the discrepancy in alignment and velocity, so that (v + c)‖b. Writing in component form gives [ vx + cx vy ] = [ bx by ] , where we have set cy = 0 due to the assumption of along-dock currents. Cross multiplying and equating the first component gives cx = bxvy − byvx by . (4.2) For each sequence, one frame is analyzed, and an average alignment vector b̂ is calculated by manually marking the body alignment vector for each individual, and then averaging over all individuals in the frame. Also, the velocity of each individual is known via the tracking algorithm, from which an average velocity v̂ for the frame is calculated. These values are used to calculate cx via Eq. (4.2). In this analysis, we assume currents to be constant throughout a given sequence (the longest of which lasts approximately 45 seconds). Analyzed data was collected over 4 days in 2008: March 1 (29 minutes), March 6 (60 minutes), March 7 (107 minutes) and March 12 (50 minutes). Sequence-specific details on currents are given in Appendix C.2. At this stage, alignment data for each event could be calculated from the cx values. However, the waste (tracer) drift in frames is a much more direct indicator of currents (i.e., less prone to measurement error), and so instead, a relationship between observed tracer velocity and cx is established via a linear least-squares fit, and this linear function is then used to compute the resultant current vector applied to velocities. In Fig. 4.13, data for tracer measurements and cx calculations are shown, with the linear fit. 115 4.3. Results 0 0.2 0.4 0.6 0.8−0.6 −0.4 −0.2 0 0.2 current vector vs. tracer shift tracer shift, BL/s cu rr e n t v ec to r Figure 4.13: Current vector x-value cx vs. waste tracer drift speed, from raw images, with linear least-squares fit. 4.3 Results 4.3.1 Basic Statistics We report basic statistics on observed surf scoter flocks in this study. Units of distance are body-lengths (BL), calculated from an average of individual scoters in images. In the analyzed sequences, surf scoters within groups had a mean nearest-neighbour distance of 1.435 ± 0.251 BL, and moved at an average speed of 2.007 BL/sec. Spacing and velocity was relatively similar across sequences, although sequences featuring the group returning from the dock to open water (S-3,S6b-c,S-12) showed a higher average speed (2.85 BL/sec, 2.45 BL/s, 2.67 BL/s, respectively). In S-3, the group was returning from a dive and had a considerably lower nearest-neighbour distance (aver- age of 1.15 ± 0.18 BL). It was generally noticed that immediately following resurfacing from a dive, aggregates were tighter than those approaching the dock. The sequence S-6a featured individuals converging into a group, and so individuals had not achieved equilibrium spacing, resulting in larger nearest- neighbour spacing (1.71 ± 0.26 BL). 116 4.3. Results Table 4.1: Basic statistics: data sequences Sequence avg. speed NND std(NND) S-1 2.150 1.538 0.1979 S-2 2.129 1.477 0.2156 S-3 2.847 1.148 0.1836 S-4 1.741 1.345 0.2135 S-5 1.985 1.323 0.1883 S-6a 1.758 1.710 0.2606 S-6bc 2.449 1.482 0.2214 S-7 1.891 1.623 0.2971 S-8 1.376 1.486 0.2335 S-9 2.006 1.5126 0.2402 S-10 1.900 1.4304 0.2366 S-11 2.123 1.3519 0.2037 S-12 2.672 1.4463 0.2426 S-13 2.080 1.3208 0.1942 S-14 2.715 n/a n/a total: 2.007 1.435 0.2506 Table 4.2: Average speed and nearest neighbor distances for the 14 separate sequences of data that were analyzed. Units are body-lengths (BL). Speed and NND was averaged over all individuals in all frames. The last column is the standard deviation of NND. We summarize basic statistics for each event, where distance is given in units of body lengths, for all individuals, and excluding edge individuals. 4.3.2 Nearest-Neighbour Distance Distributions For each individual within a group, its nearest neighbours are ordered with respect to distance, and labeled in increasing order of distance as NN1, NN2, NN3. . .. Such ordered sequences of neighbours are calculated for all individ- uals over all time frames, and frequency distributions of successive nearest 117 4.3. Results neighbours are constructed. Distributions of successive nearest neighbours are plotted in Fig. 4.14. Distribution means and standard deviations increase nearly linearly with successive neighbours, as seen in Fig. 4.15. Motivated by the regularity of distributions in Fig. 4.14, a functional form is sought to account for nearest-neighbour distributions observed in scoter flocks. In [16], a probability density function (pdf) for nearest-neighbor dis- tance is derived in the context of biological entities forming a (static) spatial pattern via competition for space (e.g., placement of trees in a forest via competition for light and nutrients, etc.) Although the context of the deriva- tion is different, the driving principles are very similar. In the case of scoter flocks, the use of space is not for resources, but instead to maintain a minimal distance from neighbors so as to avoid collisions. In [16], space usage of an individual is modeled as a Gaussian distribution, where overlaps with nearby neighbors present competition for that space, analogous to a repulsion force acting between individuals in a group. In this framework, a pdf is derived, using methods from information theory, for observed nearest-neighbor dis- tance. A good fit between the derived form of the pdf, and empirical data on spacing in scoter flocks, suggests that the spacing properties of groups in our data set are part of a more general class of spatial patterns in biology resulting from interactions among individuals under spatial competition (or equivalently, attractive and repulsive interactions). The pdf has the form q(d) = rC exp ( −β exp (−d 2/4w2) 1− exp (−d2/4w2) − ǫπd 2 ) , (4.3) where w represents the effective size of an individual, ǫ is the number density of the group (that is, number of individuals per unit of area), β is a fitting parameter, and C is a scale factor such that ∫ ∞ 0 q(d)dd = 1. (4.4) To match nearest-neighbor distributions to this derived form, we first normal- ize nearest-neighbor distributions for the observed data so that the distribu- tion curve of observed data integrates to 1. Then, ǫ is estimated by counting individuals in a fixed area (within the flock) across a number of trials, giving ǫ = 0.3358. Next, a curve-fitting algorithm is used to fit β and w, and finally C is calculated by satisfying Eq. (4.4). We note that the parameter w is 118 4.3. Results 0 1 2 3 4 50 1000 2000 3000 4000 Distributions of NN1−NN8 distances NND, in BL fre q. Figure 4.14: Nearest-neighbour distributions for the first 8 nearest neigh- bours. Successive means are also plotted (◦) above distributions. typically derived from the data, as opposed to fitting it, but because body sizes of individuals in this study are not circular, we leave this parameter to be fitted. Results of the curve-fitting give optimal parameters β = 2.44 and w = 0.9743. Because units of distance are in body lengths, and given that body width is slightly less than body-length, a value of w slightly less than 1 is reasonable. Fig. 4.16 shows the observed nearest neighbor distribution to- gether with the fit to q(d). The nearest-neighbour distribution data was also fitted to a Gaussian distribution, resulting in a poorer fit than Eq. (4.3), as the right-hand tail of the data distribution is longer than the left-hand tail. Although the functional form of 4.3 was derived only for the first nearest- neighbor, each successive neighbor distribution can be accurately fitted by varying β and w. However, there is currently no theory to explain successive distributions of interacting particles (whereas functions have been derived for distributions of successive neighbors for random points). 119 4.3. Results 0 2 4 6 81 1.5 2 2.5 3 3.5 Mean distance, NN1−8 NN number di st an ce , B L (a) Mean of successive nearest neighbours 0 2 4 6 80.2 0.3 0.4 0.5 0.6 0.7 Standard deviation of NND, NN1−8 NN number st an da rd d ev ia tio n (b) Standard deviation of successive nearest neighbours Figure 4.15: 0 1 2 30 0.5 1 1.5 2 Theoretical and observed NND distribution (NN1) distance, d, in BL q(d ) Figure 4.16: Nearest neighbour distributions of the first neighbour, with fitted probability density function q(d) (dashed) overlaid. 120 4.4. Conclusion 4.4 Conclusion Empirical data on animal groups has previously been limited to static anal- yses of large groups, or dynamic analyses of a few tens of individuals. In this work, we have extended such empirical studies to dynamic descriptions of relatively large (a few hundred) individuals. This data was gathered by photographing in time series surf scoters collectively swimming in waters near Canada Place, Vancouver, BC. Individuals moved in well-spaced two- dimensional flocks with relatively low velocities, thus overcoming many of the traditional obstacles associated with empirical studies of animal group motion. Scoter flocks foraged in predictable locations near a dock with an elevated area that allowed for oblique overhead photography. Corrections due to the angle of photography were easily made by measured features of the dock and building where the study took place. Data was gathered with a digital SLR camera firing at 3 fps, and resulting images were processed to isolate indi- viduals in MATLAB via software designed for this experiment. A tracking algorithm was developed in conjunction with conventional particle-tracking software to accurately track polarized groups of individuals, and individual trajectories were reconstructed. Methods were developed to isolate interior individuals within a group to facilitate statistical measures of the group, and also to calculate water currents, and eliminate the effects of these currents on individual motion. 14 data sequences were analyzed in total, representing a total of over 75000 reconstructed positions over 828 time frames. Basic statistics, including neighbour spacing, mean velocity, and successive neigh- bour distributions were reported. The basic results presented here are primarily limited to static aspects of the data, and do not make use of the individual trajectories that have been constructed. Related work featured in Chapter 5 of this thesis, however, analyzes dynamic aspects of the empirical data described here. 121 Bibliography [1] I. Aoki and T. Inagaki, Photographic observations on the behaviour of japanese anchovy engraulis japonica at night in the sea, Mar. Ecol. Prog. Ser. 43 (1988), 213–221. [2] M. Ballerini, N Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Gi- ardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic, Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study, Proc. Nat. Acad. Sci. 105 (2007). [3] R. Budgey, Three dimensional bird flock structture and its implications for birdstroke tolerance in aircraft, International Bird Strike Committee IBSC 24/WP 29 (1998). [4] A. Cavagna, I. Giardina, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic, Empirical investigation of starling flocks: a bench- mark study in collective animal behaviour, Anim. Behav. 76 (2008). [5] J. Cullen, E. Shaw, and H. Baldwin, Methods for measuring the three- dimensional structure of fish schools, Can. J. Zool. 71 (1993), 1494–1499. [6] D. Grunbaum, S. Viscido, and J. Parrish, Extracting interactive control algorithms from group dynamics of schooling fish, Cooperative Control LNCIS 309 (2004), 103–117. [7] T. Ikawa, H. Okabe, T. Mori, K. Urabe, and T. Ikejoshi, Order and flexibility in the motion of fish schools, Journal of Insect Behavior 7(2) (1994), 237–247. [8] P. F. Major and L. M. Dill, The three-dimensional structure of airborne bird flocks, Behavioural Ecology and Sociobiology 4 (1978), 111–122. [9] B. L. Partridge, Fish school density and volume, Marine Biology 54 (1979), 383–394. 122 Chapter 4. Bibliography [10] , The effect of school size on the structure and dynamics of min- now schools, Animal behaviour 28 (1980), 68–77. [11] , Internal dynamics and the interrelations of fish in schools, J Comp Physiol 144 (1981), 313–325. [12] B. L. Partridge, T. Pitcher, J. M. Cullen, and J. Wilson, The three- dimensional structure of fish schools, Behav Ecol Sociobiol 6 (1980), 277–288. [13] H. Pomeroy and F. Heppner, Structure of turning in airborne rock dove flocks, Auk 109(2) (1992), 256–267. [14] J.-P. Savard, D. Bordage, and A. Reed, Surf scoter (melanitta perspicil- lata). in the birds of north america no. 363 (a. poole and f. gill, eds.), The Birds of North America Inc., Philadelphia, PA, 1998. [15] S. Viscido, J. Parrish, and D. Grunbaum, Individual behaviour and emer- gent properties of fish schools: a comparison of observation and theory, Mar. Ecol. Prog. Ser. 273 (2004), 239–249. [16] G. Zou and H. Wu, Nearest-neighbor distribution of interacting biological entities, J. Theor. Biol. 172 (1995), 347–353. 123 Chapter 5 Ducks in a Row: Inferring Interaction Mechanisms from Field Data of Surf Scoters 4 5.1 Introduction Across many species, animals often aggregate into cohesive, ordered groups. To understand the mechanisms of interaction governing such groups, a rela- tively large number of hypothetical models have been developed [1, 6, 9, 10, 11, 15, 19, 21, 23]. However, the variety of different models raises an issue: given a particular collective behavior of interest, one must select the proper model, and validate the hypotheses in that model. The key to choosing be- tween competing models and validating a given choice is by comparison to empirical data. Unfortunately, there is very little empirical data with which to make this comparison, due to difficulties in gathering data on large, dy- namic animal groups. The data that has been obtained is limited either by small group sizes [8, 16, 18, 20] or an inability to track individuals through time [2], resulting in only static data (which limits comparison to a dynamic model). Model hypotheses are typically founded on a series of interaction forces: repulsion away from individuals very close, attraction to individuals far away, and alignment with nearby neighbours. The wide variation in models arises in the particular implementation of these forces. Here, an attempt is made to build a model to match empirical data by proposing and testing a series of hypotheses. 4A version of this chapter will be submitted for publication as’ Ducks in a Row: Infer- ring Interaction Mechanisms from Field Data of Surf Scoters’, R. Lukeman, Y.-X. Li and L. Edelstein-Keshet. 124 5.2. Empirical Methods 5.1.1 Specific Goals of this Work In this chapter, my aim is to draw inferences from the analysis of data col- lected by methods in Chapter 4. I aim to identify the following aspects of interaction. 1. The minimal components, from repulsion, attraction, and alignment, that are necessary to obtain simulation results comparable to observa- tions. 2. Whether there is a hierarchy of interaction rules. 3. The spatial extent of these repulsive, aligning and attractive forces. 4. Angular preferences for neighbors, and how to implement such prefer- ences in a model. 5. The relative weighting of forces that best describes observed data. 5.1.2 Using Spatial Distributions Because interactions are typically modeled as functions of distance between individuals, the spatial structure of individual behavior relative to neighbors presents a key statistic to decipher these interactions. Obtaining such data reliably on animal groups requires a significant spatial extent of the group (i.e., large group sizes, particularly in three dimensions), as otherwise edge effects dominate. Unfortunately, as group size increases, tracking individual trajectories becomes more difficult [5]. Yet, trajectories are needed to de- scribe individual velocities (and differences among mutual velocities), which in highly polarized groups provides a useful indicator of behavioral response. 5.2 Empirical Methods Reconstructing trajectories of animals within large groups in the field is com- plicated by the often three-dimensional nature and associated occlusion that occurs in imaging, high speeds of movement, and difficulty calibrating equip- ment for groups whose location is not highly predictable. To overcome these difficulties, we focused on groups of a few hundred surf scoters, collectively foraging near a dock on the waterfront of Vancouver, British Columbia. Their 125 5.3. Empirical Results collective swimming movement on the water surface during approach to the dock from open water was captured in time series via oblique overhead pho- tography (Figure 5.1.a). These foraging groups moved at relatively slow speeds, were well-spaced, and occupied predictable locations having con- venient calibration properties. The experimental setup thus allowed high quality data to be gathered, permitting reconstruction (via image processing software and particle-tracking software) of individual positions and trajecto- ries (and thus velocities) (Figure 5.1.b). Data was gathered over a two-week period in March, 2008, and 13 sep- arate sequences of group motion were reconstructed. Each sequence was comprised of between 25 and 137 frames taken at 3 fps, with number of individuals per frame ranging from a few individuals to over 200. Overall, over 75000 positions were reconstructed. Unlike [2], individual positions in a frame were completely recovered, which permitted trajectory reconstruction (Figure 5.2). Across all sequences, individuals had a mean nearest-neighbor distance of 1.45 body-lengths (BL) with standard deviation 0.2506 BL. Mean speed was 2.0 BL/s. For complete details on empirical methods, consult Chapter 4. 5.3 Empirical Results 5.3.1 Relative Location of Neighbors We first investigate spatial distributions of nearest neighbors, in two dimen- sional space around individuals. For a given focal individual i placed at the origin, the relative position of nearest neighbours j, ~xj−~xi, is plotted relative to the heading of the focal individual. By summing over neighbours of each individual in each analyzed frame, a two-dimensional distribution emerges (normalized to have maximal value 1), indicating preferred locations rela- tive to neighbors. for computational reasons, only the first 40 neighbors are plotted, a convention we carry through on similar plots both for data and simulations. In Figure 5.3, a circular region of essentially no neighbors indicates a zone of repulsion which maintains the well-spaced structure within groups. Moving out radially, a region of high density indicates preferred location of nearest neighbours (with mean nearest-neighbour spacing 1.44 BL). Further, there is a region directly in front, and directly behind a focal individual, of higher 126 5.3. Empirical Results (a) Raw image x y (b) Reconstruction of position and velocity Figure 5.1: A typical analyzed frame: input image (a) is processed, trans- formed to real space, and individual positions are reconstructed (b). Posi- tions in successive frames are linked using particle-tracking software, giving velocities (b). 127 5.3. Empirical Results 0 500 1000 15000 200 400 600 800 1000 1200 1400 0 500 1000 15000 200 400 600 800 1000 1200 1400 0 500 1000 15000 200 400 600 800 1000 1200 1400 0 500 1000 15000 200 400 600 800 1000 1200 1400 Figure 5.2: Trajectories for 4 sequences, showing surf scoters approaching the dock (y = 0) in a highly polarized fashion. Starting positions are plotted in green, final positions in red. Units are given in pixels, where 1 BL = 46 pixels. 128 5.3. Empirical Results density, indicating a preference for neighbors directly in front and behind, as opposed to alongside. We note that such a plot has an inherent circular reflective property due to mutual neighbour pairs, and so it is reasonable to suggest that the high density of neighbours to the back is a plot artefact (as opposed to an indication of interaction to the back). Moving out radially again, there is an ‘echo effect’ of neighbor density (wherein secondary neighbors appear at double the preferred distance). This results in a periodic structure that is more evident in the frontal direction than sideways, indicating a more structured positioning in the frontal direc- tion. The neighbor distribution indicates that individuals have a tendency to ‘follow the leader’, while maintaining a minimal distance in all directions. Interestingly, the group was occasionally observed to approach the dock in true single-file formation, though such aggregations were not included in dis- tributions. To further characterize the angular bias in neighbours in the radial region of preferred spacing, a mean of distribution values is calculated around an annulus centered on the maximal distribution value with radial width 0.5 BL. To account for the reflective symmetry in distribution plots, data for [0,−180◦] is combined with data for [0, 180◦], where angle 0 corresponds to the right-side of an individual, and 180◦ to the left-side. Fig. 5.4 shows more clearly what is shown in Fig. 5.3 : there is a higher neighbor density centered at 90◦, with an approximate width of 40◦. Beyond this high-density region, fluctuations in neighbor density are small, and likely due to inherent variability in the system. Fig. 5.4 shows more clearly that there is a higher density of neighbors in front. 5.3.2 Relative Deviation of Neighbors Because trajectories of individuals were constructed, we are able to calculate the spatial distribution of the relative heading of neighbours. Similar to Fig. 5.3, a focal individual i is placed at the origin, and at each relative neighbor position j, the absolute value of the difference in heading, |~vj−~vi| is recorded. Averaging over each individual in each frame, a two-dimensional distribution of heading deviation is constructed, shown in Fig. 5.5. At radial distances near the mean nearest-neighbor spacing, deviation is minimal, with slightly lower deviation sideways. However, at distances less than the preferred spacing, deviation is strongly angle-dependent: directly in front (and behind), deviation is high, while at alongside positions, deviation is 129 5.3. Empirical Results Figure 5.3: Neighbour distributions from empirical data, normalized to have maximal value 1. Neighbor positions are calculated relative to the heading of a focal individual (in the direction [0 1] in the plot, indicated by the white graphic). A frontal bias is seen in neighbor positioning. 130 5.4. Building a Model 0 50 100 1500.5 0.6 0.7 0.8 0.9 1 angular neigh. density at preferred distance (data) relative angle, in degrees de ns ity Figure 5.4: Average density of neighbours in a circular shell at the preferred distance, as a function of angle. A distinct region of higher density is centered at 90◦. low, indicating a collision avoidance mechanism whereby individuals deviate to the side if too close to a frontal neighbor, but do not strongly deviate when to close to the side of neighbors. Moving out radially, deviation increases beyond the preferred distance, due to a decrease in correlation with neighbor heading with distance. Taken together, Figures 5.3 and 5.5 suggest that the dominant distance- dependent interaction occurs in the frontal direction : individuals modify po- sition to move in line with individuals to the front, while less-than-preferred inter-individual distances to the front result in a strong repulsive reaction via a deviation to the side. In contrast, alignment with neighbors is strongest to the side. We summarize this information in Figure 5.6. 5.4 Building a Model Because real animal groups do not strictly adhere to rules in the way that simulated groups do, we cannot hope to exactly match observed trajectories of individuals within the group with a model, and further, the more aspects 131 5.4. Building a Model Figure 5.5: Spatial distribution of heading deviation, from empirical data. Radial bands are plotted at 1 BL, 2 BL, and 3 BL for reference. The region where deviation is 0 corresponds to the repulsion zone where no individuals are found. 132 5.4. Building a Model Figure 5.6: A schematic diagram partitioning local space around a focal individual into sectors of location preference, high tendency to align, and high tendency to deviate rapidly (move away from) the focal individual. Schematic was assembled according to trends shown in the data in Figs. 5.3 and 5.5. of the data one attempts to recover in a model, the more complicated the model becomes. Our goal here, in contrast, is to build as simple a model as possible to recover one particular aspect of the observed data: the spatial distribution, (and in particular the frontal bias) of neighbors. Our approach is to compare analyzed data to models that incorporate increasing levels of complexity. The structure of relative position and alignment in Figures 5.3 and 5.5 motivate a zonal modeling framework similar to [6, 10], where space around an individual is partitioned into radial bands of repulsion, alignment, and attraction. We note that the relatively fixed density across all data sequences precludes any hypothesis testing of a topological versus metric interaction, as in [2]. 5.4.1 Model Framework We formulate an acceleration analogue to the fixed-speed zonal model in [6]. n individuals with position ~xi and velocity ~vi (i = 1, . . . , n) are modeled as particles obeying equations of motion, according to ~̇xi = ~vi, (5.1) ~̇vi = ~fi,aut + ~fi,int + ~ξi, (5.2) 133 5.4. Building a Model where ~fi,int corresponds to social forces of interaction with neighbors, ~fi,aut to autonomous forces due to influences other than interaction, and ~ξi to (Gaussian) noise. Writing ~ξi = ωξ ~̂ξi, where ~̂ξi has mean 0 and standard deviation 1, strength of noise is controlled by the parameter ωξ. We set ~fi,aut = ~a − γ~vi, as in [12, 13, 19], where γ is a drag coefficient, and ~a is an intrinsic direction of travel (in the context of our data, the approach to the dock to forage). An upper limit of 10 is set on number of interacting neighbors. Model I: Repulsion Only We set ~fi,int = ωrep ~fi,rep, where ωrep is a weighting parameter and ~fi,rep = 1 nrep nrep∑ j=1 g(|~xj − ~xi|) ~xj − ~xi|~xj − ~xi| , j ∈ R, (5.3) where nrep is the number of neighbors sensed in the repulsion zone R (Fig. 5.8.a), and g(x) is as in Fig. 5.7. We set rrep = 1.44, the radius of the repulsion zone, corresponding to the mean spacing observed in the empirical data. The neighbor distribution in Fig. 5.8.b shows that repulsion alone cap- tures the well-spaced property of the observed data, with highest neighbor density corresponding to rrep. However, as one expects, no angle preference to neighbors occurs with only repulsion, as the force acts equally in all di- rections, as shown in Fig 5.8.c. Thus, Model I is inadequate. Although not able to fully capture observed neighbor distributions, clearly repulsion is an important component of the interaction. Model II: Repulsion and Attraction We add a long-range attractive force to ~fi,int, so that ~fi,int = ωrep ~fi,rep + ωatt ~fi,att, where ωatt is a weighting parameter and ~fi,att = 1 natt natt∑ j=1 g(|~xj − ~xi|) ~xj − ~xi|~xj − ~xi| , j ∈ ATT. (5.4) Here, ATT is the region of attraction, as in Fig. 5.9.a, with width ratt. To eliminate the artificial case of attraction and repulsion forces canceling, we impose a hierarchy of interaction such that if nrep > 0, then ~fi,att = 0. 134 5.4. Building a Model AL ATTR x g(x) 1 -1 rrep ral ratt Figure 5.7: Attraction-repulsion function g(x) is negative in repulsion region R, zero in alignment region AL, and positive in attraction region ATT. Mag- nitudes of attraction and repulsion are controlled by weighting parameters ωatt and ωrep. Thus, repulsion is the dominant interaction, and attraction occurs only if no neighbors are within the repulsion zone. With the inclusion of attraction, groups are tighter, shown by a less diffuse neighbor distribution in Fig. 5.9.b, yet still no angle preference exists for neighbors (Fig. 5.9.c) and so Model II is inadequate. Model III: Repulsion, Neutrality, Attraction Following [9, 20], we include a neutral zone between repulsion and attraction zones (Fig. 5.10.a). Neighbors detected within the neutral zone impart no interaction force. We maintain the hierarchy that if nrep > 0, then ~fi,att = 0. Inclusion of a neutral zone does not affect neighbor distributions (Fig. 5.10.b- c), as relative location to neighbors still arises from a balance between short- range repulsion, and cohesion imparted by attraction; the only difference is the distance at which neighbors are detected to generate cohesion. Model III is therefore inadequate. Model IV: Repulsion, Alignment, Attraction We add an alignment force force to ~fi,int, so that ~fi,int = ωrep ~fi,rep+ωatt ~fi,att+ ωal ~fi,al, where ωal is a weighting parameter and ~fi,al = 1 nal nal∑ j=1 ~vj |~vj | , j ∈ AL. (5.5) 135 5.4. Building a Model AL is an alignment zone, between zones of attraction and repulsion (i.e., having replaced the neutral zone of Model III), with width ral. nal is the number of detected neighbors in the alignment zone AL. The alignment force acts differently than attraction and repulsion in that is a propulsive force which increases the speed of an individual. To eliminate the artefact of large speed changes depending on whether nal is non-zero or not, we set ~fi,al = ~vi |~vi| if nal = 0, following [11]. Inclusion of alignment imparts a bias in neighbor location, as fewer neigh- bors are found directly to the side (Fig. 5.11.b). This occurs due to an in- terplay between alignment and noise-induced heading changes. However, the bias is gradual (Fig. 5.11.c), and does not match the frontal preference shown by the data, which has more distinct regions of higher and lower neighbor density. Exploration of parameter space did not reveal any parameter com- bination that matched the angular preference in data well, and so Model IV is deemed inadequate. Model V : Repulsion, Alignment, Attraction, Frontal Interaction In order to obtain a frontal preference for neighbors as in empirical observa- tions (Fig. 5.3-5.4), a further modification is proposed. A fourth force ~fi,front is added to the interaction force, so that ~fi,int = ωrep ~fi,rep + ωatt ~fi,att + ωal ~fi,al + ωfront ~fi,front, where ωfront is a weighting parameter. ~fi,front is an attraction/repulsion interaction (with strength as in Fig. 5.7) with a single neighbor ~xj,θ in a frontal angular region with angle θ; i.e., ~fi,front = gf(|~xj,θ − ~xi|) ~xj,θ − ~xi|~xj,θ − ~xi| , where gf(x) is given in Fig. 5.12. If no such ~xj,θ is detected, ~fi,front = 0. This force differs from other forces in that it is ‘topological’ (following the terminology of [3]) : should a neighbor exist in the frontal θ region (and within sensing range ratt), an interaction takes place with the first encoun- tered neighbor, regardless of distance. In this way, individuals within the 136 5.4. Building a Model Rep. (a) distance, in BL di st an ce , i n BL Relative Position of Neighbours −5 0 5−6 −4 −2 0 2 4 6 0 0.2 0.4 0.6 0.8 1 (b) 0 50 100 1500 0.2 0.4 0.6 0.8 1 Angular neighbor density at pref. dist. relative angle, in degrees de ns ity (c) Figure 5.8: Model I: (a) A schematic diagram representing the interaction zones, here being simply repulsion. (b) Spatial distributions of neighbors, relative to a focal individual oriented in the direction [0 1] (as indicated by the white graphic at the origin), where density has been normalized to have maximum value equal to 1. White dashed lines are superimposed at radial distances 1, 2 and 3 BL. (c) Angular density distributions at the preferred distance are plotted, following the convention of Fig. 5.4. I compare distri- butions for data (green) and the model (blue). 137 5.4. Building a Model Rep. Att. (a) distance, in BL di st an ce , i n BL Relative Position of Neighbours −5 0 5−6 −4 −2 0 2 4 6 0 0.2 0.4 0.6 0.8 1 (b) 0 50 100 1500 0.2 0.4 0.6 0.8 1 Angular neighbor density at pref. dist. relative angle, in degrees de ns ity (c) Figure 5.9: As in Fig. 5.8, but for Model II. 138 5.4. Building a Model Rep. Att.N. (a) distance, in BL di st an ce , i n BL Relative Position of Neighbours −5 0 5−6 −4 −2 0 2 4 6 0.2 0.4 0.6 0.8 1 (b) 0 50 100 1500 0.2 0.4 0.6 0.8 1 Angular neighbor density at pref. dist. relative angle, in degrees de ns ity (c) Figure 5.10: As in Fig. 5.8, but for Model III. 139 5.4. Building a Model Rep. Att.Al. (a) distance, in BL di st an ce , i n BL Relative Position of Neighbours −5 0 5−6 −4 −2 0 2 4 6 0 0.2 0.4 0.6 0.8 1 (b) 0 50 100 1500 0.2 0.4 0.6 0.8 1 Angular neighbor density at pref. dist. relative angle, in degrees de ns ity (c) Figure 5.11: As in Fig. 5.8, but for Model IV. 140 5.5. Parameters x g f (x) 1 -1 rrep ratt Figure 5.12: Attraction-repulsion function gf(x) is negative in repulsion re- gion R, and positive and constant beyond R (up to a limit of ratt). The magnitude of gf(x) is controlled by the weighting parameter ωfront. group balance forces of repulsion, alignment, and attraction in all directions with a follow-the-leader type of force with a frontal individual. Inclusion of a frontal interaction force, in combination with attraction, repulsion, and alignment, generates spatial distributions and angular bias similar to those generated from empirical data, and so we deem Model V to be adequate in capturing the spatial distribution of neighbors. Distributions are dependent on the relative contribution of forces, which we investigate via weighting parameters in the next section. Our goal is twofold; to understand the effect of parameter variation on neighbor distributions, and to obtain an optimal parameter set that best matches observed data. 5.5 Parameters Parameters in the full model are categorized into two groups: those that are fixed (either determined from data, or having negligible effect on model output), and those that are free parameters, fitted to the model. A summary of parameters is given in Table 5.2. 5.5.1 Fixed Parameters Interaction Zones Interaction zones of repulsion, alignment, attraction are controlled by radii rrep, ral, and ratt, respectively. We fix rrep = 1.45, corresponding to the mean 141 5.5. Parameters Table 5.1: Model : summary of interaction forces force hierarchy formulation ~fi,rep if nrep > 0 ~fi,rep = 1 nrep nrep∑ j=1 g(|~xj − ~xi|) ~xj − ~xi|~xj − ~xi| , j ∈ R ~fi,al if nrep = 0, nal > 0 ~fi,al = 1 nal nal∑ j=1 ~vj |~vj| , j ∈ AL, ~fi,al = ~vi |~vi| if nal = 0 ~fi,att if nrep = 0, natt > 0 ~fi,att = 1 natt natt∑ j=1 g(|~xj − ~xi|) ~xj − ~xi|~xj − ~xi| , j ∈ ATT ~fi,front if ~xj,θ exists ~fi,front = g(|~xj,θ − ~xi|) ~xj,θ − ~xi|~xj,θ − ~xi| nearest-neighbor distance in the data. In Fig. 5.5, deviation as a function of radius is minimal at the preferred distance of 1.45 BL, and then increases with radius. We estimate ral = 3, twice the distance of minimal deviation, which corresponds to a threshold mean deviation value of approximately 11◦. Although alignment is clearly not radially symmetric in Fig. 5.5, it is beyond the scope of this paper to match alignment distributions, and so we assume that the alignment zone is circular. We choose ratt = 5 as an estimation of interaction radius, though model results are not sensitive to this width, due to the upper limit on interacting neighbors. The angular width of the frontal interaction zone is chosen to be θ = 60◦, as the angular region of preferred neighbors in the data extends approximately 30◦ to the left and right of the heading of an individual (as in Fig. 5.4). Autonomous Forces Autonomous propulsion ~a is chosen to be [0 a] (where a is some positive parameter); direction choice is arbitrary, having no effect on spatial distri- butions. Because we focus on the relative weight of the various interaction forces, we fix a = 0.5 as a reference parameter. Equilibrium speed in Eq. (5.2) can be solved as v0 = |~a|+ ωal γ , (5.6) 142 5.5. Parameters fixed parameters parameter value source rrep 1.45 BL data ral 3 BL data (est.) ratt 5 BL chosen θ 60◦ data a 0.5 BLs−2 chosen (reference) γ (ωal + a)/2 s −1 data (via Eq. 5.6) free parameters parameter range optimal ωrep 1-20 10 ωal 0-2 0.5 ωatt 0.1-3 0.5 ωfront 0-1 0.1 ωξ 0.05-0.5 0.325 Table 5.2: Summary of parameters: Using Model V together with an opti- mization routine, an optimal parameter set is found that best matches simu- lated neighbor distributions to observed neighbor distributions. Parameters are dimensionless unless units are given. and so the drag coefficient γ is determined by ~a and ωal by matching equi- librium speed to mean speed calculated from empirical data (2.0 BL/s). 5.5.2 Free Parameters In order to explore the effects of relative weighting of the various interaction forces, we leave ωrep, ωatt, ωal and ωfront as free parameters. We first explore changes in distributions resulting from parameter variation, and then use an optimization algorithm to find the set of parameters which provide the best fit to the observed data. We also explore variations in the relative magnitude of noise in the system, and so let ωξ vary. 5.5.3 Parameter Effects on Radial Distributions Matching simulated distributions to observed distributions involves accu- rately representing the radial distribution of neighbors in the data, as well as 143 5.5. Parameters 0 2 4 60 0.2 0.4 0.6 0.8 1 radius, in BL de ns ity radial neighbor density (data) Figure 5.13: Average radial neighbor density, from data. the angular distribution. To first investigate radial distribution dependence on parameters, we set ωfront = 0, and set the remaining free parameters to typical values (see caption of Fig. 5.14). Then, free parameters are varied one by one, keeping the remaining parameters fixed at the typical values. In Figure 5.14, radial neighbor density is plotted over ranges of the free param- eters ωrep, ωatt, ωal and ωξ. Radial neighbor density from the data is shown for reference in Fig. 5.13. In Fig. 5.14.a, increasing ωrep from 0 lowers density within R, and en- hances periodicity in density due to well-spaced group structure. In Fig. 5.14.b, varying ωatt has little effect because at (quasi-)steady state, attraction forces occur in all directions, and tend to cancel. Self-propulsion is governed by the sum of ~a and ~fi,al; increasing ωal (Fig. 5.14.c) weights propulsion via alignment more strongly than autonomous forcing from ~a, which reduces group cohesion imparted by ~a (individuals sharing a common preferred di- rection will naturally move cohesively). Less cohesiveness in turn increases neighbor density at distances beyond preferred spacing rrep. Though it may seem counterintuitive that increasing alignment reduces cohesion, it is due to the tradeoff experienced by reducing the relative strength of ~a. It is only at considerably larger alignment radii ral (≈ 5 BL) that alignment forcing alone (i.e., with ~a = 0) leads to cohesive groups. In Fig. 5.14.d, increasing noise reduces well-spaced group structure, leading to decreased periodic structure in neighbor density. 144 5.5. Parameters 0 2 4 60 0.2 0.4 0.6 0.8 1 ω rep = 10 ω rep = 0 radius de ns ity (a) Effect of varying ωrep = 0, 1, 3, 5, 10. 0 2 4 60 0.2 0.4 0.6 0.8 1 ω att = 10 ω att = 0 radius de ns ity (b) Effect of varying ωatt = 0, 1, 5, 10. 0 2 4 60 0.2 0.4 0.6 0.8 1 ω al = 0 ω al = 10 radius de ns ity (c) Effect of varying ωal = 0, 1, 3, 5, 10. 0 2 4 60 0.2 0.4 0.6 0.8 1 ωξ = 0 ωξ = 0.4 radius de ns ity (d) Effect of varying ωξ = 0, 0.1, 0.2, 0.3, 0.4. Figure 5.14: Effect of parameter variation on radial neighbor distribution in basic model simulations. Unless varied as indicated, parameters are ωrep = 5, ωatt = 1, ωal = 1, ωξ = 0.15. 50 Simulations of 100 individuals were performed to t = 100, and average densities over all simulations (after quasi-equilibrium was reached, t = 50 to t = 100) were calculated. 145 5.6. An Optimal Parameter Set 5.6 An Optimal Parameter Set 5.6.1 Goodness-of-fit Measure The modeling goal is to recover both the observed spatial distribution of near- est neighbors, with particular emphasis on observed angular bias in neigh- bor location. To do this, we first develop a simple goodness-of-fit measure. Ideally, if the model was able to completely recover Fig. 5.3, the error be- tween the model-generated and data-generated spatial distribution of neigh- bors would be sufficient. However, real groups are much more complex than those generated by a model, and so Fig. 5.3 contains influences not included in the model developed here. As a result, to emphasize matching the angu- lar bias in neighbors, the goodness-of-fit measure used is a combination of overall spatial distribution fit, and matching angular variation in neighbor density in the radial region of preferred distance. The measure is defined as E = 〈|ρdata(x, y)− ρsim(x, y)|〉+ 〈|hdata(θ)− hsim(θ)|〉, (5.7) where ρdata (Fig. 5.3) and ρsim are the observed and simulated two-dimensional spatial distributions of neighbors, respectively, and hdata (Fig. 5.4) and hsim are the observed and simulated angular variation in neighbors. 〈·〉 denotes an average (two-dimensional and one-dimensional, respectively, in Eq. (5.7)). 5.6.2 Optimization Process Parameter space was first extensively investigated manually. From visual inspection, a reasonable set of parameters was chosen to initialize the free parameters listed in Table 5.2 in the optimization routine. Lower and upper bounds are placed on each parameter, given in Table 5.2, obtained manually from inspection of model output, requiring that simulations appear realistic. Then, a series of 30 simulations to t = 200 are performed (randomized initial conditions), with data from the last half of the simulation used to construct fsim and hsim. These distributions are used to calculate E via Eq. (5.7). A pattern search algorithm (suited for optimization of a stochastic function) is then used to update the parameter set that, after another series of sim- ulations, attempts to reduce E. This process is continued iteratively, until the parameter set that leads to the minimum value of E is found, giving the optimal set. The entire optimization process was repeated with different initial conditions to verify the computed minimum. 146 5.6. An Optimal Parameter Set 5.6.3 Optimization Results The optimization found a set of parameters that minimized the error E was found. The optimal parameter set is listed in Table 5.2. Of all influences, repulsion is weighted most strongly ωrep = 10, an order of magnitude larger than attraction (ωatt = 0.5). Interestingly, an optimal value of ωfront = 0.1 (two orders of magnitude less than repulsion) suggests that the frontal interaction tendency need not be very strong to have a clear effect on the angular preference of neighbors. Alignment strength is optimal at ωal = 0.5, such that individuals are propelled via a = 0.5 equally as much as the propulsion component due to aligning with neighbors. Spatial distribution of neighbors for the opti- mal set is shown in Fig. 5.15.a together with the corresponding distribution from data (repeated from Fig. 5.3) in Fig. 5.15.b. The essential features of the spatial distribution are captured in the simulated data, although the secondary shell of neighbors at radial distance 3 BL is more evident in sim- ulated data, likely due to irregularities in the shape and structure of real groups that cannot be accounted for by white noise in simulations. In Fig. 5.16, angular distribution of neighbors at the preferred distance is plotted for both the data and simulations with the optimal parameter set. For simula- tions upon which Figs. 5.15.a and 5.16 are based, the mean absolute error in spatial distributions is 〈|fdata(x, y)− fsim(x, y)|〉 = 0.0388, and the mean absolute error in angular distribution is 〈|hdata(θ)− hsim(θ)|〉 = 0.0201, for a total error E = 0.0589. 5.6.4 Detailed Parameter Exploration Near Optimal Set To verify that the optimal parameter set minimizes error E (locally, at least), we vary each parameter in turn in a local region about the optimal value, keeping all others fixed at the respective optimal value. For each variation, 30 simulations are run three times, and an average of the three values of E is computed, giving the error measure for a particular variation. Fig. 5.17 shows the results for varying each of the 5 free parameters about the optimal set (where parameter values are all scaled to have coinciding optimal values). It is clear that for each parameter, the optimal value minimizes error. We note that in Fig. 5.17, error values at the optimal parameter line do not coincide for each exploration due to the stochastic nature of the 147 5.6. An Optimal Parameter Set (a) Simulations (b) Data Figure 5.15: A comparison of spatial neighbor density for [a] the best-fit simulations of the most appropriate model, Model V , and [b] data (repeating Fig. 5.3). Essential features of the data are observed in simulated data, including the spatial extent of neighbors, and the frontal bias. 148 5.6. An Optimal Parameter Set -90 00 0.2 0.4 0.6 0.8 1 angular neighbor density at pref. dist. relative angle, in degrees de ns ity 90 Figure 5.16: A comparison of average angular neighbor density at the pre- ferred distance for data (green), and simulation of Model V with optimal parameters (blue). 149 5.7. Discussion 0.05 0.06 0.07 0.08 0.09 0.1 opt. par. e rr o r Parameter variation around optimal set ωξ ω att ω rep ωfront ω al Figure 5.17: Error as each parameter is varied about the optimal value. Each parameter is scaled in the plot so that optimal values coincide (blue dotted line). Parameter values explored were as follows: ωrep = 5, 7.5, 10, 12.5, 15, ωatt = 0.5, 1, 1.5, 2, ωal = 0.25, 0.5, 0.75, 1, ωfront = 0.05, 0.1, 0.15, 0.2, and ωξ = 0.275, 0.3, 0.325, 0.35, 0.375. In each case, the optimal value is a mini- mum. system (though parameters are equal at this line). 5.7 Discussion Recent advances in camera technology have made empirical studies on large animal groups much more feasible, although many challenges still exist. Much of the previous work has focused on flocks of birds flying [2, 4, 14, 18], or fish schooling [7, 8, 16, 17, 22], both of which occur naturally in three dimensions. In this study a group of surf scoters foraging collectively was analyzed because their motion of interest was limited to two dimensions (the water surface). Such groups eliminate many of the challenges associated with empirical data collection of collective motion, yet still exhibit relatively complex responses and interactions. Local interactions of individuals were investigated via spatial distribu- 150 5.7. Discussion tions of alignment and relative positioning, indicating a frontal bias in neigh- bor position, and a distinct avoidance manoeuvre in the frontal direction. The frontal bias of neighbors suggests that scoters have a follow-the-leader tendency, within a two-dimensional group structure. Interestingly, the scoters studied were observed on numerous occasions to form distinct follow-the-leader formations while approaching the dock, (see Fig. 5.18). It thus seems likely that scoters are able to modify group shape by changing interactions, although frontal tendencies that entirely govern one-dimensional formations are still present to a degree in two-dimensional formations. Surf scoters observed in this study were foraging on mussels growing on dock pilings, and so one possible explanation for the various group formations is to optimize foraging depending on the spatial proper- ties of the food source (whether it is widely dispersed, or concentrated). It has been suggested that scoters may cooperatively facilitate foraging of mus- sels because a removed mussel from the bed in turn facilitates removal of neighboring mussels due to increased purchase by the scoter. The spatial structure of deviation reinforces that interactions to the front are the dominant influence. Should an individual scoter move too close to a neighbor in front, a strong deviation is taken to avoid the neighbor. In contrast, individuals who are too close to neighbors at the side (though in- frequent, as shown in Fig. 5.3) do not show as strong an avoidance tendency via deviation. This mode of interaction is likely influenced by both the environment (swimming in water), and physiology of scoters. Turning ma- noeuvres are likely easier to perform than to modify speed while swimming, while avoiding many velocity changes would confer energetic benefits dur- ing approach to the dock for foraging. It is of future interest to develop a model that simultaneously matches both neighbor distribution and deviation distributions. Another feature of deviation distribution is (at distances beyond preferred spacing) an increase with distance from a focal individual. This fact is strong evidence that alignment forces exist. Consider the alternative explanation: that heading is determined by the common influence of the goal (in our study, approaching the dock). If this were the case, no clear correlation would be observed between deviation and distance, as each individual would be governed by the same global influence, and deviation differences would be caused by noise. The model proposed here, built step-by-step from general collective mo- tion principles, is an attempt to simulate the joint tendencies of existing in 151 5.7. Discussion Figure 5.18: An example of a follow-the-leader formation observed in surf scoters approaching a dock to forage on mussels. Here, interactions to the front are dominant. 152 5.7. Discussion a well-spaced cohesive group, and having particular interaction directly in front (as per a follow-the-leader formation). By modifying a dynamic, zonal attraction-repulsion-alignment model to include a topological frontal inter- action, groups were formed having similar spatial distribution of neighbors, with a distinct frontal bias near the preferred spacing. Although the model convincingly accounts for the observed neighbor distribution, it is possible that other models are able to do so as well (or better). However, it is unlikely that such a model is distinctly simpler than that proposed here, because the necessity of each component included was shown via the series of more com- plex models. Typical simulation studies of collective motion differ from this work in a fundamental way. Here, empirical data motivates both the implementation of the model, and the choice of parameters. In studies based on simulation, modeling choices are hypothetical, and parameters often arbitrarily chosen. As more empirical data becomes available, more focus should be placed on evaluating competing hypotheses for collective motion. 153 Bibliography [1] I. Aoki, A simulation study on the schooling mechanism in fish, Bull. Jpn. Soc. Sci. Fish. 48 (1982), 1081–1088. [2] M. Ballerini, N Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Gi- ardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic, Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study, Proc. Nat. Acad. Sci. 105 (2007). [3] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giar- dina, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic, An empirical study of large, naturally occurring starling flocks: a bench- mark in collective animal behaviour, arXiv:0802.1667 (2008). [4] R. Budgey, Three dimensional bird flock structture and its implications for birdstroke tolerance in aircraft, International Bird Strike Committee IBSC 24/WP 29 (1998). [5] A. Cavagna, I. Giardina, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic, Empirical investigation of starling flocks: a bench- mark study in collective animal behaviour, Anim. Behav. 76 (2008). [6] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and N. R. Franks, Collective memory and spatial sorting in animal groups, J. Theor. Biol. 218 (2002), 1–11. [7] J. Cullen, E. Shaw, and H. Baldwin, Methods for measuring the three- dimensional structure of fish schools, Can. J. Zool. 71 (1993), 1494–1499. [8] D. Grunbaum, S. Viscido, and J. Parrish, Extracting interactive control algorithms from group dynamics of schooling fish, Cooperative Control LNCIS 309 (2004), 103–117. 154 Chapter 5. Bibliography [9] S. Gueron, S. A. Levin, and D. I. Rubenstein, The dynamics of herds: From individual to aggregations, J theor Biol 182 (1996), 85–98. [10] A. Huth and C. Wissel, The simulation of movement of fish schools, J. Theor. Biol. 156 (1992), 365–385. [11] H. Levine, W.J. Rappel, and I. Cohen, Self-organization in systems of self-propelled particles, Physical Review E 63 (2001), 017101. [12] Y.-X. Li, R. Lukeman, and L. Edelstein-Keshet, Minimal mechanisms for school formation in self-propelled particles, Phys. D. 237(5) (2008), 699–720. [13] R. Lukeman, Y.-X. Li, and L. Edelstein-Keshet, A mathematical model for milling formations in biological aggregates, Bull. Math. Biol. 71(2) (2009), 352–382. [14] P. F. Major and L. M. Dill, The three-dimensional structure of airborne bird flocks, Behavioural Ecology and Sociobiology 4 (1978), 111–122. [15] H.-S. Niwa, Self-organizing dynamic model of fish schooling, J theor Biol 171 (1994), 123–136. [16] B. L. Partridge, Fish school density and volume, Marine Biology 54 (1979), 383–394. [17] B. L. Partridge, T. Pitcher, J. M. Cullen, and J. Wilson, The three- dimensional structure of fish schools, Behav Ecol Sociobiol 6 (1980), 277–288. [18] H. Pomeroy and F. Heppner, Structure of turning in airborne rock dove flocks, Auk 109(2) (1992), 256–267. [19] S. Sakai, A model for group structure and its behavior, Biophysics Japan 13 (1973), 82–90. [20] J. Tien, S. Levin, and D. Rubenstein, Dynamics of fish shoals: identify- ing key decision rules, Evolutionary Ecology Research 6 (2004), 555–565. [21] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Phys. Rev. Lett. 75(6) (1995), 1226–1229. 155 Chapter 5. Bibliography [22] S. Viscido, J. Parrish, and D. Grunbaum, Individual behaviour and emer- gent properties of fish schools: a comparison of observation and theory, Mar. Ecol. Prog. Ser. 273 (2004), 239–249. [23] K. Warburton and J. Lazarus, Tendency-distance models of social cohe- sion in animal groups, J. theor. Biol. 150 (1991), 473–488. 156 Chapter 6 Conclusion 6.1 Introduction In this thesis, I explored collective motion in animal groups with two ap- proaches. First, perfect schools were studied with a minimal model, and existence and stability conditions for particular formations were derived. Sec- ond, real data on motions of surf scoters swimming in polarized groups were gathered and used to test a series model hypotheses, leading to a specific model that matched observations of real scoter flocks. 6.2 Analysis of Perfect School Solutions 6.2.1 Summary In Chapters 2 and 3, a minimal model for school formation was presented. Though the model exhibited a wide array of possible solutions, I focused on a number of solutions that featured regular spacing and group geometry, and a predictable interaction structure, so-called ‘perfect schools.’ The regular- ity of these perfect school solutions was exploited to obtain conditions on model parameters that guarantee existence of solutions, and conditions that determine linear stability. Where a number of equilibria were possible, sta- bility conditions could determine the stability of each solution. For realistic interactions, a unique stable solution was shown to exist. The geometric regularity of these schools appeared in the model’s equa- tion system as a structured coefficient matrix. This structure was then used to obtain a set eigenvalues for arbitrary dimension (i.e., n individuals). In chapter 2, the matrix structure was block- and block-diagonal, and in Chap- ter 3, block-circulant. Sets of eigenvalues led to conditions on model inputs that determine stability. In Chapter 2, the stability condition was expressed in terms of the attraction- repulsion function g(x) governing interactions between individuals. I showed 157 6.2. Analysis of Perfect School Solutions that if the derivative of g(x) was positive at a given equilibrium school solu- tion, then that solution was stable. This was true for perfect schools in one dimension, and for ‘soldier formations’ in two dimensions. In Chapter 3, these techniques were used to explore another perfect school: individuals in a milling formation, following one another in a cir- cular path. Importantly, the change in group geometry changed conditions on g(x) leading to stability. Whereas explicit solutions for eigenvalues were found for perfect schools in Chapter 2, the milling formation stability analysis led to a characteristic equation in the form of a complex-coefficient quartic polynomial. This equation was studied numerically which led to an interval of values of the derivative of g(x) at steady-state which lead to stability. In contrast with results in Chapter 2, an upper limit existed on g′(x) to guarantee stability, while also, negative values of this derivative led to stable formations under certain parameters. Although the solutions studied in Chapters 2 and 3 were idealized, they led to clear relationships between model parameters and observed group be- haviour. Previously, models were typically investigated by observing model output in simulations while varying parameters, and then developing a hy- pothesis to explain the observed output. In contrast, the conclusions in Chapters 2 and 3 state explicitly the properties of the group in terms of parameters, eliminating the need for detailed parameter exploration (in the case of the solutions we study). 6.2.2 Perfect School Analysis in Context The types of solutions studied in this thesis (soldier formations and mill for- mations) are representative of two solutions that occur across many different models for collective motion [2, 3, 4, 6, 8, 13, 14]: polarized motion, and cir- cular (milling) motion. Although the solutions obtained in these models are more complex, and are truly two-dimensional, our results show that besides model parameters, the group geometry itself is an important component in determining stable patterns observed in a model. There have been a number of other attempts to obtain analytical pre- dictions of model behavior. Analogy to statistical mechanics [2, 4] led to conditions on whether groups remained well-spaced (H-stable), or became ever more dense (catastrophic), as number of individuals increased. Related results on existence of well-spaced groups were obtained using Lyapunov sta- bility [7] and assuming particular interaction forces. The analytical work in 158 6.3. Empirical Data and Modeling Chapters 2 and 3 is more general in that no particular form of interaction function is assumed, and stability results are to general perturbations, as op- posed to the particular case of group collapse (i.e. H-stable vs. catastrophic behavior). However, there are a number of limitations to this work apart from the requirement of regular geometry. In order to obtain results on stability for groups of arbitrary size, a well-defined structure of interindividual interac- tion must exist, as interactions lead to coupling of model equations, which in turn leads to structure in the coefficient matrix. Without a clear con- nection structure, these methods become unwieldy. A particular example is in extending the analysis to regular arrays of individuals, as in a lattice structure. It is a significant challenge to label individuals, and determine a connection structure, that leads to a structured matrix which can be solved. The formations we studied shared the property of being organized relative to one another in a quasi-one-dimensional manner, thus avoiding the connection and labeling challenges. Furthermore, these methods do not always general- ize to other Lagrangian models. The models of [4, 6] for example, contain a non-linear drag term which complicates stability analysis of mill forma- tions substantially. However, in spite of the challenges, there remains the potential to develop these techniques further to overcome these limitations. Additionally, analytical results obtained in this thesis have direct implication for designing stable multi-agent systems, especially in the case where motion is not subject to any constraints (such as underwater vehicles). 6.3 Empirical Data and Modeling 6.3.1 Summary In Chapters 4 and 5, I described research collecting and analyzing data on groups of surf scoters. Chapter 4 described the empirical methods used to obtain and process data, and also reported on basic statistics measured from the data. The experimental subject, flocks of scoters swimming collectively, did not possess the typical difficulties associated with gathering data on other animal groups. This was mostly due to their motion in two dimensions, but also to their predictable location and relatively slow speeds of movement. An array of techniques were developed to process data, including image detec- tion software and routines for correcting camera angle, tracking, and current 159 6.3. Empirical Data and Modeling effects. The main result of Chapter 4 was in obtaining the final data set. This data set represents a significant step forward in data collection on animal groups, being the first empirical study (to my knowledge) to capture trajectories of each individual within large groups (a few hundred individuals). Previous work has been limited to small groups (an order of magnitude less than those considered in this thesis) [5, 9, 10, 11, 12, 15], or studies that abandon tracking altogether, instead gathering static images of groups without any description of movement [1]. In Chapter 5, the data set from Chapter 4 was used to build a model to match properties of real groups. Spatial neighbor distributions were con- structed to deduce the spatial preference for neighbors in local space around an individual, and to relate deviation response to relative location. It was shown that surf scoters diligently maintain empty space around themselves, and also have a bias to positioning themselves directly behind a neighbor, in a type of follow-the-leader formation. Surf scoters also exhibited a strong tendency to deviate sideways, should they approach a frontal neighbor too closely, whereas the avoidance tendency was much weaker with neighbors at the side. A positive correlation between distance from a neighbor and heading deviation gave direct evidence for alignment forces existing locally between individuals. The spatial distribution of neighbors was then used to develop and test a model for collective motion. A model was built from generic collective motion principles of attraction, repulsion, and alignment. Showing that these forces were not sufficient to account for observations of relative neighbor locations, an additional topological interaction to the front was proposed. By combining this frontal interaction force with the other typical interaction forces, essential elements of the neighbor distribution were captured, namely the frontal preference. An error measure was devised to measure goodness- of-fit, and this measure was used to obtain an optimal set of parameters (representing relative weights of the constituent forces) that best described the empirical data. 160 6.3. Empirical Data and Modeling 6.3.2 Empirical Data and Modeling in Context Empirical Data: Strengths and Limitations We have emphasized throughout Chapters 4 and 5 that the empirical data gathered is, in a sense, the first of its kind: large scale trajectory data of indi- viduals within groups of a few hundred individuals. There have been studies of much larger groups (over 1000 starlings [1]), but the inherent complex- ity and three-dimensional structure prevented reconstructing all positions, which is a crucial requirement for generating trajectories. However, there are a number of limitations, both in our data set, and in applying these experimental methods to other groups. Foremost on the list of limitations is the temporal resolution of the data, gathered at 3 frames per second. At this resolution, tracking is much more of a challenge than at higher resolutions, especially in analyzing non-polarized groups (for example, comparison to properties of polarized groups). Also movement at scales finer than our resolution is, of course, lost. This leads to rather noisy descriptions of sensitive measurements such as acceleration. Another current limitation on the methods is in the need for manual inter- vention during image processing. Because there are situations that cannot currently be automatically detected (such as the separation of two individu- als who are very close together, appearing as one), one has to inspect each image and fix any detection errors. Continued refining of the software can alleviate much of this manual effort, however. The methods here relied on an overhead location from which to take pictures, and also calibrate the experimental setup (i.e., using the dock width and height above water level). To extend these methods to other groups not near convenient overhead locations, the camera could for example be mounted on a retractable pole that could be placed in a location of interest. I am currently working with Dr. J. Heath on extending these techniques to eider duck groups in the Canadian Arctic. Modeling of Real Groups: Strengths and Limitations In Chapter 5, the distribution of neighbors and deviation in a local neigh- borhood of an individual was used to uncover interaction structure between individuals. These statistics did indeed show clear spatial patterns from which interactions were inferred. In particular, the spatial distribution of deviation has never been shown for any group of animals, owing to the joint 161 6.4. Summary: Main Contributions need for large spatial extent and trajectory reconstruction. Empirical data was used to refute or accept a series of models to describe spatial structure in scoter flocks. The final model was successful in matching the particular quantity chosen for comparing simulation results to data, and it was clearly shown that the preceding sequence of null models were not sufficient to describe the data. However, other quantities, such as the distribution of deviation, were not accurately captured in the model. Additional complexity is likely needed to capture more aspects of the data. The presentation of a feasible model for frontal preference also does not preclude the possibility that another competing hypothesis could account for the data as well, or better than the model proposed. However, it is a strong first step in linking data with dynamic models of collective motion. Empirical Data and Modeling: Future Work There remains much work to be done in interpreting the data set of Chapter 4. In particular, attraction and repulsion responses can be measured, so that attraction-repulsion curves can be empirically derived, as opposed to being formulated in terms of hypothetical functional forms (exponential [4, 6], inverse-power [7] etc.). Also, heading correlations to particular neighbors can be calculated to attempt to deduce with which neighbors an individual interacts. Another interesting aspect of the raw data is the sequences in which individuals do not move in polarized groups, but instead, move in disordered directions, maintaining a stationary group location (such groups were not analyzed in this thesis). As these transitions from disorder to order have been studied extensively in models (e.g., [3, 8, 14]), analyzing real occurrences of such transitions would be very informative. 6.4 Summary: Main Contributions My main contributions in this thesis were: 1. Obtaining existence and stability conditions on one-dimensional and ‘soldier formation’ solutions in terms of an interaction force and its derivative; 162 6.4. Summary: Main Contributions 2. Obtaining existence and stability conditions on ‘milling formation’ so- lutions in terms of an interaction force and its derivative; 3. Development of field methods to obtain a data set of hundreds of indi- vidual surf scoters moving collectively, where each individual’s trajec- tory has been reconstructed; 4. Analysis of the field data set via spatial distributions; 5. Development of a model, built from a series of null models of increasing complexity, that matches empirical spatial distributions well. 163 Bibliography [1] M. Ballerini, N Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Gi- ardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic, Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study, Proc. Nat. Acad. Sci. 105 (2007). [2] Y. L. Chuang, M. R. D’Orsogna, D. Marthaler, A. L. Bertozzi, and L. S. Chayes, State transitions and the continuum limit for a 2d interacting, self-propelled particle system, Physica D 232 (2007), 33–47. [3] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and N. R. Franks, Collective memory and spatial sorting in animal groups, J. Theor. Biol. 218 (2002), 1–11. [4] M. R. D’Orsogna, Y. L. Chuang, A. L. Bertozzi, and L. S. Chayes, Self- propelled particles with soft-core interactions: Patterns, stability, and collapse, Physical Review Letters 96 (2006), 104302. [5] D. Grunbaum, S. Viscido, and J. Parrish, Extracting interactive control algorithms from group dynamics of schooling fish, Cooperative Control LNCIS 309 (2004), 103–117. [6] H. Levine, W.J. Rappel, and I. Cohen, Self-organization in systems of self-propelled particles, Physical Review E 63 (2001), 017101. [7] A. Mogilner, L. Edelstein-Keshet, L. Bent, and A. Spiros, Mutual in- teractions, potentials, and individual distance in a social aggregation, J Math Biol 47 (2003), 353–389. [8] H.-S. Niwa, Self-organizing dynamic model of fish schooling, J theor Biol 171 (1994), 123–136. [9] B. L. Partridge, Fish school density and volume, Marine Biology 54 (1979), 383–394. 164 Chapter 6. Bibliography [10] , The effect of school size on the structure and dynamics of min- now schools, Animal behaviour 28 (1980), 68–77. [11] B. L. Partridge, T. Pitcher, J. M. Cullen, and J. Wilson, The three- dimensional structure of fish schools, Behav Ecol Sociobiol 6 (1980), 277–288. [12] H. Pomeroy and F. Heppner, Structure of turning in airborne rock dove flocks, Auk 109(2) (1992), 256–267. [13] S. Sakai, A model for group structure and its behavior, Biophysics Japan 13 (1973), 82–90. [14] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Phys. Rev. Lett. 75(6) (1995), 1226–1229. [15] S. Viscido, J. Parrish, and D. Grunbaum, Individual behaviour and emer- gent properties of fish schools: a comparison of observation and theory, Mar. Ecol. Prog. Ser. 273 (2004), 239–249. 165 Appendix A A.1 The Stability Analysis in 1D We write (2.44)-(2.47) in matrix form as d dt δ1 δ2 ... δn ω1 ω2 ... ωn = | 0 | I | −− −− −− | −− −− −− | V | U | δ1 δ2 ... δn ω1 ω2 ... ωn , (A.1) where all the sub-matrices are n× n. In particular, 0 is the zero matrix and I is the identity matrix. V and U are given by V = −g′−(d) g′−(d) 0 0 · · · 0 0 g′+(d) −σ g′−(d) 0 · · · 0 0 0 g′+(d) −σ g′−(d) · · · 0 0 · · · · · · · · · · · · · · · · · · · · · 0 0 0 0 · · · −σ g′−(d) 0 0 0 0 · · · g′+(d) −g′+(d) , U = −h′−(0)− γ h′−(0) 0 0 · · · 0 0 h′+(0) −s h′−(0) 0 · · · 0 0 0 h′+(0) −s h′−(0) · · · 0 0 · · · · · · · · · · · · · · · · · · · · · 0 0 0 0 · · · −s h′−(0) 0 0 0 0 · · · h′+(0) −h′+(0)− γ , where s = h′+(0)+ h ′ −(0)+ γ and σ = g ′ +(d)+ g ′ −(d). Note that U and V are tridiagonal matrices. 166 A.1. The Stability Analysis in 1D The eigenvalues of J are obtained by solving Λ(λ) ≡ det[λI − J ] = det [ λI −I −V λI − U ] = 0. (A.2) We can simplify this eigenvalue equation via the result stated in Eq. (3.30). In this case, because λI and −I commute, then det[λI − J ] = det[λI − UλI − IV ] = det[λ2I − λU − V ]. (A.3) Since U and V are tridiagonal matrices, so is the matrix in (A.3). We can now write (A.2) explicitly as det χ+ η α 0 · · · 0 0 η χ α · · · 0 0 · · · · · · · · · · · · · · · · · · 0 0 0 · · · χ α 0 0 0 · · · η χ+ α = 0, (A.4) where χ = λ2 + sλ+ σ, η = −h′+(0)λ− g′+(d), and α = −h′−(0)λ− g′−(d). Thus, the eigenvalue equation is a polynomial of λ that is defined by the determinant of a tridiagonal matrix. Although the determinant can be cal- culated recursively, the complexity of the resulting formula prevents us from solving the eigenvalues easily. However, the eigenvalues are easily solvable in the case when either g+ or g− is absent, i.e., when individuals interact only with neighbours in front or behind them. Let us consider the case that particles see only what is in front of them. This is the model for schools formed by following one immediate predecessor. In this case, h′−(0) = g−(d) = 0. Thus, α = 0, χ = λ 2+(h′+(0)+γ)λ+g ′ +(d), and χ+ η = λ2 + γλ = λ(λ+ γ). The eigenvalue equation now requires only the determinant of a triangular matrix which is given by the simple formula Λ(λ) = (χ + η)χn−1. Thus, the eigenvalues are determined by Λ(λ) = (χ+ η)χn−1 = λ(λ+ γ)[λ2 + (h′+(0) + γ)λ+ g ′ +(d)] n−1 = 0. (A.5) Therefore, the eigenvalues are λ1 = 0, λ2 = −γ, and λ± = −(h′+(0) + γ)/2±√ [(h′+(0) + γ)/2]2 − g′+(d). Both λ+ and λ− have multiplicity n− 1. Similarly, if each particle can only detect the a particle behind but not a particle in front (a rather peculiar situation), η = 0 and χ = λ2 + (h′−(0) + γ)λ+ g′−(d). Λ(λ) = (χ + α)χn−1 = λ(λ+ γ)[λ2 + (h′−(0) + γ)λ+ g ′ −(d)] n−1 = 0. (A.6) 167 A.2. Stability Analysis for the Soldier Formation in 2D The eigenvalues are λ1 = 0, λ2 = −γ, and λ± = −(h′−(0)+γ)/2± √ [(h′−(0) + γ)/2]2 − g′−(d). Both λ+ and λ− have n− 1 multiplicity. We find that the roles of h′−(0) and g′−(d) in determining the stability are quite similar to the roles of h ′ +(0) and g′+(d) in the case of forward sensing. A.2 Stability Analysis for the Soldier Formation in 2D We can explicitly write the Jacobian matrix ∂ ~fi ∂~zi in Eq. (2.78) as J = ∂fix∂zix ∂fix∂ziy ∂fiy ∂zix ∂fiy ∂ziy = g+(|~zi|) |~zi| + g′ + (|~zi|)z2ix |~zi|2 − z2ixg+(|~zi|) |~zi|3 g′ + (|~zi|)zixziy |~zi|2 − g+(|~zi|)zixziy |~zi|3 g′ + (|~zi|)zixziy |~zi|2 − g+(|~zi|)zixziy |~zi|3 g+(|~zi|) |~zi| + g′ + (|~zi|)z2iy |~zi|2 − z2iyg+(|~zi|) |~zi|3 . (A.7) Evaluating J at steady state ~zi = ~z s i simplifies the matrix considerably, due to the judicious choice of coordinate system aligned with the axis of the school. Specifically, since ~zsi = (zix , ziy) s = (d , 0), (A.8) substituting zix = |~zi| = d, ziy = 0 into (A.7) yields Df ≡ ∂ ~fi ∂~zi |~zi=~zsi = [ g′+(d) 0 0 g+(d)/d ] . Next, we write Eqs. (2.81)–(2.84) in matrix form as 168 A.2. Stability Analysis for the Soldier Formation in 2D d dt ~δ1 ~δ2 ... ~δn ~ω1 ~ω2 ... ~ωn = | 0 | I | −− −− −− | −− −− −− | D | −γI | ~δ1 ~δ2 ... ~δn ~ω1 ~ω2 ... ~ωn , (A.9) where D = 0 0 0 . . . 0 0 Df −Df 0 . . . 0 0 0 Df −Df 0 . . . 0 ... 0 0 . . . 0 Df −Df , (A.10) where each 0 in the matrix represents a 2×2 matrix of zeros. D is a bidiagonal block matrix. We obtain stability conditions by evaluating the eigenvalues of the coefficient matrix in Eq. (A.9) and determining under what conditions these eigenvalues have negative real part. To find the eigenvalues of the coefficient matrix in Eq. (A.9), we form the determinant equation given by this coefficient matrix: det [ −λI I D −(λ+ γ)I ] = 0. Next, we use the block matrix property in Eq. (3.30). Noting that −λI and I commute, we have det [ −λI I D −(λ+ γ)I ] = det (λ(γ + λ)I −D) . Recalling the definition of D in (A.10), where Df = [ g′+(d) 0 0 g+(d)/d ] , 169 A.2. Stability Analysis for the Soldier Formation in 2D it is clear that D is a lower triangular matrix, and thus (λ(γ + λ)I −D) is also a lower triangular matrix. λ(γ + λ) appears in the first two diagonal entries (since each D entry is 2× 2), while the remaining 2(n− 1) diagonal entries contain n− 1 pairs λ(γ + λ)− g′+(d), λ(γ + λ)− g+(d)/d. The lower triangular structure allows us to immediately infer the eigenvalue equation, given by the product of diagonal terms, i.e., ( λ(γ + λ)2 ) ( λ(γ + λ)− g′+(d) )n−1 (λ(γ + λ)− g+(d)/d)n−1 = 0. We thus get solutions λ = 0 with multiplicity 2, λ = −γ with multiplicity 2, λ = −γ 2 ± √ γ2 − 4g′+(d) 2 , each with multiplicity n− 1, and λ = −γ 2 ± √ γ2 − 4g+(d)/d 2 each with multiplicity n− 1. 170 Appendix B B.1 Transformation to Relative Coordinates Expanding (3.18)–(3.19) gives x̂i = cosφi (xi+1 − xi) + sinφi (yi+1 − yi) , (B.1) ŷi = − sinφi (xi+1 − xi) + cosφi (yi+1 − yi) , (B.2) ûi = cosφi (ui+1 − ui) + sinφi (vi+1 − vi) , (B.3) v̂i = − sinφi (ui+1 − ui) + cosφi (vi+1 − vi) , (B.4) where ~xi = (xi, yi), and ~vi = (ui, vi). We can obtain the first model equation through time-differentiation of x̂i and ŷi, which gives d dt x̂i = −ω0 sinφi (xi+1 − xi) + cosφi (ẋi+1 − ẋi) +ω0 cosφi (yi+1 − yi) + sinφi (ẏi+1 − ẏi) , d dt ŷi = −ω0 cosφi (xi+1 − xi)− sinφi (ẋi+1 − ẋi) −ω0 sinφi (yi+1 − yi) + cosφi (ẏi+1 − ẏi) . Using (B.1)–(B.4), we rewrite the above in the more compact vector no- tation d dt ~̂xi = ~̂vi + ω0k× ~̂xi, (B.5) where k× ~̂xi = det i j kx̂i ŷi 0 0 0 1 = (ŷi,−x̂i). Similarly, we time-differentiate (B.3)–(B.4), giving ˙̂ui = −ω0 sin φi (ui+1 − ui) + cosφi (u̇i+1 − u̇i) +ω0 cosφi (vi+1 − vi) + sin φi (v̇i+1 − v̇i) , (B.6) ˙̂vi = −ω0 cos φi (ui+1 − ui)− sin φi (u̇i+1 − u̇i) −ω0 sin φi (vi+1 − vi) + cosφi (v̇i+1 − v̇i) . (B.7) 171 B.1. Transformation to Relative Coordinates This can be rewritten more compactly using (B.3)–(B.4) as ˙̂ ~vi = ω0k× ~̂vi +R(φi) ( ~̇vi+1 − ~̇vi ) = ω0k× ~̂vi +R(φi) ( ~f (~xi+2 − ~xi+1)− ~f (~xi+1 − ~xi) ) −R(φi)γ (~vi+1 − ~vi) . (B.8) Following (3.3), we write the interaction force terms in (B.8) in terms of unit direction vectors and magnitudes: R(φi) ( ~f (~xi+2 − ~xi+1)− ~f (~xi+1 − ~xi) ) = R(φi) (~xi+2 − ~xi+1) |~xi+2 − ~xi+1| g (|~xi+2 − ~xi+1|) −R(φi)(~xi+1 − ~xi)|~xi+1 − ~xi| g (|~xi+1 − ~xi|) . (B.9) We make a number of observations. First, the rotational matrix operator R is distance-preserving, and so |~xi+1 − ~xi| = |~̂xi|. Second, R(a)R(b) = R(a + b). We use this property to recast (B.9), using (3.18), as R(φi) (~xi+2 − ~xi+1) |~xi+2 − ~xi+1| g (|~xi+2 − ~xi+1|)−R(φi) (~xi+1 − ~xi) |~xi+1 − ~xi| g (|~xi+1 − ~xi|) = R (−2π n ) ~̂xi+1 |~̂xi+1| g ( |~̂xi+1| ) − ~̂xi|~̂xi| g ( |~̂xi| ) , where we have used R−1(a) = R(−a). We now write the equation system in relative coordinates: ˙̂ ~xi = ω0k× ~̂xi + ~̂vi, (B.10) ˙̂ ~vi = ω0k× ~̂vi +R (−2π n ) ~̂xi+1 |~̂xi+1| g ( |~̂xi+1| ) − ~̂xi|~̂xi| g ( |~̂xi| ) − γ~̂vi. (B.11) 172 B.1. Transformation to Relative Coordinates Next, we calculate the steady-state quantities. In equations (B.1)–(B.2), we replace the steady-state expressions for xi, xi+1, yi and yi+1 to obtain x̂si = − cosφi ( r0 cos ( θ̂i + 2π n ) − r0 cos θ̂i ) − sin φi ( r0 sin ( θ̂i + 2π n ) − r0 sin θ̂i ) , (B.12) ŷsi = sin φi ( r0 cos ( θ̂i + 2π n ) − r0 cos θ̂i ) − cosφi ( r0 sin ( θ̂i + 2π n ) − r0 sin θ̂i ) , (B.13) where θ̂i = 2πi/n + ω0t. These expressions can be simplified using trigono- metric difference formulas to x̂si = −r0 cos ( φi − θ̂i − 2π n ) + r0 cos ( φi − θ̂i ) , ŷsi = r0 sin ( φi − θ̂i − 2π n ) + r0 sin ( θ̂i − φi ) , which we can further simplify using the definition of φi in (3.21) to x̂si = −r0 cos ( −π 2 − π n ) + r0 cos ( π n − π 2 ) , ŷsi = r0 sin ( −π 2 − π n ) + r0 sin ( −π n + π 2 ) . With some simple trigonometric manipulation, we write x̂si = −r0 sin ( −π n ) + r0 sin ( π n ) , ŷsi = −r0 cos ( −π n ) + r0 cos ( π n ) , i.e., x̂si = 2r0 sin ( π n ) = d, (B.14) ŷsi = 0. (B.15) The steady-state positions in this coordinate frame are independent of the index i, and each position is a vector parallel to the x-axis with magnitude 173 B.2. Derivation of the Linearized Perturbed System d, the interindividual distance at steady state in the original coordinates. Instead of calculating ûsi and v̂ s i directly as we did with position, we can evaluate (B.10) at steady state, which implies that ~̂v s i = −ω0k× ~̂x s i , i.e., ûsi = 0, v̂si = ω0d, both of which are independent of index. B.2 Derivation of the Linearized Perturbed System Substituting (3.25)–(3.26) into (3.22)–(3.23) gives d dt ( ~̂x s i + ~δi(t) ) = ω0k× ( ~̂x s i + ~δi(t) ) + ~̂v s i + ~ξi(t), (B.16) d dt ( ~̂v s i + ~ξi(t) ) = ω0k× ( ~̂v s i + ~ξi(t) ) −R− (−2π n ) ~f ( ~̂x s i+1 + ~δi+1(t) ) −~f ( ~̂x s i + ~δi(t) ) − γ ( ~̂v s i + ~ξi(t) ) , (B.17) where i = 1, . . . , n, and n+ 1 is identified with 1 (henceforth assumed to be the case). Because k× (~a+~b) = k×~a+ k×~b, (B.16) can be separated into steady-state and perturbed components. Using ω0k× ~̂x s i + ~̂v s i = 0, (B.16) can be simplified to d dt ~δi = ω0k× ~δi + ~̂ξi, (B.18) 174 B.2. Derivation of the Linearized Perturbed System where we have dropped arguments of ~δi and ~ξi. We expand the nonlinear interaction terms in (B.17) via a linear approximation: R (−2π n ) ~f ( ~̂x s i+1 + ~δi+1(t) ) − ~f ( ~̂x s i + ~δi(t) ) ≈ R (−2π n )f (~̂xsi+1)+ ∂ ~f ∂~̂xi |(d,0)~δi+1 −~f ( ~̂x s i ) − ∂ ~f ∂~̂xi |(d,0)~δi. (B.19) Using this expansion in (B.17) and noting that ω0k× ( ~̂v s i ) +R (−2π n ) ~f ( ~̂x s i+1 ) − ~f ( ~̂x s i ) − γ ( ~̂v s i ) = 0, we simplify (B.17) as d dt ~ξi = ω0k× ~ξi +R (−2π n ) ∂ ~f ∂~̂xi |(d,0)~δi+1 − ∂ ~f ∂~̂xi |(d,0)~δi − γ~ξi. Combining with (B.18), we now write the set of equations for the perturbed system: d dt ~δi = ω0k× ~δi + ~̂ξi, (B.20) d dt ~ξi = ω0k× ~ξi +R (−2π n ) ∂ ~f ∂~̂xi |(d,0)~δi+1 − ∂ ~f ∂~̂xi |(d,0)~δi − γ~ξi. (B.21) We explicitly calculate the Jacobian matrix ∂ ~f ∂~̂xi by expanding ~f as in (B.11), i.e., ~f(~̂xi) = x̂i |~̂xi| g ( |~̂xi| ) i + ŷi |~̂xi| g ( |~̂xi| ) j ≡ f1ii+ f2ij, where |~̂xi| = √ x̂2i + ŷ 2 i . 175 B.3. Eigenvalue Equation We can now calculate the Jacobian matrix ∂f1i∂x̂i ∂f1i∂ŷi ∂f2i ∂x̂i ∂f2i ∂ŷi , via derivatives: ∂f1i ∂x̂i = g(|~xi|) |~xi| + g′(|~xi|)x2i |~xi|2 − x2i g(|~xi|) |~xi|3 , ∂f1i ∂ŷi = g′(|~xi|)xiyi |~xi|2 − g(|~xi|)xiyi |~xi|3 , ∂f2i ∂x̂i = g′(|~xi|)xiyi |~xi|2 − g(|~xi|)xiyi |~xi|3 , ∂f2i ∂ŷi = g(|~xi|) |~xi| + g′(|~xi|)y2i |~xi|2 − y2i g(|~xi|) |~xi|3 . Evaluating these matrix entries at the steady state values (x̂i, ŷi) = (d, 0) gives ∂f1i ∂x̂i |(d,0) = g(d) d + g′(d)d2 d2 − d 2g(d) d3 = g′(d), ∂f1i ∂ŷi |(d,0) = 0, ∂f2i ∂x̂i |(d,0) = 0, ∂f2i ∂ŷi |(d,0) = g(d) d . Then we can write the Jacobian matrix as D = ∂ ~f ∂~̂xi |(d,0) = [ g′(d) 0 0 g(d)/d ] . With definitions of Ω and RD from (3.27), we now can write the system in (3.28). B.3 Eigenvalue Equation The set of diagonal blocks is given by Di = Λ +D +̟ i−1RD, i = 1, . . . , n. We equivalently consider the set Di+1, i = 0, . . . , n − 1, for notational ease. 176 B.4. Derivation of Inequality (3.37) Figure B.1: A schematic diagram of particles at the threshold of mill break- ing. The angle between particles is π/2, beyond which particle i no longer senses particle i+ 1. ~vmi indicates the velocity of particle i in the absence of autonomous self-propulsion. Expanding detDi+1 = 0 gives p(λ) = ( λ(λ+ γ)− ω20 + g′(d)− exp ( 2πij n ) cos −2π n g′(d) ) ·( λ(λ+ γ)− ω20 + g(d)/d− exp ( 2πij n ) cos −2π n g(d) d ) + ( γω0 + 2λω0 + exp ( 2πij n ) sin −2π n g(d) d ) · ( γω0 + 2λω0 + exp ( 2πij n ) sin −2π n g′(d) ) = 0. (B.22) Using the expression for ω0 in (3.16), and for g(d)/d in (3.29), we reduce the parameters to γ, n, and g′(d), while the index i = 0, . . . , n−1 (corresponding to D1, . . . , Dn). B.4 Derivation of Inequality (3.37) At the steady state of the moving mill, we have ~vsi = ω0r0 − sin ( θ̂i(t) ) cos ( θ̂i(t) ) + ~a γ , ~ui = sin ( θ̂i(t) + π n ) cos ( θ̂i(t) + π n ) , 177 B.4. Derivation of Inequality (3.37) where θ̂i(t) = 2πi/n+ω0t. The threshold of ~a occurs when these two vectors are perpendicular, that is, ~vsi · ~ui = 0 (see Figure B.1). Expanding this condition gives − ω0r0 sin ( θ̂i(t) ) sin ( θ̂i(t) + π n ) − a1 γ sin ( θ̂i(t) + π n ) + ω0r0 cos ( θ̂i(t) ) cos ( θ̂i(t) + π n ) + a2 γ cos ( θ̂i(t) + π n ) = 0, (B.23) where ~a = [a1 a2] T . This equation simplifies using trigonometric identities to a1 γ sin ( θ̂i(t) + π n ) − a2 γ cos ( θ̂i(t) + π n ) = ω0r0 cos ( π n ) . Using the expressions for ω0 and r0 from the existence conditions, we further simplify the threshold condition to a1 sin ( θ̂i(t) + π n ) − a2 cos ( θ̂i(t) + π n ) (B.24) = γ ( γ sin(π n ) cos(π n ) )( g(d) cos2(π n ) γ2 sin(π n ) ) cos ( π n ) , = g(d) cos2 ( π n ) . (B.25) Mill breaking actually occurs whenever ~vsi · ~ui ≤ 0, i.e., a1 sin ( θ̂i(t) + π n ) − a2 cos ( θ̂i(t) + π n ) ≥ g(d) cos2 ( π n ) , so we must maximize the left-hand side over all possible angles to find the minimal values of a1 and a2 that lead to mill breaking. It can easily be shown using calculus that the left-hand side attains a global maximum value of |~a| = √ a21 + a 2 2. This gives our mill breaking condition in (3.37). 178 Appendix C C.1 Data Sequences: Details 14 sequences were gathered over 4 days in March, 2008. From all raw data gathered, sequences to be analyzed were chosen based on quality of data (including image quality, lack of external disturbance, having large groups within the frame). Due to limitations of the tracking software used, sequences longer than 50 frames are split and subcategorized by letters. Details of analyzed sequences are listed in Table C.1. We note that all groups analyzed showed a high degree of polarization, and all groups were two-dimensional in shape, except sequence S-14, wherein individuals approached the dock in a one-dimensional, follow-the-leader organization. These data were excluded from statistical measures. C.2 Data Analysis: Details on Currents For each data sequence, average alignment is manually calculated via a graph- ical user interface in MATLAB and compared with actual velocity (as de- termined by particle-tracking). Then, tracer velocity is calculated for each sequence, and together this information is used to recover individual align- ment, so as to account for effects of currents. Details are shown in Table C.2. 179 C.2. Data Analysis: Details on Currents Table C.1: Data sequences: details event date frames total (w/o edge) S-1 mar 1 25 4206 (2980) S-2 mar 6 42 2294 (903) S-3 mar 6 40 3920 (2739) S-4(a,b) mar 7 62 9257 (6277) S-5 mar 7 47 2833 (1533) S-6(a,b,c) mar 7 137 12423 (7665) S-7(a,b) mar 7 91 7285 (2571) S-8(a,b,c) mar 7 101 11023 (6667) S-9(a,b) mar 7 70 6765 (2995) S-10(a,b) mar 7 61 6628 (4057) S-11 mar 12 31 1250 (500) S-12 mar 12 22 2540 (1536) S-13(a,b) mar 12 55 3886 (2176) S-14 mar 12 44 959 (n/a) total: 828 75269 (42599) 180 C.2. Data Analysis: Details on Currents Table C.2: Current statistics: by sequence event frame b̂ v̂ cx (angle diff.) waste S-1 img. 10 [-0.88 -0.46] [-0.88 -0.46] -0.0058 (0.15o) 0 S-2 img. 38 [0.54 -0.83] [0.53 -0.84] 0.0176 (0.84o) 0 S-3 img. 3 [-0.27 0.96] [-0.28 0.96] 0.0147 (0.81o) 0 S-4 img. 25 [-0.47 -0.88 ] [-0.03 -0.99] -0.5020 (26.34o) 0.75 S-5 img. 39 [-0.25 -0.96 ] [ 0.10 -0.99] -0.3680 (20.85o) 0.74 S-6 img. 76 [0.24 0.97] [0.47 0.88] -0.2487 (13.95o) 0.66 S-6 img. 111 [-0.26 0.96] [0.07 0.99] -0.3410 (19.23o) - S-7 img. 16 [-0.17 0.98] [0.07 0.99] -0.2477 (14.13o) 0.67 S-8 img. 38 [ -0.76 -0.64] [-0.41 -0.91 ] -0.6858 (26.11o) 0.61 S-8 img. 82 [-0.08 -0.99] [0.34 -0.94] -0.4276 (25.21o) - S-9 img. 59 [-0.17 -0.98] [0.04 -0.99] -0.2190 (12.45o) 0.40 S-10 img. 2 [-0.52 -0.85] [-0.24 -0.97] -0.3588 (17.77o) 0.50 S-11 img. 21 [-0.29 -0.95] [-0.18 -0.97] -0.1213 (6.70o) 0.35 S-12 img. 18 [-0.01 0.99] [0.02 0.99] -0.0363 (2.08o) 0.15 S-13 img. 11 [0.12 -0.99] [0.19 -0.98] -0.0752 (4.28o) 0.14 S-14 img. 27 [0.53 -0.82] [0.60 -0.79] -0.0825 (3.97o) 0.44 181
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Modeling collective motion in animal groups : from...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Modeling collective motion in animal groups : from mathematical analysis to field data Lukeman, Ryan J. 2009
pdf
Page Metadata
Item Metadata
Title | Modeling collective motion in animal groups : from mathematical analysis to field data |
Creator |
Lukeman, Ryan J. |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | Animals moving together cohesively is a commonly observed phenomenon in biology, with bird flocks and fish schools as familiar examples. Mathematical models have been developed in order to understand the mechanisms that lead to such coordinated motion. The Lagrangian framework of modeling, wherein individuals within the group are modeled as point particles with position and velocity, permits construction of inter-individual interactions via `social forces' of attraction, repulsion and alignment. Although such models have been studied extensively via numerical simulation, analytical conclusions have been difficult to obtain, owing to the large size of the associated system of differential equations. In this thesis, I contribute to the modeling of collective motion in two ways. First, I develop a simplified model of motion and, by focusing on simple, regular solutions, am able to connect group properties to individual characteristics in a concrete manner via derivations of existence and stability conditions for a number of solution types. I show that existence of particular solutions depends on the attraction-repulsion function, while stability depends on the derivative of this function. Second, to establish validity and motivate construction of specific models for collective motion, actual data is required. I describe work gathering and analyzing dynamic data on group motion of surf scoters, a type of diving duck. This data represents, to our knowledge, the largest animal group size (by almost an order of magnitude) for which the trajectory of each group member is reconstructed. By constructing spatial distributions of neighbour density and mean deviation, I show that frontal neighbour preference and angular deviation are important features in such groups. I show that the observed spatial distribution of neighbors can be obtained in a model incorporating a topological frontal interaction, and I find an optimal parameter set to match simulated data to empirical data. |
Extent | 6516564 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-08-06 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0067464 |
URI | http://hdl.handle.net/2429/11873 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2009-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_fall_lukeman_ryan.pdf [ 6.21MB ]
- Metadata
- JSON: 24-1.0067464.json
- JSON-LD: 24-1.0067464-ld.json
- RDF/XML (Pretty): 24-1.0067464-rdf.xml
- RDF/JSON: 24-1.0067464-rdf.json
- Turtle: 24-1.0067464-turtle.txt
- N-Triples: 24-1.0067464-rdf-ntriples.txt
- Original Record: 24-1.0067464-source.json
- Full Text
- 24-1.0067464-fulltext.txt
- Citation
- 24-1.0067464.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0067464/manifest