## O

which has nontrivial solutions if and only if the determinant of the coefficient matrix, | XI + D| k|2 - A |= 0 (recall Chapter 2). So, the dispersion relation X(k2), k2 = | k |2 is given by the characteristic equation

+ dudv k4 - (dugU + dv fU + au* X * g'U)k2 + fU gU - ft gU

with solutions denoted by A.+ and X-.

We are interested in determining pattern modes which have at least one positive growth rate; that is, at least one solution has RlX > 0 so we focus on the larger of the two solutions, X+, given by

c(k2) = dudvk4 - (dugV + dv fV + auvxvgV)k2 + fVgV - fvgV. (5.54)

Recalling the discussion in Chapter 2, if RlX+ is positive (negative) then perturbations about the steady state will grow (decay). We look for solutions which are stable (A.+ < 0) to purely temporal perturbations (that is, k2 = 0), but unstable (A.+ > 0) to at least one spatial mode, which is just like the pure diffusionally driven instability situation for some nonzero k1 and k2. Mathematically, we look for a range of parameters such that

A+(0) < 0, X+(k2) > 0 for all k such that 0 < k2 < k2 < kf. (5.55)

To satisfy the first of (5.55) we must have b(0) = fV+ gV < 0, c(0) = fVV gV- fV gV > 0. (5.56)

Note that these conditions imply that is always negative so we need only focus on A.+ . Graphically, these conditions yield an inverted parabolic curve for X(k2), which has its maximum to the right of k2 = 0 and so is the most basic dispersion relation which gives diffusion-chemotaxis-driven instability discussed in detail in Chapter 2. Typical dispersion relations are shown in Figure 5.11, where by way of example we have shown how they vary with the parameter for / small enough there is no range k2, k2.

Figure 5.11. The dispersion relation, X(k2) for parameter values dv = 1.0, du = 0.3, a = 80, f = 2, / = 4, 6, 8 and 10, S = 2, p = 1 and w = 10. The curve corresponding to / = 4 is the lowest one.

We are interested in the intersection of the curve with X = 0 which gives k12 and k|. These correspond to a boundary in parameter space, which from (5.54) is given by the two solutions k2 of dudvk4 - (dug*v + dvf* + au*x*gU)k2 + f*g* - f*g* = 0. (5.57)

In general, for each set of parameter values there are two, one or zero values of k2 which satisfy (5.57). At bifurcation, where there is one value, we have k2 dug*+ dv f*v+ auvxvg* (5 58) kc =-2dd- (5.58)

(ddugv + dv fU + au*Xvg*)2 - 4dudv( fVg* - fVg*) = 0, (5.59)

where kc is the critical wavenumber. Solving the last equation for a and substituting it into (5.58) we get the crucial values a _ -(dugt + dv/*) + 2^dudv(g* f* - g*u f*) k2 = (f u*g*v - f*g*)

which give the parameter spaces for diffusion-chemotaxis-driven spatial instability in terms of the other parameters in the model system via the functions x , f and g, here defined by (5.46), and their derivatives evaluated at the steady state.

For these results the positive solutions to the quadratic were chosen. Equation (5.60) defines a critical boundary set in parameter space which separates the two regions of positive and negative X. The curve (5.59) is the bifurcation where X+ = 0 and wave patterns will neither grow nor decay. A sequence of these curves is shown in Figure 5.12 for various values of On the upper side of each curve, X+ > 0 and all unstable pattern modes will grow, that is, for all k2 < k2 < k| obtained from (5.57). We do not know at this point what patterns to expect, only that spatial patterns are possible. To go further and determine which patterns will emerge is a nonlinear problem. To date the only analytical way to determine these is by what is called a weakly nonlinear analysis which means an anlysis near where the solutions bifurcate from spatial homogeneity to spatial heterogeneity. Some of the references where this has been done for reaction diffusion equations were given in Chapter 2. Exactly the same procedures are used in these types of equations but are just a little more complex. Zhu and Murray (1995) carried out a nonlinear analysis and evaluated, analytically, the parameter spaces for both reaction diffusion and diffusion-chemotaxis systems to compare their potential for generating spatial pattern. As already mentioned, Tyson (1996) used the method to analyze the specific model equations we have discussed here.

We do not carry out the (rather complicated) nonlinear analysis here (refer to Zhu and Murray 1995 and Tyson 1996 for a full discussion). Here we just sketch the proce-

Pattcrn/no pattern boundary

Pattcrn/no pattern boundary

Figure 5.12. The boundary in (ac, fc) parameter space between regions where patterns are possible and where they are not. For (a, f) pairs above each line patterns are possible (I > 0) while below, they are not. The parameter values are duc = 0.3, dvc = 1.0, /j,c = 4 (uppermost line), to 8 (lowermost line) in steps of 0.5, Sc = 2, pc = 0.5 and wc = 5. (From Tyson 1996)

Figure 5.12. The boundary in (ac, fc) parameter space between regions where patterns are possible and where they are not. For (a, f) pairs above each line patterns are possible (I > 0) while below, they are not. The parameter values are duc = 0.3, dvc = 1.0, /j,c = 4 (uppermost line), to 8 (lowermost line) in steps of 0.5, Sc = 2, pc = 0.5 and wc = 5. (From Tyson 1996)

dure and give the results. From the linear analysis we first determine in what parameter regions patterns are possible. The nonlinear analysis gives, in effect, the type of patterns which will form. The analysis is based on the assumption that the parameters are such that we are close to the bifurcation curve in parameter space. We develop an asymptotic analysis based on one parameter (the one of specific interest), for example, the chemotaxis parameter a, which is close to its bifurcation value, ac, on the bifurcation curve but whose value moves the system into the pattern formation space as it passes through its critical value. We start with the linear solution to the boundary value problem; for example, in one dimension it involves a solution including ex(t) cos kx, where k is in the unstable range of wavenumbers. On a linear basis this solution will start to grow exponentially with time with a growth rate ek(k)t with X(k) the dispersion relation. Intuitively, by examining the undifferentiated terms in the system of equations (5.38)— (5.40), the solutions can not grow unboundedly. In a weakly nonlinear analysis we first consider the linear solution to be the solution to the linear boundary value problem with k2 = kc2; that is, the parameters are such that we are close to the bifurcation curve from no pattern to pattern. An asymptotic perturbation method is then used to study the solutions in the form where the magnitude (or amplitude) of the solution is a slowly varying function of time. The conditions under which the amplitude is bounded as t ^to are determined. This procedure then determines which of the various solution possibilities will evolve as stable solutions. We give a few more details to further explain the general procedure when we discuss the forms of the possible linear solutions to the boundary value problem.

### Linear Boundary Value Problem

A necessary prerequisite for the nonlinear analysis is solving the linear boundary value problem relevant to the nonlinear analysis. The pattern types which are possible depend on the number of different wavevectors k allowed by the boundary conditions. From our point of view we are interested in having the experimental domain tesselated by repeating patterns. In the numerical simulations a square domain was chosen for numerical simplicity and so we are interested in square or rectangular tiles giving stripes and spots. So, each square or rectangular unit is subject to periodic boundary conditions. In general, with regular tesselations, we can have squares, stripes, hexagons, and so on as we discussed in Chapter 2.

Consider the rectangular domain, S, defined by 0 < x < lx, 0 < y < ly with the sides denoted by Si:x = 0, 0 < y < ly, S2:x = lx, 0 < y < ly, S3:y = 0, 0 < x < lx, S4:y = ly, 0 < x < lx. The spatial eigenvalue problem with periodic boundary conditions is then

The possible eigensolutions of the partial differential equation (5.61) are f = Acos(k„ ■ x) + Bsin(k„ ■ x),

where the k2 are allowable eigenvectors, which we discuss below. Substituting the boundary conditions in (5.61) into the solutions (5.62) we obtain cos(k„ ■ (0, y)) = cos(k„ ■ (lx, y)) sin(k„ ■ (x, 0)) = cos(k„ ■ (x, lx)) cos(k„ ■ (0, y)) = cos(k„ ■ (ly, y)) sin(k„ ■ (x, 0)) = cos(k„ ■ (x, ly))

which after using some trigonometric identities can be written as cos(kyy) [1 - cos(k%lx^ + sin(k;^y) sin(k%lx) = 0 sin(kly) [1 - cos(k^lx^ - cosky) sin(k%lx) = 0 cos(kxx) [1 — cos(knly^ + sin(k%x) sin(kyly) = 0 sin(k%x) [1 — cos(kynly^ — cos(k£x) sin(k1ly) = 0,

where kj = (k%, kn). The general solutions of (5.63) are k1 =

0 0