Stochastic dynamics of predator-prey interactions.

The interaction between a consumer (such as, a predator or a parasitoid) and a resource (such as, a prey or a host) forms an integral motif in ecological food webs, and has been modeled since the early 20th century starting from the seminal work of Lotka and Volterra. While the Lotka-Volterra predator-prey model predicts a neutrally stable equilibrium with oscillating population densities, a density-dependent predator attack rate is known to stabilize the equilibrium. Here, we consider a stochastic formulation of the Lotka-Volterra model where the prey's reproduction rate is a random process, and the predator's attack rate depends on both the prey and predator population densities. Analysis shows that increasing the sensitivity of the attack rate to the prey density attenuates the magnitude of stochastic fluctuations in the population densities. In contrast, these fluctuations vary non-monotonically with the sensitivity of the attack rate to the predator density with an optimal level of sensitivity minimizing the magnitude of fluctuations. Interestingly, our systematic study of the predator-prey correlations reveals distinct signatures depending on the form of the density-dependent attack rate. In summary, stochastic dynamics of nonlinear Lotka-Volterra models can be harnessed to infer density-dependent mechanisms regulating predator-prey interactions. Moreover, these mechanisms can have contrasting consequences on population density fluctuations, with predator-dependent attack rates amplifying stochasticity, while prey-dependent attack rates countering to buffer fluctuations.

1 Introduction

Predator-prey dynamics has been traditionally studied using an ordinary differential equation framework starting from the seminal work of Lotka and Volterra over a century ago [1-7]. The classical Lotka-Volterra model captures the dynamics of a predator-prey system, where h(t) and p(t) are the average population densities (number of individuals per unit area) of the prey, and the predator at time t. Here r represents the prey’s growth rate and h(t) grows exponentially over time in the absence of the predator. Predators consume prey with a constant rate f that we refer to as the attack rate, and each attacked prey leads to a new predator. Finally, each predator dies at a rate γ. In addition to predator-prey systems, ecological examples of such dynamics include host-parasitoid interactions that have tremendous application in biological control of pest species [8-14]. In a typical interaction, parasitoid wasps search and attack their host insect species by laying an egg within the body of the host. The egg hatches into a juvenile parasitoid that develops within the host by eating it from the inside out. Once fully developed, the parasitoid emerges from the dead host to repeat the life cycle.

The steady-state prey and predator equilibrium densities corresponding to the Lotka-Volterra model (1) are given by respectively. It turns out that this equilibrium is neutrally stable resulting in cycling population densities with a period of (assuming perturbations around the equilibrium) [15], and such population cycles have fascinated theoretical ecologists with several interpretations/extensions [16-18]. There is a rich body of literature expanding the Lotka-Volterra model to understand how diverse processes can push the equilibrium towards stability or instability [19-28]. For example, self-limitation in the prey’s growth in the form of a carrying capacity stabilizes the equilibrium [15]. Interestingly, a wide class of two-dimensional predator-prey models with an unstable equilibrium results in a stable limit cycle [29, 30].

In this contribution, we focus on generalizing (1) to that considers a density-dependent attack rate f(h, p), where f is a continuously differentiable function in both arguments. A generalized attack rate encompasses a wide range of ecological mechanisms. At one end of the spectrum are prey-dependent attack rates that capture nonlinear functional responses. For a Type II functional response is a decreasing function of the prey density, where c1 > 0 is the attack rate at small prey densities and T is the handling time [31-33]. Basically, the total attack rate per predator f(h, p)h increases linearly with h at low prey densities, but saturates to 1/T at high prey densities. Similarly, a Type III functional response corresponds to a sigmoidal function that initially accelerates with increasing prey density and then saturates to 1/T. At the other end of the spectrum are predator-dependent attack rates. For example, a decreasing attack rate with increasing predator density implies mutual interference between predators [34-37], or aggregation of predators to a subpopulation of high-risk individuals [38-43]. In contrast, cooperation between predators is reflected in f increasing with predator density.

