GLOBAL STABILITY IN CHEMOSTAT-TYPE PLANKTON

MODELS WITH DELAYED NUTRIENT RECYCLING

Xue-Zhong He School of Mathematics and Statistics The University of Sydney Sydney, NSW 2006 Australia

and

Shigui Ruan^{∗}

Department of Mathematics, Statistics and Computing Science Dalhousie University Halifax, Nova Scotia Canada B3H 3J5

Abstract. In this paper, we consider chemostat-type plankton models in which plankton feeds on a limiting nutrient and the nutrient is supplied at a constant rate and is partially recycled after the death of plankton by bacterial decomposition. We use a distributed delay to describe nutrient recycling and a discrete delay to model the planktonic growth response to nutrient uptake. When one or both delays occur, by constructing Liapunov functionals, we obtain some sufficient conditions for the global stability of the positive equilibrium, which can be regarded as estimates of the delays for persistence of global stability.

Key words: Plankton models – Nutrient recycling – Time delay – Liapunov functionals – Global stability

Typeset by AMS-TEX

∗Research supported by the Natural Sciences and Engineering Research Council of Canada.

1. Introduction

Chemostat models have been studied by many authors in simulating the growth of planktonic communities of unicellular algae (see Smith and Waltman [19] and the references cited therein). An important difference between a chemostat situation and a lake situation is that lakes generally have a residence time of nutrient and sediments measured in years (see Powell and Richerson [15]). This implies that in models of natural systems the washout rate is very small and the regeneration of nutrient due to bacterial decomposition of the dead biomass must be considered (see Svirezhev and Logofet [20]).

The effect of nutrient recycling on stability of closed ecosystems has been studied by Nisbet and Gurney [13], Nisbet, McKinstry and Gurney [14], and Ulanowicz [21], etc. Powell and Richerson [15] and Nisbet and Gurney [13] regarded nutrient recycling as an instantaneous term, thus neglecting the time required to regenerate nutrient from dead biomass by bacterial decomposition. However, a delay in nutrient recycling is always present in a natural system and it increases when temperature decreases (see Whittaker [22]). Based on this observation, Beretta, Bischi and Solomano [1] proposed an open system with a single species feeding on a limit nutrient which is partially recycled after the death of the organisms. They used a distributed delay to model nutrient recycling and studied its effect on the stability of the positive equilibrium. Bischi [5] considered the same model and discussed the effect of the delay on resilience.

Evidence of delayed growth response has been observed in chemostat experiments with microalgae Chilamidomonas Reinhardii even when the limiting nutrient is at undetectable small concentration (Caperon [6] and Cunningham and Nisbet [7]). In a previous paper [17], one of us introduced a discrete delay in the growth response of the species to nutrient uptake in the model of Beretta, Bischi and Solomano [1]. By using the discrete delay as a bifurcation parameter, it was shown that the model undergoes a Hopf bifurcation. It was also found that both the distributed delay and the discrete delay do not change the property of persistence. Recently, Beretta and Takeuchi [2] used another distributed delay to describe the delayed growth response, namely, they considered a system of two retarded functional differential equations with two distributed delays. By constructing Liapunov functionals, they investigated the global asymptotic stability of the positive equilibrium for the case with a Lotka-Volterra coupling between nutrient and the biotic species. In two other papers, they extended the results to the cases that the biotic species has selfregulation which accounts for a finite carrying capacity of the environment (Beretta and Takeuchi [3]) and that the growth response is described by the Michaelis-Menten law (Beretta and Takeuchi [4]).

In this paper, we consider the model proposed in Ruan [17], a system of two retarded functional differential equations with a distributed delay and a discrete delay, and its special cases. By using Liapunov functionals method, we obtain some sufficient conditions for the global asymptotic stability of the positive equilibrium. We show that global stability of the positive equilibrium persists if the delays are sufficiently small and we give explicit estimates of the delays. The paper is organized as follows. In Section 2, we first consider the model with continuously delayed nutrient recycling. In Section 3, we suppose that the delay describing nutrient recycling is a discrete (constant or variable) delay. The general case with a distributed delay and a discrete delay is considered in Section 4. Finally, a discussion is carried out in Section 5.

2. Continuously Delayed Nutrient Recycling

Consider the following integrodifferential model of plankton-nutrient interaction with a distributed time delay describing nutrient recycling

\begin{matrix}Ṅ = & D(N^{0} − N) − aU(N)P + bγ ∫_{0}^{∞} f(s)P (t − s) ds \\ Ṗ = & P [−(γ + D) + cU(N)], \end{matrix}  (2.1)

where N is the limiting nutrient concentration and P is a measure of plankton concentration at any time.

In system (2.1), all parameters are positive constants. N^{0} is the input concentration of the limiting nutrient, a is the maximum uptake rate of nutrient, c (≤ a) is the maximum specific growth rate of plankton, b (0 < b < 1) is the fraction of the nutrient recycled by bacterial decomposition of the dead plankton, γ is the death rate of plankton and D is the washout rate, so γ + D represents the total loss rate of the plankton.

