A Predator-Prey Model from Collective Dynamics and Self-Propelled Particles Approach

: The deﬁnition and description of the dynamics of predator-prey system consists in one of the fundamental problems of the population biology. Since 1925, several models have been introduced. Although they are highly effective, most of them neglect certain relevant criteria such as the spatial and temporal distribution of the studied species. We introduced these criteria by coupling two models designed initially for collective dynamics. The ﬁrst is for predators mobility, a Vicsek type model and the second is a Brownian Particle (BP) model for prey. We observed, as naturally in the classical models, periodic cycles of the density of predators and prey. In this case, the period of oscillations depends relatively on the collective dynamics parameters.


Introduction
The study of population dynamics is a subject of multidisciplinary research involving ecology, biology, mathematics and physics. It's sometimes oriented specifically in the dynamics of two species: Predators and Prey. In this perspective, the model developed independently by A. J. Lotka (1925) [1] and V. Volterra (1926) [2] provides a satisfactory explanation of the observed phenomena such as the periodic variation in species density.
Since then, several variants of this model (called LV) are derived. They are discrete or continuous dynamical systems whose qualitative study most often shows the existence of a steady state and periodic solutions around it [3,4]. Some variants of these models incorporate the mobility (spatio-temporal distribution) of the studied species, they are known as reaction-diffusion systems [5]. The common point between the above-mentioned approaches is that they are Eulerian and deterministic and known as the mean-field approach. There is still some stochastic formulations studied in [6,7].
Recently, Agent Based Model (ABM) are also used to interpret predator-prey dynamics [8]. This concept was first introduced to study the collective behavior of a given species. These are the Self Propelled Particles Models (SPP) ; they are presented in two or threedimensional space by an extensive investigation [9][10][11][12][13]. Central issues here include the mechanisms by which autonomous agents interact to exhibit emergent collective behavior and the properties of the resulting behavior. On the other hand, the mobility of certain species such as bacteria or phytoplankton is not a coherent collective dynamic but rather a random walk. Individuals from these species are generally considered as Brownian Particles (BP) [14][15][16]. In sum, Self-Propelled Particles Models are introduced to study the collective dynamics of a single species but never considered for the simultaneous description of two species ; Less so for two species of Predator-Prey. If so, what about the basic properties deduced from these models such as orientation ordering, cluster formation etc.?
Our approach consists in constructing a framework interpreting a Predator-Prey's type of dynamics. Individuals are considered as material points inside the domain: a periodic box of size B × B located by their positions on the Cartesian plane and animated with some velocity.
To implement this approach, we develop in the first section two separate models to simulate mobility of prey and predators. For this purpose, we use the SPP and BP models respectively. Then, we mix the two with assumptions of living interactions to illustrate via flowchart and equations a prey-predator type dynamics. We make then a python program to incorporate the global model.
Throughout the analysis of this model, we look for a range of standard parameters for which periodic oscillations between the densities of species can be found. This is followed by an underlying analysis of the characteristics of a collective behavior.

Models and Simulation Setting
We designed separately prey and predator model describing their motions. Then we added Predator-Prey model which take in account their living interactions.

Prey's Motion
For prey motion, n p particles are randomly uniformly distributed inside the domain, their density by unit area is ρ p = n p /B 2 . r i p (t), v i p (t) represent respectively, the position and velocity of particle i at time t where the index p is set for prey. A complex type of this models using thrust force and interactions between particles is studied in [16].
In the beginning of the simulation (t = 0), we set prey's velocities from a Gaussian distribution with zero mean and σ standard deviation (i.e., each particle has a random where, M is the mass of prey, is set to 1 for simplicity ; v i p the velocity of prey i; γ is the coefficient of viscous friction, which is set by the properties of the environment and the particle geometry; R i (t) is a random force given by: R i (t) = 2γk B Tζ i , with k B , T and ζ i are respectively the Boltzmann constant, temperature and Gaussian white noise with zeromean and unit variance. T B = k B T is called Boltzmann temperature. We used Box-Muller transformation to compute the components of ζ i , it consists in choosing two elements u 1 and u 2 from a uniform distribution from ]0, 1] and mapping them to two standard, normally distributed samples: ζ i x = −2 log(u 1 ) cos(2πu 2 ), ζ i y = −2 log(u 1 ) sin(2πu 2 ). Thus, we used Verlet algorithm to perform prey's position and Störmer-Verlet method for their velocities. (Cf. [17]) r i In this system, each particle has a Brownian trajectory (random walk) and the average velocity of all particles depends closely on two parameters γ and T B (Cf. additional materials).

