GLOBAL STABILITY IN CHEMOSTAT-TYPE EQUATIONS WITH DISTRIBUTED DELAYS

Xue-Zhong He^{†}, Shigui Ruan^{‡} and Huaxing Xia^{§}

Abstract We consider a chemostat-type model in which a single species feeding on a limiting nutrient supplied at a constant rate. The model incorporates a general nutrient uptake function and two distributed (infinite) delays. The first delay models the fact that the nutrient is partially recycled after death of biomass by bacterial decomposition and the second delay indicates that the growth of the species depends on the past concentration of nutrient. By constructing appropriate Liapunov-like functionals, we obtain sufficient conditions for local and global stability of the positive equilibrium of the model. Quantitative estimates on the size of the delays for local and global stability are also obtained with the help of the Liapunov-like functionals. The technique we use in this paper may be used as well to study global stability of other types of physical models with distributed delays.

Short Title: Global Stability in Chemostat-Type Equations

Keywords: Chemostat-type equations, distributed delay, Liapunov functionals, local and global stability, nutrient recycling.

AMS Subject Classifications: 34K15, 34K20, 45J05, 92A15

Typeset by AMS-TEX

†,School of Mathematics and Statistics, University of Sydney, Sydney, NSW 2006, Australia ‡Department of Mathematics, Statistics and Computing Science, Dalhousie University, Halifax, Nova Scotia, Canada B3H 3J5. Research supported by the NSERC of Canada

§Department of Mathematics and Statistics, McMaster University, Hamilton, Ontario, Canada L8S 4K1.

1. Introduction

The effect of material (nutrient) recycling on ecosystem stability has been previously studied for closed systems (see Nisbet and Gurney [16], Nisbet, McKinstry and Gurney [17] and Ulanowicz [22]). Powell and Richerson [18] and Nisbet and Gurney [16] regarded nutrient recycling as an instantaneous process, thus neglecting the time required to regenerate nutrient from dead biomass by bacterial decomposition. However, as pointed out in Whittaker [23], a delay in nutrient recycling is always present in a natural system and it increases when temperature decreases. To simulate the growth of planktonic communities of unicellular algae in the lakes, Beretta, Bischi and Solimano [1] proposed an open system in which a single species feeds on a limiting nutrient supplied at a constant rate. They assumed that the nutrient is partially recycled after the death of the organisms and used a distributed delay to model nutrient recycling. Bischi [5] observed that the delay involved in nutrient recycling alone does not have destabilizing effect on the equilibrium.

Evidence of delayed growth response has also been observed from chemostat experiments with microalgae Chilamidomonas Reinhardii even when the limiting nutrient is at undetectable small concentration (see Caperon [7]). Following Caperon [7], Ruan [19] introduced a discrete delay to the model of Beretta, Bischi and Solimano [1] to describe the delayed growth response of the species to nutrient uptake. It is shown (see He and Ruan [11]) that the positive equilibrium is globally stable if the delays are sufficiently small. However, there is a threshold value of the discrete delay involved in the growth response, when the discrete delay is increased and passes this critical value, the equilibrium may lose its stability and a Hopf bifurcation may occur (see Ruan [19]). Recently, Beretta and Takeuchi [2-4] used an additional distributed delay to model the delayed growth response. By assuming that the response function is either a Lotka-Volterra function or a Michaelis-Menten function, they studied the global stability of the positive equilibrium.

In this paper, we consider a chemostat-type model with (distributedly) delayed growth response and (distributedly) delayed nutrient recycling, namely, a system of two retarded functional differential equations with two distributed delays. This model was first proposed and studied by Beretta and Takeuchi [2]. However, their stability results are only local. Although global stability was also considered in [2-4] for similar models with an additional instantaneous negative feedback or without delayed growth response at all, the problem is more difficult to study when delayed growth response is introduced. By constructing appropriate Liapunov-like functionals, we study both local and global stability of the positive equilibrium of the model. It turns out that the positive equilibrium can be globally asymptotically stable if the mean delays are sufficiently small, and quantitative estimates on the size of these delays can be obtained with the help of the Liapunov-like functions. Moreover, our approach to the local stability problem is slightly different than that used by Beretta and Takeuchi [2], and we improve their local stability result.

We remark that distributed (infinite) delay equations have been used in biological modeling since the work of Volterra (see Scudo and Ziegler [20]) and they are regarded to be more realistic than discrete (finite) delay equations (see Caperon [7]). The fundamental theory and some properties such as stability, existence of periodic solutions, etc. of distributed delay equations are well understood now and are discussed in the books of Burton [6], Hale and Verduyn Lunel [10], and Hino, Murakami and Naito [12]. The monographs of Cushing [8] and MacDonald [15] give excellent descriptions of distributed delay models and study the local stability, bifurcation, and periodic solutions of these models. Although global stability of some biological models with distributed delays has been studied (see Gopalsamy [9], Kuang [14], Wolkowicz, Xia and Ruan [24] and the references cited therein), in general global results for models involving distributed delays are hard to obtain. The reason probably is that there are few methods available in investigating the global stability of infinite delay equations. The most powerful and most important method is perhaps the Liapunov function(al) method. However, there is no general procedure to follow in constructing a desirable Liapunov function(al) and completely different forms of Liapunov function(al)s are used for different kinds of equations. In the present paper, we try to construct the Liapunov-like functionals step by step so that the idea and technique can be easily followed. We believe that our technique can be used as well to study global stability of some other types of physical models with infinite delays.

