An agent-based approach for modelling collective dynamics in animal groups distinguishing individual speed and orientation

1. Introduction

Collective organization and movement occur in many biological and ecological systems, from colonies of bacteria to human crowds [1–3]. In particular, many animal species form groups and are able to undergo synchronized and coordinated dynamics, which represent an evolutionary adaptation conferring certain collective benefits, including more efficient new nest exploration and colonization [4], food gathering [5], predator avoidance [6,7] and heat preservation [8], as schematically shown in figure 1.

• Download figure
• Open in new tab
• Download PowerPoint

The description of animal group movement and patterning has increased in the last decades the multidisciplinary interest of various research communities, including applied mathematicians and physicists. In this respect, a wide range of approaches are present in the theoretical and computational literature, as reviewed for instance in [1,9,10]. In particular, continuous models are characteristic of a macroscopic point of view and rely on the definition of a proper density of agents, whose dynamics are set to obey conservation equations and phenomenological assumptions for their closure [11–15]. Hydrodynamic arguments and Boltzmann-like evolution laws for statistical distributions of individual position and velocity are instead at the basis of mesoscopic kinetic approaches [16–18]. Finally, discrete models and cellular automata employ a microscopic perspective and describe the system of interest as a set of single/isolated elements, which are individually considered and behave according to individual rules (based, for instance, on energy minimization or on Newtonian laws) [19–25]. For the sake of completeness, we recall that spatial organization and collective motion of systems of animals can be also approached with lattice-based models. This is the case, for instance, in the series of works by the group of Kamimura (refer to [24,26–29]) and of the so-called swarm lattice-gas cellular automata introduced by Deutsch in [30].

The evolution of an animal group is here tackled with a particle-based approach: it belongs to the family of discrete methods being characterized by the fact that each agent is individually represented by a material point, whose velocity varies according to an ordinary differential equation. Our general model (that is not referred to a particular species) realistically distinguishes between individual speed and direction of motion. The former is in fact typically established by the purpose of movement and/or by physiological limitations, the latter results instead from the competition between different behavioural stimuli, mainly coming from the environment and from interagent interactions. In this respect, we here account for repulsion, attraction and alignment inputs that are active in non-overlapping individual neighbourhoods. Escape strategies are also included in the case of heterogeneous systems, i.e. formed by a predator chasing a group of prey. Remarkably, each of the considered behavioural stimuli is here simply defined by an orientation (unit) vector and a weight, which is correlated to an individual preference. A constraint on the sum of the weights is finally given to avoid simultaneous minimization/maximization of all individual movement traits.

Our model is not completely new: rather, it combines and modifies at various extents concepts and ingredients already present in the literature. More specifically, focusing on the similarities and differences of our mathematical framework with respect to analogous published approaches, we first start from the second-order model proposed by Cucker & Smale in [31,32] to study flocking phenomena (with the terminology ‘second-order approach’, we hereafter mean a model based on the second Newtonian law, that relates the acceleration of a particle to the acting forces). Therein, the authors only account for an alignment mechanism that forces each individual to adjust its orientation to the groupmates. In particular, the tendency to move in the same direction is set to be higher for close enough individuals, whereas it decreases in the case of pairs of distant agents. An interesting extension of this approach is discussed in [1] and deals with role differentiation within the animal population. In particular, the emergence of some instantaneous leaders (defined as individuals whose movement is not affected by others) can be observed by introducing anisotropy in the synchronization process, i.e. by assuming that each animal interacts only with the groupmates falling within a given visual cone. A further proposed model development is the addition of stochastic perturbations in individual dynamics. As commented again in [1], a noise contribution can be also included in the second-order particle model developed by D’Orsogna and coworkers which, in its original version, accounts for a self propulsion term, a friction contribution (following Rayleigh’s Law) and a classical Morse potential for attractive/repulsive pairwise interactions [33].

The analysis and classification of the collective motion of identical interacting particles is also the topic of the well-known Vicsek’s model [34]. In particular, it assumes that flocking phenomena are owing to alignment mechanisms, with uncertainties arising from a random contribution in individual velocity. This approach is then extended in [35] with the introduction of a Lennard–Jones-type potential to represent repulsive and attractive dynamics. Interestingly, each velocity contribution is therein weighted by a parameter that, at variance with our work, can independently vary in $\mathbb{R}$. Another difference with respect to our model relies on the fact that the extension of the interaction regions are significantly small so that the particles interact strictly locally.

The influence of a hunter on group collective dynamics is deeply investigated by particle-based models as well. For instance, in [20], each prey is assumed to be subjected to an increasingly linear attraction towards groupmates and to a hyperbolic repulsion both towards groupmates and towards the predator. The cohesive term and the velocity component relative to hunter avoidance are multiplied by coefficients varying in ${\mathbb{R}}^{+}$, whereas the repulsion between prey individuals has a constant unitary weight. Predator dynamics are then completely determined by an attraction term which has a hyperbolic or a more than hyperbolic law. In particular, the hunter is simultaneously attracted by all prey, i.e. it is confused according to the terminology that will be introduced later on.

In [36], predator–prey interactions are instead tackled by a second-order particle approach, eventually applicable to the case of birds attacking crabs and whales attacking small fishes. In particular, Lee and coworkers assume that prey individuals are subjected to alignment, attraction, repulsion, friction and hunter escape forces. A random component is accounted for as well. In particular, the friction term, proportional to the present individual speed, is employed to prevent prey from moving too quickly (this is not necessary in our approach as the speed is possibly limited by a threshold value). The hunter is finally set to point towards the centre of mass of the group of prey (i.e. it has a confused strategy). Differently from our approach, each behavioural trait included in their model has a distinct form.

The approaches reviewed so far share with our model the main principles underlying animal movement, i.e. repulsion, attraction, alignment and eventually predator avoidance. However, some relevant differences emerge: in fact, in all these works, the speed closely depends on the modulus of the vector establishing the direction of motion and/or each contribution taken into account has a substantially distinct mathematical law. These two aspects make our work conceptually closer to the following group of microscopic models dealing with group dynamics.

