Lotka–Volterra predator-prey model with periodically varying carrying capacity

Mean-field rate equations for stochastic dynamical reaction systems are approximate deterministic equations that aptly describe a well-mixed setup. Even though they neglect spatial correlations and temporal fluctuations, they are often useful to gain intuition on the system’s expected behavior. In Sec. III, we shall compare the results obtained with the mean-field equations with the Monte Carlo simulation data from the full stochastic model (1).

For the Lotka–Volterra predator-prey competition model (1), the classical mean-field rate equations that describe the time evolution of the mean predator and prey densities $a(t)$ and $b(t)$ read

These rate equations can be understood as representing gain / loss terms for reactions that increase / decrease the population densities. Linear stability analysis of this system shows that the system exhibits a species coexistence fixed point, and numerical integration of these equations leads to oscillatory behavior, namely neutral limit cycles. We will perform our analysis on the more generalized Lotka–Volterra model with a growth-limiting factor for the prey species (for reviews, see Refs. [51, 5]).

The original Lotka–Volterra rate equations can be generalized by including a growth-limiting factor $1-a(t)/K_{1}-b(t)/K_{2}$, where $K_{1}$ and $K_{2}$ respectively represent the (global) carrying capacities induced by prey-predator and prey-prey resource competition. For simplicity, we set $K_{1}=K_{2}$, since this does not change the qualitative behavior of the system on the mean-field level; this implies the modified set of rate equations

where $K$ denotes the (global) carrying capacity. Their mean-field character resides in the assumed factorization for the non-linear predation reaction with rate $\lambda$ of a two-point correlation function into a mere density product, which assumes statistical independence and the absence of correlations. The growth-limiting factor is used to model limited finite resources, and vanishes if $a(t)+b(t)=K$. In that case, the prey density’s temporal derivative becomes negative, indicating a strictly decreasing prey population. We remark that adding an explicit growth limiting term for the predator density is not required since the predators’ growth is determined by the prey density. Hence, if the prey species has a growth limiting factor, this will indirectly constrain the predator population abundance as well.

The stationary states of this system are given by constant solutions to (3). This results in three fixed points $(a^{*},b^{*})=\{(0,0),(0,K),(a_{0},b_{0})\}$, where

The solution $(0,0)$ represents total population extinction. At the fixed point $(0,K)$, the predator species goes extinct while the prey species fills the entire system to full capacity $K$. Finally, the solution $(a_{0},b_{0})$ with non-zero densities for both species represents predator-prey coexistence. Note that $a_{0}>0$ requires $\mu/\lambda<K$.

Next we consider the (linear) stability of these solutions, which is achieved by linearizing (3) around the three distinct stationary states. Shifting the densities by their stationary solutions $a(t)=a^{*}+\delta a(t)$, $b(t)=b^{*}+\delta b(t)$, inserting this transformation into the original rate equations, and keeping only terms linear in the small deviations $(\delta a(t),\delta b(t))$, we obtain the matrix equation $\mathbf{\dot{x}=Jx}$, where $\mathbf{x}=\bigl{(}\delta a(t)\ \delta b(t)\bigr{)}^{T}$, the dot represents the time derivative, and the Jacobian matrix $\mathbf{J}$ is explicitly given by

The dynamical behavior of the system in the vicinity of a fixed point follows from the eigenvalues $\epsilon_{\pm}$ of the Jacobian matrix at each stationary point. First let us consider the extinction fixed point $(0,0)$ with associated eigenvalues $(\epsilon_{-},\epsilon_{+})=(-\mu,\sigma)$. Both eigenvalues are real, indicating exponential behavior near the fixed point. Yet the extinction stationary point is linearly unstable in mean-field theory against prey growth, since $\epsilon_{+}=\sigma$ is positive. While this result is intuitive considering the fact that any small deviation in the prey density leads to exponential growth of the prey, we recall that in the original stochastic model, in any finite system total extinction represents the only asymptotically stable stationary absorbing state. Next, the eigenvalues for the predator extinction fixed point $(0,K)$ are $(\epsilon_{-},\epsilon_{+})=(\lambda K-\mu,-\sigma)$, which are also both real. This stationary state is only stable with respect to small perturbations if $\lambda<\lambda_{c}=\mu/K$. Finally, the two-species coexistence stationary point $(a_{0},b_{0})$ has associated eigenvalues

For $\lambda_{s}=(\mu/2K)\left(1+\sqrt{1+\sigma/\mu}\right)>\lambda>\lambda_{c}$, both eigenvalues are real and negative, and hence the stationary point is stable and small perturbations exponentially relax back towards it. If $\lambda>\lambda_{s}$, the eigenvalues acquire complex conjugate imaginary components with a negative real part, indicating that the stationary point is still stable, but the system exhibits decaying oscillations in its vicinity. Yet for $\lambda<\lambda_{c}$ the eigenvalues are both real with $\epsilon_{-}<0$ and $\epsilon_{+}>0$ assuming opposite signs. Consequently, the stationary solution $(a_{0},b_{0})$ turns into an unstable saddle point.

