Electronic Journal of Qualitative Theory of Differential Equations Proc. 8th Coll. QTDE, 2008, No. 11 1-18; http://www.math.u-szeged.hu/ejqtde/

Neutral Delay Equations from and for Population Dynamics

K.P. Hadeler^{1} Department of Mathematics and Statistics Arizona State University Tempe, AZ 85287, U.S.A.

Abstract

For a certain class of neutral differential equations it is shown that these equations can serve as population models in the sense that they can be interpreted as special cases or caricatures of the standard Gurtin-MacCamy model for a population structured by age with birth and death rate depending on the total adult population. The delayed logistic equation does not belong to this class but the blowfly equation does. These neutral delay equations can be written as forward systems of an ordinary differential equation and a shift map. There are several quite distinct ways to perform the transformation to a system, either following a method of Hale or following more closely the renewal process. Similarly to the delayed logistic equation, the neutral equation (and the blowfly equation as a special case) exhibit periodic solutions, although only for a restricted range of parameters.

subject classification: Primary 34K40, 34K17 ; Secondary 92D25

key words: neutral delay equation, population model, blowfly equation, quasilinear

1 Introduction

The classical model for a single population is the Verhulst equation

u̇ = au(1 − u/K)  (1.1)

1email: hadeler@uni-tuebingen.de This paper is in final form and no version of it is submitted for publication elsewhere.

which is based on the idea of an a priori given carrying capacity K. Another approach uses a birth rate b(u) and a death rate μ(u) depending on population size,

u̇ = b(u)u − μ(u)u.  (1.2)

In this case a carrying capacity K can be defined a posteriori as the positive solution of the equation b(u) = μ(u) (provided the solution exists and is unique). Experiments and field observations have suggested that isolated populations can show oscillatory behavior even under constant nutrient supply. Such oscillations appear to differ from predator-prey cycles and they have been explained as being caused by delays. The eminent ecologist G.E. Hutchinson ([18], [19]) thought of a delayed response to diminishing resources and suggested the equation

u̇(t) = au(t) (1 − \frac{u(t − τ)}{K})  (1.3)

as a modification of (1.1) while others, starting from experiments [31], [32], thought of a delay in maturation and arrived at what is now called the blowfly equation [35],

u̇(t) = b(u(t − τ))u(t − τ) − μ(u(t))u(t).  (1.4)

The blowfly equation has been investigated and extended in [10], [39], it has been rediscovered many times, a recent example is [1].

The equations (1.3) and (1.4) behave quite differently with respect to the onset of oscillations [12]. It seems difficult to justify Hutchinson’s equation in terms of basic biological assumptions ([28], [29]) whereas the blowfly equation can be derived from age-structure models.

Neutral delay equations have been suggested as population models occasionally [22]. But usually they did not fit well into the available theory of neutral delay equations. In general neutral delay equations have been considered to be difficult in comparison with standard delay equations.

In [3] it has been shown that there is a large class of neutral delay equations of the general form

u̇(t) = {b_{1}(u(t − τ))u(t − τ)
+b_{2}(u(t − τ)) [u̇(t − τ) + μ_{1}(u(t − τ))u(t − τ)] }
−μ_{1}(u(t))u(t)  (1.5)

which can be derived by projection from the Gurtin-MacCamy age structure model and thus can be justified in basic biological terms. The blowfly equation is contained as a special (non-neutral) case for b_{2} = 0. It turns out that this class of equations is more or less identical with a class of neutral equations that have been studied by J. Hale and coauthors [17], [4], [16], and, in particular, [15], p.35.

The equation (1.5) can be written (in many different ways) as a coupled system of an ordinary differential equation and a shift map, for instance as a system

u̇(t) = z(t) − μ_{1}(u(t))u(t)
z(t) = b_{1}(u(t − τ))u(t − τ) + b_{2}((t − τ))z(t − τ)  (1.6)

for the (adult) population size u and the rate z of entry into the adult class.