The paper is organized as follows. In Section 2 we describe the model equations. Local stability is studied in Section 3 and global stability is discussed in Section 4. Finally, a brief discussion is carried out in Section 5.

2. The Model

Let N(t) denote the limiting nutrient concentration and P (t) denote the plankton concentration at time t. Consider the following integrodifferential equations model of plankton-nutrient interaction with delayed growth response and delayed nutrient recycling:

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

where 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. Throughout, we assume that U(N) is non-negative, increasing, and vanishes when there is no nutrient, and there is a saturation effect when the nutrient is very abundant. That is, we assume that U(N) is a continuously differentiable function defined on [0, ∞) and

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

These general hypotheses are satisfied by the Michaelis-Menten function

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

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

The delay-kernels f(s) and g(s) are non-negative bounded functions defined on [0, ∞). f(s) describes the contribution of the plankton population dead in the past to the nutrient recycled and g(s) describes the delayed growth response of plankton to nutrient uptake. The presence of the distributed time delays must not affect the equilibrium values, so we normalize the kernels such that

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

As in MacDonald [15], we define the average time delays as

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

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.3)

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.4)

Throughout, we always assume that (2.3) is satisfied and T_{f} and T_{g} are finite.

System (2.1) was introduced and studied by Beretta and Takeuchi [2, 4]. It was shown there that all solutions of system (2.1) are nonnegative if the initial data chosen from BC(−∞, 0] (the Banach space of all continuous bounded functions) are nonnegative. They also discussed stability of the positive equilibrium E^{∗}. However, their stability results about system (2.1) are only local, and global stability results hold only for systems (similar to (2.1)) with an additional instantaneous negative feedback or without delayed growth response at all. The object of this paper is to improve their local stability result and investigate as well the global stability of the positive equilibrium.

3. Local Asymptotic Stability

We first study the local stability of the positive equilibrium E^{∗} = (N^{∗}, P^{∗}). Let

x_{1} = N − N^{∗}, x_{2} = P − P^{∗},

where −N^{∗} ≤ x_{1} < ∞, −P^{∗} ≤ x_{2} < ∞, and define

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

so that −U(N^{∗}) ≤ ξ(x_{1}) < 1 − U(N^{∗}) by assumption (2.2). Then the linearized equations about E^{∗} are

\begin{aligned}ẋ_{1} & = −(D + aP^{∗}U′(N^{∗}))x_{1} − aU(N^{∗})x_{2} + bγ ∫_{0}^{∞} f(s)x_{2}(t − s)ds \\ ẋ_{2} & = cP^{∗}U′(N^{∗}) ∫_{0}^{∞} g(s)x_{1}(t − s)ds. \end{aligned}  (3.2)

For convenience, we define

A = D + aP^{∗}U′(N^{∗}), B = acP^{∗}U(N^{∗})U′(N^{∗}), C = bγcP^{∗}U′(N^{∗}). (3.3)

Then B > C, and system (3.2) becomes

\begin{aligned}ẋ_{1} & = −Ax_{1} − \frac{B}{C} bγx_{2} + bγ ∫_{0}^{∞} f(s)x_{2}(t − s)ds \\ ẋ_{2} & = \frac{C}{bγ} ∫_{0}^{∞} g(s)x_{1}(t − s)ds. \end{aligned}  (3.4)

Let (x_{1}(t), x_{2}(t)) be an arbitrary solution of system (3.4). We first consider the function V_{11}(t) = x_{1}^{2}(t). It follows from (3.4) that

V̇_{11}(t) = 2x_{1}(t)ẋ_{1}(t)
= −2Ax_{1}^{2}(t) − 2bγ \frac{B}{C} x_{1}(t)x_{2}(t) + 2bγx_{1}(t)x_{2}(t)
− 2bγx_{1}(t) ∫_{0}^{t} f(s) ∫_{t−s}^{t} ẋ_{2}(u)duds + I(t)
= −2Ax_{1}^{2}(t) − 2bγ(\frac{B}{C} − 1)x_{1}(t)x_{2}(t) + I(t)
− 2Cx_{1}(t) ∫_{0}^{t} f(s) ∫_{t−s}^{t} ∫_{0}^{∞} g(v)x_{1}(u − v)dvduds
≤ −2Ax_{1}^{2}(t) − 2bγ(\frac{B}{C} − 1)x_{1}(t)x_{2}(t) + I(t)
+ C ∫_{0}^{∞} f(s) ∫_{t−s}^{t} ∫_{0}^{∞} g(v)[x_{1}^{2}(t) + x_{1}^{2}(u − v)]dvduds
= −2Ax_{1}^{2}(t) − 2bγ(\frac{B}{C} − 1)x_{1}(t)x_{2}(t) + CT_{f}x_{1}^{2}(t)
+ C ∫_{0}^{∞} f(s) ∫_{t−s}^{t} ∫_{0}^{∞} g(v)x_{1}^{2}(u − v)dvduds + I(t),  (3.5)

where

I(t) = −2bγx_{1}(t) ∫_{t}^{∞} f(s)(x_{2}(t) − x_{2}(t − s)) ds.  (3.6)