In this contribution we provide analytical conditions for having stable population dynamics in terms of the sensitivity of f to prey/predator densities, thus combining the impact of different mechanisms into a single generalized stability criterion. Furthermore, we consider a stochastic formulation of the model by allowing the prey’s growth rate to follow an Ornstein-Uhlenbeck random process that drives the deterministic predator-prey dynamics (3). While both demographic and environmental stochasticity have been previously incorporated in predator-prey models, they have primarily focused on characterizing the role of noise in driving population extinctions and facilitating coexistence [44-47]. It remains to be seen how prey- and predator-dependent attack rates impact population density fluctuations, and to address this we take a novel moment-based approach using the Linear Noise Approximation technique to derive closed-form formulas quantifying the extent of fluctuations. Systematic investigation of these formulas reveals how random fluctuations in the prey’s growth rate propagate to population densities and uncover mechanisms that amplify or attenuate these fluctuations. Moreover, our results show how simple statistical signatures (such as, the correlation between population densities) can inform density-dependent mechanisms at play in regulating population dynamics.

2 Stability analysis of the generalized Lotka-Volterra model

Setting the left-hand-side of (3) to zero, the equilibrium population densities h* and p* of the generalized Lotka-Volterra model are the solution to

We assume that the function form of f is such that (6) yields a unique non-trivial equilibrium. Before performing a local stability analysis around the equilibrium, we define two dimensionless log sensitivities where is the partial derivative of f with respect to h evaluated at the equilibrium. To be biology realistic, we assume that f(h, p)p is an increasing function of the predator density that constrains f > −1, i.e., the decrease in f with increasing p cannot be faster than 1/p.

log sensitivity of the attack ratef(h, p) to the prey density

log sensitivity of the attack ratef(h, p) to the predator density

Linearizing the right-hand-side of (3) around the equilibrium yields the following Jacobian matrix

Stability of the equilibrium requires both equilibrium of the A matrix to have strictly negative real parts, and such a matrix is referred to as a Hurwitz matrix [48, 49]. For a two-dimensional system, the equilibrium is asymptotically stable, if and only if, the determinant of the A matrix is positive and its trace is negative [48, 49]. This implies that the equilibrium obtained as the solution to (6) is asymptotically stable, if and only if, both these inequalities hold

The grey shaded region in Fig 1 shows the stability region as a function of the log sensitivities f and f, with the neutrally stable Lotka-Volterra equilibrium corresponding to f = f = 0 on the edge of stability. Moreover, the intersection of the two lines in Fig 1 reveals as a necessary condition for stability. These stability conditions for continuous-time predator-prey models are analogous counterparts to recently developed stability conditions for discrete-time predator-prey models [50, 51].

Fig 1

The grey shaded region represents the region of stability for the equilibrium of the generalized Lotka-Volterra model (3).

The stability criterion (10) is plotted in terms of the log sensitivities of the attack rate f(h, p) to the prey (f) and predator (f) densities. The yellow line within the stability region separates the regions of negative real eigenvalues of the Jacobian matrix and complex eigenvalues with negative real parts. For this plot, the prey’s growth rate is assumed to be r = 2 per unit time and γ = 1 per unit time.

It is clear from Fig 1 that as reported in previous analysis [52], a Type II functional response with f < 0 and f = 0 will lead to an unstable equilibrium. In contrast, f > 0 and f = 0 stabilizes the equilibrium. As discussed earlier, f > 0 arises in the initial phase of a Type III functional response, where the predator attack rate accelerates with increasing prey density. Interestingly, a Type II functional response (f < 0) can provide stability in a narrow range if combined with other mechanisms, such as mutual interference between predators where f < 0. Overall these results show that an attack rate that increases with prey density is sufficient to stabilize the equilibrium as long as −1 < f < r*f. Similarly, a predator-dependent attack rate with f = 0 and −1 < f < 0 is sufficient to stabilize the equilibrium. Finally, we point out that the line divides the stability region into two parts—negative real eigenvalues of the Jacobian matrix below the line, and complex eigenvalues with negative real parts above the line. With increasing r, the line shifts further to the left. This separation within the stability region is relevant in the stochastic formulation of the model, where a stable equilibrium with complex eigenvalues of the A matrix can yield signatures of oscillatory dynamics in the presence of noise.

3 Stochastic formulation of the generalized Lotka-Volterra model

Having determined the stability regions of the deterministic model (3), we next turn our attention to the stochastic formulation of this model. In this context, much prior work has studied demographic stochasticity arising at low population abundances using Lotka-Volterra and spatial predator-prey models [48, 53–59] with important consequences on optimal harvesting strategies of a fluctuating resource [60, 61]. Here, we focus on environmental stochasticity that arises through randomness in the prey’s growth rate. Towards that end, we let the prey’s growth rate r(t) evolve as per an Ornstein-Uhlenbeck (OU) process where w(t) is the Wiener process and r* is the mean level of r(t). By using an OU process we capture memory in growth-rate fluctuations, with parameters γ > 0 and σ > 0 characterizing the time-scale and magnitude of r(t) fluctuations, respectively. These growth-rate fluctuations in turn drive population-density fluctuations through the model