The function U(N) describes the nutrient uptake rate of plankton. We assume the following general hypotheses on U(N) :

(i) U(N) is non-negative, increasing and vanishes when there is no nutrient; (ii) there is a saturation effect when the nutrient is very abundant.

That is, U(N) is a continuous function defined on [0, ∞) and

U(0) = 0, \frac{dU}{dN} > 0, lim_{N→∞} U(N) = 1.

These hypotheses are satisfied by the Michaelis-Menten function (Caperon [6])

U(N) = \frac{N}{k + N} ,

where k > 0 is the half-saturation constant or Michaelis-Menten constant.

The delay-kernel f(s) is a non-negative bounded function defined on [0, ∞) and describes the contribution of the plankton population dead in the past to the nutrient recycled at time t. The presence of the distributed time delay must not affect the equilibrium values, so we normalize the kernel such that

∫_{0}^{∞} f(s)ds = 1.  (2.2)

According to MacDonald [12], we define the average time lag as

T_{f} = ∫_{0}^{∞} sf(s)ds.  (2.3)

In particular, the exponential kernel f(s) = αe^{−αs}, α > 0 is called a weak kernel and is frequently used in biological modelling (see Cushing [8] and Cunningham and Nisbet [7]). For the weak kernel defined in (2.4), the average time lag is T_{f} = 1/α.

Note that E_{0} = (N^{0}, 0) is always an equilibrium for system (2.1), and if

γ + D < c and U^{−1}(\frac{γ + D}{c} ) < N^{0},  (2.4)

system (2.1) has a positive interior equilibrium E^{∗} = (N^{∗}, P^{∗}) with

N^{∗} = U^{−1}(\frac{γ + D}{c} ), P^{∗} = \frac{D(N^{0} − N^{∗})}{aU(N^{∗}) − bγ} .  (2.5)

Model (2.1) was proposed by Beretta, Bischi and Solomano [1]. They studied the local stability of the positive equilibrium. Bischi [5] considered the effect of delay on resilience, that is, the rate at which a system returns to a stable steady state following a perturbation. His results showed that when the system is characterized by oscillating behavior, an increase of the delay can have a stabilizing effect. The results in Ruan [17] indicated that if the kernel f(s) is a weak kernel, then system (2.1) is uniformly persistent, a term used to describe long term survival of the interacting species. When the functional response is Michaelis-Menten type, Beretta and Takeuchi [4] studied the global asymptotically stability of the positive equilibrium E^{∗}.

To study the global stability of E^{∗}, we define

x_{1} = N − N^{∗},  x_{2} = ℓn(P/P^{∗})  (2.6)

and

ξ(x_{1}) = U(x_{1} + N^{∗}) − U(N^{∗}),  (2.7)

such that −N^{∗} < x_{1} < +∞, −∞ < x_{2} < +∞, −U(N^{∗}) ≤ ξ(x_{1}) < 1 − U(N^{∗}). Then

N = x_{1} + N^{∗}, P = P^{∗}e^{x_{2}}

and x_{1}ξ(x_{1}) > 0 for any x_{1} ∈ (−∞, +∞) and x_{1}ξ(x_{1}) = 0 if and only if x_{1} = 0. By the variable change (2.6), system (2.1) can be written as follows:

\begin{matrix}ẋ_{1}(t) = & − Dx_{1}(t) − aP^{∗}e^{x_{2}(t)}ξ(x_{1}(t)) − aP^{∗}U(N^{∗})[e^{x_{2}(t)} − 1] \\ + bγP^{∗} ∫_{0}^{∞} f(s)[e^{x_{2}(t−s)} − 1] ds \\ ẋ_{2}(t) = & cξ(x_{1}(t)), \end{matrix}  (2.8)

in which now the positive equilibrium has been transformed to the trivial equilibrium x_{1} = x_{2} = 0.

Now define

V_{1}(x_{1}) = ∫_{0}^{x_{1}} ξ(u) du.  (2.9)

Then V_{1}(x_{1}) ≥ 0 and, along the solutions of (2.8),

V̇_{1}(t) = ξ(x_{1}(t))ẋ_{1}(t)

= −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t)) − aP^{∗}U(N^{∗})ξ(x_{1}(t))[e^{x_{2}(t)} − 1]

+ bγP^{∗}ξ(x_{1}(t)) ∫_{0}^{∞} f(s)[e^{x_{2}(t−s)} − 1] ds = −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t)) − aP^{∗}U(N^{∗})ξ(x_{1}(t))[e^{x_{2}(t)} − 1]