For technical reasons, we assume that ∫_{0}^{∞} s^{2}f(s) ds < ∞. Then the function

V_{12}(t) = C ∫_{0}^{∞} f(s) ∫_{t−s}^{t} ∫_{r}^{t} ∫_{0}^{∞} g(v)x_{1}^{2}(u − v)dvdudrds

is well defined, and by (3.5), we have

V̇_{11}(t) + V̇_{12}(t) ≤ −2Ax_{1}^{2}(t) − 2bγ(\frac{B}{C} − 1)x_{1}(t)x_{2}(t) + CT_{f}x_{1}^{2}(t)
+ C ∫_{0}^{∞} f(s) ∫_{t−s}^{t} ∫_{0}^{∞} g(v)x_{1}^{2}(t − v)dvduds + I(t)
= −2Ax_{1}^{2}(t) − 2bγ(\frac{B}{C} − 1)x_{1}(t)x_{2}(t) + CT_{f}x_{1}^{2}(t)
+ CT_{f} ∫_{0}^{∞} g(s)x_{1}^{2}(t − s) ds + I(t).  (3.7)

We now consider the function

V_{1}(t) = V_{11}(t) + V_{12}(t) + CT_{f} ∫_{0}^{∞} g(s) ∫_{t−s}^{t} x_{1}^{2}(u)duds.

It follows from (3.7) that

V̇_{1}(t) ≤ −2(A − CT_{f})x_{1}^{2}(t) − 2bγ(\frac{B}{C} − 1)x_{1}(t)x_{2}(t) + I(t).  (3.8)

On the other hand, by the second equation of (3.4), we have

\frac{d}{dt}[x_{2} + \frac{C}{bγ} ∫_{0}^{∞} g(s) ∫_{t−s}^{t} x_{1}(u)duds] = \frac{C}{bγ} x_{1}(t).

Define

V_{2}(t) = [x_{2} + \frac{C}{bγ} ∫_{0}^{∞} ∫_{t−s}^{t} g(s)x_{1}(u)du]^{2} + (\frac{C}{bγ} )^{2} ∫_{0}^{∞} g(s) ∫_{t−s}^{t} ∫_{v}^{t} x_{1}^{2}(u)dudvds.

We find that

V̇_{2}(t) ≤ \frac{2C}{bγ} x_{1}(t)x_{2}(t) + 2(\frac{C}{bγ} )^{2}T_{g}x_{1}^{2}(t).  (3.9)

Therefore, for the function

V (t) = V_{1}(t) + \frac{(bγ)^{2}}{C} ( ∫_{t}^{∞} f(s) ds + \frac{B}{C} − 1)V_{2}(t),

we have from (3.6), (3.8) and (3.9) that

V̇ (t) ≤ − 2[A − CT_{f} − (C ∫_{t}^{∞} f(s) ds + B − C)T_{g}]x_{1}^{2}(t)
+ 2bγx_{1}(t)x_{2}(t) ∫_{t}^{∞} f(s) ds + I(t)

= − 2[A − CT_{f} − (C ∫_{t}^{∞} f(s) ds + B − C)T_{g}]x_{1}^{2}(t)
+ 2bγx_{1}(t) ∫_{t}^{∞} f(s)x_{2}(t − s) ds
≤ − 2[A − CT_{f} − (C ∫_{t}^{∞} f(s) ds + B − C)T_{g}]x_{1}^{2}(t)
+ 2bγ|x_{1}(t)|∥ϕ_{2}∥ ∫_{t}^{∞} f(s) ds
≤ − 2[A − CT_{f} − (C ∫_{t}^{∞} f(s) ds + B − C)T_{g}]x_{1}^{2}(t)
+ bγx_{1}^{2}(t) ∫_{t}^{∞} f(s) ds + bγ∥ϕ_{2}∥^{2} ∫_{t}^{∞} f(s) ds
= − 2[A − CT_{f} − (B − C)T_{g}]x_{1}^{2}(t)
+ (2CT_{g} + bγ)x_{1}^{2}(t) ∫_{t}^{∞} f(s) ds + bγ∥ϕ_{2}∥^{2} ∫_{t}^{∞} f(s) ds, (3.10)

where ϕ_{2} ∈ BC(−∞, 0] is the initial data of x_{2}(t). By using (3.10), we now can prove the following local stability result.

Theorem 3.1 Assume that ∫_{0}^{∞} s^{2}f(s) ds < ∞. If

CT_{f} + (B − C)T_{g} < A,  (3.11)

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

Proof. Let (x_{1}(t), x_{2}(t)) be an arbitrary solution of (3.4) with ϕ_{2} ∈ BC(−∞, 0] being the initial data for x_{2}(t). By (3.11), we can find ε > 0 such that

Q(ε) ≜ CT_{f} + (B − C)T_{g} + (CT_{g} + \frac{1}{2}bγ)ε < A.

Let T = T (ε) > 0 be such that ∫_{t}^{∞} f(s) ds < ε for all t ≥ T. It then follows from (3.10) that for all t ≥ T,

V̇ (t) ≤ −2(A − Q(ε))x_{1}^{2}(t) + bγ∥ϕ∥^{2} ∫_{t}^{∞} f(s) ds.

Integrating V̇ (t) from T to t ≥ T gives