This analysis demonstrates that the mean-field rate equations (3) predict a continuous active-to-absorbing state transition at $\lambda=\lambda_{c}$. The absorbing state is the predator extinction phase $(0,K)$ which is stable only for $\lambda<\lambda_{c}$. The active phase is the species coexistence phase $(a_{0},b_{0})$ which only exists and is then stable for $\lambda>\lambda_{c}$. This fixed point is a stable node for $\lambda<\lambda_{s}$ and becomes an attractive focus for $\lambda>\lambda_{s}$. The active-to-absorbing phase transition describing predator extinction is also observed in spatially extended stochastic systems. Away from criticality the system’s behavior only changes quantitatively relative to the mean-field analysis. Near the phase transition the critical exponents governing the model’s dynamical scaling laws acquire substantial corrections due to fluctuations in dimensions $d\leq d_{c}=4$. For a more thorough review of the stochastic Lotka–Volterra predator-prey model in a static environment, we refer to Refs. [50, 51, 52].

In this section, we describe results obtained from numerically integrating the coupled rate equations (3) subject to a periodically switching carrying capacity. As depicted in Fig. 1, the carrying capacity $K(t)$ is taken to be a rectangular time signal ranging between the low and high values $K_{-}$ and $K_{+}$, and with full switching period $T_{k}$ (i.e., from $K_{\mp}$ back to $K_{\mp}$). This functional variation of the carrying capacity does of course not constitute a realistic model for species interacting in nature since food resources do not change in a discontinuous manner. However, it can be argued that seasonal changes lead to a sudden carrying capacity drop / increase between winter and summer, as resource availability may seasonally vary. The following results will later be utilized to highlight the differences between the mean-field approximation and the stochastic lattice model. We remark that a full quantitative comparison between the two models is uninformative due to the fact that in the lattice model one prescribes microscopic reaction probabilities, whereas in the mean-field system one controls the effective macroscopic reaction rates. A thorough quantitative analysis would require fitting the stochastic lattice data to the mean-field results in order to extract the effective (and usually scale-dependent) macroscopic rates. Here, we are not interested in the detailed quantitative differences between the lattice and mean-field models. Rather we shall focus on the qualitative distinctions between the two models, and will specifically highlight features predicted by the mean-field equations that are not present in the stochastic lattice system.

The mean-field equations were numerically integrated by employing a fourth-order Runge–Kutta scheme with (dimensionless) time increment $\Delta t=0.01$, i.e., $t_{0}=100\,\Delta t$ sets the basic unit time scale relative to which all times and inverse rates will henceforth be measured in this section. We set the initial conditions to $\rho_{a}(0)=\rho_{b}(0)=0.5$ and $K(0)=K_{-}$, and have confirmed that our results do not depend on these chosen initial values. Figure 2 displays the resulting predator and prey densities $\rho(t)$ as functions of time. For both switching periods $T_{k}=60$ and $T_{k}=80$, we clearly observe period-doubling effects in the time traces. This is further confirmed by the Fourier transforms of these temporal evolutions shown in Fig. 3. For $T_{k}=60$, the highest Fourier peak occurs at a period $t=120$ indicating period-doubling. However, an additional smaller peak emerges at $t=240$, reflecting that the density repeats after four switching periods of the carrying capacity, suggesting even the presence of a period-quadrupling effect. Similar period-doubling is visible for $T_{k}=80$, but no period-quadrupling is discernible. Further increase in the carrying capacity period evidently eliminates period-doubling phenomena as shown in Fig. 3c. We detect the highest peak in the density Fourier transforms for $T_{k}=200$ at $t=T_{k}/3$, a harmonic of the driving period. This feature is in fact also observed in the lattice model, in contrast to the period-doubling at smaller periods $T_{k}$, for which we shall find that the internal reaction noise in the stochastic model washes away these intriguing non-linear effects.

The stationary mean-field population densities in the coexistence phase are given in Eq. (4). For an environment where $K$ periodically switches between two constant values $K_{-}$ and $K_{+}$, the long-time behavior of the system depends on these stationary densities. If the period of the oscillating environment is sufficiently long such that the system reaches the stationary state for either $K$ value, then the densities can effectively be described as oscillating between two constant values with the same period $T_{k}$ as the carrying capacity. In that case, the averages of the predator and prey densities over one period can simply be approximated by the arithmetic means $(\tilde{a},\tilde{b})$ of the two stationary values $(a_{-},b_{-})$ and $(a_{+},b_{+})$ pertaining to $K=K_{-}$ and $K=K_{+}$, respectively. Thus we obtain

