"Predator and prey" model revisited -- influence of external fluxes and noise

TL;DR

This paper revisits the predator-prey model by incorporating external fluxes and noise in kinetic coefficients, exploring how these modifications affect the dynamics of predator and prey populations.

Contribution

It introduces two natural modifications to the classical predator-prey model: external fluxes via licensing and stochastic noise in kinetic coefficients.

Findings

External fluxes influence population stability.

Noise in kinetic coefficients affects oscillation patterns.

Modified model provides more realistic ecological dynamics.

Abstract

In this paper, we suggest two ways of modifications, which seem natural: a) we will introduce sponsors/hunters of predators or/and preys with "license" for constant rate of sponsoring/hunting (regular external fluxes); b) we will introduce noise of kinetic coefficients.

Equations

$$\frac{d\bar{X}}{d\bar{t}}={\bar{k}}_1\bar{X}-{\bar{k}}_2\bar{X}\bar{Y}$$

$$\frac{d\bar{Y}}{d\bar{t}}={\bar{k}}_3\bar{X}\bar{Y}-{\bar{k}}_4\bar{Y}$$

$$\frac{1}{a}\frac{dX}{dt}=k_{1}X-k_{2}X\cdot Y$$

$$a\frac{dY}{dt}=k_{3}X\cdot Y-k_{4}Y$$

$$\frac{d\delta X}{dt}=a\cdot ((1-Y^{st})\delta X-X^{st}\cdot \delta Y)=0\cdot \delta X+(-a)\cdot \delta Y$$

$$\frac{dY}{dt}=\frac{1}{a}(Y^{st}\delta X+(X^{st}-1)\cdot \delta Y)=\frac{1}{a}\cdot \delta X+0\cdot \delta Y$$

$$\frac{dX}{dt}=X-X\cdot Y+b_x$$

$$\frac{dY}{dt}=X\cdot Y-Y+b_y$$

$$X^{st}=\frac{1-b_x-b_y}{2}+\sqrt{\frac{(1-b_x-b_y{)}^2}{4}+b_x},Y^{st}=\frac{1+b_x+b_y}{2}+$$

$$+\sqrt{\frac{(1-b_x-b_y{)}^2}{4}+b_x}$$

$$\frac{d\delta X}{dt}=(1-Y^{st})\delta X-X^{st}\cdot \delta Y$$

$$\frac{dY}{dt}=Y^{st}\delta X+(X^{st}-1)\cdot \delta Y$$

$$\det \Vert \begin{array}{cc}1-Y^{st}-\lambda & -X^{st}\\ Y^{st}& X^{st}-1-\lambda \end{array}\Vert =0\Rightarrow {\lambda }^2-(X^{st}-Y^{st})\lambda +X^{st}+$$

$$+Y^{st}-1=0,$$

$${\lambda }_1=\frac{1}{2}((X^{st}-Y^{st})+\sqrt{(X^{st}-Y^{st}{)}^2+4(1-X^{st}-Y^{st})})$$

$${\lambda }_2=\frac{1}{2}((X^{st}-Y^{st})-\sqrt{(X^{st}-Y^{st}{)}^2+4(1-X^{st}-Y^{st})})$$

$$k_1=1+\xi (t)$$

$$<\xi (t)\xi (t^{\prime })>=A^2\delta (t-t^{\prime })$$

$$k_1=1+\frac{A}{\sqrt{dt}}\sqrt{3}(2\cdot random-1),$$

$$k_1=1+\frac{A}{\sqrt{dt}}\sin (2\pi random)\sqrt{2\ln (1/random)}$$

$$\frac{dX}{dt}=(1+\xi (t))X-X\cdot Y$$

$$\frac{dY}{dt}=X\cdot Y-Y$$

$$<(\Delta X{)}^2+(\Delta Y{)}^2>\approx \alpha A^{2}t,\alpha \approx 1$$

$$\frac{dX}{dt}=(1+\xi (t))X-X\cdot Y,\frac{dY}{dt}=X\cdot Y-Y,<\xi (t)\xi (t^{\prime })>=$$

$$=A^2\delta (t-t^{\prime }),$$

$$<(\delta X{)}^2>+<(\delta Y{)}^2>=1\cdot A^{2}t$$

$$\frac{d\delta X}{dt}=\xi (t)-\delta Y$$

$$\frac{d\delta Y}{dt}=\delta X$$

$$\delta X(t=0)=0,\delta Y(t=0)=0$$

$$\frac{d<(\delta X{)}^2>}{dt}=-2<\delta X\delta Y>+2<\delta X(t)\xi (t)>$$

