#  Efficient ensemble stochastic algorithms for agent-based models with
spatial predator-prey dynamics

Giacomo Albi, Roberto Chignola and Federica Ferrarese

###### Abstract.

Experiments in predator-prey systems show the emergence of long-term cycles.
Deterministic model typically fails in capturing these behaviors, which emerge
from the microscopic interplay of individual based dynamics and stochastic
effects. However, simulating stochastic individual based models can be
extremely demanding, especially when the sample size is large. Hence, we
propose an alternative simulation approach, whose computation cost is lower
than the one of the classic stochastic algorithms. First, we describe the
agent-based model with predator-prey dynamics, and its mean-field
approximation.Then, we provide a consistency result for the novel stochastic
algorithm at the microscopic and mesoscopic scale. Finally, we perform
different numerical experiments in order to test the efficiency of the
proposed algorithm, focusing also on the analysis of the different nature of
oscillations between mean-field and stochastic simulations.

Giacomo Albi: Department of Computer Science, University of Verona, e-mail:
<EMAIL_ADDRESS>
Roberto Chignola: Department of Biotechnology, University of Verona, e-mail:
<EMAIL_ADDRESS>
Federica Ferrarese: Department of Mathematics, University of Trento, e-mail:
<EMAIL_ADDRESS>
##  1\. Introduction

The description of how biological populations interact and evolve in both
space and time is central to theoretical ecology, [14, 31, 25, 32, 28, 4]. To
this aim, several mathematical models have been proposed to approximate the
evolution of large ensembles of interacting species, starting from microscopic
stochastic agent-based dynamics, passing from mesoscopic equations, up to
macroscopic systems, such as reaction-diffusion systems. This hierarchy of
descriptions, at different scales, has been used extensively to model
interactions among species and their spatial spread beyond theoretical
ecology, and it has a profound impact in several fields such as epidemiology,
[6, 38, 2], biology, physiology and medicine, [8, 34, 33], chemical reactions
[20, 11], and in socio-economic systems opinions, markets, veichular traffic
and crowds [30, 10, 1]. In the present manuscript, we will focus on the
efficient simulation of predator-prey dynamics in population biology, such as
[28, 9, 25]. One striking feature of population biology is the emergence of
temporal cycles of species densities, [28, 29, 3]. Various mechanisms can
cause populations to oscillate, the most investigated being the cyclic
dynamics that arise from the interactions between predators and preys (see
ref. [5] and references cited therein). Traditionally, these interactions have
been studied in the context of deterministic Lotka-Volterra-type models, [28,
29, 3]. In their simplest form, these are coupled deterministic equations that
include specific birth, death, competition and predation processes, and they
can show limit cycles behaviors for appropriate choice of parameters values.
Long-term persistence of predator-prey cycles have been recently observed in a
series of well controlled microcosm experiments with freshwater organisms,
[5]. Planktonic rotifers cultured together with their prey unicellular green
algae showed oscillatory behaviors of both prey and predator densities that
persisted for up to approximately 50 cycles or $\sim$300 generations, [5].
Interestingly, dominant dynamics characterized by coherent oscillations were
interrupted by shorter episodes of irregular non-coherent oscillations, [5].
The experiments clearly demonstrate that sustained oscillations in population
dynamics do arise even in simple well-controlled ecosystems. They strongly
indicate that stochastic forces are at work to drive the reversible shift from
coherent to non-coherent oscillatory behaviours of observed populations, [5],
but their role in driving population dynamics can not be investigated within
the context of deterministic predator-prey models.

This novel experimental evidence, therefore, calls for new theoretical
explanations, and we note that these might be found within the theoretical
framework developed by McKane and Newman (see [26, 27]). They studied
individual based predator-prey models described by simple stochastic processes
of mortality, reproduction and predation, and showed that large cycles of
species densities persisted unless the number of individuals was taken to be
strictly infinite, [26, 27]. Biological cycles have been found to arise from a
novel resonance effect that include the statistical fluctuations of a given
finite-size population, an effect that fade away when different realization of
the model are averaged to recover the mean-field behavior described by the
classic deterministic predator-prey equations, [26, 27].

Finding the solution of stochastic individual-based models, however, can
easily become computationally challenging, above all when the final goal is to
investigate the dynamics of large, albeit finite, population sizes over
hundreds of species generations. Moreover, depending on the context of
application, ad-hoc methods are of paramount importance to capture the
essential features of phenomena at various scale, see for example [13, 21, 12,
35].

In the present work, to cope with the increasing complexity of the model we
rely on Monte Carlo methods, introducing a novel stochastic algorithm for the
simulation of the population dynamics with individual based models. We
investigate the efficiency of the new method on spatial predator-prey models
and we compare the computational costs, also by exploring in part the space of
parameters, to that of different stochastic algorithms such as the classic
Monte Carlo, the direct method or the $\tau-$leaping approach considered as
benchmarks [23, 15, 7].

The fundamental idea behind our work is that with classic approaches the whole
population sample must be reconstructed any time interaction between
individuals occurs with the updated number of individuals. However, we show
that there exists a maximum number of interactions that can take place before
updating the sample, and this allows us to reduce the total number of time
steps required for the whole simulation. As the consequence, the new algorithm
has a computational cost that is up to 2 orders of magnitude less than those
of the classic stochastic algorithms considered for a wide range of parameter
values. The simulations obtained with the new algorithm converge to the mean-
field solutions for increasing sample size $N$ with error $\propto
1/\sqrt{N}$, as also observed for classic algorithms, while individual
realizations for finite population sizes oscillates as observed by McKane and
Newman.

In Section 2 we present an individual-based stochastic model where predators
compete with preys and migration of individuals is allowed. We derive the
master equation and the mean-field equations which are obtained when the
number of individuals is taken to be strictly infinite. In Section 3 we
present a novel ensemble stochastic algorithm and we prove its consistency
with the classical formulation. In Section 4 we show different numerical
experiments aimed at testing the new algorithm and at comparing it with the
classical approaches. We also partly investigate the stochastic persistency of
the oscillatory behavior shown by both predator and prey dynamics. Overall our
approach can be exploited to efficiently simulate stochastic agent-based
models and thus to explore the emergence of long-lasting persistent resonant
effects in population dynamics. In A we show how to reduce the model to the
homogeneous case, assuming that individuals are not allowed to migrate. In B
we briefly report some of the classical exact and approximated stochastic
algorithms and we describe them in the homogeneous case.

##  2\. Agent-based models with predator-prey interactions

We consider an agent based stochastic model of predation between two species
with migration. We will show that for considerably large number of individuals
the evolution of this model reduces to a system of partial differential
equations.

###  2.1. Spatially heterogeneous predator-prey model

We consider a stochastic process describing the evolution of individuals
distributed in the spatial domain $\Omega=[0,L]$ with $L>0$, divided in
$M_{c}$ spatial cells. We assume that every cell $C_{\ell}$, with
$\ell=1,\ldots,M_{c}$, has $N_{c}$ components, where each component is
realized in one of the following states: predator ($s_{A,\ell}$), prey
($s_{B,\ell}$), or empty $(s_{E,\ell})$. To describe the stochastic dynamics,
we consider cell $C_{\ell}$ and its nearest cells $C_{\ell\pm 1}$, for a fixed
time $t>0$ we define different interaction events.

  * •

Competition & Birth events. We sample a component in cell $C_{\ell}$ and with
probability $q_{1}\in[0,1]$ we let it interacts with another component,
randomly chosen among the other $N_{c}-1$ components in cell $C_{\ell}$,
according to the following birth and competition rules with rates
$b^{r},p^{r}_{1},p^{r}_{2}$

(1)
$\begin{split}&s_{B,\ell}s_{E,\ell}\xrightarrow{b^{r}}s_{B,\ell}s_{B,\ell},\quad
s_{A,\ell}s_{B,\ell}\xrightarrow{p^{r}_{1}}s_{A,\ell}s_{A,\ell},\quad
s_{A,\ell}s_{B,\ell}\xrightarrow{p^{r}_{2}}s_{A,\ell}s_{E\ell},\end{split}$

where $b^{r},p^{r}_{1},p^{r}_{2}>0$ are constant parameters. We assume that if
interactions occur among two empty components then no change in the
populations sizes is accounted.

  * •

Migration event. We sample one component in cell $C_{\ell}$ and we assume that
with probability $q_{2}\in[0,1]$ such that $q_{1}+q_{2}\leq 1$ it interacts
with another component, randomly chosen among the $N_{c}$ positions in one of
the nearest cells $C_{\ell\pm 1}$. Changes in the state happen according to
migration rates $m^{r}_{1},m^{r}_{2}$

(2) $\begin{split}&s_{A,\ell}s_{E,\ell\pm
1}\xrightarrow{m^{r}_{1}}s_{E,\ell}s_{A,\ell\pm 1},\quad
s_{B,\ell}s_{E,\ell\pm 1},\xrightarrow{m^{r}_{2}}s_{E,\ell}s_{B,\ell\pm 1},\\\
&s_{A,\ell\pm 1}s_{E,\ell}\xrightarrow{m^{r}_{1}}s_{E,\ell\pm
1}s_{A,\ell},\quad s_{B,\ell\pm
1}s_{E,\ell}\xrightarrow{m^{r}_{2}}s_{E,\ell\pm 1}s_{B,\ell},\end{split}$

where $m^{r}_{1},m^{r}_{2}>0$. If interactions occur among predator/prey and
empty components or among two empty components then no changes in the
populations sizes happen.

  * •

Death event. We sample a component in cell $C_{\ell}$, and we assume that with
probability $(1-q_{1}-q_{2})$ it changes according to death rates
$d^{r}_{1},d^{r}_{2}$

(3) $\begin{split}s_{A,\ell}\xrightarrow{d^{r}_{1}}s_{E,\ell},\quad
s_{B,\ell}\xrightarrow{d^{r}_{2}}s_{E,\ell},\end{split}$

where $d_{1}^{r},d_{2}^{r}>0$. We suppose that if the selected component is
empty then no change is accounted.

We assume that at most one birth/competition, one migration and one death
event can occur at each time step. In order to introduce a master equation for
such process, we consider the vector state ${\bf x}=({\bf A},{\bf B},{\bf
E})\in\mathbb{R}^{3M_{c}}$, where ${\bf A}=(A_{1},A_{2},\ldots,A_{M_{c}})$,
${\bf B}=(B_{1},B_{2},\ldots,B_{M_{c}})$, ${\bf
E}=(E_{1},E_{2},\ldots,E_{M_{c}})$ account respectively the number of
predators, prey and empty spaces in each cell. We first define the
stoichiometry matrices $\hat{V}$, associated to single competition/birth and
death events, and $\hat{V}_{M}$, associated to migration events

$\hat{V}=\begin{bmatrix}0&1&-1\\\ 1&-1&0\\\ 0&-1&1\\\ -1&0&1\\\
0&-1&1\end{bmatrix},\quad\hat{V}_{M}=\begin{bmatrix}1&0&-1\\\ -1&0&1\\\
0&1&-1\\\ 0&-1&1\\\ \end{bmatrix}.$

Hence, the stoichiometry matrix of the full process over the lattice of
$M_{c}$ cells is defined as follows,

(4) $\tilde{V}=\begin{bmatrix}\hat{V}\otimes I\\\ \hat{V}_{M}\otimes M_{-}\\\
\hat{V}_{M}\otimes M_{+}\end{bmatrix},$

where $\otimes$ denotes the Kronecker product, and $I$ is the identity matrix
of size $M_{c}\times M_{c}$, and $M_{-},M_{+}$ are squared are square matrices
of size $M_{c}\times M_{c}$ defined as