We rewrite the mean predator density in terms of an equivalent time-averaged effective carrying capacity $K^{*}$ defined through

Comparison with the explicit result (7a) yields

which reduces to the rate-independent harmonic average ${\bar{K}}=2K_{-}K_{+}/(K_{-}+K_{+})$ for large $K_{-},K_{+}\gg 1$. Hence, in the slow-switching regime, the system can be described as oscillating around the average population densities corresponding to the constant rate-dependent effective carrying capacity $K^{*}$.

Through numerical integration of the mean-field rate equations, we tested the harmonic average hypothesis for different switching periods, and confirmed Eq. (9) in the slow-switching regime. We note that this comparison is facilitated for the mean-field model compared to the stochastic lattice system because we have exact formulas available for the stationary density values, and $K$ is not required to be an integer. Figure 4 shows the comparison of the numerically obtained population densities with periodically varying $K(t)$ with the corresponding stationary values obtained with a simple harmonic average of the carrying capacities and the rate-dependent effective carrying capacity (9). Interestingly, computing the stationary prey density from the straightforward harmonic carrying capacity average ${\bar{K}}$ yields accurate results for a large range of switchting periods, as is apparent in Fig. 4b. This is due to the fact that for a static carrying capacity, the prey density oscillates about its stationary value, and the fluctuations about it almost precisely average out, as verified in Fig. 5. Since within the mean-field framework, $b^{*}$ and hence $\tilde{b}$ do not depend on $K$, any equivalent carrying capacity would work for the prey population. The predator density also follows the harmonically averaged carrying capacity for small periods, see below; and is indeed aptly captured by the rate-dependent equivalent carrying capacity (9) for large switching periods. For intermediate periods $T_{k}$, we observe a non-monotonic crossover regime with a large resonance-like spike, see Fig. 4a. We verified that these findings do not depend on the initial conditions of the system.

The coupled mean-field rate equations (3) suggest that in the fast-switching regime, both species’ densities oscillate about values that are equal to the stationary population densities for an equivalent carrying capacity ${\bar{K}}$ that is just the harmonic average of $K_{-}$ and $K_{+}$. This follows from the fact that the prey density rate equation (3a) depends explicitly on $1/K$. Based on these observations, we construct an ansatz for the long-time behavior of both species’ densities as follows.

We first shift time according to $t\to t-NT_{k}$, where $N$ is a large integer such that at $t=NT_{k}$ the system has reached its quasi-stationary state. Hence this time axis shift defines $t=0$ to be the start of an environmental cycle in the long-time regime. If the system is thus initialized at the onset of the low carrying capacity state, i.e., at $t=0$ it just switched from $K_{+}$ to $K_{-}$, then at $t=T_{k}/2$ it will flip back from $K_{-}$ to $K_{+}$, and that cycle repeats at $t=T_{k}$. We now derive an approximate solution that describes the densities in one cycle $t\in[0,T_{k}]$. Henceforth we shall refer to the region $t\in[0,T_{k}/2]$ as $\mathcal{T_{-}}$, and the time interval $t\in[T_{k}/2,T_{k}]$ as $\mathcal{T_{+}}$.

Since the prey density exhibits a discontinuity in its first time derivative at $t=T_{k}/2$, it can be described by a piece-wise function. In the fast-switching regime, we may apply a short-time Taylor expansion for the population dynamics, and retain only the linear term. The absolute values of the prey density slope in the intervals $\mathcal{T_{-}}$ and $\mathcal{T_{+}}$ must be the same, due to the fact that the prey density is periodic, $b(T_{k})=b(0)$, and continuous at the jumps in between these two regions. In $\mathcal{T_{-}}$ the system is in the low carrying capacity state, therefore the prey density is a decreasing function of time, and its slope should be negative. For $t\in\mathcal{T_{+}}$, the prey density has a positive slope, since now the system is in the high carrying capacity state. These considerations motivate the following simple ansatz for the prey density,

which is numerically verified in Fig. 6.

The prey density is continuous at the boundary $t=T_{k}/2$, whence $b_{2}=b_{1}-\alpha T$. Moreover, in the fast-switching regime, the density variations of both species over one period of the carrying capacity should be assumed to be small relative to their average values. Consequently, $1/K(t)$ is the only significant term when averaging Eq. (3b). Its average leads to an equivalent carrying capacity that is equal to the harmonic average ${\bar{K}}$. Therefore, the system reaches a quasi-stationary state, where both densities oscillate around their stationary values for an equivalent carrying capacity ${\bar{K}}$. The temporal average of Eq. (10) needs to be $b_{0}$. Imposing this condition, we obtain

the slope constant $\alpha$ shall be determined later.

The rate equation for the predator density (3a) may now be cast into a more suggestive form,

