Group coordination in a biologically-inspired vectorial network model

Most of the mathematical models of collective behavior describe uncertainty in individual decision making through additive uniform noise. However, recent data driven studies on animal locomotion indicate that a number of animal species may be better represented by more complex forms of noise. For example, the popular zebrafish model organism has been found to exhibit a burst-and-coast swimming style with occasional fast and large changes of direction. Based on these observations, the turn rate of this small fish has been modeled as a mean reverting stochastic process with jumps. Here, we consider a new model for collective behavior inspired by the zebrafish animal model. In the vicinity of the synchronized state and for small noise intensity, we establish a closed-form expression for the group polarization and through extensive numerical simulations we validate our findings. These results are expected to aid in the analysis of zebrafish locomotion and contribute a new set of mathematical tools to study collective behavior of networked noisy dynamical systems.


INTRODUCTION
Collective behavior is characterized by self-organized states [13,24,40] that have been observed across almost every phylum of life, including fish schools, bird flocks, insect swarms, and bacteria [13,22,37].Dissecting the determinants of collective behavior is central to elucidate the underpinnings of social behavior and inform new endeavors in engineering research on neural networks [23], communications [42], and multi-vehicle robotics [36].Collective behavior is often thought to be a manifestation of local interaction of closely connected individuals [7], yet the general underlying mechanisms are not fully understood.Animal experimentation is certainly the key approach to test hypotheses on collective behavior and aid our comprehension of this fascinating phenomenon [15,26].
In recent years, zebrafish has emerged as a popular animal model for the investigation of several behavioral and neurological disorders [20,10].In particular, this freshwater fish species has been intensively used to study the effects of psychoactive compounds on individual and social behavior [12,25,28].Studies on zebrafish behavioral response often require the use of a relatively large number of fish to cope with stringent requirements of rigorous statistical analysis [17,19].In the context of the three Rs, that is, reduce, refine, and replace animal use [38], computational models of zebrafish behavior may complement established empirical methods [29,32] through a new class of in silico experiments.
Several mathematical models have been proposed to study collective behavior, such as the classical and elegant Vicsek model [41].Within this model, each individual is described as a self-propelled particle averaging its heading with neighboring individuals under the influence of a uniform noise.However, recent empirical studies have found that biological groups might exhibit more complex patterns.In [8,21], the successive steps taken by ants and caribou have been found to be correlated.Similarly, the turn rate dynamics of a flagtail fish has been modeled as a mean reverting stochastic process in [16], under the premises of the socalled persistent turning walker (PTW) model, see also [16,43].Different from these larger fish, zebrafish locomotion is characterized by episodes of fast and large changes of direction [20] (Figure 1) that compose a unique burst-and-coast swimming style [14].To model such swimming behavior, we have proposed in [29] the jump persistent turning walker model (JPTW), which consists of a stochastic mean reverting model with jumps.The model comprises four salients parameters that, upon calibration on empirical data, yield accurate predictions for zebrafish turn rate dynamics [29,32].
To describe zebrafish collective behavior, in [30] we have proposed a new form of noise based on a discrete time approximation of the JPTW.The turning rate noise takes the form of an autoregressive process that includes jumps modeled as a discrete time compound process combining a Bernoulli random variable that controls the jump frequency and a normal random variable that determines the jump intensity.With a given probability, the turn rate noise is driven by a Gaussian process, otherwise subject to an additional turn rate increment.Variants of such process are commonly used in finance to model sudden rises or drops in the price of financial assets [18]; in physics to capture the extreme deviation in particle kinematics [27]; and in biology to explain the foraging behavior of organisms [11], the self-organized vortices of microtubules [39], and the behavior of starling flocks [5].
In this work, we propose the integration of this new form of noise into the classical vectorial network model (VNM) [3,34], which has been proposed as a valid and mathematicallytractable approximation of the more complex Vicsek model in the case of uniform noise.Within this approximation, each individual can interact with any other individual, irrespective of their location.As such, the VNM should be considered as a feasible approximation of the Vicsek model under the assumption that individuals are so fast that they completely mix at each time step.This model differs from the case study considered in [30] in the mode of selection of the neighbors, which is completely random in the vein of the traditional VNM [31,35].Specifically, the approach considered in [30] is based on the notion of numerosity-constrained networks [2], where each individual consistently uses its prior heading for deciding its next heading.In this case, we instead posit that prior information is not necessarily used in the decision process, which draws a given number of indi-viduals for estimating an average heading.In the particular context of swarm robotics, individual robots may not have full, precise access to their prior state due to a variety of factors, such as external perturbations, intrinsic noise, or component failure.From a mathematical point of view, this difference results in a distinctive formulation of the problem, which warrants a separate treatment.
In the following, we first introduce the JPTW along with the proposed model of collective behavior incorporating the turn rate noise.Then, we linearize the VNM in the neighborhood of a synchronized state, in which all the individuals share a common heading and study the steady state deviation of the linearized system.We propose a closed-form expression for an order parameter of the collective behavior [34].Using numerical simulations, we verify the closed-form expression and we parametrically investigate the effects of the turn rate noise parameters on group coordination.