+ bγP^{∗}ξ(x_{1}(t))[ − ∫_{0}^{∞} f(s) ∫_{t−s}^{t} e^{x_{2}(u)}ẋ_{2}(u) du ds + [e^{x_{2}(t)} − 1]] = −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t)) − aP^{∗}U(N^{∗})ξ(x_{1}(t))[e^{x_{2}(t)} − 1]

+ bγP^{∗}ξ(x_{1}(t))[[e^{x_{2}(t)} − 1] − c ∫_{0}^{∞} f(s) ∫_{t−s}^{t} e^{x_{2}(u)}ξ(x_{1}(u)) du ds] ≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t)) − P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t)} − 1]

+ \frac{1}{2}bcγP^{∗} ∫_{0}^{∞} f(s) ∫_{t−s}^{t} e^{x_{2}(u)}[ξ^{2}(x_{1}(t)) + ξ^{2}(x_{1}(u))] du ds

= −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t)) − P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t)} − 1]

+ \frac{1}{2}bcγP^{∗}[( ∫_{0}^{∞} f(s) ∫_{t−s}^{t} e^{x_{2}(u)} du ds)ξ^{2}(x_{1}(t))

+ ∫_{0}^{∞} f(s) ∫_{t−s}^{t} e^{x_{2}(u)}ξ^{2}(x_{1}(u)) du ds].  (2.10)

Let x(t) = (x_{1}(t), x_{2}(t)) and

V_{2}(x(t)) = \frac{1}{2}bcγP^{∗} ∫_{0}^{∞} f(s) ∫_{t−s}^{t} ∫_{v}^{t} e^{x_{2}(u)}ξ^{2}(x_{1}(u)) du dv ds.  (2.11)

From (2.10) and (2.11), we have

\frac{d}{dt}(V_{1} + V_{2})(x(t))

≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t))

− P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t)} − 1]
+ \frac{1}{2}bcγP^{∗}[( ∫_{0}^{∞} f(s) ∫_{t−s}^{t} e^{x_{2}(u)} du ds) + T_{f}e^{x_{2}(t)}]ξ^{2}(x_{1}(t)).  (2.12)

On the other hand, from the second equation of (2.1) one can see that Ṗ (t) ≥ −(γ + D)P (t), which implies that

ℓn\frac{P (t)}{P (s)} ≥ −(γ + D)(t − s).

Thus, we have

P (s) ≤ e^{(γ+D)(t−s)}P (t)  (2.13)

and

P^{∗} ∫_{0}^{∞} f(s) ∫_{t−s}^{t} e^{x_{2}(u)} du ds
= ∫_{0}^{∞} f(s) ∫_{t−s}^{t} P (u) du ds
≤ ( ∫_{0}^{∞} f(s) ∫_{t−s}^{t} e^{(γ+D)(t−u)} du ds)P (t)
= \frac{1}{γ + D} [ ∫_{0}^{∞} f(s)(e^{(γ+D)s} − 1)ds]P (t).  (2.14)

By (2.12) and (2.14), we obtain

\frac{d}{dt}(V_{1} + V_{2})(x(t))

≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t))

− P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t)} − 1]
+ \frac{1}{2}bcγ[\frac{1}{γ + D} ∫_{0}^{∞} f(s)(e^{(γ+D)s} − 1) ds + T_{f}]P (t)ξ^{2}(x_{1}(t)).  (2.15)

Let

V_{3}(x_{2}) = ∫_{0}^{x_{2}}[e^{s} − 1] ds.  (2.16)

Then V_{3}(x) ≥ 0 for x ∈ (−∞, +∞) and, along the solutions of (2.8),

V̇_{3}(x_{2}(t)) = [e^{x_{2}(t)} − 1]ẋ_{2}(t) = cξ(x_{1}(t))[e^{x_{2}(t)} − 1].  (2.17)

One can choose α > 0 such that

αc = P^{∗}[aU(N^{∗}) − bγ].

Then, for a Liapunov functional V (x)(t) defined by

V (x)(t) = V_{1}(x_{1})(t) + V_{2}(x)(t) + αV_{3}(x_{2})(t),  (2.18)

we have from (2.18), (2.15), and (2.17) that

V̇ (x)(t) ≤ − Dx_{1}ξ(x_{1})
+ P ξ^{2}(x_{1})[ − a + \frac{1}{2}bcγ(\frac{1}{γ + D} ∫_{0}^{∞} f(s)[e^{(γ+D)s} − 1] ds + T_{f})]. (2.19)

Denote

T_{f}^{∗} = \frac{1}{γ + D} ∫_{0}^{∞} f(s)[e^{(γ+D)s} − 1] ds < ∞.  (2.20)

By applying the standard global stability theorem on the method of Liapunov functional (cf Theorem 2.5.1 in Kuang [11]), we can summarize the above analysis as the following theorem.

Theorem 2.1. Assume that (a) the inequalities in (2.4) hold;

(b) T_{f} < ∞, T_{f}^{∗} < ∞ and

