Stability in predator – prey models and discretization of a modified Volterra – Lotka model

Abstract

We consider n  2 populations of animals that are living in mutual predator – prey relations or are pairwise neutral to each other. We assume the temporal development of the population densities to be described by a system of differential equations which has an equilibrium state solution. We derive sufficient conditions for this equilibrium state to be stable by Lyapunov's method. The results supplement those published elsewhere.

Further we consider a modification of the Volterra – Lotka model which admits an asymptotically stable steady state solution. This model is discretized in two ways and we investigate how small the time step size has to be chosen in order to guarantee that the steady state solution is an attractive fixed point of the discretized model. This investigation is connected with the determination of the model parameters from given data.

Keywords:

• Predator–prey model
• Equilibrium state
• Stability
• Discretization

1. Introduction

The oldest mathematical model to describe the predator – prey relation of two populations in a closed habitat is that of Volterra and Lotka. It consists of two differential equations of the form

where x(t) and y(t) denote the density of the predator and prey population at the time t, respectively. The coefficients satisfy the conditions

and are chosen such that the predator population decreases exponentially in the absence of prey and the prey population increases exponentially in the absence of predators. The conditions c ₁₂ > 0 and c ₂₁ < 0 determine that the encountering of predators and prey has a positive and negative effect on the speed of growth of the corresponding population, respectively.

Obviously, the above system has a steady state solution which is given by

This solution can be shown to be stable in the sense of Lyapunov (see, for instance 1).

In Section 2 we consider a generalization of the Volterra – Lotka model to the case of n  2 populations living in mutual predator – prey relations or being pairwise neutral to each other. This model is described by the system (2.1) and it is assumed that this system admits a steady state solution. In Section 2 we derive sufficient conditions for this steady state solution to be stable in the sense of Lyapunov. The results in this section supplement those in 2.

So far the predator – prey models are purely theoretical models.

The question now arises as to how these models can be adapted to measured data in order to find out whether these reflect predator – prey behaviour. Since the data are normally given at discrete points of time, it is reasonable to replace the time-continuous model by a time-discrete model that can be considered as an approximation of the time-continuous model. The simplest way of doing this is to replace the time derivatives [xdot](t) and [ydot](t) by difference quotients

with some given step size h > 0. In this way we obtain the following system of difference equations

This system has

as a unique fixed point solution. This fixed point, however, is repelling for every choice of h and hence not stable (see 1). So the time-discrete model is not a good approximation of the time-continuous model, which has the fixed point solution as a stable steady state solution.

In Section 3 we replace the Volterra – Lotka model by the system (3.1) where the additional coefficient c ₂₂ is non-positive. If c ₂₂ < 0 and the condition (3.3) is satisfied, then the system (3.1) has a steady state solution which is asymptotically stable (see 1) and the discrete system (3.4) has this solution as an attractive fixed point, if the step size is sufficiently small (also see 1). So the time-discrete system (3.4) can be considered as a good approximation of the system (3.1), if the time step size is sufficiently small. In addition, estimates for the step size are derived which guarantee that the fixed point of the time-discrete system (3.4) is attractive.

In Section 3 we also consider a second discretization of the system (3.1) which is given by (3.8) and also has the steady state solution of (3.1) as a fixed point which is attractive, if the time step size is sufficiently small. Again estimates for the step size are derived which guarantee the attractiveness of the fixed point of system (3.8).

In Section 4 we try to adapt the two systems (3.4) and (3.8) to realistic data which are taken from 3 with the aid of the least squares approximation method. It turns out that this is only possible for the second discretization of system (3.1) for which we can show that the sufficient conditions for the attractiveness of the fixed point are satisfied.

2 The model and sufficient conditions for stability of equilibrium states

In 2 we have investigated a general predator – prey model with respect to stability and asymptotic stability. In this model we consider n  2 populations X_i, i = 1, … , n, of animals that are living in mutual predator – prey relations or are pairwise neutral to each other. Let us denote by x_i (t) the density of population X_i at time t. We assume the temporal development of these densities to be described by the following system of differential equations