So, what does one achieve by replacing one infinite-dimensional system (the Gurtin-MacCamy system) by another, i.e., (1.5) and (1.6)? First of all, one can identify those delay equations which allow a viable interpretation as population models. In particular, one finds a class of neutral equations with this property. The derivation tells further, which biological mechanism causes the neutral term, and that the neutral term always comes together with another term (see the b_{2} term in (1.5)). This is to say, one cannot get a viable model from just putting a neutral delay term into, say, the blowfly equation. Second, systems of the form (1.6) are easy to handle for numerical analysis, they depend on few parameters and they are amenable to the theory of delay equations, see, e.g., [41]. Second, our approach distinguishes a somewhat larger class of neutral delay equations (which we call quasi-linear) for which a reduction to a system of a shift and an ordinary differential equations is possible.

How does our approach fit into the general theory of population models with age (or size) structure? The Gurtin-MacCamy system has been analyzed in great detail in [36], [42]. At present there are several rather independent lines of research using operator semi-groups which still reflect the two classical approaches, by integral equations (Sharpe-Lotka) and by partial differential equations (McKendrick). The Volterra integral equations approach ([9], [6], using adjoint semigroups, [7] [8]) states that “every reasonable structured population model (not only age structured) can be written as a system of Volterra integral equations (delay equations) of the form x(t) = F (x_{t}) or x(t) = F_{1}(x_{t}, y_{t}), ẏ(t) = F_{2}(x_{t}, y_{t})” (x_{t} is the usual t-segment of the function x). Of course (1.6) has this form. However, our goal is not further generalization but getting delay equations. The other approach is based on operator semigroup theory following [36], [42]. The recent paper [26] based on integrated semigroups includes an excellent review of earlier work, These authors also discuss the relation between neutral delay equations and population models exhibited in [3].

The paper is organized as follows. In section 2 we start from the Gurtin-MacCamy model, assume that the coefficients are of a simple but highly singular type, and derive a system of coupled delay equations. In section 3 we discuss the blow-fly equation as a special case. In section 4 we derive a neutral equation as a special case and we present several equivalent systems of a shift map coupled to an ordinary differential equation. In section 5 we study a general class of neutral equations which we call “quasi-linear” and we establish the connection to the theory in [15] [16]. In section 6 we address forward continuation and in section 7 linear stability.

2 The Gurtin-MacCamy model

Linear models for mortality have been designed by Euler and Daniel Bernoulli in the eighteenth century. Sharpe and Lotka included a birth law and designed a model in the form of a Volterra integral equation or renewal equation. McKendrick 1926 cast this model in the form of a hyperbolic differential equation. Feller 1941 proved for the Volterra equation a renewal theorem, i.e., convergence to the unique persistent (exponential) solution. The Sharpe-Lotka-McKendrick model is linear and hence does not show a saturation effect. Gurtin and MacCamy [11] modified the model by letting the birth and death rates depend on total population size, We consider the following system of Gurtin-MacCamy type

\frac{∂u(t, a)}{∂t} + \frac{∂u(t, a)}{∂a} + μ(a, W (t))u(t, a) = 0

u(t, 0) = ∫_{0}^{∞} b(a, W (t))u(t, a)da

W (t) = ∫_{0}^{∞} ρ(a)u(t, a)da  (2.1)

where W is the weighted population size. It is well known that in some cases this model can be “reduced” to a system of ordinary differential equations. If the functions μ, b and ρ ≡ 1 do not depend on age a then the total population size U(t) = ∫_{0}^{∞} u(t, a)da satisfies U̇ = b(U)U − μ(U)U. More generally, if the functions b(a, W), μ(a, W) are exponential polynomials of the variable a, with W as a parameter, and ρ ≡ 1, then certain functionals P_{i}(t) = ∫_{0}^{∞} p_{i}(a)u(t, a)da satisfy a system of ordinary differential equations, see [23] [24], [2] [30], [3]. If the exponential polynomials have many terms then the system of ordinary differential equations gets cumbersome, and this “chain trickery approach” is of little use, even more so, as the positive cone of the o.d.e. system is not the image of the positive cone of the Gurtin-MacCamy system under the reducing projection, see [3]. A similar reduction works for delay equations with exponentially distributed delays.

Here we apply another reduction principle which has been used by several authors [2] [37] [38] [40] (for a general view on reducible systems see [14]). We assume that the coefficients are piece-wise constant functions of a. We even allow that some of these functions have delta peaks. In general it does make sense to study a hyperbolic system with delta peaks in the coefficients but in the case of the Gurtin-MacCamy system these coefficients can be justified if the system is reformulated in terms of renewal equations. We assume a single jump or peak at the age a = τ where τ can be seen as the length of the juvenile period. Hence the coefficients are