which indicates that the extrema of $a(t)$ occur at times when $b(t)=b_{0}$. Using the ansatz (11), this happens at $t=T_{k}/4$ and $t=3T_{k}/4$. Equation (12) can then be integrated to solve for the predator density,

where $A$ and $A^{\prime}$ are integration constants. Since the predator density is required to be continuous at $t=T_{k}/2$, one arrives at the relation $A^{\prime}=A\,e^{\lambda\alpha T_{k}^{2}/4}$, which yields the approximate predator density solution

The average of the predator density over one cycle of environmental switching then becomes

Under the assumption of fast environmental switching, $T_{k}$ should be the smallest time scale in the system, and the explicit form of Eq. (14) suggests that the fast-switching regime is quantitatively delineated by $T_{k}\sqrt{\alpha\lambda}\ll 1$. The still undetermined parameter is the (initial) slope of the prey density $\alpha=|\dot{b}|_{0}$. To zeroth order in $\alpha$, either immediately from Eq. (13) or by expanding Eq. (14) in $T_{k}\sqrt{\alpha\lambda}$, gives the simple result

Since this average must equal the stationary value of predator density for a harmonically averaged carrying capacity, we may fix the integration constant

to leading order in an expansion in powers of $T_{k}\sqrt{\alpha\lambda}$.

The left-hand side of Eq. (3b) equals the constant slope of the prey density under the fast-switching approximation. Since $t<T_{k}$, we also have $t\sqrt{\alpha\lambda}\ll 1$, and with $a(t)=A+O(T_{k}\sqrt{\alpha\lambda})$ one has $b(t)=b_{0}+O(T_{k}\alpha)$. Upon inserting these asymptotic values into (3b) for $t\in\mathcal{T_{-}}$, we arrive at

and thus inserting Eq. (16) we obtain

These approximations fully characterize the long-time quasi-stationary state in the fast-switching regime.

In order to test this approximate solution, we computed the root mean-square error between our ansatz and the result of numerically integrating the mean-field equations. This error was then divided by the actual density average as obtained from numerical integration to obtain a dimensionless error measure $\bar{\sigma}$. In Fig. 7, this relative error $\bar{\sigma}$ is plotted against the dimensionless carrying capacity period $T_{k}\sqrt{\lambda\alpha}$. As expected, our approximation yields small relative errors for $T_{k}\sqrt{\lambda\alpha}\ll 1$. Interestingly, the asymptotic expansion seems to work even up to values $T_{k}\sqrt{\lambda\alpha}=2$ with relative errors less than $10\%$.

In this section, we employ a lattice model to numerically simulate the stochastic Lotka–Volterra predator-prey system (1), which allows us to investigate spatial structures and reaction-induced spatio-temporal correlations. Utilizing a stochastic lattice model allows us to investigate resonance effects on correlations, and their relation to the intrinsic spatio-temporal patterns of our system. Direct comparison with the mean-field rate equation approximation delineates the latter’s validity range, thus providing information when stochastic fluctuations and correlation effects may be ignored without losing pertinent qualitative features. The stochasticity of the system is implemented through an individual-based Monte Carlo algorithm. We implement the model on a two-dimensional square lattice with periodic boundary conditions (i.e., a toroidal simulation domain), where each lattice site holds information about the number of individuals of each species at that location. The initial configuration of the system is set up as a disordered state where each individual is placed at a randomly selected lattice site. We employ the following notation: $n_{a}(x,y;t)$: number of predator individuals at site $(x,y)$ and at time $t$; $n_{b}(x,y;t)$: number of prey individuals at site $(x,y)$ and at time $t$; $N_{a}(t)$: total number of predator individuals across the entire lattice at time $t$; $N_{b}(t)$: total number of prey individuals across the entire lattice at time $t$; $\bar{n}_{i}(t)$: $N_{i}(t)/L^{2}$, where $L$ is the linear lattice size, denotes the average species $i\in(a,b)$ density. Time is simulated via Monte Carlo steps (MCS), such that at each Monte Carlo time step

1. 1.

   a random location on the lattice $(x,y)$ is picked;

2. 2.

   a random neighboring site is selected from the von-Neumann neighborhood (four nearest neighbors) $(x_{\rm new},y_{\rm new})$;

3. 3.

   if $(x,y)$ contains a predator individual, we attempt $n_{b}(x_{\rm new},y_{\rm new};t)$ predation reactions as follows:

   • •

     generate a uniformly distributed random number $r$;

   • •

     if $r<\lambda$, decrease the number of prey at $(x_{\rm new},y_{\rm new})$ by $1$ and increase the number of predators at $(x_{\rm new},y_{\rm new})$ by $1$;

4. 4.

   next attempt a death reaction for the predator as described below:

   • •

     generate a uniformly distributed random number $r$.

   • •

     if $r<\mu$, decrease the number of predators at $(x,y)$ by $1$;