Predator's Motion
To describe predator's motion in the domain, we used a variant of Vicsek model and added a repulsion rule (Cf. [9,12,18]). So, n q predator-particles are distributed inside the domain and move with a constant speed v o q ; Their density is ρ q = n q /B 2 and r i q (t), v i q (t) represent respectively position and velocity vectors for i th particle at time t where subscript q denote predator. For initial conditions, we set to each particle, a uniform random position r i q (0) and velocity r i q (0) (Cf. [16]). Then, in any time step τ, the position r i q (t) changes according to Equation r i q (t + τ) = r i q (t) + τv i q (t), and the velocity v i q depends on the position of the i th particle which is either in the area of repulsion (zor) or alignment (zoa). These areas are defined by circles whose radii are noted R r for zor and R a for zoa with R r R a . The first behavioral zone, represented by a circle of smaller radius R r , is responsible for the maintenance of a minimum distance between neighboring particles, while the second one with a grater radius R a is used for aligning the particle velocities ( Figure 1a). In the same figure, we put flowchart illustrating movement choices of a predator i whose rules are described as follows: (a) Repulsion rule: If there are n r neighbors around i in zor, the velocity v i where v o q is a constant speed for all particles; r ij q = r j q − r i q is the vector in the direction of neighbor j and R(ξ i (t)) is the rotation matrix. ξ i (t) is a random angular noise obtained from [−2πη; 2πη] by Gaussian-distribution, where η ∈ [0; 1] represents the noise strength. The standard deviation is then around 2.π 3 η. (b) Alignment rule: If there are no neighbors in the zone of repulsion, the individual i responds to others which are in the zone of alignment with: where: b is the cosine of the angle measured between v i q (t) and r ij q ; It's evaluated to check if particle i see particle j ; a ∈ [−1; 0] has a fixed value which is the cosine of half angle of perception α, the blind angle should be 360 • − α. Thus, if a = −1, central particle i interact with all particles inside zone of alignment. The desired direction of particle i in with n a is the number of particles in the zone of alignment and φ(t) the polar order (or polarization) given by: φ(t) = 1 n a n a ∑ i=1 v i q (t) ; and λ is the degree of polarization.
(c) Absolute rule: If repulsion happens, the alignment is neglected .
(d) Alternative rule: If there are no individuals in the both zones, the particle i keeps its previous direction thus, v i q (t + τ) = v i q (t)v o q . Finally, the parameters involved in this model for predators are shown in Table 2. Thus, in the next section, we will model the vital interactions simulating a prey-predator dynamic.

The Predator-Prey Model
In addition to previous models describing movement of prey and predator (Sections 2.1 and 2.2), we interested in this part on their dynamic of evolution as two species living in a common environment (the domain). To describe this dynamic, we introduced a Predator-Prey model by coupling those two models through underline assumptions.
Firstly, we treated prey as living agents with the principle idea that any prey has some energy state which varies according to its position towards predators. It features life states of all individuals: survival, reproduction and mortality. Then, at t = 0, we set n p (0) the initial size of prey population and an initial energy E i p (0) to each one. It's a random value gotten from Gaussian distribution with zero mean and σ E p standard deviation (i.e., E i p (0) → ℵ(0, σ E p )). Therefore, after each time step, i th prey survive when it's outside the predatory area of all predators ( Equation ??). It gains a survival bonus, an energy E s p which can increase or decrease its next state of energy, this modeling the fact that each prey can survive by engendering or not another one or just die. ∀ j |r Consequently, if the energy state reach some certain reproduction threshold E r p , the individual will be able to reproduce by cloning. Then a new prey k born inside a birth zone (zob p ) of radius R b p ; he inherits from his ancestor the similar movement abilities E i p (t + τ) ≥ E r p ⇒ i th prey clone. Subsequently, its mortality is due to predatory or deficiency of energy, it happens when the energy become less than a fixed mortality threshold E d p , E i p (t) ≤ E d p ⇒ i th prey die. With the same order of ideas, we approached in second instance the dynamics of predators. Thus, for an initial size of population n q (0), we assumed that predator will need energy to survive, it's also a normal distribution E i q (0) with zero mean and σ E q standard deviation. Therefore, in each time step, it uses E m q amount of energy for metabolism (2). It's deriving energy come from predation of preys present in it's zone of predatory zop (a circle of radius R p ). Thus, its state of energy decrease by the metabolism cost E m q or increase by the predatory bonus E a q , Following this energy update, predator dies when it's state of energy decrease below the mortality threshold E d q . Otherwise, if its energy level reaches more than the reproduction threshold E r q , a new predator k born with a random normal energy E k q (t + τ) → ℵ(0, σ E q ) and random uniform velocity inside the zone of birth (zob q ) of radius R b q . This is the cloning reproduction principle that we used for prey E i q (t) ≤ E d q ⇒ i th predator die, E i q (t + τ) ≥ E r q ⇒ i th predator clone. Finally, after recap on Table 3, the parameters involved in the model described above, we designed a python program for simulations. The files are available on this link where "main.py" generates the data and "anim.py" simulates interactions between the species.