MODELING COLLECTIVE BEHAVIOR WITH TURN RATE NOISE
In this section, we recall the JPTW model introduced in [29] followed by its discrete time approximation that is used to model the turn rate noise.The new form of noise is then utilized to explain the VNM with turn rate noise.

The jump persistent turning walker
The JPTW proposed in [29,30] is a mean reverting stochastic differential equation (SDE) that incorporates jumps to capture the fast and large change of directions observed in zebrafish swimming.The SDE of the turn rate process ωt(rad s −1 ), where the subscript t (s) identifies the continuous time, takes the form where dt is an infinitesimal time increment; β (s −1 ) is the relaxation rate controlling the gradual return of turn rate to zero after a change of direction; σ (rad s −1 ) quantifies the variability of the turn rate process, and Wt is a standard Wiener process, defined such that dWt are independent and identically distributed (i.i.d.) Gaussian processes with zero mean and variance dt; Jt is the jump term defined as a compound Poisson process Jt = νt j=1 Zj, where Zj 's are i.i.d.Gaussian random variables with zero mean and standard deviation γ and νt is a counting process.For time r ≤ t, νt −νr is a Poisson random variable with parameter λ(t − r).The intensity of the jumps is γ (rad s −1 ) and their frequency is λ (s −1 ).
For convenience, the JPTW hypothesizes that at most a single jump is observed in a small time step ∆t [29].This simplification is used to approximate the Poisson random variable ν(k∆t) − ν((k − 1)∆t) with a Bernoulli random variable with parameter λ∆t [6].Such a quantity measures the probability of a jump in a single time increment.The discrete time process in Eq. ( 2) is used to formulate the likelihood probability density function (pdf) of the turn rate with parameters (β, σ, γ, λ), that is, where m = ω(k − 1)e −β∆t and φ(ω; m, v) is the pdf of a Gaussian random variable with mean m and variance v.In the absence of jumps, that is λ = 0, the right hand side of Eq. ( 3) corresponds to the pdf of the PTW.For λ = 0, jumps are present with probability λ∆t, and the pdf of the turn rate is augmented through the second summand in the right hand side in Eq. ( 3).
Model parameters are identified by maximizing the loglikelihood function T k=2 log f ) over a given time series {ω(k)} T k=1 of zebrafish turn rate of T points.The JPTW calibrated on experimental data of zebrafish swimming has been shown to fit the empirical distribution of zebrafish turn rate using quantile quantile plots [29,32].Moreover, using a likelihood ratio test, the model has also been found to provide a better prediction of zebrafish turn rate dynamics compared to the PTW.