5. 5.

   if $(x,y)$ contains a prey individual, attempt a reproduction reaction as follows:

   • •

     generate a uniformly distributed random number $r$;

   • •

     if $r<\sigma$ and $n_{a}(x_{\rm new},y_{\rm new};t)\\ +n_{b}(x_{\rm new},y_{\rm new};t)<K$, increase the number of prey at site $(x_{\rm new},y_{\rm new})$ by $1$;

6. 6.

   if $(x,y)$ is empty ($n_{a}(x,y;t)+n_{b}(x,y;t)=0$), return to step 1;

7. 7.

   the above steps are repeated $N_{a}(t)+N_{b}(t)$ times.

This implementation ensures that at each Monte Carlo time step, on average, all individuals in the lattice attempt a reaction. We utilize random updates (i.e., picking new lattice sites at random) rather than systematic sequential updates (going over each lattice site in a specific sequence) in order to avoid introducing any bias in how reactions occur in the system.

A choice now has to be made in how to precisely manage the population after switching from the high to the low carrying capacity, because there will likely be an excess number of individuals at some lattice sites. We have considered two implementations to deal with this issue: In the first variant, we randomly removed any excess individuals to immediately reach the allowed low carrying capacity value $K_{-}$. While this implementation leads to interesting period-doubling behavior, we deemed it to be unrealistic. In the second implementation, we left the excess particles on site, but restricted further prey reproduction at lattice locations with more individuals than permitted. Therefore, we allow the system to intrinsically relax to a configuration without excess individuals, since eventually any superfluous predators would be forced to perish, and any excess prey would be devoured by predators. This intrinsic relaxation introduces a time scale set by the internal response time of the system, which is in turn determined by the reaction rates.

The stochastic lattice system was simulated over multiple runs, thus averaging both over ensembles of different initial conditions and distinct temporal histories; $\langle\ldots\rangle$ denotes the resulting (double) ensemble averages. We measured the average spatial species densities $\rho_{i}(t)=\langle\bar{n}_{i}(t)\rangle$, and computed the (connected) auto-correlation functions at fixed positions, $C_{ij}(t,t_{0})=\langle\bar{n}_{i}(t)\bar{n}_{j}(t_{0})\rangle-\langle\bar{n}_{i}(t)\rangle\langle\bar{n}_{j}(t_{0})\rangle$. The static correlations as functions of spatial distance $|x-x_{0}|$ were extracted using the definition $C_{ij}(x,x_{0};y_{0},t_{0})=\langle n_{i}(x,y_{0},t_{0})n_{j}(x_{0},y_{0},t_{0})\rangle-\langle n_{i}(x,y_{0},t_{0})\rangle\langle n_{j}(x_{0},y_{0},t_{0})\rangle$. In the long-time regime, the system should be isotropic at length scales large compared with the lattice constant, so that static correlations along the $x$ or $y$ directions will become identical. We also assume that the system is homogeneous at those scales, and hence that the auto-correlations should be independent of the reference positions $(x_{0},y_{0})$. Consequently we determine the auto-correlations using the densities averaged over lattice sites, which improves our statistics. Both these assumption were confirmed via explicit simulations. Furthermore, we evaluated the static correlations at a specific, sufficiently late fixed time step $t_{0}$, but again checked that all correlations are invariant under discrete time translation $t_{0}\to t_{0}+T_{k}$ with the environmental switching period $T_{k}$.

The various parameters of the system are the three reaction rates $(\sigma,\mu,\lambda)$, the low and high carrying capacity values $(K_{-},K_{+})$, and the period of the oscillating environment $T_{k}$. However, we can eliminate one of these parameters by rescaling the units of time. In our Monte Carlo simulations, which are not intended to match any specific experimental or observational data, we chose to always fix $\sigma=\mu$, because these parameters represent the rates for the linear prey reproduction and predator death reactions, and we are predominantly interested in the behavior of the system as the non-linear coupling $\lambda$ is varied. Thus our varying control parameters consist of the set $(\lambda,K_{-},K_{+},T_{k})$. For fixed $\sigma$ and $\mu$, the critical predation rate $\lambda_{c}$ only depends on the carrying capacity. Therefore, for the remainder of this paper we shall implicitly assume a fixed value for $\sigma=\mu$ and indicate the critical threshold as $\lambda_{c}(K)$.