Results And Discussion
We performed computer simulations to analyze the predator-prey model. Unless explicitly stated, we kept all parameters standard as indicated in Table 1. Standard parameters provide a regime where the density of species oscillates over time,a pseudopereriodic regime for all trials done. However, by varying certain parameters, one can have other types of regime like the simultaneous extinction of two species or the survival of one and the extinction of the other.
We started to run ones the model by focusing mainly on three characteristics time in [ 3 × 10 3 , 8 × 10 3 ]. The choice of characteristic time is with respect to predator density as: (ρ q min ← t min ), (ρ q mid ← t mid ), (ρ q max ← t max ) where subscripts min, mid, and max denote respectively the minimum, middle and maximum densities. Not that in this notation t min might be greater than t max etc., because subscripts here are only rigorously attached to density values. These three cases typically represent the different phases of predators collective dynamics. In this sense, we illustrated the spatial distribution of species on Figure 2A snapshots (a-c); i.e., the locations inside the domain of predators (blue) and preys (orange) at specific time (t min , t mid , or t max ). We noticed that prey density is less important where that of predators is elevated and vice versa, such a fact gives the intuition to cyclic evolution of species which is represented in Figure 2A(d). For hole experience, we recorded the range of densities for predators and prey respectively in [0.089, 0.25] and [0.29, 0.89]. The corresponding characteristics densities are (ρ q min , ρ q mid , ρ q max ) = (0.089, 0.21, 0.25). Hence, we observed that density variations have an impact on the collective dynamics notably on the particle orientation ordering and cluster formation groups ; this results is also found in many study of Viscek-type models.
For this, we investigated particle orientation ordering by evaluating polar order parameter Φ(t), which quantifies the alignment of predator particles to the average instantaneous velocity vector Φ(t) = 1 ; this parameter goes from zero in the case of a disordered phase to one in the case of a completely ordered phase. We have shown on Figure 2A(e) the evolution of polar order parameter Φ(t) where particles seem less orderly at low density. Anyway, it's not obvious to find a correlation between the variation of particles density and their polar order parameter because the appearance of newborns disturbs alignment of particles, and we observe emergence of a spontaneous local phase transition. To more understand the evolution of predator-prey densities (ρ q (t), ρ p (t)) and the equivalent predator polar ordering (Φ(t)), we run 18 trials with total time L t = 5 × 10 4 . Results are like a same i.e., the density of predators ρ q (t) oscillate in [0.013; 0.33] and that of preys ρ p (t) in [0.12; 1.32]. It results as on the first trial that ρ q (t) increase when ρ p (t) decrease and vice versa. Thus, there is an equilibrium point (ρ q (t),ρ p (t)) and limit cycle around it. We also highlighted that the curve of ρ q is smoother than the curve of ρ p , this is due to the height impact of randomness in prey dynamic assumptions. This randomness effect appears also on polar order parameter Φ(t) which is bounded in [0.076, 0.67]. Therefore, predators could be completely disordered some time, and they cannot reach a degree of order more than 0.67. This indicates that the presence of preys has an impact on predator collective dynamic as disordering factor.
Since only the evolution of the predator density does not follow a Gaussian distribution, we interested on fluctuation in the number of predator particles in the system: ∆n q (l) = n q (l) 2 − n q (l) 2 where n q (l) denotes the number of predators in a box of linear size l . ∆n q (l) is in general proportional to n q (l) β . β = 0.5 correspond to normal fluctuation while β0.5 is the giant fluctuation. In our system, the number fluctuation of predator particle is estimated for three density case, for minimum and middle densities, the critical exponent is β 0.8 ; this corresponds to giant fluctuation. For maximum density, the fluctuation is normal because of β 0.5 (Figure 3a). We examined finally cluster distribution of predator particle. Formally, the cluster dynamics of the SPP system can be described by deriving a master equation for the evolution of the probability p(m), where m = m 1 , m 2 , ..., m N with m 1 being the number of isolated particles, m 2 the number of two-particle clusters, m 3 the number of three-particle clusters etc. Cluster in our model is defined as a group of particles with a distance between neighbors smaller or equal to the radius of alignment zone R a , i.e., particles interacting directly or via neighboring agents are included into one cluster ( [13]).  ρ q ≤ ρ q min ,ρ q min ≤ ρ q ≤ ρ q max ,ρ q ≥ ρ q max . ρ q min = 0.16 and ρ q max = 0.21. L t = 10 5 and the rest of parameters are standard.

Conclusions
In summary, we built a predator-prey model from collective dynamics and selfpropelled particles approach. Despite the plethora of assumptions, the stability of the system is possible with the admissible parameters which gives a quasi-periodic regime.
The main results related to SPP & BP models are kept with the particularity that the densities are dynamic here. Let us note however the difficulty in the design of the algorithm caused by the large size of parameters.