γ[T_{f}^{∗} + T_{f}] < \frac{2a}{bc}.  (2.21)

Then the positive equilibrium E^{∗} of system (2.1) is globally asymptotically stable.

Remark 2.2 From the second equation of system (2.1), we also have

Ṗ ≤ P [c − (γ + D)].

Denote γ_{1} = c − (γ + D). Then γ_{1} > 0 and

P (s) ≤ P (t)e^{γ_{1}(s−t)}.

Using this in the derivation of (2.14), we can obtain a condition slightly different from (2.21), that is,

γ[T_{f}^{∗} + T_{f}] < \frac{2a}{bc},  (2.22)

where T_{f}^{∗} = \frac{1}{γ_{1}} ∫_{0}^{∞} f(s)[1 − e^{−γ_{1}s}] ds.

Remark 2.3 Global stability of some models related to system (2.1) has been studied by Beretta and Takeuchi [2, 3, 4]. In Section 4 of [2], they discussed the global stability of system (1.1) with a Lotka-Volterra coupling between nutrient and the plankton; in Section 4 of [3], they considered the global stability of system (1.1) with an additional negative term −dP in the second equation, i.e. they supposed that the planktonic species has self-regulation, but their results cannot be applied to the case when d = 0; in Section 4 of [4], they studied the global stability of system (1.1) with the response function being of the Michaelis-Menten form. Also, they derived global stability conditions for solutions with initial conditions in Q_{H} = {ϕ : ∥ϕ} < H} and depending on K, the upper bound of the solutions. In our case, such restrictions on initial conditions are not required.

3. Discretely Delayed Nutrient Recycling

3.1 Constant Delay. In this section, we suppose that the delay describing nutrient recycling is discrete. We first consider the following equations

\begin{matrix}Ṅ = & D(N^{o} − N) − aU(N)P + bγP (t − τ) \\ Ṗ = & P [−(γ + D) + cU(N)], \end{matrix}  (3.1)

where τ > 0 is a constatnt. Notice that if the inequalities in (2.4) hold, then the positive equilibrium E^{∗} given by (2.5) exits. By using the transformation (2.6), system (3.1) can be written as follows:

\begin{aligned}ẋ_{1}(t) & = −Dx_{1}(t) − aP^{∗}e^{x_{2}(t)}ξ(x_{1}(t)) − aP^{∗}U(N^{∗})[e^{x_{2}(t)} − 1] \\ & + bγP^{∗}[e^{x_{2}(t−τ)} − 1] \\ ẋ_{2}(t) & = cξ(x_{1}(t)). \end{aligned}  (3.2)

Define V_{1}(x_{1}) as in (2.9), then along the solutions of (3.2) we have

V̇_{1}(x_{1}(t)) = ξ(x_{1}(t))ẋ_{1}(t)

≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t)) − P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t)} − 1]

+ \frac{1}{2}bcγP^{∗}[( ∫_{t−τ}^{t} e^{x_{2}(u)} du)ξ^{2}(x_{1}(t))

+ ∫_{t−τ}^{t} e^{x_{2}(u)}ξ^{2}(x_{1}(u)) du].  (3.3)

Now define

W_{2}(x(t)) = \frac{1}{2}bcγP^{∗} ∫_{t−τ}^{t} ∫_{v}^{t} e^{x_{2}(u)}ξ^{2}(x_{1}(u)) du dv.  (3.4)

From (3.3) and (3.4), it follows that

\frac{d}{dt}(V_{1} + W_{2})(x(t))

≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t)) − P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t)} − 1]

+ \frac{1}{2}bcγP^{∗}[( ∫_{t−τ}^{t} e^{x_{2}(u)} du ds) + τe^{x_{2}(t)}]ξ^{2}(x_{1}(t)).  (3.5)

Using (2.13) and following the argument in (2.14), we obtain

P^{∗} ∫_{t−τ}^{t} e^{x_{2}(u)} du ≤ \frac{1}{γ + D} [e^{(γ+D)τ} − 1]P (t).  (3.6)

Hence, (3.5) becomes

\frac{d}{dt}(V_{1} + W_{2})(x(t))

≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t)) − P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t)} − 1]

+ \frac{1}{2}bcγ[\frac{1}{γ + D} (e^{(γ+D)τ} − 1) + τ]P (t)ξ^{2}(x_{1}(t)).  (3.7)

Let V_{3} be defined by (2.16). We can choose α such that αc = P^{∗}[aU(N^{∗}) − bγ]. Define a Liapunov functional as follows

W (x)(t) = V_{1}(x)(t) + W_{2}(x)(t) + αV_{3}(x)(t).  (3.8)

Then from (3.7), (3.8), and (2.17), we have

Ẇ (x)(t) ≤ − Dx_{1}ξ(x_{1})
+ P ξ^{2}(x_{1})[ − a + \frac{1}{2}bcγ(\frac{1}{γ + D} [e^{(γ+D)τ} − 1] + τ)].  (3.9)