Snapshots of a single simulation run of a representative system at different times are depicted in Fig. 8. According to our setup, the system is initially in a random configuration, so the predators consume the prey available in their neighborhood. At $t=46$, the predators have devoured most of their prey, and their number decays over time: The system is in the predator extinction phase for $K=1$. Therefore, without the external periodic environmental variation, the predators would eventually go extinct. However, we see that at $t=70$, after the carrying capacity has jumped at $t=T_{k}/2=50$ from $K=1$ to $K=10$, the prey are permitted to reproduce more abundantly. The prey population increase induces spreading waves of predators, in turn causing an enhancement of the predator density, until by $t=88$ the latter almost fill the entire lattice. When the carrying capacity drops back to $K=1$ at $T_{k}=100$, the predator density starts decaying again over time towards the point of extinction until the carrying capacity is once more reset and the whole process (stochastically) repeats.

Figure 9 shows the long-time behavior of the density for two different values of the predation rate $\lambda$. For $\lambda=0.1$ (a), the system oscillates between the predator extinction phase, approached when $K=1$, and the two-species coexistence phase, when $K=10$. In contrast, for $\lambda=0.275$ (b) the system resides in the species coexistence phase at both $K$ values. Both population time traces show stable oscillations with the switching period $T_{k}$, as expected for a dynamical system driven by a periodic external force. This is further confirmed by the Fourier transform plots displayed in Fig. 10. The prey density becomes non-smooth at points where the carrying capacity switches from $K_{+}$ to $K_{-}$ or vice versa, while the predator density remains smooth at those points. This is indicative of the fact that only the prey density explicitly depends on the carrying capacity, while the predator density depends on $K$ through its coupling to the prey species. In Fig. 8a, we see that even though the system is in the predator extinction phase when $K=1$, the $A$ species are still able to maintain a non-zero population density through the periodic environmental variation. Indeed, we observe that the key difference between the runs for $\lambda=0.1$ and $\lambda=0.275$ resides in the amplitude of the oscillations, which drops significantly when the predation rate increases. This is a general feature of the static Lotka–Volterra model. However, the amplitude of the oscillation in Fig. 8a is even higher than would be attained in a static system with fixed $K=10$: Driving the system away from reaching the absorbing state causes the densities to overshoot their stationary state values for $K=10$. While a static system would go extinct for low values of the predation rate, the periodic temporal variation of the carrying capacity allows both species to coexist in this situation.

In fact, the population oscillations become most prominent if the carrying capacity effectively switches the system between the predators’ absorbing and active phases. To demonstrate that this is a generic feature of our model, we plot the maximum density values reached in the simulations in the long-time limit in Fig. 11. For $\lambda>6\sigma$, we observe predator extinction; as was noted in Ref. [13], the system may, depending on the initial conditions, evolve into one of the two absorbing states for large predation rates. We interpret this extinction transition to be caused by stochastic fluctuations in our finite simulation system: As the predation rate becomes large, stochastic fluctuations are increasingly likely to drive the simulation towards the absorbing predator extinction state. For smaller predation rates, the asymptotic predator density decreases with growing $\lambda$. In the two-species coexistence region, the simulation results for the systems with periodically varying environment exhibit markedly larger oscillation amplitudes for both predatator and prey populations. This enhancement of the maximum population density in a periodically varying environment relative to the static case is responsible for sustaining species coexistence in an extended region of parameter space. Moreover, the extinction transition at high predation rate is moved to larger values of $\lambda/\sigma$ for the simulation runs with periodically varying carrying capacities compared to systems with fixed environment.

The predator-prey density phase space plots are constructed in Fig. 12 by simulating the system for multiple predation rate values. We see that for each $\lambda$ the system fluctuates around a closed orbit. Upon increasing the predation rate $\lambda$, the radius of this closed orbit becomes smaller, while the influence of stochastic fluctuations become more apparent. For $\lambda=0.8$, the orbit approaches $\rho_{A}=0$ which means that the predator population is close to extinction. Raising the predation rate further to $\lambda=0.9$, the system reaches the (finite system size) absorbing state with vanishing predator density, see Fig. 11.

We next carefully investigate how the system behaves in the two opposite limits of fast and slow environmental switching, relative to the intrinsic period of the Lotka–Volterra population oscillations. Figures 13(a) and (b) show the both populations’ densities for $T_{k}=10$ (fast switching) and $T_{k}=460$ (slow switching). The time-averaged behavior of the density in the fast-switching regime resembles a system with a constant effective equivalent carrying capacity $K^{*}$ that should be related to $K_{-}$ and $K_{+}$. In the slow-switching regime the system is given sufficient time to approach a (quasi-)stationary state when $K_{+}=10$. The prey density then reaches very high values, and the system is slowly driven to predator extinction; however, it would take many cycles of the changing environment for this absorbing state to be attained. As the switching period $T_{k}$ is increased, the predator population may only survive for a few cycles; eventually, when $T_{k}$ is set too large, it will go extinct before the prey food resources become abundant again.