where

and, for i ≠ j,

We assume that there is an equilibrium state [xbar] = ([xbar] ₁,…,[xbar] _n ) of equation (2.1) that is a solution of

with

Then equation (2.1) can be rewritten in the form

Let us assume that, for every initial state x ⁰ ∊ IR ⁿ with for i = 1, … , n, there is exactly one solution x(t) = (x ₁(t), … , x_n (t)) of equation (2.3) with

For every such solution we define

and obtain from equation (2.3) the system

This system has (0, … , 0) ∊ IR ⁿ as its equilibrium state. Further it follows from the above assumption that, for every initial state u ⁰ ∊ IR ⁿ , there is exactly one solution u(t) = (u ₁(t), … , u_n (t)) of equation (2.4) for t  0 with

In this section we want to derive sufficient conditions for the equilibrium state (with (2.2)) of equation (2.1) to be stable, which is equivalent to (0, … , 0) ∊ IR ⁿ being a stable equilibrium state of equation (2.4).

For this purpose we select a non-empty subset J of (1, … , n) and assume that

Then we define a Lyapunov function by

and conclude

Further we obtain

and with

it follows that

if

By Satz 1.12 in 1 we therefore conclude that (0, … , 0) ∊IR ⁿ is a stable equilibrium state of equation (2.4), if the conditions (2.5) and (2.6) are satisfied.

Now we consider the following special case: For some m∊{1, … , n−1} let J = {m + 1, … , n}. Further let

The condition (2.5) is satisfied and condition (2.6) reads

The case where

in which the second condition of (2.7) is satisfied has been considered in 2.

We end with the case n = 3, m = 2 in which, instead of (2.8), we assume that

Then the conditions (2.7) are equivalent to

Equation (2.2) reads

The unique solutions are given by

As necessary and sufficient conditions for [xbar] ₁ > 0, [xbar] ₂ > 0, [xbar] ₃ > 0 we obtain

3 Discretization of a modified Volterra – Lotka model

We consider a predator – prey model which is described by the following system of differential equations

Here x(t) and y(t) denote the density of the predator and prey population at time t, respectively. The coefficients satisfy the conditions

First we consider the case c ₂₂ = 0 in which (3.1) are the differential equations of the well-known Volterra – Lotka model. In this case the system (3.1) has the steady state solution

In 1 we have shown that this solution is stable.

If c ₂₂ < 0, then (3.1) has the steady state solution

if

In 1 we have shown that this solution is asymptotically stable.

Now let h > 0 be a given step size. Then we discretize (3.1) by replacing [xdot](t) and [ydot](t) by the difference quotients

In this way we obtain the following system of difference equations:

This system has

as a unique fixed point solution, if condition (3.3) is satisfied.

In the case c ₂₂ = 0 (where (3.3) is satisfied) this fixed point, however, is repelling for every choice of h > 0 and hence not stable (see 1). If we replace the system (3.4) with c ₂₂ = 0 by

then this also has the unique fixed point solution (3.5) with c ₂₂ = 0 which is stable in the sense that the linearization of (3.6) around the fixed point has (0,0) as stable fixed point, if

(see 2).

In the case c ₂₂ < 0 we have shown in 1 that the fixed point solution (3.5) of the system (3.4) is an attractor, if the step size h > 0 is chosen sufficiently small. The question is how small h has to be chosen. In order to answer this question we consider the Jacobi matrix of the right-hand side of (3.4) at the fixed point (3.5) which is given by

As eigenvalues of J([xbar], [ybar]) we obtain

We distinguish three cases:

1. Then it follows that
   and −1 < λ₁ = λ₂, if and only if If this condition is satisfied, the fixed point (3.5) is an attractor.

2. Then it follows that
   if and only if
   If this condition is satisfied, the fixed point (3.5) is an attractor.

3. Then it follows that λ₂ < λ₁ < 1 and
   If this condition is satisfied, the fixed point (3.5) is an attractor.