$M_{-}=\begin{bmatrix}0&0&\ldots&0&0\\\ 1&-1&0&&\vdots\\\
0&1&\ddots&\ddots&0\\\ \vdots&\ddots&\ddots&-1&0\\\
0&\ldots&0&1&-1\end{bmatrix},\quad M_{+}=\begin{bmatrix}-1&1&0&\ldots&0\\\
0&-1&1&\ddots&\vdots\\\ 0&0&\ddots&\ddots&0\\\ \vdots&\ddots&\ddots&-1&0\\\
0&\ldots&0&0&0\end{bmatrix}.$

Here, each row of the stoichiometry matrix (4) represents the changes in the
populations sizes due to the occurrence of the events described in (1)-(2)-(3)
and each column represents the predators, preys and empty components in each
cell. Then, for any $\ell=1,\ldots,M_{c}$, we write the associated transition
rates as follows

(5) $\begin{split}&\pi_{{\bf v}_{\ell_{0}}}({\bf
x})=2b^{r}q_{1}\frac{B_{\ell}}{N_{c}}\frac{E_{\ell}}{N_{c}-1},\qquad\pi_{{\bf
v}_{\ell_{1}}}({\bf
x})=2p_{1}^{r}q_{1}\frac{A_{\ell}}{N_{c}}\frac{B_{\ell}}{N_{c}-1},\cr&\pi_{{\bf
v}_{\ell_{2}}}({\bf
x})=2p_{2}^{r}q_{1}\frac{A_{\ell}}{N_{c}}\frac{B_{\ell}}{N_{c}-1},\qquad\pi_{{\bf
v}_{\ell_{3}}}({\bf
x})=d_{1}^{r}(1-q_{1}-q_{2})\frac{A_{\ell}}{N_{c}},\cr&\pi_{{\bf
v}_{\ell_{4}}}({\bf
x})=d_{2}^{r}(1-q_{1}-q_{2})\frac{B_{\ell}}{N_{c}},\quad\pi_{{\bf
v}_{\ell_{5}}}({\bf
x})=m_{1}^{r}q_{2}\frac{A_{\ell}}{N_{c}}\frac{E_{\ell-1}}{N_{c}},\cr&\pi_{{\bf
v}_{\ell_{6}}}({\bf
x})=m_{1}^{r}q_{2}\frac{A_{\ell-1}}{N_{c}}\frac{E_{\ell}}{N_{c}},\quad\qquad\pi_{{\bf
v}_{\ell_{7}}}({\bf
x})=m_{2}^{r}q_{2}\frac{B_{\ell}}{N_{c}}\frac{E_{\ell-1}}{N_{c}},\cr&\pi_{{\bf
v}_{\ell_{8}}}({\bf
x})=m_{2}^{r}q_{2}\frac{B_{\ell-1}}{N_{c}}\frac{E_{\ell}}{N_{c}},\quad\qquad\pi_{{\bf
v}_{\ell_{9}}}({\bf
x})=m_{1}^{r}q_{2}\frac{A_{\ell}}{N_{c}}\frac{E_{\ell+1}}{N_{c}},\cr&\pi_{{\bf
v}_{\ell_{10}}}({\bf
x})=m_{1}^{r}q_{2}\frac{A_{\ell+1}}{N_{c}}\frac{E_{\ell}}{N_{c}},\quad\qquad\pi_{{\bf
v}_{\ell_{11}}}({\bf
x})=m_{2}^{r}q_{2}\frac{B_{\ell}}{N_{c}}\frac{E_{\ell+1}}{N_{c}},\cr&\pi_{{\bf
v}_{\ell_{12}}}({\bf
x})=m_{2}^{r}q_{2}\frac{B_{\ell+1}}{N_{c}}\frac{E_{\ell}}{N_{c}},\end{split}$

where we define the operators $\pi_{{\bf v}_{\ell_{j}}}({\bf x})=\pi({\bf
x}+{\bf v}_{\ell_{j}}|{\bf x})$ for any $j=0,\ldots,M$. Here, $M+1=13$ is the
total number of events described in (1)-(2)-(3) and ${\bf v}_{\ell_{j}}$ is
the $\ell_{j}$-th row of the stoichiometry matrix $\tilde{V}$ defined in (4),
$\ell_{j}=\ell+jM_{c}$ for $j=0,\ldots,M$. We assume that migrations events
are not allowed at the boundaries. Hence we introduce
$A_{0},E_{0},B_{0},A_{M_{c}+1},E_{M_{c}+1},B_{M_{c}+1}$ to be equal to zero,
in order to properly define the transition rates associated to the boundary
cells. The density $P({\bf x},t)$ describing the probability of the state
${\bf x}$ evolves according to the master equation as follows

(6) $\frac{dP({\bf x},t)}{dt}=\sum_{\ell_{j}\in\mathcal{J}}\Biggl{[}\pi({\bf
x}|{\bf x}-{\bf v}_{\ell_{j}})P({\bf x}-{\bf v}_{\ell_{j}},t)-\pi({\bf x}+{\bf
v}_{\ell_{j}}|{\bf x})P({\bf x},t)\Biggr{]},$

where

(7) $\mathcal{J}=\\{\ell+jM_{c}|j=0,\ldots,M\text{ and
}\ell=1,\ldots,M_{c}\\}.$

######  Remark 1.

  * •

We assume that each spatial cell contains at most $N_{c}$ components that can
be either preys or predators or empty spaces. Hence, even if predators and
preys are restricted to move only to empty neighbor sites, the model can
include regions occupied both by predators and preys. In the case of a lattice
consisting of a single spatial cell, we refer to spatial homogeneous model,
where migration events are neglected and maximal density per cell is
$N_{c}=N$. We report in A the corresponding agent-based dynamics, and the
associated master equation (6).

  * •

We observe that the vector state ${\bf x}$ can be dimensionaly reduced if we
assume to remove the empty components, considering the relation ${\bf
E}=N_{c}-{\bf A}-{\bf B}$. However, we prefer to give a more general
description, to allow for further generalization, such as models that may
include empty sites with different nature.

###  2.2. Mean-field approximation

In order to recover the mean-field behavior of the stochastic process, we
introduce the empirical densities $f_{\ell}^{N_{c}}(t),g_{\ell}^{N_{c}}(t)$,
for $\ell=1,\ldots,M_{c}$, as the averaged quantities

(8) $\displaystyle f_{\ell}^{N_{c}}(t)=\frac{\left\langle
A_{\ell}\right\rangle}{N_{c}}=\frac{1}{N_{c}}\sum_{A_{\ell}=0}^{N_{c}}\sum_{B_{\ell}=0}^{N_{c}}\sum_{E_{\ell}=0}^{N_{c}}A_{\ell}P({\bf
x},t),$ $\displaystyle g_{\ell}^{N_{c}}(t)=\frac{\left\langle
B_{\ell}\right\rangle}{N_{c}}=\frac{1}{N_{c}}\sum_{A_{\ell}=0}^{N_{c}}\sum_{B_{\ell}=0}^{N_{c}}\sum_{E_{\ell}=0}^{N_{c}}B_{\ell}P({\bf
x},t),$

where $\langle\cdot\rangle$ denotes the expected value, and the empirical
density for the empty component can be recovered in each cell computing
$1-f_{\ell}^{N_{c}}(t)-g_{\ell}^{N_{c}}(t)$.

Hence, multiplying the master equation by $A_{\ell}$ and by $B_{\ell}$, and
summing over all the values of $A_{\ell}$, $B_{\ell}$ and $E_{\ell}$ for any
$\ell=1\,\ldots,M_{c}$ leads to the following Proposition.

######  Proposition 1.

By standard assumptions of the mean-field limit we assume that for any
$\ell=1\,\ldots,M_{c}$, $\left\langle
A_{\ell}B_{\ell}\right\rangle=\left\langle A_{\ell}\right\rangle\left\langle
B_{\ell}\right\rangle$, $\left\langle A_{\ell}^{2}\right\rangle=\left\langle
A_{\ell}\right\rangle^{2}$, $\left\langle
B_{\ell}^{2}\right\rangle=\left\langle B_{\ell}\right\rangle^{2}$,
$\left\langle A_{\ell\pm 1}B_{\ell}\right\rangle=\left\langle A_{\ell\pm
1}\right\rangle\left\langle B_{\ell}\right\rangle$ and $\left\langle
A_{\ell}B_{\ell\pm 1}\right\rangle=\left\langle
A_{\ell}\right\rangle\left\langle B_{\ell\pm 1}\right\rangle$ then the time
evolution of the empirical population densities $f_{\ell}^{N_{c}}$,
$g_{\ell}^{N_{c}}$ is given by

(9)
$\begin{split}\frac{df_{\ell}^{N_{c}}}{d\tau}=&2\tilde{p^{r}_{1}}\frac{\langle
A_{\ell}\rangle}{N_{c}}\frac{\langle
B_{\ell}\rangle}{N_{c}-1}-\tilde{d^{r}_{1}}\frac{\langle
A_{\ell}\rangle}{N_{c}}\cr&\qquad\qquad\qquad+\tilde{m}^{r}_{1}\left(\frac{\Delta_{\epsilon}\langle
A_{\ell}\rangle}{N_{c}}+\frac{\langle
A_{\ell}\rangle}{N_{c}}\frac{\Delta_{\epsilon}\langle
B_{\ell}\rangle}{N_{c}}-\frac{\langle
B_{\ell}\rangle}{N_{c}}\frac{\Delta_{\epsilon}\langle
A_{\ell}\rangle}{N_{c}}\right),\\\
\frac{dg_{\ell}^{N_{c}}}{d\tau}=&r\frac{\langle
B_{\ell}\rangle}{N_{c}}\left(1-\frac{\langle
B_{\ell}\rangle}{qN_{c}}\right)-\alpha\frac{\langle
A_{\ell}\rangle}{N_{c}}\frac{\langle
B_{\ell}\rangle}{N_{c}}\cr&\qquad\qquad\qquad+\tilde{m}^{r}_{2}\left(\frac{\Delta_{\epsilon}\langle
B_{\ell}\rangle}{N_{c}}+\frac{\langle
B_{\ell}\rangle}{N_{c}}\frac{\Delta_{\epsilon}\langle
A_{\ell}\rangle}{N_{c}}-\frac{\langle
A_{\ell}\rangle}{N_{c}}\frac{\Delta_{\epsilon}\langle
B_{\ell}\rangle}{N_{c}}\right),\end{split}$

where

$\Delta_{\epsilon}h_{\ell}=\sum_{s\in\\{\ell-1,\ell+1\\}}\frac{h_{s}-h_{\ell}}{\epsilon^{2}},$

for any function $h$ is the discrete Laplace operator, $\epsilon$ is the
lattice spacing, $\tau=t/N_{c}$, the parameters are defined according to the
following scaling

(10)
$\begin{split}&\tilde{b}^{r}=b^{r}q_{1}\quad\tilde{p_{1}}^{r}=p^{r}_{1}q_{1}\quad\tilde{p_{2}}^{r}=p^{r}_{2}q_{1}\quad\tilde{d^{r}_{1}}=(1-q_{1}-q_{2})d^{r}_{1}\\\
&\tilde{d_{2}}^{r}=(1-q_{1}-q_{2})d^{r}_{2}\quad\tilde{m}^{r}_{1}=q_{2}\epsilon^{-2}m^{r}_{1}\quad\tilde{m}^{r}_{2}=q_{2}\epsilon^{-2}m^{r}_{2},\end{split}$

and

