7 Stable oscillations of a predator-prey probabilistic
cellular automaton: a mean-field approach.
Tânia Tomé and Kelly C de Carvalho
Instituto de F́ısica, Universidade de São Paulo
Caixa Postal 66318
05315-970 São Paulo, São Paulo, Brazil
E-mail: ttome@if.usp.br
Abstract.
We analyze a probabilistic cellular automaton describing the dynamics of
coexistence of a predator-prey system. The individuals of each species are localized
over the sites of a lattice and the local stochastic updating rules are inspired on
the processes of the Lotka-Volterra model. Two levels of mean-field approximations
are set up. The simple approximation is equivalent to an extended patch model,
a simple metapopulation model with patches colonized by prey, patches colonized
by predators and empty patches. This approximation is capable of describing the
limited available space for species occupancy. The pair approximation is moreover
able to describe two types of coexistence of prey and predators: one where population
densities are constant in time and another displaying self-sustained time-oscillations
of the population densities. The oscillations are associated with limit cycles and arise
through a Hopf bifurcation. They are stable against changes in the initial conditions
and, in this sense, they differ from the Lotka-Volterra cycles which depend on initial
conditions. In this respect, the present model is biologically more realistic than the
Lotka-Volterra model.
PACS numbers: 87.23.Cc, 05.65.+b, 05.70.Ln, 02.50.Ga
http://arxiv.org/abs/0704.0512v1
Stable oscillations of a predator-prey probabilistic cellular automaton 2
1. Introduction
The simplest model exhibiting time-oscillations in a two-component system is the model
proposed independently by Lotka [1, 2, 3] and by Volterra [4]. In this model the
individuals of two species are dispersed over an assumed homogeneous space. It is
implicitly assumed in this approach that any individual can interact with any other one
with equal intensity implying that their positions are not taken into account. The time
evolution of the densities of the two species in the Lotka-Volterra model is given by a set
of two ordinary differential equations [5, 6, 7, 8] and is set up in analogy with the laws
of mass-action. Depending on the level of description wanted, the approach based on
mass-action laws, contained on the Lotka-Volterra model, suffices. However, there are
situations in which the coexistence takes place in a spatially heterogeneous habitat such
that the population densities can be very low in some regions. In this case we need to
proceed beyond the mass-law equations and consider the space structure of the habitat.
In other words, it becomes necessary to analyze the coexistence by taking explicitly into
account spatial structured models.
In fact, the role of space in the description of population biology problems has been
recognized by several authors in the last years [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26]. In a very clear manner, Durrett and Levin [11] have pointed
out that the modelling of population dynamics systems which are spatially distributed
by interacting particle systems [11, 27, 28, 29] is the appropriate theoretical approach
that is able to give the more complete description of the problem. We include in this
approach probabilistic cellular automata (PCA) [29, 30, 31], which will concern us here.
We refer to interacting particle systems and PCA as stochastic lattice models. They are
both Markovian processes defined by discrete stochastic variables residing on the sites
of a lattice; the former being a continuous time process and the latter a discrete time
process.
In the present work we study the coexistence and the emergence of stable self-
sustained oscillations in a predator-prey system by considering a PCA previously studied
by numerical simulations [24, 26]. This PCA is defined by local rules, similar to the
ones of the contact process [27], that are capable of describing the interaction between
prey and predator. Here, we focus on the analysis of the PCA by means of dynamic
mean-field approximations [10, 28, 29, 30, 32, 33]. In this approach the equations for
the time evolution of correlations of various orders are truncated at a certain level and
high order correlations of sites are written in terms of small order correlations. The
simplest approximation is the one in which all correlations are written in terms of one-
site correlation, called simple approximation. In a more sophisticated approximation,
called pair approximation [10, 29], any correlation is written in terms of one-site and
two-site correlations.
The simple mean-field approximation is capable of predicting coexistence of
individuals in a stationary state where the densities of each species, and of empty
sites, are constant. However, it is not capable of predicting possible time oscillating
Stable oscillations of a predator-prey probabilistic cellular automaton 3
behavior of the population densities and we have proceed to the next order of mean-
field approximation. The simple approximation, on the other hand, can be placed in an
explicit correspondence with a patch model [7, 12, 34] where unoccupied patches can be
colonized by prey and patches occupied by prey can be colonized by predators that in
turn may become extinct. In this approximation the PCA can be seen as an extended
version of the Lotka-Volterra model which includes an extra logistic term related to the
empty sites.
The pair-mean field approximation is able to predict possible time oscillating
behavior of the population densities that are self-sustained and are attained thorough
Hopf bifurcations. This is in contrast with the Lotka-Volterra model which presents
no stable oscillations but exhibits instead infinite cycles that are associated to different
initial conditions. However, from the biological point of view, one does not expect that a
small variation in the initial densities of prey and predator result in different amplitudes
of oscillations. Within our approach, a PCA treated in the pair-approximation, the
oscillations are associated with limit cycles what mean to say that they are stable
against the changes in the initial conditions. According to our point of view, the pair-
approximation, in which the correlation between neighboring sites are treated exactly,
provides a basic description of the predator-prey spatial interactions. For this reason,
we will refer to the PCA in this approximation as a quasi-spatial-structured model.
2. Model
2.1. Probabilistic cellular automaton
We consider interacting particles living on the sites of a lattice and evolving in time
according to Markovian local rules. The lattice is the geometrical object that plays the
role of the spatial region occupied by particles, in a general case, or by individuals of each
species in the present case. The lattice sites are the possible locations for the individuals.
Each site can be either empty or occupied by one individual of different species and a
stochastic variable ηi is introduced to describe the state of each site at a given instant
of time. The state of the entire system is denoted by η = (η1, . . . , ηi, . . . , ηN) where
N is the total number of sites. The transition between the states is governed by the
interactions between neighbor sites in the lattice and by a synchronous dynamics.
The probability P (ℓ)(η) of configuration η at time step ℓ evolves according to the
Markov chain equation
P (ℓ+1)(η) =
W (η|η′)P (ℓ)(η′), (1)
where the summation is over all the microscopic configurations of the system, and
W (η|η′) is the conditional transition probability from state η′ at time ℓ to state η at
time ℓ + 1. This transition probability does not depend on time and contains all the
information about the dynamics of the system. Taking into account that all the sites
are simultaneously updated, which is the fundamental property of a PCA, the transition
Stable oscillations of a predator-prey probabilistic cellular automaton 4
Figure 1. Transitions of the predator-prey model. The three states are: prey or
herbivorous (H or 1), predator (P or 2) and empty or vegetable (V or 0). The allowed
transitions obey the cyclic order shown.
probability can be factorized and written in the form [29, 30]
W (η|η′) =
wi(ηi|η
′), (2)
where wi(ηi|η
′) is the conditional transition probability that site i takes the state ηi
given that the whole system is in state η′. Being a probability distribution, the quantity
wi(ηi|η
′) must satisfy the following properties: wi(ηi, η
′) ≥ 0 and
wi(ηi|η
′) = 1. (3)
The average of any state function F (η) is evaluated by
〈F (η)〉ℓ =
F (η)P (ℓ)(η). (4)
The time evolution equation for 〈F (η)〉 is obtained from definition (4) and equation (1).
For example, we can derive the equations for the time evolution of densities and two-site
correlations.
2.2. Predator-prey probabilistic cellular automaton
To model a predator-prey system by a PCA, the stochastic variable ηi associated to site
i will represent the occupancy of the site by one prey, or the occupancy by one predator
or the vacancy (a site devoid of any individual). The variable ηi is assumed to take
the value 0, 1, or 2, according to whether the site is empty (V), occupied by a prey
individual (H) or by a predator (P), respectively. That is,
0, empty (V),
1, prey (H),
2, predator (P),
which defines a three-state per site PCA.
The stochastic rules, embodied in the transition rate wi(ηi|η
′), are set up according
to the following assumptions. (a) The space is homogeneous, which means to say that
no region of the space will be privileged against the others, that is, in principle the
individuals have the same conditions of surveillance in any space region. (b) The space
is isotropic, which means to say that there is no preferential direction in this space for
any interaction. (c) The allowed transitions between states are only the ones that obey
Stable oscillations of a predator-prey probabilistic cellular automaton 5
the cyclic order shown in figure 1. Prey can only born in empty sites; prey can give place
to a predator, in a process where a prey individual dies and a predator is instantaneously
born; finally a predator can die leaving an empty site. The empty sites are places where
prey can proliferate and can be seen as the resource for prey surveillance. The death of
predators complete this cycle, reintegrating to the system the resources for prey.
The predator-prey PCA has three parameters: a, the probability of birth of prey, b,
the probability of birth of predator and death of prey, and c, the probability of predator
death. Two of the process are catalytic: the occupancy of a site by prey or by a predator
is conditioned, respectively, to the existence of prey or predator in the neighborhood of
the site. The third reaction, where predator dies, is spontaneous, that is, it occurs, with
probability c, independently of the neighbors of the site. We assume that
a+ b+ c = 1, (6)
with 0 ≤ a, b, c ≤ 1.
The transition probabilities of the predator-prey PCA are described in what follows:
(a) If a site i is empty, ηi = 0, and there is at least one prey in its first neighborhood
there is a favorable condition for the birth of a new prey. The probability of site i being
occupied in next time step by a prey is proportional to the parameter a and to the
number of prey that are in the first neighborhood of the empty site.
(b) If a site is occupied by a prey, ηi = 1, and there is at least one predator in its
first neighborhood then the site has a probability of being occupied by a new predator
in the next instant of time. In this process the prey dies instantaneously. The transition
probability is proportional to the parameter b and the number of predators in first
neighborhood of the site.
(c) If site i is occupied by a predator, ηi = 2, it dies with probability c.
The transition probabilities associated to the three processes above mentioned can
be summarized as follows:
wi(0|η) = cδ(ηi, 2) + [1− fi(η)]δ(ηi, 0), (7)
wi(1|η) = fi(η)δ(ηi, 0) + [1 − gi(η)]δ(ηi, 1), (8)
wi(2|η) = gi(η)δ(ηi, 1) + (1− c)δ(ηi, 2), (9)
where
fi(η) =
δ(ηk, 1), gi(η) =
δ(ηk, 2), (10)
and the summation is over the four nearest neighbors of site i in a regular square lattice.
The notation δ(x, y) stands for the Kronecker delta function. These stochastic local
rules, when inserted in equation (2), define the dynamics of the PCA for a predator-
prey system.
The present stochastic dynamics predicts the existence of states, called absorbing
states, in which the system becomes trapped. Once the system has entered such a state
Stable oscillations of a predator-prey probabilistic cellular automaton 6
it cannot escape from it anymore remaining there forever. There are two absorbing
states. One of them is the empty lattice. Since the predator death is spontaneous, a
configuration where just predators are present is not stationary. This situation happens
whenever the prey have been extinct. In this case the predator cannot reproduce
anymore and also get extinct, leaving the entire lattice with empty sites. The other
absorbing state is the lattice full of prey. This situation occurs if there are few predators
and they become extinct. The remaining prey will then reproduce without predation
filling up the whole lattice. The existence of absorbing stationary states is an evidence
of the irreversible character of the model or, in other words, of the lack of detailed
balance [29]. However, the most interesting states, the ones that we are concerned with
in the present study, are the active states characterized by the coexistence of prey and
predators.
2.3. Time evolution equations for state functions
We start by defining the densities, which are the one-site correlations, and the two-site
correlations. These quantities will be useful in our mean-field analysis to be developed
below. The density of prey, predator, and empty sites at time step ℓ are defined thought
the expressions
i (1) = 〈δ(ηi, 1)〉ℓ, (11)
i (2) = 〈δ(ηi, 2)〉ℓ, (12)
i (0) = 〈δ(ηi, 0)〉ℓ. (13)
The evolution equations for the above densities are obtained from their definitions as
state functions, as given by equation (4), and by using the evolution equation for P (ℓ)(η),
given by equation (1). The resulting equations can be formally written as
(ℓ+1)
i (1) = 〈wi(1|η)〉ℓ, (14)
(ℓ+1)
i (2) = 〈wi(2|η)〉ℓ, (15)
(ℓ+1)
i (0) = 〈wi(0|η)〉ℓ, (16)
where the transition probabilities for this model are given in equations (7), (8) and (9).
The correlation between a prey localized at site i and a predator localized at site j
at time step ℓ is defined by
ij (1, 2) = 〈δ(ηi, 1)δ(ηj, 2)〉ℓ. (17)
The other two-site correlations are defined similarly. The time evolution equation for
the correlation of two neighbor sites i and j, one being occupied by a prey and the other
by a predator, is given by
(ℓ+1)
ij (1, 2) = 〈wi(1|η)wj(2|η)〉ℓ. (18)
The other two-site evolution equations are given by similar formal expressions. We
can also derive equations for three-site correlations. Since we are interested here on
Stable oscillations of a predator-prey probabilistic cellular automaton 7
approximations in which only the one-site and two-site correlations should be treated
exactly, the above equations suffice.
We call the attention to the fact that equation (18) includes the product of two
transition probabilities. This is a consequence of the synchronous update of the PCA
which allows that both neighboring sites i and j have their states changed at same
time step. This situation does not occur when we consider a continuous time one-
site dynamics. Therefore, although local interaction in the present PCA and in the
continuous time model considered in reference [10] are the same, the predator-prey
system evolves according to different global dynamics which leads to different time
evolution equations for the densities and the correlations.
The exact evolution equations for the one-site correlations are
P ′j (1) =
Pji(01)−
Pji(12) + Pj(1), (19)
P ′j (2) =
Pji(12) + (1− c)Pj(2), (20)
where the summation in j is over the ζ nearest neighbors of site i. To simplify notation
we are using unprimed and primed quantities to refer to quantities taken at time ℓ and
ℓ+ 1, respectively.
The exact evolution equations for the correlations of two nearest neighbor sites j
and k are
P ′jk(01) =
n(6=j)
Pjkn(001)−
i(6=k)
Pijkn(1001)
+ (1−
Pjk(01)−
n(6=j)
Pjkn(012)
i(6=k)
Pijk(101)−
n(6=j)
Pijkn(1012)
)Pjk(21)−
n(6=j)
Pjkn(212)
n(6=j)
Pjkn(201), (21)
P ′jk(12) =
n(6=j)
Pjkn(012) +
i(6=k)
Pijkn(1012)
n(6=j)
Pjkn(112)−
i(6=k)
Pijkn(2112)
+ (1− c)
)Pjk(12)−
i(6=k)
Pijk(212)
(1− c)
i(6=k)
Pijk(102), (22)
Stable oscillations of a predator-prey probabilistic cellular automaton 8
P ′jk(02) =
n(6=j)
)Pjkn(012)−
i(6=k)
Pijkn(1012)
+ (1− c)
cPjk(22) + Pjk(02)−
i(6=k)
Pijk(102)
Pjk(21) +
n(6=j)
Pjkn(212)
 , (23)
