## D

In dimensional terms D has units cm2s-1 (or km2yr-1 or whatever scale we are interested in) and f '(0) has units s-1 since it is the linear birth rate (for u small f (u) ~ f '(0)u) which together give Lc in centimetres. Thus if the domain size L is less than the critical size Lc, u ^ 0 as t ^^ and no spatial structure evolves. The larger the diffusion coefficient the larger is the critical domain size; this is in keeping with the observation that as D increases so also does the flux out of the region.

The scenario for spatial structure in a growing domain is that as the domain grows and L just passes the bifurcation length Lc, u = 0 becomes unstable and the first mode a1 exp f/(0) — D(r)t sin n x ~L

starts to grow with time. Eventually the nonlinear effects come into play and u(x, t) tends to a steady state spatially inhomogeneous solution U(x), which, from (2.77), is determined by

where the prime denotes differentiation with respect to x. Because f (U) is nonlinear we cannot, in general, get an explicit solution for U.

From the spatial symmetry in (2.77) and (2.81) — setting x ^—x leaves the equations unchanged—we expect the solutions to be symmetric in x about the midpoint x = L/2. Since u = 0 at the boundaries we assume the midpoint is the maximum, um say, where U' = 0; it is helpful now to refer to Figure 2.19(a). If we multiply (2.81) by U' and integrate with respect to x from 0 to L we get

2 J0

since U = um when U' = 0. It is convenient to change the origin to L/2 so that

U'(0) = 0 and U (0) = um; that is, set x ^ x — L/2. Then

which integrates to give

which gives the solution U (x ) implicitly; typical solutions are illustrated schematically in Figure 2.19(b). The boundary conditions u = 0 at X = ±L/2 and the last equation give um

L = (2D)X/2 \ [ F (Um ) - F (w)]-1/2d w ^ Um = Um (L ). (2.84)

Jo um

L = (2D)X/2 \ [ F (Um ) - F (w)]-1/2d w ^ Um = Um (L ). (2.84) Figure 2.19. (a) Typical steady state pattern in the population u governed by (2.77) when the domain length L > Lc, the critical size for instability in the zero steady state. Note the symmetry about L/2. (b) Schematic steady state solution with the origin at the symmetry point where u = um and ux = 0.

Figure 2.19. (a) Typical steady state pattern in the population u governed by (2.77) when the domain length L > Lc, the critical size for instability in the zero steady state. Note the symmetry about L/2. (b) Schematic steady state solution with the origin at the symmetry point where u = um and ux = 0.

We thus obtain, implicitly, um as a function of L. The actual determination of the dependence of um on L has to be carried out numerically. Note the singularity in the integrand when w = um, but because of the square root it is integrable. Typically um increases with L as illustrated in Figure 2.19(b).

Spatial Patterning of the Spruce Budworm

Now consider the model for the spruce budworm, the dynamics for which we derived in Chapter 1, Volume I. Here, using (1.8) for f (u), (2.77) becomes

where the positive parameters r and q relate to the dimensionless quantities associated with the dimensional parameters in the model defined by (1.7) in Chapter 1, Volume I; q is proportional to the carrying capacity and r is directly proportional to the linear birth rate and inversely proportional to the intensity of predation. The population dynamics f (u) is sketched in Figure 2.20(a) when the parameters are in the parameter domain giving three positive steady states ui, u2 and u3, the first and third being linearly stable and the second unstable. With F(u) defined by (2.82) and substituted into (2.84) we have u m as a function of the domain size L. This was evaluated numerically by Ludwig et al. (1979) the form of which is shown in Figure 2.20(b); there is another critical length, L0 say, such that for L > L0 more than one solution exists. We analyse this phenomenon below.