$r=2\tilde{b}^{r}-\tilde{d_{2}}^{r},\qquad
q=1-\frac{\tilde{d_{2}}^{r}}{2\tilde{b}^{r}},\qquad\alpha=2(\tilde{p}^{r}_{1}+\tilde{p}^{r}_{2}+\tilde{b}^{r}).$

Finally the mean-field behavior is recovered for $N_{c}\gg 1$ as the result of
the following Theorem.

######  Theorem 1.

Consider the discrete mean-field model (9) for $N_{c}$ individuals. Then
taking the limit for $N_{c}\to\infty$, $\epsilon\to 0$, the mean-field
equations for the densities $f(x,t)$ and $g(x,t)$,

(11)
$\displaystyle{\partial_{\tau}f}=2\tilde{p}_{1}fg-\tilde{d}_{1}f+\tilde{m}_{1}(f\Delta
g+(1-g)\Delta f),$
$\displaystyle{\partial_{\tau}g}=rg\left(1-\frac{g}{q}\right)-\alpha
fg+\tilde{m}_{2}(g\Delta f+(1-f)\Delta g),$

where each entry of the vectors $f({\bf x},t)$ and $g({\bf x},t)$ represents
the predators and preys densities in each cell $C_{\ell}$,
$\ell=1,\ldots,M_{c}$ and

$\Delta h(x,t)=\lim_{\epsilon\to 0}{\Delta_{\epsilon}}h$

for any function h.

We refer to [26, 27, 37], for a detailed proof of the mean-field limit.

######  Remark 2.

Note that under condition $(1-q_{1}-q_{2})d^{r}_{2}<(2q_{1}b^{r})$ model (11)
corresponds to the well known Lotka-Volterra equations with logistic growth
term and diffusion, [28].

##  3\. Efficient ensemble stochastic algorithms

The dynamics described by model (11) can be properly simulated with stochastic
algorithms. Exact stochastic algorithms predict which is the next firing event
and at which time it will fire, as for example in [15]. However, a direct
implementation is often prohibitively expensive since the final goal is the
simulation of a large stochastic process. On the other hand, classic
approximated algorithms, such as the Monte Carlo algorithm and the
$\tau$-leaping method, can speed up the simulations but their efficiency and
accuracy are strongly related to the choice of the time step that can vary in
time or be inversely proportional to the sample size. The main idea of the
procedure is to fix a constant time step, and to allow multiple events to
happen at the same time. Hence, we describe the associated stochastic process
for the spatial heterogeneous predator-prey dynamics and show its consistency
with the previous description. Second we describe the ensemble stochastic
algorithm in the spatial homogeneous and heterogeneous setting.

###  3.1. Predator-prey model with ensemble interactions

At each time step $t\geq 0$ consider two fractions $q_{1},q_{2}\in[0,1]$ such
that $q_{1}+q_{2}\leq 1$. Assume that the following events can occur.

  * •

Competition & Birth. Sample $\lfloor q_{1}N_{c}\rfloor$ components in each
cell $C_{\ell}$, $\ell=1,\ldots,M_{c}$ assuming that each one interacts with
another component, randomly chosen among the $N_{c}$ positions in cell
$C_{\ell}$ without repetition, according to the rules described in (1).

  * •

Migration. Sample other $\lfloor q_{2}N_{c}\rfloor$ components in any cell
$C_{\ell}$, different from the ones previously selected, and let each one to
interact with another state, randomly chosen in cell $C_{\ell\pm 1}$ without
repetition, according to the migration rules defined in (2).

  * •

Death. The remaining $N_{c}-\lfloor q_{1}N_{c}\rfloor-\lfloor
q_{2}N_{c}\rfloor$ components in any cell $C_{\ell}$ change according to death
rules defined in equation (3).

Hence, in any infinitesimal time interval $[t,t+d\tau_{N_{c}}]$, where
$d\tau_{N_{c}}:=dtN_{c}$, a transition can occur from the state with ${\bf
x}=({\bf A},{\bf B},{\bf E})$ individuals to the states with ${\bf x}+k{\bf
v}_{\ell_{j}}$ individuals for any $k=1,\ldots,N_{c}$, where ${\bf
v}_{\ell_{j}}$ is the $\ell_{j}$-th row of the stoichiometry matrix
$\tilde{V}$ defined in (4) for any $\ell=1,\ldots,M_{c}$, $j=0,\ldots,M$. For
$k=1,\ldots N$, the transition rates write as follows

(12) $\pi({\bf x}+k{\bf v}_{\ell_{j}}|{\bf x})=\prod_{i=1}^{k}\pi^{i}_{{\bf
v}_{\ell_{j}}}({\bf x})$

where

(13) $\begin{split}&\pi_{{\bf v}_{\ell_{0}}}^{i}({\bf
x})=2b^{r}q_{1}\frac{\tilde{B}_{\ell}^{i}}{\tilde{N_{c}}^{i}}\frac{\tilde{E}_{\ell}^{i}}{\tilde{N}_{c}^{i}},~{}~{}\quad\qquad\pi_{{\bf
v}_{\ell_{1}}}^{i}({\bf
x})=2p_{1}^{r}q_{1}\frac{\tilde{A}_{\ell}^{i}}{\tilde{N}_{c}^{i}}\frac{\tilde{B}_{\ell}^{i}}{\tilde{N}_{c}^{i}},\\\
&\pi_{{\bf v}_{\ell_{2}}}^{i}({\bf
x})=2p_{2}^{r}q_{1}\frac{\tilde{A}_{\ell}^{i}}{\tilde{N}_{c}^{i}}\frac{\tilde{B}_{\ell}^{i}}{\tilde{N}_{c}^{i}},~{}~{}\quad\qquad\pi_{{\bf
v}_{\ell_{3}}}^{i}({\bf
x})=d_{1}^{r}(1-q_{1}-q_{2})\frac{\tilde{A}_{\ell}^{i}}{\tilde{N}_{c}^{i}},\\\
&\pi_{{\bf v}_{\ell_{4}}}^{i}({\bf
x})=d_{2}^{r}(1-q_{1}-q_{2})\frac{\tilde{B}_{\ell}^{i}}{\tilde{N}_{c}^{i}},\quad\pi_{{\bf
v}_{\ell_{5}}}^{i}({\bf
x})=m_{1}^{r}q_{2}\frac{\tilde{A}_{\ell}^{i}}{\tilde{N}_{c}^{i}}\frac{\tilde{E}_{\ell-1}^{i}}{\tilde{N}_{c}^{i}},\\\
&\pi_{{\bf v}_{\ell_{6}}}^{i}({\bf
x})=m_{1}^{r}q_{2}\frac{\tilde{A}_{\ell-1}^{i}}{\tilde{N}_{c}^{i}}\frac{\tilde{E}_{\ell}^{i}}{\tilde{N}_{c}^{i}},\qquad\quad\pi_{{\bf
v}_{\ell_{7}}}^{i}({\bf
x})=m_{2}^{r}q_{2}\frac{\tilde{B}_{\ell}^{i}}{\tilde{N}_{c}^{i}}\frac{\tilde{E}_{\ell-1}^{i}}{\tilde{N}_{c}^{i}},\\\
&\pi_{{\bf v}_{\ell_{8}}}^{i}({\bf
x})=m_{2}^{r}q_{2}\frac{\tilde{B}_{\ell-1}^{i}}{\tilde{N}_{c}^{i}}\frac{\tilde{E}_{\ell}^{i}}{\tilde{N}_{c}^{i}},\qquad\quad\pi_{{\bf
v}_{\ell_{9}}}^{i}({\bf
x})=m_{1}^{r}q_{2}\frac{\tilde{A}_{\ell}^{i}}{\tilde{N}_{c}^{i}}\frac{\tilde{E}_{\ell+1}^{i}}{\tilde{N}_{c}^{i}},\\\
&\pi_{{\bf v}_{\ell_{10}}}^{i}({\bf
x})=m_{1}^{r}q_{2}\frac{\tilde{A}_{\ell+1}^{i}}{\tilde{N}_{c}^{i}}\frac{\tilde{E}_{\ell}^{i}}{\tilde{N}_{c}^{i}},\qquad\quad\pi_{{\bf
v}_{\ell_{11}}}^{i}({\bf
x})=m_{2}^{r}q_{2}\frac{\tilde{B}_{\ell}^{i}}{\tilde{N}_{c}^{i}}\frac{\tilde{E}_{\ell+1}^{i}}{\tilde{N}_{c}^{i}},\\\
&\pi_{{\bf v}_{\ell_{12}}}^{i}({\bf
x})=m_{2}^{r}q_{2}\frac{\tilde{B}_{\ell+1}^{i}}{\tilde{N}_{c}^{i}}\frac{\tilde{E}_{\ell}^{i}}{\tilde{N}_{c}^{i}}.\end{split}$

Here we define the operators $\pi_{{\bf v}_{\ell_{j}}}^{i}({\bf x})=\pi({\bf
x}+i{\bf v}_{\ell_{j}}|{\bf x}+(i-1){\bf v}_{\ell_{j}})$ for any
$i=1,\ldots,N_{c}$, we use the following notation $\tilde{Y}^{i}=Y-i+1$ for
any state $Y$ in the set $\\{A_{\ell},B_{\ell},E_{\ell},A_{\ell\pm
1},B_{\ell\pm 1},E_{\ell\pm 1},N_{c}\\}$ and since migration outside the
boundaries is not allowed, we introduce
$A_{0},E_{0},B_{0},A_{M_{c}+1},E_{M_{c}+1},B_{M_{c}+1}$ to be equal to zero.
The master equation associated to this process is derived in the following
Proposition.

######  Proposition 2 (Consistency).

The master equation can be written as

(14) $\frac{dP({\bf x},s)}{ds}=\sum_{\ell_{j}\in\mathcal{J}}\Biggl{[}\pi({\bf
x}|{\bf x}-{\bf v}_{\ell_{j}})P({\bf x}-{\bf v}_{\ell_{j}},s)-\pi({\bf x}+{\bf
v}_{\ell_{j}}|{\bf x})P({\bf x},s)\Biggr{]},$

where $s=tN_{c}$, the set $\mathcal{J}$ is defined as in (7) and
$|\mathcal{J}|=M_{c}(M+1)$ is the number of events that can occur.

###### Proof.

The probability to be in the state with ${\bf x}$ individuals at time
$t+\tau_{N_{c}}$ is given by two contributions:

  * •

the probability of staying in the state ${\bf x}$ at time $t+\tau_{N_{c}}$
that is

(15) $1-\sum_{\ell_{j}\in\mathcal{J}}\sum_{k=1}^{N_{c}}\pi({\bf x}+k{\bf
v}_{\ell_{j}}|{\bf x})d\tau_{N_{c}}^{k};$

  * •

the probability of moving from the state with ${\bf x}$ individuals to the
state with ${\bf x}+k{\bf v}_{\ell_{j}}$ individuals in the time interval
$[t,t+\tau_{N_{c}}]$, for any $k=1,\ldots,N_{c}$ and $\ell_{j}\in\mathcal{J}$
that is

(16) $\sum_{\ell_{j}\in\mathcal{J}}\sum_{k=1}^{N_{c}}\pi({\bf x}|{\bf x}-k{\bf
v}_{\ell_{j}})d\tau_{N_{c}}^{k}.$