x_{1}^{2}(t) + \frac{(bγ)^{2}}{C^{2}} (B − C)V_{2}(t) + 2(A − Q(ε)) ∫_{T}^{t} x_{1}^{2}(s) ds
≤ V (t) + 2(A − Q(ε)) ∫_{T}^{t} x_{1}^{2}(s) ds
≤ V (T) + bγ∥ϕ_{2}∥^{2} ∫_{T}^{t} ∫_{s}^{∞} f(u) duds
≤ V (T) + bγ∥ϕ_{2}∥^{2} ∫_{0}^{∞} sf(s) ds
= V (T) + bγ∥ϕ_{2}∥^{2}T_{f} < ∞.

Therefore, x_{1}(t) and x_{2}(t) are bounded, and x_{1}^{2}(t) ∈ L_{1}[0, ∞). By the Mean Value Theorem and the equations in (3.4), x_{1}(t), x_{2}(t) and their derivative functions are thus uniformly continuous on [0, ∞). Applying Barbǎlat Lemma (see Lemmas 1.2.2 and 1.2.3 in Gopalsamy [9]), we conclude that (x_{1}(t), ẋ_{1}(t)) → 0 as t → ∞. Therefore, from the first equation of (3.4), we must have

lim_{t→∞} [ − \frac{B}{C} x_{2}(t) + ∫_{0}^{∞} f(s)x_{2}(t − s) ds] = 0.  (3.12)

Let α = lim inf_{t→∞} x_{2}(t), β = lim sup_{t→∞} x_{2}(t) and {t_{m}} ↑ ∞ be a sequence such that x_{2}(t_{m}) → β as m → ∞. Then β < ∞, and from (3.12) we obtain

\frac{B}{C} β = lim_{m→∞} ∫_{0}^{∞} f(s)x_{2}(t_{m} − s) ds ≤ β.

Since B > C, this implies that β ≤ 0. A similar argument shows that α ≥ 0. Therefore α = β = 0, and (x_{1}(t), x_{2}(t)) → (0, 0) as t → ∞ for every solution (x_{1}(t), x_{2}(t)) of system (3.4).

Note that the characteristic equation of (4.3) is

Δ(λ) = λ^{2} + Aλ + G(λ)(B − CF (λ)) = 0,

where

F (λ) = ∫_{0}^{∞} f(s)e^{−λs}ds, G(λ) = ∫_{0}^{∞} g(s)e^{−λs}ds.

Since B > C and every solution of (3.4) approaches zero as t → ∞, Δ(λ) has no roots with Re(λ) ≥ 0. Therefore, all roots of Δ(λ) have negative real parts and E^{∗} is thus locally asymptotically stable. This completes the proof. □

Remark 3.2 Beretta and Takeuchi [2, 4] observed that system (3.4) has the same characteristic equation as the following system has

\begin{aligned}ẏ_{1} & = y_{2} \\ ẏ_{2} & = −Ay_{2} − ∫_{0}^{∞} Bg(s)y_{1}(t − s)ds + C ∫_{0}^{∞} ∫_{0}^{∞} g(s − v)f(v)y_{1}(t − s)dvds. \end{aligned}

(3.13)

They then constructed a Liapunov functional for (3.13) and showed that the sufficient condition for the local stability of E^{∗} is (Theorem 2 in [4])

CT_{f} + (B + C)T_{g} < A.

Our condition (3.11) improves the above condition.

4. Global Asymptotic Stability

To study the global stability of E^{∗}, we make the change of variables:

x_{1} = N − N^{∗}, x_{2} = ln(P/P^{∗}),  (4.1)

and define ξ(x_{1}) as in (3.1). Then

N = x_{1} + N^{∗}, P = P^{∗} exp(x_{2})  (4.2)

and x_{1}ξ(x_{1}) > 0 for any x_{1} ∈ [−N^{∗}, +∞), x_{1}ξ(x_{1}) = 0 if and only if x_{1} = 0. Using (4.1), we rewrite system (2.1) as follows:

ẋ_{1}(t) = − Dx_{1}(t) − aP^{∗} exp(x_{2}(t))ξ(x_{1}(t)) − P^{∗}G(exp(x_{2}(t)) − 1)
− bγP^{∗} ∫_{0}^{t} ∫_{t−s}^{t} f(s) exp(x_{2}(u))ẋ_{2}(u)du ds + J(t)
ẋ_{2}(t) =c ∫_{0}^{∞} g(s)ξ(x_{1}(t − s))ds,

(4.3)

where G = aU(N^{∗}) − bγ > 0 and

J(t) = −bγP^{∗} ∫_{t}^{∞} f(s)[exp(x_{2}(t)) − exp(x_{2}(t − s))] ds.  (4.4)

Let (x_{1}(t), x_{2}(t)) be an arbitrary solution of system (4.3). We consider the function

V_{11}(t) = ∫_{0}^{x_{1}(t)} ξ(s) ds.

Then upon using (4.3), we obtain

V̇_{11}(t) = −Dx_{1}(t)ξ(x_{1}(t)) − aP (t)ξ^{2}(x_{1}(t)) − P^{∗}G(exp(x_{2}(t)) − 1)ξ(x_{1}(t))