From an ecological viewpoint we would like to know the critical domain size L0 when the maximum population can be in the outbreak regime; that is, um > u2 in Figure 2.20(a). This is determined from numerical integration of (2.84) and is shown in Figure 2.20(b). When L > L0 we see from Figure 2.20(b) that there are three possible solutions with different um. The ones with u in the refuge and outbreak regimes are Figure 2.20. (a) Typical dynamics f (u; r, q) for the spruce budworm as defined by (2.85). (b) The maximum population um as a function of the domain size L. For um < u 1 the population is in the refuge range, whereas um > u2 for L > Lq, which is in the outbreak regime.

Figure 2.20. (a) Typical dynamics f (u; r, q) for the spruce budworm as defined by (2.85). (b) The maximum population um as a function of the domain size L. For um < u 1 the population is in the refuge range, whereas um > u2 for L > Lq, which is in the outbreak regime.

stable and the other, the middle one, is unstable. Which solution is obtained depends on the initial conditions. Later we shall consider possible ecological uses of this model in the control of the budworm. Before doing so we describe a useful technique for determining approximate values for L 0 analytically.

Analytical Method for Determining Critical Domain Sizes and Maximum Populations

The numerical evaluation of um (L) when there are three possible um 's for a given L is not completely trivial. Since the critical domain size L0, which sustains an outbreak, is one of the important and useful quantities we require for practical applications, we now derive an ad hoc analytical method for obtaining it by exploiting an idea described by Lions (1982).

The steady state problem is defined by (2.81). Let us rescale the problem so that the domain is x e (0,1) by setting x ^ x /L so that the equivalent U (x) is now determined from

From Figure 2.19 the solution looks qualitatively like a sine. With the rescaling so that x e (0,1) the solution is thus qualitatively like sin(nx). This means that U" ~ —n U and so the last equation implies

2 2 Dn2U

We are interested in the value of L such that the last equation has three roots for U: this corresponds to the situation in Figure 2.20(b) when L > L0. Thus all we need do to determine an approximate L0 is simply to plot the last equation as in Figure 2.21 and determine the value L such that three solutions exist.

For a fixed dispersal coefficient D we see how the solutions U vary with L .As L increases from L ~ 0 the first critical L, Lc, is given when the straight line Dn2U/L2 intersects f (U), that is, when Dn2/L2 = f'(00), as given by (2.80). As L increases further we can determine the critical L0 when Dn2U/L0 is tangent to the curve f (U), at P in Figure 2.21. It is simply a matter of determining L which gives a double positive root of

It is left as an exercise (Exercise 7) to determine L0 as a function of r, q and D when f (U) is given by (2.85). For any given L the procedure also determines, approximately, the maximum U. From Figure 2.21 we clearly obtain by this procedure a figure similar to that in Figure 2.20(b). This simple procedure is quite general for determining critical domain sizes, both for structure bifurcating from the zero steady state and for domains which can sustain larger populations arising from population dynamics with multiple positive steady states. Figure 2.21. Approximate analytic procedure for determining the critical domain sizes Lc and Lo which can sustain respectively a refuge and an outbreak in the species population where the dynamics is described by f (U). Lc is the value of L when Dn2U/L2 is tangent to f (U) at U = 0. Lo is given by the value of L when DnU/L2 is just tangent to F(U) at P.

Figure 2.21. Approximate analytic procedure for determining the critical domain sizes Lc and Lo which can sustain respectively a refuge and an outbreak in the species population where the dynamics is described by f (U). Lc is the value of L when Dn2U/L2 is tangent to f (U) at U = 0. Lo is given by the value of L when DnU/L2 is just tangent to F(U) at P.

2.8 Spatial Patterns in Scalar Population Interaction Diffusion Equations with Convection: Ecological Control Strategies

In practical applications of such models the domains of interest are usually two-dimensional and so we must consider (2.76). Also, with insect pests in mind, the exterior region is not generally completely hostile, so u = o on the boundaries is too restrictive a condition. Here we briefly consider a one- and two-dimensional problem in which the exterior domain is not completely hostile and there is a prevailing wind. This is common in many insect dispersal situations and can modify the spatial distribution of the population in a major way.