Indeed, recall that for any $i=1,\ldots,N_{c}$ a transition from the state
with ${\bf x}+i{\bf v}_{\ell_{j}}$ individuals to the state with ${\bf
x}+(i-1){\bf v}_{\ell_{j}}$ individuals in the infinitesimal time interval
$[t,t+\tau_{N_{c}}]$ occurs with probability $\pi_{\ell_{j}}^{i}({\bf
x})d\tau_{N_{c}}$ and hence the probability of moving from the state with
${\bf x}$ to the state with ${\bf x}+k{\bf v}_{\ell_{j}}$ individuals, for any
$k=1,\ldots,N_{c}$, is

$\prod_{i=1}^{k}\left\\{\pi_{{\bf v}_{\ell_{j}}}^{i}({\bf
x})d\tau_{N_{c}}\right\\}=\pi({\bf x}+k{\bf v}_{\ell_{j}}|{\bf
x})d\tau_{N_{c}}^{k},$

that is exactly the expression that appears in equation (15). Similarly we can
derive equation (16). Hence,

(17) $\begin{split}P({\bf
x},t+d\tau_{N_{c}})=&\sum_{\ell_{j}\in\mathcal{J}}\sum_{k=1}^{N_{c}}\Big{[}P({\bf
x}-k{\bf v}_{\ell_{j}},t)\pi({\bf x}|{\bf x}-k{\bf
v}_{\ell_{j}})d\tau_{N_{c}}^{k}\Big{]}\\\ &+P({\bf
x},t)\Biggl{(}1-\sum_{\ell_{j}\in\mathcal{J}}\sum_{k=1}^{N_{c}}\Big{[}\pi({\bf
x}+k{\bf v}_{j}|{\bf x})d\tau_{N_{c}}^{k}\Big{]}\Biggr{)}.\end{split}$

Rewrite explicitly the term for $k=1$ in equation (17) to have

(18) $\begin{split}&P({\bf x},t+d\tau_{N_{c}})-P({\bf
x},t)=d\tau_{N_{c}}\sum_{\ell_{j}\in\mathcal{J}}\Big{[}P({\bf x}-{\bf
v}_{j},t)\pi({\bf x}|{\bf x}-{\bf v}_{j})-P({\bf x},t)\pi({\bf x}+{\bf
v}_{j}|{\bf x})\Big{]}\\\
&+d\tau_{N_{c}}^{k}\sum_{\ell_{j}\in\mathcal{J}}\sum_{k=2}^{N_{c}}\Big{[}P({\bf
x}-k{\bf v}_{j},t)\pi({\bf x}|{\bf x}-k{\bf v}_{j})-P({\bf x},t)\pi({\bf
x}+k{\bf v}_{j}|{\bf x})\Big{]}.\end{split}$

Dividing both sides of equation (18) by $d\tau_{N_{c}}$ and letting
$d\tau_{N_{c}}\to 0$ we obtain the consistency with the master equation (6),
where the time variable is scaled by a factor $N_{c}$, i.e. $s=tN_{c}$. ∎

######  Remark 3.

We observe that thanks to the consistency result of Proposition 2, the
corresponding mean-field approximation is equivalent to (11). According to
(11) the time scale of the mean-field dynamics is such that

$\tau=s/{N_{c}}=N_{c}t/N_{c}=t,$

then the time scales of the ensemble agent-based dynamics and of the mean-
field dynamics are equivalent.

###  3.2. Efficient Monte-Carlo methods

The idea of the efficient Monte Carlo algorithm is to allow multiple events to
occur at the same time step that is supposed to be fixed and constant for any
choice of $N_{c}$. Therefore, as we will see in the numerical experiments, its
accuracy is comparable with the one of classic approximated algorithms but its
efficiency is improved. In the following we will focus first on the spatial
homogeneous case, where a single cell is accounted ($M_{c}=1$) and $N_{c}=N$,
second we will consider the full dynamics with spatial interaction. We refer
to A for further details on the homogeneous case.

#### Spatially homogeneous case

Divide a priori the time interval considering a constant time step $\tau$ that
is independent of the sample size. Introduce a parameter $\mu\in[0,1]$ and at
each time $t$ assume that simultaneously $\lfloor\mu N\rfloor$ individuals
interact two by two. If the two individuals are

  * •

a prey and an empty space: with probability $b^{r}\tau$ a new prey born and
occupies the empty space;

  * •

a predator and a prey: with probability $p_{1}^{r}\tau$ then a new predator
born and with probability $p_{2}^{r}\tau$ the prey dies and a new empty space
is added to the system.

Assume that the remaining $N-\lfloor\mu N\rfloor$ individuals are subjected to
death events that happen with probability $d_{1}^{r}\tau$ in the predators
population and with probability $d_{2}^{r}\tau$ in the preys one. Algorithm 1
defines the details of the efficient Monte Carlo algorithm.

######  Algorithm 1 (Efficient ensemble stochastic algorithm - homogeneous).

  1. 1.

Define the sample, the initial time $t=0$, the final time $T$, the time step
$\tau$ and a parameter $\mu\in[0,1]$ .

  2. 2.

while $t<T$

    1. (a)

Birth and competition events happen between $\lfloor\mu N\rfloor$ individuals
selected two by two

$s_{B}s_{E}\xrightarrow{b}s_{B}s_{B},\quad
s_{A}s_{B}\xrightarrow{p_{1}}s_{A}s_{A},\quad
s_{A}s_{B}\xrightarrow{p_{2}}s_{A}s_{E},$

according to probability $b=b^{r}\tau$, $p_{1}^{r}\tau$ and $p_{2}^{r}\tau$,
respectively.

    2. (b)

The remaining $N-\lfloor\mu N\rfloor$ individuals are subjected to death
events

$s_{A}\xrightarrow{d_{1}}s_{E},\quad s_{B}\xrightarrow{d_{2}}s_{E},$

according to probability $d_{1}=d_{1}^{r}\tau$ and $d_{2}=d_{2}^{r}\tau$,
respectively.

    3. (c)

Update the sample.

    4. (d)

Set $t\leftarrow t+\tau$.

repeat

#### Spatially heterogeneous case

We extend Algorithm 1 to the spatially heterogeneous case. Hence, the sample
is divided in cells and populations are subjected also to migration events.
Migrations can occur either between cell $C_{\ell}$ and cell $C_{\ell+1}$ or
between cell $C_{\ell}$ and cell $C_{\ell-1}$ for any $\ell=1,\ldots,M_{c}$.
The spatially heterogeneous Algorithm 1 reads as follows

######  Algorithm 2 (Efficient ensemble stochastic algorithm).

  1. 1.

Define the sample, the initial time $t=0$, the final time $T$, the time step
$\tau$ and two parameters $q_{1},q_{2}\in[0,1]$ such that $q_{1}+q_{2}\leq 1$.

  2. 2.

while $t<T$

    1. (a)

In each cell $C_{\ell}$, $\lfloor q_{1}N_{c}\rfloor$ individuals interact two
by two and are subjected to birth and competition events (1) that occur with
probabilities $b=b^{r}\tau$, $p_{1}=p_{1}^{r}\tau$ and $p_{2}=p_{2}^{r}\tau$;

    2. (b)

In each cell $C_{\ell}$, $\lfloor q_{2}N_{c}/2\rfloor$ individuals interact
with $\lfloor q_{2}N_{c}/2\rfloor$ individuals sampled in cell $C_{\ell+1}$
and other $\lfloor q_{2}N_{c}/2\rfloor$ interact with $\lfloor
q_{2}N_{c}/2\rfloor$ individuals sampled in cell $C_{\ell-1}$ according to
(2), and migrate with probability $m_{1}=m_{1}^{r}\tau$ for predators and
$m_{2}=m_{2}^{r}\tau$ for preys;

    3. (c)

In each cell $C_{\ell}$, the remaining $N_{c}-\lfloor
q_{1}N_{c}\rfloor-\lfloor q_{2}N_{c}\rfloor$ individuals are subjected to
death events (3) according to probabilities $d_{1}=d_{1}^{r}\tau$ and
$d_{2}=d_{2}^{r}\tau$;

    4. (d)

Update the sample;

    5. (e)

Set $t\leftarrow t+\tau$.

##  4\. Numerical experiments

###  4.1. Test 1: Validation

In this section we present a comparison between the numerical solutions of the
mean-field equations and the stochastic simulations. In both the homogeneous
and heterogeneous case, stochastic simulations have been performed with the
efficient version of the Monte Carlo algorithm presented in Section 3. In the
homogeneous case, the numerical solutions of the mean-field equations are
computed using the Matlab function $\mathtt{ode45}$, [24], while in the
heterogeneous case with a combination of finite difference methods and
numerical methods for ODEs assuming periodic boundary conditions, [22]. The
parameters choice for all the tests in Section 4.1 is specified in Table 1.
The sample size $N$ and the total number of individuals $N_{c}$ in any cell
change in any test and will be defined later.

Table 1. Model parameters for the different scenarios. | $b^{r}$ | $d_{1}^{r}$ | $d_{2}^{r}$ | $p_{1}^{r}$ | $p_{2}^{r}$ | $m_{1}^{r}$ | $m_{2}^{r}$ | $\mu$ | $q_{1}$ | $q_{2}$
---|---|---|---|---|---|---|---|---|---|---
Homogeneous | 0.1 | 0.1 | 0 | 0.25 | 0.05 | - | - | 0.5 | - | -
Heterogeneous | 0.1 | 0.1 | 0 | 0.25 | 0.05 | 0.5 | 0.5 | - | 0.3 | 0.3

Figure 1 shows a comparison between stochastic and mean-field solutions in the
homogeneous case assuming the sample size to be $N=1000$ and the initial
predators and preys density to be $A_{0}=N/4$ and $B_{0}=N/2$ respectively.
Note that the solution of the mean-field equations (11), without spatial
dependence, converges to the stable equilibrium

(19)
$f^{*}=\frac{2\tilde{b}^{r}\tilde{p}^{r}_{1}-\tilde{b}^{r}\tilde{d}^{r}_{1}-\tilde{p}^{r}_{1}\tilde{d}^{r}_{2}}{2\tilde{p}^{r}_{1}(\tilde{p}^{r}_{1}+\tilde{p}^{r}_{2}+\tilde{b}^{r})},\qquad
g^{*}=\frac{\tilde{d}^{r}_{1}}{2\tilde{p}^{r}_{1}},$

that for the parameters choice of Table 1 takes the values $f^{*}=g^{*}=0.2$.

Figure 2 shows that the error of the efficient Monte Carlo algorithm is
proportional to $1/\sqrt{N}$, as the one of classic algorithms in both the
predators (left) and preys (right) cases. The errors $e_{f}^{N}$ and
$e_{g}^{N}$ are defined as

(20) $\begin{split}e_{f}^{N}=\|f(t)-f^{N}(t)\|_{\infty},\qquad
e_{g}^{N}=\|g(t)-g^{N}(t)\|_{\infty},\end{split}$

where $f(t)$, $g(t)$ denote the mean-field solutions and $f^{N}(t)$,
$g^{N}(t)$ denote the stochastic solutions obtained with the efficient Monte
Carlo algorithm at time $t$ for a certain value of $N$. Recall that $f(t)$,
$f^{N}(t)$ refer to the predators population while $g(t)$, $g^{N}(t)$ to the
preys one.

Figure 1. Homogeneous predator-prey model: simulation of the processes
described in (26)-(27) with the efficient Monte Carlo algorithm and solutions
of the spatial homogeneous mean-field equations (11) for $N=1000$. Markers
have been added just to indicate different lines.

Figure 2. Homogeneous predator-prey model: efficient Monte Carlo algorithm’s
error computed as in equations (20) for $N=[10,\ldots,10^{4}]$. The error is
proportional to $1/\sqrt{N}$ as the one of the classic Monte Carlo algorithm.
Markers correspond to the values $e_{f}^{N}$, $e_{g}^{N}$ for a fixed $N$.