Therefore, we have the following global stability theorem.

Theorem 3.1. Assume that (a) the inequalities in (2.4) hold;

(b) τ > 0 is finite such that

γ[\frac{1}{γ + D} [e^{(γ+D)τ} − 1] + τ] < \frac{2a}{bc}.  (3.10)

Then the positive equilibrium E^{∗} of (3.1) is globally asymptotically stable.

Remark 3.2 As in Reamrk 2.2, condition (3.10) in Theorem 3.1 can be replaced by

γ(\frac{[1 − e^{−γ_{1}τ}]}{γ_{1}} + τ) < \frac{2a}{bc}.

3.2. Variable delay. Consider the following equations with a variable delay

\begin{matrix}Ṅ = & D(N^{o} − N) − aU(N)P + bγP (t − τ(t)) \\ Ṗ = & P [−(γ + D) + cU(N)], \end{matrix}  (3.11)

where τ(t) is differentiable and satisfies

0 ≤ τ(t) ≤ T, 0 ≤ τ̇(t) < 1. (3.12)

Denote σ(t) = t − τ(t), then its inverse σ^{−1}(t) exists and denote q(t) = σ^{−1}(t) − t. Under the transformation (2.6), we can rewrite (3.11) as follows:

\begin{aligned}ẋ_{1}(t) & = −Dx_{1}(t) − aP^{∗}e^{x_{2}(t)}ξ(x_{1}(t)) − aP^{∗}U(N^{∗})[e^{x_{2}(t)} − 1] \\ & + bγP^{∗}[e^{x_{2}(t−τ(t))} − 1] \\ ẋ_{2}(t) & = cξ(x_{1}(t)). \end{aligned}  (3.13)

Define V_{1}(x_{1}) as in (2.9), then along the solutions of (3.13), we have

V̇_{1}(x_{1}(t)) = ξ(x_{1}(t))ẋ_{1}(t)

≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t)) − P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t)} − 1]

+ \frac{1}{2}bcγP^{∗}[( ∫_{t−τ(t)}^{t} e^{x_{2}(u)} du)ξ^{2}(x_{1}(t))

+ ∫_{t−τ(t)}^{t} e^{x_{2}(u)}ξ^{2}(x_{1}(u)) du].  (3.14)

Now define

W_{3}(x(t)) = \frac{1}{2}bcγP^{∗} ∫_{t}^{σ^{−1}(t)} ∫_{σ(s)}^{t} e^{x_{2}(u)}ξ^{2}(x_{1}(u)) du dv.  (3.15)

Then by (3.14) and (3.15), we obtain

\frac{d}{dt}(V_{1} + W_{3})(x(t))

≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t))

− P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t)} − 1]
+ \frac{1}{2}bcγP^{∗}[( ∫_{t−τ(t)}^{t} e^{x_{2}(u)} du ds) + q(t)e^{x_{2}(t)}]ξ^{2}(x_{1}(t)).  (3.16)

Using (2.13) and following the argument in (2.14), we have

P^{∗} ∫_{t−τ(t)}^{t} e^{x_{2}(u)} du ≤ \frac{1}{γ + D} [e^{(γ+D)τ(t)} − 1]P (t).  (3.17)

Hence, it follows from (3.16) that

\frac{d}{dt}(V_{1} + W_{2})(x(t))

≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t))

− P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t)} − 1]
+ \frac{1}{2}bcγ[\frac{1}{γ + D} (e^{(γ+D)τ(t)} − 1) + q(t)]P (t)ξ^{2}(x_{1}(t)).  (3.18)

Let V_{3} be defined by (2.16). One can select α such that αc = P^{∗}[aU(N^{∗}) − bγ]. Define

V (x)(t) = V_{1}(x)(t) + W_{3}(x)(t) + αV_{3}(x)(t)  (3.19)

as a Liapunove functional. Based on the above analysis, we have the following theorem:

Theorem 3.3. Assume that (a) the inequalities in (2.4) hold;

(b) the delay τ(t) satisfies (3.12) and

γ[\frac{1}{γ + D} [e^{(γ+D)τ(t)} − 1] + q(t)] < \frac{2a}{bc},  (3.20)

where q(t) = σ^{−1}(t) − t and σ(t) = t − τ(t). Then the positive equilibrium E^{∗} of (3.11) is globally asymptotically stable.

4. Delayed Growth Response

In this section, we suppose that there is a discrete delay in the growth response of plankton to nutrient uptake, that is, we consider the following system:

\begin{aligned}Ṅ & = D(N^{o} − N) − aU(N)P + bγ ∫_{0}^{∞} f(s)P (t − s) ds \\ Ṗ & = P [−(γ + D) + cU(N(t − τ))], \end{aligned}  (4.1)

where τ ≥ 0 is a constant. Notice that system (2.1) is a special case of system (4.1) when τ = 0. If the inequalities in (2.4) hold, then there is a positive equilibrium E^{∗} given by (2.5).