Suppose, for algebraic simplicity, that the two-dimensional domain is a rectangular region B defined byo < x < a, o < y < b having area A. The completely hostile problem is then given by id 2 u d 2 u \

Following the same procedure as in the last section for u small we get the solution of the linearised problem to be i(x, y, t) = O-mn exp-

So the critical domain size, which involves both a and b, is given by any combination of a and b such that a a2b2 Dn2

Since

we get an inequality estimate for spatial patterning to exist, namely,

Estimates for general two-dimensional domains were obtained by Murray and Sperb (l983). Clearly the mathematical problem is that of finding the smallest eigenvalue for the spatial domain considered.

In all the scalar models considered above the spatial patterns obtained have only a single maximum. With completely hostile boundary conditions these are the only type of patterns that can be generated. With two-species reaction diffusion systems, however, we saw that more diverse patterns could be generated. It is natural to ask whether there are ways in which similar multi-peak patterns could be obtained with single-species models in a one-dimensional context. We now show how such patterns could occur.

Suppose there is a constant prevailing wind w which contributes a convective flux (w ■ V)u to the conservation equation for the population u(r, t). Also suppose that the exterior environment is not completely hostile in which case appropriate boundary conditions are

where n is the unit normal to the domain boundary d B. The parameter h is a measure of the hostility: h = to implies a completely hostile exterior, whereas h = 0 implies a closed environment, that is, zero flux boundaries. We briefly consider the latter case later. The mathematical problem is thus ut + (w ■ V)u = f (u) + DV2 u, (2.92)

with boundary conditions (2.9l) and given initial distribution u(r, 0). Here we consider the one-dimensional problem and follow the analysis of Murray and Sperb (l983), who also deal with the two-dimensional analogue and more general aspects of such problems.

The problem we briefly consider is the one-dimensional system which defines the steady state spatially inhomogeneous solutions U(x). From (2.9l) and (2.92), since

(n ■ V)u + hu = 0 ^ ux + hu = 0, x = L; ux — hu = 0, x = 0, where w1 is the x-component of the wind w, the mathematical problem for U (x) is

DU" - W1U' + f (U) = 0, U'(0) - hU(0) = 0, U'(L) + hU(L) = 0.

We study the problem using phase plane analysis by setting

dV w1 V - f (U) U'= V, DV' = wiV - f (U) ^ — = 1 DVf , (2.94)

and we look for phase plane trajectories which, from the boundary conditions in (2.93), join any point on one of the following lines to any point on the other line,

The phase plane situation is illustrated in Figures 2.22(a) and (b) as we now show.

Refer first to Figure 2.22(a). From (2.94) we get the sign of dV/dU at any point (U, V). On the curve V = f (U)/wi, dV/dU = 0 with dV/dU positive and negative when (U, V) lies respectively above (if V > 0) and below it. So, if we start on the boundary line V = hU at say, P, the trajectory will qualitatively be like T1 Figure 2.22. (a) With h sufficiently large the possible trajectories from V = hU to V = -hU admit solution trajectories like Ti and T2 with only a single maximum Um. (b) For small enough h it is possible to have more complex patterns as indicated by the specimen trajectory T. (c) A typical solution U (x) for a phase trajectory like T in (b).

Figure 2.22. (a) With h sufficiently large the possible trajectories from V = hU to V = -hU admit solution trajectories like Ti and T2 with only a single maximum Um. (b) For small enough h it is possible to have more complex patterns as indicated by the specimen trajectory T. (c) A typical solution U (x) for a phase trajectory like T in (b).

since dV/dU < 0 everywhere on it. If we start at S, say, although the trajectory starts with dV/dU < 0, it intersects the dV/dU = 0 line and passes through to the region where dV/dU > 0 and so the trajectory turns up. The trajectories T2, T3 and T4 are all possible scenarios depending on the parameters and where the solution trajectory starts. T3 and T4 are not solution trajectories satisfying (2.94) since they do not terminate on the boundary curve V = -hU. T1 and T2 are allowable solution paths and each has a single maximum U where the trajectory crosses the V = 0 axis.

