Oscillations and patterns in interacting populations of two species
Interacting populations often create complicated spatiotemporal behavior, and understanding it is a basic problem in the dynamics of spatial systems. We study the two-species case by simulations of a host–parasitoid model. In the case of co-existence, there are spatial patterns leading to noise-sustained oscillations. We introduce a new measure for the patterns, and explain the oscillations as a consequence of a timescale separation and noise. They are linked together with the patterns by letting the spreading rates depend on instantaneous population densities. Applications are discussed.
pacs:87.23.Cc 02.50.Ey 82.20.-w 87.18.Hf
A fundamental aim in studying population dynamics is to understand species interactions. A paradigmatic, still interesting, case is that of two species (1), where one feeds or lives from the other. These may be predators and their prey or parasitoids living on the expense of hosts. If the interaction is strong and if the parasitoid or predator is specialized, its growth will have a delayed negative feedback as its resource is diminished. In such cases, oscillations are an inherent feature (2). Many surface reactions also have similar dynamics (3).
The classical ways of looking at such systems assume fully stirred populations. The encounters between individuals are assumed to be proportional to the product of their densities, analogously to the mass action principle in chemistry (1). This assumption is also at the heart of the mean-field -like Lotka–Volterra equations. In general, however, spreading and interaction are restricted in space. In this case, correlated structures arise and the assumption about complete mixing no longer holds (4).
If the parasite abundance is small, any feedback effect is weak. Population sizes then show no oscillations, and the predating species is locally concentrated in a cluster-like arrangement. This has been theoretically studied in Ref. (5). With strong feedback spatiotemporal patterns emerge in a multitude of forms. These include disordered flame-like patterns (6); (7); (8), and ripple-like spatiotemporal ones (9). There is also a large body of work on similar patterning in individual-based models (e.g. (10)), in statistical physics (11); (12), and in calcium concentration oscillations in living cells (13). A paradigmatic pattern-forming system is the complex Ginzburg-Landau equation (CGLE) (14), which exhibits spiral-like geometries.
Voles in Northern Britain (15), mussels in the Wadden sea (16), and lemmings in Northern Europe (17), are good examples of empirical observations of such patterning. These involve either predation or being predated. Spatial structures weaken the interactions since species tend to be aggregated within themselves. They also provide the prey a refuge since around the prey there are less predators. Therefore, spatial inhomogeneity can stabilize the dynamics and promote coexistence (7); (4); (18).
Here, we analyze the full spatiotemporal dynamics of two interacting species using instantaneous configurations and time series of the population densities. We introduce a measure for the level of patterning in such systems. When the patterns form, one observes persistent oscillations in the population densities. We show that the underlying dynamics follows a particular logic: it originates in the response of the rates to changes in instantaneous densities, and the emerging system proves different from the limit cycle in Lotka–Volterra systems, or recent developments where three-species models have been mapped (19); (20) to CGLE. The present mechanism works by the interplay of oscillatory transients to a stable fixed point and stochasticity. The response of the interaction rates is due to spatial correlations. The mechanism is novel, and does not in fact need any non-linearities to work, which will become evident below based on simulations and effective equations describing them.
We study a host–parasitoid model in discrete time and space. It is inspired by (21), but has a wider interaction range as in the incidence function models of metapopulation dynamics ((22)). The model describes annual host–parasitoid dynamics on a two-dimensional square lattice (23). At each time step, a site can be either empty (in state ) or populated by a host without (state ) or with parasitoids (state ). Transitions between the states are cyclic, , neglecting possible spontaneous deaths of non-parasitized hosts assumed to be rare, for simplicity. Although this means that hosts live forever if there are no parasitoids, this is not a serious restriction; in coexistence it boils down to assuming faster extinction for the parasitoids than for the hosts. The model can also be described as an SIR model with rebirth.
At each site transition probabilities depend on the surrounding populations through the connectivity
of site with respect to species ( or ). Here if the state of is , and 0 else. The kernel has an exponential decay with the scale and is normalized by . Dispersal lengths are chosen since these are biologically motivated (22) and lead to generalizability.
In a timestep, the transition takes place with probability and with probability (in the parameter range of interest, practically always). Parasitized hosts may die () with probability irrespective of the surroundings. Note that they do not reproduce. Periodic boundary conditions and parallel updates are used. There are two absorbing states, an empty lattice () and one full of hosts. Fig. 1 (a) shows an example with coexistence.
One finds moving regions predominantly in one of the three states. To quantify these patterns, we define the dominance regions (Fig. 1b) as follows. By smoothing one obtains continuous densities . For each site, and are compared to the space–time averages and . The densities are positive, lying in the first quadrant of , divided into three regions shown in Fig. 1c. The site at time is then defined to belong to a domain according to the region (, , or ) containing . In essence, the regions coarse-grain on a scale . In this regime, they are insensitive to changes in .
The domains are separated by walls, joining at triple points, vortices (6). A vortex has a sign (), if one encounters the domains in the order following a small cycle around the vortex counter-clockwise (clockwise). Pairs of vortices of opposite signs are created and annihilated together. The domains rotate around the vortices, which are relatively stable. In other words, the species invade the appropriate neighboring domains so that the walls rotate around the vortices. Similar structures have been identified earlier in related systems (e.g. (6); (11)). In three dimensions, the vortices generalize to strings (6).
First, consider static measures such as the domain wall length from source to sink vortex. It has an exponential distribution, whose mean is drawn for different parameters in Fig. 1d. Its ratio to its counterpart in uncorrelated random arrangements with the same densities is shown. Patterns and oscillations lead to walls with lattice units (l.u.). This is more than ten times larger than the smoothing width and also several times that in the random arrangements ( l.u.). The coarse-graining gives a measure of patterns distinguishing between uncorrelated and patterned states.
Next, turn to the spatially averaged densities and . Fig. 2a shows them as a function of time in three cases: (i) a non-patterned state with a small parasitoid population, (ii) a state with patterns, and (iii) a small subsystem () out of a large system () with patterns. In the patterned systems, there is a high-frequency oscillation matching the angular velocity of single vortices and a slow variation of the amplitude. Below, we explain both as a consequence of a timescale separation, connect them to the patterns, and explain why the oscillations do not conform to the usual limit cycles.
To build a description of the dynamics of the model in a novel fashion using aggregated variables, consider Poincaré maps (Fig. 3). For large enough systems and are unique functions of and up to noise. Also by attractor reconstruction (24) we find that the full system with degrees of freedom coarse-grains into a two-dimensional one. The points in the maps lie close to a two–dimensional surface, and for large enough (with many patterns in the system) even on the tangential plane through the average . Based on these numerical observations, the dynamics is linear in and :
In other words, the observations imply that even though the dynamics is expected to be non-linear based on the MF approximation, in a large system with many patterns the possible non-linearities self-average out.
It is then useful to consider the expansion around the average – a linear iterative map. In oscillatory cases, its eigenvalues form a conjugated pair . These are associated with two timescales, the period and the decay rate of the amplitude. Their ratio
tells whether dynamics is oscillatory or “just noisy”. indicates patterned systems and oscillatory dynamics. Fig. 4 shows the typical behaviors of the two kinds of dynamics observed depending on the presence or absence of patterns, and whether one adds noise (as additional Gaussian uncorrelated noise terms on the RHS of Eq. (2)) to the coarse-grained dynamical system to mimic the finite- simulations of the full spatial system. Note that in all cases the fixed point is attractive.
So far we have given separately a temporal and a spatial diagnosis of the pattern dynamics. Next we link these together. For small, they equal the spreading probabilities, and the dynamics can be written as and , where the influence of the connectivities is expressed by the correlation functions . A corresponding non-spatial approximation is
This is an approximation of the usual MF form with the interaction parameters and generalized to be arbitrary functions of the instantaneous densities. They can be non-linear and they do not have to conform to the standard MF nor to any ad-hoc approximations (25). By an expansion of Eqs. (4) around the fixed point one arrives at Eq. (2) with
where , , and their derivatives are evaluated at the fixed point. The derivatives are necessary for consistency. The matrix elements and the densities and are measured from the simulations. Since and are parameters, omitting the derivatives would make Eqs. (5) overdetermined and thus unsatisfiable. By keeping them, there are four equations and four unknowns to be solved.
The effect of the nonzero derivatives is best illustrated by a phase diagram, Fig. 2. In the spatially extended system, there are three qualitative phases: the extinction of the parasitoids, non-oscillatory and oscillatory coexistence. Except for the extinction, the boundaries are not sharp. Instead, there is a transition zone, defined via the timescale ratio (Eq. (3)) as the region where . In MF, there is also a fourth phase, absent here: oscillatory coexistence in a limit cycle. The phase structure resembles that in earlier work on a related model with only nearest-neighbour spreading (26); (27); (28). There, as well, oscillatory and non-oscillatory phases are recovered, the latter identified as a limit cycle using the pair approximation. Based on our findings, it could also be noise-sustained.
Let us now compare the explained mechanism with recent approaches. A possibility is to make the angular velocity of the oscillation either amplitude- or phase-dependent ((29); (30)). However, Eq. (2) does not allow for either dependence. Another one is to map the population model (19); (20) to CGLE (14). An unstable fixed point is necessary for the mapping, yielding a limit cycle.
To conclude, we have studied spatiotemporal dynamics of a model of two kinds of interacting particles, in biological terms hosts and parasitoids. A large parasitoid population creates patterns and noisy oscillations of population sizes. We have introduced a new measure for the patterns, and explained the noisy oscillation as a consequence of a time-scale separation. In other words, even with a limit cycle at the well-mixed limit, the spatial case has stable dynamics with long-lived oscillatory transients. This is due to spatial correlations making the spreading rates functions of the instantaneous population densities. Since the type of oscillation determines its properties (e.g. the fluctuating amplitude), which in turn affect vulnerability to extinction, the distinction is important. The connection offers a shortcut to study the effect of, e.g. environment: it could be related directly to the matrix elements in (2), in contrast to a full form of the interactions, lightening the analysis. We expect that the observation of patterns and oscillations arising from local dynamics and self-averaging - in finite systems since the noise amplitude depends on system size - will find other applications beyond the biology-inspired model. They are not restricted to only two species or two-dimensional systems, since the analysis can be carried out also for more complicated cases. There is no restriction to cyclic dynamics either, nor to discrete-time systems since continuous-time ones can be handled by considering snapshots taken at regular intervals. Further examples of applications include chemical reactions on surfaces (3), and metapopulations on disordered and scalefree landscapes.
Acknowledgments. Ilkka Hanski and Otso Ovaskainen are thanked for stimulating discussions and Lasse Laurson for assistance. This work was supported by the Academy of Finland through the Center of Excellence program (M.A. and M.P.) and Deutsche Forschungsgemeinschaft via SFB 611 (M.R.). The authors thank the Lorentz center for kind hospitality.
- A.J. Lotka, J. Am. Chem. Soc. (1920).
- P. De Bach and H.S. Smith, Ecology 22, 363 (1941).
- V. P. Zhdanov, Surf. Sci. Rep. 45, 231 (2002).
- R.V. Solé and J. Bascompte, Self-Organization in Complex Ecosystems, Princeton University Press, Princeton (2006).
- O. Ovaskainen and S. Cornell, Proc. Natl. Acad. Sci. 103, 12781 (2006).
- K.-i. Tainaka, Phys. Rev. E 50, 3401 (1994).
- F. Zhang, Z. Li, and C. Hui, Ecological Modelling 193, 721 (2006).
- M. Mobilia and I. T. Georgiev and U. C. Täuber, J. Stat. Phys. 182, 447 (2007).
- T. Nguyen-Huu, C. Lett, J.C. Poggiale, and P. Auger, Ecological Modelling 197, 290 (2006).
- A. M. de Roos, E. McCauley, and W. G. Wilson, Proc. Biol. Sci. 246, 117 (1991).
- G. Szabó and A. Szolnoki, Phys. Rev. E 65, 036115 (2002).
- K. Wood, C. Van den Broeck, R. Kawai, and K. Lindenberg, Phys. Rev. Lett. 96, 145701 (2006).
- D. E. Clapham, Cell 80, 259 (1998).
- I. S. Aranson and L. Kramer, Rev. Mod. Phys. 74, 99 (2002).
- X. Lambin, D.A. Elston, S.J. Petty, and J.L. MacKinnon, Proc. R. Soc. Lond. B 265, 1491 (1998).
- J. van de Koppel, M. Rietkerk, N. Dankers, and P.M.J. Herrman, Am. Nat. 165, E66 (2005).
- I. Hanski, H. Henttonen, E. Korpimäki, L. Oksanen, and P. Turchin, Ecology 82, 1505 (2001).
- C. Briggs and M.F. Hoopes, Theor. Pop. Biol. 65, 299 (2004).
- T. Reichenbach, M. Mobilia, and E. Frey, Nature 448, 1046 (2007).
- T. Reichenbach, M. Mobilia, and E. Frey, e-print, arXiv:0801.1798.
- M.P. Hassell, H.N. Comins, and R.M. May, Nature 353, 255 (1991).
- I. Hanski, Nature 396, 41 (1998).
- One can also study in a similar way a realistic example of a metapopulation landscape, the Åland archipelago on the Baltic Sea (manuscript in preparation).
- N.H. Packard, J.P. Crutchfield, J.D. Farmer, and R.S. Shaw, Phys. Rev. Lett. 45, 712 (1980).
- M. Pascual, P. Mazzega, and S.A. Levin, Ecology 82, 2357 (2001).
- J.E. Satulovsky and T. Tomé, Phys Rev. E 49, 5073 (1994).
- T. Tomé and K.C. de Carvalho, J. Phys. A Math. Gen. 40, 12901 (2007).
- E. Arashiro, A. L. Rodrigues, M. J. de Oliveira, and T. Tomé, Phys. Rev. E 77, 061909 (2008).
- R. Abta, M. Schiffer, and N. M. Shnerb, Phys. Rev. Lett. 98, 098104 (2007).
- R. Abta and N. M. Shnerb, Phys. Rev. E 75, 051914 (2007).