In 1987, Reynolds published a well-known agent-based approach to simulate flocking phenomena of bird-like individuals (named with the acronym BOIDs) [37]. Therein, each particle is moved by attraction, repulsion and alignment stimuli. In particular, each behavioural trait is mathematically given by an orientation unit vector multiplied by a weight (as in our case). However, in the BOID model, the weights are allowed to independently vary, i.e. no constraint on their relationship (e.g. on their sum) is established. An agent may therefore simultaneously maximize or minimize all its behavioural stimuli which, according to us, is somehow unrealistic. Furthermore, in the approach by Reynolds, the modulus of the resultant velocity vector gives the agent speed, if a threshold values is not exceeded. In this perspective, speed and orientation are not decoupled. The BOID model is then extended in [38] by including an additional escape rule with the same mathematical structure as the other three classical contributions. Further, Delgado-Mata and coworkers multiply each velocity component by a factor that accounts for the emotional state of the animals (e.g. fear). Also in this case, the velocity vector of each agent establishes its speed.

In the work by Couzin and colleagues [21,39], the velocity of each individual is instead defined by a constant speed and an orientation unit vector, exactly as in our approach. However, the direction of movement is determined in a slightly different way. In fact, a generic agent, characterized by an anisotropic visual cone, is set to be subjected only to repulsion if it detects another mate within its avoidance region. Otherwise, it is set to be subjected to both attraction and alignment towards the groupmates falling within the corresponding interaction areas: in this case the resultant orientation vector is given by the mean of the two contributions. Perturbations are employed by modifying the direction of movement at a randomly chosen extent: however, an agent is not allowed to undergo a complete change in its orientation, owing to a constraint on the turn-around angle. At variance with our approach, Couzin and colleagues therefore do not weight individual behavioural preferences; rather, repulsion is the primary stimulus, whereas alignment and attraction equally affect animal dynamics.

Couzin and another group of coworkers perform and describe in [40] an interesting experiment, where a predator bluegill sunfish is allowed to hunt mobile virtual targets in a controlled environment. In particular, the prey agents are assumed to move following almost the set of rules defined in [21,39]. Each agent is in fact subjected either to repulsion or to the balance between three traits: it can align, ignore and be attracted by the groupmates. In this respect, the resultant direction of motion of a virtual prey, in the absence of repulsion, is given by a weighted sum of the three contributions, where the weights, as in our case, sum up to one. However, differently from our work, such three behavioural inputs are simultaneously present within the same interaction region (for instance, in our case the alignment and the attraction areas do not overlap).

The approaches by Couzin and colleagues are extended in [41] to describe escape and foraging manoeuvres in a selfish herd. In particular, all agents are also set to move away from the position of the predator and/or towards the location of a food source. Further, each velocity contribution is therein multiplied by a weight, that is allowed to freely and independently vary: this is in analogy to the Reynolds’ BOID model but differently from our approach. Therefore, also in this case, implausible situations may in principle occur: e.g. a prey may simultaneously maximize the competing stimuli of predator avoidance and food search. Finally, in the work by Wood & Ackland, the predator points to the closest prey, i.e. it is non-confused.

Summing up, with respect to similar published approaches, the proposed model has the following distinct features: (i) individual speed and direction of movement are clearly independent; (ii) the description of animal dynamics does not require the introduction of complex laws of motion and/or of algorithmic rules; rather all behavioural stimuli are defined by the same simple mathematical form (i.e. by an orientation vector and a weight); (iii) only a minimal set of parameters, which have a direct empirical meaning as well, is necessarily introduced; (iv) each agent behaves according to a hierarchy of behavioural preferences (i.e. it is not allowed to maximize/minimize the response to all movement inputs); and (v) the simplicity of our method allows us to provide not only numerical but also analytical insights.

Structure of the work. The rest of the paper is then organized as it follows. In §2, we will present the mathematical framework. Section 3 will show the model ability to reproduce and classify collective movement and patterning of an animal population emerging from different individual behavioural preferences. Section 4 will then focus on a coupled analytical/numerical study of the equilibrium configurations of our animal system when subjected only to attractive and repulsive interactions. A possible model extension accounting for the presence of a predator which impacts on group dynamics will be proposed and investigated in §5. Finally, some hints for further applications and future developments of our approach will be given in the conclusive §6.

2. Mathematical model

A generic group of animals is described by a particle system within the d-dimensional space ${\mathbb{R}}^d$. In this perspective, the ith agent (with i = 1, …, n, n being the total number of individuals) is represented by a material point with position x_i(t). Individual dynamics are then described by the following first-order model:

which can be derived from a generic second-order Newtonian approach under the assumption of overdamped force-velocity response (a consistent hypothesis for living entities, see [2,42,43] for comments). Equation (2.1) can be derived from a more general first-order model (see appendix A) and effectively decouples the magnitude of the velocity of the ith animal, given by the scalar $v_i(t)\in {\mathbb{R}}_{+}$ (possibly, $v_i(t)\in [0,v_{\max }]$ to account for physiological limitations), from its direction, given by the unit vector ${\displaystyle \frac{{\mathit{\omega }}_i(t)}{|{\mathit{\omega }}_i(t)|}\in B_1^d\subset {\mathbb{R}}^d}$ (where $B_1^d$ denotes the d-dimensional ball, centred at the origin, with unitary radius). As a relevant feature of the proposed model, these two quantities are assumed to be independent, because they have distinct physical meanings. For instance, the individual speed v_i(t) can be determined by the intrinsic status of the agent (e.g. calm or in a hurry) as well as by the purpose of its movement and by physiological/physical limitations. On the other hand, the velocity orientation of the particle, ω_i(t)/|ω_i(t)|, can be formally correlated to the competition between different behavioural stimuli (triggered, for example, by environmental signals, own preference or interactions with the surrounding agents). In this respect, the vector ω_i(t) can be determined by a weighted sum of proper contributions, i.e. where ${\mathcal{J}}_i$ is the set of behavioural inputs that influence the dynamics of i while the unit vectors ${\mathbf{w}}_i^{\hspace{0.17em}j}\in B_1^d\subset {\mathbb{R}}^d$ define the corresponding orientations. Finally, the coefficients ${\alpha }_i^j(t)$ are weights that describe the relative importance of each stimulus with respect to the others; they may also evolve in time owing to variations of internal or external conditions, e.g. visibility. In particular, to have comparable effects on agent dynamics and to avoid simultaneous maximization/minimazion of all behavioural traits, we assume that for all i = 1, …, n. In this respect, relation (2.2) is a linear combination of the directional cues that affect the movement of the ith agent.