We now have to relate the corresponding domain length L to these solution trajectories. To be specific let us focus on the trajectory T2. Denote the part of the solution with V > 0 by V+(U) and that with V < 0 by V-(U). If we now integrate the first equation in (2.94) from Uq to UQ, that is, the U-values at either end of the T2 trajectory, we get the corresponding length of the domain for the solution represented by T2 as r Um r U Q

UQ Um

So, for each allowable solution trajectory we can obtain the corresponding size of the solution domain. The qualitative form of the solution U (x) as a function of x can be deduced from the phase trajectory since we know U and U' everywhere on it and from the last equation we can calculate the domain size. With the situation represented by Figure 2.22(a) there can only be a single maximum in U (x). Because of the wind convection term, however, there is no longer the solution symmetry of the solutions as in the last Section 2.7.

Now suppose the exterior hostility decreases, that is, h in (2.95) decreases, so that the boundary lines are now as illustrated in Figure 2.22(b). Proceeding in the same way as for the solution trajectories in Figure 2.22(a) we see that it is possible for a solution to exist corresponding to the trajectory T. On sketching the corresponding solution U(x) we see that here there are two maxima in the domain: see Figure 2.22(c). In this situation, however, we are in fact patching several possible solutions together. Referring to Figure 2.22(b) we see that a possible solution is represented by that part of the trajectory T from A to B1. It has a single maximum and a domain length L1 given by the equivalent of (2.96). So if we restrict the domain size to be L1 this is the relevant solution. However, if we allow a larger L the continuation from B1 to B2 is now possible and so the trajectory AB1 B2 corresponds to a solution of (2.93). Increasing L further we can include the rest of the trajectory to B3. It is thus possible to have multi-humped solutions if the domain is large enough. The length L corresponding to the solution path T is obtained in exactly the same way as above, using the equivalent of (2.96).

So, for small enough values of h it is possible to have more and more structure as the trajectory winds round the point u2 in the (U, V) phase plane. For such solutions to exist, of course, it is essential that w1 = 0. If w1 = 0 the solutions are symmetric about the U-axis and so no spiral solutions are possible. Thus, a prevailing wind is essential for complex patterning. It also affects the critical domain size for patterns to exist. General results and further analysis are given by Murray and Sperb (1983).

### An Insect Pest Control Strategy

Consider now the problem of insect pest control. The forest budworm problem is very much a two-dimensional spatial problem. As we pointed out in Chapter 1, Volume I, Section 1.2, a good control strategy would be to maintain the population at a refuge level. As we also showed in Section 1.2, Volume I it would be strategically advantageous if the dynamics parameters r and q in (2.85) could be changed so that only a single positive steady state exists. This is not really ecologically feasible. With the more realistic spatial problem, however, we have a further possible means of keeping the pest levels within the refuge range by ensuring that their spatial domains are of a size that does not permit populations in the outbreak regime. The arguments go through for two-dimensional domains, but for illustrative purposes let us consider first the one-dimensional situation.

Refer to Figure 2.20(b). If the spatial region were divided up into regions with size L < L0, that is, so that the maximum um was always less than u 1, the refuge population level, we would have achieved our goal. So, a possible strategy is to spray the region in strips so that the non-sprayed regions impose an effective L < L0 as in Figure 2.23(a): the solid vertical lines separating the sprayed regions are the boundaries to a completely hostile exterior.

Of course it is not practical to destroy all pests that stray out of the unsprayed region, so a more realistic model is that with boundary conditions (2.91) where some insects can survive outside the untreated domain. The key mathematical problem to be solved then is the determination of the critical width of the insect 'break' Lb. This must be such that the contributions from neighbouring untreated areas do not contribute a sufficient number of insects, which diffuse through the break, to initiate an outbreak in the neighbouring patches even though L < L0, the critical size in isolation. A qualitative population distribution would typically be as shown by the dashed line in Figure 2.23(a).