$$\frac{d<(\delta Y{)}^2>}{dt}=2<\delta X\delta Y>$$

$$\frac{d}{dt}(<(\delta X{)}^2>+<(\delta Y{)}^2>)=2<\delta X(t)\xi (t)>$$

$$\frac{d}{dt}\left(\begin{array}{c}\delta X\\ \delta Y\end{array}\right)=\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\left(\begin{array}{c}\delta X\\ \delta Y\end{array}\right)+\left(\begin{array}{c}\xi \\ 0\end{array}\right)\equiv \hat{M}\left(\begin{array}{c}\delta X\\ \delta Y\end{array}\right)+\hat{f},$$

$$\delta X(t=0)=0,\delta Y(t=0)=0$$

$$\left(\begin{array}{c}\delta X(t)\\ \delta Y(t)\end{array}\right)={\int }_0^t\exp ((t-t^{\prime })\hat{M})\left(\begin{array}{c}\xi (t^{\prime })\\ 0\end{array}\right)d{t}^{\prime }$$

$$<\left(\begin{array}{c}\delta X(t)\\ \delta Y(t)\end{array}\right)\xi (t)>=\left(\begin{array}{c}<\delta X(t)\xi (t)>\\ <\delta Y(t)\xi (t)>\end{array}\right)={\int }_0^t\exp ((t-t^{\prime })\hat{M})\cdot$$

$$\cdot <\left(\begin{array}{c}\xi (t^{\prime })\\ 0\end{array}\right)\xi (t)>d{t}^{\prime }={\int }_0^t\exp ((t-t^{\prime })\hat{M})\left(\begin{array}{c}{A}^2\delta (t-t^{\prime })\\ 0\end{array}\right)d{t}^{\prime }=$$

Full text

**MSC 2010. **34A34, 34F05, 34D08

**Key words. **Lotka-Volterra model, nonlinear differential equations, stability analysis, noise

1 Introduction

Predator-prey model, introduced by Lotka and Volterra, is a basic synergetic model demonstrating oscillatory behavior of the nonlinear biological, chemical or economic systems [1, 2, 3]. It is governed by the simple set of two nonlinear equations taking into account the natural birthrate $\overline{k}_{1}$ of preys, natural death rate $\overline{k}_{4}$ of predators, as well as ‘‘collisions’’ of preys with predators, unlucky for preys (rate $\overline{k}_{2}$) and lucky for feeding the new generations of predators (rate $\bar{k}_{3}$ ):

[TABLE]

(X - number of preys, Y - number of predators, ). If the ‘‘CREATOR’’ of this ecosystem, choosing the initial numbers of both species,‘‘misses’’ the stationary numbers $\bar{X}^{st}=\bar{k}_{4}/\bar{k}_{3},\bar{Y}^{st}=\bar{k}_{1}/\bar{k}_{2}$, system demonstrates the oscillatory behavior, all oscillations proceeding around the mentioned stationary point. Transition to non-dimensional variables, $X=\bar{X}/\bar{X}^{st},Y=\bar{Y}/\bar{Y}^{st},t=\sqrt{\bar{k}_{1}\bar{k}_{4}}\bar{t}$ ,provides the set of two universal equations:

[TABLE]

with $a=\sqrt{\bar{k}_{1}/\bar{k}_{4}},{k}_{1}=1,{k}_{2}=1,{k}_{3}=1,{k}_{4}=1$. This model is very idealized and almost closed (excluding unlimited feed for preys). It is a main reason why the phase trajectories in the standard Lotka-Volterra model neither converge nor diverge, but just oscillate. Mathematically, it means that the linearized set of equations for deviations,

[TABLE]

provides the purely imaginary (with zero real part) roots of the characteristic equation: $\det\begin{Vmatrix}-\lambda&-a\\ 1/a&-\lambda\end{Vmatrix}=0\Rightarrow\lambda^{2}=-1\Rightarrow\lambda=\pm{i}$, which means oscillations without divergence or convergence.

Everywhere below we will limit ourselves to the particular, most symmetric case $a=1$.

There exist a lot of modifications and generalizations of predator-prey model [4, 5, 6, 7, 8], including non-homogeneity of the system and account of diffusion, several types of preys or/and predators, noise (fluctuations) of preys or/and predators numbers. In this paper, we suggest two more ways of modifications, which seem natural:

We will introduce sponsors/hunters of predators or/and preys with ‘‘license’’ for constant rate of sponsoring/hunting (regular external fluxes). 2. 2.