Before describing mathematical tools for quantifying statistical moments of population densities, we point out that our approach of incorporating environmental stochasticity is different to other works that have considered either seasonal deterministic variations in r(t) [48, 62] or have added memoryless Brownian noise terms to the deterministic population dynamics [44, 63, 64]. The main reason for choosing an OU process to emulate environmental stochasticity is that it has a certain timescale of fluctuations. In that sense, a given environment can persist for some time, and have an impact on the population dynamics. We believe this is a much more ecologically relevant way to introduce environmental noise in contrast to white noise that has a flat frequency spectrum, and this is indeed a novelty of the paper.

To obtain the time evolution of the statistical moments of r(t), h(t) and p(t) corresponding to the nonlinear stochastic dynamical system (11) and (12) we use the following result. For any continuously differentiable function ψ(r, h, p), its expected value evolves as and moment dynamics is obtained by simply using a monomial in (13) [65]. Throughout the manuscript we use 〈 〉 to denote the expected value operation. For example, taking m1 = 1 or 2 with m2 = m3 = 0 yields the time evolution of the first two statistical moments of r(t) that result in the following steady-state mean and variance respectively. Similarly, to derive the mean dynamics of the prey’s population density we use m2 = 1, m1 = m3 = 0 to obtain where the right-hand-side now consists of higher-order moments. This problem of unclosed moment dynamics, where the time evolution of lower-order moments depends on higher-order moments has been well described for nonlinear stochastic systems, and often arises in the modeling of biochemical and ecological processes [14, 66–82]. Typically, different closure schemes are employed to approximate moment dynamics and we use one such approach known as the Linear Noise Approximation (LNA) [83-87]. In essence, assuming a stable equilibrium (h*, p*) in the deterministic formulation as given by (replacing r by r* in (6)) then for small fluctuations in r(t), h(t), p(t) around their respective equilibriums, the model nonlinearities can be linearized as

Moment dynamics is then derived after replacing these linear approximations in place of their nonlinear terms in (13) resulting in a closed system—the time derivative of a second-order moment now only depends on moments of order up to two. More specifically, if we collect all the first and second-order moments within the vector then it’s time evolution is given by a linear time-invariant system for some and matrix A. Solving this linear system at steady-state quantifies the magnitude of fluctuations in the population densities.

4 Quantifying random fluctuations in population densities

In the previous section, we described a LNA-based approach to quantify the statistical moments of population densities. Here we present some of the key results and insights from the analysis of moments. We quantify noise in the random processes r(t), h(t), p(t) using the square of their respective coefficient of variations

Solving (21) in Wolfram Mathematica yields the following analytical expressions for the noise in the prey and the predator population densities (normalized by the noise in the prey’s growth rate) respectively, and they depend on four dimensionless parameters—the sensitivity of the attack rate to the prey density (f), the sensitivity of the attack rate to the predator density (f), the prey’s average growth rate and the time-scale of fluctuations in r(t) normalized by the predator’s death rate

Recall from (10) that the stability of the deterministic equilibrium constraints f and f in the stability region of Fig 1 which ensures positivity of noise levels. Moreover, as one gets closer to the Lotka-Volterra model (f → 0 and f → 0), the system approaches the stability boundary leading to ( and ) in the LNA framework of noise derivation.

For an attack rate that only depends on the prey’s density (f = 0), (23) reduces to and fluctuations in population densities monotonically decrease with increasing dependence of the attack rate on the prey’s density (Fig 2). A closer look at (23) reveals that while the noise in the predator’s abundance approached a non-zero limit illustrating noise propagation of growth-rate fluctuations to the predator density via the prey.

Fig 2

Noise in the fluctuations of population densities as determined by (25) plotted for increasing f for f = 0.

Realizations of the prey’s growth rate, population densities of the prey and the predator as obtained by simulating (11) and (12) using are shown for f = 0.5 (left) and f = 10 (right). Each realization is normalized to its initial value so that it starts at a value of one. Note that for large values of f fluctuations in the predator density remain pronounced even though fluctuations in the prey density are minimal. For this plot, , , .