The two-dimensional analogue is clear but the solution of the optimisation problem is more complicated. First the critical domain A0 which can sustain an insect pest Figure 2.23. (a) A possible control strategy to contain the insect pest in a refuge rather than an outbreak environment. Strips—insect 'breaks'—are sprayed to maintain an effective domain size L < L0, the critical size for an outbreak. The broken line is a more typical situation in practice. (b) Equivalent two-dimensional analogue where A > Ac is a typical domain which can sustain a pest refuge population but which is not sufficient to sustain an outbreak; that is, A < A0.

Figure 2.23. (a) A possible control strategy to contain the insect pest in a refuge rather than an outbreak environment. Strips—insect 'breaks'—are sprayed to maintain an effective domain size L < L0, the critical size for an outbreak. The broken line is a more typical situation in practice. (b) Equivalent two-dimensional analogue where A > Ac is a typical domain which can sustain a pest refuge population but which is not sufficient to sustain an outbreak; that is, A < A0.

outbreak has to be determined for boundary conditions (2.91). Then the width of the sprayed strips has to be determined. It is not a trivial problem to solve, but certainly a possible one. A preliminary investigation of these problems was carried out by Ben-Yu et al. (1986).

Although we have concentrated on the budworm problem the techniques and control strategies are equally applicable to other insect pests. The field of insect dispersal presents some very important ecological problems, such as the control of killer bees now sweeping up through the western United States (see, for example, Taylor 1977) and locust plagues in Africa. Levin (see, for example, 1981a,b) has made realistic and practical studies of these and other problems associated with spatially heterogeneous ecological models. The concept of a break control strategy to prevent the spatial spread of a disease epidemic will be discussed in some detail later in Chapter 13 when we discuss the spatial spread of rabies; it is now in use.

In an interesting ecological application of diffusion-driven instability Hastings et al. (1997) investigated an outbreak of the western tussock moth which had been hypothesised to be the result of a predator-prey interaction between the moths and parasitoids. The model they analyse qualitatively is a two-species system in which the prey do not move. We saw in Chapter 13, Section 13.7, Volume I how the effect of having a percentage of the prey sessile gave rise to counterintuitive results. Hastings et al. (1997) also obtain counterintuitive results by considering a quite general system with typical predator-prey interactions in which the prey does not diffuse. Their new analysis is very simple but highly illuminating. Their approach is reminiscent of that in Section 1.6 on excitable waves for determining where the trajectory for the wave in the phase plane is closed, thereby determining the wavespeed among other things, and revolves around the existence of three possible steady states at each spatial position, two of which are stable, with the possibility of a jump, or rather a steep singular region joining them. They then show that, counterintuitively, the spatial distribution of the prey will have its highest density at the edge of the outbreak domain. This phenomenon has been observed in western tussock moth outbreaks. The role of theory in predicting counterintuitive behaviours and subsequent experimental or observational confirmation is particularly important when, as is frequently the case, the plethora of observational facts is bewildering rather than illuminating. The article by Kareiva (1990) is particularly relevant to this relation of theory to data.

2.9 Nonexistence of Spatial Patterns in Reaction Diffusion Systems: General and Particular Results

The scalar one-dimensional reaction diffusion system with zero flux is typically of the form ut = f (u) + uxx, x e (0,1), t > 0

Intuitively the only stable solution is the spatially homogeneous one u = u0, the steady state solution of f (u) = 0: if there is more than one stable solution of f (u) = 0

which will obtain depends on the initial conditions. It can be proved that any spatially nonuniform steady state solution is unstable (the analysis is given in the first edition of this book: it involves estimating eigenvalues). This result does not carry over completely to scalar equations in more than one space dimension as has been shown by Matano (l979) in the case where f (u) has two linearly stable steady states. The spatial patterns that can be obtained, however, depend on specific domain boundaries, non-convex to be specific. For example, a dumbbell shaped domain with a sufficiently narrow neck is an example. The pattern depends on the difficulty of diffusionally transporting enough flux of material through the neck to effect a change from one steady state to another so as to achieve homogeneity.