We will introduce noise of kinetic coefficients.

We will see that incorporation of adding or deleting predators or/and preys substantially broadens the spectrum of possible regimes: (I) system may remain eternally oscillating without convergence neither divergence, as in the classic LV-model, (II) system can be stable and converge to the steady-state limit, (III) system can be metastable, converging to steady-state from the initial positions in some critical vicinity of stationary solution and diverging from the positions outside this critical region, (IV) totally unstable, always diverging system.

We will also see that the noise of kinetic coefficients, in average, leads to divergence, but the time law for the growth of mean squared distance from the stationary solution is peculiar and resembles Brownian motion.

2 Influence of regular adding/hunting

The basic equations for LV-system (Lotka-Volterra system) with external fluxes are:

[TABLE]

We will start from some simple examples. Regular adding preys without predators being touched ($b_{x}>0,b_{y}=0$) stabilizes the system (Fig.1a), regular hunting preys without predators being touched ($b_{x}<0,b_{y}=0$) destabilizes the system (Fig.1b):

In cases of adding/hunting only predators without preys being touched, we have analogic situation (Fig.2), but adding of predators stabilizes the system only for $0<b_{y}<1$.

Now we will consider the general case. Equations for stationary points give the formal solution

[TABLE]

Of course, only positive solutions ($X^{st}>0,Y^{st}>0$) should be further considered. Moreover, these solutions should be at least locally stable. To find the criteria of local stability, we linearize Eq. (2) in the vicinity of stationary solution determined by Eq. (2):

[TABLE]

Then the local stability is determined by the real parts of the roots of characteristic equation

[TABLE]

[TABLE]

The result of steady-state analysis is summarized at Fig.3.

In general, we should distinguish globally stable and locally stable (metastable) stationary solutions. It will be done elsewhere.

Two examples are demonstrated at Fig. 4(a, b). Namely, we demonstrate the convergence region (green) at X-Y plane around the metastable stationary points corresponding to following cases: (a)$b_{x}=1.1,b_{y}=-1.09,X^{st}=1.65,Y^{st}=1.66$, (b)$b_{x}=-0.1,b_{y}=0.3,X^{st}=0.64,Y^{st}=0.84$

3 Noise of kinetic coefficients

Influence of noise on the behavior of the LV system is not a new problem. Yet, as far as we know, only noise of X and Y had been explored - random adding or deleting of small number of preys or predators[4, 7]. Additionally of this, in our paper we explore the random fluctuations of kinetic coefficients $k_{1},k_{2},k_{3},k_{4}$ (for example, fluctuation of birthrate due to rainy days).

We will start from the Langevine noise of reduced birthrate without memory and with fixed amplitude $A$:

[TABLE]

In case of numeric modeling, one should introduce this noise in such a way that the change of the time step should not change the impact of noise. Our suggestion was a step-wise probability distribution:

[TABLE]

(Mean square of random function $\sqrt{3}(2\cdot{random}-1)$ is equal to 1). Alternatively, one may use the Gaussian distribution

[TABLE]

.

In more details, the introduction of the noise of kinetic coefficients is discussed in [8, 9] for the case of atomic migration. We started from stationary point ($X=1,Y=1$) as an initial condition. Typical phase trajectory as a numeric solution of the set

[TABLE]

is shown at Fig.5.

Then, we took ensemble of $M=100$ LV-systems, originating at the same stationary point (1,1), and found the mean squared displacement from this point as a function of time. The results for different time-steps and the same amplitude are shown at Fig.6.

One can see that dependencies for different time-steps are close to each other and can be approximated as

[TABLE]

If the initial point differs from the stationary point, then the initial mean squared displacement at first decreases and then follows the same time law - see Fig.7.

Analogic dependencies were obtained by numeric simulations for the noise of $k_{2},k_{3},k_{4}$.

Most probably parabolic dependence is related to ‘‘zero stability’’ of LV system, so that we observe something like random walk in XY space.

We can prove Eq. (3.4) analytically and find exact value of $\alpha$ at least for small deviations from steady-state using linearized kinetic equations.

Theorem 1**.**

Linearization (first-order approximation) of the Lotka-Volterra model in the vicinity of steady-state with the Langevin noise of the birthrate of preys and without external fluxes,

[TABLE]

provides the following parabolic law for the sum of dispersions for the ensemble of LV-systems:

[TABLE]

Proof.

Consider the first order approximation for the deviations from steady-state: $X=1+\delta{X},Y=1+\delta{Y},|\delta{X}|<<1,|\delta{Y}|<<1,$