The vectorial network model with turn rate noise
In the VNM, an individual i = 1, . . ., N with heading angle θi(k) at time step k = 0, 1, . .., interacts with a subset of K randomly selected neighbors [3,34].In a discrete time setting, the orientation updating rule of individual i is given by where the operator Arg(•) returns the angle of a vector; ve Iθ j (k) is the average orientation that is used by individual i to update its heading angle, with I being the imaginary unit and v being a constant; ω1, . . ., ωN ∈ [−π, π) are i.i.d.uniform random variables; η is a nonnegative scaling factor of the noise intensity.The updating rule in Eq. ( 2) implies an absence of temporal correlation [30] and longrange interactions [4].
To account for the distinctive behavior of zebrafish, we replace the uniform noise in Eq. ( 4) with the turn rate noise in Eq. ( 2) by setting δ = 1, that is, where for simplicity, we have omitted the dependence of εi and τi on the time step k and we have introduced α = 1 − e −β∆t , with e −β∆t being the one-step correlation coefficient of ωi.The turn rate noise is constructed such that a Gaussian noise determines the individual orientation with probability 1−λ, and an additional turn rate increment with intensity γ is observed with probability λ.The turn rate noise defined above can be shaped to a desired distribution by adjusting the parameters α, σ, λ and γ as illustrated in Figure 2. The VNM with turn rate noise can be written as: where for simplicity, we have set v = 1.To analyze group coordination, we use the well-known polarization [34] defined by: where E[•] indicates expectation.The polarization is equal to 1 when all individuals have the same heading.This is the case when noise is absent from the system.On the other hand, the polarization approaches 0 as coordination is lost.
In the case of uniform additive noise, group coordination has been investigated for various group size N , number of connected neighbors K, and noise intensity η [3,34].A linear approximation of the problem has also been proposed in [35] and in [31] by considering a fraction of individuals as group leaders.
For small networks, independent of the number of connected neighbors K, the polarization has been shown to reach a plateau for large noise intensity [35].For large networks, a phase transition from order to complete disorder has been found for large noise intensity [3,34].Also, as the number of connected neighbors increases, the noise intensity required to drive the system to complete disorder has been shown to increase [35].

STEADY STATE DEVIATION
We introduce ζi(k) = θi(k) − θ0, for i = 1, . . ., N where θ0 defines a synchronized state, to linearize Eq. ( 6) as follows: We introduce the compact notation ζ(k) = [ζ1(k), . . ., ζN (k)] T and ω(k) = [ω1(k), . . ., ωN (k)] T , where T indicates matrix transposition.Using this formalism, the linear stochastic system in Eq. ( 8) can be compactly rewritten as where ε = [ε1, . . ., εN ] T and τ = [τ1, . . ., τN ] T are N × 1 vectors.Matrix W in Eq. ( 9) defines a random variable with rows that are i.i.d.vectors with K randomly selected entries taking value 1/K, while any other entry equals 0. A variant of the problem defined in Eq. ( 9) with uniform additive noise has been addressed in [35], where a closed-form expression for the polarization was first derived.Following a similar line of arguments, we compute the steady state deviation of the VNM with turn rate noise toward establishing a new closed-form expression for the polarization.
We multiply both sides of Eq. (9a) by the projection matrix N to obtain a stochastic system for disagreement dynamics, that is, where ξ(k) = Rζ(k) is the disagreement vector.For any time step k, we define the expected value of the disagreement norm as δ(k) = E ξ(k) 2 .Similar to [35,31], we investigate the mean square behavior of the system in Eq. ( 10) through the second moment matrix Ξ(k) = E ξ(k)ξ(k) T whose trace corresponds to the disagreement norm δ(k).
When it exists, we define the steady state mean square deviation as δ∞ = lim k→∞ δ(k).In the absence of noise, that is, η = 0, Eq. (9a) is mean square stable and δ∞ = 0 independently of the initial condition.In such case, the disagreement dynamics is said to be asymptotically mean square stable.
Similar to [30], we evaluate the second moment matrix Ξ(k) using a well-known property of the Kronecker algebra [9], that is, vec (ABC) = (C T ⊗ A)vec (B), where A, B, C are matrices with appropriate dimensions, ⊗ is the Kronecker product, and vec(•) indicates vectorization of a matrix by stacking its columns.Thus, we obtain the following recursive relation for the second moment matrix: which is identical to the one computed in [35].
Mirroring the line of arguments of [35], the expected values of W and W ⊗ W are given by E where 1N is the N -dimensional column vector of all ones.Using these expressions, it can be verified that where 0N is the N -dimensional column vector of all zeros.
We iterate the expression in Eq. ( 11) from time 0 to time k, and we use the independence between W and ω(s) for s ≤ k to obtain vec (Ξ(k + 1)) = G k+1 vec (Ξ(0)) ( 12) where, without loss of generality, the initial value of the turn rate vector is set to zero, that is, ω(0) = 0N .
Accounting also for the independence between ε and τ , the variance of the random variable εi + γτi is computed as 1 + λγ 2 .Similar to [30], from the relation in Eq. ( 10b) and again the independence between ε and τ , the evolution of the second moment matrix is computed as By replacing Eq. ( 13) into Eq.( 12) and iterating over k, we find In [35], we have computed the spectral radius of G as follows The two series in Eq. ( 14) converge if and only if ra < (1 − α) 2 < 1 and the limit is equal to The steady state deviation is evaluated by right multiplying ( 16) by vec (IN ) T to obtain This expression can be evaluated in closed-form in two different ways, using the fact that [35], or considering Proposition 2 of [1].In either cases, we can derive