We now explore the equivalent static environment hypothesis in the fast-switching regime in more detail. The mean-field rate equations suggest that for very short periods this equivalent carrying capacity equals the harmonic average ${\bar{K}}$ of $K_{+}$ and $K_{-}$, since Eqs. (3) only explicitly depend on $1/K$. For longer periods, the mean-field model predicts that the dynamics becomes effectively equivalent to a quasi-static system with a rate-dependent equivalent carrying capacity $K^{*}$, Eq. (9). One should expect the slow-switching equivalent carrying capacity in the stochastic lattice model to display a similar dependence on the microscopic reaction probabilities. As mentioned earlier, their precise relationship with macroscopic reaction rates such as $\sigma$ and $\lambda$ is however subtle and difficult to capture quantitatively, which poses a problem for stringently testing Eq. (9) for the stochastic lattice model. Yet for large $K_{-}$ and $K_{+}$, $K^{*}$ reduces approximately to the harmonic average ${\bar{K}}$, independent of the reaction rates. Hence we focus on testing the equivalent static environment hypothesis mainly with this effective carrying capacity.

To this end, we first present Monte Carlo simulation data for our system with $K_{-}=2$ and $K_{+}=6$, hence ${\bar{K}}=3$, obtained for a series of different switching periods $T_{k}$, measured relative to the intrinsic population oscillation period $\tau$ at fixed $\bar{K}$. For comparison, we also display simulations with fixed carrying capacity $K=3$, and display the resulting population densities in Fig. 14. We find that the predator density in the oscillating environment does not behave as if the environment were static with a harmonically averaged carrying capacity ${\bar{K}}$, with a discrepancy in the predator density of at least $18.4\%$. In contrast, the time-averaged prey density $\rho_{\infty}$ matches with the static equivalent ${\bar{K}}$ value for $T_{k}\approx 2.2\tau$. Yet for faster switching rates, we observe worse agreement with a discrepancy of up to $5.45\%$. For periods $T_{k}>2.2\tau$, $\rho_{\infty}$ increases monotonically with $T_{k}/\tau$, deviating further from the average prey density for the static equivalent ${\bar{K}}$. For larger $T_{k}/\tau$, the discrepancy between the harmonically averaged and the oscillating environments becomes more enhanced, although deviations remain less than $10\%$. Hence we conclude that our prey density data for an oscillating environment can be satisfactorily described by an equivalent constant environment for a wide range of oscillation periods. We note that both time-averaged population densities exhibit resonance-like extrema at $T_{k}\approx\tau$, owing to the environment switching just after the predators and prey have reached their maximum and minimum population counts, respectively, following their intrinsic Lotka–Volterra oscillations. As the period of the environment increases, more of these population oscillations may occur before the carrying capacity is reset, and integrating over one cycle of the environmental switching effectively averages over multiple periods of the intrinsic oscillations.

In Fig. 15 we repeat this numerical investigation for $K_{-}=4$, $K_{+}=12$, thus ${\bar{K}}=6$. The time-averaged prey density $\rho_{\infty}$ for the oscillating environment agrees well with the corresponding value for the static equivalent environment for all switching periods. However, the predator density for low periods does not match the harmonic mean hypothesis. For periods $T_{k}>\tau$, we see that the predator density with the oscillating environment approaches $\rho_{\infty}$ for the static equivalent environment. This suggests that the harmonically averaged carrying capacity works well to describe the mean predator population density for large $K_{-}$ and $K_{+}$ values, and for large environment oscillation periods, such that the system reaches the stationary state before switching occurs. In conclusion, stochastic fluctuations may change the form of the general equivalent static carrying capacity (9), yet it can still be approximated by the harmonic average for large carrying capacities.

Our simulation results indicate that the functional dependence of the prey density on the carrying capacity can be well approximated as $b\sim 1/K$ for a large range of environmental switching periods. However, the predator density exhibits a more complicated dependence on the carrying capacity values and $T_{k}$; it can only be approximated by $a\sim 1/K$ for large $K_{-}$ and $K_{+}$ and for $T_{k}\gg\tau$. In the latter limit, the system reaches its quasi-stationary state before the environment switches, which for the used parameter values corresponds to a stable node with non-oscillatory kinetics; consequently, there is little variation with $T_{k}$. Generally we observe that the long-time behavior of both population densities depends on the carrying capacity period in a non-monotonic manner.

The predator-prey pursuit and evasion waves characteristic of the stochastic spatial Lotka–Volterra model are more prominent in systems with high reaction rates. Therefore, we study the ensuing correlations for $\sigma=\mu=0.5$, $\lambda=0.1$, and leave $K_{-}=1$, $K_{+}=10$. For these parameters the system resides deep in the predator extinction absorbing phase when $K=K_{-}$, and in the active two-species coexistence phase for $K=K_{+}$. The behavior of the system for environmental switching period $T_{k}=10$ is exemplified by the simulation snapshots depicted in Fig. 16. The predators are initially almost driven to extinction, but due to the switching environment the prey population increases until it fills most of the lattice. We observe that at $t=23$ there remain only a few surviving predators which become localized sources for spreading waves. At $t=28$, the prey may proliferate in the interior of the fronts as well, causing the population waves to spread both outwards and inwards, until they eventually collide and interfere with each other as seen at $t=50$. Starting from $t=101$, the lattice exhibits a global density oscillation, and it becomes difficult to discern the original locations of the wavefront sources.