[TABLE]

We consider an ensemble of LV-systems with identical initial condition

[TABLE]

Let us multiply Eq. (3.7a) by $2\delta{X}$ and average over ensemble:

[TABLE]

Let us multiply Eq. (3.7b) by $2\delta{Y}$ and average over ensemble:

[TABLE]

Now we add Eqs (3.8a)+(3.8b). It gives

[TABLE]

To make the set (3.8a, 3.8b) self-consistent, one should find the value of $<\delta{X}(t)\xi(t)>$ .

For this, first of all, we find the formal solution of set (3.7a, 3.7b) which we reformulate in the matrix form:

[TABLE]

Solution of Eq. (3.10) is:

[TABLE]

Then we can find two mean values $<\delta{X}(t)\xi(t)>$ and $<\delta{Y}(t)\xi(t)>$.

[TABLE]

(Note that the factor $1/2$ appeared because the argument of delta function $\delta(t-t^{\prime})$ is equal zero not within the integration interval but instead at it’s boundary.) Thus,

[TABLE]

Theorem 1 is proved. ∎

Now we proceed to the noise of birthrate of predators.

Theorem 2**.**

Linearization of the Lotka-Volterra model with the Langevin noise of the birthrate of predators and without external fluxes,

[TABLE]

provides the parabolic law (3.6) for the sum of dispersions.

Proof.

Again, consider the first order approximation for the deviations from steady-state:$X=1+\delta{X},Y=1+\delta{Y},|\delta{X}|<<1,|\delta{Y}|<<1,$

[TABLE]

We again consider an ensemble of LV-systems with identical initial condition

[TABLE]

Let us multiply Eq. (3.15a) by $2\delta{X}$ and average over ensemble:

[TABLE]

Let us multiply Eq. (3.15b) by $2\delta{Y}$ and average over ensemble:

[TABLE]

Now we add Eqs (3.16a)+(3.16b). It gives

[TABLE]

Now, one should find the value of $<\delta{Y(t)\xi(t)>}$.

For this, first of all, we find the formal solution of set (3.15a, 3.15b) which we reformulate in the matrix form:

[TABLE]

Solution of Eq. (3.18) is:

[TABLE]

Then we can find two mean values$<\delta{X}(t)\xi(t)>$ and $<\delta Y(t)\xi(t)>$ .

[TABLE]

Thus,

[TABLE]

[TABLE]

Theorem is proved. ∎

Now we proceed to the noise of extinction rate of preys (caught by predators).

Theorem 3**.**

Linearization of the Lotka-Volterra model with the Langevin noise of extinction rate of preys and without external fluxes,

[TABLE]

provides the parabolic law (3.6) for the sum of dispersions.

Proof.

Again, consider the first order approximation for the deviations from steady-state: $X=1+\delta X,\,\,Y=1+\delta Y,\,\,\,\left|{\delta X}\right|<<1,\,\,\left|{\delta Y}\right|<<1,$

[TABLE]

Eqs. (3.23a, 3.23b) are identical with eq. (3.7a, 3.7b) except sign before noise term. All characteristics (probabilities and mean values) of Langevin noise are symmetrical in respect to sign. Therefore Theorem 3 directly follows from Theorem 1. ∎

Now we proceed to the noise of reproduction rate of predators due to eating of preys.

Theorem 4**.**

Linearization of the Lotka-Volterra model with the Langevin noise of reproduction rate of predators due to eating of preys and without external fluxes,

[TABLE]

provides the parabolic law (3.6) for the sum of dispersions.

Proof.

Again, consider the first order approximation for the deviations from steady-state: $X=1+\delta X,\,\,Y=1+\delta Y,\,\,\,\left|{\delta X}\right|<<1,\,\,\left|{\delta Y}\right|<<1\,,$

[TABLE]

Eqs. (3.25a, 3.25b) are identical with eq. (3.15a, 3.15b) except sign before noise term. All characteristics (probabilities and mean values) of Langevin noise are symmetrical in respect to sign. Therefore, Theorem 4 directly follows from Theorem 2. ∎

Now we proceed to the noise in the LV-model with external fluxes – adding noise to various terms in eqs.(2). We will start with the simplest cases – when the external flux is present but it’s mean value is zero.

Theorem 5**.**

Linearization of the Lotka-Volterra model with the Langevin noise of external flux of preys,

[TABLE]

provides the parabolic law (3.6) for the sum of dispersions.

Proof.