Given the generical mathematical structure, we now recall that the majority of animal groups undergoes collective motion mainly guided by interindividual social interactions. In particular, three fundamental behavioural rules are classically identified, as explained in [37,44–46] and reviewed in [1] in the case of swarming populations. First, a short-range repulsion describes the tendency of each animal to maintain a minimum comfort space within the group: it is here implemented by a directional unit vector that allows the ith agent to move away from close enough neighbours:

where $d_{\mathrm{rep}}$ being the extension of the avoidance region.

A middle-range alignment mechanism then reproduces the aim of each agent to synchronize its movement, in terms of orientation, with the groupmates falling within a given surrounding area of radius $d_{\mathrm{align}}$, i.e.,

In equation (2.6), $\mathrm{M}({\mathit{\omega }}_j(t)/|{\mathit{\omega }}_i{(t)|)}_{\hspace{0.17em}j\in {\mathcal{N}}_i^{\mathrm{align}}(t)}$ denotes in fact the mean of the actual velocity directions evaluated over the set of particles A long-range attraction finally enters the picture if the ith animal falls too far apart from the rest of the group and tries to reach again its mates: in mathematical terms, we have indeed with The extension of the attraction region, $d_{\mathrm{attr}}$, can be also interpreted as the visual distance of the animal of interest as done for instance in [47]. A schematic representation of the different interparticle interaction areas is proposed in figure 2a.

• Download figure
• Open in new tab
• Download PowerPoint

Putting the introduced velocity components in equation (2.2), ω_i results specified as

being ${\mathcal{J}}_i=\{\text{repulsion, alignment, attraction}\}$ and ${\alpha }_i^{\mathrm{rep}}(t)+{\alpha }_i^{\mathrm{align}}(t)+{\alpha }_i^{\mathrm{attr}}(t)=1$ for any i = 1, …, n. As written, the proposed approach includes a minimal set of parameters: individual speeds and extensions of interparticle interaction regions, that can be easily quantified for a given animal species, as well as the intensities of the behavioural inputs, that need at least a hierarchical estimate but have clear biological meaning. No other strictly technical coefficients need to be introduced.

We finally propose a proper dimensionless form of the model, in order to facilitate its application to any group of animals, regardless of their characteristic measures. For this purpose, we scale lengths with the repulsion radius $d_{\mathrm{rep}}$, velocities with the characteristic speed of the agents of interest, say v, and times with $d_{\mathrm{rep}}/v$. Trivial calculations allow then to rewrite equations (2.1) as

where the dimensionless variables are overlined. In particular, the extensions of the interaction regions scale as ${\overline{d}}_{\mathrm{rep}}=1$, ${\overline{d}}_{\mathrm{align}}=d_{\mathrm{align}}/d_{\mathrm{rep}}$ and ${\overline{d}}_{\mathrm{attr}}=d_{\mathrm{attr}}/d_{\mathrm{rep}}$, so that ${\overline{d}}_{\mathrm{align}}$ and ${\overline{d}}_{\mathrm{attr}}$ measure the relative extension of the corresponding areas with respect to the radius of the avoidance neighbourhood.

3. Simulation details and results

The proposed model is employed in a planar setting, i.e. d = 2. In all forthcoming simulations, the animal population of interest is formed by n = 100 individuals. In particular, as shown in figure 2c, they are initially arranged in an almost round configuration with diameter equal to 6 (in non-dimensional terms) and randomly assigned positions. The initial velocity direction of each generic ith agent, i.e. ${\overline{\mathit{\omega }}}_i(\overline{t}=0)$ is randomly established as well. Given that the proposed scaling results in a unitary dimensionless repulsion radius ${\overline{d}}_{\mathrm{rep}}=1$, we set the linear extension of the alignment and the attraction regions equal to ${\overline{d}}_{\mathrm{align}}=8$ and ${\overline{d}}_{\mathrm{attr}}=30$, respectively. The ratios between the depth of the different interacting areas are consistent with the biologically plausible ranges explored by Couzin et al. [39] and also employed by Wood & Ackland [41]. In this respect, as also commented in the conclusive section, the application of the model to a specific animal species would lead to a more proper parameter setting. We furthermore assume that all individuals have a constant and common speed equal to the characteristic value v: this implies a non-dimensional quantity ${\overline{v}}_i(\overline{t})=1$, for any i and $\overline{t}$. The weights αs are finally hypothesized to be independent from both time and animal individuality, i.e. ${\alpha }_i^{(\cdot )}(\overline{t})={\alpha }^{(\cdot )}$ for any i and $\overline{t}$, being ( · ) ∈ {rep, align, attr}. In particular, the resulting set of permitted triplets (i.e. those satisfying constraint (2.3)) identifies the simplex $\mathcal{S}$, represented by the triangle in figure 3a.

• Download figure
• Open in new tab
• Download PowerPoint

We now analyse how individual behavioural preferences affect group dynamics. In more detail, our strategy is to vary the values of the α-coefficients and to numerically explore the large-time states of the animal system. In this perspective, we introduce the following classification for the simulation outcomes.

Definition 3.1.

(i) The animal group has a time-asymptotic well-spaced/non collapsed configuration if the following condition is satisfied:

$${\overline{d}}_{\min }=\underset{\overline{t}\to \mathrm{\infty }}{lim}\underset{\genfrac{}{}{0ex}{}{(i,j)\in \{1,\dots ,n\}\times \{1,\dots n\}}{i\ne j}}{min}|{\overline{\mathbf{x}}}_i(\overline{t})-{\overline{\mathbf{x}}}_j(\overline{t})|\in ({\overline{d}}_{\mathrm{rep}}-ϵ;{\overline{d}}_{\mathrm{rep}}+ϵ).$$3.1

(ii) The animal group has a time-asymptotic synchronized movement if the following condition is satisfied:

$$\sigma =\underset{\overline{t}\to \mathrm{\infty }}{lim}{n}^{-\frac{1}{2}}\sqrt{\sum _{i=1}^n{|\frac{{\overline{\mathit{\omega }}}_i(\overline{t})}{|{\overline{\mathit{\omega }}}_i(\overline{t})|}-\mathrm{M}{\left(\frac{{\overline{\mathit{\omega }}}_j(\overline{t})}{|{\overline{\mathit{\omega }}}_j(\overline{t})|}\right)}_{\genfrac{}{}{0ex}{}{\hspace{0.17em}j=1,\dots ,n}{j\ne i}}|}^2}<{\sigma }_{\mathrm{c}},$$3.2