System (4.1) was considered in Ruan [17]. It was showed that the solutions are bounded, the discrete delay could distabilize the positive equilibrium E^{∗} and Hopf bifurcations may occur. In the following, we will give estimates of the delay such that stability of the positive equilibrium persists.

Let x(t) = (x_{1}(t), x_{2}(t)) be defined by (2.6) and ξ(x_{1}) be defined by (2.7). Then it follows from (4.1) that

\begin{aligned}ẋ_{1}(t) & = −Dx_{1}(t) − aP^{∗}e^{x_{2}(t)}ξ(x_{1}(t)) − aU(N^{∗})P^{∗}[e^{x_{2}(t)} − 1] \\ & + bγP^{∗} ∫_{0}^{∞} f(s)[e^{x_{2}(t−s)} − 1] ds \\ ẋ_{2}(t) & = cξ(x_{1}(t − τ)). \end{aligned}  (4.2)

Define

V_{1}(x_{1}(t)) = ∫_{0}^{x_{1}(t)} ξ(s) ds, V_{2}(x_{2}(t)) = ∫_{0}^{x_{2}(t)}[e^{s} − 1] ds.  (4.3)

Then along the solutions of system (4.2) we have

V̇_{1}(x_{1}(t)) = −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t))
− P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t)} − 1]
+ bγP^{∗}ξ(x_{1}(t)) ∫_{0}^{∞} f(s) ∫_{t−s}^{t} e^{x_{2}(u)}ẋ_{2}(u) du ds
= −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t))
− P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t+τ)} − 1]
+ P^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t))[e^{x_{2}(t+τ)} − e^{x_{2}(t)}]
+ bcγP^{∗}ξ(x_{1}(t)) ∫_{0}^{∞} f(s) ∫_{t−s}^{t} e^{x_{2}(u)}ξ(x_{1}(u − τ)) du ds  (4.4)

and

V̇_{2}(x_{2}(t + τ)) = cξ(x_{1}(t))[e^{x_{2}(t+τ)} − 1].  (4.5)

Choose α = P^{∗}[aU(N^{∗}) − bγ]/c and

V_{3}(x(t)) = V_{1}(x_{1}(t)) + αV_{2}(x_{2}(t + τ)).  (4.6)

Then it follows that

V̇_{3} = −Dx_{1}(t)ξ(x_{1}(t)) − aP^{∗}e^{x_{2}(t)}ξ^{2}(x_{1}(t))

+ cP^{∗}[aU(N^{∗}) − bγ]ξ(x_{1}(t)) ∫_{t}^{t+τ} e^{x_{2}(u)}ξ(x_{1}(u − τ)) du + bcγP^{∗}ξ(x_{1}(t)) ∫_{0}^{∞} f(s) ∫_{t−s}^{t} e^{x_{2}(u)}ξ(x_{1}(u − τ)) du ds

≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP (t)ξ^{2}(x_{1}(t))

+ \frac{cP^{∗}[aU(N^{∗}) − bγ]}{2} ∫_{t}^{t+τ} e^{x_{2}(u)}[ξ^{2}(x_{1}(t)) + ξ^{2}(x_{1}(u − τ))] du

+ \frac{bcγP^{∗}}{2} ∫_{0}^{∞} f(s) ∫_{t−s}^{t} e^{x_{2}(u)}[ξ^{2}(x_{1}(t)) + ξ^{2}(x_{1}(u − τ))] du ds.  (4.7)

Now define a Liapunov functional as follows

V (x(t)) = V_{3}(x(t)) + V_{4}(x(t)),  (4.8)

in which

V_{4}(x(t)) =\frac{1}{2}cP^{∗}[aU(N^{∗}) − bγ] ∫_{t−τ}^{t} ∫_{v}^{t} e^{x_{2}(u+τ)}ξ^{2}(x_{1}(u)) dv du

+ \frac{1}{2}bγc[ ∫_{0}^{∞} f(s) ∫_{t−s}^{t} ∫_{v}^{t} P (u)ξ^{2}(x_{1}(u − τ)) dv du ds

+ T_{f} ∫_{t−τ}^{t} P (s + τ)ξ^{2}(x_{1}(s)) ds].

Thus we have

\begin{aligned}V̇ & ≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP (t)ξ^{2}(x_{1}(t)) \\ & + \frac{1}{2}cP^{∗}[aU(N^{∗}) − bγ][ ∫_{t}^{t+τ} e^{x_{2}(u)} du + τe^{x_{2}(t+τ)}]ξ^{2}(x_{1}(t)) \end{aligned}

+ \frac{1}{2}bcγ[( ∫_{0}^{∞} f(s) ∫_{t−s}^{t} P (u) du ds) + T_{f}P (t + τ)]ξ^{2}(x_{1}(t)).  (4.9)