where the summation in i is over the nearest neighbors of j and the summation in n is
over the nearest neighbors of k.
3. Mean-field approximation
3.1. One and two site approximations
The evolution equation for a density in any interacting particle system which evolves in
time according to local interaction rules always contains terms related to the correlations
between neighbor sites in a lattice. The evolution equations for the correlations of two
neighbor sites includes the correlation of clusters of three or more sites in the lattice and
so on. In this way we can have an infinite set of coupled equations for the correlations
which is equivalent to the evolution equation for the probability P (ℓ)(η), described in
equation (1) for the automaton. The scope of the dynamic mean-field approximation
consists in the truncation of this infinite set of coupled equations [30, 31, 32, 33].
The lowest order dynamic mean-field approximation is the one where the probability
of a given cluster is written as the product of the probabilities of each site. That is, all
the correlations between sites in the cluster are neglected. For example, let us consider
the cluster constituted by a center (C) site and its first neighboring sites to the north
(N), south (S), east (E) and west (W) as shown in figure 2.
Within the one-site approximation the probability P (N,E,W, S, C) corresponding
to the cluster shown in figure 2 is approximated by
P (N,E,W, S, C) = P (N)P (E)P (W )P (S)P (C), (24)
where P (X), X = N,E,W, S, C are the one-site probabilities corresponding to each
site. For some stochastic dynamics models this approximation is able to give qualitative
results that are in agreement with the expected results.
In order to get a better approximation we must include fluctuations. The
simplest mean-field approximation that includes correlations is the pair-mean field
approximation. This approximation is better explained by taking again, as an example,
the cluster constituted by a center site which and its four nearest neighbors, shown
above. Within the pair-approximation the conditional probability P (N,E,W, S |C) is
approximated by
P (N,E,W, S |C) = P (N |C)P (E, |C)P (W |C)P (S |C), (25)
Stable oscillations of a predator-prey probabilistic cellular automaton 9
Figure 2. A site (C) of the square lattice and its four nearest neighbor sites (N, E,
W, S).
that is, the conditional probability P (N,E,W, S |C) is written in terms of the product
of the conditional probabilities P (X|C), X = N,E,W, S. Now using the definition of
conditional probability we have
P (N,E,W, S, C)
P (C)
P (N,C)
P (C)
P (E,C)
P (C)
P (W,C)
P (C)
P (S, C)
P (C)
, (26)
P (N,E,W, S, C) =
P (N,C)P (E,C)P (W,C)P (S, C)
[P (C)]3
. (27)
We see that the resulting probability is written as a function of two-site correlations
P (X,C), and the one-site correlation P (C).
3.2. Patch model
The simple mean-field approximation of the predator-prey PCA describes exactly
the same properties of an extended Levins patch model [7, 34]. That is, the PCA
with local rules similar to the contact process becomes, in the simple mean-field
approximation, analogous to the Levins model for metapopulation with empty patches,
patches colonized by prey and patches colonized by predators.
In the one-site mean-field approximation we consider that the probability of any
cluster of sites can be written as the product of the probabilities of each site, as in
equation (24). Using this approach, and writing x = Pi(1), y = Pi(2), and z = Pi(0) it
can be seen that the set of equations can be reduced to the following two-dimensional
map [26]
x′ = x+ axz − bxy, (28)
which is an evolution equation for prey density x, and
y′ = y + bxy − cy, (29)
which is an evolution equation for predator density yℓ. Notice that
z = 1− x− y. (30)
The fixed point of this map are those that represent the stationary solutions x′ = x
and y′ = y, and they correspond to the three following solutions x1 = 0, y1 = 0, and
x2 = 1, y2 = 0, and x3 = a/b, y3 = (1 − c/b)/(1 + b/a). The first solution corresponds
to an absorbing states where both species have been extinct. The second corresponds
Stable oscillations of a predator-prey probabilistic cellular automaton 10
−0.5 0 0.5
Figure 3. Phase diagram of the patch model. The continuous line represents
the transition, c1(p), between the prey absorbing (A) state and the active species
coexistence (C) state. The dashed line separates the two asymptotic time behavior of
the active state.
to an absorbing state where predators have extinct. The third solution corresponds to
an active state where prey and predator coexist.
Due to the constraint (6), the parameters a, b and c are not all independent and
only two can be chosen as independent. For this reason it is convenient to introduce the
following parametrization [10]
− p, b =
+ p, (31)
and consider p and c as the independent variables. The parameter p is such that
−1/2 ≤ p ≤ 1/2 and 0 ≤ c ≤ 1 as before. This parametrization will useful in the
determination of the different phases displayed by the model.
A linear stability analysis reveals that solution the (x1, y1) is a hyperbolic saddle
point for any set of the parameters a, b and c and so it is always unstable. The empty
absorbing state will never be reached. A linear stability analysis also shows that the
solution (x2, y2) is a stable node in the following region of the phase diagram c > c1
where
c1(p) =
(1 + 2p). (32)
The active solution is stable in the region c < c1 and is attained in two ways: by an
asymptotic stable focus, where the successive interactions of the map show damped
oscillations; or trough an asymptotic stable node. In the phase diagram of figure 3 we
show the transition line between the prey absorbing state and the active state given by
c = c1.
In figure 4 it is shown the behavior of the densities against the parameter c, the
probability of predators death, for the special case p = 0.2. In terms of phase transitions
what happens is that in the phase diagram there is a transition line separating the
absorbing prey phase and the active phase which is characterized by constant and
nonzero densities of prey and predator.
Stable oscillations of a predator-prey probabilistic cellular automaton 11
0 0.1 0.2 0.3 0.4 0.5 0.6
predator
Figure 4. Densities of predator and prey as functions of the parameter c for p = 0.2,
for the patch model.
We may conclude that the mean-field approximation for the predator-prey
probabilistic cellular automaton with rules (7), (8), and (9) is capable to show, under a
robust set of control parameters, that prey and predators can coexist without extinction.
However the map defined by equations (28) and (29) is not able to describe self-sustained
oscillations of species population densities.
3.3. Quasi-spatial model
In order to find if oscillations in the species populations can be described within a mean-
field approach we consider a more sophisticated approximation, the pair-approximation,
where correlations of two neighbor sites are included in the time evolution equations for
the densities. This is the lowest order mean-field approximation which takes into account
the spatial localization of neighboring individuals.
In this analysis we will maintain the correlations of one site and the correlations
of two-sites in the equations. Correlations of three and four neighbor sites will be
approximated by means of equation (27). With these approximation the model is
described by the following set of five coupled equations
x′ = au− bv + x, (33)
y′ = bv + (1− c)y, (34)
u′ = αa[
] + [(1− βa)− αa
][u− αb
+ αac
+ c[(1− βb)v − αb
], (35)
v′ = αb[βa
] + αa(1− c)
+ αb[
] + (1− c)[(1− βb)v − αb
], (36)
Stable oscillations of a predator-prey probabilistic cellular automaton 12
−0.5 0 0.5
Figure 5. Phase diagram of the quasi-spatial model. The upper continuous
line represents the transition, c1(p), between the prey absorbing (A) state and the
nonoscillating coexistence (CNO) state. The lower continuous line represents the
transition, c2(p), between the nonoscillating coexistence and the oscillating (COS)
coexistence state. The dashed line separates the two asymptotic time behavior of
the nonoscillating coexistence state.
w′ = αb[(1− βa)
] + (1− c)[w − αa
+ c[βbv + αb
] + c(1− c)s, (37)
where α and β are numerical fractions defined by α = (ζ − 1)/ζ and β = 1/ζ where
ζ is the coordination number of the lattice. For the present case of a square lattice,
ζ = 4 so that α = 3/4 and β = 1/4. We are using the following notation: u = P (0, 1),
v = P (1, 2), and w = P (0, 2) and also r = P (1, 1), q = P (0, 0) and s = P (2, 2). The
last three correlations are not independent but are related to others by
r = x− u− v, (38)
q = z − u− w, (39)
s = y − v − w. (40)
We used the properties P (1, 0) = P (0, 1), P (1, 2) = P (2, 1) and P (2, 0) = P (0, 2), that
follows from the assumption that space is isotropic and homogeneous.
We have analyzed numerically the five-dimensional map, described by the set of
equations (33), (34), (35), (36) and (37), and we have obtained four types of solutions.
Two solutions are trivial and are given by x = y = u = v = w = 0 and x = 1,
y = u = v = w = 0. They correspond to the empty and prey absorbing states,
respectively. The empty absorbing state, where both species have been extinct is
an unstable solution and never occurs. However, the prey absorbing state is one of
the possible stable stationary solutions and is stable above the critical transition line
Stable oscillations of a predator-prey probabilistic cellular automaton 13
0 0.1 0.2 0.3 0.4 0.5 0.6
predator
Figure 6. Densities of predator and prey as functions of c for the quasi-spatial model,
for p = 0.2.
c = c1(p) shown in figure 5. Below this line it becomes unstable giving rise to the active
state.
The other solutions correspond to the active states where both prey and predator
coexist. These solutions are of two kinds: a stationary solution where there is a
coexistence of the two species with densities constant in time, which we call the
nonoscillating (NO) active state; and another solution where both population densities
oscillate in time. This solution corresponds to a self-sustained oscillation of the predator-
prey system and will be called the oscillating (O) active state. In the phase diagram of
figure 5 there is a line c = c2(p) that separates the NO and O active phases. Figure 6
shows the behavior of the densities as a function of c for p = 0.
3.4. Oscillatory behavior
In figure 7 we show an example of self-sustained oscillations of the densities of prey
and predators as functions of time. The oscillating solutions are attained from the
nonoscillating solutions by a Hopf bifurcation. The fixed point associated to this solution
is an unstable center which produces a stable limit cycle as trajectories in the phase-
space of the predator density versus prey density, as can be seen in figure 7. Notice
that the oscillations are not damped and have a well defined period which is the same
for the prey density and for the predator density, which implies that the oscillations are
coupled. A maximum of predators always follow a maximum of prey. This means that
the abundance of prey is a condition that favors the increase in the number of predators.
As the predator number increases the prey population decays. The evanescence of prey
is followed by a decrease in the predator number, giving conditions for the increase of
prey population until the cycle starts again.
A well defined oscillatory behavior is found for many biological population, the
most famous being the one related to the time oscillations of the population of lynx
and snowshoe hare in Canada for which data were collected for a long period of time
Stable oscillations of a predator-prey probabilistic cellular automaton 14
1000 1200 1400 1600 1800 2000
predator
0 0.05 0.1 0.15
Prey population density
Figure 7. (a) Densities of predator and prey as functions of time and (b) density of
predator versus density of prey, for the quasi-spatial model, for p = 0 and c = 0.016.
[7, 8]. If the hare population cycles are mainly governed by the lynx cycle then the
oscillations shown by the present model reproduces qualitatively some of the features of
this predator-prey dynamics.
Next we analyze the behavior of the frequency and amplitude of oscillations. Fixing
the parameter p and varying the parameter c, we verify that in all the oscillating region
the frequency of oscillation is proportional to parameter c,
ω ∼ c, (41)
as can be seen in figure 8. Low frequencies are associated to low values of c; what means
that, for small values of c, the greater the predator lifetime the greater will be period of
the oscillation. As to the amplitude A of the oscillations, we have verified, that fixing
the value of p and varying the parameter c, it increases as c decreases. Our results show
that,
A ∼ (c− c2)
1/2, (42)
as expected for a Hopf bifurcation and shown in figure 8. The transition line c = c2
from the oscillating phase to the nonoscillating phase can either be obtained by using the
criterion given by equation (42) or by analyzing the eigenvalues associated to the map
given by the set of equations (33), (34), (35), (36) and (37). This last criterion means
to find the points of phase diagram such that the real part of the dominant complex
eigenvalue equals 1.
4. Discussion and conclusion
The main result coming from the pair mean-field approximation applied to the predator-
prey PCA is that it is possible to describe coexistence and self-sustained time oscillations.
Moreover, these are stable oscillations. Given a set of parameters, just one limit cycle
is achieved, no matter what the initial conditions are. This property is essential in
describing a biological system since a small variation in the initial condition can not
Stable oscillations of a predator-prey probabilistic cellular automaton 15
0 0.005 0.01 0.015 0.02
0 0.005 0.01 0.015 0.02
Figure 8. (a) Frequency of oscillations ω versus the parameter c. The frequency
vanishes linearly as one approaches c = 0. (b) Amplitude A of oscillations versus c
near the Hopf bifurcation point c2 = 0.019. The quantity A
2 vanishes linearly when
c → c2 in accordance with a Hopf bifurcation.
modify the amplitude, frequency and mean value of the time oscillation densities of a
predator-prey system. Similar results were obtained from a continuous time version of
the present model [10]. Although the simple mean-field equations are essentially the
same in both versions this is not the case concerning the pair mean-field approximation.
The time evolutions of the pair correlations for the PCA, presented here, depend on
higher order correlations (up to fourth) when compared to the ones of the continuous
version (up to third).
The model studied here is a spatial structured model with individuals residing in
sites of a lattice and described by discrete dynamic variables. When we perform simple
mean-field approximation we neglect all the correlations of sites in the lattice. But we
take into account that there are limited resources for the surveillance of each species.
For example in the time evolution equation for the density of prey we have an explicit
term relative to reaction of birth of prey which is the product of the density of prey x by
the density of empty sites z = (1−x−y). This coincides with an extended patch model
approach for predator-prey systems. The presence of this term is what differs the simple
mean-field equations from the Lotka-Volterra equations. However, taking into account
the limitation of space and resources the simple mean-field equations are not sufficient
to get self-sustained oscillations although able to describe damped time oscillations of
population densities.
To get self-sustained time oscillations we had to proceed to the next level of
approximation in which a pair of nearest neighbor sites is treated exactly. This
approximation can be seen as representing a pair of nearest neighbor sites immersed
in a mean field produced by the rest of the lattice. The most important feature being
the fact that the two sites of this pair can be seen as localized in space. The set of
five equations which results from the pair approximation for the PCA is indeed able
to produce self-sustained oscillations of population densities. It presents an important
Stable oscillations of a predator-prey probabilistic cellular automaton 16
property that the Lotka-Volterra model lacks, namely, the oscillating solutions are stable
and are unique for a given set of the control parameters.
Acknowledgements
The authors have been supported by the Brazilian agency CNPq.
References
[1] Lotka A 1920 J. Am. Chem. Soc. 42 1595
[2] Lotka A 1920 Proc. Nat. Acad. of Sciences USA 6 410
[3] Lotka A 1924 Elements of Mathematical Biology (new York: Dover)
[4] Volterra V 1931 Leçons sur la Théorie Mathématique de la Lutte pour la Vie Paris: Gauthier-
Villars)
[5] Haken H 1976 Synergetics, An Introduction (Berlin: Springer)
[6] Renshaw E 1991 Modelling Biological Populations in Space and Time (Cambridge: Cambridge
University Press)
[7] Hastings A 1997 Population Biology: Concepts and Models (New York: Springer)
[8] Ricklefs R E and Miller G L 2000 Ecology (New York: Freeman)
[9] Tainaka K 1989 Phys. Rev. Lett. 63 2688
[10] Satulovsky J and Tomé T 1994 Phys. Rev. E 49 5073
[11] Durrett R and Levin S 1994 Theor. Popul. Biol. 46 363
[12] Hanski I and Gilpin M E (eds.) 1997 Metapopulation Biology: Ecology, Genetic and Evolution
(San Diego: Academic Press)
[13] Satulovsky J and Tomé T 1997 J. Math. Biol. 35 344
[14] Tilman D and Kareiva P 1997 Spatial Ecology: the Role of Space in Population Dynamics and
Interactions (Princeton: Princeton University Press)
[15] Fracheburg L and Krapvisky P 1998 J. Phys. A 31 L287
[16] Liu Y C, Durrett R and Milgroom M 2000 Ecol. Model. 127 291
[17] Antal T, Droz M, Lipowsky A and Odor G 2001 Phys. Rev. E 64 036118
[18] Ovaskanien O, Sato K, Bascompte J and Hanski I 2002 J. Theor. Biol. 215 95
[19] Aguiar M A M, Sayama H, Baranger M and Bar-Yam Y 2003 Braz. J. Phys. 33 514
[20] de Carvalho K C and Tomé T 2004 Mod. Phys. Lett. B 18 873
[21] Nakagiri N and Tainaka K 2004 Ecol. Model. 174 103
[22] Szabó G 2005 J. Phys. A 38 6689
[23] Stauffer D, Kunwar A and Chowdhury D 2005 Physica A 352 202
[24] de Carvalho K C and Tomé T 2006 Int. Mod. Phys. C 17 1647
[25] Mobilia M, Georgiev I T and Tauber U C 2006 Phys. Rev. E 73 040903
[26] Arashiro E and Tomé T 2007 J. Phys. A 40 887
[27] Liggett T M 1985 Interacting Particle Systems (New York: Springer)
[28] Marro J and Dickman R 1999 Nonequilibrium Phase Transitions (Cambridge: Cambridge
University Press)
[29] Tomé T and de Oliveira M J 2001 Dinâmica Estocástica e Irreversibilidade (São Paulo: Editora
da Universidade de São Paulo)
[30] Tomé T 1994 Physica A 212 99
[31] Tomé T, Arashiro E, Drugowich de Feĺıcio J R and de Oliveira M J, 2003 Braz. J. Phys. 33 458
[32] Dickman R 1986 Phys. Rev. A 34 4246
[33] Tomé T and Drugowich de Feĺıcio J R 1996 Phys. Rev. E 53 3976
[34] Levins R 1969 Bull. Entomol. Soc. Am. 15 237
	Introduction
	Model
	Probabilistic cellular automaton
	Predator-prey probabilistic cellular automaton
	Time evolution equations for state functions
	Mean-field approximation
	One and two site approximations
	Patch model
	Quasi-spatial model
	Oscillatory behavior
	Discussion and conclusion

  We analyze a probabilistic cellular automaton describing the dynamics of
coexistence of a predator-prey system. The individuals of each species are
localized over the sites of a lattice and the local stochastic updating rules
are inspired on the processes of the Lotka-Volterra model. Two levels of
mean-field approximations are set up. The simple approximation is equivalent to
an extended patch model, a simple metapopulation model with patches colonized
by prey, patches colonized by predators and empty patches. This approximation
is capable of describing the limited available space for species occupancy. The
pair approximation is moreover able to describe two types of coexistence of
prey and predators: one where population densities are constant in time and
another displaying self-sustained time-oscillations of the population
densities. The oscillations are associated with limit cycles and arise through
a Hopf bifurcation. They are stable against changes in the initial conditions
and, in this sense, they differ from the Lotka-Volterra cycles which depend on
initial conditions. In this respect, the present model is biologically more
realistic than the Lotka-Volterra model.

Introduction
The simplest model exhibiting time-oscillations in a two-component system is the model
proposed independently by Lotka [1, 2, 3] and by Volterra [4]. In this model the
individuals of two species are dispersed over an assumed homogeneous space. It is
implicitly assumed in this approach that any individual can interact with any other one
with equal intensity implying that their positions are not taken into account. The time
evolution of the densities of the two species in the Lotka-Volterra model is given by a set
of two ordinary differential equations [5, 6, 7, 8] and is set up in analogy with the laws
of mass-action. Depending on the level of description wanted, the approach based on
mass-action laws, contained on the Lotka-Volterra model, suffices. However, there are
situations in which the coexistence takes place in a spatially heterogeneous habitat such
that the population densities can be very low in some regions. In this case we need to
proceed beyond the mass-law equations and consider the space structure of the habitat.
In other words, it becomes necessary to analyze the coexistence by taking explicitly into
account spatial structured models.
In fact, the role of space in the description of population biology problems has been
recognized by several authors in the last years [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26]. In a very clear manner, Durrett and Levin [11] have pointed
out that the modelling of population dynamics systems which are spatially distributed
by interacting particle systems [11, 27, 28, 29] is the appropriate theoretical approach
that is able to give the more complete description of the problem. We include in this
approach probabilistic cellular automata (PCA) [29, 30, 31], which will concern us here.
We refer to interacting particle systems and PCA as stochastic lattice models. They are
both Markovian processes defined by discrete stochastic variables residing on the sites
of a lattice; the former being a continuous time process and the latter a discrete time
process.
In the present work we study the coexistence and the emergence of stable self-
sustained oscillations in a predator-prey system by considering a PCA previously studied
by numerical simulations [24, 26]. This PCA is defined by local rules, similar to the
ones of the contact process [27], that are capable of describing the interaction between
prey and predator. Here, we focus on the analysis of the PCA by means of dynamic
mean-field approximations [10, 28, 29, 30, 32, 33]. In this approach the equations for
the time evolution of correlations of various orders are truncated at a certain level and
high order correlations of sites are written in terms of small order correlations. The
simplest approximation is the one in which all correlations are written in terms of one-
site correlation, called simple approximation. In a more sophisticated approximation,
called pair approximation [10, 29], any correlation is written in terms of one-site and
two-site correlations.
The simple mean-field approximation is capable of predicting coexistence of
individuals in a stationary state where the densities of each species, and of empty
sites, are constant. However, it is not capable of predicting possible time oscillating
Stable oscillations of a predator-prey probabilistic cellular automaton 3
behavior of the population densities and we have proceed to the next order of mean-
field approximation. The simple approximation, on the other hand, can be placed in an
explicit correspondence with a patch model [7, 12, 34] where unoccupied patches can be
colonized by prey and patches occupied by prey can be colonized by predators that in
turn may become extinct. In this approximation the PCA can be seen as an extended
version of the Lotka-Volterra model which includes an extra logistic term related to the
empty sites.
The pair-mean field approximation is able to predict possible time oscillating
behavior of the population densities that are self-sustained and are attained thorough
Hopf bifurcations. This is in contrast with the Lotka-Volterra model which presents
no stable oscillations but exhibits instead infinite cycles that are associated to different
initial conditions. However, from the biological point of view, one does not expect that a
small variation in the initial densities of prey and predator result in different amplitudes
of oscillations. Within our approach, a PCA treated in the pair-approximation, the
oscillations are associated with limit cycles what mean to say that they are stable
against the changes in the initial conditions. According to our point of view, the pair-
approximation, in which the correlation between neighboring sites are treated exactly,
provides a basic description of the predator-prey spatial interactions. For this reason,
we will refer to the PCA in this approximation as a quasi-spatial-structured model.
2. Model
2.1. Probabilistic cellular automaton
We consider interacting particles living on the sites of a lattice and evolving in time
according to Markovian local rules. The lattice is the geometrical object that plays the
role of the spatial region occupied by particles, in a general case, or by individuals of each
species in the present case. The lattice sites are the possible locations for the individuals.
Each site can be either empty or occupied by one individual of different species and a
stochastic variable ηi is introduced to describe the state of each site at a given instant
of time. The state of the entire system is denoted by η = (η1, . . . , ηi, . . . , ηN) where
N is the total number of sites. The transition between the states is governed by the
interactions between neighbor sites in the lattice and by a synchronous dynamics.
The probability P (ℓ)(η) of configuration η at time step ℓ evolves according to the
Markov chain equation
P (ℓ+1)(η) =
W (η|η′)P (ℓ)(η′), (1)
where the summation is over all the microscopic configurations of the system, and
W (η|η′) is the conditional transition probability from state η′ at time ℓ to state η at
time ℓ + 1. This transition probability does not depend on time and contains all the
information about the dynamics of the system. Taking into account that all the sites
are simultaneously updated, which is the fundamental property of a PCA, the transition
Stable oscillations of a predator-prey probabilistic cellular automaton 4
Figure 1. Transitions of the predator-prey model. The three states are: prey or
herbivorous (H or 1), predator (P or 2) and empty or vegetable (V or 0). The allowed
transitions obey the cyclic order shown.
probability can be factorized and written in the form [29, 30]
W (η|η′) =
wi(ηi|η
′), (2)
where wi(ηi|η
′) is the conditional transition probability that site i takes the state ηi
given that the whole system is in state η′. Being a probability distribution, the quantity
wi(ηi|η
′) must satisfy the following properties: wi(ηi, η
′) ≥ 0 and
wi(ηi|η
′) = 1. (3)
The average of any state function F (η) is evaluated by
〈F (η)〉ℓ =
F (η)P (ℓ)(η). (4)
The time evolution equation for 〈F (η)〉 is obtained from definition (4) and equation (1).
For example, we can derive the equations for the time evolution of densities and two-site
correlations.
2.2. Predator-prey probabilistic cellular automaton
To model a predator-prey system by a PCA, the stochastic variable ηi associated to site
i will represent the occupancy of the site by one prey, or the occupancy by one predator
or the vacancy (a site devoid of any individual). The variable ηi is assumed to take
the value 0, 1, or 2, according to whether the site is empty (V), occupied by a prey
individual (H) or by a predator (P), respectively. That is,
0, empty (V),
1, prey (H),
2, predator (P),
which defines a three-state per site PCA.
The stochastic rules, embodied in the transition rate wi(ηi|η
′), are set up according
to the following assumptions. (a) The space is homogeneous, which means to say that
no region of the space will be privileged against the others, that is, in principle the
individuals have the same conditions of surveillance in any space region. (b) The space
is isotropic, which means to say that there is no preferential direction in this space for
any interaction. (c) The allowed transitions between states are only the ones that obey
Stable oscillations of a predator-prey probabilistic cellular automaton 5
the cyclic order shown in figure 1. Prey can only born in empty sites; prey can give place
to a predator, in a process where a prey individual dies and a predator is instantaneously
born; finally a predator can die leaving an empty site. The empty sites are places where
prey can proliferate and can be seen as the resource for prey surveillance. The death of
predators complete this cycle, reintegrating to the system the resources for prey.
The predator-prey PCA has three parameters: a, the probability of birth of prey, b,
the probability of birth of predator and death of prey, and c, the probability of predator
death. Two of the process are catalytic: the occupancy of a site by prey or by a predator
is conditioned, respectively, to the existence of prey or predator in the neighborhood of
the site. The third reaction, where predator dies, is spontaneous, that is, it occurs, with
probability c, independently of the neighbors of the site. We assume that
a+ b+ c = 1, (6)
with 0 ≤ a, b, c ≤ 1.
The transition probabilities of the predator-prey PCA are described in what follows:
(a) If a site i is empty, ηi = 0, and there is at least one prey in its first neighborhood
there is a favorable condition for the birth of a new prey. The probability of site i being
occupied in next time step by a prey is proportional to the parameter a and to the
number of prey that are in the first neighborhood of the empty site.
(b) If a site is occupied by a prey, ηi = 1, and there is at least one predator in its
first neighborhood then the site has a probability of being occupied by a new predator
in the next instant of time. In this process the prey dies instantaneously. The transition
probability is proportional to the parameter b and the number of predators in first
neighborhood of the site.
(c) If site i is occupied by a predator, ηi = 2, it dies with probability c.
The transition probabilities associated to the three processes above mentioned can
be summarized as follows:
wi(0|η) = cδ(ηi, 2) + [1− fi(η)]δ(ηi, 0), (7)
wi(1|η) = fi(η)δ(ηi, 0) + [1 − gi(η)]δ(ηi, 1), (8)
wi(2|η) = gi(η)δ(ηi, 1) + (1− c)δ(ηi, 2), (9)
where
fi(η) =
δ(ηk, 1), gi(η) =
δ(ηk, 2), (10)
and the summation is over the four nearest neighbors of site i in a regular square lattice.
The notation δ(x, y) stands for the Kronecker delta function. These stochastic local
rules, when inserted in equation (2), define the dynamics of the PCA for a predator-
prey system.
The present stochastic dynamics predicts the existence of states, called absorbing
states, in which the system becomes trapped. Once the system has entered such a state
Stable oscillations of a predator-prey probabilistic cellular automaton 6
it cannot escape from it anymore remaining there forever. There are two absorbing
states. One of them is the empty lattice. Since the predator death is spontaneous, a
configuration where just predators are present is not stationary. This situation happens
whenever the prey have been extinct. In this case the predator cannot reproduce
anymore and also get extinct, leaving the entire lattice with empty sites. The other
absorbing state is the lattice full of prey. This situation occurs if there are few predators
and they become extinct. The remaining prey will then reproduce without predation
filling up the whole lattice. The existence of absorbing stationary states is an evidence
of the irreversible character of the model or, in other words, of the lack of detailed
balance [29]. However, the most interesting states, the ones that we are concerned with
in the present study, are the active states characterized by the coexistence of prey and
predators.
2.3. Time evolution equations for state functions
We start by defining the densities, which are the one-site correlations, and the two-site
correlations. These quantities will be useful in our mean-field analysis to be developed
below. The density of prey, predator, and empty sites at time step ℓ are defined thought
the expressions
i (1) = 〈δ(ηi, 1)〉ℓ, (11)
i (2) = 〈δ(ηi, 2)〉ℓ, (12)
i (0) = 〈δ(ηi, 0)〉ℓ. (13)
The evolution equations for the above densities are obtained from their definitions as
state functions, as given by equation (4), and by using the evolution equation for P (ℓ)(η),
given by equation (1). The resulting equations can be formally written as
(ℓ+1)
i (1) = 〈wi(1|η)〉ℓ, (14)
(ℓ+1)
i (2) = 〈wi(2|η)〉ℓ, (15)
(ℓ+1)
i (0) = 〈wi(0|η)〉ℓ, (16)
where the transition probabilities for this model are given in equations (7), (8) and (9).
The correlation between a prey localized at site i and a predator localized at site j
at time step ℓ is defined by
ij (1, 2) = 〈δ(ηi, 1)δ(ηj, 2)〉ℓ. (17)
The other two-site correlations are defined similarly. The time evolution equation for
the correlation of two neighbor sites i and j, one being occupied by a prey and the other
by a predator, is given by
(ℓ+1)
ij (1, 2) = 〈wi(1|η)wj(2|η)〉ℓ. (18)
The other two-site evolution equations are given by similar formal expressions. We
can also derive equations for three-site correlations. Since we are interested here on
Stable oscillations of a predator-prey probabilistic cellular automaton 7
approximations in which only the one-site and two-site correlations should be treated
exactly, the above equations suffice.
We call the attention to the fact that equation (18) includes the product of two
transition probabilities. This is a consequence of the synchronous update of the PCA
which allows that both neighboring sites i and j have their states changed at same
time step. This situation does not occur when we consider a continuous time one-
site dynamics. Therefore, although local interaction in the present PCA and in the
continuous time model considered in reference [10] are the same, the predator-prey
system evolves according to different global dynamics which leads to different time
evolution equations for the densities and the correlations.
The exact evolution equations for the one-site correlations are
P ′j (1) =
Pji(01)−
Pji(12) + Pj(1), (19)
P ′j (2) =
Pji(12) + (1− c)Pj(2), (20)
where the summation in j is over the ζ nearest neighbors of site i. To simplify notation
we are using unprimed and primed quantities to refer to quantities taken at time ℓ and
ℓ+ 1, respectively.
The exact evolution equations for the correlations of two nearest neighbor sites j
and k are
P ′jk(01) =
n(6=j)
Pjkn(001)−
i(6=k)
Pijkn(1001)
+ (1−
Pjk(01)−
n(6=j)
Pjkn(012)
i(6=k)
Pijk(101)−
n(6=j)
Pijkn(1012)
)Pjk(21)−
n(6=j)
Pjkn(212)
n(6=j)
Pjkn(201), (21)
P ′jk(12) =
n(6=j)
Pjkn(012) +
i(6=k)
Pijkn(1012)
n(6=j)
Pjkn(112)−
i(6=k)
Pijkn(2112)
+ (1− c)
)Pjk(12)−
i(6=k)
Pijk(212)
(1− c)
i(6=k)
Pijk(102), (22)
Stable oscillations of a predator-prey probabilistic cellular automaton 8
P ′jk(02) =
n(6=j)
)Pjkn(012)−
i(6=k)
Pijkn(1012)
+ (1− c)
cPjk(22) + Pjk(02)−
i(6=k)
Pijk(102)
Pjk(21) +
n(6=j)
Pjkn(212)
 , (23)