The associated temporal auto- and static correlation functions are displayed in Fig. 17. The auto-correlation functions exhibit damped oscillations with a peak period $2T_{k}=20$, twice the switching period of the carrying capacity. This is due to the fact that the two-point correlation function contains a product of particle densities, and the square of sinoidal functions may be decomposed into sine functions with doubled period. Note that the auto-correlations decay to zero after approximately $40$ time steps. The on-site population restrictions induce anti-correlations between individuals of the same species; the cross-correlation function $C_{ab}$ becomes positive after some time has elapsed., indicating that surviving predators follow the prey with some time delay. The static correlation functions rapidly decay to zero, demonstrating that the spatial correlation lengths are small, on the scale of a few lattice spacings.

Figure 18 shows simulation snapshots for the system parameters, but with a larger switching period $T_{k}=30$. In this run, only one predator patch has survived by $t=39$. Subsequently it serves as a source for a spreading population wave that later interferes with itself owing to the periodic boundary conditions of the lattice. At $t=56$ the wave starts spreading in both directions until at $t=72$, when the system returns to the low carrying capacity regime, and the prey in the interior of the front are not allowed to reproduce further. Even after a long time period at $t=264$, there is only a single density oscillation center that is sourced by the sole predator patch that had survived at $t=39$.

In Fig. 19a we plot the corresponding auto-correlation functions, which exhibit a much slower decay compared to Fig. 17(a) for $T_{k}=10$. This suggests that a carrying capacity period of $T_{k}=30$ causes a resonance effect, which indeed becomes apparent in the simulation movies, as in this case the switching happens approximately when the waves travel back to the location of the source. The resonance sustains the spatial and temporal correlations and thereby stabilizes the travelling waves, leading to a sustained asynchrony that promotes species coexistence. The Fourier transform again confirms that the auto-correlation functions oscillate with a period $2T_{k}$. Since the carrying capacity switching period is $T_{k}=30$, and it is initialized with $K(t=0)=K_{-}$, the behavior of the system at different $t_{0}$ values can be described as follows: For $t_{0}=990$, the system has just switched from $K(t)=K_{+}$ to $K_{-}$; at $t_{0}=1000$, it still resides at carrying capacity $K_{-}$; for $t_{0}=1005$, the system has just switched from $K_{-}$ back to $K_{+}$; at $t_{0}=1015$, the carrying capacity is still $K_{+}$. The static correlation functions, shown in Figs. 19(b,c,d), exhibit similar behavior for $t_{0}=990$ and $t_{0}=1015$, and for $t_{0}=1000$ and $t_{0}=1005$, respectively, which suggests a common delay time for the correlations. At $t_{0}=990$, the system is in the state with $K(t)=K_{-}$, while at $t_{0}=1015$, $K(t)=K_{+}$, about to switch to $K_{-}$; and similarly at $t_{0}=1005$ and $t_{0}=1000$. Compared with the system with faster switching period $T_{k}=10$, the static correlations decay over a larger distance, in agreement with the movies and snapshots which show wider wavefronts. The predator-prey cross-correlation function $C_{ab}(x)$ displays maxima at positive values for $t_{0}=1000$ and $t_{0}=1005$, when the carrying capacity is low and few individuals are present per site. Conversely at $t_{0}=990$ and $t_{0}=1015$, when the population densities are large, the only positive peak occurs at $x=0$, due to the fact that predators tend to be on the same site as prey for large $K$. For low carrying capacities, the predators cannot reside on the same locations as the prey, so instead they are most likely to be in the close prey neighborhood.

Finally, we further investigate the properties of our system by applying an asymmetric square signal for the switching carrying capacity, such that $K=K_{-}$ for $T_{-}$ time steps, and then $K=K_{+}$ for the subsequent time interval of length $T_{+}$, where $T_{-}\neq T_{+}$. The total switching period of the carrying capacity then is $T_{k}=T_{-}+T_{+}$. Simulating such a stochastic lattice system reveals a period-doubling effect for an intermediate range of $T_{+}/T_{-}$ ratios, as shown in Fig 20. For either too small or too large time interval ratios, no period-doubling effect could be observed. The origin of this intriguing period-doubling effect appears to be that prey particles are not able to reproduce quickly enough while the system has attained the high carrying capacity $K_{+}$. Hence, it takes the system two cycles of the oscillating environment for the prey density to reach its peak value.