- Email: [email protected]
A population model with time-dependent delay
Mathl. Comput. Modelling Vol. 26, No. 12, pp. 13-30, 1997 Copyright@1997 Elsevier Science Ltd Printed in Great Britain. All rights reserved 08957177/97 $17.00 + 0.00
PII: SO8957177(97)00237-9
A Population Model with Time-Dependent Delay M. V. BARTUCCELLIAND S. A. GOURLEY* Department of Mathematical and Computing Sciences University of Surrey, Guildford, Surrey GU2 5XH, United Kingdom Qn.bartuccelli>~s.gourley>@mcs.surrey.ac.uk (Received July 1997; accepted Augwt 1997)
Abstract-we present analytical and computational results concerning the linear stability and instability of the uniform steady-state solution of a system of reaction-diffusion equations where a parameter in the kinetic terms is periodic in time. Under suitable sssumptions the system is equivalent to a scalar equation with a periodically varying delay. Such a varying delay can model the seasonal fluctuations to the regeneration time of a resource. We study the effect such a varying delay can have on the stability of the spatially uniform steady-state. Analytical results reveal that instability can set in if the delays are large, while computational methods of analysing the stability equations reveal the precise shape of the instability boundary. The nonlinear stability of the uniform state is also examined using ladder methods.
Keywords-Reaction-diffusion
equations, Time delay, Stability, Parameter domains.
1.
INTRODUCTION
Oneof the most commonly studied problems in developmental biology is that of pattern formation arising from an initially uniform state. This problem has been widely known in the literature as the Turing bifurcation since the seminal work of Turing in 1952 [l]. He studied reactiondiffusion systems and discovered the mechanism by which a uniform state can become unstable under the effect of diffusion. Since Turing’s discovery, many authors have applied and generalized his idea to a number of systems arising in ecology and embryology, and a good review is given in Murray [2]. I n recent years, many authors have extended Turing’s ideas by considering the effect of allowing the kinetic and the diffusion parameters to vary in space, in order to take into account spatial inhomogeneities. Ecological reaction-diffusion models which incorporate spatial variation in the kinetic terms were considered by Shigesada [3] and Cantrell and Cosner [4]. The effect of spatial variation in the diffusivities was considered by Maini et al. [5] and by Benson et al. [6,7]. Another extension, having more relevance to our work in the present paper is the possibility of parameters varying with time, and we shall briefly summarise those works that particularly motivated us here. Timm and Okubo [8] analysed the planktonic predator-prey model of Levin and Segel [9] with a temporal sinusoidal variation introduced into the diffusivity
of one of the species. Their particular type of temporal variation indicated a stabilising elfect and Gourley et al. [lo] proved a more general result to this effect for two-species reaction-diffusion models with small amplitude variations in the diffisivities of both species. Sherratt [ll] analysed *Author to whom all correspondence should be addressed.
13
14
M.V.
BARUCCELLI
AND S. A. GOURLEY
the case when the diffusion coefficients take the form of a periodic step function which alternates between two constant values but need not have small amplitude. In fact, he showed that large amplitude oscillations of this form can have a destabilising effect. This work was developed further in 1121. Though not having such direct relevance to the present paper, there are numerous other papers concentrating on the effects of periodic forcing in ecological modelling, for example, [13-151. Existence theory is treated in (16,171, while criteria for permanence (long term coexistence) are presented in 1181. Introducing time dependence into the parameters of population models allows us to incorporate daily, seasonal, or annual fluctuations in the values of those parameters and the importance of these effects is clear. However, in population models where delay is incorporated, for whatever reason, there naturally also arises the question of whether the delay itself should vary with time. For example, delay is sometimes incorporated to take into account the regeneration time of a resource, but there is no reason why such a regeneration time should always be the same-it is quite likely to be subject to seasonal fluctuations. In the present paper, we shall study the system it = Au + f&v), (1.1)
vt = ~Av + @(t)(u - v),
in which u is the population density of some species, v is a weighted average of past values of u (see below), 13(t) is a temporally periodic parameter related to the delay in the system (in a sense to be discussed shortly), and the constant ,u satisfies 0 < ~15 1. If 8 is constant or piecewise constant then the second equation can be solved to give v in terms of u: v=
t +CQexp ss
-((’ - ~““‘*~” - ‘))gexp
--oo -co
[4np(t - s)]1’2
[-@(t
-
s)]u(y,
s) dy &
(1.2)
justifying our assertion that v is delayed 2~. For more general periodic functions f?(t), (1.2) is only an approximation but it becomes increasingly accurate when the period of 8 is large. The type of application we have in mind involves the period being large so this is not too much of a drawback. The model we are studying is thus a single species reaction-diffusion population model in which f(u, v) depends on the current density ~(5, t) and also on the spatio-temporal convolution v, given by (1.2), which takes account in a particular way of past densities and nonlocal spatial effects. It follows closely the type of model introduced by Britton 1191and Gourley and Britton [20] who argued that whenever delay effects are incorporated (for whatever reason) into a diffusion model, it is necessary to take into account the fact that because of the diffusion the individuals will not have been at the same point in space at times in the past. To account for this a detailed derivation is provided in [19], using random walk arguments, for the use of the double convolution (1.2) (with @ = 1) as an expression for the lagged density, and in that paper it is also observed that such an expression satisfies a differenti~ equation, namely, the second equation of (1.1) (with p = 1). In the present paper, we will in fact allow ,Qto be any constant satisfying 0 5 p < 1. We then have a choice of either taking ,u = 1 or we can take the limiting case /I --) 0, in which case expression (1.2) can easily be seen to degenerate to a purely temporal convolution of the kind frequently seen in ordinary differential equation models with delay 1211. Turning now to the significance of B, note from (1.2) that if 0 is large the delay in the system is weak, in the sense that although all values of u in the past contribute to the spat&temporal average v, only recent values can have a strong influence. Decreasing the value of 6 increases the influence given to values of u far in the past. It is quite reasonable, therefore, to think of the quantity l/e as being a quantitative measure of the delay in the system. In the present paper, we are allowing 8 to be a periodic function of t, alternating between two constant values. Thus,
Population
Model
15
the delay will do likewise and we can regard this as a primitive means of modelling situations where the delay in the system is a periodic function of time. For example, delay is often used to model the regeneration time of resources and this is quite likely to be longer in winter than in summer. Numerous other possible reasons for incorporating delay in population models are described in [21]. As far as time-varying delays are concerned, the use of a piecewise constant periodic delay is perhaps the simplest case that allows progress to be made analytically. As the notation implies, f(~, v) allows dependence on both current density u and lagged density v given by (1.2). Any function f that is biologically realistic and sufficiently smooth is allowable. An obvious possibility is the logistic case f(u, v) = ~(1 - v). An interesting generalisation of this is given by f(~, v) = ~(1 + a~ - (1 + a)~), where a is an aggregation parameter. This particular case will be considered in Section 2.2. In Section 2, we present the linearised stability analysis for the spatially uniform steady-state solution of our model. Section 3 then addresses the question of nonlinear stability using ladder theory. Then, in Section 4, we study the linearised stability criteria numerically to augment our analytical results and finally, in Section 5, we outline what conclusions may be drawn.
2. LINEAR
STABILITY
ANALYSIS
We assume that the system (1.1) has a spatially uniform steady-state solution (u, v) E (u*, v*) such that (u*, v*) > 0. We are interested in analyzing the linear stability of this uniform state. Notice that the second equation of the system implies u* = v’, so the steady state satisfies f(~*, u*) = 0. Setting 21 = u* + ii, v = 2L*+ 6, substituting
into (1.1) and then dropping the tildes, we obtain the linearised system ut = AU + au + bv,
(2.1)
vt = j.~Av+ B(t)(u - v).
In (2.1), a and b are the partial derivatives of f(u,v) with respect to the variables u and v, respectively, evaluated at the uniform state. The limiting case of zero delay is found by letting 0 --) 00 and from (1.2) it is easily seen that in this limit v + u (see, also, [20]). Then the first equation of (2.1) becomes ut = Au f (a + b)u.
(2.2)
Hence, if a + b < 0, the uniform state is stable, with and without diffusion. We will henceforth assume that a + b < 0, since we are interested in instability for the full system that is delay induced, rather than instability that is present even without delay. We shall now analyse the stability of the uniform state for the case when 8 is constant. Then we shall investigate a particular time-dependent 8 in the form of a periodic step function. Seeking solutions of (2.1) (with ~9constant) of the form (u, v) = (cl, ~2)exp (st + ik . x) , we find that
s-a+K
-8
,+0&K)
(::)
= (:)’
where K = k2 = k . k. For nontrivial solutions, we need the determinant zero, so s2+(tl+pK-a+K)s+(B+pK)(-a+K)-M=O.
of the matrix to be (2.3)
16
M. V. BARUCCELLI AND S.A. GOURLEY
Instability of the uniform state can set in either through a steady-state bifurcation (s = 0), or a Hopf bifurcation. For a given wavenumber k, a steady-state bifurcation may occur when (4 + /.&K)(-a -I- K) - be = 0,
(2.4
0 = PK(-a + K) (a -I-b) - K ’
(2.5)
i.e., when
Since a + b < 0 and B > 0, we see imm~iately that steady-state bi~rcatio~ will not occur if a < 0. To seek a Hopf bifurcation we set s = iw, with w real and positive. Such a ~furcation may occur when 0-a+(p+l)K=O, (2.6) provided that (Q $ pK)(-a
+ K) - b@ > 0.
(2.7)
Hence, again one can see that if a < 0 there are no Hopf bifurcations. Thus, if a c 0, the uniform state is stable regardless of the value of 0. Now let a + b < 0, but a > 0. It is convenient here to think of l/e as a bifurcation parameter, since it represents delay. Small delay usually implies stability and, in this particular case, graphical arguments help to show that the uniform state is stable to all su~ciently small ~~urbations when
(2.8) As I/0 is increased instability may set in either through a Hopf bifurcation (when 8 = a), or a steady-state bifurcation, depending on the values of the other parameters (for example, if ,G is sufficiently small we expect a Hopf bi~rcation). Further technical conditions need to be checked to be certain that bifurcation actually occurs and a useful reference here is the article by Crandall and Rabinowitz [22). The bifurcating solutions could be constructed using methods similar to those in [20], in which the noncompactness of the spatial domain is dealt with by restricting attention to spatially periodic solutions. In the case of a steady-state bifurcation, the critical wa~n~ber at the bi~~ation can be cumulated and is the value of K E (0, a) at which the relevant minimum in (2.8) is attained. This completes the analysis for the case when 6’ is constant and, as we have seen, it is quite straightforward to enumerate the various possibilities. We are now going to introduce temporal variation in 19and study the case when e(t) is a periodic step function with period T, so that
e(t) =
(2.9)
where n is any integer. This choice of e(t) follows the work of Sherratt Ill] and will enable us to make some analytical progress towards determining conditions for stability/instability of the uniform state. Introducing e(t) defined above into system (2.1), we obtain ut = Au + au + bu, vt = PAW + e(t)(~ - vf.
(2.10)
To study stability of the uniform state to perturbations of wave vector k, we set (u,v) (U(t), V(t)) exp(ik . x) and substitute into (2.10) to obtain u, = -KU v, = -pKV
=
-t aU + bV,
+ e(t)(u - V),
(2.11)
Population Model
17
where K = k2 = k a k. Denoting w = (V, V)T, we can write this in the form
dw=Afw
’
dt
d”=A-w dt
’
where
The method for solving systems like (2.11) is Floquet theory, and we summarise this theory very briefly. For more details see, for example, [23]. Floquet theory applies to linear ordinary differential equations with periodic coefficients. The solutions of this class of equation will generally be a periodic function times an exponential. The general technique is to solve the system of linear differential equations for any set of complete linearly independent solutions, for 0 5 t 5 T. This yields a so-called propagator matrix E such that the general solution can be expressed as ~(t + mT) = E”“y(t), where y(t) is the general solution of the system of order n, andm=1,2,3 ,.... Note that the propagator E can be expressed in terms of any fundamental matrix Q(t) as E = @-l(O)@(T). The eigenvalues of E are known as the Floquet multipliers of the system. We shall denote them by pi (i = 1,2,. . . , n) for a system of order n. It can be shown that the characteristic multipliers do not depend on the particular choice of the fundamental matrix a(T), and therefore, are intrinsic properties of the system. Furthermore, to determine the stability of the solutions of the problem, one needs to compute the eigenvalues of the propagator matrix E. From y(t + mT) = E”‘y(t), it follows that the uniform state is stable if and only if all the eigenvalues of E have modulus less than one. Hence, if any of the eigenvalues is larger than one in magnitude then the solution is unstable. We shall use this theory to determine the Floquet multipliers, which in general are complex numbers, and to give information about the stability of the solutions of (2.11). Following the notation used in [ll], we denote the eigenvalues of A* by Xt2 and the eigenvectors by z$. We further define Zf to be the matrix whose columns are the eigenvectors of A*. After some algebra, we find that
zf,
Xf=i(a-B*-(p+l)K+P’), (2.12) X:=l(a-@-(p+l)K-P’),
and that
-b
4
;(Qi-P*)
’
;(Q&+Pk)
where l’* = 44b6’* + (Q*)2,
(2.13)
Q*=a+8*+(~-l)K. Defining n*(t)
=
0
exp(Xft)
(
o
exp(Xtt)
)
’
it follows that any fundamental matrices of (2.11) have the form Z+A+(t)C+ and Z-A-(t)C-, We take C+ = where C+ and C- are matrices whose entries are constants of integration. (A+(T/2))-l. Because we want continuity at t = T/2, we require that C- = (A-(T/2))-1 (Z-)-‘Z+. Hence, we obtain E = (@(O))-%(T)
= (2+)-l
T
z
Z-A0
(Z-)-‘Z+A+
0
;
M. V. BARUCCELLIAND S. A. GOURLEY
18
The last matrix A+(T/2) in the above actually appears at the front in the computation of E, but we have written it at the back for convenience in the subsequent calculation. Matrix multiplication is not, of course, commutative, but what is important are the e~ge~vu4~esof E and we are using the result that if A and B are any square matrices, then AB and BA have the same eigenvalues. The uniform state will be stable (to perturbations of wavenumber k) if and only if the two eigenvalues of E have modulus less than 1. After some algebra, we can see that the matrix E has the form E = exp (-lY’/4)fi‘,
where (2.14)
r=-2(X,+X;t)+~++P-=8++8---2[a-(~+l)K] and B=exp
(P+ + P-)T [
4
x
1
[(z+W (: exp(+))Pm’z+(; _(~))].
It is straightforward to see that det(@) = 1. Hence, the characteristic equation of B has the form (2.15)
&++1=0, where B, which is the trace of k is given by i3={(l-exp(e))
2b (e+ + e-) + Q+Q-
(
(I-exp(F))
+ (l+exp(T))
P+P-
>
(2.16)
(l+exp(~))}expi(P+~P-)T’41,
with Q*, P* given by (2.13). Analytical results are difficult to obtain in general, but it is possible to make progress by considering various limiting cases. We start by considering the case when Bf and 8- are very close together. Then we see what can be said if they are widely different. We then treat the case when 9+ and 13~are both very large, and after that the case when both are very small. These two cases correspond to small and large delays, respectively. Then we shall prove that if a + b < 0, a < 0, and b > 0, the uniform state is stable for all values of Bf and 8-. Note that the condition for stabidity of the uniform state, as a solution of the reaction-diffusion system is that for every K>O,weneed lexp(F)
Q/Cl,
(2.17)
where fi satisfies (2.15). Note also that the coefficient B in (2.15) need not be a real number, since one or both of P+ or P- can be purely imaginary. 2.1. Small Amplitude
Periodic
Variation
To keep the algebra as simple as possible we will, in this section and the next, consider stability of the uniform state to spatially uniform (K = 0) perturbation only. This is equivalent to studying its stability as a solution of the corresponding ordinary differential equation system, i.e., system (1.1) without diffusion terms. We assume, as always, that a + b < 0. Recall that in the case of constant 8, instability of the uniform state is possible if a > 0 and, in this case, that (2.8) gives the value of 8 at which instability actually sets in. We will assume in this section that a > 0 (note this implies b < 0). Since K = 0, the instability sets in through a Hopf bifurcation when 9 = a, with stability when 0 > a.
Population Model
19
Turning now to the case of varying 8, our aim is to get some insight into how this periodic variation changes things. Analytical progress can be made by considering the effect of a small amplitude periodic pe~urbation whose mean value is a, so we set 8+ = a + E and B- = a - E, whereO
with a similar expression for P”, and that
’ sin’ fT,“i;;T;L”-“;jr o2(a + b)
TbmsinTdM~&&~(u+ b)
.
(2.18)
Since I’ = 0, condition (2.17) becomes Ifi] < 1, where fi2 - Bb f 1 = 0. For “most” values of the period T of B(t), the leading order term in the above expression is less than two in absolute value and this implies, given s is small, that also IEI 5 2. The roots 3 are then a complex conjugate pair on the unit circle. So this implies the system usually remains marginally stable under the effect of a small amplitude periodic perturbation. The values Tn:=: d&,
n=l,2,3
(2.19)
,..,,
need to be considered in more detail because the leading order term in (2.18) then has modulus 2. Note that these values are closely related to the period of the Hopf bifurcating solution in the constant B case, which can be shown to be
The cases n even and n odd have to be looked at separately as follows. CASE 1: n ODD. In this case, cosT,dm Tn dj
= 0, and sin2(1/2)
= -1, sinT,dm
= 1. Thus,
BZ-2-E2
b <-2, u2(a + b)
since a + b < 0 and b < 0.
The roots fi are therefore real, and are located just on either side of -1. Thus, the uniform state becomes unstable. CASE 2: n EVEN. In this case, the 0(s2) term in (2.18) is zero, so higher order computations are needed. We have done some numerical simulations, however, and these suggest the uniform state is unstable, but much more weakly so than in the case of n odd. 2.2. Large
Amplitude
Periodic
Variation
In this section, we again assume K = 0 and study stability to spatially uniform perturbations only. What we will do is show that under certain conditions a large ~plitude periodic variation in 6(t) may have a stabilising effect. We shall illustrate this phenomenon with reference to the population equation
t s -00
6e-e@-S)u(s)
ds
, >
20
M. V. BARUCCELLIAND S. A. GOURLEY
with a, 8 > 0 and its equivalent system dtt z
=%6(1+at4 -(Ifo)u), (2.21)
dv
-& = L9(u- v). This system (in the constant 6 case) is one of numerous special cases of a more general equation studied in [19j. The parameter a measures the tendency of the species to aggregate or group together as a defense against predators. The quantity a used elsewhere in this paper as a coefficient of the linesrised system is equal to the a of this particular model. The linear&&ion of (2.21) about (u, v) = (1,l) is in fact dU z
= aU
-
(1 i- u)u, (2.22)
dv -& = @(u-v),
so, in this model, b = -(l -t a). In the constant 0 case, a Hopf bifurcation occurs at 8 = a. We now introduce variation in 8 according to (2.9) as usual, but not necessarily with mean value a this time. If 8’ is large and 8- is small, say, this guarantees Pi and P- are both real (and also that I? > 0). Now let the period 2’ -+ 00, then B
2
O(1) exp
(P+ + P-) T 4 ( > s
The two roots of pz - Bfi + 1 = 0 are, therefore, (~ymptotic~ly) B and l/B. When ji = l/B, it is quite obvious that 1exp (-I’T/4)$I < 1. Fbr the root fi = B, the stability condition becomes ]O(l)exp(~)exp[(P++~-)T]/cl,
forTlarge,
or I‘ > P+ + P-. To simpli~ this, we shall exploit the assumed largeness of 8+ by using an asymptotic expression available from the next section (expression (2.23)). Also, if we let 8- 40, then P- + a. The stability condition reduces down to a < 1, which is perfectly possible biologically and means that the tendency of the species to aggregate is not too large. Since large 0 means small delay and vice versa, what we have done in this section is to demonstrate that if one delay is very small, the other is large and also the period is very large then stabilisation is a possibility, depending on the values of the other parameters. Actually this result did not surprise us. The differential equation y
=
u(t) (1 - u (t - o(t))),
with a(t) =
has recently been studied by Schley and Gourley [24]. In this equation, the delay one, alternating between zero and 2~s. We have shown that if crs is related to the T = 4~70,the constant state u E 1 is stable for any 00 > 0, and in some respects this the one of the present section. For smaller periods, however, the stability depends of 00.
is a discrete period T by result is like on the value
Population Model
21
2.3. The Case 19+and 8” Large We now return to the general linearised system of the paper, system (2.11) with (2.9), and we allow perturbations of arbitrary wavenumber so that K 2 0. When 8+ and 6- are large we can show that P+-8++(2b+a+&-l)K)+O
& (
with a similar result for P-.
(2.23)
, >
Other relevant facts, when 0+ and 8- are large and that r > 0,
P+, and P- (and hence, B) are real, and B N exp[(P+ + P-)T/4] >> 1. Again the two roots of fi2 - Bji + 1 = 0 are asymptotically B and l/B. Much as in the previous section, the stability condition reduces to
r>p++p-. From (2.14) and (2.23), it follows that we require 4(a + b) - 4K < 0 for all K 2 0, and this is true because a + b < 0 by hypothesis. Therefore, for 8+ and 0- large (i.e., when the delays are small), we have stability of the uniform state. Small delay usually does imply stability, and in this section, we have shown this to be the case when a periodic delay oscillates between two small values.
2.4. The Case 8+ and 8- Small When 8+ and 8- are small, we can see that 2b
P*-la+(p-l)KI+
+
a + (cl - 1)K
Ia + (p
_
l)KI
‘* +’ (e*2) ’
>
(2.24)
and that Bpexp
1 [
(P++P-V 4
-(P++P-)T
+exp
4
I
= 2 cash The roots of fi2 - Bji + 1 = 0 are, therefore,
1*
f(P++P-)T
ji = exp [
4
The condition for stability is exp[-I’T/4] p < 1. This requires -I? f (P+ + P-) and P- are positive, l? > P+ + P-. Using (2.14) and (2.24), this gives 0++6--2[a++l)K]
>2la+(p-l)KI+
2bl+;;($!,K)
< 0 or, since P+
(0++8-).
(2.25)
It is convenient to consider the two cases a + (II - l)K < 0 and a + (p - l)K > 0 separately. CASE 1: a -t (/.J- l)K < 0. In this case, (2.25) reads, after some algebra
(e+ + e-)
a+b+b---1)K a + (p - 1)K
>
_zpK
1
Since a + (cl - l)K < 0, a + b < 0, and 0 5 p 5 1, the above inequality is trivially satisfied, since the square bracketed term is positive. CASE 2: a + (p - l)K < 0. Note that for this to hold it is necessary that a > 0, and this in turn implies b < 0 (since a + b < 0). If these conditions are satisfied then a + (/J - l)K is positive for all wavenumbers K if /J = 1, and if p < 1 it is positive for wavenumbers satisfying
M. V. BAAUCCELIJ ANDS. A. GOURLEY
22
K E [O,(u/l - CL)).It turns out in this case that instability is possible under a further condition which we shall obtain. Reversing the inequality in (2.25), we obtain the condition for instability that e++e--2(a-(p+l)K) or after
2b + a + (fi - 1)K a + (/J - l)K
<2(a+(/J-l)K)$
(e+ + 8-) ,
(2.26)
some algebra that Q(K) := 2(K - a) (a + (p - 1)K) < b (8+ +
@-),
(2.27)
for some K E [O,(a/l - p)), or [O,00) if p = 1. Note that the right-hand side of (2.27) is negative, since b < 0. In both of the two cases /.J= 1 and p< 1, simple graphical arguments will show that we require Q(0) < b(B+ + P), and this holds whenever
e++ewhich we can regard as automatically
2
< $fq,
satisfied, by the smallness of 8*.
Thus, if 8+ and B- are small (i.e., if the delays are large), then our conclusions are as follows. If a < 0, the uniform state is stable (because Case 1 applies for every wavenumber). If a > 0, then Case 2 applies, either for every wavenumber or for a finite range of wavenumbe~, and in either case there are wavenumbers which give instability of the uniform state. 2.5. TheCasea+bO We are going to prove that under the above assumptions, the uniform state is stable (to arbitrary small perturbations) for all positive values of 8+ and 8-. Our particular method of proof actually depends on the coefficient B in the equation fi2 - Sfi + 1 = 0 being real (for all K 2 0). The condition b 2 0, among other things, ensures this. We would actually conjecture, though, that the uniform state is in fact stable under the conditions a + b < 0 and a < 0 only. First, note that b 2 0 implies, from (2.13), that Pf and P- (and hence, B) are real. The inequality of Lemma 1 (below) helps to show, among other things, that in fact B > 0.Now,the roots of ,G2- Bfi f 1 = 0 may be either real or complex. The latter p~sibility is trivial to deal with since the roots are then a complex conjugate pair fir and &, and their product is fi& = 1. Hence, [,&ii= 1 and so trivially 1exp (-I’T/4)fiJ < 1, since a < 0 implies r > 0. The case when the roots of b2 - B/i+ 1 = 0 are real is more awkward. First, note that if the roots are real then they are both positive, since B > 0.This means that for stability it is enough to show that each root is less than exp @‘T/4). A simple graphical ~gument shows that if we set f(P) = p2 - B@ + 1 then it is sufficient to show that the inequalities (2.28)
G fi=exp@‘T/4)
(2.29)
> 0,
both hold. To proceed we now need the following lemma. LEMMA
1. Ifa+b
then 2b (e+ I
+ e-1 + Q+Qpip-
PROOF. After some algebraic manipulation, we see
b(a+b+(~-l)K)IO, which is clearly true under the stated assumptions.
< I
-
I ’
that the above inequality is equivalent to
Population Model
23
In order to prove (2.28), we proceed as follows. Now ~(exp(~))
=exp(a>
-~exp(~)
+l,
(2.30)
and B={(l--exp(T))
(I-exp(T))
+ (1 +exp (q))
(1 +exp
(2b(8f+~~J+Q+Q-) (l$2))}
exp[(P+
~P-)*~41,
Using Lemma 1, B 5 f exp (P++:S)T]
{(l-exp(T))
[
+ (l.+exp(T))
(l-exp(T))
(l+exp(F))}.
(232) .
We point out here that Lemma 1 can also be used to confirm that B > 0.Now let us introduce ,
We find, using (2.32), that f (exp (3))
2 z2 - i XYZ[(1+;)(1++(1-$)(I-$)I+1 =,r2--2 =
(
(z -xy)
Cry+;
)
+1
(*-$ ><
The second bracket in the last line is trivially positive, since z,y > 1 and z > 1. The first will be positive if t > xy, that is, if exp @‘T/4) > exp (P+ + P-)T/4 or I’ > P+ + P-. This is equivalent to e++e--z[a-(/.A+1)K]
> da+
and a sufficient condition to ensure this is that the inequalities 6+ - (a - (p + l)K) > 44bB+ + (Q+)2, e- -(a--(p+l)K)
> da,
should both hold, which they do if and only if 4(a + b)8f + 4K (a/A- e* - /AK) < 0. The latter is certainly true, because a < 0 and a + b < 0. This completes the proof that
f (em (rT/4)) > 0.
24
M. V. BARUCCELLIAND S. A. GOURLEY
It remains to show that %I P=exp(r~,l~ > 0. Using the same variables z, y, z as above, we get
>2z-,,(I+--&),
by(2.32) but .z > zy
>XY--=
1
X2Y2 - 1
XY
XY
2 0,
as desired. For cases when B can be complex there is not much we can do analytically and we have to resort to analysing the stability criteria by numerical inv~tigation of a parameter space of interest. The results of some such numerical investigations are presented in Section 4. When the result of the present section is applied to the particular model studied in Section 2.2 (equation (2.20)), the conclusion is that if a < -1, then the steady-state u z 1 of that model is stable for all values of 8” and I!?“. Britton [19] proved stability in the case of constant 8 under the condition a < -l/2. In the next section, we shall consider the question of the nonlinear stability of the uniform steady-state solution of the system (1.1).
3. NONLINEAR
STABILITY
ANALYSIS
In this section, we exploit some of the ladder methods developed by Bartuccelli et al. [25281 to prove some nonlinear stability results for the original system, which we restate here for convenience: it = Au + f(u, w), (3.1) vt = ~Av + O(t)(u - v). In order to make the ladder methods available we shall assume, in this section, that we are working on the finite spatial domain x f [O,Lid, with periodic boundary conditions, where d is the spatial dimension. Also, for the time being, O(t) is any positive differentiable function and not necessarily periodic. The analysis outlined below will hold for any biologically meaningful nonlinear term, but to keep the algebra simple we restrict ourselves to the particular case vzt = Aa+u(l
-cru-flu),
2rt = ~Av + O(t)(u - v),
(3.2)
where CYand p are positive constants. We are interested in analyzing the nonlinear stability of the uniform state. To use the ladder methods effectively, it is very desirable to centre the system on the uniform steady&ate (2t*, u*), so we again introduce u =
us+ ii,
v =
u* + 6,
where in this system u* = l/(cy + p), substituting into (3.2) and dropping the tildes (but not the higher-order terms this time), we obtain the nonlinear system ut=Au-azb-bv-~~2-/3~v, vt = /.JAv + O(t)(u - v),
(3.3)
where a = &/(a! -I-@) and b = /3/(a + /3). Note that a + b = 1, and so we are in the same hypotheses of the linear analysis of the previous section, since a and b in (3.3) are equivalent to
Population Model
25
-a and -b in the linear system (2.1). The u and v in (3.3) are not positive functions, of course, since they are really fi and 6. However, solutions of the original system (3.2) are undoubtedly positive if the initial data is, by the maximum principle. In view of this, we can say that solutions of (3.3) satisfy u* + % > 0 and u* + v 2 0. This will be very important later. We introduce the positive definite quantities Hn = 8
J
(Vn~)2 dx,
Jn = b
J
(V’%Q2 dx,
Fn=H,+J,.
(3.4)
If we differentiate H, with respect to time and use (3.3), and then integrate by parts on the Laplacian term (the integrated terms vanish by virtue of the periodic boundary conditions), we obtain iI!& = -H,,+l
- aH, - be / (V%) (V’%) dx J
Y
-
a~J
(vnu2) dx - pe
(vu)
J (vu)(VW)dz + f
(3.5) j+.
Doing the same calculation for J, gives ;j_
=
J
-pJ,+l + b8
(V%) (V”u) dx - SJ,.
(3.6)
Adding the & and the JR, we get a differential equality (ladder) for the F,. Now it turns out that for investigating stability of the uniform state it is not necessary to analyse the ladder for arbitrary nonnegative integers 7~2 0, but it is sufficient to study the bottom rung of the ladder, which is given by n = 0. We shall then obtain what it is known ss nonlinear stability in the L2 norm. Adding the H, and the J,, differential equalities for n = 0 gives (3.7) The above first-order nonlinear ordinary differenti~ equation contains all the information we need to analyse the stability of the zero solution of (3.3), and hence, of the uniform state (u*, u*) of the original system (3.2). Naturally the stability will depend, as in the linear stability analysis already done, on the functional expression of the kinetic parameter 19(t). First, we note that the fifth and sixth terms in the right-hand side of (3.7) are in general of indefinite sign. But, as previously noted, we know that U* + u I 0 and U* + v > 0. This gives --u < U* and -v < u*, and we use this information to obtain the estimates
J u3
-CUB
dx 5 CXU* Ho
and
-PO
J
u2vdx I /3u*Ho.
Employing these estimates in (3.7) yields a differential inequality
;po <-qJ1 Note that the HO term accounts a “linearised” inequality
-HI
18 - aHo - OJ, + Ho + z sHr,.
for the nonlinear terms in (3.3). By dropping it we can obtain
;po-pJ1 5
- HI - aHo - 6Jo + f ;Ho,
(W
from which results can be found for the linearised system. Let us now list various consequences of both the full and the linear&d inequalities (3.8) and (3.9). The idea is to try to reconstruct Fo in the right-hand side of (3.8), and then find
M. V. BARUCCELLIAND S. A. GOURLEY
26
conditions on 19(t) and the other parameters which will guarantee that FO + 0 as t ---f 00. This will give nonlinear stability, in the L2 norm, of the zero solution of (3.3). First, note that, by Poincare’s inequality, Hi 2 I~-~I-&,, and Jr 2 K2Jo,. Hence, (3.8) becomes
“@ <2
Ho -
O-
($ +3)Jo 5 -&3)Fo,
where 4(t) = min
$+a-1
-$$),
(S+@(t))].
K Hence, if (3.10)
linl~f 4(t) > 0, then FO --f 0 as t -+ 00. There are several consequences of this in the following.
(1) Provided
B and 8 are bounded for all t > 0, we can ensure (3.10) holds by taking L sufficiently small. Thus, if L is sufficiently small, the uniform state is stable. This is in accordance with the general observation in reaction diffusion systems that spatial patterns cannot emerge if the spatial domain is too small. (2) 8(t), a nonincreasing function of time (in particular, constant). Under this assumption the term (I/2)(~/~)~ e is nonp~itive, and therefore, from (3.10), it follows that Fe(t) + 0 as t ---f cm provided that min K
;+a-1=;-b)‘(s+e(t))]
>o.
(3.11)
(3) B(t)
large. In this csse, ~1/2)(~/~)~0 is small and so & < 0, provided again (3.11) holds. Therefore, we have stability also for large values of 8, i.e., for small values of the delay.
Notice that the argument in Cases (l)-(3) a b ove, also holds for nonperiodic functions 8 = O(t), imposed.
so long they satisfy the restrictions
(4) Assume now that e(t) is a periodic function of time, and sssume its period is large; let us introduce the new time variable i = t/T < 1; this will enable us to analyse, in particular, the stability of our uniform state for large values of the period of O(t). The Fo ladder now becomes 14 f&j L 2’ (-/A& - Hi - aHi, - SJ, + Ho) + 2 5Ho. (3.12) Therefore, if 2’ is sufllciently large the (1/2)(~/~)Ho function 4(t) defined above becomes
term can be neglected,
and the
$*(t)=min[(++a_1).($+B(t))]. Therefore, provided (3.11) (5) If the period of our kinetic This can be seen from (4) and so Fo increase in time,
is satisfied, we again obtain that the uniform state is stable. parameter B(t) is small we can show that we have instability. above: when T is small the 6/0He term is the dominant one that is, we have instability.
This concludes the nonlinear analytical analysis of the system.
Population
(a) p=
Model
1, a = 4, b = -5,
27
T = 4.
4.0
3.0
z-
2.0
1.0
1.0
2.0
(b) /A = 1, a = 1, b = -2,
3.0
4.0
T = 4.
Figure 1. Diagram showing typical regions of stability and instability in the paramfor the zero solution of the eter plane (r+, T-), where r+ = l/e+, T- = l/e-, linearised system (2.10) with 0(t) defined by (2.9). In the light regions, the zero solution is stable to perturbations of arbitrary wavenumber while in the dark regions there is at least one unstable mode.
4. NUMERICAL
COMPUTATIONS
It has become apparent by now, that in order to make analytical progress with the type of system we are studying in this paper, it is almost invariably necessary to consider various limiting
M. V. BARUCCELLI AND S.
28
A. GOURLEY
5.0
4.0
3.0
z2.0
1.0
0.0
I
=
1, a = 0.5,b = -0.52,
P =
a = 0.5, b = -0.52,
(c) P
T = 4.
5.0
4.0
3.0 2-
2.0
1.0
0.0 2+
(4
0,
‘1’ = 4.
Figure 1. (cont.)
csses and use ~ymptotics. The only exception to this was in Section 2.5, but the gumptions made there are difficult to interpret and the methods used in that section are unlikeIy to generate the best possible result. This leaves the question of what the stability regions really look like in parameter space. We wrote a simple computer program to determine regions of stability and instability in the parameter space (T+,T-), where T+ = l/6+ and T- = l/19-, so the r’s may be interpreted as
Population Model
29
delays. At each point on a grid of parameter values, the program simply checks the stability condition 1exp (-P/4) fi] < 1 for each of a large number of (equally spaced) wavenumber values up to a suitably large wavenumber. If instability is found for any of these wavenumbers, the point is coloured black. If all wavenumbers have been checked and no instability found, the point is left white and the computer moves to the next point on the grid. Some stability diagrams generated by this simple code are shown in Figure 1. The diagrams are all symmetric about the line T+ = T- as we should expect. In all four diagrams, one can see that we have stability for small values of the parameters r+ and T-, as expected (small delays give stability). Two of the diagrams also bring out rather clearly the phenomenon noted in Section 2.2, for one specific system, that if one delay is very large and the other very small, then the result can be stability. Our diagrams suggest that this phenomenon may not just be a peculiarity of that particular system and also that it may occur even when the period T is not particularly large. Recall that in Section 2.2, we could only demonstrate it analytically for very large values of T. Various other experiments suggest a general result that if a > 0, there is always an instability region in the parameter plane we have considered. We were not able to get instability for any negative value of a, we looked at but have no analytical result to this effect. Recall that to establish the stability result of Section 2.5 it was necessary to assume not only that a + b < 0 and a < 0, but also that b > 0. The latter assumption was necessary for our particular method of proof to work but we believe the result may hold true even without this assumption. If this is indeed the case it means that in the a + b < 0 and a < 0 situation the periodic variation does not change stability, since these conditions certainly do guarantee stability in the constant 19case. 5. CONCLUSION In this paper, we have described a simple means of incorporating a time-dependent delay into a reaction-diffusion population model, while still keeping the model amenable to analytic treatment. The periodic delay can be regarded as a means of modelling periodic fluctuations to the regeneration time of a resource and the periodic step function is the simplest case which can be dealt with analytically. We have analysed how the stability of the spatially uniform steadystate is affected by the temporal variation of the parameter 13,which is the reciprocal of the delay in the system. We have found that if a > 0, then a small amplitude periodic variation in e(t) will be destabilising if a resonance effect occurs, that is if the period of e(t) matches that of the Hopf bifurcating periodic solution that we would expect to see if t9 were constant. We have also seen that if conditions are right, it is possible for a large amplitude variation to be stabilising. However, it would be imprudent to draw too general a conclusion on this latter point from what we have done here and further investigation of the large amplitude case is essential. The observations made above would seem to be at variance with those of Timm and Okubo [8], Gourley et al. [lo), and Sherratt [ll], all of whose results suggest that small amplitude perturbations are stabilising while larger ones may be destabilising. In all of these papers, however, it was the diffusion coefficients, and not the kinetic parameters, that were oscillating in time. Furthermore, these authors were considering the effect of temporally periodic perturbations to systems which in the unperturbed case are at bifurcation to a nonhomogeneous steady solution, rather than systems undergoing a Hopf bifurcation to a solution periodic in time. We have also found that if 8+ and C9- are both large then the uniform state is stable, as we expected since this is the small delay case. If they are small the uniform state may be unstable, depending on the values of the other parameters. Furthermore, we obtained sufficient conditions for stability which are independent of 8+ and 8- . Numerically computed stability diagrams confirm the analytical analysis as well as revealing the precise shape of the stability boundary. Finally, we performed a nonlinear stability analysis of the uniform state of a particular system by applying ladder methods. The use of these methods has shown that, as expected, the stability or
M. V. BARIJCCELLIAND S. A. GOURLEY
30
otherwise of the uniform state depends on the size of the spatial domain. The ladder approach also seems to work well for a much broader class of time-dependent functions O(t) than does the linearised
analysis.
REFERENCES 1. A.M. Turing, The chemical basis of morphogenesis, Phil. Runs. R. Sot. Lon. B237, 37-72, (1952). 2. J.D. Murray, Mathematical Biology, Springer, Berlin, (1989). 3. N. Shigesada, Spatial distribution of rapidly dispersing animals in heterogeneous environments, In Lecture Notes in Biomathematics, Volume 54, (Edited by S.A. Levin and T.G. Hallam), pp. 478-491, Springer, Berlin, (1984). 4. R.S. Cantrell and C. Cosner, The effects of spatial heterogeneity in population dynamics, J. Math. Biol. 29, 315-338, (1991). 5. P.K. Maini, D.L. Benson and J.A. Sherratt, Pattern formation in reaction-diffusion models with spatially inhomogeneous diffusion coefficients, IMA J. Math. Appl. Med. Biol. 9, 197-213, (1992). 6. D.L. Benson, J.A. Sherratt and P.K. Maini, Diffusion driven instability in an inhomogeneous domain, Bull. Math. Biol. 55, 365-384, (1993). 7. D.L. Benson, P.K. Maini and J.A. Sherratt, Pattern formation in reaction-diffusion models with spatially inhomogeneous diffusion coefficients, Math. Comp. Modelling 17 (12), 29-34, (1993). 8. U. Timm and A. Okubo, Diffusion-driven instability in a predator-prey system with time varying diffusivities, J. Math. Biol. 31, 307-320, (1992). 9. S.A. Levin and L.A. Segel, Hypothesis for origin of planktonic patchiness, Nature 259, 259, (1976). 10. S.A. Gourley, N.F. Britton, M.A.J. Chaplain and H.M. Byrne, Mechanisms for stabilisation and destabilisation of systems of reaction-diffusion equations, J. Math. Biol. 34, 857-877, (1996). 11. J.A. Sherratt, Turing bifurcations with a temporally varying diffusion coefficient, J. Math. ,Eiol. 33, 295-308, (1995). 12. J.A. Sherratt, Diffusion-driven instability in oscillating environments, Euro. J. Appl. Math. 6, 355-372, (1995). 13. P. Hess and H. Weinberger, Convergence to spatial-temporal clines in the Fisher equation with time periodic fitness, J. Math. Biol. 28, 83-98, (1990). 14. G.J. Butler and H.I. Freedman, Periodic solutions of a predator-prey system with periodic coefficients, Math. Biosci. 55, 27-38, (1981). 15. J.M. Cushing, Oscillatory population growth in periodic environments, Theor. Pop. Biol. 30, 289-308, (1986). 16. P. Hess, Periodic-Parabolic Boundary Value Problems and Positivity, Volume 247, Pitman Research Notes in Mathematics, Longman, (1991). 17. J. Lopez-Gomez, Positive periodic solutions of Lotka-Volterra reaction-diffusion systems, Diff. and Znt. Eqns. 5, 55-72, (1992). 18. X.-Q. Zhao and V. Hutson, Permanence in Kolmogorov periodic predator-prey models with diffusion, Nonlinear Analysis 23, 651-668, (1994). 19. N.F. Britton, Spatial structures and periodic travelling waves in an integro-differential reaction-diffusion population model, SIAM. J. Appl. Math. 50, 1663-1688, (1990). 20. S.A. Gourley and N.F. Britton, A predator prey reaction diffusion system with nonlocal effects, J. Math. Biol. 34, 297-333, (1996). 21. J.M. Cushing, Integrodifferential equations and delay models in population dynamics, In Lecture Notes in Biomathematics, Volume 20, Springer-Verlag, Berlin, (1977). 22. M.G. Crandall and P.H. Habinowitz, Mathematical theory of bifurcation, In Bifurcntion Phenomena in Mathematical Physics and Related Topics, (Edited by C. Bardos and D. Be&s), pp. 3-46, Ffeidel, Dordrecht, (1980). 23. W. Kaplan, Operational Methods for Linear Systems, Addison-Wesley, Heading, MA, (1962). 24. D. Schley and S.A. Gourley, On a scalar differential equation with a temporally oscillating piecewise constant delay, (In preparation). 25. M.V. Bartuccelli, CR. Doering, J.D. Gibbon and S.J.A. Malham, Length scales in solutions of the NavierStokes equations, Nonlinear-ity 6, 549-568, (1993). 26. M.V. Bartuccelli, J.D. Gibbon and M. Oliver, Length scales in solutions of the complex Ginzburg-Landau equation, Physica D 89, 267-286, (1996). 27. S.A. Gourley and M.V. Bartuccelli, Length scales in solutions of a scalar reaction-diffusion equation with delay, Phys. Lett. A. 202, 79-87, (1995). 28. M.V. Bartuccelli, S.A. Gourley and C.J. Woolcock, Ladder theorems and length scales in solutions of a generalized diffusion model, Physica Scripta (to appear).
PII: SO8957177(97)00237-9
A Population Model with Time-Dependent Delay M. V. BARTUCCELLIAND S. A. GOURLEY* Department of Mathematical and Computing Sciences University of Surrey, Guildford, Surrey GU2 5XH, United Kingdom Qn.bartuccelli>~s.gourley>@mcs.surrey.ac.uk (Received July 1997; accepted Augwt 1997)
Abstract-we present analytical and computational results concerning the linear stability and instability of the uniform steady-state solution of a system of reaction-diffusion equations where a parameter in the kinetic terms is periodic in time. Under suitable sssumptions the system is equivalent to a scalar equation with a periodically varying delay. Such a varying delay can model the seasonal fluctuations to the regeneration time of a resource. We study the effect such a varying delay can have on the stability of the spatially uniform steady-state. Analytical results reveal that instability can set in if the delays are large, while computational methods of analysing the stability equations reveal the precise shape of the instability boundary. The nonlinear stability of the uniform state is also examined using ladder methods.
Keywords-Reaction-diffusion
equations, Time delay, Stability, Parameter domains.
1.
INTRODUCTION
Oneof the most commonly studied problems in developmental biology is that of pattern formation arising from an initially uniform state. This problem has been widely known in the literature as the Turing bifurcation since the seminal work of Turing in 1952 [l]. He studied reactiondiffusion systems and discovered the mechanism by which a uniform state can become unstable under the effect of diffusion. Since Turing’s discovery, many authors have applied and generalized his idea to a number of systems arising in ecology and embryology, and a good review is given in Murray [2]. I n recent years, many authors have extended Turing’s ideas by considering the effect of allowing the kinetic and the diffusion parameters to vary in space, in order to take into account spatial inhomogeneities. Ecological reaction-diffusion models which incorporate spatial variation in the kinetic terms were considered by Shigesada [3] and Cantrell and Cosner [4]. The effect of spatial variation in the diffusivities was considered by Maini et al. [5] and by Benson et al. [6,7]. Another extension, having more relevance to our work in the present paper is the possibility of parameters varying with time, and we shall briefly summarise those works that particularly motivated us here. Timm and Okubo [8] analysed the planktonic predator-prey model of Levin and Segel [9] with a temporal sinusoidal variation introduced into the diffusivity
of one of the species. Their particular type of temporal variation indicated a stabilising elfect and Gourley et al. [lo] proved a more general result to this effect for two-species reaction-diffusion models with small amplitude variations in the diffisivities of both species. Sherratt [ll] analysed *Author to whom all correspondence should be addressed.
13
14
M.V.
BARUCCELLI
AND S. A. GOURLEY
the case when the diffusion coefficients take the form of a periodic step function which alternates between two constant values but need not have small amplitude. In fact, he showed that large amplitude oscillations of this form can have a destabilising effect. This work was developed further in 1121. Though not having such direct relevance to the present paper, there are numerous other papers concentrating on the effects of periodic forcing in ecological modelling, for example, [13-151. Existence theory is treated in (16,171, while criteria for permanence (long term coexistence) are presented in 1181. Introducing time dependence into the parameters of population models allows us to incorporate daily, seasonal, or annual fluctuations in the values of those parameters and the importance of these effects is clear. However, in population models where delay is incorporated, for whatever reason, there naturally also arises the question of whether the delay itself should vary with time. For example, delay is sometimes incorporated to take into account the regeneration time of a resource, but there is no reason why such a regeneration time should always be the same-it is quite likely to be subject to seasonal fluctuations. In the present paper, we shall study the system it = Au + f&v), (1.1)
vt = ~Av + @(t)(u - v),
in which u is the population density of some species, v is a weighted average of past values of u (see below), 13(t) is a temporally periodic parameter related to the delay in the system (in a sense to be discussed shortly), and the constant ,u satisfies 0 < ~15 1. If 8 is constant or piecewise constant then the second equation can be solved to give v in terms of u: v=
t +CQexp ss
-((’ - ~““‘*~” - ‘))gexp
--oo -co
[4np(t - s)]1’2
[-@(t
-
s)]u(y,
s) dy &
(1.2)
justifying our assertion that v is delayed 2~. For more general periodic functions f?(t), (1.2) is only an approximation but it becomes increasingly accurate when the period of 8 is large. The type of application we have in mind involves the period being large so this is not too much of a drawback. The model we are studying is thus a single species reaction-diffusion population model in which f(u, v) depends on the current density ~(5, t) and also on the spatio-temporal convolution v, given by (1.2), which takes account in a particular way of past densities and nonlocal spatial effects. It follows closely the type of model introduced by Britton 1191and Gourley and Britton [20] who argued that whenever delay effects are incorporated (for whatever reason) into a diffusion model, it is necessary to take into account the fact that because of the diffusion the individuals will not have been at the same point in space at times in the past. To account for this a detailed derivation is provided in [19], using random walk arguments, for the use of the double convolution (1.2) (with @ = 1) as an expression for the lagged density, and in that paper it is also observed that such an expression satisfies a differenti~ equation, namely, the second equation of (1.1) (with p = 1). In the present paper, we will in fact allow ,Qto be any constant satisfying 0 5 p < 1. We then have a choice of either taking ,u = 1 or we can take the limiting case /I --) 0, in which case expression (1.2) can easily be seen to degenerate to a purely temporal convolution of the kind frequently seen in ordinary differential equation models with delay 1211. Turning now to the significance of B, note from (1.2) that if 0 is large the delay in the system is weak, in the sense that although all values of u in the past contribute to the spat&temporal average v, only recent values can have a strong influence. Decreasing the value of 6 increases the influence given to values of u far in the past. It is quite reasonable, therefore, to think of the quantity l/e as being a quantitative measure of the delay in the system. In the present paper, we are allowing 8 to be a periodic function of t, alternating between two constant values. Thus,
Population
Model
15
the delay will do likewise and we can regard this as a primitive means of modelling situations where the delay in the system is a periodic function of time. For example, delay is often used to model the regeneration time of resources and this is quite likely to be longer in winter than in summer. Numerous other possible reasons for incorporating delay in population models are described in [21]. As far as time-varying delays are concerned, the use of a piecewise constant periodic delay is perhaps the simplest case that allows progress to be made analytically. As the notation implies, f(~, v) allows dependence on both current density u and lagged density v given by (1.2). Any function f that is biologically realistic and sufficiently smooth is allowable. An obvious possibility is the logistic case f(u, v) = ~(1 - v). An interesting generalisation of this is given by f(~, v) = ~(1 + a~ - (1 + a)~), where a is an aggregation parameter. This particular case will be considered in Section 2.2. In Section 2, we present the linearised stability analysis for the spatially uniform steady-state solution of our model. Section 3 then addresses the question of nonlinear stability using ladder theory. Then, in Section 4, we study the linearised stability criteria numerically to augment our analytical results and finally, in Section 5, we outline what conclusions may be drawn.
2. LINEAR
STABILITY
ANALYSIS
We assume that the system (1.1) has a spatially uniform steady-state solution (u, v) E (u*, v*) such that (u*, v*) > 0. We are interested in analyzing the linear stability of this uniform state. Notice that the second equation of the system implies u* = v’, so the steady state satisfies f(~*, u*) = 0. Setting 21 = u* + ii, v = 2L*+ 6, substituting
into (1.1) and then dropping the tildes, we obtain the linearised system ut = AU + au + bv,
(2.1)
vt = j.~Av+ B(t)(u - v).
In (2.1), a and b are the partial derivatives of f(u,v) with respect to the variables u and v, respectively, evaluated at the uniform state. The limiting case of zero delay is found by letting 0 --) 00 and from (1.2) it is easily seen that in this limit v + u (see, also, [20]). Then the first equation of (2.1) becomes ut = Au f (a + b)u.
(2.2)
Hence, if a + b < 0, the uniform state is stable, with and without diffusion. We will henceforth assume that a + b < 0, since we are interested in instability for the full system that is delay induced, rather than instability that is present even without delay. We shall now analyse the stability of the uniform state for the case when 8 is constant. Then we shall investigate a particular time-dependent 8 in the form of a periodic step function. Seeking solutions of (2.1) (with ~9constant) of the form (u, v) = (cl, ~2)exp (st + ik . x) , we find that
s-a+K
-8
,+0&K)
(::)
= (:)’
where K = k2 = k . k. For nontrivial solutions, we need the determinant zero, so s2+(tl+pK-a+K)s+(B+pK)(-a+K)-M=O.
of the matrix to be (2.3)
16
M. V. BARUCCELLI AND S.A. GOURLEY
Instability of the uniform state can set in either through a steady-state bifurcation (s = 0), or a Hopf bifurcation. For a given wavenumber k, a steady-state bifurcation may occur when (4 + /.&K)(-a -I- K) - be = 0,
(2.4
0 = PK(-a + K) (a -I-b) - K ’
(2.5)
i.e., when
Since a + b < 0 and B > 0, we see imm~iately that steady-state bi~rcatio~ will not occur if a < 0. To seek a Hopf bifurcation we set s = iw, with w real and positive. Such a ~furcation may occur when 0-a+(p+l)K=O, (2.6) provided that (Q $ pK)(-a
+ K) - b@ > 0.
(2.7)
Hence, again one can see that if a < 0 there are no Hopf bifurcations. Thus, if a c 0, the uniform state is stable regardless of the value of 0. Now let a + b < 0, but a > 0. It is convenient here to think of l/e as a bifurcation parameter, since it represents delay. Small delay usually implies stability and, in this particular case, graphical arguments help to show that the uniform state is stable to all su~ciently small ~~urbations when
(2.8) As I/0 is increased instability may set in either through a Hopf bifurcation (when 8 = a), or a steady-state bifurcation, depending on the values of the other parameters (for example, if ,G is sufficiently small we expect a Hopf bi~rcation). Further technical conditions need to be checked to be certain that bifurcation actually occurs and a useful reference here is the article by Crandall and Rabinowitz [22). The bifurcating solutions could be constructed using methods similar to those in [20], in which the noncompactness of the spatial domain is dealt with by restricting attention to spatially periodic solutions. In the case of a steady-state bifurcation, the critical wa~n~ber at the bi~~ation can be cumulated and is the value of K E (0, a) at which the relevant minimum in (2.8) is attained. This completes the analysis for the case when 6’ is constant and, as we have seen, it is quite straightforward to enumerate the various possibilities. We are now going to introduce temporal variation in 19and study the case when e(t) is a periodic step function with period T, so that
e(t) =
(2.9)
where n is any integer. This choice of e(t) follows the work of Sherratt Ill] and will enable us to make some analytical progress towards determining conditions for stability/instability of the uniform state. Introducing e(t) defined above into system (2.1), we obtain ut = Au + au + bu, vt = PAW + e(t)(~ - vf.
(2.10)
To study stability of the uniform state to perturbations of wave vector k, we set (u,v) (U(t), V(t)) exp(ik . x) and substitute into (2.10) to obtain u, = -KU v, = -pKV
=
-t aU + bV,
+ e(t)(u - V),
(2.11)
Population Model
17
where K = k2 = k a k. Denoting w = (V, V)T, we can write this in the form
dw=Afw
’
dt
d”=A-w dt
’
where
The method for solving systems like (2.11) is Floquet theory, and we summarise this theory very briefly. For more details see, for example, [23]. Floquet theory applies to linear ordinary differential equations with periodic coefficients. The solutions of this class of equation will generally be a periodic function times an exponential. The general technique is to solve the system of linear differential equations for any set of complete linearly independent solutions, for 0 5 t 5 T. This yields a so-called propagator matrix E such that the general solution can be expressed as ~(t + mT) = E”“y(t), where y(t) is the general solution of the system of order n, andm=1,2,3 ,.... Note that the propagator E can be expressed in terms of any fundamental matrix Q(t) as E = @-l(O)@(T). The eigenvalues of E are known as the Floquet multipliers of the system. We shall denote them by pi (i = 1,2,. . . , n) for a system of order n. It can be shown that the characteristic multipliers do not depend on the particular choice of the fundamental matrix a(T), and therefore, are intrinsic properties of the system. Furthermore, to determine the stability of the solutions of the problem, one needs to compute the eigenvalues of the propagator matrix E. From y(t + mT) = E”‘y(t), it follows that the uniform state is stable if and only if all the eigenvalues of E have modulus less than one. Hence, if any of the eigenvalues is larger than one in magnitude then the solution is unstable. We shall use this theory to determine the Floquet multipliers, which in general are complex numbers, and to give information about the stability of the solutions of (2.11). Following the notation used in [ll], we denote the eigenvalues of A* by Xt2 and the eigenvectors by z$. We further define Zf to be the matrix whose columns are the eigenvectors of A*. After some algebra, we find that
zf,
Xf=i(a-B*-(p+l)K+P’), (2.12) X:=l(a-@-(p+l)K-P’),
and that
-b
4
;(Qi-P*)
’
;(Q&+Pk)
where l’* = 44b6’* + (Q*)2,
(2.13)
Q*=a+8*+(~-l)K. Defining n*(t)
=
0
exp(Xft)
(
o
exp(Xtt)
)
’
it follows that any fundamental matrices of (2.11) have the form Z+A+(t)C+ and Z-A-(t)C-, We take C+ = where C+ and C- are matrices whose entries are constants of integration. (A+(T/2))-l. Because we want continuity at t = T/2, we require that C- = (A-(T/2))-1 (Z-)-‘Z+. Hence, we obtain E = (@(O))-%(T)
= (2+)-l
T
z
Z-A0
(Z-)-‘Z+A+
0
;
M. V. BARUCCELLIAND S. A. GOURLEY
18
The last matrix A+(T/2) in the above actually appears at the front in the computation of E, but we have written it at the back for convenience in the subsequent calculation. Matrix multiplication is not, of course, commutative, but what is important are the e~ge~vu4~esof E and we are using the result that if A and B are any square matrices, then AB and BA have the same eigenvalues. The uniform state will be stable (to perturbations of wavenumber k) if and only if the two eigenvalues of E have modulus less than 1. After some algebra, we can see that the matrix E has the form E = exp (-lY’/4)fi‘,
where (2.14)
r=-2(X,+X;t)+~++P-=8++8---2[a-(~+l)K] and B=exp
(P+ + P-)T [
4
x
1
[(z+W (: exp(+))Pm’z+(; _(~))].
It is straightforward to see that det(@) = 1. Hence, the characteristic equation of B has the form (2.15)
&++1=0, where B, which is the trace of k is given by i3={(l-exp(e))
2b (e+ + e-) + Q+Q-
(
(I-exp(F))
+ (l+exp(T))
P+P-
>
(2.16)
(l+exp(~))}expi(P+~P-)T’41,
with Q*, P* given by (2.13). Analytical results are difficult to obtain in general, but it is possible to make progress by considering various limiting cases. We start by considering the case when Bf and 8- are very close together. Then we see what can be said if they are widely different. We then treat the case when 9+ and 13~are both very large, and after that the case when both are very small. These two cases correspond to small and large delays, respectively. Then we shall prove that if a + b < 0, a < 0, and b > 0, the uniform state is stable for all values of Bf and 8-. Note that the condition for stabidity of the uniform state, as a solution of the reaction-diffusion system is that for every K>O,weneed lexp(F)
Q/Cl,
(2.17)
where fi satisfies (2.15). Note also that the coefficient B in (2.15) need not be a real number, since one or both of P+ or P- can be purely imaginary. 2.1. Small Amplitude
Periodic
Variation
To keep the algebra as simple as possible we will, in this section and the next, consider stability of the uniform state to spatially uniform (K = 0) perturbation only. This is equivalent to studying its stability as a solution of the corresponding ordinary differential equation system, i.e., system (1.1) without diffusion terms. We assume, as always, that a + b < 0. Recall that in the case of constant 8, instability of the uniform state is possible if a > 0 and, in this case, that (2.8) gives the value of 8 at which instability actually sets in. We will assume in this section that a > 0 (note this implies b < 0). Since K = 0, the instability sets in through a Hopf bifurcation when 9 = a, with stability when 0 > a.
Population Model
19
Turning now to the case of varying 8, our aim is to get some insight into how this periodic variation changes things. Analytical progress can be made by considering the effect of a small amplitude periodic pe~urbation whose mean value is a, so we set 8+ = a + E and B- = a - E, whereO
with a similar expression for P”, and that
’ sin’ fT,“i;;T;L”-“;jr o2(a + b)
TbmsinTdM~&&~(u+ b)
.
(2.18)
Since I’ = 0, condition (2.17) becomes Ifi] < 1, where fi2 - Bb f 1 = 0. For “most” values of the period T of B(t), the leading order term in the above expression is less than two in absolute value and this implies, given s is small, that also IEI 5 2. The roots 3 are then a complex conjugate pair on the unit circle. So this implies the system usually remains marginally stable under the effect of a small amplitude periodic perturbation. The values Tn:=: d&,
n=l,2,3
(2.19)
,..,,
need to be considered in more detail because the leading order term in (2.18) then has modulus 2. Note that these values are closely related to the period of the Hopf bifurcating solution in the constant B case, which can be shown to be
The cases n even and n odd have to be looked at separately as follows. CASE 1: n ODD. In this case, cosT,dm Tn dj
= 0, and sin2(1/2)
= -1, sinT,dm
= 1. Thus,
BZ-2-E2
b <-2, u2(a + b)
since a + b < 0 and b < 0.
The roots fi are therefore real, and are located just on either side of -1. Thus, the uniform state becomes unstable. CASE 2: n EVEN. In this case, the 0(s2) term in (2.18) is zero, so higher order computations are needed. We have done some numerical simulations, however, and these suggest the uniform state is unstable, but much more weakly so than in the case of n odd. 2.2. Large
Amplitude
Periodic
Variation
In this section, we again assume K = 0 and study stability to spatially uniform perturbations only. What we will do is show that under certain conditions a large ~plitude periodic variation in 6(t) may have a stabilising effect. We shall illustrate this phenomenon with reference to the population equation
t s -00
6e-e@-S)u(s)
ds
, >
20
M. V. BARUCCELLIAND S. A. GOURLEY
with a, 8 > 0 and its equivalent system dtt z
=%6(1+at4 -(Ifo)u), (2.21)
dv
-& = L9(u- v). This system (in the constant 6 case) is one of numerous special cases of a more general equation studied in [19j. The parameter a measures the tendency of the species to aggregate or group together as a defense against predators. The quantity a used elsewhere in this paper as a coefficient of the linesrised system is equal to the a of this particular model. The linear&&ion of (2.21) about (u, v) = (1,l) is in fact dU z
= aU
-
(1 i- u)u, (2.22)
dv -& = @(u-v),
so, in this model, b = -(l -t a). In the constant 0 case, a Hopf bifurcation occurs at 8 = a. We now introduce variation in 8 according to (2.9) as usual, but not necessarily with mean value a this time. If 8’ is large and 8- is small, say, this guarantees Pi and P- are both real (and also that I? > 0). Now let the period 2’ -+ 00, then B
2
O(1) exp
(P+ + P-) T 4 ( > s
The two roots of pz - Bfi + 1 = 0 are, therefore, (~ymptotic~ly) B and l/B. When ji = l/B, it is quite obvious that 1exp (-I’T/4)$I < 1. Fbr the root fi = B, the stability condition becomes ]O(l)exp(~)exp[(P++~-)T]/cl,
forTlarge,
or I‘ > P+ + P-. To simpli~ this, we shall exploit the assumed largeness of 8+ by using an asymptotic expression available from the next section (expression (2.23)). Also, if we let 8- 40, then P- + a. The stability condition reduces down to a < 1, which is perfectly possible biologically and means that the tendency of the species to aggregate is not too large. Since large 0 means small delay and vice versa, what we have done in this section is to demonstrate that if one delay is very small, the other is large and also the period is very large then stabilisation is a possibility, depending on the values of the other parameters. Actually this result did not surprise us. The differential equation y
=
u(t) (1 - u (t - o(t))),
with a(t) =
has recently been studied by Schley and Gourley [24]. In this equation, the delay one, alternating between zero and 2~s. We have shown that if crs is related to the T = 4~70,the constant state u E 1 is stable for any 00 > 0, and in some respects this the one of the present section. For smaller periods, however, the stability depends of 00.
is a discrete period T by result is like on the value
Population Model
21
2.3. The Case 19+and 8” Large We now return to the general linearised system of the paper, system (2.11) with (2.9), and we allow perturbations of arbitrary wavenumber so that K 2 0. When 8+ and 6- are large we can show that P+-8++(2b+a+&-l)K)+O
& (
with a similar result for P-.
(2.23)
, >
Other relevant facts, when 0+ and 8- are large and that r > 0,
P+, and P- (and hence, B) are real, and B N exp[(P+ + P-)T/4] >> 1. Again the two roots of fi2 - Bji + 1 = 0 are asymptotically B and l/B. Much as in the previous section, the stability condition reduces to
r>p++p-. From (2.14) and (2.23), it follows that we require 4(a + b) - 4K < 0 for all K 2 0, and this is true because a + b < 0 by hypothesis. Therefore, for 8+ and 0- large (i.e., when the delays are small), we have stability of the uniform state. Small delay usually does imply stability, and in this section, we have shown this to be the case when a periodic delay oscillates between two small values.
2.4. The Case 8+ and 8- Small When 8+ and 8- are small, we can see that 2b
P*-la+(p-l)KI+
+
a + (cl - 1)K
Ia + (p
_
l)KI
‘* +’ (e*2) ’
>
(2.24)
and that Bpexp
1 [
(P++P-V 4
-(P++P-)T
+exp
4
I
= 2 cash The roots of fi2 - Bji + 1 = 0 are, therefore,
1*
f(P++P-)T
ji = exp [
4
The condition for stability is exp[-I’T/4] p < 1. This requires -I? f (P+ + P-) and P- are positive, l? > P+ + P-. Using (2.14) and (2.24), this gives 0++6--2[a++l)K]
>2la+(p-l)KI+
2bl+;;($!,K)
< 0 or, since P+
(0++8-).
(2.25)
It is convenient to consider the two cases a + (II - l)K < 0 and a + (p - l)K > 0 separately. CASE 1: a -t (/.J- l)K < 0. In this case, (2.25) reads, after some algebra
(e+ + e-)
a+b+b---1)K a + (p - 1)K
>
_zpK
1
Since a + (cl - l)K < 0, a + b < 0, and 0 5 p 5 1, the above inequality is trivially satisfied, since the square bracketed term is positive. CASE 2: a + (p - l)K < 0. Note that for this to hold it is necessary that a > 0, and this in turn implies b < 0 (since a + b < 0). If these conditions are satisfied then a + (/J - l)K is positive for all wavenumbers K if /J = 1, and if p < 1 it is positive for wavenumbers satisfying
M. V. BAAUCCELIJ ANDS. A. GOURLEY
22
K E [O,(u/l - CL)).It turns out in this case that instability is possible under a further condition which we shall obtain. Reversing the inequality in (2.25), we obtain the condition for instability that e++e--2(a-(p+l)K) or after
2b + a + (fi - 1)K a + (/J - l)K
<2(a+(/J-l)K)$
(e+ + 8-) ,
(2.26)
some algebra that Q(K) := 2(K - a) (a + (p - 1)K) < b (8+ +
@-),
(2.27)
for some K E [O,(a/l - p)), or [O,00) if p = 1. Note that the right-hand side of (2.27) is negative, since b < 0. In both of the two cases /.J= 1 and p< 1, simple graphical arguments will show that we require Q(0) < b(B+ + P), and this holds whenever
e++ewhich we can regard as automatically
2
< $fq,
satisfied, by the smallness of 8*.
Thus, if 8+ and B- are small (i.e., if the delays are large), then our conclusions are as follows. If a < 0, the uniform state is stable (because Case 1 applies for every wavenumber). If a > 0, then Case 2 applies, either for every wavenumber or for a finite range of wavenumbe~, and in either case there are wavenumbers which give instability of the uniform state. 2.5. TheCasea+b
G fi=exp@‘T/4)
(2.29)
> 0,
both hold. To proceed we now need the following lemma. LEMMA
1. Ifa+b
then 2b (e+ I
+ e-1 + Q+Qpip-
PROOF. After some algebraic manipulation, we see
b(a+b+(~-l)K)IO, which is clearly true under the stated assumptions.
< I
-
I ’
that the above inequality is equivalent to
Population Model
23
In order to prove (2.28), we proceed as follows. Now ~(exp(~))
=exp(a>
-~exp(~)
+l,
(2.30)
and B={(l--exp(T))
(I-exp(T))
+ (1 +exp (q))
(1 +exp
(2b(8f+~~J+Q+Q-) (l$2))}
exp[(P+
~P-)*~41,
Using Lemma 1, B 5 f exp (P++:S)T]
{(l-exp(T))
[
+ (l.+exp(T))
(l-exp(T))
(l+exp(F))}.
(232) .
We point out here that Lemma 1 can also be used to confirm that B > 0.Now let us introduce ,
We find, using (2.32), that f (exp (3))
2 z2 - i XYZ[(1+;)(1++(1-$)(I-$)I+1 =,r2--2 =
(
(z -xy)
Cry+;
)
+1
(*-$ ><
The second bracket in the last line is trivially positive, since z,y > 1 and z > 1. The first will be positive if t > xy, that is, if exp @‘T/4) > exp (P+ + P-)T/4 or I’ > P+ + P-. This is equivalent to e++e--z[a-(/.A+1)K]
> da+
and a sufficient condition to ensure this is that the inequalities 6+ - (a - (p + l)K) > 44bB+ + (Q+)2, e- -(a--(p+l)K)
> da,
should both hold, which they do if and only if 4(a + b)8f + 4K (a/A- e* - /AK) < 0. The latter is certainly true, because a < 0 and a + b < 0. This completes the proof that
f (em (rT/4)) > 0.
24
M. V. BARUCCELLIAND S. A. GOURLEY
It remains to show that %I P=exp(r~,l~ > 0. Using the same variables z, y, z as above, we get
>2z-,,(I+--&),
by(2.32) but .z > zy
>XY--=
1
X2Y2 - 1
XY
XY
2 0,
as desired. For cases when B can be complex there is not much we can do analytically and we have to resort to analysing the stability criteria by numerical inv~tigation of a parameter space of interest. The results of some such numerical investigations are presented in Section 4. When the result of the present section is applied to the particular model studied in Section 2.2 (equation (2.20)), the conclusion is that if a < -1, then the steady-state u z 1 of that model is stable for all values of 8” and I!?“. Britton [19] proved stability in the case of constant 8 under the condition a < -l/2. In the next section, we shall consider the question of the nonlinear stability of the uniform steady-state solution of the system (1.1).
3. NONLINEAR
STABILITY
ANALYSIS
In this section, we exploit some of the ladder methods developed by Bartuccelli et al. [25281 to prove some nonlinear stability results for the original system, which we restate here for convenience: it = Au + f(u, w), (3.1) vt = ~Av + O(t)(u - v). In order to make the ladder methods available we shall assume, in this section, that we are working on the finite spatial domain x f [O,Lid, with periodic boundary conditions, where d is the spatial dimension. Also, for the time being, O(t) is any positive differentiable function and not necessarily periodic. The analysis outlined below will hold for any biologically meaningful nonlinear term, but to keep the algebra simple we restrict ourselves to the particular case vzt = Aa+u(l
-cru-flu),
2rt = ~Av + O(t)(u - v),
(3.2)
where CYand p are positive constants. We are interested in analyzing the nonlinear stability of the uniform state. To use the ladder methods effectively, it is very desirable to centre the system on the uniform steady&ate (2t*, u*), so we again introduce u =
us+ ii,
v =
u* + 6,
where in this system u* = l/(cy + p), substituting into (3.2) and dropping the tildes (but not the higher-order terms this time), we obtain the nonlinear system ut=Au-azb-bv-~~2-/3~v, vt = /.JAv + O(t)(u - v),
(3.3)
where a = &/(a! -I-@) and b = /3/(a + /3). Note that a + b = 1, and so we are in the same hypotheses of the linear analysis of the previous section, since a and b in (3.3) are equivalent to
Population Model
25
-a and -b in the linear system (2.1). The u and v in (3.3) are not positive functions, of course, since they are really fi and 6. However, solutions of the original system (3.2) are undoubtedly positive if the initial data is, by the maximum principle. In view of this, we can say that solutions of (3.3) satisfy u* + % > 0 and u* + v 2 0. This will be very important later. We introduce the positive definite quantities Hn = 8
J
(Vn~)2 dx,
Jn = b
J
(V’%Q2 dx,
Fn=H,+J,.
(3.4)
If we differentiate H, with respect to time and use (3.3), and then integrate by parts on the Laplacian term (the integrated terms vanish by virtue of the periodic boundary conditions), we obtain iI!& = -H,,+l
- aH, - be / (V%) (V’%) dx J
Y
-
a~J
(vnu2) dx - pe
(vu)
J (vu)(VW)dz + f
(3.5) j+.
Doing the same calculation for J, gives ;j_
=
J
-pJ,+l + b8
(V%) (V”u) dx - SJ,.
(3.6)
Adding the & and the JR, we get a differential equality (ladder) for the F,. Now it turns out that for investigating stability of the uniform state it is not necessary to analyse the ladder for arbitrary nonnegative integers 7~2 0, but it is sufficient to study the bottom rung of the ladder, which is given by n = 0. We shall then obtain what it is known ss nonlinear stability in the L2 norm. Adding the H, and the J,, differential equalities for n = 0 gives (3.7) The above first-order nonlinear ordinary differenti~ equation contains all the information we need to analyse the stability of the zero solution of (3.3), and hence, of the uniform state (u*, u*) of the original system (3.2). Naturally the stability will depend, as in the linear stability analysis already done, on the functional expression of the kinetic parameter 19(t). First, we note that the fifth and sixth terms in the right-hand side of (3.7) are in general of indefinite sign. But, as previously noted, we know that U* + u I 0 and U* + v > 0. This gives --u < U* and -v < u*, and we use this information to obtain the estimates
J u3
-CUB
dx 5 CXU* Ho
and
-PO
J
u2vdx I /3u*Ho.
Employing these estimates in (3.7) yields a differential inequality
;po <-qJ1 Note that the HO term accounts a “linearised” inequality
-HI
18 - aHo - OJ, + Ho + z sHr,.
for the nonlinear terms in (3.3). By dropping it we can obtain
;po-pJ1 5
- HI - aHo - 6Jo + f ;Ho,
(W
from which results can be found for the linearised system. Let us now list various consequences of both the full and the linear&d inequalities (3.8) and (3.9). The idea is to try to reconstruct Fo in the right-hand side of (3.8), and then find
M. V. BARUCCELLIAND S. A. GOURLEY
26
conditions on 19(t) and the other parameters which will guarantee that FO + 0 as t ---f 00. This will give nonlinear stability, in the L2 norm, of the zero solution of (3.3). First, note that, by Poincare’s inequality, Hi 2 I~-~I-&,, and Jr 2 K2Jo,. Hence, (3.8) becomes
“@ <2
Ho -
O-
($ +3)Jo 5 -&3)Fo,
where 4(t) = min
$+a-1
-$$),
(S+@(t))].
K Hence, if (3.10)
linl~f 4(t) > 0, then FO --f 0 as t -+ 00. There are several consequences of this in the following.
(1) Provided
B and 8 are bounded for all t > 0, we can ensure (3.10) holds by taking L sufficiently small. Thus, if L is sufficiently small, the uniform state is stable. This is in accordance with the general observation in reaction diffusion systems that spatial patterns cannot emerge if the spatial domain is too small. (2) 8(t), a nonincreasing function of time (in particular, constant). Under this assumption the term (I/2)(~/~)~ e is nonp~itive, and therefore, from (3.10), it follows that Fe(t) + 0 as t ---f cm provided that min K
;+a-1=;-b)‘(s+e(t))]
>o.
(3.11)
(3) B(t)
large. In this csse, ~1/2)(~/~)~0 is small and so & < 0, provided again (3.11) holds. Therefore, we have stability also for large values of 8, i.e., for small values of the delay.
Notice that the argument in Cases (l)-(3) a b ove, also holds for nonperiodic functions 8 = O(t), imposed.
so long they satisfy the restrictions
(4) Assume now that e(t) is a periodic function of time, and sssume its period is large; let us introduce the new time variable i = t/T < 1; this will enable us to analyse, in particular, the stability of our uniform state for large values of the period of O(t). The Fo ladder now becomes 14 f&j L 2’ (-/A& - Hi - aHi, - SJ, + Ho) + 2 5Ho. (3.12) Therefore, if 2’ is sufllciently large the (1/2)(~/~)Ho function 4(t) defined above becomes
term can be neglected,
and the
$*(t)=min[(++a_1).($+B(t))]. Therefore, provided (3.11) (5) If the period of our kinetic This can be seen from (4) and so Fo increase in time,
is satisfied, we again obtain that the uniform state is stable. parameter B(t) is small we can show that we have instability. above: when T is small the 6/0He term is the dominant one that is, we have instability.
This concludes the nonlinear analytical analysis of the system.
Population
(a) p=
Model
1, a = 4, b = -5,
27
T = 4.
4.0
3.0
z-
2.0
1.0
1.0
2.0
(b) /A = 1, a = 1, b = -2,
3.0
4.0
T = 4.
Figure 1. Diagram showing typical regions of stability and instability in the paramfor the zero solution of the eter plane (r+, T-), where r+ = l/e+, T- = l/e-, linearised system (2.10) with 0(t) defined by (2.9). In the light regions, the zero solution is stable to perturbations of arbitrary wavenumber while in the dark regions there is at least one unstable mode.
4. NUMERICAL
COMPUTATIONS
It has become apparent by now, that in order to make analytical progress with the type of system we are studying in this paper, it is almost invariably necessary to consider various limiting
M. V. BARUCCELLI AND S.
28
A. GOURLEY
5.0
4.0
3.0
z2.0
1.0
0.0
I
=
1, a = 0.5,b = -0.52,
P =
a = 0.5, b = -0.52,
(c) P
T = 4.
5.0
4.0
3.0 2-
2.0
1.0
0.0 2+
(4
0,
‘1’ = 4.
Figure 1. (cont.)
csses and use ~ymptotics. The only exception to this was in Section 2.5, but the gumptions made there are difficult to interpret and the methods used in that section are unlikeIy to generate the best possible result. This leaves the question of what the stability regions really look like in parameter space. We wrote a simple computer program to determine regions of stability and instability in the parameter space (T+,T-), where T+ = l/6+ and T- = l/19-, so the r’s may be interpreted as
Population Model
29
delays. At each point on a grid of parameter values, the program simply checks the stability condition 1exp (-P/4) fi] < 1 for each of a large number of (equally spaced) wavenumber values up to a suitably large wavenumber. If instability is found for any of these wavenumbers, the point is coloured black. If all wavenumbers have been checked and no instability found, the point is left white and the computer moves to the next point on the grid. Some stability diagrams generated by this simple code are shown in Figure 1. The diagrams are all symmetric about the line T+ = T- as we should expect. In all four diagrams, one can see that we have stability for small values of the parameters r+ and T-, as expected (small delays give stability). Two of the diagrams also bring out rather clearly the phenomenon noted in Section 2.2, for one specific system, that if one delay is very large and the other very small, then the result can be stability. Our diagrams suggest that this phenomenon may not just be a peculiarity of that particular system and also that it may occur even when the period T is not particularly large. Recall that in Section 2.2, we could only demonstrate it analytically for very large values of T. Various other experiments suggest a general result that if a > 0, there is always an instability region in the parameter plane we have considered. We were not able to get instability for any negative value of a, we looked at but have no analytical result to this effect. Recall that to establish the stability result of Section 2.5 it was necessary to assume not only that a + b < 0 and a < 0, but also that b > 0. The latter assumption was necessary for our particular method of proof to work but we believe the result may hold true even without this assumption. If this is indeed the case it means that in the a + b < 0 and a < 0 situation the periodic variation does not change stability, since these conditions certainly do guarantee stability in the constant 19case. 5. CONCLUSION In this paper, we have described a simple means of incorporating a time-dependent delay into a reaction-diffusion population model, while still keeping the model amenable to analytic treatment. The periodic delay can be regarded as a means of modelling periodic fluctuations to the regeneration time of a resource and the periodic step function is the simplest case which can be dealt with analytically. We have analysed how the stability of the spatially uniform steadystate is affected by the temporal variation of the parameter 13,which is the reciprocal of the delay in the system. We have found that if a > 0, then a small amplitude periodic variation in e(t) will be destabilising if a resonance effect occurs, that is if the period of e(t) matches that of the Hopf bifurcating periodic solution that we would expect to see if t9 were constant. We have also seen that if conditions are right, it is possible for a large amplitude variation to be stabilising. However, it would be imprudent to draw too general a conclusion on this latter point from what we have done here and further investigation of the large amplitude case is essential. The observations made above would seem to be at variance with those of Timm and Okubo [8], Gourley et al. [lo), and Sherratt [ll], all of whose results suggest that small amplitude perturbations are stabilising while larger ones may be destabilising. In all of these papers, however, it was the diffusion coefficients, and not the kinetic parameters, that were oscillating in time. Furthermore, these authors were considering the effect of temporally periodic perturbations to systems which in the unperturbed case are at bifurcation to a nonhomogeneous steady solution, rather than systems undergoing a Hopf bifurcation to a solution periodic in time. We have also found that if 8+ and C9- are both large then the uniform state is stable, as we expected since this is the small delay case. If they are small the uniform state may be unstable, depending on the values of the other parameters. Furthermore, we obtained sufficient conditions for stability which are independent of 8+ and 8- . Numerically computed stability diagrams confirm the analytical analysis as well as revealing the precise shape of the stability boundary. Finally, we performed a nonlinear stability analysis of the uniform state of a particular system by applying ladder methods. The use of these methods has shown that, as expected, the stability or
M. V. BARIJCCELLIAND S. A. GOURLEY
30
otherwise of the uniform state depends on the size of the spatial domain. The ladder approach also seems to work well for a much broader class of time-dependent functions O(t) than does the linearised
analysis.
REFERENCES 1. A.M. Turing, The chemical basis of morphogenesis, Phil. Runs. R. Sot. Lon. B237, 37-72, (1952). 2. J.D. Murray, Mathematical Biology, Springer, Berlin, (1989). 3. N. Shigesada, Spatial distribution of rapidly dispersing animals in heterogeneous environments, In Lecture Notes in Biomathematics, Volume 54, (Edited by S.A. Levin and T.G. Hallam), pp. 478-491, Springer, Berlin, (1984). 4. R.S. Cantrell and C. Cosner, The effects of spatial heterogeneity in population dynamics, J. Math. Biol. 29, 315-338, (1991). 5. P.K. Maini, D.L. Benson and J.A. Sherratt, Pattern formation in reaction-diffusion models with spatially inhomogeneous diffusion coefficients, IMA J. Math. Appl. Med. Biol. 9, 197-213, (1992). 6. D.L. Benson, J.A. Sherratt and P.K. Maini, Diffusion driven instability in an inhomogeneous domain, Bull. Math. Biol. 55, 365-384, (1993). 7. D.L. Benson, P.K. Maini and J.A. Sherratt, Pattern formation in reaction-diffusion models with spatially inhomogeneous diffusion coefficients, Math. Comp. Modelling 17 (12), 29-34, (1993). 8. U. Timm and A. Okubo, Diffusion-driven instability in a predator-prey system with time varying diffusivities, J. Math. Biol. 31, 307-320, (1992). 9. S.A. Levin and L.A. Segel, Hypothesis for origin of planktonic patchiness, Nature 259, 259, (1976). 10. S.A. Gourley, N.F. Britton, M.A.J. Chaplain and H.M. Byrne, Mechanisms for stabilisation and destabilisation of systems of reaction-diffusion equations, J. Math. Biol. 34, 857-877, (1996). 11. J.A. Sherratt, Turing bifurcations with a temporally varying diffusion coefficient, J. Math. ,Eiol. 33, 295-308, (1995). 12. J.A. Sherratt, Diffusion-driven instability in oscillating environments, Euro. J. Appl. Math. 6, 355-372, (1995). 13. P. Hess and H. Weinberger, Convergence to spatial-temporal clines in the Fisher equation with time periodic fitness, J. Math. Biol. 28, 83-98, (1990). 14. G.J. Butler and H.I. Freedman, Periodic solutions of a predator-prey system with periodic coefficients, Math. Biosci. 55, 27-38, (1981). 15. J.M. Cushing, Oscillatory population growth in periodic environments, Theor. Pop. Biol. 30, 289-308, (1986). 16. P. Hess, Periodic-Parabolic Boundary Value Problems and Positivity, Volume 247, Pitman Research Notes in Mathematics, Longman, (1991). 17. J. Lopez-Gomez, Positive periodic solutions of Lotka-Volterra reaction-diffusion systems, Diff. and Znt. Eqns. 5, 55-72, (1992). 18. X.-Q. Zhao and V. Hutson, Permanence in Kolmogorov periodic predator-prey models with diffusion, Nonlinear Analysis 23, 651-668, (1994). 19. N.F. Britton, Spatial structures and periodic travelling waves in an integro-differential reaction-diffusion population model, SIAM. J. Appl. Math. 50, 1663-1688, (1990). 20. S.A. Gourley and N.F. Britton, A predator prey reaction diffusion system with nonlocal effects, J. Math. Biol. 34, 297-333, (1996). 21. J.M. Cushing, Integrodifferential equations and delay models in population dynamics, In Lecture Notes in Biomathematics, Volume 20, Springer-Verlag, Berlin, (1977). 22. M.G. Crandall and P.H. Habinowitz, Mathematical theory of bifurcation, In Bifurcntion Phenomena in Mathematical Physics and Related Topics, (Edited by C. Bardos and D. Be&s), pp. 3-46, Ffeidel, Dordrecht, (1980). 23. W. Kaplan, Operational Methods for Linear Systems, Addison-Wesley, Heading, MA, (1962). 24. D. Schley and S.A. Gourley, On a scalar differential equation with a temporally oscillating piecewise constant delay, (In preparation). 25. M.V. Bartuccelli, CR. Doering, J.D. Gibbon and S.J.A. Malham, Length scales in solutions of the NavierStokes equations, Nonlinear-ity 6, 549-568, (1993). 26. M.V. Bartuccelli, J.D. Gibbon and M. Oliver, Length scales in solutions of the complex Ginzburg-Landau equation, Physica D 89, 267-286, (1996). 27. S.A. Gourley and M.V. Bartuccelli, Length scales in solutions of a scalar reaction-diffusion equation with delay, Phys. Lett. A. 202, 79-87, (1995). 28. M.V. Bartuccelli, S.A. Gourley and C.J. Woolcock, Ladder theorems and length scales in solutions of a generalized diffusion model, Physica Scripta (to appear).