## Brazilian Journal of Physics

*On-line version* ISSN 1678-4448

### Braz. J. Phys. vol.37 no.2a São Paulo June 2007

#### http://dx.doi.org/10.1590/S0103-97332007000300020

**REGULAR ARTICLES**

**Anisotropic probabilistic cellular automaton for a predator-prey system**

**Kelly C. de Carvalho; Tânia Tomé**

Instituto de Física, Universidade de São Paulo Caixa Postal 66318 05315-970 São Paulo, SP, Brazil

**ABSTRACT**

We consider a probabilistic cellular automaton to analyze the stochastic dynamics of a predator-prey system. The local rules are Markovian and are based in the Lotka-Volterra model. The individuals of each species reside on the sites of a lattice and interact with an unsymmetrical neighborhood. We look for the effect of the space anisotropy in the characterization of the oscillations of the species population densities. Our study of the probabilistic cellular automaton is based on simple and pair mean-field approximations and explicitly takes into account spatial anisotropy.

** Keywords:** Cellular automata; Predator-prey systems; Spatial models; Spatio-temporal patterns; Mean-field approach

**I. INTRODUCTION **

In the last years a particular great effort has been done in order to understand the role of space given by a spatial structure and local interactions in the characterization of the dynamics of competing biological species systems [1-18]. In this context it has been studied irreversible stochastic lattice models [19-21] with the purpose of mimic predator-prey systems with Markovian local rules based in the Lotka-Volterra model [22,23]. One of the problems that has been object of study is the connection between the time oscillations of population densities and spatial pattern distribution of the individuals of each species.

Here we study the coexistence of a two-species system by considering a probabilistic cellular automaton (PCA) which is a modified version of the automaton devised in [10,14,17,18]. The model, to be called anisotropic predator-prey PCA possess local rules that are similar to the ones of the cellular automaton proposed in [24] and was introduced by us in order to explore the effect of spatial anisotropy in the temporal oscillations.

We report dynamic mean-field approximations which take into account the spatial dependence of densities and pair correlations of sites. In the next section we present the model. In Sec. 2 we show the equations for the time evolution of species densities and show the spatial dependence of these equations. In Sec. 3 and 4 we perform dynamic simple and pair mean-field approximations. Last section summarize the model, method and results.

**II. MODEL**

The physical space occupied by the species is represented by a regular lattice of *N* sites in which each site can be in one of three states. At each site of the lattice we attach a stochastic variable h_{i} that takes the values 0, 1 and 2 according whether the site *i* is empty, or occupied by a prey or occupied by a predator, respectively. The state of the system can be represented by the vector h = (h_{1},h_{2}, ...,h* _{N}*). The time evolution equation for (h), the probability of state h at time , is given by

