Competitive Brownian and Levy walkers
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r Competitive Brownian and L´evy walkers
E. Heinsalu,
1, 2
E. Hern´andez-Garc´ıa, and C. L´opez IFISC, Instituto de F´ısica Interdisciplinar y Sistemas Complejos (CSIC-UIB),Campus Universitat de les Illes Balears, E-07122 Palma de Mallorca, Spain National Institute of Chemical Physics and Biophysics, R¨avala 10, Tallinn 15042, Estonia (Dated: November 18, 2011)Population dynamics of individuals undergoing birth and death and diffusing by short or longranged twodimensional spatial excursions (Gaussian jumps or L´evy flights) is studied. Competitiveinteractions are considered in a global case, in which birth and death rates are influenced by allindividuals in the system, and in a nonlocal but finite-range case in which interaction affects indi-viduals in a neighborhood (we also address the noninteracting case). In the global case one singleor few-cluster configurations are achieved with the spatial distribution of the bugs tied to the typeof diffusion. In the L´evy case long tails appear for some properties characterizing the shape anddynamics of clusters. Under non-local finite-range interactions periodic patterns appear with peri-odicity set by the interaction range. This length acts as a cut-off limiting the influence of the longL´evy jumps, so that spatial configurations under the two types of diffusion become more similar. Bydividing initially everyone into different families and following their descent it is possible to showthat mixing of families and their competition is greatly influenced by the spatial dynamics.
PACS numbers: 05.40.-a, 05.40.Fb, 87.18.Hf, 87.23.Cc
I. INTRODUCTION
Birth and death are the most relevant processes in de-termining the dynamics of biological populations whichin the context of statistical physics can be modeled us-ing interacting particle models where particle number ischanging in time. As it is understood by now, birth anddeath processes are also responsible for clustering mech-anisms in systems where random-walking individuals un-dergo reproduction and death. As a result, aggregationof organisms can occur even in simple models where birthand death processes are combined with spatial diffusion.In fact, in the most simple
Brownian bug model , whereparticles reproduce and die with the same probability andundergo Brownian motion [1–3], clustering of particleswas observed. In this model the clustering is producedsimply by the reproductive correlations (the offspring isborn at the same location of the parent) and by the irre-versibility of the death process. Birth and death modelsof moving individuals are the pertinent framework to cap-ture properties of biological systems such as planktonicpopulations [2], or patterns in amoebae [4] and bacteria[5].Taking into account another central ingredient that ispresent in ecological systems, namely, the competitionwith other individuals in the neighborhood for resources,the formation of periodic spatial structures was observedin Refs. [6–8]. In these nonlocally interacting Brownianbug models it was assumed that the reproduction proba-bility depends on the number of other organisms in theneighborhood. In Ref. [9] nonlocally interacting L´evybugs , i.e., reproducing and dying organisms that undergoL´evy flights, were studied. This type of motion is relevantto model cell migration [10], biological searching strate-gies [11, 12], bacteria dynamics [13], or pattern formationof mussels [14]. In Ref. [9] it was shown that the forma- tion of a periodic pattern is robust with respect to thetype of spatial motion that the particles perform. Theperiodic arrangement of clusters in these nonlocally in-teracting bug models is a consequence of the competitiveinteraction and has a spatial scale determined by the in-teraction range [6]. However, a deeper analysis of the dif-ferences and similarities between the Brownian and L´evycases is still missing. In particular, as shown in [7, 15],this analysis can be very conveniently performed by con-sidering the limit of the interaction distance reaching thesystem size (global interaction), since a unique clusterappears which helps to understand and characterize thecluster properties, and the fluctuations of the populationsize.In the present paper we report on differences betweenthe systems of Brownian and L´evy bugs, in the situationsof global and non-local interactions, as well as in the non-interacting case. In addition, results on the dependenceof population on diffusion, and mixing of families of par-ticles are presented for the finite-range interaction case.The paper is organized as follows: in Sec. II we describethe models to be analyzed. In Sec. III the noninteractingbug systems are studied. The infinite competition rangewhere each particle is competing with all the others is an-alyzed in Sec. IV. Finally, the nonlocally interacting (i.e.with a finite interaction range) models are investigatedin Sec. V.
II. MODEL AND NUMERICAL ALGORITHM
We consider a system consisting initially of N point-like particles, which could be thought as being biologicalorganisms or bugs, placed randomly in a two-dimensional L × L square domain with periodic boundary condi-tions. Except when explicitly stated, we take L = 1,so that lengths are measured in units of system size.Bugs diffuse, reproduce at rate r ib , and die at rate r id ; i = 1 , . . . , N , and N ≡ N ( t ) is the number of bugs inthe system at time t . The numerical algorithm used toevolve the system follows the one suggested in Ref. [16].The following sequence of steps is repeated until the finalsimulation time is reached:We first compute the random time τ after which thenext demographic event (birth or death) will occur. Forthis we need to determine the total birth and death rates, B tot = N X i =1 r ib , D tot = N X i =1 r id , (1)and compute also the total rate R tot = B tot + D tot = N X i =1 ( r ib + r id ) . (2)For the random times τ we choose an exponential prob-ability density with the complementary cumulative dis-tribution p ( τ ) = exp( − τ / ˜ τ ) (3)so that values of τ could be generated from τ = − ˜ τ ln( ξ ),where ξ is a uniform random number on (0 ,
1) [17]. Thecharacteristic time or time-scale parameter ˜ τ = h τ i isdetermined by the total rate:˜ τ = R − . (4)After the random time τ , an individual i , chosenamong all the N ( t ) bugs, either reproduces or disappears.With probability B tot /R tot the event is reproduction andwith probability D tot /R tot it is death. The probability ofchoosing a particular individual i is weighted proportion-ally to its contribution to the corresponding total rate.In the case of reproduction, the new bug is located atthe same position ( x i , y i ) as the parent individual i . Fi-nally, all the bugs perform a jump of random length ℓ ina random direction characterized by an angle uniformlydistributed on (0 , π ) ( ℓ and the direction of the jump areindependent for each particle). The new present time is t + τ , bugs are relabeled with indices i = 1 , , ..., N ( t + τ ),and the process is repeated.When bugs undergo normal diffusion (Brownian bugs),a Gaussian jump-length probability density function isused, ϕ ( ℓ ) = 2˜ ℓ √ π exp (cid:18) − ℓ ℓ (cid:19) , l ≥ h ℓ i = ˜ ℓ ; ˜ ℓ is the space-scale pa-rameter. Since we draw the angle specifying the directionof the jump from the interval (0 , π ), we restrict ℓ in Eq.(5) to have positive sign. The random jump length ℓ canbe computed from ℓ = ˜ ℓ ξ G , where ξ G is sampled from the standard Gaussian distribution with average 0 and stan-dard deviation 1, and neglecting the sign. Note that therandom walk defined in this way is not exactly the sameas the one in which the walker performs jumps extractedfrom a two-dimensional Gaussian distribution, but it alsoleads to normal diffusion and allows a more direct com-parison with the L´evy case. The corresponding diffusioncoefficient can be estimated as κ = h ℓ i / (2 h τ i ) = ˜ ℓ / (2˜ τ ) . (6)As we choose to fix the value of κ , and the demographicrates, then the space-scale parameter is determined by˜ ℓ = √ κ ˜ τ = p κ/R tot . (7)In order to simulate the system where the bugs un-dergo superdiffusive L´evy flights (L´evy bugs) one canuse a symmetric L´evy-type probability density functionfor the jump size ( ℓ ≥ ϕ µ ( ℓ ) ≈ ˜ ℓ µ | ℓ | − µ − , ℓ → ∞ ( ℓ ≫ ˜ ℓ ) (8)with the L´evy index 0 < µ <
2. For all L´evy-typeprobability density functions with µ < h ℓ i = ∞ , leading to the occurrence ofextremely long jumps, and typical trajectories are self-similar, showing at all scales clusters of shorter jumpsinterspersed with long excursions. For 0 < a < µ < h ℓ a i are finite. For the L´evy index inthe range 1 < µ < h ℓ i is finite. The com-plementary cumulative distribution corresponding to (8)behaves as P µ ( ℓ ) ≈ µ − ( ℓ/ ˜ ℓ ) − µ , ℓ → ∞ . (9)As a simple form of complementary cumulative distri-bution function which behaves asymptotically as (9), weuse P µ ( ℓ ) = (1 + b /µ ℓ/ ˜ ℓ ) − µ , (10)with b = [Γ(1 − µ/ µ/ / Γ( µ ), and ℓ ≥
0. As before,the direction of the jump is assigned by drawing an ran-dom angle on (0 , π ). The particular expression for b ischosen for consistency with previous work [9]. It gives tothe tail of the jump distribution the same prefactor as forthe L´evy-stable distribution [20], but any other positivevalue of b should lead to the same results reported here.One can generate a random step-length ℓ by inverting(10): ℓ = ˜ ℓ ( ξ − /µ − b /µ . (11)with ξ being a uniform random variable on the unit in-terval. Now the diffusion coefficient (6) is infinite, butone can define a generalized diffusion coefficient in termsof the scales ˜ ℓ and ˜ τ as [18, 19] κ µ = ˜ ℓ µ / (2˜ τ ) . (12)Therefore, in the case of the L´evy flights, when fixing thevalue of κ µ , the space-scale parameter is:˜ ℓ = (2 κ µ ˜ τ ) /µ = (2 κ µ /R tot ) /µ . (13)As we consider the bugs to be point-like, the spatialdynamics does not include any interaction between them.The interaction is instead taken into account throughreproduction and death rates, which we assume to beaffected by competitive interactions.If the birth and death rates of a bug are influenced bythe number of other bugs within a certain radius R , onetalks about a nonlocal interaction of finite range. In thepresent paper we assume that the birth and death ratesof the i -th individual depend linearly on the number ofneighbors in the interaction range [6], r ib = max (cid:0) , r b − αN iR (cid:1) , (14) r id = max (cid:0) , r d + βN iR (cid:1) . (15)Here N iR is the number of bugs which are at a distancesmaller than R from bug i , the parameters r b and r d are the zero-neighbor birth and death rates, and the pa-rameters α and β determine how r ib and r id depend onthe neighborhood. For positive values of α and β , themore neighbors an individual has within the radius R ,the smaller is the probability for reproduction and thelarger is the probability that the bug does not survive,which could arise from competition for resources. Thefunction max() enforces the positivity of the rates. Sincewe take R < L/ R ≪ L ), the periodic boundaryconditions are straightforwardly implemented and bugsare never counted twice.If the birth and death rates of a bug are instead influ-enced by all the other individuals in the system, i.e., r ib ≡ r b = max { , r b − α [ N ( t ) − } , (16) r id ≡ r d = max { , r d + β [ N ( t ) − } , (17)then one talks about global interaction . This is formallyequivalent to Eqs. (14) and (15) with R sufficiently largefor the interaction domain to include the whole system,but taking care of counting each bug only once, so that N iR = N ( t ) − r ib ≡ r b = r b , r id ≡ r d = r d , (18)which is equivalent to α = β = 0 in Eqs. (14), (15), bugsare noninteracting .In the following we discuss the Brownian and L´evy bugsystems when individuals do not influence each other andwhen inter-particle interaction occurs, either global ornonlocal. Although we formally maintain the parameter β in Eqs. (15) and (17), in our numerical examples werestrict to β = 0. III. SIMPLE BUG MODELS WITH NOINTERACTIONA. Noninteracting Brownian bugs
The simple Brownian bug model with no interaction,i.e., when the birth and death rates of the individualsare given by Eq. (18), has been studied and discussed invarious works [1, 2]. The ensemble average of the totalpopulation size follows h N ( t ) i = N exp [∆( t − t )] , (19)independently of the diffusivity of the bugs; it only de-pends on the difference ∆ = r b − r d . If the birth rate islarger than the death rate, ∆ >
0, the total populationgenerally explodes exponentially, though there is a finiteprobability for extinction that depends on the initial sizeof the population and decreases with increasing ∆. Ifthe death rate is larger than the birth rate, ∆ < h N ( t ) i = N and the average lifetime is infinite. However, in singlerealizations the fluctuations in the number of individualsare huge leading to fast extinction in some runs. In fact,there exists a typical lifetime proportional to N , definedas the time for which the fluctuations become as large asthe mean value, beyond which the population is extinctwith probability close to 1 [1].As a surprising effect, in the systems where the nonin-teracting Brownian bugs undergo death and reproductionwith equal probabilities, spatial clustering of the bugswas observed in single realizations [1–3]. A typical timeevolution of such a system is illustrated in Fig. 1a. Wenote that in all figures presenting the spatial configura-tions of the bugs, we have divided the individuals accord-ing to their initial position into nine groups characterizedby different colors as in Fig. 1 at time t = 0; if reproduc-tion takes place, the newborn bug assumes the same coloras the parent. From Fig. 1a one can see that many smallclusters form some time after starting from a uniforminitial distribution. The occurrence of the clustering isrelated to the fact that in the case of reproduction thenew bug is located at the same position as the parent.Due to the fluctuations and irreversibility of death thenumber of clusters decreases in time, until there will bea single cluster consisting of individuals coming from asingle ancestor. There are constant, spontaneous, short-lived break-offs from the main cluster, which are alwayslocated near it. The center of mass of such a cluster un-dergoes a motion similar to that of a single bug [1] andits linear width fluctuates with a typical value propor-tional to √ N [1]. Furthermore, the larger the diffusioncoefficient κ , the wider is the cluster (notice that whensimulating the system numerically, if the diffusivity be-comes so large that the jump lengths become comparableto the system size, one needs to take a larger simulation t = 0 (a) y t = 100 t = 1000 t = 10000 (b) x y x x x FIG. 1: (Color) Simple bug models with no interaction; spatial configuration of bugs at different times t : (a) Brownian bugs with κ = 10 − and (b) L´evy bugs with κ µ = 10 − and µ = 1. Reproduction and death occur with equal probability, r b = r d = 0 . N = 1000. Bugs are colored with the color their ancestors bear in the panel at t = 0. box). Finally, due to the fluctuations also the last clusterdisappears. B. Noninteracting L´evy bugs
In the case of noninteracting L´evy bugs, the numberof individuals still follows Eq. (19), independently of theL´evy index µ , and also the cluster formation observed inthe case of Brownian bugs takes place. Now, however, asbugs can perform long jumps, there are also small clus-ters continuously appearing and disappearing far fromthe main clusters (Fig. 1b). The smaller the value of µ the more anomalous the system, i.e., the larger is theprobability for long jumps and therefore there are moreflash-clusters. When the number of clusters has alreadydecreased to one, due to the long jumps and fluctuationof the number of individuals, new clusters that are placedfar from the central cluster can appear in the system alsofor some time and often the disappearance of the maincluster takes place whereas another new central clusterappears somewhere else. As a result the center of massundergoes anomalous diffusion as single bugs do. Thevalue of the L´evy index µ influences also the linear sizeof the main clusters: the smaller is µ , the more compactare the clusters, although also more particle jumps tolong distances occur. The influence of the value of κ µ is similar as in the case of Brownian bugs, i.e., a largervalue of the anomalous diffusion coefficient results in alarger linear size of the clusters. IV. GLOBAL INTERACTIONA. Formation of a cluster
Let us now investigate the behavior of the Brownianand L´evy bug systems in the case of global interaction,i.e., birth and death rates of the individuals are given byEqs. (16), (17). The time evolutions of the globally in-teracting Brownian and L´evy bug systems are illustratedby Fig. 2a and 2b, respectively. In both systems we startfrom N = 500 bugs uniformly distributed in the sim-ulation area and choose for the parameters characteriz-ing death and birth rates the following values: r b = 1, r d = 0 . α = 0 . β = 0. As in the noninteractingcase, the final state of the dynamics is complete extinc-tion, since there is always a nonvanishing probability fora fluctuation strong enough to produce that. However, ifthe number of bugs in the system is large this happensat very long times [21]. Then, there is a long-lived qua-sistable state for which the average number of individuals h N ( t ) i can be estimated from the condition that deathand birth are equally probable, r ib = r id . From there, h N ( t ) i = ∆ α + β + 1 , (20)where ∆ = r b − r d . We have restricted to parametervalues so that the max functions in Eqs. (16-17) do notoperate. Since we have chosen N > h N ( t ) i = 46 in Fig.2, death is more probable at small times and the numberof bugs decreases rapidly. Approximately at time t = 30the number of individuals has reached the value at which t = 0 (a) y t = 20 t = 50 t = 100 t = 500 (b) x y x x x x FIG. 2: (Color) Globally interacting bug models; spatial configuration of bugs at different times t : (a) Brownian bugs with κ = 10 − and (b) L´evy bugs with κ µ = 10 − and µ = 1. The parameters in the reproduction and death rates are: r b = 1, r d = 0 . α = 0 . β = 0. Bugs are colored with the color their ancestors bear in the panel at t = 0. death and birth become in average equally probable andafter this time particle number fluctuates around thatvalue; parameters of the birth and the death rates canbe chosen so that these fluctuations are weak. At thistime small clusters start to form due to the reproduc-tive pair correlations. As in the case of noninteractingbugs, fluctuations and irreversibility of death makes thenumber of clusters to decrease in time, although now thefluctuations of the particle content of the different clus-ters are correlated to keep the total number close to thevalue given by Eq. (20) and the process is faster. Finallya single cluster consisting of bugs coming from the sameancestor remains (as stated before, it will also disappearat very long times due to finite-size fluctuation effects)though there are also now spontaneous short-lived break-offs from the central cluster as in the case of noninteract-ing bugs. The center of mass of such a cluster is movingin space and its linear size is a fluctuating quantity. Theclustering of the globally interacting bugs was quantita-tively discussed in Ref. [7] for the one-dimensional Brow-nian bug system. B. Fluctuations of the number of bugs
As indicated by Eq. (20), for given values of α and β , the average number of individuals in the system withglobal interaction depends solely on the difference ∆ = r b − r d . It is independent of the concrete values of r b and r d , as well as of the value of κ or κ µ and µ ; in fact, itdoes not even depend on whether the system consists ofBrownian or L´evy bugs. Nevertheless, fluctuations of thenumber of bugs do indeed depend on the values of r b and r d , even if the difference ∆ , and thus the average num-ber of bugs, has the same value. To illustrate this, let us calculate from the simulations time series the probabilitydistribution of the number of individuals in the globallyinteracting Brownian and L´evy bug systems. As can beseen from Fig. 3, for a given value of ∆ , larger valuesof r b and r d lead to larger fluctuations. This is a sim-ple consequence of the Poisson character of the birth anddeath events for which fluctuations in each of the instan-taneous rates are proportional to the mean rates. Whenthe distributions are narrow, they are close to Gaussian.For larger rates particle number distribution gets broaderimplying that there is an enhanced probability that par-ticle number becomes zero at some moment, after whichbugs become extinct (remember that what is in fact plotin Fig. 3 is the numerical particle number distribution inthe long-lived metastable state before extinction). Forthe present case with ∆ = 0 . α = 0 . β = 0,rate values above the ones shown in Fig. 3 (i.e. r b > r d > .
1) lead to observable extinction after some tenthsof thousands of steps. An ecological implication of thiscould be the following: one can think of two colonies oforganisms of the same type, having both the same equi-librium size determined for example by the size of theliving area. Now if in one of the systems the populationhas no enemies and the natural death rate is low, butin the other the death rate is higher due to the presenceof a predator, then the latter system will more proba-bly go to extinction sooner due to the presence of largerfluctuations.
C. The average cluster shape, cluster width, andcenter of mass motion
Let us keep in the following α = 0 . β = 0 and r b = 1, r d = 0 . p N (a)(b) (c) FIG. 3: (Color online) Probability distribution of the numberof bugs in globally interacting bug systems. The results arenumerically obtained from the time series of the particle num-ber in the very long-lived state before the fluctuations lead thesystem to the extinction. For all the curves α = 0 . β = 0,and the rate difference is ∆ = r b − r d = 0 .
9, but the rates r b , r d assume different values: (a) r b = 1, r d = 0 .
1; (b) r b = 1 . r d = 0 .
6; (c) r b = 2, r d = 1 .
1. The overlappingcurves correspond to Brownian and L´evy bug systems; thedistributions do not depend on the type of diffusion nor onthe values of κ , κ µ or µ in this globally interacting case. Fig. 2 and in Fig. 3 for curve (a)] and study the behaviorof the cluster formed in the case of a system with globalinteraction defined by Eqs. (16)-(17). As mentioned, evenafter the transition from an uniform distribution of bugsto a single cluster (and before eventual extinction at largetimes), at some moment the system can consist actuallyof more than one cluster. In such cases we define that allthe individuals in the system belong to the same cluster,even though in the L´evy case the distance between thebugs (subclusters) can be rather large. In order to avoidthe boundary effects as much as possible, in Figs. 4-7 thelinear size of the simulation area was taken as L = 1000and to have enough statistics simulations were run until t = 5 × .Let us start by analyzing the average shape of the clus-ter. The average cluster, ρ ( x, y ), is obtained setting ateach time the origin in the center of mass of the cluster(distances under the periodic boundary conditions arecomputed under a minimum distance convention) andaveraging over a long time (after the transition from uni-form distribution to one single cluster but before long-time extinction). A one-dimensional cut of it (say across x for y = 0, i.e., ρ ( x ) ≡ ρ ( x, y = 0)) is shown in Figs. 4and 5. For the case of Brownian bugs, the tail of the aver-age cluster is approximately exponential. A pair distribu-tion function, which should be related but not identical tothe average cluster discussed here, was analytically calcu-lated in Ref. [16] for a globally interacting Brownian bugmodel of our type but in which total extinction was for-bidden. This quantity also displayed a fast decaying tail. In the case of L´evy bugs the tail of ρ ( x ) follows insteada power law, ρ ( x ) ∼ x − (2+ µ ) , see Fig. 4b, arising fromthe long jumps. Note that, in the present case of circularsymmetry, the relation ρ ( x, y ) dx dy = R ( r )(2 π ) − dr dθ of ρ ( x, y ) with the radial density of the average cluster, R ( r ), where r and θ are the polar coordinates centeredat the cluster center, implies ρ ( x ) = R ( r = | x | )(2 π | x | ) − ,so that the asymptotic behavior of the radial density is R ( r ) ∼ r − (1+ µ ) . This is the same asymptotic behavioras the individual radial jumps in (8) and it is also theasymptotic tail of the probability of displacement fromthe original position of nonreproducing bugs moving byL´evy flights [19]. We note also that, for κ = κ µ , thecentral part of ρ ( x ) is narrower and higher in the L´evythan in the Brownian bug system, and it is narrower andhigher the smaller the value of µ (see Fig. 4a). This isa somehow counterintuitive effect of the L´evy motion onclusters, already commented in the noninteracting case:increasing anomalous diffusion (smaller µ ) induces largerjumps and longer tails, but at small scales it acts as mak-ing the cluster more compact.The influence of the diffusivity is similar in both sys-tems: the larger is the value of κ or κ µ the more spreadis the average cluster (see Fig. 5). For the Brownian one-dimensional case it was shown in Ref. [7] that clusterwidth is essentially the distance associated to the Brow-nian walk during the lifetime of a bug and their descen-dants. Thus, the width increases as κ / . In the L´evycase, defining the distance associated to the walk is moresubtle, since higher moments of displacements diverge.But the behavior of lower ones and dimensional analysisindicate that typical displacements during a lifetime scaleas κ /µµ , and then this should determine the width of ρ ( x )or R ( r ) (i.e., the spatial dependence should occur onlythrough the combinations [ xκ − /µµ ] or [ rκ − /µµ ]). Impos-ing additionally that the total number of bugs in the aver-age cluster in this global interaction case does not dependon particle motion or distribution, and it is thus indepen-dent on the value of κ µ we have R ( r ) = κ − /µµ F ( rκ − /µµ ),or ρ ( x ) = 1 κ /µµ G xκ /µµ ! , (21)with G ( u ) = F ( u ) /u . The analogous scaling form for theaverage cluster in the Gaussian diffusion case is ρ ( x ) = 1 κ G (cid:16) xκ / (cid:17) . (22)The insets in Fig. 5 show the validity of these scalingforms.In addition to the average cluster shape, which givesinformation on the cluster width, it is also interesting tostudy the fluctuations of the cluster width in time. Wecharacterize cluster width at each instant of time by thestandard deviation of the bug positions with respect tothe center of mass of the cluster at that time. Then, aprobability density π ( σ ) is constructed from the values -5 -3 -1 -0.5 -0.25 0 0.25 0.5 ρ x (a) LB, µ = 0.511.5 BB -1 ρ x (b) µ = 0.511.5 FIG. 4: (Color online) (a) ρ ( x ), the cross-section of the two-dimensional particle density of the average cluster in semi-log scale; comparison between the Brownian and L´evy bugsystems; κ = 10 − , κ µ = 10 − , r b = 1, r d = 0 . α = 0 . β = 0. (b) The tails of ρ ( x ) in log-log scale in the case ofthe L´evy bug systems for different values of µ . Solid linescorrespond to fitting curves ∝ x − (2+ µ ) . of σ at different times. Figure 6 shows that in the caseof globally interacting Brownian bugs the distribution of σ is short-tailed. In the case of globally interacting L´evybugs, in contrast, the distribution for σ is characterizedby tails with a power law decay with exponent − (1 + µ ).This means that in the latter case the cluster width canundergo arbitrarily large fluctuations in time. We notethat the tails in π ( σ ) decay with the same exponent as theradial density R ( r ) of the average cluster, thus suggest-ing that the tails of the average cluster are produced bythe large fluctuations in the width of the instantaneousclusters (which in fact include splitting events).The individual motion of bugs drives the behavior ofthe center of mass of the system. Figure 7 depicts theprobability density p (∆ CM ) of the jump lengths ∆ CM performed by the center of mass each time the bug mo-tion step is executed in the globally interacting Brownianand L´evy bug systems. For Brownian bugs it is short- FIG. 5: (Color online) ρ ( x ), the cross-section of the two-dimensional particle density of the average cluster in semi-log scale for different values of diffusivity: (a) Brownian bugsand (b) L´evy bugs with µ = 1. Other parameters are as inFigs. 2 and 4. The insets check the correctness of the scalingforms (22) (with G ( u ) = κρ and u = x/κ / ) and (21) (with G ( u ) = κ /µµ ρ and u = x/κ /µ ). tailed. In fact, from the arguments in Ref. [7], the cen-ter of mass motion of the cluster for globally interactingBrownian bugs is characterized by a Brownian processwith the same diffusion coefficient as the individual bugs.In the case of L´evy bugs the probability density of thejump lengths of the center of mass is described by a dis-tribution with a power-law tail with exponent − (1 + µ ),i.e., the center of mass of the cluster formed in the caseof globally interacting L´evy bugs undergoes jumps thatfollow asymptotically the same law as the single bugs,Eq. (8), and as the radial tails of the average cluster.This reflects the fact commented previously that, due tothe long jumps of the L´evy bugs, additional clusters farfrom the main one appear from time to time, stronglydisplacing the center of mass of the system. Due to thefluctuations it is even possible that the cluster that usedto be the main cluster disappears and a new main clus-ter forms somewhere else. As a result the center of massmotion undergoes the same type of superdiffusion as theindividual bugs of the system.Extending the arguments for the Brownian bugs [7] -7 -5 -3 -1 -2 -1 π ( σ ) σ LB, µ = 0.511.5 BB FIG. 6: (Color online) Probability density π ( σ ) of the stan-dard deviation σ of the bug positions with respect to thecenter of mass of the cluster in the Brownian and L´evy bugsystems; κ = 10 − , κ µ = 10 − . The distribution is obtainedaveraging over a long time. The curves corresponding to theL´evy bug systems are well fitted by ∝ σ − (1+ µ ) (not shown).Other parameters are as in Figs. 2, 4, and 5. -8 -6 -4 -2 -5 -4 -3 -2 -1 p ∆ CM LB , µ = 0.511.5 BB FIG. 7: (Color online) Probability density p (∆ CM ) of thejump lengths of the center of mass in the Brownian and L´evybug systems. Same parameter values as in Fig. 6. The tailsof the curves corresponding to L´evy bugs are well fitted by ∝ ∆ − (1+ µ ) CM (not shown). (which were themselves adapted from the ones in [1]) tothe L´evy case one can heuristically show that the dis-tributions of σ and ∆ CM are related. To this aim onemakes the approximation that the number of bugs in thesystem is constant, say N , instead of being constant onaverage. The center of mass receives a positive contribu-tion from the new bugs appearing (at location ~x i ) due tothe reproduction between diffusion steps (say at time t i ),a negative contribution from the bugs disappearing dur-ing that time (say from position ~x j at time t j ), and the direct contribution from the L´evy jumps ~ℓ k of all bugspresent at the diffusion step: ~ ∆ CM = 1 N X i ∈ B ~x i ( t i ) − N X j ∈ D ~x j ( t j ) + 1 N N X k =1 ~ℓ k . (23) B and D denote the sets of bugs that have been born ordead, respectively, between diffusion steps. The two firstterms can combined in a single one ~S ≈ N − P np =1 ~σ p by considering that the two sets have approximately thesame number of individuals, n . ~σ p = ~x i − ~x j is the dis-placement between a pair of these bugs, one just bornand the other just disappeared, sampled inside the samecluster. Then the modulus of each σ p should be of the or-der of the cluster width σ , which fluctuates in time withprobability tails ruled by an exponent − (1 + µ ). Thiscontribution in Eq. (23) gives the motion of the center ofmass due to the birth and death processes. The contri-bution from the individual particle jumps is in the lastterm in (23), which is a sum of L´evy jumps of parameter µ so that the tails of the probability density are char-acterized by a decay with the same exponent − (1 + µ ).These heuristic arguments imply that the modulus ∆ CM will also be distributed with long tails characterized byan exponent − (1 + µ ), as observed. V. NONLOCAL INTERACTIONA. Formation of a periodic pattern
In Refs. [6, 8, 9] on the nonlocally interacting Brow-nian and L´evy bugs it was assumed that the birth anddeath rates of the i -th individual are given by Eqs. (14),(15). In the case of Brownian bugs, for small enoughdiffusion coefficient and large enough ∆ , the occurrenceof a periodic pattern consisting of clusters that are ar-ranged in a hexagonal lattice was observed (see Fig. 8a-c) [6, 8]. For large values of the diffusion coefficient suchperiodic pattern is replaced by a more homogeneous dis-tribution of bugs (Fig. 8d). In the case of L´evy bugs,since the diffusion coefficient (6) is infinite, one couldexpect that the spatial distribution will not reveal a pe-riodic pattern; however, as shown in Ref. [9], for properparameters periodic cluster arrangements do indeed oc-cur (see Fig. 9). The reason for the divergence of thediffusion coefficient in the L´evy case is in the statisticalweight of large jumps. These large jumps have some in-fluence on the characteristics of the pattern formed, butthe relevant structure is ruled mainly by the interactionsbetween individuals. In the L´evy bug system however, atvariance with the Brownian case, even at small values of κ µ there are many solitary bugs appearing for short timeperiods in the space between the periodically arrangedclusters due to the large jumps [9], c.f. Figs. 8a and 9a.However, the periodicity of the pattern is of the orderof R (the interaction range) in both systems, being only (a) y (b)(c) x y (d) x FIG. 8: (Color) Interacting Brownian bug model with R =0 . r b = 1, r d = 0 . α = 0 . β = 0. Spatial configura-tion of bugs at time 45000 in systems with different diffusioncoefficients: (a) κ = 10 − , (b) κ = 2 × − , (c) κ = 4 × − ,(d) κ = 10 − . The initial configuration of bugs is the sameas in Figs. 1 and 2 at time t = 0. slightly influenced by κ or κ µ and µ , as demonstrated inRefs. [6, 9] through a mean-field theory calculation.In Ref. [9] also the two-dimensional particle density ofthe average cluster, obtained by setting the origin at thecenter of mass of each cluster forming the periodic pat-tern and averaging over all the clusters in the simulationarea and over time, was studied. In both, Brownian andL´evy bug systems the central part of the average cluster,where most of the individuals are concentrated, was wellfitted by a Gaussian function, but the way the particledensity decreases when moving away from the center ofmass of the cluster is rather different. In the Browniancase a Gaussian decay provides a good approximation,whereas in the L´evy case it is close to exponential. Thecomparison with the systems with global interaction, dis-cussed in Sec. IV C, reveals therefore that the interactionrange R turns the exponential decay into Gaussian andthe power law decay into exponential.For a given value of diffusion coefficient, there existsa critical value of ∆ below which the system gets ex-tinct, independently of α [6]. Above this critical value,for every α the increase of ∆ results in the increase ofthe average number of bugs, but the pattern formation isnot much influenced. The latter is, however, true solelyif ∆ increases through the increase of r b and the deathrate is low. Namely, as in the case of global interac-tion discussed in Sec. IV B, an increase of the death rate,though accompanied by a compensating increase of birthrate, leads to larger fluctuations in the particle number.In numerical simulations we have observed that the largerare the fluctuations in the number of bugs, the more dif- (a) y (b)(c) x y (d) x FIG. 9: (Color) Interacting L´evy bug model with R = 0 . r b = 1, r d = 0 . α = 0 . β = 0 (same parametersas in Fig. 8 for Brownian bugs). The spatial configurationof bugs at time 45000 in systems with different generalizeddiffusion coefficients and anomalous exponent: (a) κ µ = 10 − , µ = 1; (b) κ µ = 10 − , µ = 1; (c) κ µ = 10 − , µ = 1 .
5; (d) κ µ = 5 × − , µ = 1 .
5. The initial configuration of bugs isthe same as in Figs. 1 and 2 at time t = 0. ficult is the formation of the periodic pattern, and finallythe individuals do not arrange in the periodic pattern butin random clusters (see also Ref. [16]). This effect mayin fact make difficult to observe the periodic clusteringphenomenon in real competitive biological systems.In the following we keep for the parameters in the birthand death rate the same values as in the case of globalinteraction, i.e., r b = 1, r d = 0 . α = 0 . β = 0.For these parameter values the number of bugs fluctu-ates only weakly around the mean value. Differently fromthe case of global interaction, now the average number ofbugs in the system is influenced not only by the birth anddeath rates, but also by the diffusion, i.e., in the case ofBrownian bugs by κ and in the case of L´evy bugs by κ µ and µ , see Fig. 10. The smaller is κ , κ µ or µ , the largerthe particle number. At the same time Figs. 8 and 9 in-dicate that by decreasing κ or κ µ the linear width of theclusters becomes smaller, the particle density in the clus-ters higher, and the density between the clusters lower(c.f. Sec. IV C and see also Ref. [9]). Somehow counter-intuitively, the effect of decreasing µ seems to have thesame effects, as commented above for the noninteractingand global cases. Furthermore, the value of κ or κ µ and µ seems to weakly influence the number of clusters in thesystem: In Figs. and 8 and 9 smaller values lead to largernumber of clusters. This observation is not explained bythe linear instability analysis of [6, 9].0 (a) × -5 × -5 × -5 × -5 -4 〈 N 〉 κ (b) 〈 N 〉 µ κ µ = 10 -5 × -5 -4 FIG. 10: (Color online) a): Average number of bugs versusdiffusion coefficient in the system with Gaussian jumps. b):Average number of bugs versus anomalous exponent µ in thesystem with L´evy jumps for various values of the anomalousdiffusion coefficient. Other parameters as in Figs. 8 and 9. B. Mixing of different families
It is interesting to study the evolution of the systemalso regarding the disappearance or survival of the dif-ferent groups, by dividing initially the bugs into differ-ent families and following their descent. In the case ofnonlocally interacting Brownian bugs, a very low diffu-sion coefficient leads to the situation in which after clus-ter formation the inter-cluster travel is very rare becausethe individuals are not capable to make the jumps fromone cluster to another one, and it is also very unlikelyto arrive to the next cluster doing a multistep randomwalk because death is very probable between the clus-ters. Therefore, in the case of very low diffusion differentfamilies would remain inside their initial clusters. If oneassumes that initially each individual represents a dif-ferent family, then only inter-cluster competition occursand the final number of families is equal to the numberof clusters. If instead initially individuals are assigned to families according to large areas of initial positions(larger than typical cluster size as done in Figs. 1 and2 at time t = 0), there is no family competition inter-nal to the clusters, most families survive and the clusterscoming from different families occupy approximately theterritory of the ancestors even after a long time, as can beseen from Fig. 8a. In that case, the travel of a cluster toa new territory away from the other clusters of the samefamily can take place due to the diffusion of the clusteras a whole during the clusters arrangement into the pe-riodic pattern. For larger values of κ the inter-clustertravel is possible which leads to the conquering of newterritories, i.e., bugs can be found in a region where theirancestors were not from, Fig. 8b. The effect is larger forlarger κ and leads to the disappearance of some families,as can be seen from Fig. 8c. Finally, for increased dif-fusion, intra-cluster competition will force all survivingbugs to be from a single family (in fact, from a single an-cestor); which family (ancestor) wins is a random event.The process is faster for larger diffusion. Increasing thediffusivity further even the periodic pattern disappears,Fig. 8d.Figure 9 illustrates the family mixing for nonlocally in-teracting L´evy bugs. In this case the inter-cluster travel-ing is supported by the long jumps. Differently from thecase of Brownian bugs, now the individuals can reachnot only the next neighboring cluster but also clustersfar away. Consequently, a cluster originally consisting ofbugs coming from one ancestor can after some time turninto a cluster consisting of bugs coming from differentfamilies placed initially far away. Thus, intra-cluster bugcompetition becomes soon competition between families,and even if the diffusivity of the bugs is very low, at theend the L´evy bug system would consist of individualscoming from one or just a few ancestors. As in the caseof Brownian bugs, the process of the disappearance offamilies is faster the greater is the generalized diffusioncoefficient.Besides the diffusion of a cluster as a whole during theformation of the periodic pattern and the conquering ofnew territories through the migration to and survival inanother cluster, the mixing of clusters from different fam-ilies can take place also due to the appearance of a newcluster if in the periodic pattern there is a dislocation.In the case of Brownian bugs the new cluster is formedthrough the splitting of an old cluster. In the case ofL´evy bugs, instead, the new cluster can appear also farfrom the original territory. VI. CONCLUSIONS AND OUTLOOK
We have presented some detailed properties of interact-ing particle systems in which the individuals are Brown-ian or L´evy random walkers which interact in a compet-itive manner. We have seen strong differences betweenthe globally and the finite-range nonlocally interactingsystems. In the systems with global interaction the spa-1tial distribution of the bugs becomes tied to the type ofdiffusion, Brownian or L´evy. Typical configurations con-sist of a single or a few clusters for both types of motion.For the L´evy bug systems long tails appear in the meancluster shape and in probability distributions of clusterwidth and of jumps of the center of mass. For Brownianbug systems these quantities appear to be much shorterranged. This is qualitatively also the situation in the non-interacting case, although then the effects of the particle-number fluctuations are much stronger. Under non-localfinite-range interactions the situation is rather different.First, single cluster configurations are generally replacedby periodic patterns with periodicity set by the interac-tion range R . Motion of individual clusters is severelyrestricted by the presence of the neighboring clusters. Inaddition, the natural spatial cut-off introduced by the in-teraction range R seems to limit the influence of the longL´evy jumps, so that measures of spatial cluster shapedo not generally exhibit power laws, making spatial con-figurations under both types of diffusion more similar.Mixing of families and their competition is neverthelessgreatly influenced by the type of motion. This suggeststhat it would be interesting to consider the influence ofdifferent types of diffusion into competitive genetic mix-ing processes [22]. Obtaining analytic understanding in this type of inter-acting systems is difficult, but at least the nature of theinstability leading to pattern formation and its relevantspatial scale have been clarified by using partial integro-differential equation descriptions of the mean field type[6, 9], which are useful in broader contexts [23–25]. How-ever, from previous work in the Gaussian case [7, 8, 26], itis known that quantities such as cluster width and struc-ture or transition thresholds strongly depend on particle-number fluctuations. Thus, obtaining additional resultsfrom differential equation approaches would need the in-clusion of effective multiplicative noise terms [5] or focuson statistical quantities such as pair correlation functions[2, 16, 27]. Acknowledgments
This work has been supported by the targeted financ-ing project SF0690030s09, Estonian Science Foundationthrough grant no. 7466, by the Balearic Government(E.H.), and by Spanish MICINN and FEDER throughproject FISICOS (FIS2007-60327). [1] Y.-C. Zhang, M. Serva, and M. Polikarpov, J. Stat. Phys. , 849 (1990).[2] W. R. Young, A. J. Roberts, and G. Stuhne, Nature ,328 (2001).[3] J. Felsenstein, The American Naturalist , 359 (1975).[4] B. Houchmandzadeh, Phys. Rev. Lett. , 078103(2008).[5] F. Ramos, C. L´opez, E. Hern´andez-Garc´ıa, and M. A.Mu˜noz, Phys. Rev. E , 021102 (2008).[6] E. Hern´andez-Garc´ıa and C. L´opez, Phys. Rev. E ,016216 (2004).[7] E. Hern´andez-Garcia and C. L´opez, J. Phys.: Condens.Matter , S4263 (2005).[8] C. L´opez and E. Hern´andez-Garcia, Physica D , 223(2004).[9] E. Heinsalu, E. Hern´andez-Garcia, and C. L´opez, Euro-phys. Lett. , 40011 (2010), Erratum: , 69902 (2011).[10] P. Dieterich, R. Klages, R. Preuss, and A. Schwab, Proc.Natl. Acad. Sci. USA , 459 (2008).[11] D. Sims, E. Southall, N. Humphries, G. Hays, C. Brad-shaw, J. Pitchford, A. James, M. Ahmed, A. Brierley,M. Hindell, et al., Nature , 1098 (2008).[12] F. Bartumeus, M. G. E. D. Luz, G. M. Viswanathan, andJ. Catal´an, Ecology , 3078 (2005).[13] M. Levandowsky, B. S. White, and F. L. Schuster, ActaProtozool. , 237 (1997).[14] M. de Jager, F. J. Weissing, P. Herman, B. A. Nolet, andJ. van de Koppel, Science , 1551 (2011).[15] E. Brigatti, V. Schwammle, and M. A. Neto, Phys. Rev. E , 021914 (2008).[16] D. A. Birch and W. R. Young, Theoretical PopulationBiology , 26 (2006).[17] W. H. Press, B. P. Flannery, S. A. Teulolsky, and W. Vet-terling, Numerical Recipes in C: The Art of ScientificComputing (Cambridge University Press, Cambridge, 2edition, 1992).[18] R. Klages, G. Radons, and I. M. Sokolov,
AnomalousTransport: Foundations and Applications (Wiley-VCH,2008).[19] R. Metzler and J. Klafter, Phys. Rep. , 1 (2000).[20] J. P. Nolan,
Stable Distributions - Modelsfor Heavy Tailed Data (Birkhauser, Boston,2012), to appear. Chapter 1 online athttp://academic2.american.edu/ ∼ jpnolan.[21] C. R. Doering, K. V. Sargsyan, and L. M. Sander, Mul-tiscale Modeling & Simulation , 283 (2005).[22] H. Sayama, M. A. M. de Aguiar, Y. Bar-Yam, andM. Baranger, Phys. Rev. E , 051919 (2002).[23] M. A. Fuentes, M. N. Kuperman, and V. M. Kenkre,Phys. Rev. Lett. , 158104 (2003).[24] M. G. Clerc, D. Escaff, and V. M. Kenkre, Phys. Rev. E , 056217 (2005).[25] Y. E. Maruvka and N. M. Shnerb, Phys. Rev. E ,011903 (2006).[26] E. Hern´andez-Garc´ıa and C. L´opez, Physica A , 95(2005).[27] B. Houchmandzadeh, Phys. Rev. E80