Let us now focus on the heterogeneous model and consider first the one
dimensional case. Assume that the dynamics evolves in an interval area of land
$[0,L_{x}]$, $L_{x}>0$, divided in $M_{c}$ cells. Figure 3 shows three
snapshots taken at time $t=5$, $t=50$ and $t=100$ in which the stochastic and
mean-field solutions are compared. At time $t=5$, $A_{0}=N_{c}/4$ predators
and $B_{0}=N_{c}/2$ preys are concentrated in the same central cells. At time
$t=50$, preys migrate in regions where the predators concentration is lower
and predators decrease their size and start to migrate to reach the regions
occupied by preys. At time $t=100$, preys continue their migration increasing
their size in the regions where the predators concentration is lower.

Figure 3. Heterogeneous one dimensional predator-prey model: simulation of the
processes described in (1)-(2)-(3) with the efficient Monte Carlo algorithm
and solutions of the mean-field equations (11) for $N_{c}=1000$. This figure
shows three snapshots taken at time $t=5$ (left), $t=50$ (centre), $t=100$
(right). Markers have been added just to indicate different lines.

Figure 4 shows the asymptotic behavior of the two populations at the mean-
field and stochastic level. Predators and preys migrate in the whole available
space reaching in each cell the value given by the equilibrium (19).

Figure 4. Heterogeneous one dimensional predator-prey model: asymptotic
behavior (at time $t=500$) of predators (on the left) and preys (on the right)
populations at the mean-field and stochastic level for $N_{c}=1000$. Markers
have been added just to indicate different lines.

Figure 5 shows that the error of the efficient Monte Carlo algorithm as a
function of $N_{c}$ is proportional to $1/\sqrt{N_{c}}$, as the one of classic
algorithms. The errors are computed as

(21)
$\begin{split}&e_{f}^{N_{c}}=\left\langle\max_{t}|f(x,t)-f^{N_{c}}(x,t)|\right\rangle_{x},\,\,\,e_{g}^{N_{c}}=\left\langle\max_{t}|g(x,t)-g^{N_{c}}(x,t)|\right\rangle_{x},\end{split}$

where $\left\langle\cdot\right\rangle$ denotes the expected value with respect
to $x$.

Figure 5. Heterogeneous one dimensional predator-prey model: efficient Monte
Carlo algorithm’s error computed as in equations (21) for
$N_{c}=[10,\ldots,10^{4}]$. The error is proportional to $1/\sqrt{N_{c}}$, as
the one of the classic Monte Carlo algorithm. Markers correspond to the values
$e_{f}^{N_{c}}$, $e_{g}^{N_{c}}$ for a fixed $N_{c}$.

Let us now focus on the two dimensional case. Assume that the dynamics evolves
in a square area of land $[0,L_{x}]\times[0,L_{y}]$, $L_{x},L_{y}>0$ divided
in $C_{\ell_{x}\ell_{y}}$ cells, for $\ell_{x}=1,\ldots,M_{c}^{x}$ and
$\ell_{y}=1,\ldots,M_{c}^{y}$. Populations can born, compete, die and migrate
in one of the four nearest cells. The error computed as in equation (21)
assuming $x=(x,y)$ is still proportional to $1/\sqrt{N_{c}}$, as shown in
Figure 6.

Figure 6. Heterogeneous two dimensional predator-prey model: efficient Monte
Carlo algorithm’s error computed as in equations (21) for
$N_{c}=[10,\ldots,10^{4}]$. The error is proportional to $1/\sqrt{N_{c}}$, as
the one of the classic Monte Carlo algorithm. Markers correspond to the values
$e_{f}^{N_{c}}$, $e_{g}^{N_{c}}$ for a fixed $N_{c}$.

Figure 7 and Figure 8 show three snapshots describing the time evolution of
preys and predators densities at the mean-field and stochastic level in the
heterogeneous two dimensional case for $N_{c}=1000$. At time $t=5$,
$B_{0}=N_{c}/2$ preys are concentrated in the central cells and surrounded by
$A_{0}=N_{c}/4$ predators. At time $t=100$ predators migrate in the central
cells while preys reduce their size and start to migrate in the regions where
the predators concentration is lower. At time $t=150$ preys are still
migrating and increasing their size. Predators on the contrary are reducing
their size and migrating to reach the regions in which preys are mainly
concentrated.

Figure 7. Heterogeneous two dimensional predator-prey model (preys
population): simulation of the processes described in (1)-(2)-(3) with the
efficient Monte Carlo algorithm and solutions of the mean-field equations (11)
in the two dimensional case for $N_{c}=1000$. This figure shows three
snapshots taken at time $t=5$ (top), $t=100$ (middle), $t=150$ (bottom). On
the left, mean-field solutions and on the right, stochastic simulations.

Figure 8. Heterogeneous two dimensional predator-prey model (predators
population): simulation of the processes described in (1)-(2)-(3) with the
efficient Monte Carlo algorithm and solutions of the mean-field equations (11)
in the two dimensional case for $N_{c}=1000$. This figure shows three
snapshots taken at time $t=5$ (top), $t=100$ (middle), $t=150$ (bottom). On
the left, mean-field solutions and on the right, stochastic simulations.

In Figure 9 we see for simplicity just the asymptotic behavior of the
predators population. One can show that both populations in long time migrate
in the whole available space reaching in each cell the value given by the
equilibrium (19).

Figure 9. Heterogeneous two dimensional predator-prey model: asymptotic
behaviour (at time $t=500$) of predators population at the mean-field (on the
right) and stochastic (on the left) level for $N_{c}=1000$.

######  Remark 4.

Note that the mean-field solutions of equations (11) present a damped behavior
in time while the stochastic solutions obtained with the efficient Monte Carlo
algorithm have a persistent behavior in time. One can prove that the
stochastic persistency is due to a resonant effect, [27]. We refer to section
4.4 for a dedicated test to analyzed this behavior.

###  4.2. Test 2: Computational cost

Let us consider the homogeneous case and assume to fix the parameters
$\mu=0.5$, $b^{r}=1$, $d_{1}^{r}=d_{2}^{r}=0.3$, and to let the competition
parameters $p^{r}_{1},p^{r}_{2}$ to vary between $0.1$ and $0.9$. Figure 10
shows that the computational cost of the efficient Monte Carlo algorithm is
lower than the one of the other algorithms.

Figure 10. Homogeneous predator-prey model: computational cost of efficient
and classic Monte Carlo algorithms, direct method and $\tau$-leaping method as
the competition parameters vary for $N$ fixed. The dynamics in (26)-(27) is
simulated for $N=500$, $\mu=0.5$, $b^{r}=1$, $d_{1}^{r}=d_{2}^{r}=0.3$,
$p_{1}^{r}=p_{2}^{r}=0.1,\ldots,0.9$. Markers represent the computational
costs relative to the parameters choice indicated.

Figure 11 shows that the efficient Monte Carlo algorithm, the direct method
and the $\tau$-leaping algorithm have a computational complexity of order $N$
in time while the one of the classic Monte Carlo algorithm is of order
$N^{2}$.

Figure 11. Homogeneous predator-prey model: computational complexity of
efficient Monte Carlo algorithm, direct method and $\tau$-leaping method (on
the left), and classic Monte Carlo algorithms (on the right) as $N$ varies.The
dynamics in (26)-(27) is simulated for $N=500,\ldots,4000$, $\mu=0.5$,
$b^{r}=1$, $d_{1}^{r}=d_{2}^{r}=0.3$, $p_{1}^{r}=p_{2}^{r}=0.5$. Markers
represent the computational costs relative to the parameters choice indicated.

Figure 12 shows the comparison between the computational costs of the
stochastic algorithms in the one dimensional heterogeneous case. The dynamics
in (1)-(2)-(3) is simulated for $N_{c}=100$, $q_{1}=0.3$, $q_{2}=0.3$ and
$b^{r}=0.1$, $d_{1}^{r}=0.1$, $d_{2}^{r}=0$, $p_{1}^{r}=0.25$,
$p_{2}^{r}=0.05$, letting the migration parameters $m_{1}^{r}$, $m_{2}^{r}$ to
vary between $0.1$ and $0.9$. The computational cost of the efficient Monte
Carlo algorithm is lower than the one of classic algorithms.

Figure 12. Heterogeneous one dimensional predator-prey model: computational
cost of efficient and classic Monte Carlo algorithms, direct method and
$\tau$-leaping method as the migration parameters vary, for $N_{c}$ fixed. The
dynamics in (1)-(2)-(3) is simulated for $N_{c}=100$, $q_{1}=0.3$,
$q_{2}=0.3$, $b^{r}=0.1$, $d_{1}^{r}=0.1$, $d_{2}^{r}=0$, $p_{1}^{r}=0.25$,
$p_{2}^{r}=0.05$, $m_{1}^{r}=m_{2}^{r}=0.1,\ldots,0.9$. Markers represent the
computational costs relative to the parameters choice indicated.

Figure 13 shows the computational cost of the stochastic algorithms for
different values of the migration rates as $N_{c}$ varies. On the left the
computational costs for fixed migration rates $m_{1}^{r}=m_{2}^{r}=0.1$ and on
the right for $m_{1}^{r}=m_{2}^{r}=0.9$. Note that the computational cost of
the efficient Monte Carlo algorithm is always lower than the one of the
classic algorithms.

Figure 13. Heterogeneous one dimensional predator-prey model: computational
cost of efficient and classic Monte Carlo algorithms, direct method and
$\tau$-leaping method as $N_{c}$ varies. The dynamics in (1)-(2)-(3) is
simulated for $N_{c}=100,\ldots,1000$, $q_{1}=0.3$, $q_{2}=0.3$, $b^{r}=0.1$,
$d_{1}^{r}=0.1$, $d_{2}^{r}=0$, $p_{1}^{r}=0.25$, $p_{2}^{r}=0.05$,
$m_{1}^{r}=m_{2}^{r}=0.1$(on the left), $m_{1}^{r}=m_{2}^{r}=0.9$ (on the
right). Markers represent the computational costs relative to the parameters
choice indicated.

In Figure 14 a comparison between the computational costs of the stochastic
algorithms in the two dimensional heterogeneous case. The dynamics is
simulated for $N_{c}=50$ fixed, $q_{1}=0.3$, $q_{2}=0.3$, $b^{r}=0.1$,
$d_{1}^{r}=0.1$ $d_{2}^{r}=0$, $p_{1}^{r}=0.25$, $p_{2}^{r}=0.05$, letting the
migration parameters to vary between $0.1$ and $0.9$. Also in this case the
efficient Monte Carlo algorithm has a computational cost lower than the one of
classic algorithms.

Figure 14. Heterogeneous two dimensional predator-prey model: computational
cost of efficient and classic Monte Carlo algorithms, direct method and
$\tau$-leaping method as the migration parameters vary, for $N_{c}$ fixed. The
dynamics in (1)-(2)-(3) is simulated for $N_{c}=50$, $q_{1}=0.3$, $q_{2}=0.3$,
$b^{r}=0.1$, $d_{1}^{r}=0.1$, $d_{2}^{r}=0$, $p_{1}^{r}=0.25$,
$p_{2}^{r}=0.05$, $m_{1}^{r}=m_{2}^{r}=0.1,\ldots,0.9$. Markers represent the
computational costs relative to the parameters choice indicated.

###  4.3. Test 3: Accuracy & performances

In the following we will compare the accuracy of the direct method with the
one of the efficient Monte Carlo algorithm and of the $\tau$-leaping method.
We first focus on the homogeneous case and we define the errors as