μ(a) = μ_{0}(W) + (μ_{1}(W) − μ_{0}(W))H_{τ}(a)
b(a) = b_{1}(W)H_{τ}(a) + b_{2}(W)δ_{τ}(a)
ρ(a) = \{ \begin{matrix}α & a < τ \\ β & a ≥ τ \end{matrix}  (2.2)

where δ_{τ}(a) is the delta peak at a = τ and H_{τ}(a) is the Heaviside function, H_{τ}(a) = 0 for a < τ and H_{τ}(a) = 1 for a ≥ τ. Hence the mortality jumps from μ_{0} to μ_{1}; the fertility jumps from 0 to b_{1} with a delta peak at a = τ. In [3] some juvenile fertility and a delta peak in the mortality have also been allowed. Now introduce the variables

V (t) = ∫_{0}^{τ} u(t, a)da, U(t) = ∫_{τ}^{∞} u(t, a)da  (2.3)

which represent the total juvenile population and the total adult population. To ease the presentation here we have exchanged the roles of U and V with respect to [3].

Suppose u(t, a) is a solution of the system (2.1) with coefficients (2.2). For 0 < t < τ the variables V and U satisfy a non-autonomous system of ordinary differential equations whereby the coefficients depend on the history, i.e., on the values of the initial data u(0, a) = u_{0}(a) in the interval 0 ≤ a ≤ τ. For t > τ the variables U and V satisfy an autonomous system of neutral delay differential equations

V̇ (t) = b_{1}(W (t))U(t) − μ_{0}(W (t))V (t)
+(b_{2}(W (t)) − 1){b_{1}(W (t))U(t)
+b_{2}(W (t)) [U̇(t − τ) + μ_{1}(W (t − τ))U(t − τ)] }
× exp{− ∫_{0}^{τ} μ_{0}(W (t − τ + σ))dσ}
U̇(t) = {b_{1}(W (t − τ))U(t − τ)
+b_{2}(W (t − τ)) [U̇(t − τ) + μ_{1}(W (t − τ))U(t − τ)] }
× exp{− ∫_{0}^{τ} μ_{0}(W (t − τ + σ))dσ} − μ_{1}(W (t))U(t)
W (t) = αV (t) + βU(t).  (2.4)

For the initial system we refer to [3]. Here we assume that the population has evolved at least for a time interval of length τ such that (2.4) applies.

Suppose that W depends only on the adult population, i.e., α = 0, β = 1. Then the system (2.4) decomposes and the equation for U does not depend on V. We get a scalar neutral differential delay equation

U̇(t) = {b_{1}(U(t − τ))U(t − τ)
+b_{2}(U(t − τ)) [U̇(t − τ) + μ_{1}(U(t − τ))U(t − τ)] }
× exp{− ∫_{0}^{τ} μ_{0}(U(t − τ + σ))dσ} − μ_{1}(U(t))U(t) (2.5)

with a point delay and a (nonlinear) distributed delay. The distributed delay describes the effect of the adult population on juvenile mortality. Of course cannibalism is known in some fish and reptile species but for most invertebrate and even vertebrate species such effect does not exist. The system becomes much simpler if μ_{0} is a constant.

3 The blowfly equation

Assume that μ_{0} is a constant and that b_{2} ≡ 0. Then the equation (2.5) reduces to

u̇(t) = b_{1}(u(t − τ))u(t − τ)e^{−μ_{0}τ} − μ_{1}(u(t))u(t).  (3.1)

This is the blowfly equation (1.4) with b(u) = b_{1}(u)e^{−μ_{0}τ} and μ(u) = μ_{1}(u). The formula (3.1) has a very natural interpretation. The equation (3.1) is essentially an equation of the form (1.2) applied to the adult population only. Adults produce offspring with a birth rate b_{1} but this offspring enters the adult population only after time τ and with a discount factor e^{−μ_{0}τ} that accounts for juvenile mortality.

Hence we can state: The blowfly equation describes the dynamics of the Gurtin-MacCamy system for the special case where the juvenile mortality is a constant and the adult mortality and the adult fertility depend only on the adult population density.

In [35] the coefficients have been adapted to data, the special case in [1] has perhaps the simplest coefficients: the birth rate decreases like 1/u and the death rate increases like 1 + u.

Finally we comment on the further history of Hutchinson’s and the blowfly equations. Hutchinson’s equation can be carried into Wright’s equation ẋ(t) = −αx(t − τ)(1 + x(t)), and then to ẋ(t) = f(x(t − τ)) with f(x) = −α(exp{x} − 1). Motivated by applications from biology, an instantaneous feedback term has been introduced into such equation (see [25], [13]) leading to the equation ẋ(t) = f(x(t − τ)) − νx(t), and then parameters have been rescaled to achieve the form ϵẋ(t) = f(x(t − 1)) − x(t) that suggests comparison with the discrete map x(t) = f(x(t − 1)). In spite of seminal papers by Nussbaum [33] [34] on the existence of periodic solutions and the considerable progress in understanding the dynamics of this equation, see the overview [41], and the extension to state-depending delays [27], the basic equation (1.3) remains a somewhat dubious biological model.

4 The neutral delay equation

Now consider the equation (2.5) with constant μ_{0} and general b_{2}. We get a neutral delay differential equation

U̇(t) = {b_{1}(U(t − τ))U(t − τ)

+b_{2}(U(t − τ)) [U̇(t − τ) + μ_{1}(U(t − τ))U(t − τ)] }e^{−μ_{0}τ}
−μ_{1}(U(t))U(t).  (4.1)

So far our approach has shown three things: There is a large class of neutral delay equations which can be seen as population models. This class does not contain Hutchinson’s equation. This class is still special within the set of all nonlinear neutral delay equations.

For a solution U(t) of (4.1) define the function

Z(t) = U̇(t) + μ_{1}(U(t))U(t).  (4.2)

With the variable Z we can rewrite the neutral equation as a system

U̇(t) = Z(t) − μ_{1}(U(t))U(t)
Z(t) = [b_{1}(U(t − τ))U(t − τ) + b_{2}(U(t − τ))Z(t − τ)]e^{−μ_{0}τ}. (4.3)

The variable Z can be interpreted as the number of entries per time into the adult class. Hence the first equation of (4.3), which is an ordinary differential equation, says that surviving juveniles enter the adult class and then die with death rate μ_{1}(U). The second equation is essentially a shift map which yields new values for Z in terms of the history of U and Z. This system is perhaps the most transparent formulation of the neutral equation in biological terms and also well suited for numerical simulations.

We know that a second order differential equation can be written as a first order system in various ways. A neutral delay equation is somewhat similar to a second order equation because of the leading term U̇(t−τ). Hence there are various other ways to write (4.1) as a system of an ordinary differential equation and a shift map. For example, introduce a new variable

B_{2}(U) = ∫_{0}^{U} b_{2}(V)dV e^{−μ_{0}τ}.  (4.4)

Then the equation (4.1) is equivalent with the system

U(t) = X(t) + B_{2}(U(t − τ))
Ẋ(t) = −μ_{1}(X(t) + B_{2}(U(t − τ))(X(t) + B_{2}(U(t − τ)). (4.5)

Here U is again the adult population whereas X is the adult population without those individuals which have been produced by the delta peak. Of course, from a biological point of view, the system (4.3) is more transparent than (4.5). However, from a mathematical point of view, the system (4.5) fits better into the general pattern of quasilinear equations which will be defined in the next section.

5 Quasilinear neutral delay equations

Now we are leaving biological models and we enter the discussion of general neutral delay equations. Therefore, and in accordance with section 1, we use lower case letters for the dependent variable. The equation (4.1) has the general form

u̇(t) = g(t, u(t − τ))u̇(t − τ) + h(t, u(t), u(t − τ))  (5.1)

where in the particular case the coefficient functions g and h are given by

g(t, v) = b_{2}(v)e^{−μ_{0}τ},
h(t, u, v) = (b_{1}(v) + b_{2}(v)μ_{1}(v))e^{−μ_{0}τ} − μ_{1}(u)u.  (5.2)

The distinctive feature of equation (5.1) is that the term u̇(t−τ) (the “highest derivative”) occurs linearly. For this reason we call this neutral equation a “quasi-linear” neutral equation, similar to the notion “quasi-linear” for parabolic and elliptic partial differential equations.

In comparison with equation (5.1) the equation (4.1) is special insofar as the function h splits into a sum of an instant term and a delayed term, i.e., (4.1) has the form

u̇(t) = g(t, u(t − τ))u̇(t − τ) + q(t, u(t)) + h_{0}(t, u(t − τ)).  (5.3)

In [15], p.35, equations of the form (5.1) have been distinguished as a special class of neutral delay equations for which there is a coherent theory [4], [17]. We think it important to stress that there is a close connection between this class of “easy” neutral delay equations and population models. The theory of such equations has been extended to equations with distributed delays [16].

Here we present the nucleus of this theory and relate it to the systems studied in section 4. Similar transitions from neutral equations to systems have been used in [20] [21], and in [43], p. 55.

Assume that in equation (5.1) the functions g and h are continuously differentiable. We introduce the primitive of the coefficient function g as

G(t, u) = ∫_{0}^{u} g(t, z)dz.  (5.4)

Let

G_{t}(t, u) = ∫_{0}^{u} g_{t}(t, z)dz, G_{u}(t, u) = g(t, u)  (5.5)

denote the partial derivatives of the function G with respect to the first and second variable. Define the function

F (t, v, z) = h(t, v + G(t, z), z) − G_{t}(t, z).  (5.6)

The following proposition establishes a connection between the scalar equation and a system. The connection is local. We do not impose conditions which would ensure global existence of solutions.

Proposition 5.1 For every solution u of the equation (5.1) there is a function v such that the pair (u, v) is a solution of the system of an ordinary differential equation and a shift map

u(t) = v(t) + G(t, u(t − τ))
v̇(t) = F (t, v(t), u(t − τ)).  (5.7)

Conversely, if a solution of the system (5.7) is given with a differentiable function u then this function is a solution of equation (5.1).

Proof: We write the equation (5.1) as

\frac{d}{dt}(u(t) − G(t, u(t − τ))) = h(t, u(t), u(t − τ)) − G_{t}(t, u(t − τ)). (5.8)

Then we introduce a new dependent variable v by

v(t) = u(t) − G(t, u(t − τ)).  (5.9)

Then (5.8) becomes

v̇(t) = h(t, u(t), u(t − τ)) − G_{t}(t, u(t − τ)).  (5.10)

We use the definition (5.9) as a dynamic equation for u,

u(t) = v(t) + G(t, u(t − τ))  (5.11)

and get a system

u(t) = v(t) + G(t, u(t − τ))
v̇(t) = h(t, u(t), u(t − τ)) − G_{t}(t, u(t − τ)).  (5.12)

In the second equation we replace u(t) from the first equation and get

u(t) = v(t) + G(t, u(t − τ))
v̇(t) = h(t, v(t) + G(t, u(t − τ)), u(t − τ)) − G_{t}(t, u(t − τ)). (5.13)

Finally use the definition (5.6) and get (5.7).

Now suppose a solution of (5.7) is given with differentiable u. Differentiate the first equation,

u̇(t) = v̇(t) + G_{t}(t, u(t − τ)) + G_{u}(t, u(t − τ))u̇(t − τ),

replace v̇(t) from the second equation, use the definition of F as given by (5.6) and again use the first equation of (5.7) to replace v(t).

□

Next we look at the special equation (5.3).

Proposition 5.2 The equation (5.3) is equivalent with the system

u(t) = v(t) + G(t, u(t − τ))  (5.14)
v̇(t) = q(t, v(t) + G(t, u(t − τ))) + h_{0}(t, u(t − τ)) − G_{t}(t, u(t − τ))

and also with the system

u̇(t) = z(t) − q(t, u(t))  (5.15)
z(t) = g(t, u(t − τ))[z(t − τ) + q(t − τ, u(t − τ))] + h_{0}(t, u(t − τ)).

Proof: (5.14) follows by applying Proposition 5.1 to (5.3) while (5.15) is obtained by applying the first method of section 4 to (5.3).

□

Hence we have two different systems, both consisting of an ordinary differential equation and a shift map, for the neutral differential equation (5.3).

The common idea of both approaches is the following. Because of the term u̇(t − τ) the neutral equation has the character of a “second order equation”. In general, to carry a second order equation into a system one introduces the first derivative as a new variable. Here the “first derivative” is either the difference v(t) = u(t) − G(u(t − τ)) or the modified derivative z(t) = u̇(t) + q(t, u(t)).

In order to better understand these different systems we take a look at the linear non-autonomous equation.

Linear non-autonomous equations

Consider a general linear neutral equation

u̇(t) = a(t)u̇(t − τ) + b(t)u(t − τ) + c(t)u(t) + f(t)  (5.16)

with continuously differentiable a, continuous b, c, and a constant delay τ > 0. The system (5.14) becomes

u(t) = v(t) + a(t)u(t − τ)
v̇(t) = c(t)v(t) + [b(t) − ȧ(t) + a(t)c(t)]u(t − τ) + f(t) (5.17)

whereas the system (5.15) becomes

u̇(t) = z(t) + c(t)u(t)
z(t) = a(t)z(t − τ) + [b(t) + a(t)c(t − τ)]u(t − τ) + f(t). (5.18)

For (5.18) continuity of the coefficient a suffices.

6 Forward continuation

We assume, for simplicity, that the function h in (5.1) has a global Lipschitz constant with respect to the second variable u(t).

For the equation (5.1) we have the equivalent system (5.7) and for the special equation (5.3) we have even two equivalent systems (5.14) and (5.15). In each of these systems a shift map is joined to an ordinary differential equation. Here we look at suitable state spaces and forward continuation.

For the system (5.7) a natural choice of a state variable is (u_{t}, v(t)) ∈ C[−τ, 0] × IR. If (ϕ, ψ) is an initial data in this space then we can find a unique forward continuation (u_{t}, v(t)) with (u_{0}, v(0)) = (ϕ, ψ) in the following manner. First solve the differential equation for v with the initial condition v(0) = ψ in the interval 0 ≤ t ≤ τ. Then use the first equation (v(t) is now known) to determine u in the interval 0 < t ≤ τ.

Proposition 6.1 For the system (5.14) the u component of the solution is continuous at t = 0. The function v is continuous at t = 0 if and only if

ϕ(0) = ψ + G(0, ϕ(−τ)).  (6.1)

The u component is continuously differentiable at t = 0 if and only (6.1) is satisfied and also ϕ ∈ C^{1}[−τ, 0] and

ϕ̇(0) = g(0, ϕ(−τ))ϕ̇(−τ) + h(0, ϕ(0), ϕ(−τ)).  (6.2)

Proof: Continuity at t = 0 says u(0) = v(0)+G(0, u(−τ)). Use u(0) = ϕ(0) and v(0) = ψ to get (6.1). Continuity of u̇ at t = 0 says

u̇(0) = v̇(0) + G_{t}(0, u(−τ)) + g(0, u(−τ))u̇(−τ)).

Replace v̇(0) (the G_{t} term cancels) and get (6.2).

□

The condition (6.2) just says that the initial data satisfy the neutral differential equation (5.1) which makes sense. Continuity and differentiability, respectively, at t > 0 implies the same property at t = τ, 2τ, . . ..

Now consider the system (5.15). In this case a natural choice of a state is (u(t − τ), z_{t}) ∈ IR × C[−τ, 0]. If (u(−τ), z_{0}) = (ψ, ϕ) is given, then use the first equation to compute u(t) for −τ ≤ t ≤ 0 and then use the second equation to compute z(t) for 0 < t ≤ τ. Then z is not necessarily continuous at t = 0.

Proposition 6.2 For the system (5.15) the function u is continuous at t = 0 by construction. The function z is continuous at t = 0 if and only if

ϕ(0) = g(0, ψ)[ϕ(−τ) − q(−τ, ψ)] + h_{0}(0, ψ).  (6.3)

If this condition is satisfied then u̇ is also continuous at t = 0.

Again, continuity at t = 0 implies continuity at t = τ, 2τ, . . ..

In particular, for the population model (4.3) we find the compatibility condition

ϕ(0) = [b_{1}(ψ)ψ + b_{2}(ψ)ϕ(−τ)]e^{−μ_{0}τ}.  (6.4)

Finally look at the linear equation (5.16). Then the conditions (6.1)(6.2) become

ϕ(0) = ψ − a(0)ϕ(−τ)  (6.5)
ϕ̇(0) = a(0)ϕ̇(−τ) + b(0)ϕ(−τ) + c(0)ϕ(0)  (6.6)

while (6.3) becomes simply

ϕ(0) = a(0)ϕ(−τ) + a(0)c(−τ)ψ + b(0)ψ.  (6.7)

7 Stability

Hutchinson’s goal as well as that of the authors [35] of the blowfly equation has been to explain population oscillations in terms of delays. Hutchinson’s equation has stable non-constant periodic solutions for aτ > π/2. The question of the largest lower bound on aτ for such solutions to exist seems of mainly academic interest (there are no such solutions for aτ < 3/2). Hence Hutchinson’s model suggests that oscillations are most likely to occur once the delay is sufficiently large. In [5] the effect of maturation periods on stability has been studied.

The blowfly equation and the more general neutral equation send the message that not always a large delay causes instability. In [12] it has been shown that the constant solution ū > 0 of the blowfly equation (1.4) becomes unstable and gives rise to periodic oscillations for some τ > 0 if and only if the inequality

\frac{d}{du} [(b_{1}(u) + μ_{1}(u))u] |_{u=ū} < 0  (7.1)

holds. Since b_{1}(u) is a decreasing function and μ_{1} is an increasing function, this inequality can only be satisfied if b′_{1}(ū) is strongly negative. According to the blowfly model, stable periodic oscillation occur if the birth rate is strongly decreasing with increasing population size and if the delay is large.

The corresponding inequality for the neutral equation (1.5) reads

\frac{d}{du} [(b_{1}(u) + (b_{2}(ū) + 1)μ_{1}(u))u] |_{u=ū} < 0.  (7.2)

Hence we have a similar situation as in the blowfly equation. Notice that the equilibrium ū in both equations is not the same. The term b_{2}μ_{1} may be increasing or decreasing. Hence there is no simple statement as to whether the neutral equation is more stable than the blowfly equation or conversely.

References

[1] Arino, J., Wang, Lin, Wolkowicz, G.S.K., An alternative formulation for a delayed logistic equation. J. Theor. Biol. 241, 109-119 (2006)

[2] Blythe, S.P., Nisbet, R.M., Gurney, W.S.C., The dynamics of population models with distributed maturation period. Theor. Population Biol. 25, 289-311 (1984)

[3] Bocharov, G., Hadeler, K.P., Structured population models, conservation laws, and delay equations. J. Diff. Equ. 168, 212-237 (2000)

[4] Cruz, M.A., Hale, J.K., Existence, uniqueness and continuous dependence for hereditary systems. Annali Mat. Pura Appl. 85, 63-82 (1970)

[5] Cushing, J.M., Stability and maturation periods in age structured populations. In: S. Busenberg and K. Cooke, eds.) Differential Equations and Applications in Ecology, Epidemics, and Population Problems, pp 163-181, Academic Press, New York, 1981

[6] Diekmann, O., Getto, P., Boundedness, global existence and continuous dependence for nonlinear dynamical systems describing physiologically structured populations. J. Diff. Equ. 215, 268-219 (2005)

[7] Diekmann, O., Getto, Ph., Gyllenberg, M., Stability and bifurcation analysis of Volterra functional equations in the light of suns and stars. SIAM J. Math. Anal. 39, 1023-1069 (2007)

[8] Diekmann, O., Gyllenberg, M., Abstract delay equations inspired by population dynamics. In: H. Amann et al. (eds) Functional Analysis and Evolution Equations. Birkhäuser, pp 187-200 (2007)

[9] Diekmann, O., Gyllenberg, M., Huang, H., Kirkilionis, M., Metz, J.A.J., Thieme, H.R., On the formulation and analysis of general deterministic structured models. II. Nonlinear theory. J. Math. Biol. 43, 157-189 (2001)

[10] Gurney, W.S.C., Blythe, S.P., Nisbet, R.M., Nicholson’s blowflies revisited. Nature 287, 17-21 (1980)

[11] Gurtin, M., MacCamy, R., Nonlinear age-dependent population dynamics. Arch. Rat. Mech. Anal. 54, 281-300 (1974)

[12] Hadeler, K.P., Bocharov, G., Where to put delays in population models, in particular in the neutral case. Canadian Applied Mathematics Quarterly (CAMQ) 11, 159-173 (2003)

[13] Hadeler, K.P., Tomiuk, J., Periodic solutions of difference-differential equations. Arch. Rat. Mech. Anal. 65, 87-95 (1977)

[14] Hadeler, K.P., Walcher, S., Reducible ordinary differential equations. J. Nonlinear Sci. 16, 583-613 (2006)

[15] Hale, Jack, Theory of Functional Differential Equations. Springer Verlag 1977

[16] Hale, J.K., Verduyn Lunel, S.M., Introduction to Functional Differential Equations. Springer Appl. Math. Series 99, Springer Verlag 1993

[17] Hale, J.K., Meyer, K.R., A class of functional equations of neutral type. Mem. Amer. Mat. Soc. 76 (1967)

[18] Hutchinson, G.E., Circular causal systems in ecology. Ann. New York Acad. Sci. 50, 221-246 (1948)

[19] Hutchinson, G.E., An Introduction to Population Ecology. Yale University Press 1978

[20] Krisztin, T., Wu, J., Monotone semiflows generated by neutral equations with different delays in neutral and retarded parts. Acta Math. Univ. Comen., New Ser. 63, 207-220 (1994)

[21] Krisztin, T., Wu, J., Asymptotic periodicity, monotonicity, and oscillation of solutions of scalar neutral functional differential equations. J. Math. Anal. Appl. 199, 502-525 (1996)

[22] Kuang, Yang, Delay Differential Equations with Applications to Population Dynamics. Academic Press, Boston 1993 [23] MacDonald, M., Time Lags in Biological Models. Lecture Notes in

Biomathematics 27, Springer Verlag 1978

[24] MacDonald, Biological Delay Systems. Cambridge University Press 1989

[25] Mackey, M., Glass, L., Oscillations and chaos in physiological control systems. Science 197, 287-289 (1977)

[26] Magal, P., Ruan, S., On integrated semigroups and age structured models in L^{p} spaces. Differential and Integral Equations 20, 197-243 (2007)

[27] Mallet-Paret, J., Nussbaum, R.D., Boundary layer phenomena for differential-delay equations with state dependent time lags I. Arch. Rat. Math. Anal. 120, 99-146 (1992)

[28] May, R.M., Stability and Complexity of Ecosytems. Princeton University Press 1974

[29] Maynard Smith, J., Models in Ecology. Cambridge University Press 1953

[30] Metz, J.A.J., Diekmann, O., Exact finite-dimensional representations for physiologically structured populations I. The abstract formulation of linear chaintrickery. In: Differ. Equ. with Appl. to Physics and Engineering, (Leibnitz 1989) Lect. Notes in Pure and Appl. Math. 139, 269-289, Dekker, N.Y. 1991.

[31] Nicholson, A.J., An outline of the dynamics of animal populations. J. of Zoology 2, 9-65 (1954)

[32] Nicholson, A.J., The selfadjustment of populations to change. Cold Spring Harbor Symp. on Quantitative Biol. 22, 153-173 (1957)

[33] Nussbaum, R.D., Periodic solutions of some nonlinear autonomous functional differential equations II. J. Diff. Equ. 14, 360-394 (1973)

[34] Nussbaum, R.D., A global bifurcation theorem with applications to functional differential equations. J. Funct. Anal. 19, 319-338 (1975)

[35] Perez, J.F., Malta, C.P., Coutinho, F.A.B., Qualitative analysis of oscillations in isolated populations of flies. J. Theor. Biol. 71, 505-514 (1978)

[36] Prüß, J., Stability anlysisi for equilibria in age-specific population dynamics. Nonlinear Anal. 7, 1291-1313 (1983)

[37] Smith, H.L., Reduction of structured population models to thresholdtype delay equations and functional differential equations: a case study. Math. Biosc. 113, 1-23 (1993)

[38] Smith, H.L., Equivalent dynamics for a structured population model and related functional differential equation. Rocky Mtn. J. Math. 25, 491-499 (1995)

[39] So, J.W.-H., Yu, J.S., Global attractivity and uniform persistence in Nicholson’s blowflies. Diff. Equ. Dyn. Syst. 2, 11-18 (1994)

[40] Sulsky, D., Vance, R.R., Newman, W.I., Time delays in age-structured populations. J. Theor. Biol. 141, 403-422 (1989)

[41] Walther, H.-O., The 2-dimensional attractor of x′(t) = −μx(t)+f(x(t− 1)). Memoirs of the AMS No. 544, Providence 1995

[42] Webb, G., Theory of Nonlinear Age-Dependent Population Dynamics. Marcel Dekker, New York 1985

[43] Wu, Jianhong, Theory and Applications of Partial Functional Differential Equations. Appl. Math. Series 119, Springer 1996

(Received August 16, 2007)