where the summation in i is over the nearest neighbors of j and the summation in n is
over the nearest neighbors of k.
3. Mean-field approximation
3.1. One and two site approximations
The evolution equation for a density in any interacting particle system which evolves in
time according to local interaction rules always contains terms related to the correlations
between neighbor sites in a lattice. The evolution equations for the correlations of two
neighbor sites includes the correlation of clusters of three or more sites in the lattice and
so on. In this way we can have an infinite set of coupled equations for the correlations
which is equivalent to the evolution equation for the probability P (ℓ)(η), described in
equation (1) for the automaton. The scope of the dynamic mean-field approximation
consists in the truncation of this infinite set of coupled equations [30, 31, 32, 33].
The lowest order dynamic mean-field approximation is the one where the probability
of a given cluster is written as the product of the probabilities of each site. That is, all
the correlations between sites in the cluster are neglected. For example, let us consider
the cluster constituted by a center (C) site and its first neighboring sites to the north
(N), south (S), east (E) and west (W) as shown in figure 2.
Within the one-site approximation the probability P (N,E,W, S, C) corresponding
to the cluster shown in figure 2 is approximated by
P (N,E,W, S, C) = P (N)P (E)P (W )P (S)P (C), (24)
where P (X), X = N,E,W, S, C are the one-site probabilities corresponding to each
site. For some stochastic dynamics models this approximation is able to give qualitative
results that are in agreement with the expected results.
In order to get a better approximation we must include fluctuations. The
simplest mean-field approximation that includes correlations is the pair-mean field
approximation. This approximation is better explained by taking again, as an example,
the cluster constituted by a center site which and its four nearest neighbors, shown
above. Within the pair-approximation the conditional probability P (N,E,W, S |C) is
approximated by
P (N,E,W, S |C) = P (N |C)P (E, |C)P (W |C)P (S |C), (25)
Stable oscillations of a predator-prey probabilistic cellular automaton 9
Figure 2. A site (C) of the square lattice and its four nearest neighbor sites (N, E,
W, S).
that is, the conditional probability P (N,E,W, S |C) is written in terms of the product
of the conditional probabilities P (X|C), X = N,E,W, S. Now using the definition of
conditional probability we have
P (N,E,W, S, C)
P (C)
P (N,C)
P (C)
P (E,C)
P (C)
P (W,C)
P (C)
P (S, C)
P (C)
, (26)
P (N,E,W, S, C) =
P (N,C)P (E,C)P (W,C)P (S, C)
[P (C)]3
. (27)
We see that the resulting probability is written as a function of two-site correlations
P (X,C), and the one-site correlation P (C).
3.2. Patch model
The simple mean-field approximation of the predator-prey PCA describes exactly
the same properties of an extended Levins patch model [7, 34]. That is, the PCA
with local rules similar to the contact process becomes, in the simple mean-field
approximation, analogous to the Levins model for metapopulation with empty patches,
patches colonized by prey and patches colonized by predators.
In the one-site mean-field approximation we consider that the probability of any
cluster of sites can be written as the product of the probabilities of each site, as in
equation (24). Using this approach, and writing x = Pi(1), y = Pi(2), and z = Pi(0) it
can be seen that the set of equations can be reduced to the following two-dimensional
map [26]
x′ = x+ axz − bxy, (28)
which is an evolution equation for prey density x, and
y′ = y + bxy − cy, (29)
which is an evolution equation for predator density yℓ. Notice that
z = 1− x− y. (30)
The fixed point of this map are those that represent the stationary solutions x′ = x
and y′ = y, and they correspond to the three following solutions x1 = 0, y1 = 0, and
x2 = 1, y2 = 0, and x3 = a/b, y3 = (1 − c/b)/(1 + b/a). The first solution corresponds
to an absorbing states where both species have been extinct. The second corresponds
Stable oscillations of a predator-prey probabilistic cellular automaton 10
−0.5 0 0.5
Figure 3. Phase diagram of the patch model. The continuous line represents
the transition, c1(p), between the prey absorbing (A) state and the active species
coexistence (C) state. The dashed line separates the two asymptotic time behavior of
the active state.
to an absorbing state where predators have extinct. The third solution corresponds to
an active state where prey and predator coexist.
Due to the constraint (6), the parameters a, b and c are not all independent and
only two can be chosen as independent. For this reason it is convenient to introduce the
following parametrization [10]
− p, b =
+ p, (31)
and consider p and c as the independent variables. The parameter p is such that
−1/2 ≤ p ≤ 1/2 and 0 ≤ c ≤ 1 as before. This parametrization will useful in the
determination of the different phases displayed by the model.
A linear stability analysis reveals that solution the (x1, y1) is a hyperbolic saddle
point for any set of the parameters a, b and c and so it is always unstable. The empty
absorbing state will never be reached. A linear stability analysis also shows that the
solution (x2, y2) is a stable node in the following region of the phase diagram c > c1
where
c1(p) =
(1 + 2p). (32)
The active solution is stable in the region c < c1 and is attained in two ways: by an
asymptotic stable focus, where the successive interactions of the map show damped
oscillations; or trough an asymptotic stable node. In the phase diagram of figure 3 we
show the transition line between the prey absorbing state and the active state given by
c = c1.
In figure 4 it is shown the behavior of the densities against the parameter c, the
probability of predators death, for the special case p = 0.2. In terms of phase transitions
what happens is that in the phase diagram there is a transition line separating the
absorbing prey phase and the active phase which is characterized by constant and
nonzero densities of prey and predator.
Stable oscillations of a predator-prey probabilistic cellular automaton 11
0 0.1 0.2 0.3 0.4 0.5 0.6
predator
Figure 4. Densities of predator and prey as functions of the parameter c for p = 0.2,
for the patch model.
We may conclude that the mean-field approximation for the predator-prey
probabilistic cellular automaton with rules (7), (8), and (9) is capable to show, under a
robust set of control parameters, that prey and predators can coexist without extinction.
However the map defined by equations (28) and (29) is not able to describe self-sustained
oscillations of species population densities.
3.3. Quasi-spatial model
In order to find if oscillations in the species populations can be described within a mean-
field approach we consider a more sophisticated approximation, the pair-approximation,
where correlations of two neighbor sites are included in the time evolution equations for
the densities. This is the lowest order mean-field approximation which takes into account
the spatial localization of neighboring individuals.
In this analysis we will maintain the correlations of one site and the correlations
of two-sites in the equations. Correlations of three and four neighbor sites will be
approximated by means of equation (27). With these approximation the model is
described by the following set of five coupled equations
x′ = au− bv + x, (33)
y′ = bv + (1− c)y, (34)
u′ = αa[
] + [(1− βa)− αa
][u− αb
+ αac
+ c[(1− βb)v − αb
], (35)
v′ = αb[βa
] + αa(1− c)
+ αb[
] + (1− c)[(1− βb)v − αb
], (36)
Stable oscillations of a predator-prey probabilistic cellular automaton 12
−0.5 0 0.5
Figure 5. Phase diagram of the quasi-spatial model. The upper continuous
line represents the transition, c1(p), between the prey absorbing (A) state and the
nonoscillating coexistence (CNO) state. The lower continuous line represents the
transition, c2(p), between the nonoscillating coexistence and the oscillating (COS)
coexistence state. The dashed line separates the two asymptotic time behavior of
the nonoscillating coexistence state.
w′ = αb[(1− βa)
] + (1− c)[w − αa
+ c[βbv + αb
] + c(1− c)s, (37)
where α and β are numerical fractions defined by α = (ζ − 1)/ζ and β = 1/ζ where
ζ is the coordination number of the lattice. For the present case of a square lattice,
ζ = 4 so that α = 3/4 and β = 1/4. We are using the following notation: u = P (0, 1),
v = P (1, 2), and w = P (0, 2) and also r = P (1, 1), q = P (0, 0) and s = P (2, 2). The
last three correlations are not independent but are related to others by
r = x− u− v, (38)
q = z − u− w, (39)
s = y − v − w. (40)
We used the properties P (1, 0) = P (0, 1), P (1, 2) = P (2, 1) and P (2, 0) = P (0, 2), that
follows from the assumption that space is isotropic and homogeneous.
We have analyzed numerically the five-dimensional map, described by the set of
equations (33), (34), (35), (36) and (37), and we have obtained four types of solutions.
Two solutions are trivial and are given by x = y = u = v = w = 0 and x = 1,
y = u = v = w = 0. They correspond to the empty and prey absorbing states,
respectively. The empty absorbing state, where both species have been extinct is
an unstable solution and never occurs. However, the prey absorbing state is one of
the possible stable stationary solutions and is stable above the critical transition line
Stable oscillations of a predator-prey probabilistic cellular automaton 13
0 0.1 0.2 0.3 0.4 0.5 0.6
predator
Figure 6. Densities of predator and prey as functions of c for the quasi-spatial model,
for p = 0.2.
c = c1(p) shown in figure 5. Below this line it becomes unstable giving rise to the active
state.
The other solutions correspond to the active states where both prey and predator
coexist. These solutions are of two kinds: a stationary solution where there is a
coexistence of the two species with densities constant in time, which we call the
nonoscillating (NO) active state; and another solution where both population densities
oscillate in time. This solution corresponds to a self-sustained oscillation of the predator-
prey system and will be called the oscillating (O) active state. In the phase diagram of
figure 5 there is a line c = c2(p) that separates the NO and O active phases. Figure 6
shows the behavior of the densities as a function of c for p = 0.
3.4. Oscillatory behavior
In figure 7 we show an example of self-sustained oscillations of the densities of prey
and predators as functions of time. The oscillating solutions are attained from the
nonoscillating solutions by a Hopf bifurcation. The fixed point associated to this solution
is an unstable center which produces a stable limit cycle as trajectories in the phase-
space of the predator density versus prey density, as can be seen in figure 7. Notice
that the oscillations are not damped and have a well defined period which is the same
for the prey density and for the predator density, which implies that the oscillations are
coupled. A maximum of predators always follow a maximum of prey. This means that
the abundance of prey is a condition that favors the increase in the number of predators.
As the predator number increases the prey population decays. The evanescence of prey
is followed by a decrease in the predator number, giving conditions for the increase of
prey population until the cycle starts again.
A well defined oscillatory behavior is found for many biological population, the
most famous being the one related to the time oscillations of the population of lynx
and snowshoe hare in Canada for which data were collected for a long period of time
Stable oscillations of a predator-prey probabilistic cellular automaton 14
1000 1200 1400 1600 1800 2000
predator
0 0.05 0.1 0.15
Prey population density
Figure 7. (a) Densities of predator and prey as functions of time and (b) density of
predator versus density of prey, for the quasi-spatial model, for p = 0 and c = 0.016.
[7, 8]. If the hare population cycles are mainly governed by the lynx cycle then the
oscillations shown by the present model reproduces qualitatively some of the features of
this predator-prey dynamics.
Next we analyze the behavior of the frequency and amplitude of oscillations. Fixing
the parameter p and varying the parameter c, we verify that in all the oscillating region
the frequency of oscillation is proportional to parameter c,
ω ∼ c, (41)
as can be seen in figure 8. Low frequencies are associated to low values of c; what means
that, for small values of c, the greater the predator lifetime the greater will be period of
the oscillation. As to the amplitude A of the oscillations, we have verified, that fixing
the value of p and varying the parameter c, it increases as c decreases. Our results show
that,
A ∼ (c− c2)
1/2, (42)
as expected for a Hopf bifurcation and shown in figure 8. The transition line c = c2
from the oscillating phase to the nonoscillating phase can either be obtained by using the
criterion given by equation (42) or by analyzing the eigenvalues associated to the map
given by the set of equations (33), (34), (35), (36) and (37). This last criterion means
to find the points of phase diagram such that the real part of the dominant complex
eigenvalue equals 1.
4. Discussion and conclusion
The main result coming from the pair mean-field approximation applied to the predator-
prey PCA is that it is possible to describe coexistence and self-sustained time oscillations.
Moreover, these are stable oscillations. Given a set of parameters, just one limit cycle
is achieved, no matter what the initial conditions are. This property is essential in
describing a biological system since a small variation in the initial condition can not
Stable oscillations of a predator-prey probabilistic cellular automaton 15
0 0.005 0.01 0.015 0.02
0 0.005 0.01 0.015 0.02
Figure 8. (a) Frequency of oscillations ω versus the parameter c. The frequency
vanishes linearly as one approaches c = 0. (b) Amplitude A of oscillations versus c
near the Hopf bifurcation point c2 = 0.019. The quantity A
2 vanishes linearly when
c → c2 in accordance with a Hopf bifurcation.
modify the amplitude, frequency and mean value of the time oscillation densities of a
predator-prey system. Similar results were obtained from a continuous time version of
the present model [10]. Although the simple mean-field equations are essentially the
same in both versions this is not the case concerning the pair mean-field approximation.
The time evolutions of the pair correlations for the PCA, presented here, depend on
higher order correlations (up to fourth) when compared to the ones of the continuous
version (up to third).
The model studied here is a spatial structured model with individuals residing in
sites of a lattice and described by discrete dynamic variables. When we perform simple
mean-field approximation we neglect all the correlations of sites in the lattice. But we
take into account that there are limited resources for the surveillance of each species.
For example in the time evolution equation for the density of prey we have an explicit
term relative to reaction of birth of prey which is the product of the density of prey x by
the density of empty sites z = (1−x−y). This coincides with an extended patch model
approach for predator-prey systems. The presence of this term is what differs the simple
mean-field equations from the Lotka-Volterra equations. However, taking into account
the limitation of space and resources the simple mean-field equations are not sufficient
to get self-sustained oscillations although able to describe damped time oscillations of
population densities.
To get self-sustained time oscillations we had to proceed to the next level of
approximation in which a pair of nearest neighbor sites is treated exactly. This
approximation can be seen as representing a pair of nearest neighbor sites immersed
in a mean field produced by the rest of the lattice. The most important feature being
the fact that the two sites of this pair can be seen as localized in space. The set of
five equations which results from the pair approximation for the PCA is indeed able
to produce self-sustained oscillations of population densities. It presents an important
Stable oscillations of a predator-prey probabilistic cellular automaton 16
property that the Lotka-Volterra model lacks, namely, the oscillating solutions are stable
and are unique for a given set of the control parameters.
Acknowledgements
The authors have been supported by the Brazilian agency CNPq.
References
[1] Lotka A 1920 J. Am. Chem. Soc. 42 1595
[2] Lotka A 1920 Proc. Nat. Acad. of Sciences USA 6 410
[3] Lotka A 1924 Elements of Mathematical Biology (new York: Dover)
[4] Volterra V 1931 Leçons sur la Théorie Mathématique de la Lutte pour la Vie Paris: Gauthier-
Villars)
[5] Haken H 1976 Synergetics, An Introduction (Berlin: Springer)
[6] Renshaw E 1991 Modelling Biological Populations in Space and Time (Cambridge: Cambridge
University Press)
[7] Hastings A 1997 Population Biology: Concepts and Models (New York: Springer)
[8] Ricklefs R E and Miller G L 2000 Ecology (New York: Freeman)
[9] Tainaka K 1989 Phys. Rev. Lett. 63 2688
[10] Satulovsky J and Tomé T 1994 Phys. Rev. E 49 5073
[11] Durrett R and Levin S 1994 Theor. Popul. Biol. 46 363
[12] Hanski I and Gilpin M E (eds.) 1997 Metapopulation Biology: Ecology, Genetic and Evolution
(San Diego: Academic Press)
[13] Satulovsky J and Tomé T 1997 J. Math. Biol. 35 344
[14] Tilman D and Kareiva P 1997 Spatial Ecology: the Role of Space in Population Dynamics and
Interactions (Princeton: Princeton University Press)
[15] Fracheburg L and Krapvisky P 1998 J. Phys. A 31 L287
[16] Liu Y C, Durrett R and Milgroom M 2000 Ecol. Model. 127 291
[17] Antal T, Droz M, Lipowsky A and Odor G 2001 Phys. Rev. E 64 036118
[18] Ovaskanien O, Sato K, Bascompte J and Hanski I 2002 J. Theor. Biol. 215 95
[19] Aguiar M A M, Sayama H, Baranger M and Bar-Yam Y 2003 Braz. J. Phys. 33 514
[20] de Carvalho K C and Tomé T 2004 Mod. Phys. Lett. B 18 873
[21] Nakagiri N and Tainaka K 2004 Ecol. Model. 174 103
[22] Szabó G 2005 J. Phys. A 38 6689
[23] Stauffer D, Kunwar A and Chowdhury D 2005 Physica A 352 202
[24] de Carvalho K C and Tomé T 2006 Int. Mod. Phys. C 17 1647
[25] Mobilia M, Georgiev I T and Tauber U C 2006 Phys. Rev. E 73 040903
[26] Arashiro E and Tomé T 2007 J. Phys. A 40 887
[27] Liggett T M 1985 Interacting Particle Systems (New York: Springer)
[28] Marro J and Dickman R 1999 Nonequilibrium Phase Transitions (Cambridge: Cambridge
University Press)
[29] Tomé T and de Oliveira M J 2001 Dinâmica Estocástica e Irreversibilidade (São Paulo: Editora
da Universidade de São Paulo)
[30] Tomé T 1994 Physica A 212 99
[31] Tomé T, Arashiro E, Drugowich de Feĺıcio J R and de Oliveira M J, 2003 Braz. J. Phys. 33 458
[32] Dickman R 1986 Phys. Rev. A 34 4246
[33] Tomé T and Drugowich de Feĺıcio J R 1996 Phys. Rev. E 53 3976
[34] Levins R 1969 Bull. Entomol. Soc. Am. 15 237
	Introduction
	Model
	Probabilistic cellular automaton
	Predator-prey probabilistic cellular automaton
	Time evolution equations for state functions
	Mean-field approximation
	One and two site approximations
	Patch model
	Quasi-spatial model
	Oscillatory behavior
	Discussion and conclusion