On the other hand, from the second equation of system (4.1), we have

Ṗ (t) ≥ −(γ + D)P

and hence

P (s) ≤ P (t)e^{(γ+D)(t−s)}.  (4.10)

It follows from (4.10) that

P^{∗}[ ∫_{t}^{t+τ} e^{x_{2}(s)} ds + τe^{x_{2}(t+τ)}]
= ∫_{t}^{t+τ} P (s) ds + τP (t + τ)

≤ P (t)[ ∫_{t}^{t+τ} e^{(γ+D)(t−s)} ds + τe^{−(γ+D)τ} ]

= [\frac{1}{γ + D} (1 − e^{−(γ+D)τ} ) + τe^{−(γ+D)τ} ]P (t).  (4.11)

By using (2.14), (4.11) and Remark 3.2, we obtain

\begin{aligned}V̇ & ≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP (t)ξ^{2}(x_{1}(t)) \\ & + \{\frac{1}{2}c[aU(N^{∗}) − bγ][\frac{1}{γ + D} (1 − e^{−(γ+D)τ} ) + τe^{−(γ+D)τ} ] \end{aligned}

+ \frac{bcγ}{2} [T̃_{f} + T_{f}e^{−(γ+D)τ}]\}P (t)ξ^{2}(x_{1}(t)),  (4.12)

where

T̃_{f} = min\{\frac{1}{γ + D} ∫_{0}^{∞} f(s)[e^{(γ+D)s} − 1]ds, \frac{1}{γ_{1}} ∫_{0}^{∞} f(s)[1 − e^{−γ_{1}s}] ds\}

with γ_{1} = c − (γ + D) > 0.

We have the following global stability theorem for system (4.1).

Theorem 4.1. Assume that (a) the inequalities in (2.4) holds;

(b) the delay τ > 0 and the kernel f satisfy

T_{τ} + T_{τ,f} < a,  (4.13)

where

T_{τ} =\frac{c}{2}[aU(N^{∗}) − bγ][\frac{1}{γ + D} (1 − e^{−(γ+D)τ} ) + τe^{−(γ+D)τ} ],
T_{τ,f} =\frac{1}{2}bcγ[T̃_{f} + T_{f}e^{−(γ+D)τ}].

Then the positive equilibrium E^{∗} of (4.1) is globally asymptotically stable.

Remark 4.2 Obviously, when τ = 0, Theorem 4.1 reduces to Theorem 2.1 and Remark 2.2.

A special case of system (4.1) is that when the kernel f is a Delta function, that is, the nutrient recycling is instantaneous. The model now has the form

\begin{aligned}Ṅ & = D(N^{o} − N) − aU(N)P + bγP \\ Ṗ & = P [−(γ + D) + cU(N(t − τ))]. \end{aligned}  (4.14)

System (4.14) was briefly discussed in Ruan [17]. It was found that the delay has an effect on stability of the positive equilibrium. However, by Theorem 4.1, we have the following result on the global stability of the positive equilibrium of system (4.14).

Theorem 4.3. Assume that (a) the inequalities in (2.4) holds;

(b) the delay τ > 0 satisfies

T_{τ} < a.  (4.15)

Then the positive equilibrium E^{∗} of (4.14) is globally asymptotically stable.

Remark 4.4 If τ = 0, (4.15) is automatically satisfied. This is Theorem 3.2 of Beretta, Bischi and Solimano [1] or Theorem 2.7 of Ruan [16], that is, if the positive equilibrium exists, it is globally asymptotically stable. Theorem 4.3 indicates that if the delay τ is sufficient small, then the positive equilibrium E^{∗} is still globally stable.

5. Discussion

In this paper, we have discussed global stability of the positive equilibrium of planktonnutrient models with delayed nutrient recycling. We supposed that there is a delay in the growth response of plankton to nutrient uptake. We also considered the case where the delay describing nutrient recycling is a constant or variable discrete delay.

In the first model, system (2.1), we supposed that the limit nutrient is partially recycled after the death of plankton due to bacterial decomposition and used a distributed delay to model nutrient recycling. Recall that the sufficient conditions we obtained for global stability of the positive equilibrium E^{∗} are as follows:

γ + D < c,  U^{−1}(\frac{γ + D}{c} ) < N^{0}  (2.4)

and

γ[T_{f}^{∗} + T_{f}] < \frac{2a}{bc},  (2.21)

where T_{f} is the average time delay defined in (2.3) and T_{f}^{∗} is defined by (2.20) and is somehow similar to the average time delay. The first inequality in (2.4) indicates that c, the maximum specific growth rate of plankton must be greater than γ + D, the total loss rate of plankton and the second inequality in (2.4) means that there must be enough nutrient input concentration. Notice that the inequalities in (2.4) are required throughtout the paper to gurantee the existence of the positive equilibrium and when the delays are absent, they are exactly the global stability conditions. Condition (2.21) demonstrates that the positive equilibrium is still globally asymptotically stable if the average time delay involved in nutrient recycling is relatively small, in other words, if the process of recycling is relatively short.