(22) $E_{f}^{N}=\|f^{N}_{DM}(t)-f^{N}\|_{\infty},\quad
E_{g}^{N}=\|g^{N}_{DM}(t)-g^{N}\|_{\infty},$

where $f^{N}_{DM}(t)$, $g^{N}_{DM}(t)$ denote the simulations obtained with
the direct method and $f^{N}(t)$, $g^{N}(t)$ denote the simulations obtained
either with the efficient Monte Carlo or with the $\tau$-leaping algorithms
for the predators and preys populations, respectively. Figure 15 shows the
errors $E_{f}^{N}$, $E_{g}^{N}$ as $N$ varies and for fixed birth, competition
and death parameters (see Table 1). Note that the error related to the
simulations obtained with the efficient Monte Carlo algorithm is lower than
the one related to the simulated $\tau$-leaping solutions for almost every
$N$. Note also that both errors decrease as the sample size increases reaching
a constant value when $N$ is large enough.

Figure 15. Homogeneous predator-prey model: accuracy of the efficient Monte
Carlo and of the $\tau$-leaping algorithms with respect to the direct method
as $N$ varies. The dynamics in (26)-(27) is simulated for
$N=10^{2},\ldots,5\times 10^{6}$, $\mu=0.5$, $b^{r}=0.1$, $d_{1}^{r}=0.1$,
$d_{2}^{r}=0$, $p_{1}^{r}=0.25$, $p_{2}^{r}=0.05$. Markers correspond to the
values $E_{f}^{N},E_{g}^{N}$ computed as in equations (22) for a fixed $N$.

The same results can be obtained in the heterogeneous case with the parameters
choice reported in Table 1, as show in Figure 16. Here, for simplicity we
focus on the one dimensional case and we define the errors as

(23)
$\begin{split}&E_{f}^{N_{c}}=\left\langle\max_{t}|f^{N_{c}}(x,t)-f^{N_{c}}(x,t)|\right\rangle_{x},\quad
E_{g}^{N_{c}}=\left\langle\max_{t}|g^{N_{c}}(x,t)-g^{N_{c}}(x,t)|\right\rangle_{x},\end{split}$

where $\left\langle\cdot\right\rangle$ denotes the expected value with respect
to $x$.

Figure 16. Heterogeneous predator-prey model: accuracy of the efficient Monte
Carlo and of the $\tau$-leaping algorithms with respect to the direct method
as $N_{c}$ varies. The dynamics in (1)-(2)-(3) is simulated for
$N_{c}=10^{2},\ldots,10^{5}$, $q_{1}=q_{2}=0.3$, $b^{r}=0.1$, $d_{1}^{r}=0.1$,
$d_{2}^{r}=0$, $p_{1}^{r}=0.25$, $p_{2}^{r}=0.05$, $m_{1}=0.5$, $m_{2}=0.5$.
Markers correspond to the values $E_{f}(N_{c}),E_{g}(N_{c})$ computed as in
equations (23) for a fixed $N_{c}$.

Figure 17 shows the trade off between the error $E_{f}$ and the computational
costs of the efficient Monte Carlo algorithm and of the $\tau$-leaping method
in both the homogeneous and heterogeneous case for the predators dynamics.
Note that both errors decrease as the computational time increases, that is as
the sample size increases. Note also that in the homogeneous case the
computational cost is lower but the error is higher than in the heterogeneous
case. Here, for simplicity we focus on the predators dynamics but the same
results can be obtained if we focus on the preys dynamics.

Figure 17. Homogeneous and heterogeneous one dimensional predator-prey model:
trade-off of error $E_{f}^{N_{c}}$ computed as in equation (22)-(23) vs
computational cost of the efficient Monte Carlo and of the $\tau$-leaping
algorithms as the computational cost varies. On the left the trade-off in the
homogeneous case, on the right in the heterogeneous case. Markers correspond
to the values $E_{f}^{N_{c}}$ in relation to the computational costs.

###  4.4. Test 4: Stochastic persistency

We have observed in Figure 1 that both mean-field solutions and stochastic
simulations present an oscillatory behavior. However, the nature of the
oscillations seems to be different. Following the idea in [27], we can prove
that oscillations in predators and preys stochastic simulations should be of
order $1/\sqrt{N}$, as the error of the Monte Carlo method, but indeed are
amplified by a large factor due to noise resonant effect.

In order to analyze the nature of oscillations we introduce the following
linearized model corresponding to a perturbation of the mean-field dynamics
(9) in absence of spatial diffusion, see also Appendix A,

(24) $\frac{d}{dt}\begin{bmatrix}{f^{N}}\\\
{g^{N}}\end{bmatrix}=\Psi\begin{bmatrix}f^{N}\\\ g^{N}\end{bmatrix}+\Phi\xi,$

where $\Psi$ is the Jacobian matrix of the mean-field equations evaluated at
the equilibrium $(f^{*},g^{*})$, defined in (19), $\Phi$ corresponds to the
following matrix

$\Phi=\begin{bmatrix}0&\sqrt{2\tilde{p}_{1}^{r}f^{*}g^{*}}&0&-\sqrt{\tilde{d}_{1}^{r}f^{*}}\\\
\sqrt{2\tilde{b}^{r}g^{*}(1-f^{*}-g^{*})}&-\sqrt{2\tilde{p}_{1}^{r}f^{*}g^{*}}&-\sqrt{2\tilde{p}_{2}^{r}f^{*}g^{*}+\tilde{d}_{2}^{r}g^{*}}&0\end{bmatrix},$

associated to the white noise vector
$\xi=[\xi_{1},\xi_{2},\xi_{3},\xi_{4}]^{T}$, with $\xi_{i}$ i.i.d distributed
random numbers with zero mean and variance $1/\sqrt{N}$, [17],[36]. We focus
on the predator density, we compute the Fourier transform of the mean-field
density $f(t)$, i.e. $\tilde{f}(\omega)$, and we derive the power spectrum

(25)
$P(\omega)=\left\langle|\tilde{f}(\omega)|^{2}\right\rangle=\frac{\Theta+\Lambda\omega^{2}}{(\omega^{2}-\Omega_{0}^{2})^{2}+\Gamma^{2}\omega^{2}},$

where $\left\langle\cdot\right\rangle$ is the expected value w.r.t. $\xi$. The
parameters $\Theta,\Lambda,\Omega_{0},\Gamma$ are defined as follows