Again, consider the first order approximation for the deviations from steady-state: $X=1+\delta X,\,\,Y=1+\delta Y,\,\,\,\left|{\delta X}\right|<<1,\,\,\left|{\delta Y}\right|<<1\,,$

[TABLE]

Eqs.(3.27a, 3.27b) are identical with (3.7a, 3.7b). Therefore, Theorem 5 directly follows from Theorem 1. ∎

Theorem 6**.**

Linearization of the Lotka-Volterra model with the Langevin noise of external flux of predators,

[TABLE]

provides the parabolic law (3.6) for the sum of dispersions.

Proof.

Again, consider the first order approximation for the deviations from steady-state: $X=1+\delta X,\,\,Y=1+\delta Y,\,\,\,\left|{\delta X}\right|<<1,\,\,\left|{\delta Y}\right|<<1\,$

[TABLE]

Eqs.(3.29a, 3.29b) are identical with (3.25a, 3.25b). Therefore, Theorem 6 directly follows from Theorem 4. ∎

Now we proceed to the noise in LV-systems with non-zero mean external fluxes. Of course, in this case the steady-state reduced values are not equal to 1 anymore, but instead determined by eq. (2). In this paper we consider only one case - the Lottka-Volterra model with the Langevin noise of the birthrate

[TABLE]

Linearization for first-order deviations from steady-state ($X={X^{st}}+\delta X,\,\,Y={Y^{st}}+\delta Y,\,\,\,\left|{\delta X}\right|<<{X^{st}},\,\,\left|{\delta Y}\right|<<{Y^{st}}\,$), gives

[TABLE]

We again consider an ensemble with identical initial condition $\delta X\left({t=0}\right)=0,\,\,\delta Y\left({t=0}\right)=0$ Let us multiply eq. (3.31a) by $2\delta X$ and average over ensemble:

[TABLE]

Let us multiply eq. (3.31b) by $2\delta Y$ and average over ensemble:

[TABLE]

Let us add the product of eq. (3.31a) with $2\delta Y$ and the product of eq. (3.31b) with $2\delta X$ and average over ensemble:

[TABLE]

To make the set (3.32a, 3.32b, 3.32c) self-consistent, one should find the values of $<\delta X\left(t\right)\xi\left(t\right)>$, $<\delta Y\left(t\right)\xi\left(t\right)>$.

For this, first of all, we find the formal solution of set (3.7a, 3.7b) which we reformulate in the matrix form:

[TABLE]

Solution is:

[TABLE]

Then we can find two mean values $<\delta X\left(t\right)\xi\left(t\right)>$ and $<\delta Y\left(t\right)\xi\left(t\right)>$.

[TABLE]

Thus, the set of equations (3.32a, 3.32b, 3.32c) becomes self-consistent:

[TABLE]

The formal solution (in matrix form) of eq. (3.34) is:

[TABLE]

where

[TABLE]

If the steady-state solution satisfies the local stability criterion (${\mathop{\rm Im}\nolimits}{\lambda_{1}}<0,\,{\mathop{\rm Im}\nolimits}{\lambda_{2}}<0$ in eqs. (2.4)), then the set (3.34) should describe the competition between dynamical tendency of attraction to the metastable steady-state and stochastic tendency of migration from this steady state. It is physically evident that sooner or later the stochastic (noise) will bring the system beyond the limits of metastability (beyond the convergence region). It is also evident that the “Mean Time To Failure” (MTTF) of the LV-system should depend on the noise amplitude. This problem will be considered elsewhere. (See also [10].)

Conclusion

Regular hunting of preys without troubling predators destabilizes the LV-system. 2. 2.

Regular hunting of predators without troubling preys destabilizes the LV-system. 3. 3.

Regular adding of preys always stabilizes the LV-system. 4. 4.

Regular adding of predators stabilizes the LV-system only if adding rate less than some threshold. 5. 5.

Diagram of regimes for LV-system with regular external fluxes is described by Fig.3. 6. 6.

Simultaneous adding of preys and limited hunting of predators may leave the system metastable but only within some critical region of initial parameters. 7. 7.

In numeric modeling, to make the impact of the noise of kinetic parameters independent on the choice of time step, the random perturbation kinetic coefficients should be proportional to noise amplitude and inversely proportional to the square root of the time step. 8. 8.

Mean square distance from stationary point (at least in case of small deviations from steady-state) increases proportionally to time with proportionality coefficient proportional to the noise amplitude:$<(\Delta{X})^{2}+(\Delta{Y})^{2}>={A}^{2}t$ 9. 9.

Competition between noise and metastability under non-zero external fluxes (including MTTF) will be analyzed elsewhere.