Our analysis further shows that when , then . In contrast, when , then for small value of f, and beyond a critical value f. One observation from the equilibrium analysis in (18) is that for a prey-dependent attack rate, h* is independent of r* implying that if the prey’s growth rate is chosen from a static distribution then it will not create any fluctuations in the preys’ density. This can be seen by considering the OU process in the limit γ → 0, σ → 0 keeping in (16) fixed. This limit corresponds to the scenario where, for each stochastic realization of the OU process a random initial condition is chosen with mean r* and variance , and this value remains the same over time. In this limit,

However, as p* is linearly dependent on r*

While fluctuations in population densities monotonically decrease with increasing f, the impact of a predator-dependent attack rate is quite different with noise levels varying non-monotonically with f (Fig 3). This effect can be understood in terms of the stability region in Fig 1, where for a given f > 0 stability requires −1 − f < f < r* f, and increasing f on either side puts the system closer to the stability boundary amplifying random fluctuations. This results in a scenario where fluctuations in population densities are minimized at an intermediate value of f.

Fig 3

Noise in the fluctuations of prey density as determined by (23) plotted as a function of f.

Noise is minimized at an intermediate value of f and stochastic realizations of the prey and predator densities are shown for three different values of f. For this plot, , f = 0.5, σ = 1 and . Note from Fig 1 that as f increases from −1 to 0.75 (for a fixed value of f = 0.5), one goes from real to complex eigenvalues that is reflected in the pronounced oscillatory dynamics of population densities as seen in the simulations on the right.

We next investigate the predator-prey Pearson’s correlation coefficient

The moment dynamics (21) results in the following closed-form expression highlighting an interesting result—while a negative dependence of the attack rate (f < 0) on the predator density leads to positive predator-prey correlations, a positive dependence f > 0 leads to negative correlations (Fig 4). Moreover, predator-prey density fluctuations are predicted to be uncorrelated for a prey-dependent attack rate. To understand this result, one can derive from (18) the following log sensitivities of the equilibrium densities to r*

Fig 4

Pearson’s correlation coefficient between the predator and the prey population densities as predicted by (31) for varying levels of f with f < 0 (f > 0) driving a positive (negative) correlation.

Sample trajectory paths are shown for f = 1 and f = −1. Other parameters taken as , f = 0.6 and .

Thus, when f > 0, the prey’s equilibrium density decreases with increasing r*, while the predator’s equilibrium density always increases with r*. These opposing responses of equilibrium densities intuitively explain the negative correlation seen for f > 0. In contrast, when f < 0, both equilibrium densities increase with r* and manifest in a positive correlation in the stochastic model. Recent work in host-parasitoid discrete-time models with a random host reproduction rate has also identified contrasting correlations depending on the mechanism stabilizing the population dynamics [88].

5 Conclusion

In summary, we have developed a novel stability criterion for a generalized Lotka-Volterra model with a density-dependent attack rate. (Fig 1). These results reveal that a Type II functional response can stabilize the equilibrium if combined with mechanisms involving predator inefficiency that puts the system in the grey shaded region corresponding to f < 0 and f > 0. Moreover, stability arises quite robustly for a Type III functional response when f > 0 as long as f is small enough to be in between −1 − f and rf.

Our stochastic analysis exploiting the linear noise approximation results in novel analytical formulas for and , providing insights into the propagation of growth rate fluctuations via the nonlinear dynamics to impact population density fluctuations. It is important to point out that these LNA-derived formulas are exact only in the small-noise limit (i.e., and ), and become less reliable as noise levels increase close to the stability boundary. Through these formulas, one can grasp important qualitative trends, such as, density fluctuations monotonically decrease with increasing f for a prey-dependent attack rate (Fig 2). However, predator-dependent attack rates can amplify stochasticity as f gets close to the stability boundary on either side of the x-axis in Fig 1. The impact of these sensitivities on the magnitude of density fluctuations have direct consequences on population extinctions. Finally, we have shown that population density correlations may contain signatures on stabilizing mechanisms at play with no correlation in predator-prey densities implying a prey-dependent attack rate (Fig 4), a negative correlation implying predator cooperation (f > 0), and a positive correlation implying mutual interference between predators (f < 0). Future work will expand these results to consider demographic stochasticity by explicitly modeling probabilistic birth-death events. It will also be interesting to consider competition between two or more consumers, and also look at the apparent competition between different prey species attacked by a common predator [89, 90]. Along these lines, new results have recently been developed in the discrete-time framework [50, 91].