− bγP^{∗}ξ(x_{1}(t)) ∫_{0}^{t} ∫_{t−s}^{t} f(s) exp(x_{2}(u))ẋ_{2}(u)duds + ξ(x_{1}(t))J(t)

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

− P^{∗}G(exp(x_{2}(t)) − 1)ξ(x_{1}(t)) + ξ(x_{1}(t))J(t)

− bcγP^{∗}ξ(x_{1}(t)) ∫_{0}^{t} ∫_{t−s}^{t} f(s) exp(x_{2}(u)) ∫_{0}^{∞} g(v)ξ(x_{1}(u − v))dvduds

≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP (t)ξ^{2}(x_{1}(t)) − P^{∗}G(exp(x_{2}(t)) − 1)ξ(x_{1}(t))

+ \frac{1}{2}bcγ( ∫_{0}^{∞} ∫_{t−s}^{t} f(s)P (u)duds)ξ^{2}(x_{1}(t)) + ξ(x_{1}(t))J(t)

+ \frac{1}{2}bcγ ∫_{0}^{∞} ∫_{t−s}^{t} f(s)P (u) ∫_{0}^{∞} g(v)ξ^{2}(x_{1}(u − v))dvduds.  (4.5)

Let us now consider the following two functions

V_{12}(t) = \frac{1}{2}bcγ ∫_{0}^{∞} f(s) ∫_{t−s}^{t} ∫_{w}^{t} P (u) ∫_{0}^{∞} g(v)ξ^{2}(x_{1}(u − v))dvdudwds,
V_{13}(t) = \frac{1}{2}bcγT_{f} ∫_{0}^{∞} g(s) ∫_{t−s}^{t} P (u + s)ξ^{2}(x_{1}(u))duds.

In what follows, we assume that ∫_{0}^{∞} s^{2}f(s) ds < ∞ so that the function V_{12}(t) above exists. Then for V_{1}(t) = V_{11}(t)+V_{12}(t)+V_{13}(t), we obtain from (4.5) that

\begin{aligned}V̇_{1}(t) & ≤ −Dx_{1}(t)ξ(x_{1}(t)) − aP (t)ξ^{2}(x_{1}(t)) − P^{∗}G(exp(x_{2}(t)) − 1)ξ(x_{1}(t)) \\ & + \frac{1}{2}bcγξ^{2}(x_{1}(t)) ∫_{0}^{∞} ∫_{t−s}^{t} f(s)P (u)duds \end{aligned}

+ \frac{1}{2}bcγT_{f}ξ^{2}(x_{1}(t)) ∫_{0}^{∞} g(s)P (t + s) ds + ξ(x_{1}(t))J(t).  (4.6)

Notice that from the second equation of system (2.1), we have Ṗ (t) ≥ −(γ + D)P (t), for all t > 0. Thus

P (s) ≤ P (t) exp[(γ + D)(t − s)] for t, s > 0.  (4.7)

This implies that

∫_{0}^{∞} ∫_{t−s}^{t} f(s)P (u) duds = ∫_{0}^{t} ∫_{t−s}^{t} f(s)P (u) duds + ∫_{t}^{∞} ∫_{0}^{t} f(s)P (u) duds + K(t)

≤ P (t) ∫_{0}^{t} ∫_{t−s}^{t} f(s) exp[(γ + D)(t − u)]duds

+ P (t) ∫_{t}^{∞} ∫_{0}^{t} f(s) exp[(γ + D)(t − u)]duds + K(t)

≤ \frac{P (t)}{γ + D} ∫_{0}^{∞} f(s)(exp[(γ + D)s] − 1) ds + K(t)

= T_{f}^{∗}P (t) + K(t),  (4.8)

where

K(t) = ∫_{t}^{∞} f(s) ∫_{t−s}^{0} P (u) duds,  (4.9)
T_{f}^{∗} = \frac{1}{γ + D} ∫_{0}^{∞} f(s)(exp[(γ + D)s] − 1) ds.  (4.10)

Similarly, using (4.7), we have

∫_{0}^{∞} g(s)P (t + s)ds ≤ P (t) ∫_{0}^{∞} g(s) exp[−(γ + D)s]ds
= (1 − (γ + D)T_{g}^{∗})P (t),  (4.11)

where

T_{g}^{∗} = \frac{1}{γ + D} (1 − ∫_{0}^{∞} g(s) exp[−(γ + D)s]ds).  (4.12)

If we assume that T_{f}^{∗} in (4.10) is finite, then (4.6), together with (4.8) and (4.11), implies that

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

− P^{∗}G(exp(x_{2}(t)) − 1)ξ(x_{1}(t))

+ \frac{1}{2}bcγ(T_{f}^{∗} + T_{f} − (γ + D)T_{f}T_{g}^{∗})P (t)ξ^{2}(x_{1}(t))

+ ξ(x_{1}(t))J(t) + \frac{1}{2}bcγξ^{2}(x_{1}(t))K(t).  (4.13)

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

\frac{d}{dt}[x_{2}(t) + c ∫_{0}^{∞} g(s) ∫_{t−s}^{t} ξ(x_{1}(u))duds] = cξ(x_{1}(t)).

Let

y(t) = x_{2}(t) + c ∫_{0}^{∞} g(s) ∫_{t−s}^{t} ξ(x_{1}(u))duds.

We define

V_{21}(t) = ∫_{0}^{y(t)}[exp(s) − 1]ds.