where the sum is over the 3^{N} configurations of the system and *W*(h| h^{'}) is the transition probability from a state h^{'} to state h, given that at the previous time step the system was in state h^{'}. Since we are considering probabilistic cellular automata, all the sites are updated simultaneously. In this case we have

where *w _{i}*(h

_{i}| h

^{'}) is the conditional transition probability per site.

**A. The anisotropic predator-prey PCA**

The stochastic rules, embodied in the transition rate *w _{i}*(h

_{i}| h

^{'}), are set up in order that the allowed transitions between states are only the ones that obey the cyclic order shown in Fig. 1. Prey can only born in empty sites; prey can give place to a predator in a process where the prey dies and the predator is instantaneously born; finally 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 anisotropic 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 with 0 __<__ *a,b,c* __<__ 1.

The transition probabilities of the anisotropic predator-prey PCA are described in what follows:

(a) If a site *i* is empty, h* _{i}* = 0, and if at north or east there is at least one prey, then it can be occupied by a prey, h

*= 1, in the next time step, with a probability proportional to the parameter*

_{i}*a*and to the number of prey

*n*at north and east of site

_{a}*i*.

(b) If a site is occupied by a prey, h* _{i}* = 1, then the site has a probability of being occupied by a new predator, h

*= 2, in the next time step if there are prey at north or east. In this process the prey dies instantaneously. The transition probability is proportional to the parameter*

_{i}*b*and the number of predators at north and east of site

*i*.

(c) If site i is occupied by a predator, h* _{i}* = 2, it dies spontaneously with probability

*c*.

The anisotropic cellular automaton is a variation of the automaton introduced in [10,14,17,18]. Here each site of a regular square lattice interacts with its first neighbors only at two preferential directions. This anisotropic neighborhood consists of the northern and eastern neighbors of each site as shown in Fig. 2. The set of transition probabilities per site is given by

and

where

and the sum is over the neighbor sites localized at east and north of site *i* and correspond to the number of neighbors of site *i* occupied by prey individuals and predators individuals, respectively.

We have considered this probabilistic cellular automaton with the purpose of verifying the effect of anisotropy in the properties of the time oscillations of the predator-prey system. The rules considered are in some sense inspired in the north-east-center (NEC) cellular automaton [24], which also consider an unsymmetrical neighborhood of northern and eastern sites.

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 it cannot escape from it anymore remaining there forever. There are two absorbing states. One of them is the empty lattice. The other absorbing state is the one in which the lattice is 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 [21]. 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.

**B. Evolution equation for state functions**

The densities of prey, predator and empty sites are defined as

respectively. The lower index *i* is used to denote the site and the upper index stands for the time. The pair correlation of two neighbor sites *i* and *j* one being occupied by prey and the other by predator is defined as

In this definition the sites *i* and *j* are two neighboring horizontal sites, the site *j* being at the east of *i*. For two neighboring sites placed vertically we use the notation

where *k* denotes the neighbor of *i* to the north.

The evolution equation for the densities can be obtained by using the rules of the automaton and the evolution equation for the probability (h). Using equations (2), (3), (4) and (5), we obtain the time evolution equation for the density of prey at site *i*,

where we have used prime and unprimed to denote quantities at time +1 and , respectively. Again, denotes a pair correlation of a site i which is empty and its neighbor at the north j which is occupied by a prey and *P _{ik}*(0 1) denotes the pair correlation of the site

*i*which is is empty and its neighbor to the east

*k*which is occupied by a prey. Note that site

*i*is the site to be updated. The time evolution equation for the density of predators at site i is given by

The evolution equations for correlations of two neighbor sites are given by equations which involves clusters of two, three and four neighbors to the north and east and they are too cumbersome. We will write them in the pair approximation in the next subsection.

We will assume in the present analysis that the densities and the pair correlations are homogeneous so that the correlations of sites placed horizontally or vertically are the same. In other words the systems exhibits a specular symmetry along the southwest-northeast line (see Fig. 2). Therefore,

and

**III. SIMPLE MEAN-FIELD APPROXIMATION**

We will first analyze the equations (12) and (13) for the densities of prey and predators by means of a simple mean field approximation [21,26,27]. In this approximation we write the probability of a given cluster of sites as the product of the probabilities of each site. In our analysis we maintain the space dependence of each probability. This is necessary since the rules that define the model are not isotropic. We will use the notations

where the index *n* stands for a site at the *n* layer. A layer is defined as the sites belonging to a line perpendicular to the southwest-northeast axis, as shown in Fig. 2. Using this approximation and considering equations (12) and (13) we obtain

and

where *z _{n}* = 1-

*x*. Due to this relation there is no need for the equation for the density of empty sites.

_{n}-y_{n}The analysis of the above set of equations show that the stable solutions are of two types: a prey absorbing state, where the all lattice is full of prey; and an active solution with prey and predators densities constant in time and space. So, we have just homogeneous and nonoscillating solutions when we treat the anisotropic PCA by means of simple mean-field approximation. At this level of mean-field approximation the space anisotropy does not play any role in determining the kind of coexistence of species.

**IV. PAIR-MEAN FIELD APPROXIMATION**

Now we consider the evolution equations for one-site correlations and two site (pair) correlations. The evolution equations for the pair correlations depend on clusters of three and four sites. Following the approach of the pair approximation [25-28] we write the clusters of three and four sites as products of pair correlations and one-site correlations. However, we must be careful in doing this approximation because we are maintaining the spatial anisotropy dependence of the one and two site correlations. For example, the correlation of a site at layer *n* in state 0, a site at layer *n*+1 in state 0 and a site at layer *n*+2 in state 1, is written in the pair approximation as

There are more complex clusters that appear in the evolution equation for pair correlations. For example, the following cluster of four sites is approximated by

Observe that *P _{nn}*

_{+1}(10) and

*P*

_{nn}_{+1}(01) are two independent variables.

Now we use the notation

Using the pair approximation we derive the following equations for pair correlations,

The equations for the densities are

and

Therefore, we get a set of eight equations that we have iterated to find the solutions. We considered periodic boundary conditions and an initial condition where half of the lattice with a set of densities and the other half with another set of densities.

For *a = b* and great values of *c* the system attains a stationary state which is the absorbing prey state. Decreasing *c* there is a transition to an active state, coexistence of species, and we found very interesting solutions which are characterized by travelling waves of the densities of each species. In Fig. 3 we show, for a particular set of the parameters *a = b* and *c* = 0.05, the prey and predators oscillations as a function of space and for a fixed instant of time. We can see that the behavior of the space oscillations are very complex and very different from a sinusoidal wave.

We remark that our assumption concerning the spatial dependence of densities and pair correlations underlying our mean-field approach could not easily be conceived a priori. It was set up by considering the type of local dynamics with unsymmetrical rules and from the results of Monte Carlo simulations for the present PCA [10,29]. To further clarify this point we show in Fig. 4 a snapshot generated by simulation of the present PCA on a square lattice where each site can be occupied by prey, predators or can be empty. Periodic boundary conditions were used and the system evolved in time according to the rules defined by Eqs. (3), (4) and (5). A synchronous update was used. We see that the distribution of individuals of each species is not homogeneous but exhibits a pattern composed by layers of prey and predators. These layers are displayed in a manner perpendicular to the southwest-northeast direction, the direction along which the oscillations occur as can be seen in Fig. 4. We have found that the spatial oscillations do not occur independently of the time oscillations, since the spatial layers of prey and predators are not static but are like fronts moving along the southwest-northeast direction. The spatial pattern oscillations are intimately associated to local time oscillations of the species [10,29]. These features were incorporated in the present mean-field analysis of the model.

**V. CONCLUSIONS**

We have considered a predator-prey probabilistic cellular automaton with anisotropic local stochastic rules. We studied this model by means of dynamic mean-field approximation at two orders: simple mean-field approximation and pair mean-field approximation. Due to the anisotropy these approximations only can be performed if we consider the spatial dependence of probabilities. The simple mean-field approximation for this automaton just provides the prey absorbing state and active spatial homogeneous solutions which are also constant in time. The pair mean-field approximation gives much more rich results and show that the active states characterized by time oscillations of species densities are inhomogeneous in space. Therefore, we have spatiotemporal patterns of coexistence. This is in accordance with previous Monte Carlo simulations [10] where we have found local time oscillations connected to inhomogeneous spatial distributions of species individuals.

[1] K. Tainaka, Phys. Rev. Lett. **63**, 2688 (1989). [ Links ]

[2] R. Durrett and S. Levin, Theor. Popul. Biol. **46**, 363 (1994). [ Links ]

[3] A. Hastings, *Population biology: concepts and models*, Springer, New York, 1997. [ Links ]

[4] J. Satulovsky and T. Tomé, Phys. Rev. E bf 49, 5073 (1994). [ Links ]

[5] D. Tilman and P. Kareiva, *Spatial ecology: the role of space in population dynamics and interactions*, Priceton University Press, Princeton, 1997. [ Links ]

[6] J. Satulovsky and T. Tomé, J. Math. Biol. **35**, 344 (1997). [ Links ]

[7] R. Durrett and S. Levin, J. Theor. Biol. **205**, 201 (2000). [ Links ]

[8] T. Antal and M. Droz, Phys. Rev. E **63**, 056119 (2001). [ Links ]

[9] M. A. M. de Aguiar, H. Sayama, M. Baranger, and Y. Bar-Yam, Braz. J. Phys. **33**, 514 (2003). [ Links ]

[10] K. C. de Carvalho and T. Tomé, Mod. Phys. Lett. B **18**, 873 (2004). [ Links ]

[11] A. F. Rozenfeld and E. V. Albano, Phys. Lett. A **332**, 361 (2004). [ Links ]

[12] D. Chowdhury and D. Stauffer, J. Biosc. **30**, 277 (2005). [ Links ]

[13] G. Szabó, J. Phys. A **38**, 6689 (2005). [ Links ]

[14] K. C. de Carvalho and T. Tomé, Int. J. Mod. Phys. C **17**, 1647 (2006). [ Links ]

[15] M. Mobilia, I. T. Giorgiev, and U. C. Täuber, Phys. Rev. E **73**, 040903 (2006). [ Links ]

[16] S. Morita and K. Tainaka, Popul. Ecol. **48**, 99 (2006). [ Links ]

[17] E. Arashiro and T. Tomé, J. Phys. A **40**, 887 (2007) [ Links ]

[18] T. Tomé and K. C. de Carvalho, Stable oscillations in a predator-prey system: a mean-field approach, arXiv:0704.0512v1 cond-mat stat-mech.

[19] T. M. Liggett, *Interacting Particle Systems*, Springer, New York, 1985. [ Links ]

[20] J. Marro and R. Dickman, *Nonequilibrium Phase Transitions*, Cambridge University Press, Cambridge, 1999. [ Links ]

[21] T. Tomé and M. J. de Oliveira, *Dinâmica Estocástica e Irreversibilidade*, Editora da Universidade de São Paulo, São Paulo, 2001. [ Links ]

[22] A. Lotka, J. Am. Chem. Soc. **42**, 1595 (1920); [ Links ]Proc. Nat. Acad. of Sciences USA **6**, 410 (1920); * [ Links ]Elements of Mathematical Biology*, Dover, New York (1924). [ Links ]

[23] V. Volterra, *Leçons sur la Theorie Mathématique de la Lutte pour la Vie*, Gauthier-Villars, Paris, 1931. [ Links ]

[24] C. H. Bennett and G. Grinstein, Phys. Rev. Lett. **55**, 657 (1985). [ Links ]

[25] R. Dickman, Phys. Rev. A bf 34, 4246 (1986). [ Links ]

[26] T. Tomé, Physica A **212**, 99 (1994). [ Links ]

[27] T. Tomé and J. R. Drugowich de Felício, Phys. Rev. E **53**, 3976 (1996). [ Links ]

[28] T. Tomé, E. Arashiro, J. R. Drugowich de Felício, and M. J. de Oliveira, Braz. J. Phys. **33**, 458 (2003). [ Links ]

[29] T. Tomé and K. C. de Carvalho, in preparation.

Received on 13 April, 2007