We saw in Sections 2.3 and 2.4 how reaction diffusion systems with zero flux boundary conditions could generate a rich spectrum of spatial patterns if the parameters and kinetics satisfied appropriate conditions: crucially the diffusion coefficients of the reactants had to be different. Here we show that for general multi-species systems patterning can be destroyed if the diffusion is sufficiently large. This is intuitively what we might expect, but it is not obvious if the diffusion coefficients are unequal. This we now prove. The analysis, as we show, gives another condition involving the kinetic relaxation time of the mechanism which is certainly not immediately obvious.

Before discussing the multi-species multi-dimensional theory it is pedagogically helpful to consider first the general one-dimensional two-species reaction diffusion system

Vt = g(u, V) + D2Vxx with zero flux boundary conditions and initial conditions

V(x, 0) = Vo(x), where uQ(x) and vQ(x) are zero on x = 0, l. Define an energy integral E by

This is, except for the l/2, the heterogeneity function introduced in (2.52). Differentiate E with respect to t to get dE _ r dt Jo dE _ I (UxUxt + Vx Vxt) dx o and substitute from (2.98), on differentiating with respect to x, to get, on integrating by parts, dE dt

Because of the zero flux conditions the integrated terms are zero. Now define the quantities d and m by d = min(Di, D2), m = max (f2 + f2 + g2u + gf) , (2.101)

where maxu,v means the maximum for u and v taking all possible solution values. If we want we could define m by some norm involving the derivatives of f and g; it is not crucial for our result. From the equation for dE /dt, with these definitions, we then have dE /*1 2 2 f1 2 2

— <—d (uxx + v2xx) dx + 4m (ui + vt) dx dt Jo Jo (2.102)

Jo Jo with a similar inequality for v; see Appendix A for a derivation of (2.103).

From the inequality (2.102) we now see that if the minimum diffusion coefficient d, from (2.101), is large enough so that (4m — 2n2d) < 0 then dE/dt < 0, which implies that E ^ 0 as t ^ rn since E (t) > 0. This implies, with the definition of E from (2.100), that ux ^ 0 and vx ^ 0 which implies spatial homogeneity in the solutions u and v as t ^ to. The result is not precise since there are many appropriate choices for m; (2.101) is just one example. The purpose of the result is simply to show that it is possible for diffusion to dampen all spatial heterogeneities. We comment briefly on the biological implication of this result below.

We now prove the analogous result for general reaction diffusion systems. Consider ut = f(u) + DV2 u, (2.104)

where u, with components Ui, i = 1, 2,... , n, is the vector of concentrations or populations, and D is a diagonal matrix of the positive diffusion coefficients Di, i = 1, 2,... ,n and f is the nonlinear kinetics. The results we prove are also valid for a diffusion matrix with certain cross-diffusion terms, but for simplicity here we only deal with (2.104). Zero flux boundary and initial conditions for u are

where n is the unit outward normal to d B, the boundary of the domain B. As before we assume that all solutions u are bounded for all t > 0. Practically this is effectively assured if a confined set exists for the reaction kinetics.

We now generalise the previous analysis; it helps to refer to the equivalent steps in the above. Define the energy E (t ) by

where the norm

Let d be the smallest eigenvalue of the matrix D, which in the case of a diagonal matrix is simply the smallest diffusion coefficient of all the species. Now define m = u ||Vuf(u) ||, (2.107)

max where u takes on all possible solution values and Vu is the gradient operator with respect to u.

Differentiating E (t) in (2.106), using integration by parts, the boundary conditions (2.105) and the original system (2.104) we get, with (a, b) denoting the inner product of a and b, dE

J (V u, DV2 u) d r - J (V2 u, DV2u) d r + / (V u, Vuf -Vu) d r JdB JB JB 