Then it follows that

V̇_{21}(t) = cξ(x_{1}(t)){ exp [x_{2}(t) + c ∫_{0}^{∞} g(s) ∫_{t−s}^{t} ξ(x_{1}(u))duds] − 1}
= cξ(x_{1}(t))[ exp(x_{2}(t)) − 1]

+ cξ(x_{1}(t)){ exp [x_{2}(t) + c ∫_{0}^{∞} g(s) ∫_{t−s}^{t} ξ(x_{1}(u))duds] − exp(x_{2}(t))}

= cξ(x_{1}(t))[ exp(x_{2}(t)) − 1]
+ c exp(x_{2}(t))ξ(x_{1}(t)){ exp [c ∫_{0}^{∞} g(s) ∫_{t−s}^{t} ξ(x_{1}(u))duds] − 1}
= cξ(x_{1}(t))[ exp(x_{2}(t)) − 1]
+ c^{2} exp(x_{2}(t))ξ(x_{1}(t)) exp(α(t)) ∫_{0}^{∞} g(s) ∫_{t−s}^{t} ξ(x_{1}(u))duds

for some function α(t) between 0 and c ∫_{0}^{∞} g(s) ∫_{t−s}^{t} ξ(x_{1}(u))duds. Since |ξ(x_{1})| ≤ 1, we have

|α(t)| ≤ c ∫_{0}^{∞} g(s) ∫_{t−s}^{t} |ξ(x_{1}(u))|duds ≤ cT_{g}.

Therefore,

V̇_{21}(t) ≤ cξ(x_{1}(t))[ exp(x_{2}(t)) − 1]

+ \frac{1}{2}c^{2} exp(x_{2}(t)) exp(α(t)) ∫_{0}^{∞} ∫_{t−s}^{t} g(s)[ξ^{2}(x_{1}(t)) + ξ^{2}(x_{1}(u))]duds

= cξ(x_{1}(t))[ exp(x_{2}(t)) − 1]

+ \frac{1}{2}c^{2} exp(x_{2}(t)) exp(α(t))[T_{g}ξ^{2}(x_{1}(t)) + ∫_{0}^{∞} ∫_{t−s}^{t} g(s)ξ^{2}(x_{1}(u))duds]

≤ cξ(x_{1}(t))[ exp(x_{2}(t)) − 1]

+ \frac{1}{2}c^{2}T_{g} exp(x_{2}(t)) exp(cT_{g})ξ^{2}(x_{1}(t))

+ \frac{1}{2}c^{2} exp(x_{2}(t)) exp(cT_{g}) ∫_{0}^{∞} ∫_{t−s}^{t} g(s)ξ^{2}(x_{1}(u))duds.  (4.14)

We now assume that ∫_{0}^{∞} s^{2}g(s) ds < ∞ and define