EFFECT OF THE TURN RATE NOISE ON GROUP COORDINATION
In this section, we study the effects of the turn rate noise on group coordination based on a new closed-form expression for the polarization for small noise intensity, which is complemented and validated by extensive numerical simulations.

Linear approximation of the polarization
To examine the effect of turn rate noise on polarization in the VNM, we linearize Eq. ( 7) in the neighborhood of the synchronized state by substituting θi(k) = θ0 +ζi(k), using a Taylor series expansion of cos(•) and sin(•) in the expression for e I(•) , and considering a first order expansion of 1 + (•), to find When 0 < ra < (1−α) 2 < 1, we replace δ∞ by its expression in Eq. ( 18) to obtain Using the closed-form expression in Eq. ( 20), it is evident that for small noise perturbations, the polarization is a decreasing quadratic function of the noise scaling factor η.  Similar to [30], for small noise, the polarization is also a decreasing function of the jump frequency and the jump intensity.Given the inequality K(N−1) N(K−1)+1 < 1, the net role of the turn rate noise parameters is independent of the group size N and the number of connected neighbors K. Specifically, the presence of jumps induces an additional variance of λγ 2 that contributes to further reducing polarization.The noise autoregressive coefficient has a quadratic effect on group polarization with an extremum observed at α = 1.
Beyond its utility in studying the effect of the turn rate noise on the polarization, the closed-form expression in Eq. ( 20) also enables the study of the behavior of the system in the thermodynamic limit.When N → ∞, the order parameter can be approximated as Pol In addition, for K ≫ 1, the polarization reduces to Pol ≃ 1 − η 2 1+λγ 2 1−(1−α) 2 , and is only a function of the turn rate noise parameters.In such case, as the frequency of jumps and their intensity become substantial, a very small noise intensity is required to maintain the system in an ordered state, that is, η ≪ 1−(1−α) 2 1+λγ 2 .