The results in Section 3 show that if the delay describing nutrient recycling is a constant or variable discrete delay, one can still have global asymptotic stability of the positive equilibrium as long as the delay is sufficient small.

In Section 4, we introduced a discrete delay in the planktonic growth response term. It has been found (see Ruan [17]) that this delay has a distabilizing effect on stability of the positive equilibrium. In Theorem 4.2, we showed that if the inequalities in (2.4) hold and if

T_{τ} + T_{τ,f} < a,  (4.13)

then the positive equilibrium E^{∗} of (4.1) is globally asymptotically stable, where T_{τ} and T_{τ,f} involve the discrete delay τ and the distributed delay kernel f . Inequality (4.13) thus gives us an estimate for the delays: when the delays are sufficiently small such that (4.13) holds, then the equilibrium is globally stable; when the delays are increased and (4.13) is invalid, the equilibrium may lose its stability and Hopf bifurcation may occur

(see Ruan [17]).

The results in this paper indicate that in order to have global asymptotic stability,

delays must be relatively small. These results justify the common belief that “small delays

are negligible in some modelling processes as far as stabilities are concerned”.

It should be pointed out that our method can be used to analyze models with two or

more various delays appearing in the first equation. Also, it would be interesting to study

the global asymptotic stability of the chemostat-type competition models (see Freedman

and Xu [9]) or predator-prey models (see Ruan and Wolkowicz [18]). We leave this for

future study.

References

[1] E. Beretta, G. I. Bischi and F. Solimano, Stability in chemostat equations with delayed nutrient recycling, J. Math. Biol. 28(1990), 99-111

[2] E. Beretta and Y. Takeuchi, Qualitative properties of chemostat equations with time delays: boundedness, local and global asymptotic stability, Differential Equations and Dynamical Systems 2(1994), 19-40.

[3] E. Beretta and Y. Takeuchi, Qualitative properties of chemostat equations with time delays II, Differential Equations and Dynamical Systems 2(1994), 263-288.

[4] E. Beretta and Y. Takeuchi, Global stability for chemostat equations with delayed nutrient recycling, Nonlinear World, 1(1994), 291-306.

[5] G. I. Bischi, Effects of time lags on transient characteristics of a nutrient cycling model, Math. Biosci. 109(1992), 151-175.

[6] J. Caperon, Time lag in population growth response of isochrysis galbana to a variable nitrate environment, Ecology, 50(1969), 188-192.

[7] A. Cunningham and R. M. Nisbet, Time lag and co-operativity in the transient growth dynamics of microalge, J. Theor. Biol. 84 (1980), 189-203.

[8] J. M. Cushing, Integrodifferential Equations and Delay Models in Population Dynamics, Springer-Verlag, Heidelberg, 1977.

[9] H.I. Freedman and Y. Xu, Models of competition in the chemostat with instantaneous and delayed nutrient recycling, J. Math. Biol. 31 (1993), 513-527.

[10] K. Gopalsamy, Stability and Oscillations in Dealy Differential Equations of Population Dynamics, Kluwer Academic Publishers, Dordrecht, 1992.

[11] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics. Academic Press, New York, 1993.

[12] N. MacDonald, Time Lags in Biological Models, Springer-Verlag, Heidelberg, 1978.

[13] R. M. Nisbet and W. S. C. Gurney, Model of material cycling in a closed ecosystem, Nature 264 (1976), 633-635.

[14] R. M. Nisbet, J. McKinstry and W. S. C. Gurney, A strategic model of material cycling in a closed ecosystem, Math. Biosci. 64 (1983), 99-113.

[15] T. Powell and P. J. Richerson, Temporal variation, spatial heterogeneity and competition for resource in plankton system: a theoretical model, Am. Nat. 125(1985), 431-464.

[16] S. Ruan, Persistence and coexistence in zooplankton-phytoplankton-nutrient models with instantaneous nutrient recycling, J. Math. Biol. 31 (1993), 633-654.

[17] S. Ruan, The effect of delays on stability and persistence in plankton models, Nonlinear Analysis 24 (1995), 575-585.

[18] S. Ruan and G. Wolkowicz, Persistence in plankton models with delayed nutrient recycling, Canad. Appl. Math. Quart. 3 (1995), 219-235.

[19] H. L. Smith and P. Waltman, The Theory of the Chemostat, Cambridge University Press, 1994.

[20] Y. M. Svirezhev and D. O. Logofet, Stability of Biological Communities, MIR, Moscow, 1983.

[21] R. E. Ulanowicz, Mass and energy flow in closed ecosystems, J. Theor. Biol. 34 (1972), 239-253.

[22] R. H. Whittaker, Communities and Ecosystems, Macmillan, New York, 1975.