$\displaystyle\Theta$
$\displaystyle=\left\langle[\Psi_{12}(\Phi_{21}\xi_{1}+\Phi_{22}\xi_{2}+\Phi_{23}\xi_{3})-\Psi_{22}(\Phi_{12}\xi_{2}+\Phi_{14}\xi_{4})]^{2}\right\rangle$
$\displaystyle\Lambda$
$\displaystyle=\left\langle[(\Phi_{12}\xi_{2}+\Phi_{14}\xi_{4}]^{2}\right\rangle,\quad\Omega_{0}^{2}=\Psi_{12}|\Psi_{21}|,\quad\Gamma=|\Psi_{22}|.$

Note that $P(\omega)$ resembles the one of a damped harmonic oscillator, [19].
As shown in Figure 18 on the left, the solution $f(t)$ of the system of
equations (24) without the influence of noise behaves as an under-damped
harmonic oscillator. The damped term is $\Gamma<1$ and oscillations decrease
exponentially as $\tilde{A}\exp\left(-\frac{\Gamma}{2}t\right)$ where
$\tilde{A}$ is the maximum oscillations amplitude. This behavior recalls the
one of the mean-field solutions. If we introduce the noise term, we see in
Figure 18 on the right that the oscillations amplitude is influenced both by
damped and resonant effects and the oscillations become sustained as the one
of the stochastic simulations shown in Section 4. Hence the resonant effect is
a consequence of the white noise that is not an external factor but it is due
to the stochasticity of birth, death and competitions events.

Figure 18. Predators density without noise (on the left) and with noise (on
the right) as solution of equations (24).

In Figure 19 the power spectrum computed as in equation (25) and the one
obtained averaging the results of 500 simulations of the processes described
in (26)-(27) computed with the efficient Monte Carlo algorithm and with the
direct method for different values of $N$ and for the parameters choice
reported in the first line of Table 1. Note that the two simulated power
spectra agree with the one computed as in equation (25), especially when $N$
is large. Note also that the amplitude of the oscillations decreases as $N$
increases.

Figure 19. Homogeneous case: power spectrum computed as in equation (25) in
red, averaging 500 simulations of the processes described in (26)-(27)
computed with the efficient Monte Carlo algorithm in blue and with the direct
method in black. On the left $N=10^{3}$ and on the right $N=10^{5}$. Markers
have been added just to indicate different lines.

The same results can be obtained in the heterogeneous case in each cell
$C_{\ell}$, $\ell=1,\ldots,M_{c}$, as the time $t$ varies. In Figure 20, we
see an example of persistent oscillatory behavior in stochastic simulations
compared with the damped one of the mean-field solutions in cell $C_{\ell}$,
$\ell=25,50$ for $N_{c}=100$, $M_{c}=100$.

Figure 20. Heterogeneous case: time evolution of the predators density in cell
$C_{\ell}$, $\ell=25,50$. Parameters: $N_{c}=100$, $b^{r}=0.1$,
$d_{1}^{r}=0.1$, $d_{2}^{r}=0$, $p_{1}^{r}=0.25$, $p_{2}^{r}=0.05$,
$m_{1}^{r}=m_{2}^{r}=0.5$, cell $C_{\ell}$, $\ell=25,50$. Markers have been
added just to indicate different lines.

In Figure 21 on the left the power spectrum computed averaging the results of
500 simulations of the processes described in (1)-(2)-(3) in cell $C_{\ell}$,
$\ell=25,50$ with the efficient Monte Carlo algorithm for the parameters
choice reported in the second line of Table 1 for $N_{c}=100$ and $T=2000$. On
the right, the average between the power spectra computed in cell $C_{\ell}$,
$\ell=1,\ldots,M_{c}$ obtained with the efficient Monte Carlo algorithm and
with the direct method. Note that also in the heterogeneous case the power
spectrum recalls the one of a damped harmonic oscillator. Note also that the
two average power spectra agree and so that the oscillations that characterize
the simulations obtained with the two algorithms are of the same nature.

Figure 21. Heterogeneous case: power spectrum. On the left, power spectrum
computed averaging 500 simulations of the processes described in (1)-(2)-(3)
in cell $C_{\ell}$, $\ell=25,50$ obtained with the efficient Monte Carlo
algorithm. On the right the average between the power spectra computed in cell
$C_{\ell}$, $\ell=1,\ldots,M_{c}$ obtained with the efficient Monte Carlo
algorithm (in blue) and with the direct method (in magenta). Markers have been
added just to indicate different lines.

##  5\. Conclusion

Starting from the microscopic level it is possible to derive the master
equation that describes the time evolution of the probability of being in a
certain state and the mean-field equations that provide an average description
of the behavior of the system which is valid only when the sample size $N$ is
arbitrary large. Stochastic models instead describe the individuals behavior
of interacting species in a more realistic context. However, classic
stochastic algorithms have high computational costs, especially when the
sample size increases. In this paper, we have presented a valid efficient
stochastic algorithm that, compared to other algorithms, has lower
computational costs,and equivalent, or even better performances. First we have
proved its consistency providing a reformulation of the model and showing that
the mean-field equations are still valid. Then with different numerical
experiments we have proved that the efficient Monte Carlo algorithm produces
solutions that oscillate around the mean-field ones and that the error between
stochastic simulations and mean-field solutions as a function of $N$ is still
proportional to $1/\sqrt{N}$. We have also shown which is the main advantage
of this new method proving that its computational cost is lower than the one
of classical algorithms for a huge set of parameters but its accuracy is
preserved. In the end, we have compared the oscillatory behavior of the
stochastic and mean-field solutions underling the differences in their
oscillatory behavior. Due to the flexible formulation of the ensemble
stochastic algorithm it can be used not only to simulate the microscopic
dynamics of interacting biological entities, but further generalized to other
agent-based models, with different interaction dynamics and different
macroscopic limiting models.

##  Appendix A Spatial homogeneous predator-prey model

We report here the agent-based model associated to spatial homogeneous
dynamics.

#### Agent-based dynamics

At each instant of time we select one individual and with probability
$\mu\in[0,1]$ we assume that an interaction occurs with another individual,
randomly chosen among the remaining $N-1$, according to the following rules

(26) $\begin{split}&s_{B}s_{E}\xrightarrow{b^{r}}s_{B}s_{B},\quad
s_{A}s_{B}\xrightarrow{p^{r}_{1}}s_{A}s_{A},\quad
s_{A}s_{B}\xrightarrow{p^{r}_{2}}s_{A}s_{E},\end{split}$

where $b^{r},p^{r}_{1},p^{r}_{2}$ are birth, and competition rates, and when
interaction occurs among two empty components no change is accounted.
Otherwise, with probability $1-\mu$, we assume that the sampled individual
changes according to death events with rates $d^{r}_{1},d^{r}_{2}$, as follows

(27) $\begin{split}s_{A}\xrightarrow{d^{r}_{1}}s_{E},\quad
s_{B}\xrightarrow{d^{r}_{2}}s_{E}.\end{split}$

We suppose that if the selected component is empty then no changes in the
populations sizes happen. We assume that at most one birth/competition and one
death event can occur at each time step.

We focus on the evolution of the total number of predators and preys
considering the transition rates from state ${\bf x}=(A,B,E)$ to the states
${\bf x}+{\bf v}_{j}$, where ${\bf v}_{j}$ denotes the $j$-th row of the
stoichiometry matrix

(28) $V=\begin{bmatrix}0&1&-1\cr 1&-1&0\cr 0&-1&1\cr-1&0&1\cr
0&-1&1\end{bmatrix},$

whose columns represent predators, preys and empty components respectively and
whose rows represent the changes in the populations sizes due to the birth,
competition and death events described in (26)-(27). For such predator-prey
dynamics the transition rates write as follows

(29) $\begin{split}&\pi_{{\bf v}_{1}}({\bf x})=2\mu
b^{r}\frac{B}{N}\frac{E}{N-1},\quad\pi_{{\bf v}_{2}}({\bf x})=2\mu
p_{1}^{r}\frac{A}{N}\frac{B}{N-1},\quad\pi_{{\bf v}_{3}}({\bf x})=2\mu
p_{2}^{r}\frac{A}{N}\frac{B}{N-1},\cr&\pi_{{\bf v}_{4}}({\bf
x})=(1-\mu)d_{1}^{r}\frac{A}{N},\quad\pi_{{\bf v}_{5}}({\bf
x})=(1-\mu)d_{2}^{r}\frac{B}{N},\end{split}$

where we define the operators $\pi_{{\bf v}_{j}}({\bf x})=\pi({\bf{x}}+{\bf
v}_{j}|{\bf{x}})$ for any $j=1,\ldots,M$. Here, $M=5$ denotes the total number
of events described in (26)-(27). In Figure 22 we depicted for $A$ predators,
$B$ preys and $E$ empty states all possible states ${\bf x}+{\bf v}_{j}$ after
one step of process (26)-(27). Neighbor states ${\bf x}+{\bf v}_{j}$ are
connected to state ${\bf x}$ either with outgoing arrows (red) for loss
events, or with incoming arrows (blue) for gain events.

$(A,B,E)$$(A,B+1,E)$$(A,B-1,E+1)$$(A+1,B,E-1)$$(A-1,B,E+1)$$(A+1,B-1,E)$$(A-1,B+1,E)$
Figure 22. Neighbor states are connected to state $(A,B,E)$ either with
outgoing arrows (red) for loss events, or with incoming arrows (blue) for gain
events.

Given the probability $P({\bf x},t)$ to be in the state ${\bf x}$ with $A$
predators, $B$ preys and $E$ empty spaces at time $t$, its time evolution is
governed by the associated master equation as follows

(30) $\frac{dP({\bf x},t)}{dt}=\sum_{j=1}^{M}\Biggl{[}(\mathcal{E}^{{\bf
v}_{j}}-1)\left(P({\bf x},t)\pi({\bf x}+{\bf v}_{j}|{\bf x})\right)\Biggr{]},$

where the step operator $\mathcal{E}^{{\bf v}}$ acts on an arbitrary function
$f$ according to

(31) $\mathcal{E}^{{\bf v}_{j}}\left(f({\bf x})\right)=f({\bf x}-{\bf
v}_{j}).$

#### Mean-field type approximations

In analogy to Section 2, the spatial homogeneous mean-field approximation (11)
can be derived from the master equation (30) introducing the empirical
densities $f^{N},g^{N}$ as in (8). A different approximation is obtained
considering a Taylor expansion of (31) up to second order. In this framework,
the master equation (30) is approximated by the following Fokker-Planck
equation

(32) $\displaystyle\frac{dP({\bf
x},t)}{dt}=-\sum_{i=1}^{3}\partial_{x_{i}}\Biggl{(}\sum_{j=1}^{M-1}A_{ij}({\bf
x})P({\bf
x},t)\Biggr{)}+\frac{1}{2}\sum_{i=1}^{3}\partial_{x_{i}x_{i}}\Biggl{(}\sum_{j=1}^{M-1}D_{ij}({\bf
x})P({\bf x},t)\Biggr{)}+$
$\displaystyle\qquad\qquad\quad\sum_{i=1}^{3}\sum_{k=1\atop k\neq
i}^{3}\partial_{x_{i}x_{k}}\Biggl{(}\sum_{j=1}^{M-1}C_{ikj}({\bf x})P({\bf
x},t)\Biggr{)},$

where

(33) $\begin{split}&A_{ij}({\bf x})={\bf v}_{j}^{i}\Pi_{{\bf v}_{j}}({\bf
x}),\quad D_{ij}({\bf x})=({\bf v}_{j}^{i})^{2}\Pi_{{\bf v}_{j}}({\bf
x}),\quad C_{ikj}({\bf x})={\bf v}_{j}^{i}{\bf v}_{j}^{k}\Pi_{{\bf
v}_{j}}({\bf x}).\end{split}$

Here ${\bf v}_{j}^{i}$ denotes the $(i,j)$ component of the stoichiometry
matrix $V$ and $\Pi_{{\bf v}_{l}}({\bf x})=\pi_{{\bf v}_{l}}({\bf x})$ for
$l\in\\{1,2,4\\}$ and $\Pi_{{\bf v}_{3}}({\bf x})=\pi_{{\bf v}_{3}}({\bf
x})+\pi_{{\bf v}_{5}}({\bf x})$. As stated in [17], [36], [16], equation (32)
is the Fokker-Planck equation associated to the Langevin equations

(34) $\frac{dx_{i}}{dt}=\sum_{j=1}^{M-1}\left(A_{ij}({\bf
x})+\sqrt{D_{ij}({\bf x})}\xi_{j}\right),$

for $i=1,\ldots,3$, where $\xi_{j}$, $j=1,\ldots,M-1$ are normally distributed
random numbers with zero mean and variance equal to $1/\sqrt{N}$ . Hence,
multiplying both sides of equation (34) by $P({\bf x},t)/N$ and summing over
all the values of $x_{i}$, $i=1,\ldots,3$ we can write

(35) $\displaystyle\frac{df^{N}}{d\tau}$ $\displaystyle=(\beta
g^{N}-\delta)f^{N}+\sqrt{2\tilde{p}_{1}^{r}f^{N}g^{N}}\xi_{2}-\sqrt{\tilde{d}_{1}^{r}f^{N}}\xi_{4},$
$\displaystyle\frac{dg^{N}}{d\tau}$
$\displaystyle=rg^{N}\left(1-\frac{g^{N}}{K}\right)-\alpha
f^{N}g^{N}+\sqrt{2\tilde{b}^{r}g^{N}(1-f^{N}-g^{N})}\xi_{1}$
$\displaystyle\qquad\qquad\qquad-\sqrt{2\tilde{p}_{2}^{r}f^{N}g^{N}+\tilde{d}_{2}^{r}}\xi_{3}-\sqrt{2\tilde{p}_{1}^{r}f^{N}g^{N}}\xi_{2},$

where we have substituted the value $A_{ij}$ and $D_{ij}$ defined in (33).
Time is scaled according to $\tau=t/N$, and parameters are defined as
$\beta=2\mu p^{r}_{1}$, $\delta=(1-\mu)d^{r}_{1}$,
$\alpha=2\mu(p^{r}_{1}+p^{r}_{2}+b^{r})$,$r=2\mu{b^{r}}-(1-\mu){d^{r}_{2}}$,
and $K=1-\frac{(1-\mu){d^{r}_{2}}}{2\mu{b^{r}}}$. Note that in the limit
$N\to\infty$ the Langevin equations in (35) collapse to the spatial
homogeneous mean-field equations (11).

##  Appendix B Monte Carlo algorithms for agent-based dynamics

We report some of the classic stochastic algorithms for agent-based models
with predator-prey interaction without spatial interaction. We refer to [15,
18, 9] for classical approaches, and to [23] for a comprehensive collection on
stochastic algorithms for agent-based dynamics.

#### Direct method

The main idea of the direct method is to estimate which is the next firing
event $j$ and at which time $\tau$ it will occur from the probability density
function

(36) $p(\tau,j|{\bf x},t)=a_{j}({\bf x})e^{-a_{0}({\bf x})\tau}.$

Here $a_{j}({\bf x})=N\pi({\bf x}+{\bf v}_{j}|{\bf x})$ represents the
propensity of the event $j$ and $a_{0}$ is the sum of all the propensities. At
each time $t$ suppose to split the probability density function defined in
equation (36) into the product of two probability functions, one for the
firing time $\tau$ and the other for the event $j$ that occurs at time
$t+\tau$, as follows

(37) $p(\tau,j|{\bf x},t)=p_{1}(\tau|{\bf x},t)p_{2}(j|\tau,{\bf x},t),$

where

(38) $p_{1}(\tau|{\bf x},t)=a_{0}({\bf x})e^{-a_{0}({\bf x})\tau},\quad
p_{2}(j|\tau,{\bf x},t)=\frac{a_{j}({\bf x})}{a_{0}({\bf x})}.$

Consider two uniformly distributed random numbers $r_{1},r_{2}$ $\sim U(0,1)$,
the index of the next firing event $j$ can be computed as the smallest index
such that

(39) $\sum_{k=1}^{j}a_{k}({\bf x})\geq r_{1}a_{0}({\bf x}).$

The time in which the next event happens can be computed from the first
equation in (38) as

(40) $\tau=\frac{1}{a_{0}({\bf x})}\ln(r_{2}^{-1}).$

Algorithm 3 outlines the details of the direct method.

######  Algorithm 3 (Direct method).

  1. 1.

Define the stoichiometry matrix $V$, the initial state ${\bf x}$, the initial
time $t=0$ and the final time $T$.

  2. 2.

while $t<T$

    1. (a)

Compute the propensities $a_{j}$, $j=1,\ldots,M$.

    2. (b)

Consider two uniformly distributed random numbers $r_{1},r_{2}$ $\sim U(0,1)$.

    3. (c)

Select the index of the next firing event as in equation 39.

    4. (d)

Compute the time in which the next event fires as in equation 40.

    5. (e)

Set ${\bf x}={\bf x}+{\bf v}_{j}$ where ${\bf v}_{j}$ is the $j$-th row of the
stoichiometry matrix $V$ defined as in equation 28.

    6. (f)

Set $t\leftarrow t+\tau$.

repeat

#### Classic Monte Carlo algorithm

The classic Monte Carlo method is an approximated stochastic algorithm. It
allows multiple events to occur at the same time step that is supposed to be
inversely proportional to the sample size $N$. Its computational costs
increase as the sample size increases and can be comparable with the ones of
exact algorithms. The main idea of the classic Monte Carlo algorithm is the
following. Suppose to divide a priori the time interval considering a constant
time step $\tau/N$. At each time $t$ select two components and assume that
with probability $\mu\in[0,1]$ a birth or a competition event can occur. If
the two selected components are

  * •

a prey and an empty space: with probability $b_{r}\tau$ a new prey born and
occupies the empty space;

  * •

a predator and a prey: with probability $p_{1}^{r}\tau$ a new predator born
and with probability $p_{2}^{r}\tau$ the prey dies and a new empty space is
added to the system.

Extract another component and assume that with probability $1-\mu$ a death
event can occur. If the selected component is occupied by a predator then with
probability $d_{1}^{r}\tau$ the predator dies and if it is occupied by a prey
then with probability $d_{2}^{r}\tau$ the prey dies. Algorithm 4 defines the
details of the classic Monte Carlo algorithm.

######  Algorithm 4 (Classic Monte Carlo algorithm).

  1. 1.

Define the sample, the initial time $t=0$, the final time $T$, the time step
as $\tau/N$ and a parameter $\mu\in[0,1]$.

  2. 2.

while $t<T$

    1. (a)

Select two components inside the sample.

    2. (b)

Assume that with probability $\mu$ birth and competition events happen

$BE\xrightarrow{b}BB,\quad AB\xrightarrow{p_{1}}AA,\quad
AB\xrightarrow{p_{2}}AE,$

where $b=b^{r}\tau$, $p_{1}^{r}\tau$ and $p_{2}^{r}\tau$.

    3. (c)

Select another component in the sample.

    4. (d)

Assume that with probability $1-\mu$ death events happen

$A\xrightarrow{d_{1}}E,\quad B\xrightarrow{d_{2}}E,$

with $d_{1}=d_{1}^{r}\tau$ and $d_{2}=d_{2}^{r}\tau$.

    5. (e)

Update the sample.

    6. (f)

Set $t\leftarrow t+\tau/N$.

repeat

####  $\tau$-leaping method

The $\tau$-leaping method is another approximated stochastic algorithm. Its
main idea is to assume that multiple events can happen in the same time
interval that is suppose to vary in time. Therefore this algorithm, as the
Monte Carlo and the direct method, has to cope with high computational costs.
In the following we will describe more in details how the $\tau$-leaping
algorithm works. Assume that more than one event can fire simultaneously at
any time step. The probability that $k_{j}$ events fires in the time interval
$(t,t+\tau)$ follows a Poisson distribution of parameter $a_{j}({\bf x})\tau$
where $a_{j}({\bf x})$ is the propensity of event $j$ and $\tau$ the time
step. The number of firing events $k_{j}$ is generated by sampling its
corresponding Poisson distribution, $Poi(a_{j}({\bf x})\tau)$. The state can
be updated according to the following rule

(41) ${\bf x}(t+\tau)={\bf x}(t)+\sum_{j=1}^{M}k_{j}{\bf v}_{j}={\bf
x}+\sum_{j=1}^{M}Poi(a_{j}({\bf x})\tau){\bf v}_{j}.$

The time step $\tau$ is updated according to a leap selection:

  * •

choose an initial time step $\tau$;

  * •

let $0<\epsilon<1$ and until $\Delta a_{j}({\bf x}):=|a_{j}({\bf
x}(t+\tau))-a_{j}({\bf x}(\tau))|>\epsilon$ set $\tau=\tau/2$.

Algorithm 5 outlines the details of the $\tau$-leaping method.

######  Algorithm 5 ($\tau$ \- leaping method).

  1. 1.

Define the stoichiometry matrix $V$, the initial state ${\bf x}$, the initial
time $t=0$ and the final time $T$, the time step $\tau$.

  2. 2.

while $t<T$

    1. (a)

Compute the propensities $a_{j}$, $j=1,\ldots,M$.

    2. (b)

Consider $M$ uniformly distributed random numbers $k_{1},\ldots,k_{M}$ $\sim
Poi(a_{j}({\bf x})\tau)$.

    3. (c)

Update the state as in equation (41).

    4. (d)

Choose $\tau$ with leap selection.

    5. (e)

Set $t\leftarrow t+\tau$.

repeat

## Acknowledgment

GA and FF would like to thank the Italian Ministry of Instruction, University
and Research (MIUR) to support this research with funds coming from PRIN
Project 2017 (No. 2017KKJP4X entitled “Innovative numerical methods for
evolutionary partial differential equations and applications”).

## References

  * [1] G. Albi, N. Bellomo, L. Fermo, S.-Y. Ha, J. Kim, L. Pareschi, D. Poyato, and J. Soler.  Vehicular traffic, crowds, and swarms: From kinetic theory and multiscale methods to applications and research perspectives.  Mathematical Models and Methods in Applied Sciences, 29(10):1901–2005, 2019.
  * [2] L. J. Allen, B. M. Bolker, Y. Lou, and A. L. Nevai.  Asymptotic profiles of the steady states for an sis epidemic reaction-diffusion model.  Discrete & Continuous Dynamical Systems, 21(1):1, 2008.
  * [3] A. Berryman.  Population cycles: the case for trophic interactions.  Oxford University Press, 2002.
  * [4] A. J. Black and A. J. McKane.  Stochastic formulation of ecological models and their applications.  Trends in ecology & evolution, 27(6):337–345, 2012.
  * [5] B. Blasius, L. Rudolf, G. Withoff, U. Gaedke, and G. F. Fussmann.  Long-term cyclic persistence in an experimental predator-prey system.  Nature, 577:226–230, 2020.
  * [6] W. Boscheri, G. Dimarco, and L. Pareschi.  Modeling and simulating the spatial spread of an epidemic through multiscale kinetic transport equations.  Mathematical Models and Methods in Applied Sciences, pages 1–39, 2021.
  * [7] R. E. Caflisch.  Monte Carlo and quasi-Monte Carlo methods.  Acta numerica, 7:1–49, 1998.
  * [8] J. A. Carrillo, Y. Huang, and M. Schmidtchen.  Zoology of a nonlocal cross-diffusion model for two species.  SIAM Journal on Applied Mathematics, 78(2):1078–1104, 2018.
  * [9] S. Chen and U. C. Täuber.  Non-equilibrium relaxation in a stochastic lattice Lotka–Volterra model.  Physical biology, 13(2):025005, 2016.
  * [10] E. Cristiani, B. Piccoli, and A. Tosin.  Multiscale modeling of pedestrian dynamics, volume 12.  Springer, 2014.
  * [11] A. Duncan, R. Erban, and K. Zygalakis.  Hybrid framework for the simulation of stochastic chemical kinetics.  Journal of Computational Physics, 326:398–419, 2016.
  * [12] S. Engblom, L. Ferm, A. Hellander, and P. Lötstedt.  Simulation of stochastic reaction-diffusion processes on unstructured meshes.  SIAM Journal on Scientific Computing, 31(3):1774–1797, 2009.
  * [13] M. B. Flegg, S. Hellander, and R. Erban.  Convergence of methods for coupling of microscopic and mesoscopic reaction–diffusion simulations.  Journal of computational physics, 289:1–17, 2015.
  * [14] J. Friedman and J. Gore.  Ecological systems biology: The dynamics of interacting populations.  Current Opinion in Systems Biology, 1:114–121, 2017.
  * [15] D. T. Gillespie.  A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.  Journal of computational physics, 22(4):403–434, 1976.
  * [16] D. T. Gillespie.  The multivariate langevin and fokker–planck equations.  American Journal of Physics, 64(10):1246–1257, 1996.
  * [17] D. T. Gillespie.  The chemical Langevin equation.  The Journal of Chemical Physics, 113(1):297–306, 2000.
  * [18] D. T. Gillespie.  Stochastic simulation of chemical kinetics.  Annu. Rev. Phys. Chem., 58:35–55, 2007.
  * [19] B. Hauer, J. Maciejko, and J. Davis.  Nonlinear power spectral densities for the harmonic oscillator.  Annals of Physics, 361:148–183, 2015.
  * [20] A. Hellander and P. Lötstedt.  Hybrid method for the chemical master equation.  Journal of Computational Physics, 227(1):100–122, 2007.
  * [21] S. Hellander, A. Hellander, and L. Petzold.  Mesoscopic-microscopic spatial stochastic simulation with automatic system partitioning.  The Journal of chemical physics, 147(23):234101.
  * [22] W. Hundsdorfer and J. G. Verwer.  Numerical solution of time-dependent advection-diffusion-reaction equations, volume 33.  Springer Science & Business Media, 2013.
  * [23] L. Marchetti, C. Priami, and V. H. Thanh.  Simulation algorithms for computational systems biology.  Springer, 2017.
  * [24] The Mathworks, Inc., Natick, Massachusetts.  MATLAB version 9.3.0.713579 (R2017b), 2017.
  * [25] S. J. Maynard.  Models in ecology.  Cambridge UP, Cambridge, 1974.
  * [26] A. J. McKane and T. J. Newman.  Stochastic models in population biology and their deterministic analogs.  Physical Review E, 70(4):041902, 2004.
  * [27] A. J. McKane and T. J. Newman.  Predator-prey cycles from resonant amplification of demographic stochasticity.  Physical review letters, 94(21):218102, 2005.
  * [28] J. D. Murray.  Mathematical biology: I. An introduction, volume 17.  Springer Science & Business Media, 2007.
  * [29] R. M. Nisbet and W. Gurney.  Modelling fluctuating populations: reprint of first Edition (1982).  2003\.
  * [30] L. Pareschi and G. Toscani.  Interacting multiagent systems: kinetic equations and Monte Carlo methods.  OUP Oxford, 2013.
  * [31] E. C. Pielou et al.  An introduction to mathematical ecology.  An introduction to mathematical ecology., 1969.
  * [32] E. Renshaw.  Modelling biological populations in space and time, volume 11.  Cambridge University Press, 1993.
  * [33] T. Roose, S. J. Chapman, and P. K. Maini.  Mathematical models of avascular tumor growth.  SIAM review, 49(2):179–208, 2007.
  * [34] N. Shigesada.  Spatial distribution of dispersing animals.  Journal of mathematical biology, 9(1):85–96, 1980.
  * [35] F. Spill, P. Guerrero, T. Alarcon, P. K. Maini, and H. Byrne.  Hybrid approaches for multiple-species stochastic reaction–diffusion models.  Journal of computational physics, 299:429–445, 2015.
  * [36] R. Toral and P. Colet.  Stochastic numerical methods: an introduction for students and scientists.  John Wiley & Sons, 2014.
  * [37] N. G. Van Kampen.  Stochastic processes in physics and chemistry, volume 1.  Elsevier, 1992.
  * [38] G. Webb.  A reaction-diffusion model for a deterministic diffusive epidemic.  Journal of Mathematical Analysis and Applications, 84(1):150–161, 1981.