where the overlines identify, as usual, dimensionless measures. Further, the notation $\mathrm{M}{(\cdot )}_{\genfrac{}{}{0ex}{}{\hspace{0.17em}j=1,\dots ,n}{j\ne i}}$ indicates the mean of the velocity directions over the entire animal group.

These quantities are particularly suitable to describe the phenomenology of a set of interacting particles, in terms of both patterning and motion. In particular, condition (3.1) assures that there is no large-time individual overlap. In fact, if the asymptotic minimal interagent distance falls within the range $({\overline{d}}_{\mathrm{rep}}-ϵ;{\overline{d}}_{\mathrm{rep}}+ϵ)$, $ϵ$ = 0.05 being introduced to account for numerical errors in a multiparticle system, then all pairs of animals have an equal or a larger spacing, i.e. they do not collapse. Condition (3.2) instead assures that all individuals have almost the same direction of movement. In fact, preliminary simulations have shown that standard deviations σ of the distribution of the agent orientations smaller than the critical value ${\sigma }_{\mathrm{c}}=2$ imply a fully aligned animal movement. Similar classifications have been used in the literature to differentiate the evolutions of a wide range of multiagent populations, see for instance the case of swarming birds in [1] and references therein.

Our numerical study allows us to indeed subdivide the permitted parameter space $\mathcal{S}$ in four disjointed regions, according to definition 3.1. In particular, the subregion ${\mathcal{S}}_1$ is characterized by triplets of coefficients resulting in the asymptotic collapse of the cluster of particles, which also maintain an uncorrelated movement. The group of animals is instead observed to have a large-time synchronized collective motion, however coupled with an asymptotic collapse, when the intensities of the behavioural stimuli fall within ${\mathcal{S}}_2$. On the other hand, α^(·)-triplets belonging to the subregion ${\mathcal{S}}_3$ allow avoidance of unrealistic animal overlapping but the agent locomotion remains completely individual (i.e. no common direction emerges). ${\mathcal{S}}_4$ is finally characterized by the triplets of α^(·)s that result in an asymptotic non collapsed group configuration and collective alignment, i.e. all agents finally move in the same direction keeping a comfort spacing.

To further support the above dissertation, for each subregion ${\mathcal{S}}_i$ (i = 1, …, 4) of the permitted parameter space we show in figure 3b a representative asymptotic configuration of the system and in figure 7a the mean (with the corresponding standard deviation) of the quantities introduced in Definition 3.1, evaluated over a set of 10 numerical realizations resulting from randomly chosen α^(·)-triplets.

By reviewing the simulation results summarized in figure 3a, we can also notice that below a given value of the alignment parameter ${\alpha }^{\mathrm{align}}$ (i.e. ≈ 0.2), a synchronized movement of the animal group is not obtained regardless of the intensity of the other stimuli. On the opposite, a directionally collective movement is captured for any pair $({\alpha }^{\mathrm{rep}},{\alpha }^{\mathrm{attr}})$ if ${\alpha }^{\mathrm{align}}$ is larger than ≈ 0.6. Analogously, too low values of the weight ${\alpha }^{\mathrm{rep}}$ (i.e. ≤ 0.1) result in a particle collapse, independently from the magnitude of the other pair of α-parameters. On the other hand, the condition ${\alpha }^{\mathrm{rep}}>0.95$ is sufficient to have well-spaced large-time system configurations regardless of the values given to the other pair of coefficients. It is also interesting to notice that a complete balance between the three weights (i.e. ${\alpha }^{\mathrm{rep}}={\alpha }^{\mathrm{align}}={\alpha }^{\mathrm{attr}}$) does not result in a well-spaced system nor in a synchronized movement. Region ${\mathcal{S}}_4$ is rather characterized by a significant hierarchy of the intensity of the behavioural stimuli, with the alignment coefficient that substantially overcomes the repulsive one, which is in turn one order of magnitude larger than the attraction parameter. In this respect, it is finally useful to remark that (so far) we have not included the presence of elements that may cause dispersion of single animals. Therefore, a small attraction stimulus is sufficient to avoid group scattering.

For the sake of clarity, we remark that the numerical simulations proposed in this section run until $\overline{t}=100$: preliminary simulations have in fact shown us that the values of the quantities introduced in Definition 3.1 undergo only negligible oscillations if evaluated at longer times. Therefore, the particle configurations observed at $\overline{t}=100$ represent a consistent numerical approximation of the asymptotic system behaviour.

Comparable classifications of groups dynamics, done in the case of similar models, are obtained varying the extension of the interaction regions [39,41] or the parameters defining the magnitude and the shape of the individual velocity components [1,35].

4. Analysis of the equilibrium configurations of an animal group subjected only to attraction and repulsion

In the previous section, we have classified the collective dynamics of the animal group a posteriori, i.e. as resulting from numerical simulations. However, analytical arguments concerning the H-stability properties of agent-based systems allow us to derive parametric constraints that are a priori sufficient to avoid large-time individual collapse. In this respect, let us first neglect the alignment velocity contribution and focus only on attractive–repulsive stimuli to assure the existence of asymptotic equilibrium configurations. Under this assumption, the individual velocity vector ω_i(t), defined in equation (2.2), can be specified and rewritten in the following form:

where ${\alpha }_i^{\mathrm{rep}}+{\alpha }_i^{\mathrm{attr}}=1$ (cf. equation (2.3)) and with and It has been shown that the H-stability characteristics of the potential associated with the the interaction kernel $k^{\mathrm{int}}$, which will be hereafter denoted by $K^{\mathrm{int}}$, determine the large-time spatial organization of the particle system. In particular, if $K^{\mathrm{int}}$ is H-stable, the minimal interparticle distance at the equilibrium is bounded below by a finite fixed positive value (regardless of the total amount of particles n). In other words, if $K^{\mathrm{int}}$ is H-stable, the individuals do not asymptotically overlap. In this respect, a criterion to detect the H-stability of a scalar potential has been given in [48], and recently applied to the case of cell aggregates [49] and honeybee swarms [47] (see also [1,50] for further comments). Accordingly, we can state the following.

Theorem 4.1.

The interaction potential$K^{\mathrm{int}}$related to the pairwise interaction kernel$k^{\mathrm{int}}$defined in equations (4.2)–(4.4) is H-stable if the following parametric relation holds:

$$\begin{array}{rl}& \frac{d_{\mathrm{attr}}^3-d_{\mathrm{align}}^3}{d_{\mathrm{rep}}^3+d_{\mathrm{attr}}^3-d_{\mathrm{align}}^3}<{\alpha }_{\mathrm{rep}}<1,or\hspace{0.17em}equivalently,\\ & 0<{\alpha }^{\mathrm{attr}}<\frac{d_{\mathrm{rep}}^3}{d_{\mathrm{rep}}^3+d_{\mathrm{attr}}^3-d_{\mathrm{align}}^3}.\end{array}$$4.5

Proof.

Following the calculation introduced in [49], we have first to derive the potential $K^{\mathrm{int}}\hspace{0.17em}:\hspace{0.17em}\mathbb{R}\to \mathbb{R}$ associated to the kernel $k^{\mathrm{int}}$, i.e.

$$K^{\mathrm{int}}(r)=\{\begin{array}{ll}-{\alpha }^{\mathrm{rep}}\hspace{0.17em}r+C_1,& \text{if }0<r\le d_{\mathrm{rep}};\\ C_2,& \text{if }{d}_{\mathrm{rep}}<r\le d_{\mathrm{align}};\\ {\alpha }^{\mathrm{attr}}\hspace{0.17em}r+C_3,& \text{if }{d}_{\mathrm{align}}<r\le d_{\mathrm{attr}};\\ C_4,& \text{if }r>d_{\mathrm{attr}}.\end{array}$$

where r : = |x_j(t) − x_i(t)|. The constants of integration C₁, C₂, C₃, $C_4\in \mathbb{R}$ have then to be estimated in order to guarantee (i) the continuity of the potential $K^{\mathrm{int}}(r)$ and (ii) the fact that it is essentially negligible for large enough interparticle distances (i.e. ${\displaystyle \underset{r\to \mathrm{\infty }}{lim}{K}^{\mathrm{int}}(r)=0}$). Both conditions are in fact necessary hypotheses for Theorem 4.1, as clearly stated in [48]. By simple algebraic calculations, we have that C₄ can be taken equal to 0 and the other constants result

$$\begin{array}{ll}{C}_1& ={\alpha }^{\mathrm{attr}}{d}_{\mathrm{align}}-{\alpha }^{\mathrm{attr}}{d}_{\mathrm{attr}}+{\alpha }^{\mathrm{rep}}{d}_{\mathrm{rep}};\\ C_2& ={\alpha }^{\mathrm{attr}}{d}_{\mathrm{align}}-{\alpha }^{\mathrm{attr}}{d}_{\mathrm{attr}};\\ C_3& =-{\alpha }^{\mathrm{attr}}{d}_{\mathrm{attr}},\\ C_4& =0,\end{array}$$

so that the interaction potential rewrites as

$$K^{\mathrm{int}}(r)=\{\begin{array}{ll}-{\alpha }^{\mathrm{rep}}r+{\alpha }^{\mathrm{attr}}{d}_{\mathrm{align}}-{\alpha }^{\mathrm{attr}}{d}_{\mathrm{attr}}+{\alpha }^{\mathrm{rep}}{d}_{\mathrm{rep}},& \text{if }0<r\le d_{\mathrm{rep}};\\ {\alpha }^{\mathrm{attr}}{d}_{\mathrm{align}}-{\alpha }^{\mathrm{attr}}{d}_{\mathrm{attr}},& \text{if }{d}_{\mathrm{rep}}<r\le d_{\mathrm{align}};\\ {\alpha }^{\mathrm{attr}}r-{\alpha }^{\mathrm{attr}}{d}_{\mathrm{attr}},& \text{if }{d}_{\mathrm{align}}<r\le d_{\mathrm{attr}};\\ 0,& \text{if }r>d_{\mathrm{attr}}.\end{array}$$

Recalling the definition 1.1 in [48], we can affirm that $K^{\mathrm{int}}$ is H-stable when

$${\int }_0^{+\mathrm{\infty }}{K}^{\mathrm{int}}(r)\hspace{0.17em}r\hspace{0.17em}\text{d}r=\frac{1}{2}\underset{r\to 0}{lim}{K}^{\mathrm{int}}(r)>0.$$

In this respect, with some simple algebraic calculations, we obtain

$${\int }_0^{+\mathrm{\infty }}{K}^{\mathrm{int}}(r)\hspace{0.17em}r\hspace{0.17em}\mathrm{d}r=\frac{{\alpha }^{\mathrm{rep}}{d}_{\mathrm{rep}}^3+{\alpha }^{\mathrm{attr}}(d_{\mathrm{align}}^3-d_{\mathrm{attr}}^3)}{6}\hspace{0.17em}>\hspace{0.17em}0,$$

which is non negative if

$$\frac{{\alpha }^{\mathrm{rep}}}{{\alpha }^{\mathrm{attr}}}>\frac{d_{\mathrm{attr}}^3-d_{\mathrm{align}}^3}{d_{\mathrm{rep}}^3}.$$4.6

Finally, recalling that ${\alpha }_i^{\mathrm{rep}}+{\alpha }_i^{\mathrm{attr}}=1$, we can rewrite equation (4.6) only in terms of the repulsive (respectively, the adhesive) coefficient ${\alpha }_i^{\mathrm{rep}}$ (respectively, ${\alpha }_i^{\mathrm{attr}}$):

$$\begin{array}{rl}& \frac{d_{\mathrm{attr}}^3-d_{\mathrm{align}}^3}{d_{\mathrm{rep}}^3+d_{\mathrm{attr}}^3-d_{\mathrm{align}}^3}<{\alpha }_{\mathrm{rep}}<1,\mathrm{or}\hspace{0.17em}\mathrm{equivalently},\\ & 0<{\alpha }^{\mathrm{attr}}<\frac{d_{\mathrm{rep}}^3}{d_{\mathrm{rep}}^3+d_{\mathrm{attr}}^3-d_{\mathrm{align}}^3},\end{array}$$

which is exactly the thesis of the Theorem. □

Condition (4.5) allows us to subdivide the segment s defining the permitted parameter space $({\alpha }_i^{\mathrm{rep}},{\alpha }_i^{\mathrm{attr}})$ into two parts, each resulting in a distinct stability characterization of the interaction kernel $k^{\mathrm{int}}$, see figure 4. In particular, pairs of α-values falling within the dashed part s₂ of the segment satisfy equation (4.5) and assure an asymptotic non-collapsed configuration of the system. For the sake of completeness, the non-dimensional counterparts of constraints (4.5) reads as:

• Download figure
• Open in new tab
• Download PowerPoint

The above analytical considerations can be supported by numerical results. In this respect, we perform a pair of representative computational tests that reproduce the evolution of a homogeneous population of n = 100 component individuals, whose dynamics are only regulated by adhesive and repulsive stimuli. The values of the individual speed and of the extension of the interaction regions, as well as the initial condition of the system, are exactly the same as those employed in the previous §3. As shown in figure 4, if the weights α^(·)s do not satisfy the H-stability condition (i.e. fall within s₁), then the particle cluster asymptotically collapses, as confirmed by the corresponding ${\overline{d}}_{\min }=5\times {10}^{-4}$. On the other hand, if the values of the two coefficients satisfies the H-stability constraint (i.e. fall within s₂), then the animal group stabilizes in a well-spaced configuration characterized by the absence of individual overlap. Interestingly, in this case, the large-time minimal interparticle distance ${\overline{d}}_{\min }=0.9639$ falls within the range $({\overline{d}}_{\mathrm{rep}}-ϵ;{\overline{d}}_{\mathrm{rep}}+ϵ)$, introduced in Definition 3.1 (recalling that ${\overline{d}}_{\mathrm{rep}}=1$). We finally remark that the numerical realizations here proposed run until $\overline{t}=2000$, i.e. when a steady equilibrium configuration is reached.

The coupling between the analytical and the numerical results proposed in this section allows therefore to conclude that, in the absence of the alignment velocity component, an asymptotic well-spaced configuration of the system is a priori assured by H-stable repulsive/attractive kernels, i.e. by pairs of weights $({\alpha }^{\mathrm{rep}},{\alpha }^{\mathrm{attr}})$ satisfying condition (4.7). It is finally useful to remark that, as far as we know, a similar study has been never done in the case of interaction potentials characterized by constraints and interdependence between the characteristic parameters.

5. Prey–predator dynamics: model and results

We now incorporate heterogeneity in the animal system by introducing the presence of a single predator, labelled by the identification number i = 100 (so that the overall number of agents remains unaltered). In this respect, the remaining set of individuals hereafter constitutes a group of prey chased by the hunter. In nature, we can observe a wide variety of predator–prey interactions: however, a common characteristic is the emergence of empty space around the predator, owing to the obvious tendency of the other individuals to move away, as commented again in [20] and references therein and shown in figure 1c,d. This behaviour can be implemented in our model by the inclusion of a proper avoidance term in the prey dynamical rules: in this respect, we have that ${\mathcal{J}}_i=\{$repulsion, alignment, attraction, escape} and

for i ∈ {1, …, 99}. In particular, the directional velocity contribution ${\overline{\mathbf{w}}}_i^{\mathrm{escape}}$ reads as which is actually active if $|{\overline{\mathbf{x}}}_i(\overline{t})-{\overline{\mathbf{x}}}_{100}(\overline{t})|\le {\overline{d}}_{\mathrm{escape}}$, being ${\overline{d}}_{\mathrm{escape}}$ set equal to ${\overline{d}}_{\mathrm{attr}}$, see again figure 2a. The term in equation (5.2) is indeed a long-range repulsive contribution because it enters the picture as soon as the predator falls within the visual region of the ith prey. We in fact recall that the extension of the attraction region can be interpreted as an individual gaze depth. The weights of the behavioural stimuli included in equation (5.1) have finally to satisfy constraint (2.3), i.e. for any particle i = 1, …, 99 and time $\overline{t}$.

The dynamics of the predator are described by a first-order model as well, under the assumption that it is only subjected to the hunting stimulus. In mathematical terms, ${\mathcal{J}}_{100}=\{$predation} while, in the usual non-dimensional form, equations (2.1) and (2.2) can be rewritten as

where obviously ${\alpha }_{100}^{\mathrm{pred}}(\overline{t})=1$ for any $\overline{t}$. To specify the velocity contribution in equation (5.4), we can observe that, in general, the predator is attracted and consequently oriented by the group of prey: however, different hunting strategies can be identified. For instance, the confused predator does not have the ability to identify and attack a single individual within a set of agents [51,52]. In this respect, its chasing direction points towards the centre of mass of the population of prey, where all of them are equally targeted (if seen). From a modelling perspective, this amounts in defining: with where ${\overline{d}}_{\mathrm{pred}}$ is the non-dimensional extension of the predation region, that can be assimilated to the hunter visual depth, see also figure 2b.

In other cases (e.g. in other species of animals), the predator is not confused by the presence of multiple prey, it being able to target a specific individual within the group and to change strategy accordingly. For instance, it can chase the agent which, at a given time, is the closest to its position, i.e.

with

(a) Simulation details and numerical results

Multispecies dynamics are located within the planar space. The values of the non-dimensional parameters characterizing the behavioural rules of the group of prey remain unaltered with respect to the previous section, given their independency from both time and individual specificity, i.e. ${\overline{v}}_i=1$, ${\overline{d}}_{\mathrm{rep}}=1$, ${\overline{d}}_{\mathrm{align}}=8$ and ${\overline{d}}_{\mathrm{escape}}={\overline{d}}_{\mathrm{attr}}=30$ for any i = 1, …, 99 and $\overline{t}$. For the sake of simplicity, we also assume that the predator has the same physiological characteristics as its prey: we therefore fix ${\overline{d}}_{\mathrm{pred}}={\overline{d}}_{\mathrm{attr}}=30$ and ${\overline{v}}_{100}(\overline{t})=1$ for any $\overline{t}$.

The group of prey is initially arranged within the same round area of diameter equal to 6 used in the previous set of simulations, with random positions and velocity directions. The predator is initially placed at a significant distance from the cluster of individuals, however with at least a prey within its hunting region, i.e. ${\mathcal{N}}_{100}(\overline{t}=0)\ne \mathrm{\varnothing }$. Furthermore, the initial orientation of the hunter points to the centre of mass of the group of targets (figure 2c). The predator is finally assumed to catch a target when their relative non-dimensional distance drops below 10⁻².

The dynamics of the heterogenous system are then numerically studied upon variations of (i) predator hunting strategy and (ii) hierarchy of prey behavioural preferences which are quantified, as seen, by the weights α^(·), being $(\cdot )\in \{\text{rep, align, attr, escape}\}$. In particular, we hereafter focus on representative cases characterized by the fact that one (or more) prey behavioural stimuli significantly overcome the others. With the terminology ‘significantly overcome’, we arbitrarily mean that the corresponding α-value(s) is (are) at least twofold higher than the others. The resulting simulation outcomes are then analysed qualitatively (i.e. in terms of prey escape strategies) and quantitatively (i.e. in terms of time needed by the predator to eventually reach a target individual, hereafter defined as ${\overline{t}}_{\mathrm{p}}$). A numerical realization is finally stopped either as soon as the predator catches a prey or at $\overline{t}=2000$, i.e. at a sufficiently long time to have a clear idea of the system dynamics.

As shown in figure 5, when the prey agents are mainly subjected to the escape stimulus (i.e. ${\alpha }^{\mathrm{escape}}\gg {\alpha }^{\mathrm{rep}}=$${\alpha }^{\mathrm{align}}={\alpha }^{\mathrm{attr}}$), they are able to avoid the attack of both confused and not confused predators. In particular, they quickly move away from the approaching hunters without wasting time to organize in a well-spaced configuration or to align towards a preferred direction. In this respect, the group of prey randomly dissociates in disorganized colonies of different sizes with the predators falling within the empty space in between them and therefore being unable to reach any target, as shown also by the insets i1 and i2 in figure 5.

• Download figure
• Open in new tab
• Download PowerPoint

Differentiated phenomenologies instead emerge if the prey group is in an attraction regime (i.e.${\alpha }^{\mathrm{attr}}\gg {\alpha }^{\mathrm{rep}}={\alpha }^{\mathrm{align}}={\alpha }^{\mathrm{escape}}$). The confused predator falls and oscillates at the centre of a ring of radially escaping individuals, which still maintain a collective connection (see the left-middle panel in figure 5). The confused hunter is in fact unable to choose an optimal direction of attack, lying at the centre of mass of the set of prey locations. Such an evasion pattern has been numerically captured and analytically constructed and justified in a similar model [20]. The prey agents perform almost the same strategy in the case of a not confused predator, organizing in a half-moon of escaping agents. However, as shown in the centre panel in figure 5, the tendency to remain in an almost compact configuration delays their evasive manoeuvres: the hunter has therefore enough time to catch at least one target. In this respect, figure 7b(i) shows that the predation time ${\overline{t}}_{\mathrm{p}}$ decreases upon increments in the attraction stimulus of the prey (with respect to the other behavioural inputs), until stabilizing around a non-dimensional threshold value just below 4. Half-moon evasive patterns in the case of a predator pointing the centre of the target cluster has also been obtained in the work by Lee and colleagues [36] which will be reviewed in more details in the conclusive section.

A too large repulsion stimulus (i.e. ${\alpha }^{\mathrm{rep}}\gg {\alpha }^{\mathrm{align}}=$${\alpha }^{\mathrm{attr}}={\alpha }^{\mathrm{escape}}$) has a negative effect in the case of not confused predator: in order to maintain a comfort space, some prey may in fact turn back and go in the direction of the hunter, which is therefore facilitated in its purpose (see the left-bottom panel in figure 5). Increasing differences between prey preference for interagent avoidance and the other behavioural stimuli reduce the time needed by the predator to reach a target, as captured by the central plot in figure 7b. Also in this case a threshold value for ${\overline{t}}_{\mathrm{p}}$ emerges, which is close to 3. The confused predator is instead not able to take advantage of this prey phenomenology, as captured by the inset i3 in figure 5.

The prey group safely moves away from both types of predator when they are able to synchronize their movement, regardless of their spacing and escape stimulus (i.e. ${\alpha }^{\mathrm{align}}\gg {\alpha }^{\mathrm{rep}}={\alpha }^{\mathrm{attr}}={\alpha }^{\mathrm{escape}}$), see top panels in figure 6. In fact, as soon as some individuals within the population are able to perceive the presence and the location of the predator, they start to evade in the opposite direction: such information is then quickly transmitted to the rest of group and allows all animals to behave accordingly. Such a collective phenomenology represents an example of coordinated defence mechanism.

• Download figure
• Open in new tab
• Download PowerPoint

A coupling between sufficiently high escape and attraction stimuli (i.e. ${\alpha }^{\mathrm{escape}}={\alpha }^{\mathrm{attr}}\gg {\alpha }^{\mathrm{rep}}={\alpha }^{\mathrm{align}}$) is detrimental for prey in the case of a not confused predator. The tendency to remain somewhat compact in fact delays the evasive strategy: in particular, some individuals fall behind the rest of the group and are unable to avoid the attack of the focused hunter. In such a parameter regime, increments of the difference between the pair of more relevant stimuli and the others lead to decrements of the predation time ${\overline{t}}_{\mathrm{p}}$, which finally stabilizes around a limit value close to 20, see figure 7b(iii). This threshold is much larger than the corresponding quantities obtained in the previous parameter setting: this is probably a consequence of the fact that in this case the prey individuals are also subjected to a substantially high preservation stimulus. Such a hierarchy of prey preferred behaviour does not instead represent a sufficient advantage for the confused predator, which is still unable to point out and chase any single agent.

• Download figure
• Open in new tab
• Download PowerPoint

The co-presence of significant alignment and escape stimuli (i.e.${\alpha }^{\mathrm{escape}}={\alpha }^{\mathrm{align}}\gg {\alpha }^{\mathrm{rep}}={\alpha }^{\mathrm{attr}}$) is instead sufficient for the set of prey to avoid both confused and not confused predators. In particular, the groupmates organize in more or less compact clusters, each of them having a preferred direction of evasion. In this respect, the confused hunter undergoes an almost uneffective Brownian motion, being unable to choose and follow a single subgroup of prey (see the inset i4 in figure 6). On the other hand, the not confused predator, after an initial crawling, opts to follow a single cluster of targets but it can not reach the component individuals, which are already too far from its position (see the inset i5 in figure 6).

Dispersion in small aggregates is also the prey’s successful strategy when mainly guided by repulsion and escape tendencies (i.e. ${\alpha }^{\mathrm{escape}}={\alpha }^{\mathrm{rep}}\gg {\alpha }^{\mathrm{align}}={\alpha }^{\mathrm{attr}}$). It is however useful to notice that, in this case, prey agents belonging to the same (well-spaced) cluster can undergo uncorrelated movement, as captured by the insets i6 and i7 in figure 6.

Summing up, we can conclude that confusion completely drops predator ability to hunt its targets, regardless of their strategy, as also experimentally confirmed in [44,52,53]. From a prey perspective, alignment represents a significant strategic advantage. On the other hand, excessive grouping stimuli are detrimental because they make it easier for the predator to spot and attack prey, which proceed slowly to remain in visual-contact with the mates. The tendency to maintain a comfort spacing is instead negative if not accompained by a similar or a larger (in term of intensity) escape strategy. We finally remark, for the sake of completeness, that prey and predator have the same speed. Of course, a different hypothesis in this respect may result in variations of the simulation results.

6. Conclusive remarks

Populations of intelligent living entities are mainly characterized by the fact that the component agents are not passively prone to external forces but rather undergo active decision-based processes according to individual behavioural preferences and mutual interactions.

In this respect, the description of the collective organization and motion of animal groups has become of increasing interest also for the modelling community and treated with different techniques (as specified in the Introduction).

This topic has been here addressed with a microscopic particle approach able to distinguish speed and orientation. In particular, the latter has been assumed to be the result of competing stimuli, each of them simply weighted by a coefficient (α^(·)) that defines a sort of individual preference. The sum of the α-coefficients has been fixed equal to one in order to account for a balance between individual movement traits while avoiding their simultaneous minimization/maximization.

In the case of a homogeneous population of animals, a numerical analysis of the model has then allowed us to point out the hierarchy of behavioural inputs (i.e. attraction, repulsion and movement synchronization) eventually resulting in a large-time non-collapsed configuration of the system and/or in fully aligned dynamics. An analytical study on the H-stability properties of the animal system subjected only to the interaction velocity contributions has then been provided: it has been able to give a condition on the attraction/repulsion coefficients that a priori assures an asymptotic well-spaced particle pattern. Such a mathematical analysis has also been supported by a proper series of simulations: its application to the case of interaction kernels/potentials characterized by parametric constraints (i.e. here given by the unitary sum of the α-weights) is a quite novel aspect of this work.

Our approach has then been enriched to include the effect on collective dynamics of the presence of predators. In this respect, a confused hunter, being unable to single out and chase a preferred target, has been shown to constantly fail in its purpose. Such simulation results are in agreement with experimental studies [51–53] and have been also captured by similar mathematical approaches dealing with schools of fishes [25] and swarming prey [20,54–56]. On the opposite side, the search of an aligned movement has been shown to represent an escape advantage for the prey group in the case of the not confused predator. Evasive manoeuvres are instead delayed and impeded by individual preference towards cohesion. Such model outcomes are in partial disagreement with the experimental results obtained by the authors in the previously cited work [40]. Their empirical evidence in fact shows that the predation risk is reduced in the case of prey exhibiting a balance between attraction and orientation. A possible explanation of this discrepancy is owing to the fact that in our model the virtual prey individuals react to the presence of the predator thereby giving rise to somehow different dynamics. In [40], it is also observed that the predator preferentially attacks small groups of targets and individuals located at the edge of the prey aggregate, an aspect not captured by our approach.

It is however useful to remark that the model outcomes relative to predator-prey dynamics have here been obtained under the assumption that the hunter has the same physical characteristics as its targets. For instance, a common speed implies that a predator is not able to catch a prey if both follow a similar path starting from different positions. Interactions within heterogenous systems have been recently analysed in several other contexts, i.e. not strictly related to the animal world, including crowd dynamics, in particular in the perspective of groups with leaders [2], and cell migration, e.g. in the case of the coexistence of different cell lineages and phenotypes [57,58].

We here remark that, in principle, the computational results presented in the previous parts of the paper depend on the initial condition of the particle system. However, as it is possible to observe by the error bars in the plots of figure 7, the simulation outcomes obtained by keeping constant the initial density of agents (i.e. the number of animals and diameter of the round region where they are initially located) while randomly varying individual position and orientation do not significantly differ in almost all cases.

Keeping fixed its basic structure, our model can be easily extended to include a possibly large spectrum of individual behavioural rules: this would simply amount in adding the proper contribution in equation (2.2). It would be also possible to deal with a variability in the weights α, which would introduce adaptation aspects. For instance, a dispersed population of individuals may first try to group and then synchronize their movement. This would amount to give a time-dependent decrement to ${\alpha }^{\mathrm{attr}}$ (after an initially high estimate) and a simultaneous increment to ${\alpha }^{\mathrm{align}}$ (and to ${\alpha }^{\mathrm{rep}}$ to avoid particle overlap). Variations in the values of the α-parameters may be also related to specific system dynamics. For instance, the prey’s escape stimulus, quantified by coefficient ${\alpha }^{\mathrm{escape}}$, may increase when the predator is sufficiently close, while it may decrease becoming negligible when the hunter is substantially far or eventually absent. In this respect, interesting data can be obtained from experimental works. For example, Herbert and co-workers in [59] demonstrate that, in the case of shoaling fishes, the presence of a hunter increases the cohesion within the group, while decreasing dispersion and interindividual spacing. Further, it seems that the tendency and the intensity of the alignment processes are not affected. To employ these results in our model, we would have to keep fixed ${\alpha }^{\mathrm{align}}$, while simultaneously decrease ${\alpha }^{\mathrm{rep}}$ and increase ${\alpha }^{\mathrm{attr}}$.

Of course, in nature, predators of different species may have more or less sophisticated hunting strategies, as they can have the ability to target a single prey according to a given characteristic (e.g. age or physical limitations) and not only according to its distance. This aspect can be introduced in the proposed modelling framework by defining a state variable that, for each agent, describes the characteristic of interest, thereby introducing a differentiation within the set of prey. Finally, the application to specific groups of animals is straightforward, because it would only require the inclusion of proper empirical information on individual speed, extension of the interaction regions and relative behavioural preferences (i.e. hierarchies of α values).

Data accessibility

This article has no additional data.

Authors' contributions

All authors equally contributed to the paper.

Competing interests

We declare we have no competing interests.

Funding

This research was partially supported by the Italian Ministry of Education, University and Research (MIUR) through the ‘Dipartimenti di Eccellenza’ Programme (2018-2022)–Department of Mathematical Sciences ‘G. L. Lagrange’, Politecnico di Torino (CUP: E11G18000350001). Both authors are members of GNFM (Gruppo Nazionale per la Fisica Matematica) of INdAM (Istituto Nazionale di Alta Matematica), Italy.

Footnotes

One contribution of 14 to a theme issue ‘Multi-scale analysis and modelling of collective migration in biological systems’.