V_{22}(t) = \frac{1}{2}c^{2} exp(cT_{g}) ∫_{0}^{∞} g(s) ∫_{t−s}^{t} exp(x_{2}(v + s)) ∫_{v}^{t} ξ^{2}(x_{1}(u)dudvds.

Then V_{22}(t) exists and for V_{2}(t) = V_{21}(t) + V_{22}(t), we have from (4.14) that

V̇_{2}(t) ≤ cξ(x_{1}(t))[ exp(x_{2}(t)) − 1]

+ \frac{1}{2}c^{2}T_{g} exp(cT_{g}) exp(x_{2}(t))ξ^{2}(x_{1}(t))

+ \frac{1}{2}c^{2} exp(cT_{g})ξ^{2}(x_{1}(t)) ∫_{0}^{∞} g(s) ∫_{t−s}^{t} exp(x_{2}(u + s))duds.

Note that from (4.2) and (4.7) we have

∫_{t−s}^{t} exp(x_{2}(u + s))du = ∫_{t}^{t+s} exp(x_{2}(u))du

≤ \frac{1}{γ + D} (1 − exp [ − (γ + D)s]) exp(x_{2}(t)).

Therefore,

V̇_{2}(t) ≤ cξ(x_{1}(t))[ exp(x_{2}(t)) − 1]
+ \frac{1}{2}c^{2} exp(cT_{g})(T_{g} + T_{g}^{∗}) exp(x_{2}(t))ξ^{2}(x_{1}(t)).  (4.15)

We finally define the following function

V (t) = V_{1}(t) + \frac{P^{∗}}{c} (bγ ∫_{t}^{∞} f(s) ds + G)V_{2}(t).

It then follows from (4.13) and (4.15) that

V̇ (t) ≤ −Dx_{1}(t)ξ(x_{1}(t)) − \frac{1}{2}(2a − L)P (t)ξ^{2}(x_{1}(t))

+ \frac{1}{2}bcγP^{∗} exp(cT_{g})(T_{g} + T_{g}^{∗}) exp(x_{2}(t))ξ^{2}(x_{1}(t)) ∫_{t}^{∞} f(s) ds

+ bγP^{∗}ξ(x_{1}(t)) ∫_{t}^{∞} f(s)[exp(x_{2}(t − s)) − 1] ds

+ \frac{1}{2}bcγξ^{2}(x_{1}(t))K(t),  (4.16)

where

L = bcγ(T_{f}^{∗} + T_{f} − (γ + D)T_{f}T_{g}^{∗}) + cG exp(cT_{g})(T_{g} + T_{g}^{∗}).

Notice that from (4.2), P (u) = P^{∗} exp(ϕ_{2}(u)) for u ≤ 0 and exp(x_{2}(t − s)) = exp(ϕ_{2}(t − s)) for t ≤ s, where ϕ_{2} ∈ BC(−∞, 0] is the initial data for x_{2}(t). By (4.9), we have

|K(t)| ≤ P^{∗} exp(∥ϕ_{2}∥) ∫_{t}^{∞} sf(s) ds,
|∫_{t}^{∞} f(s)[exp(x_{2}(t − s)) − 1]ds| ≤ [exp(∥ϕ_{2}∥) − 1] ∫_{t}^{∞} f(s) ds.

Since |ξ(x_{1})| ≤ 1, using the above inequalities, we obtain from (4.16) that

V̇ (t) ≤ − Dx_{1}(t)ξ(x_{1}(t)) − \frac{1}{2}(2a − L)P (t)ξ^{2}(x_{1}(t))

+ \frac{1}{2}bcγ exp(cT_{g})(T_{g} + T_{g}^{∗})P (t)ξ^{2}(x_{1}(t)) ∫_{t}^{∞} f(s) ds

+ bγP^{∗}[exp(∥ϕ_{2}∥) − 1] ∫_{t}^{∞} f(s) ds

+ \frac{1}{2}bcγP^{∗} exp(∥ϕ_{2}∥) ∫_{t}^{∞} sf(s) ds.  (4.17)

The above analysis now leads to the following global stability result.

Theorem 4.1 Let G = aU(N^{∗}) − bγ and T_{f}^{∗} and T_{g}^{∗} be the constants defined in (4.10) and (4.12), respectively. Assume that ∫_{0}^{∞} s^{2}g(s)ds < ∞ and

L ≜ bcγ(T_{f}^{∗} + T_{f} − (γ + D)T_{f}T_{g}^{∗}) + cG exp(cT_{g})(T_{g} + T_{g}^{∗}) < 2a.  (4.18)

Then every positive solution (N(t), P (t)) of system (2.1) approaches E^{∗} as t → ∞. If, in addition,

\frac{1}{2}CT_{f} − c(B − C)T_{g}^{2} ≤ D,  (4.19)

where B and C are the constants defined in (3.3), then the positive equilibrium E^{∗} is globally asymptotically stable.

Proof. Let (x_{1}(t), x_{2}(t)) be an arbitrary solution of system (4.3) with ϕ_{2} ∈ BC(−∞, 0] being the initial data of x_{2}(t). We choose ε > 0 such that

L(ε) ≜ L + \frac{1}{2}bcγ exp(cT_{g})(T_{g} + T_{g}^{∗})ε < 2a,

and find T = T (ε) > 1 such that ∫_{t}^{∞} f(s) ds < ε for all t ≥ T. Notice that L < 2a in (4.18) implies T_{f}^{∗} < ∞. Thus V (t) is well defined and (4.17) holds. Therefore, for all t ≥ T, we have

V̇ (t) ≤ − Dx_{1}(t)ξ(x_{1}(t)) − \frac{1}{2}(2a − L(ε))P (t)ξ^{2}(x_{1}(t))

+ \frac{1}{2}bγP^{∗}[(c + 2) exp(∥ϕ_{1}∥) − 2] ∫_{t}^{∞} sf(s) ds

≤ − Dx_{1}(t)ξ(x_{1}(t)) + M ∫_{t}^{∞} sf(s) ds,

where M = \frac{1}{2}bγP^{∗}[(c + 2) exp(∥ϕ_{1}∥) − 2] > 0. Integrating V̇ (t) from T to t ≥ T now gives

V_{11}(t) + \frac{P^{∗}G}{C}V_{2}(t) + D ∫_{T}^{t} x_{1}(s)ξ(x_{1}(s)) ds

≤ V (t) + D ∫_{T}^{t} x_{1}(s)ξ(x_{1}(s)) ds

≤ V (T) + M ∫_{T}^{t} ∫_{s}^{∞} uf(u) duds

≤ V (T) + M ∫_{0}^{∞} s^{2}f(s) ds

≤ V (T) + 2(γ + D)MT_{f}^{∗} < ∞.

This implies that x_{1}(t) and x_{2}(t) are bounded, and x_{1}(t)ξ(x_{1}(t)) ∈ L_{1}[0, ∞). Since ξ is uniformly continuous on [0, ∞), it follows from the equations in system (4.3) and the boundedness of x_{1}(t) and x_{2}(t) that x_{1}(t)ξ(x_{1}(t)) is also uniformly continuous. Thus, by Barbǎlat Lemma (see Lemma 1.2.2 in [9]), x_{1}(t)ξ(x_{1}(t)) → 0 as t → ∞. This leads to lim_{t→∞} x_{1}(t) = 0 , so lim_{t→∞} N(t) = N^{∗} by (4.2). Note further that ẋ_{1}(t) is also uniformly continuous on [0, ∞) by the equations in (4.3). Applying Barbǎlat Lemma (see Lemma 1.2.3 in [9]) once again gives lim_{t→∞} ẋ_{1}(t) = lim_{t→∞} Ṅ(t) = 0. Now taking the limit ( t → ∞ ) on both sides of (2.1), we obtain

lim_{t→∞} (aU(N^{∗})P (t) − bγ ∫_{0}^{∞} f(s)P (t − s) ds) = D(N^{0} − N^{∗}).  (4.20)

Let α = lim inf_{t→∞} P (t), β = lim sup_{t→∞} P (t) and {t_{m}} ↑ ∞ be a sequence such that lim_{m→∞} P (t_{m}) = β. Since β < ∞, it follows from (4.20) that

\begin{aligned}aU(N^{∗})β & = D(N^{0} − N^{∗}) + lim_{m→∞} bγ ∫_{0}^{∞} f(s)P (t_{m} − s) ds \\ & ≤ D(N^{0} − N^{∗}) + bγβ. \end{aligned}

Thus,

β ≤ \frac{D(N^{0} − N^{∗})}{aU(N^{∗}) − bγ} = P^{∗}.

Similarly, we can show that α ≥ P^{∗}. Therefore, lim_{t→∞} P (t) = P^{∗}. This proves the global attractivity of E^{∗}.

To see the asymptotic stability of E^{∗} under the additional condition (4.19), we first note that T_{f} < T_{f}^{∗} and T_{g} < T_{g}^{∗} < 1/(γ + D). Keeping (4.18) in mind, we find that the assumption L < 2a implies

\begin{aligned}bcγT_{f} + cGT_{g} & < a + \frac{1}{2}bcγ(γ + D)T_{f}T_{g}^{∗} − cG(exp(cT_{g}) − 1)T_{g} \\ & < a + \frac{1}{2}bcγT_{f} − c^{2}GT_{g}^{2}. \end{aligned}

Recall (3.3) and G = aU(N^{∗}) − bγ. The above inequality, together with the assumption (4.19), gives

\begin{aligned}CT_{f} + (B − C)T_{g} & < aP^{∗}U′(N^{∗}) + \frac{1}{2}CT_{f} − c(B − C)T_{g}^{2} \\ & ≤ aP^{∗}U′(N^{∗}) + D = A. \end{aligned}

Thus (3.11) holds. Note that ∫_{0}^{∞} s^{2}f(s) ds < 2(γ + D)T_{f}^{∗} < ∞. By Theorem 3.1, E^{∗} is locally asymptotically stable. This, together with global attractivity, implies that E^{∗} globally asymptotically stable. □

5. Discussion

In this paper, we have considered a chemostat-type plankton model with nutrient recycling. We assumed that there are a (distributed) delay in the growth response of plankton to nutrient uptake and a (distributed) delay in nutrient recycling. We have obtained some sufficient conditions for both local and global stability of the positive equilibrium by constructing appropriate Liapunov-like functionals. It is known that the delay in the growth response of the populations to nutrient uptake can cause oscillations in population density (see Caperon [7]and Ruan [19]), our results indicate that one can still have global stability of the positive equilibrium if the delays are sufficiently small, and explicit estimates (see (3.11), (4.18) and (4.19)) on the size of these delays for local and global stability can also be obtained.

We should pointed out that model (2.1) and other related models have been studied by Beretta and Takeuchi [2-4] and Kolmanovskii, Torelli and Vermiglio [13].

For local stability, following the arguments of Kolmanovskii, Torelli and Vermiglio [13], Beretta and Takeuchi [2, 4] chose a Liapunov functional for the equivalent system (3.13) and obtained some sufficient conditions. Our local stability conditions improve their conditions. For global stability, Beretta and Takeuchi [2, 4] only considered the special cases of model (2.1) when either the second delay does not appear or the growth response function is a Lotka-Volterra or a Michaelis-Menten function. Our global stability results hold for the general case when both distributed delays present and the growth response function is a general function satisfying (2.2).

References

[1] E. Beretta, G. I. Bischi and F. Solimano, Stability in chemostat equations with delayed nutrient recycling, J. Math. Biol., 28(1990), pp. 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), pp. 19-40.

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

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

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