Analytical and numerical results
We further investigate the effect of the turn rate noise on group coordination through numerical simulations for small (N = 10) and large (N = 5, 000) networks.The simulations for the small group are executed for 10, 000 time steps, and the polarization is computed by averaging the last 5, 000 time steps.The simulations for the large group (N = 5, 000) are performed for 2, 000 time steps while averaging the last 1, 000 time steps to compute the polarization.The initial headings of the individuals are randomly set prior to running the simulations.
Figure 3 presents the polarization as a function of the jump frequency and jump intensity.In agreement with theoretical expectations, the polarization is a decreasing function of both these parameters.Compared to the jump amplitude, the jump frequency tends to have a greater impact on group coordination.Displaying the polarization as a function of the relaxation rate and the jump amplitude, see Figure 4 helps illustrating that polarization decreases with the jump amplitude.Moreover, the polarization shows an inverted U-shape as α varies.This confirms the findings obtained through our new closed-form expression, where for small noise values, polarization is related to an inverted quadratic function in α.Also, different from the case studied in [30], where an individual always includes itself in its K connected neighbors, the turn rate noise parameters more severely affect the VNM with turn rate noise.
The effect of the noise scaling factor η on the polarization for fixed noise parameters α = 1.71, λ = 0.1, and γ = 1 is presented in Figure 5.For η ≤ 0.45, our prediction in Eq. ( 20) is in very good agreement with numerical simulations for small (Figure 5(a)) and large networks (Figure 5(b)).It is also observed that the polarization is a decreasing function of the noise scaling factor and reaches a plateau estimated at around 0.28 for the small networks in Figure 5(a).Such a plateau is independent of the number of connected neighbors and cannot be captured by the closed-form expression in Eq. ( 20), which holds only for small values of η.When the noise scaling factor increases, for large networks shown in Figure 5(b), the system exhibits a phase transition from a state where individuals share the same direction to a state where their heading is completely random.Similar to the traditional VNM [35], polarization is a decreasing function of the number of connected neighbors as opposed the case considered in [30] where such a conclusion holds only η ≥ 0.66.

CONCLUSIONS
In this work, we have analyzed the VNM with a new form of noise inspired by the zebrafish model organism [29].This new form of noise is modeled as a discrete time autoregressive process with jumps, which is derived from a discrete time approximation of the mean reverting stochastic jump diffusion process proposed in [29].In the model of collective behavior considered herein, the mode of selection of connected neighbors is random and identical to the traditional VNM with uniform noise.This process of network construction fundamentally differs from the instance studied in [30], where an individual always includes itself in its set of connected neighbors.
Group coordination within the VNM with turn rate noise has been thoroughly investigated through a novel closedform solution that is valid for small perturbations from a synchronized state and independent numerical simulations.The number of connected neighbors is set as follows: K = 2 in green, K = 4 in blue, and K = 8 in red.Dots are numerical simulations, and solid lines are analytical predictions from Eq. (20).For the small network, the plateau is reached at Pol ≃ 0.28.
The analytical solution provides good predictions of the polarization for small noise intensities, where the polarization monotonically decreases with the square of the noise scaling factor.Similar to variants of the VNM studied in [30,31,34,35], when the noise intensity increases, a phase transition is observed for large group size while a plateau in the polarization is found for small groups.Compared to our previous results in [30] for the VNM−ω, group coordination in the VNM with turn rate noise is more affected by the noise, whereby we observe a complex interplay between the parameters shaping the turn rate noise and the group coordination.
Our modeling framework proposed here is expected to aid in the analysis of zebrafish sociality and provide a mathematical foundation for the study of collective behavior of other species, characterized by complex individual dynamics.The analysis could also be helpful in the study of engineering systems, such as swarm robotics, in which information may be intermittently shared across units under the effects of noise.

Figure 1 :
Figure 1: Zebrafish burst-and-coast swimming in successive video frames samples at 30 frames per second.(a) A single tail flick can propel the fish forward and (b) fish can change its orientation in just a few milliseconds.

Figure 5 :
Figure 5: Polarization as a function of the noise scaling factor η for small (a) N = 10 and large (b) N = 5, 000 networks with α = 1.71, λ = 0.1 and γ = 1.The number of connected neighbors is set as follows: K = 2 in green, K = 4 in blue, and K = 8 in red.Dots are numerical simulations, and solid lines are analytical predictions from Eq.(20).For the small network, the plateau is reached at Pol ≃ 0.28.