If we replace the system (3.4) by

this system also has the fixed point solution (3.5), if the condition (3.3) is satisfied. The Jacobi matrix of the right-hand side of (3.8) at the fixed point (3.5) is given by

and has the eigenvalues

If we put

then a > b, b > 0 (if c ₂₂ < 0, which we assume) and we obtain

We again distinguish three cases:

1. (a + b)² = 2b.

   This implies b ² < (a + b)² = 2b, hence b < 2 or which is equivalent to

   Further it follows that λ_1,2 = 1 – a – b < 1 and λ₁ = λ₂ > – 1, if and only if 2 or b < 2. Therefore the fixed point (3.5) is an attractor, if the condition (3.9) is satisfied.

2. (a + b)² < 2b⇔1 – 2a > (1 – a – b)². Then it follows that |λ₁|² = |λ₂|² = 1 – 2a < 1. If we put

   then (a + b)² <2b turns out to be equivalent to which is only possible if , in which case the last inequality is equivalent to
   Therefore, if , then condition (3.10) is sufficient for the fixed point (3.5) to be an attractor. If , case (2) cannot occur.

3. (a + b)² > 2b. Then it follows that λ₂ < λ₁ <1 and −1 < λ₂, if and only if

   If condition (3.11) is satisfied, then the fixed point (3.5) is an attractor.

4 Determination of the model parameters from data

Let us assume that we are given data (x(i · h),y(i · h)) for some h > 0 and i = 0, … , N.

First we consider the system (3.4) with c ₂₂ < 0. On using c ₁ = −c ₁₂[ybar] the first equation can be rewritten in the form

The second equation is written in the form

In order to determine c ₁₂ we choose some [ybar], say

and minimize

In order to determine c ₂, c ₂₁, c ₂₂ we minimize

Let the data be given by the following table (see 3):

Table

ti	0	1	2	3	4	5	6	7	8	h = 1
x(ti )	0:300	0:148	0:144	0:248	0:171	0:284	0:243	0:191	0:166
y(ti )	0:688	0:851	0:852	0:707	0:788	0:624	0:679	0:761	0:781

If we choose [ybar] = 0.75 the minimization of (4.1) leads to c ₁₂ = 3.23. The minimization of (4.2), however, leads to c ₂ = −1.62, c ₂₁ = 3.6 and c ₂₂ = 1.16 which violate the condition (3.2).

If we replace (4.2) by

then we obtain by minimizing (4.3) the values

From (3.5) we obtain [xbar] = .21.

Further we have α = 0.9, β = 0.969 which implies and . Therefore condition (3.10) is satisfied which implies that (0.21, 0.75) is an attractive fixed point of (3.8).

The following table is also taken from 3:

Table

(ti )	0	1	2	3	4	5	6	7	8	h = 1
x(ti )	0:114	0:043	0:086	0:139	0:154	0:195	0:209	0:146	0:079
y(ti )	0:885	0:973	0:911	0:827	0:789	0:739	0:698	0:804	0:855

If we choose [ybar] = .831, we obtain c ₁₂ = 1.28 by minimizing (4.1). Minimizing (4.3) leads to

which gives [xbar] = 0.129.

Further we obtain α = 0.82269, β = 0.2205115 which implies

and

Therefore condition (3.11) is satisfied which implies that (0.129, 0.831) is an attractive fixed point of (3.8).

In the first example the condition (a + b)² <2b is satisfied which implies a ² <2b which in turn is equivalent to

Under this condition we have seen for the system (3.4) that |λ₁| = |λ₂| < 1, if and only if

In the first example we have

Hence the condition (4.4) is violated for h = 1.

This can be considered as an explanation for the fact that the discretization via (3.4) does not lead to acceptible parameters.

This can be supported by the following remark: We have seen above that the condition (a + b)² < 2b is equivalent to

Now let us assume that which is satisfied in the first example. Then it follows that

Therefore the condition (4.4) is stronger than the condition (4.5) for |λ₁| = |λ₂| < 1 to hold true, if .