[6] T. A. Burton, Volterra Integral and Differential Equations, Academic Press, New York, 1983.

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

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

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

[10] J. K. Hale and S. M. Verduyn Lunel, Introduction to Functional Differential Equations, Appl. Math. Sci. 99, Springer-Verlag, New York, 1993.

[11] X.-Z. He and S. Ruan, Global stability in chemostat-type plankton models with delayed nutrient recycling, Research Report 96-8, The School of Math. & Stat., University of Sydney, 1996.

[12] Y. Hino, S. Murakami and T. Naito, Functional Differential Equations with Infinite Delay, Lecture Notes in Math. 1473, Springer-Verlag, New York, 1991.

[13] V. B. Kolmanovskii, L. Torelli and R. Vermiglio, Stability of some test equations with delay, SIAM J. Math. Anal., 25(1994), pp. 948-961.

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

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

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

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

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

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

[20] F. M. Scudo and J. R. Ziegler, The Golden Age of Theoretical Ecology: 1923-1940, Lecture Notes in Biomathematics 22, Springer-Verlag, Berlin, 1978.

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

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

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

[24] G. S. K. Wolkowicz, H. Xia and S. Ruan, Competition in the chemostat: a distributed delay model and its global asymptotic behavior, SIAM J. Appl. Math., to appear.