Boundary Effects on Population Dynamics in Stochastic Lattice Lotka-Volterra Models
BBoundary Effects on Population Dynamics inStochastic Lattice Lotka–Volterra Models
Bassel Heiba, Sheng Chen, and Uwe C. T¨auber
Department of Physics (MC 0435) and Center for Soft Matter and Biological Physics,Robeson Hall, 850 West Campus Drive, Virginia Tech, Blacksburg, Virginia 24061, USA
Abstract
We investigate spatially inhomogeneous versions of the stochastic Lotka–Volter-ra model for predator-prey competition and coexistence by means of MonteCarlo simulations on a two-dimensional lattice with periodic boundary condi-tions. To study boundary effects for this paradigmatic population dynamicssystem, we employ a simulation domain split into two patches: Upon settingthe predation rates at two distinct values, one half of the system resides in an ab-sorbing state where only the prey survives, while the other half attains a stablecoexistence state wherein both species remain active. At the domain bound-ary, we observe a marked enhancement of the predator population density. Thepredator correlation length displays a minimum at the boundary, before reach-ing its asymptotic constant value deep in the active region. The frequency of thepopulation oscillations appears only very weakly affected by the existence of twodistinct domains, in contrast to their attenuation rate, which assumes its largestvalue there. We also observe that boundary effects become less prominent as thesystem is successively divided into subdomains in a checkerboard pattern, withtwo different reaction rates assigned to neighboring patches. When the domainsize becomes reduced to the scale of the correlation length, the mean populationdensities attain values that are very similar to those in a disordered system withrandomly assigned reaction rates drawn from a bimodal distribution.
Email address: [email protected], [email protected], [email protected] (Bassel Heiba, ShengChen, and Uwe C. T¨auber)
Preprint submitted to Physica A November 13, 2018 a r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug eywords: stochastic population dynamics, Lotka–Volterra predator-preycompetition, spatial inhomogeneities, boundary effects, quenched disorder
1. Introduction
Due to its wide range of applications and relative simplicity, variants of theLotka–Volterra predator-prey competition model represent paradigmatic sys-tems to study the emergence of biodiversity in ecology, noise-induced patternformation in population dynamics and (bio-)chemical reactions, and phase tran-sitions in far-from-equilibrium systems. In the classical deterministic Lotka–Volterra model [1, 2], two coupled mean-field rate equations describe the popu-lation dynamics of a two-species predator-prey system, whose solutions displayperiodic non-linear oscillations fully determined by the system’s initial state. Yetthe original mean-field Lotka–Volterra rate equations do not incorporate demo-graphic fluctuations and internal noise induced by the stochastic reproductionand predation reactions in coupled ecosystems encountered in nature. In a se-ries of analytical [3]–[6] and numerical simulation studies [7]–[18], the populationdynamics of several stochastic spatially extended lattice Lotka–Volterra modelvariants was found to substantially differ from the mean-field rate equationpredictions due to stochasticity and the emergence of strong spatio-temporalcorrelations: Both predator and prey populations oscillate erratically, and donot return to their initial densities; the oscillations are moreover damped andasymptotically reach a quasi-stationary state with both population densities fi-nite and constant on one- or two-dimensional square lattices [7], whereas damp-ing appears absent or is very weak in three dimensions [9]. Very similar dynam-ical properties are observed in other two-dimensional model variants, includ-ing a predator-prey system with added prey food supply and cover [12], andimplementations on a triangular lattice [14]. In a non-spatial setting, the per-sistent non-linear oscillations can be understood through resonantly amplifieddemographic fluctuations [19]. Local carrying capacity restrictions, representinglimited resources in nature, can be implemented in lattice simulations by con-2training the number of particles on each site [8, 10, 11, 15, 16, 17]. These localoccupation number restrictions cause the emergence of a predator extinctionthreshold and an absorbing phase, wherein the predator species ultimately dis-appears while the prey proliferate through the entire system. Upon tuning thereaction rates, one thus encounters a continuous active-to-absorbing state non-equilibrium phase transition whose universal features turn out to be governedby the directed percolation universality class [4, 5, 10, 11, 12, 15, 16, 18].Biologically more relevant models should include spatial rate variability toaccount for environmental disorder. The population dynamics in a patch sur-rounded by a hostile foe [20, 21, 22] is well represented by Fisher’s model [23],which includes diffusive spreading as well as a reaction term capturing interac-tions between individuals and with the environment. For the stochastic Lotka–Volterra model, the influence of environmental rate variability on the populationdensities, transient oscillations, spatial correlations, and invasion fronts was in-vestigated by assigning random reaction rates to different lattice sites [24, 25].Spatial variability in the predation rate results in more localized activity patches,a remarkable increase in the asymptotic population densities, and acceleratedfront propagation. These studies assumed full environmental disorder, as therewas no correlation at all between the reaction rates on neighboring sites.In a more realistic setting, the system should consist of several domainswith the environment fairly uniform within each patch, but differing markedlybetween the domains, e.g., representing different topographies or vegetationstates. In our simulations, we split the system into several patches and assigndifferent reaction rates to neighboring regions. By tuning the rate parameters,we can force some domains to be in an absorbing state, where the predatorsgo extinct, or alternatively in an active state for which both species coexist atnon-zero densities. One would expect the influence of the boundary between theactive and absorbing regions to only extend over a distance on the scale of thecharacteristic correlation length in the system. In this work, we study the localpopulation densities, correlation length, as well as the local oscillation frequencyand attenuation, as functions of the distance from the domain boundary. As we3uccessively divide the system further in a checkerboard pattern so that eachpatch decreases in size, the population dynamics features quantitatively tendtowards those of a randomly disordered model with reaction rates assigned tothe lattice sites from a bimodal distribution.
2. Model Description and Background
The deterministic classical Lotka–Volterra model [1, 2] is a set of two couplednon-linear dynamical rate equations that on a mean-field approximation levelcapture the following kinetic reactions of two species, respectively identified withpredators A and prey B : A µ → ∅ , A + B λ → A + A , B σ → B + B . (1)In these stochastic processes, µ corresponds to the spontaneous predator deathrate, while σ denotes the prey reproduction rate. Finally, λ is the predation ratewhich describes the non-linear reaction through which the predator and preyspecies interact with each other. The simplified Lotka–Volterra model thus as-sumes that the prey population grows exponentially in the absence of predators,but becomes diminished with growing predator population. In the presence ofthe prey, the predator population will increase with the prey population, butis subject to exponential decay once all prey are gone. The configuration withvanishing predator number represents an absorbing state for this system, sincethere exists no stochastic reaction process that would allow recovery from it. Forcompleteness, we mention that the total population extinction state of courserepresents another absorbing state. We also remark that one could add inde-pendent predator reproduction A → A + A (with rate σ (cid:48) ) and prey death B → ∅ (rate µ (cid:48) ) processes to the standard Lotka–Volterra kinetics (1). Yet this wouldinduce no qualitative changes as long as σ (cid:48) < µ and σ > µ (cid:48) ; one then simplyneeds to replace µ with the rate difference µ − σ (cid:48) , and σ with σ − µ (cid:48) .The associated rate equations, subject to mean-field mass action factoriza-tion for the non-linear predation process, and valid under well-mixed conditions4or spatially homogeneous time-dependent particle densities a ( t ) and b ( t ), read˙ a ( t ) = λ (cid:48) a ( t ) b ( t ) − µ (cid:48) a ( t ) , ˙ b ( t ) = − λ (cid:48) a ( t ) b ( t ) + σ (cid:48) b ( t ) , (2)with continuum reaction rates λ (cid:48) , µ (cid:48) , and σ (cid:48) . Since there exists a conservedfirst integral (Lyapunov function) K = λ (cid:48) ( a + b ) − σ (cid:48) ln a − µ (cid:48) ln b = const. forthis deterministic dynamics, the solutions to eqs. (2) are strictly periodic non-linear oscillations that precisely return to the system’s initial state. Althoughpopular in the fields of ecology and biology, the Lotka–Volterra model is alsooften criticized for being too simplistic and mathematically unstable. This is dueto several simplifying and likely unrealistic assumptions: First, the prey alwayshave a sufficient amount of food available, whence its depletion is neglected,and the prey’s nourishment source is not explicitly represented in the model.Second, the only source of food for the predator species is the prey, and itsconsumption is a necessary requirement for the predators’ reproduction. Third,there is no specified limit on the prey intake for the predators. Fourth, the rateof change of either population is directly proportional to its size. Finally, duringthe temporal evolution any environmental influence is assumed fixed in time,and crucial concepts such as trait inheritance, mutations, or natural selectionplay no role.In our research, we use Monte Carlo simulations for the stochastic Lotka–Volterra model based on the reactions (1) performed on a two-dimensionalsquare lattice with 512 ×
512 sites and periodic boundary conditions to fullyaccount for emerging spatial structures and internal reaction noise. We notethat we have also performed simulations on two-dimensional square latticeswith 256 ×
256 and 128 ×
128 sites; aside from overall noisier data, as one wouldexpect, we obtain no noticeable quantitative differences. Given that the cor-relation lengths ξ measured below are much smaller than these system sizes,this is not surprising. Due to our limited computational resources, we have notattempted runs on even larger systems. In the following, all listed Monte Carlodata and extracted quantitative results refer to 512 ×
512 square lattices. Also,for the reaction processes, we only consider the four nearest-neighbor sites, and5ave not extended interactions to larger distances. In our model, we imple-ment occupation number limitations or finite local carrying capacities; i.e., thenumber of particles on any lattice site is restricted to be either 0, if the siteis empty, or 1, if it is occupied by a predator or a prey individual. We shallexamine the population densities of each species, given by their total particlenumber divided by the number of lattice sites, and aim to quantify the ensuingoscillations and through characteristic observables that include their frequencyand attenuation, as well as typical population cluster sizes as determined bytheir spatial correlation length.The simulation algorithm for the death, reproduction, and mutual interac-tions of the prey and predator particles proceeds as follows [26, 7, 9]: For eachiteration, an occupied site is randomly selected and then one of its four adja-cent sites is picked at random. If the two selected sites contain a predator anda prey particle, a random number x ∈ [0 ,
1] is generated; if x < λ , the preyindividual is removed and a newly generated predator takes its place. Similarly,if the occupant is a predator, a random number x ∈ [0 ,
1] is generated, andthe particle is removed if x < µ . Yet if the initially selected occupant is a preyparticle and the chosen neighbor site empty, a random number x ∈ [0 ,
1] isgenerated; if x < σ a new prey individual is added to this site. One MonteCarlo Step (MCS) is considered completed when on average all particles haveparticipated in the reactions once.The variables that can be tuned in our simulations are: the system size L ,the initial predator density ρ A (0), the initial prey density ρ B (0), the predatordeath rate µ , the prey reproduction rate σ , the predation rate λ , and the numberof Monte Carlo steps. We chose the linear system size L = 512. Naturally onemust avoid starting the simulations from one of the absorbing states. For anynon-zero initial predator and prey density, the population numbers and particledistribution at the outset of the simulation runs influence the system merelyfor a limited time, and the final (quasi-)stationary state of the system is onlydetermined by the three reaction rates [18]. In our simulations, the rates µ and σ are kept constant for simplicity, while λ is considered to be the only relevant6 .0 0.1 0.2 0.3 0.4 0.5 ρ A ( t ) ρ B ( t ) λ =0 . λ =0 . λ =0 . Figure 1: Monte Carlo simulation trajectories for a stochastic Lotka–Volterra model on a512 ×
512 square lattice with periodic boundary conditions and restricted site occupancy shownin the predator ρ A ( t ) versus prey density ρ B ( t ) phase plane ( ρ A ( t ) + ρ B ( t ) ≤
1) with initialvalues ρ A (0) = 0 . ρ B (0), fixed reaction probabilities µ = 0 . σ = 1 .
0, and differentpredation efficiencies: (i) λ = 0 . λ = 0 .
18 (greenstars): direct exponential relaxation to the quasi-stationary state just above the extinctionthreshold in the predator-prey coexistence phase; (iii) λ = 0 . control parameter. The dynamical properties are generically determined by theratio of the reaction rates; the subsequent results apply also for different setsof µ and σ with appropriately altered predation rate λ . Since we only havetwo species, predators and prey, 0 ≤ ρ A + ρ B ≤ λ , the predators will gradually starveto death, and the remaining prey will finally occupy the whole system. On theother hand, when λ is large, there is a finite probability (in any finite lattice) thatall prey individuals would be devoured; subsequently the predators would dieout as well because of starvation. In fact, the absorbing extinction state is theonly truly stable state in a finite population with the stochastic dynamics (1).7owever, in sufficiently large systems, quasi-stable states in which both speciessurvive with relatively constant population densities during the entire simulationduration are indeed observed in certain regions of parameter space. Fig. 1 showsthree single-run simulation results, plotting the prey population density ρ B ( t )versus that of the predators ρ A ( t ) with the reaction probabilities µ = 0 . σ = 1 . λ as the only control parameter. We chose the initial population densities as ρ A (0) = 0 . ρ B (0) with the particles randomly distributed among the latticesites. With λ = 0 . ρ B → λ to 0 .
18 (green stars), just above the predator extinctionthreshold, the system relaxes exponentially to a quasi-stationary state withnon-zero densities for both species. For λ = 0 . λ c = 0 . ± .
3. Boundary Effects at a Coexistence / Predator Extinction Interface
Natural environments vary in space and boundaries are formed between dif-ferent regions, yielding often quite sharp interfaces, e.g., between river and land,desert and forest, etc. At the boundaries of such spatially inhomogeneous sys-tems, interesting phenomena may arise. In order to study boundary effectson simple predator-prey population dynamics, we split our simulation domaininto two equally large pieces with one half residing in the predator extinctionstate, and the other half in the two-species coexistence phase. We use a two-dimensional lattice with 512 ×
512 sites with periodic boundary conditions, andindex the columns with integers in the interval [0 , µ = 0 .
125 and σ = 1 . a) (b) (c) Figure 2: Snapshots of the spatial particle distribution on a 512 ×
512 lattice (with periodicboundary conditions) that is split into equally large predator extinction (left) and speciescoexistence (right) regions: prey are indicated in green, predators in red, white spaces inwhite. (a) Random initial distribution with densities ρ A = 0 . ρ B ; (b) state of the systemafter 1000 MCS, when it has reached a quasi-stationary state with uniform rates µ = 0 . σ = 1 .
0, while λ l = 0 . , λ r = 0 . , ×
50 area at the boundary. on all sites, we assign λ l = 0 . < λ c on columns [0 , λ r = 0 . > λ c for the columns on the “right”half with indices [256 , ρ A = 0 . ρ B . After the system has evolved for 1000 MCS,a quasi-steady state is obtained as shown in Fig. 2(b), and the close-up nearthe boundary (c). The predators are able to penetrate into the “left” absorbingregion by less than 10 columns, and no predator individuals are encounteredfar away from the active-absorbing interface. On the right half, we observea predator-prey coexistence state with the prey particles forming clusters sur-rounded by predators and predation reactions occurring at their perimeters.Since only the predator species is subject to the extinction transition intoan absorbing state, while the prey can survive throughout the entire simulationdomain, we concentrate on boundary effects affecting the predator population.We measure the column densities of predators ρ A ( n ), defined as the number of9redators on column n divided by L = 512, and record their averages from 1000independent simulation runs as a function of column index n . As shown in theinset of Fig. 3(a), ρ A ( n ) decreases to 0 deep inside the absorbing half of thesystem, and reaches a positive constant 0 . ± .
001 within the active region.The main graph focuses on the boundary region, where we observe a markedpredator density peak right at the interface (column 256). The predator densityenhancement at the boundary is obviously due to the net intrusion flow ofspecies A from the active subdomain with high predation rate into the predatorextinction region with abundant food in the form of the near uniformly spreadprey population. We also ran simulations for other predation rate pairs such as λ l = 0 . λ r = 0 . ρ A appeared on column 257 in that situationinstead of at n = 256).Fig. 3(b) shows the exponential decay of the predator column density ρ A ( n )as function of the distance | − n | from the boundary (located at n = 255)towards the “left”, absorbing side. A simple linear regression gives the inversecharacteristic decay length k = − . ρ A ( n ) neither fits exponential nor algebraic decay. Instead, ρ A reaches the asymptotic constant value 0 . ± .
001 deep in the coexistenceregion through an apparent stretched exponential form ρ A ( n ) ∼ e − ( n − l +0 .
195 with stretching exponent l ≈ . ξ , obtained from the equal-time correlation func-tion C ( x ), to characterize the spatial extent of these clusters. For species α, β = A, B , the (connected) correlation functions are defined as C αβ ( x ) = (cid:104) n α ( x ) n β (0) (cid:105)−(cid:104) n α ( x ) (cid:105) (cid:104) n β (0) (cid:105) , where n α ( x ) = 0 , α at site x [17]. For x = 0 and α = β , in a spatially homoge-neous system it is simply given by the density (cid:104) n A (cid:105) : C αα (0) = (cid:104) n A (cid:105) (1 − (cid:104) n A (cid:105) ).For | x | > (cid:104) n α ( x ) n β (0) (cid:105) is computed as follows: First choose a site, and then10
50 260 270 n ρ A (a) 0 5 10 15 20 l n ( ρ A ) (b) 0.0 0.4 0.8 1.2 log ( n − l o g ( − l n ( ρ A − . )) (c) n Figure 3: After the split system with rates µ = 0 . σ = 1 . λ l = 0 . , λ r = 0 . , ρ A ( n ) as functionof column index n ∈ [251 , L = 512 columns (data averaged over 1000 independent runs); (b) exponentialdecay of ρ A ( n ) from the boundary (located at n = 255) into the absorbing region with n ∈ [235 , k = − . ρ A decays to a positive constant value 0 . ± .
001 deep in the right coexistence region. Theblue dots display log ( − ln( ρ A − . ( n − l = 0 .
348 is obtained from linear data regression. a second site at distance x away from the first one. n α ( x ) n β (0) equals 1 only ifthe first site is occupied by an individual of species β , and the second one by anparticle of species α , otherwise the result is 0. One then averages over all sites.Here, we compute the predator correlations C AA ( x, n ) on a given column n , i.e., we only take the mean in the above procedure over the L = 512 siteson that column. The main panel in Fig. 4(a) shows the predator correlationfunction C AA ( x ) on column n = 274 with x ∈ [0 , C AA ( x ) graduallydecreases to zero. The inset presents the same data in a logarithmic scale,demonstrating exponential decay according to C ( x, n ) ∼ e − x/ξ ( n ) . Since thestatistical errors grow at large distances x , we only use the initial data pointsup to x = 6 for the analysis. Linear regression of ln( C AA ( x, n = 274)) over11 x C ( x ) (a) 255 260 265 270 275 n ξ (b) x l n ( C ( x )) Figure 4: (a) Main panel: the predator correlation function C AA ( x, n ) on column n = 274(data averaged over 1000 independent simulation runs). Inset: ln( C AA ( x, n )); the red straightline indicates a simple linear regression of the data points with x ∈ [1 , ξ ( n = 274) ≈ .
2; the error bars indicate the standard deviation. (b)Correlation length ξ ( n ) versus column number n , with ξ ( n ) defined as the negative reciprocalof the slope of ln( C AA ( x, n )). x ∈ [1 ,
6] gives ξ ( n = 274) ≈ .
2, indicated as red square in Fig. 4(b). Inthe same manner, we obtain the characteristic correlation lengths ξ ( n ) for eachcolumn n as shown in Fig. 4(b), starting at the interface at n = 256. Weobserve ξ ( n ) to increase by about a factor of four within the first ten columnsaway from the boundary, and then saturate at the bulk value ξ ≈ .
2. Near theabsorbing region, the predator clusters are thus much smaller, owing to the netflux of predators across the boundary into the extinction domain. These valuesof ξ are measured after the entire system has reached its (quasi-)steady stateafter 1000 MCS, and would not change for longer simulations run times. Wenote that the relationship between the correlation length ξ and the predationrate λ is manifestly not linear, i.e., a very large value of λ does not imply hugepredator clusters. We surmise that the cluster size remains finite even in thatscenario, and the predators would penetrate into the “left” absorbing region fora finite number of columns only. For sufficiently large domain size, the systemshould thus remain spatially inhomogeneous even for very high predation rates12
20 40 60 80 100 120 140 160 t [MCS] ρ A ( t ) column 256column 274 Figure 5: The temporal evolution of the average predator column densities ρ A ( t, n ) (averagedover 1000 independent runs) on columns n = 256 (blue dots) and n = 274 (red plus marks),with initial predator density ρ A (0) = 0 . µ = 0 . σ = 1 . λ l = 0 . n ∈ [0 , λ r = 0 . n ∈ [256 , λ . Finally, the dependence of the typical cluster size ξ ( n ) on column index n correlates inversely with the column density plotted in Fig. 3: High localdensity corresponds to small cluster size and vice versa. We note that theproduct ρ A ( n ) ξ ( n ) is however not simply constant across different columns;rather it is minimal near the boundary (at n = 256), then increases away fromthe interface, and ultimately reaches a fixed value within 10 columns inside theactive region.Spatially homogeneous stochastic Lotka–Volterra systems display dampedpopulation oscillations in the predator-prey coexistence phase after being ini-tialized with random species distribution, see, for example, the (red triangle)trajectory in Fig. 1 for predation rate λ = 0 .
4. We next explore the boundaryeffects on these population oscillations near the active-absorbing interface. Weprepare the system with the same parameters as mentioned above so that its“left” half is in the absorbing state while the “right” side sustains species coex-istence. The initial population densities are again set to ρ A (0) = 0 . ρ B (0),with the particles randomly distributed on the lattice. We then measure the col-umn predator densities as a function of time (MCS). Fig. 5 displays the temporalevolution of ρ A ( t, n ) on columns n = 256 and n = 274. We observe the oscil-13 .000 0.628 1.256 ω [ MCS − ] f ( ω ) (a)
256 264 272 280 288 n t c (b) Figure 6: (a) Fourier transform amplitude f A ( ω, n ) of the predator column density timeevolution on columns n = 256 (blue dots), 258 (green triangles up), 262 (red triangles down),266 (cyan squares), 270 (magenta stars), and 274 (black plus marks), with rates µ = 0 . σ = 1 .
0, and λ l = 0 . n ∈ [0 , λ r = 0 . n ∈ [256 , t c ( n ) on columns near the active-absorbing boundary, inferred from the peakwidths in (a), with the error bars representing the standard deviation. lations on the column closest to the interface to be strongly damped, whereasdeeper inside the active region the population oscillations are more persistentand subject to much weaker attenuation. Both column densities asymptoticallyreach the expected quasi-steady state values.In order to determine the dependence of the local oscillation frequencieson the distance from the active-absorbing interface, we compute the Fouriertransform amplitude f A ( ω, n ) = | (cid:82) e − iωt ρ A ( t, n ) dt | of the column density timeseries data by means of the fast Fourier transform algorithm for n ∈ [256 , ρ A ( t, n ) ∼ e − t/t c ( n ) cos(2 πt/T ( n )), we may then identify the peak position of f A ( ω, n )with the characteristic oscillation frequency 2 π/T ( n ), and the peak half-widthat half maximum with the attenuation rate or inverse relaxation time 1 /t c ( n ).We find that the oscillation frequencies are constant except for the column atthe boundary ( n = 256), which shows a very slight enhancement. We concludethat the presence of the extinction region does not markedly affect the frequency14 a) (b) (c) Figure 7: Snapshots of the distribution of predator (red) and prey (green) particles after thesystem has evolved for 1000 MCS with rates µ = 0 . σ = 1 .
0, and λ switched alternatinglybetween the values 0 . . ×
512 system is periodically divided into successively smallersquare patches with lengths 256 (a), 64 (b), and 16 (c), respectively. The square subdomainsdominantly colored in green reside in the extinction state ( λ = 0 . λ = 0 . of the population oscillations in the active regime. In contrast, the attenuationrate increases by a factor of three within about 20 columns in the vicinity ofthe interface, as demonstrated in Fig. 6(b). Beyond n ≈
278 in the coexistenceregion, the relaxation time assumes its constant bulk value.
4. Checkerboard Division of the System
To further explore boundary (and finite-size) effects in spatially inhomoge-neous Lotka–Volterra systems, we proceed by successively dividing the simula-tion domain into subdomains in a checkerboard pattern, setting the predationrate to two distinct values in neighboring patches, and thus preparing themalternatingly in either the active coexistence or absorbing predator extinctionstates. Fig. 7(a) shows a case when the system is split into four subregions with σ = 1 . µ = 0 . λ = 0 . . × ×
512 sites is respectivelysplit into 8 × ×
32 square patches: If for a given box λ is set to 0 . λ = 0 . λ switching between 0 . . .
1, but stays roughly the same on the subdomains where λ = 0 .
8. The predator density in contrast increases in both the active andabsorbing regions as the subdivision proceeds. We have also confirmed thatthese changes in the total population densities naturally become less significantif the two different predation rate values are chosen closer to each other.In Fig. 8, we plot the total (summed over all subdomains) predator and preypopulation densities ρ in the simulation domain split into N × N checkerboardpatches, as functions of log N . Here, N = 1 corresponds to the situation stud-ied in section 3, where the system was divided into two rectangular subdomains.The other values of N = 512 , , , , , , , /N . The mean population density ρ shown for eachdata point represents an average of 1000 independent simulation runs; the asso-ciated statistical error was very small, with a standard deviation of order 10 − .As apparent in the data, the overall population predator density ρ A monotoni-cally increases with growing number N of subdivisions, while the prey density ρ B decreases.We also performed the analogous sequence of measurements for other pairsof predation rate values. For instance, with checkerboard subdomains with λ = 0 . . N ρ ρ A ( N ) ρ Ar ρ B ( N ) ρ Br Figure 8: Total population densities for predators (red triangles) and prey (green squares)versus number of checkerboard-patterned subdivisions N of the simulation domain, after thesystem has evolved for 1000 MCS and reached a quasi-stationary state, with reaction rates µ = 0 . σ = 1 .
0, and λ alternatingly switched between 0 . .
8. For comparison,the graph also shows the total quasi-steady state population densities for predators (blackplus) and prey (blue star) in a system with randomly assigned predation, drawn with equalprobability from a bimodal distribution with values λ = 0 . . the population density changes with increasing N are less pronounced than inFig. 8, and ρ A , ρ B acquire maximum and minimum values at N = 256 ratherthan 512. The origin of this slight shift can be traced to the fact that thepredator correlation length is of order one lattice constant at the boundary ofthe λ = 0 . / . . / . λ = 0 . . N = 512 system, for which these two predation ratesare alternatingly assigned to the lattice sites in a periodic regular manner. Wewould expect the population densities in these two distinct systems to reachequal values if the associated correlation lengths at the boundaries were largecompared to the lattice constant, which is however not the case here.17 . Conclusion In this work, we have focused on studying boundary effects in a stochasticLotka–Volterra predator-prey competition model on a two-dimensional lattice,by means of detailed Monte Carlo simulations. We first considered a systemsplit into two equally large parts with distinct non-linear predation rates, suchthat one domain is set to be in the predator extinction state, while the other oneresides in the two-species coexistence phase. We have primarily addressed theinfluence of such an absorbing-active separation on both populations’ densityoscillations as function of the distance from the boundary.We find a remarkable peak in the column density oscillation amplitude ofthe predator population, as shown in Fig. 3(a), which reflects its net steadyinflux towards the absorbing region. Correspondingly, the predator correlationlength that characterizes the typical cluster size reaches a minimum value at theboundary, see Fig. 4(b). The population oscillation frequency there shows onlysmall deviations from its bulk value, while the attenuation rate is locally stronglyenhanced, see Fig. 6(b), inducing overdamped relaxation kinetics. Overall, theecosystem remains stable.Furthermore, upon splitting the system successively into more pieces in acheckerboard fashion, the observed boundary effects become less significant,and as demonstrated in Fig. 8, the overall population densities acquire valuesthat are close to those in a disordered system with randomly assigned predationrates drawn from a bimodal distribution.
ReferencesReferences [1] Lotka A J 1920 Undamped oscillations derived from the law of mass action
J. Am. Chem. Soc. Nature
558 183] Matsuda H, Ogita N, Sasaki A and Sat¯o K 1992 Statistical mechanics ofpopulation – the lattice Lotka–Volterra model
Progr. Theor. Phys. Phys. Rev. E Phys. Rev. E SIAM Review J. Chem. Phys.
Physica A Phys.Rev. E Physica A Physica A Phys.Rev. E Phys. Rev. E Phys. Rev. E Phys. Rev. E J.Stat. Phys.
J. Phys.Cond. Matt. Phys. Biol. Phys. Rev. Lett. Biometrika Statistical Mechanics of Biocomplexity
Phys. Rev. E Ann. Eugen. Phys. Rev. Lett.
Phys. Rev. Lett.
Phys. Rev. Lett.56