7  Spatially-structured populations

We now want to keep track of how individuals are distributed over space. So instead of describing a population by its total number of individuals \(N(t)\) we describe it by a population density \(u(x,t)\), where \(u(x,t)dx\) is the number of individuals in the interval \([x,x+dx]\) at time \(t\). We can recover the total number of individuals by integrating the density over the entire spatial domain.

7.1 Derivation of PDE

In the non-spatial case we described the time evolution of the population by the equation \[ \frac{dN}{dt} = f(N), \tag{7.1}\] where \(f(N)\) is the net growth rate that encodes the difference between birth and death rates for the population as a whole. In the spatial case we have to consider the net growth rate of the population in a small interval \([x,x+dx]\). This is again given by a function \(f\) that encodes the births and deaths, but now we also have movement of individuals into or out of the interval. We write the rate of change of the number of individuals in the interval as \[ \frac{\partial}{\partial t}\int_{x_0}^{x_0+\Delta x}u(x,t)dx = \int_{x_0}^{x_0+\Delta x}f(u(x,t),x)dx +J(x_0)-J(x_0+\Delta x). \tag{7.2}\] The flux \(J(x)\) is defined as the net rate at which individuals move through point \(x\) from left to right. If more individuals move from right to left than from left to right, the flux is negative. We have indicated explicitly that \(f\) can depend on both the population density \(u\) at \(x\) as well as the position \(x\) itself. This is because the birth and death rates can depend on the local environment, for example the availability of resources or the presence of predators. The flux \(J\) can also depend on the local population density and the position. We did not indicate this explicitly to simplify the notation.

We now use the integral mean value theorem that states that for a continuous function \(g(x)\) the integral over \(g(x)\) over an interval \([a,b]\) is equal to \(g(\xi)(b-a)\) for some \(\xi\in[a,b]\). We apply this to the integrals in Eq. 7.2 to get \[ \partial_t u(\xi_1,t)\Delta x = f(u(\xi_2,t),\xi_2)\Delta x +J(x_0)-J(x_0+\Delta x) \tag{7.3}\] for some \(\xi_1,\xi_2\in[x_0,x_0+\Delta x]\). We now divide both sides by \(\Delta x\), \[ \partial_t u(\xi_1,t) = f(u(\xi_2,t),\xi_2) +\frac{J(x_0)-J(x_0+\Delta x)}{\Delta x} \tag{7.4}\] and take the limit \(\Delta x\to 0\), where \(\xi_1\to x\) and \(\xi_2\to x\) and the difference quotient becomes the derivative, to get \[ \partial_t u(x,t) = f(u(x,t),x)-\partial_x J(x). \tag{7.5}\]

As discussed above, the flux \(J\) can depend on \(x\) both directly and through the population density \(u(x,t)\). An important example is the case where the flux is proportional to the gradient of the population density, \[ J(x) = -D\partial_x u(x), \tag{7.6}\] where \(D\) is the diffusion coefficient. This models random motion of the individuals. Each individual is equally likely to move right or left, but if there are more individuals on the left and less on the right, then the result is a net movement to the right. That is why the flux has the opposite sign of the gradient of the density.

This random motion gives us the reaction-diffusion equation \[ \partial_t u=f(u)+D\,\partial_x^2 u. \tag{7.7}\] If there is no local population dynamics, the equation simplifies to the diffusion equation, also known as the heat equation. The diffusion term has the effect of smoothing out spatial inhomogeneities in the population density. Therefore the heat equation by itself is rather boring. However, as we will see later, the reaction term in combination with the diffusion term can lead to the formation of spatial patterns, such as travelling waves or stationary patterns.

7.2 Fishing model with diffusion

We now consider a model for a population of fish that is subject to fishing and that moves around randomly. We want to model a marine protected area where fishing is limited to avoid a collapse of the fish population. Figure 7.1 shows the spatial setup. The \(x\)-axis runs perpendicular to the shore, which is at \(x=0\). A marine protected area runs out upto a distance \(L\) from the shore. We assume that the marine protected area has infinite extension in the \(y\)-direction, so that we can ignore the \(y\)-coordinate and model the fish population as a function of \(x\) only.

Beyond \(x=L\) there is no management of the fishing activity. We assume that the fishers are so efficient that they catch all fish that venture outside the protected area. So we set \(u(x,t)=0\) for \(x\geq L\). This is of course an idealisation, but it allows us to focus on the dynamics within the protected area.

Figure 7.1: Sketch of the marine protected area

We assume that in the absence of fishing the population has a logistic growth rate. The population density \(u(x,t)\) satisfies the equation \[ \partial_t u = r u \left(1-\frac{u}{K}\right) - E u + D \partial_x^2 u, \tag{7.8}\] where \(r\) is the intrinsic growth rate, \(K\) is the carrying capacity, \(E\) is the fishing rate and \(D\) is the diffusion coefficient. The first term on the right-hand side describes the logistic growth, the second term the fishing and the third term the random motion of the fish.

In the protected area, the fishing rate \(E\) must certainly be below the intrinsic growth rate \(r\), otherwise the population is guaranteed to die out. However we must also take into account that some fish will be lost when they randomly move out of the protected area and are then fished immediately. We expect this extra loss to be proportionally larger when the protected area is smaller. So when the marine protected area is planned, there will be a trade-off between the width \(L\) of the protected area and the allowed rate of fishing \(E\). Our aim is to describe this trade-off precisely by finding the condition on \(E\) and \(L\) such that the population does not collapse.

To achieve our aim, we will need to solve the PDE. But before solving a PDE we always need to be clear about the boundary conditions on the solution. We have already discussed that we assume that \(u(x,t)=0\) for \(x\geq L\), i.e., we have a Dirichlet boundary condition at \(x=L\). At \(x=0\) we impose a no-flux boundary condition, which means that there are no fish moving from the sea onto the shore or the other way around. The no-flux boundary condition is the Neumann boundary condition \(\partial_x u(0,t) = 0\).

The PDE in Eq. 7.81 is a nonlinear PDE, so we cannot solve it analytically in general. However, we can still make progress by realising that before the population goes extinct, it will get small. We can then linearise the PDE around the extinction state \(u=0\) and solve the linearised equation. If the extinction steady state is linearly unstable, then the population will not go extinct.

Linearising the PDE Eq. 7.81 around the extinction state \(u=0\) means that we neglect the term \(u^2\) in the logistic growth term. We then have the linearised equation \[ \partial_t u = (r - E)u + D \partial_x^2 u. \tag{7.9}\] We can now solve this linear PDE with the method of separation of variables. We make the Ansatz \(u(x,t) = X(x)T(t)\) and plug this into Eq. 7.9. We get \[ XT'=(r-E)XT+DX''T \tag{7.10}\] and after dividing by \(XT\) we get \[ \frac{T'}{T} = r-E + D \frac{X''}{X}. \tag{7.11}\] The left-hand side depends only on \(t\) and the right-hand side only on \(x\). Therefore both sides must be equal to a constant, which we call \(\gamma\). We then have the two ODEs \[ T' = \gamma T, ~~~ X'' = \frac{\gamma-r+E}{D}X =-\rho^2 X, \tag{7.12}\] where, in order to save writing, we have introduced \[ \rho^2 = (r-E-\gamma)/D. \tag{7.13}\] The solution of the time ODE is \[ T(t) = T(0)e^{\gamma t} \tag{7.14}\] and the solution of the spatial ODE is \[ X(x) = A \cos(\rho x) + B \sin(\rho x), \tag{7.15}\] where \(A\) and \(B\) are constants to be determined by the boundary conditions.

Imposing the no-flux boundary condition \(\partial_xu(0,t)\) gives \[ X'(0) =-A\rho\sin(\rho 0)+B\rho\cos(\rho 0)=B\rho= 0. \tag{7.16}\] So either \(\rho=0\) or \(B=0\). The case \(\rho=0\) gives the trivial solution \(X(x)=A\), which is not interesting. So we take \(B=0\) and get \(X(x) = A\cos(\rho x)\).

The Dirichlet boundary condition at \(x=L\) gives \(X(L) = A\cos(\rho L) = 0\). So either \(A=0\) or \(\cos(\rho L) = 0\). The former gives rise to the zero solution, so we want the latter, which requires \(\rho L\) to be an odd multiple of \(\pi/2\), i.e., \[ \rho=\rho_n = \frac{(2n+1)\pi}{2L} \tag{7.17}\] for \(n=0,1,2,...\). To each of these \(\rho_n\) there corresponds, according to Eq. 7.13, a \[ \gamma_n = r-E-D\rho_n^2 \tag{7.18}\] and a solution \[ u_n(x,t) = X_n(x)T_n(t) = A_ne^{\gamma_n t}\cos(\rho_n x) \tag{7.19}\] The general solution is then a linear combination of these solutions, \[ u(x,t) = \sum_{n=0}^\infty A_ne^{\gamma_n t}\cos(\rho_n x). \tag{7.20}\]

The condition for the population not to go extinct is that the extinction steady state is linearly unstable. This means that there is at least one \(\gamma_n\) with positive real part. We see from Eq. 7.18 that \(\gamma_n\) increases with \(n\). So the condition for the population not to go extinct is that \(\gamma_0\) is positive. This gives the condition \[ \gamma_0=r-E-D\rho_0^2=r-E-\frac{D\pi^2}{4L^2}>0. \tag{7.21}\] We can now solve this inequality for \(L\) to get the condition on the width of the marine protected area that ensures the population does not collapse. \[ L>\frac{\pi}{2}\sqrt{\frac{D}{r-E}}. \tag{7.22}\] Alternatively we can solve the inequality for \(E\) to get the condition on the fishing rate that ensures the population does not collapse. \[ E<r-\frac{D\pi^2}{4L^2}. \tag{7.23}\] This can now be used to inform policy decisions on the width of the marine protected area and the allowed fishing rate.

7.3 Invasion waves in SIR model

We now want to use our spatial modelling skills to study the spread of an infectious disease through space. For concreteness, think of a situation where rabies has infected foxes in Dover and we want to get an idea of how soon we will have rabid foxes in York. We model the fox population by the SIR model, where \(S(x,t)\) is the density of susceptible foxes, \(I(x,t)\) is the density of infected foxes and \(R(x,t)\) is the density of removed foxes (which for rabies unfortunately means dead foxes). So where in Section 6.1 we had ordinary differential equations for the total number of susceptible, infected and removed foxes, we now have partial differential equations for the densities of susceptible, infected and removed foxes in space.

While foxes are usually very territorial and stay in their own territory, infected foxes become a bit insane and move around randomly. This is modelled by a diffusion term in the equation for the infected foxes. The equations are \[\begin{split} \partial_t S&=-\beta \frac{SI}{N},\\ \partial_t I&=\beta \frac{SI}{N}-\gamma I+D\partial_x^2I,\\ \partial_t R&=\gamma I, \end{split} \tag{7.24}\] where \(N=S+I+R\) is the total fox population density, \(\beta\) is the infection rate, \(\gamma\) is the recovery rate and \(D\) is the diffusion coefficient. For simplicity we study the movement in the \(x\) direction only. The above equations are valid for all \(x\). We will later concentrate on the equations for \(I\) and \(S\) only, as \(R\) can then be calculated from \(N=S+I+R\).

Before analysing these equations further, let us think about what we expect to happen. The infected foxes in Dover will infect more and more susceptible foxes. The infected foxes will then move around and infect susceptible foxes in neighbouring territories. This will lead to a wave of infection spreading out from Dover. The wave will move at a speed that depends on the infection rate, the recovery rate and the diffusion coefficient. The wave will have a front where the infected foxes are and a tail where the epidemic has died out.

To make the further analysis easier, we will non-dimensionalise the equations. We introduce the non-dimensional variables \[ u = S/N, ~~~ v = I/N,~~~ \tilde{t} = \beta t, \tag{7.25}\] Then the equations become \[\begin{split} \partial_{\tilde{t}} u&=-uv\\ \partial_{\tilde{t}} v&=uv-\frac{\gamma}{\beta}v+\frac{D}{\beta}\partial_x^2v. \end{split} \tag{7.26}\] We now also introduce the non-dimensional spatial variable \(\tilde{x} = \sqrt{\frac{\beta}{D}}x\) and the non-dimensional parameter \(r = \frac{\gamma}{\beta}\). We now drop the tildes to avoid clutter and write the equations as \[\begin{split} \partial_t u&=-uv\\ \partial_t v&=uv-rv+\partial_x^2v. \end{split} \tag{7.27}\]

7.3.1 Travelling wave Ansatz and boundary conditions

We are looking for a solution describing the spread of the infection from Dover to York. We make the Ansatz that the solution is a travelling wave, i.e., that it is of the form \[ u(x,t)=U(z),~~~v(x,t)=V(z) \tag{7.28}\] with \(z=x-ct\) for some \(c>0\). This means that the wave is a right-moving wave. We plug this into Eq. 7.27 and get the system of ODEs \[\begin{split} -cU'&=-UV,\\ -cV'&=UV-rV+V''. \end{split} \tag{7.29}\] Next let us think about the boundary conditions. We expect that ahead of the wave all the individuals are still susceptible and there are no infecteds yet. So to the far right at \(z=infty\) we have \(u(\infty)=1\) (corresponding to S=N) and \(v(\infty)=0\). Behind the wave the epidemic will have run its course and so \(u(-\infty)=S_\infty/N=:a\) and \(v(-\infty)=0\). In other words, the travelling wave interpolates between the state before an epidemic (at \(z=\infty\)) and the state after an epidemic (at \(z=-\infty\)). We also have that the solution becomes flat as \(z\to\pm\infty\), so in particular \(V'(\pm\infty)=0\).

We now massage the equations a bit to get them into a form that we can integrate more easily. From the first equation we see that \(UV=cU'\) and also that \(V=cU'/U=c(\log U)'\). We plug this into the second equation and get \[ -cV'=cU'-rc(\log U)' + V''. \tag{7.30}\] Now that each term in the equation is a total derivative, we can integrate it by just removing the differentiations: \[ -cV=cU-rc\log U + V'+A, \tag{7.31}\] where \(A\) is a constant of integration. We use the boundary condition at \(z\to\infty\) to determine \(A\): \[ 0=c-rc\log 1 + 0 + A \Rightarrow A=-c. \tag{7.32}\] At \(z\to-\infty\) we get \[ 0=ca-rc\log a + 0 - c \Rightarrow a-1=r\log a. \tag{7.33}\] This equation is equivalent to Eq. 6.14 for \(S_\infty\) in the case where \(S_0=N\). Again this transcendental equation can only be solved numerically.

7.3.2 Wave speed

Next we want to learn about the wave speed \(c\). We do this by linearising the equation around the leading edge of the wave, where \(V\) is very small and \(U\) is close to 1, i.e, \(U=1-\epsilon\) for small \(\epsilon\). Substituting this into the second equation in 7.29 we get \[ -cV' = (1-\epsilon)V - r V + V''. \tag{7.34}\] Because both \(\epsilon\) and \(V\) are small, we can neglect the product \(\epsilon V\) and get \[ -cV' = (1 - r) V + V''. \tag{7.35}\] This is a linear ODE with constant coefficients and can thus be solved with the Ansatz \(V(z) = e^{-\lambda z}\). Substituting this into the ODE and dividing by \(e^{-\lambda z}\) gives \[ \lambda^2 -c\lambda+ 1-r=0. \tag{7.36}\] The solution of this quadratic equation is \[ \lambda = \frac{c\pm\sqrt{c^2-4(1-r)}}{2}. \tag{7.37}\] We need \(\lambda\) to be real so that our solution correctly describes the exponential growth at the start of an epidemic, so we need the discriminant to be non-negative, i.e., \[ c^2-4(1-r)\geq 0 \tag{7.38}\] and thus we get a lower bound on the wave speed: \[ c\geq 2\sqrt{1-r}. \tag{7.39}\]

7.4 Turing instabilities

We now come to a very intriguing question: how can the complicated spatial patterns arise that one can observe in nature. One might have thought that to explain complicated patterns one would need complicated models that encode all that complexity. However, it turns out that patterns can arise from very simple, translation-invariant models. The key idea is that the spatially uniform steady state in a translation-invariant model can be unstable to small perturbations, which then grow and form spatial patterns. This is called a Turing instability.

We consider a system of two species, described by densities \(u(x,t)\) and \(v(x,t)\), that interact with each other locally but also move around randomly. The system is described by the reaction-diffusion equations \[\begin{split} \frac{\partial u}{\partial t}&=f(u,v)+D_1\frac{\partial^2u}{\partial x^2},\\ \frac{\partial v}{\partial t}&=g(u,v)+D_2\frac{\partial^2v}{\partial x^2}, \end{split} \tag{7.40}\] where \(f\) and \(g\) are functions that describe the local dynamics of the two species. The diffusion terms model random motion of the individuals.

For values \((u^*,v^*)\) of the densities that satisfy \(f(u^*,v^*)=g(u^*,v^*)=0\), the spatially homogeneous steady state \(u(x,t)=u^*\), \(v(x,t)=v^*\) is a solution of the system. Assume that in the absence of diffusion, this steady state is stable. Thus without diffusion, any perturbations to the spatially homogeneous steady state will decay. One would then expect that the spatially homogeneous steady state is stable also in the presence of diffusion, given that diffusion has the tendency to spread individuals out away from regions of higher concentration. The result would be the absence of any spatial structure in the solution. However, we will now see that this is not always the case. Rather counter-intuitively, the random motion of the individuals can destabilise the spatially homogeneous steady state and give rise to spatial patterns.

7.4.1 Deriving conditions for Turing instabilities

Mathematics is a great tool for dealing with counter-intuitive phenomena. We just have to do the maths and see what it tells us. So we now derive the equations that describe the time-evolution of small perturbations and solve them and see whether the perturbations grow or not. So we set \[ u(x,t)=u^*+\xi(x,t),~~~v(x,t)=v^*+\eta(x,t) \tag{7.41}\] where \(\xi(x,t)\) and \(\eta(x,t)\) are small perturbations. If we plug this into Eq. 7.40 we get equations for the time evolution of the small perturbations: \[\begin{split} \frac{\partial \xi}{\partial t}&=f(u^*+\xi(x,t),v^*+\eta(x,t))\xi+D_1\frac{\partial^2\xi}{\partial x^2},\\ \frac{\partial \eta}{\partial t}&=g(u^*+\xi(x,t),v^*+\eta(x,t))\xi+D_2\frac{\partial^2\eta}{\partial x^2}. \end{split} \tag{7.42}\] We now use Taylor expansions to linearise the equations. We expand \(f\) and \(g\) around \((u^*,v^*)\) to first order and drop all terms that are higher order in \(\xi\) and/or \(\eta\). This gives the linear PDEs \[\begin{split} \frac{\partial \xi}{\partial t}&=\frac{\partial f}{\partial u}(u^*,v^*)\xi+\frac{\partial f}{\partial v}(u^*,v^*)\eta+D_1\frac{\partial^2\xi}{\partial x^2}\\ &=a_{11}\xi+a_{12}\eta+D_1\frac{\partial^2\xi}{\partial x^2},\\ \frac{\partial \eta}{\partial t}&=\frac{\partial g}{\partial u}(u^*,v^*)\xi+\frac{\partial g}{\partial v}(u^*,v^*)\eta+D_2\frac{\partial^2\eta}{\partial x^2}\\ &=a_{21}\xi+a_{22}\eta+D_2\frac{\partial^2\eta}{\partial x^2}. \end{split} \tag{7.43}\] Those partial derivatives of \(f\) and \(g\) with respect to \(u\) and \(v\) are familiar to us from the stability analysis in non-spatial models as the entries of the Jacobian matrix and we are using the same shorthand notation for them here.

We now make the harmonic wave Ansatz \[\begin{split} \xi(x,t)&=B_1e^{\sigma_k t}\sin(kx+\alpha),\\ \eta(x,t)&=B_2e^{\sigma_k t}\sin(kx+\alpha). \end{split} \tag{7.44}\] During your studies you will have seen this easy way to solve linear PDEs with constant coefficients already several times. Plugging this Ansatz into Eq. 7.43 and dividing by \(e^{\sigma_k t}\sin(kx+\alpha)\) we get \[\begin{split} \sigma_k B_1&=a_{11}B_1+a_{12}B_2-k^2D_1B_1\\ \sigma_k B_2&=a_{21}B_1+a_{22}B_2k^2D_2B_2. \end{split} \tag{7.45}\] We can write this in matrix form as \[ \sigma_k\begin{pmatrix}B_1\\B_2\end{pmatrix}=A(k)\begin{pmatrix}B_1\\B_2\end{pmatrix}. \tag{7.46}\] where \[ A(k)=\begin{pmatrix}a_{11}-D_1k^2&a_{12}\\a_{21}&a_{22}-D_2k^2\end{pmatrix}. \tag{7.47}\] This looks very similar to what we get when we do the stability analysis of a non-spatial model, except that now we have a whole family of matrices \(A(k)\), one for each \(k\). The condition for the spatially homogeneous steady state to be stable is that all the eigenvalues \(\sigma_k\) of all the matrices \(A(k)\) are negative. The spatially homogeneous steady state to be unstable is that there is at least one \(k\) for which the real part of \(\sigma_k\) is positive.

As we know from studying the stability of non-spatial models of two interacting species, an easy way to determine the stability is to look at the signs of the trace and the determinant of the matrix \(A(k)\). For the trace we find \[ \text{tr}(A(k))=a_{11}-k^2D_1+a_{22}-k^2D_2=\text{tr}(A(0))-k^2(D_1+D_2). \tag{7.48}\] We are interested in the case where the steady state is stable in the absence of diffusion and is destabilised by the diffusion. Stability in the absence of diffusion requires \(\text{tr}(A(0))<0\). We see that in that case also \(\text{tr}(A(k))<0\) for all \(k\) because the \(k\)-dependent term \(-k^2(D_1+D_2)\) is negative. So the trace is not yet giving us a hint that there might be an instability. So next we look at the determinant.

For the determinant we find \[\begin{split} \det(A(k))&=(a_{11}-D_1k^2)(a_{22}-D_2k^2)-a_{12}a_{21}\\ &=D_1D_2k^4-(D_1a_{22}+D_2a_{11})k^2+\det(A(0)) \end{split} \tag{7.49}\] Stability in the absence of diffusion tells us that \(\det(A(0))>0\). We see that \(\det(A(k))\) is a quadratic polynomial in \(k^2\) with positive leading coefficient. Its graph is a parabola. This means that the determinant is positive for small \(k\) and for large \(k\) but there may be an intermediate region where the parabola dips below the axis and the determinant is negative in that region. Perturbations with wave numbers in that region would grow over time, making the homogeneous steady state unstable and leading to spatial structure. If that region exists, it is bounded by the values \(k_\pm\) for which the determinant is zero. Using the quadratic formula we find that \[ k_\pm^2=\frac{D_1a_{22}+D_2a_{11}\pm\sqrt{(D_1a_{22}+D_2a_{11})^2-4D_1D_2\det(A(0))}}{2D_1D_2}. \tag{7.50}\] The unstable interval exists if \(k_\pm^2\) are real and positive. They are real if the discriminant under the square root is non-negative. This gives the condition \[ (D_1a_{22}+D_2a_{11})^2>4D_1D_2\det(A(0)). \tag{7.51}\] They are positive if the numerator is positive. This gives the condition \[ D_1a_{22}+D_2a_{11}>0. \tag{7.52}\] The condition for the spatially homogeneous steady state to be unstable is that both Eq. 7.51 and Eq. 7.52 are satisfied.

7.4.2 Understanding the conditions for Turing instabilities

We have now derived the conditions for Turing instabilities in terms of the entries of the Jacobian matrix \(A=A(0)\) of the local dynamics of the two species and the diffusion coefficients: \[\begin{split} a)\qquad&a_{11}+a_{22}<0,\\ b)\qquad&a_{11}a_{22}-a_{12}a_{21}>0,\\ c)\qquad&D_1a_{22}+D_2a_{11}>0,\\ d)\qquad&(a_{11}+a_{22})^2>4(a_{11}a_{22}-a_{12}a_{21}). \end{split} \tag{7.53}\] Conditions a) and b) are the conditions for the spatially homogeneous steady state to be stable in the absence of diffusion. Condition c) and d) are the condition for the spatially homogeneous steady state to be unstable in the presence of diffusion.

We will first look at what these conditions tell us about the sign structure of the Jacobian, which will introduce the concept of activator and inhibitor species. We make three observations:

  1. We observe that condition a) implies that at least one of \(a_{11}\) and \(a_{22}\) is negative. We are free to choose the numbering of our species, so, without loss of generality, let’s say that \(a_{22}<0\). This means that species 2 inhibits its own growth above the steady state. We call species 2 an inhibitor species.

  2. From condition c) we see that if \(a_{22}<0\) then \(a_{11}>0\). This means that species 1 activates its own growth above the steady state. We call species 1 an activator species.

  3. From condition b) we now observe that, because \(a_{11}a_{22}<0\), we must have that \(a_{12}a_{21}<0\).

Taking these three observations together tells us that there are only two possibilities for the sign structure of the Jacobian matrix: \[ \text{sign}(A)=\begin{pmatrix}+&-\\+&-\end{pmatrix}~~~\text{or}~~~\text{sign}(A)=\begin{pmatrix}+&+\\-&-\end{pmatrix}. \tag{7.54}\] If the signs in the Jacobian are different from both of these possibilites you can immediately rule out a Turing instability.

Next we look at what our conditions a) to d) tell us about the magnitude of the diffusion coefficients. From condition a) we know that \(a_{11}<-a_{22}\) and hence \(-a_{11}/a_{22}<1\). From condition c) we see that \(D_1/D_2<-a_{11}/a_{22}\). Putting these together gives \(D_1/D_2<1\) and hence \(D_1<D_2\): The inhibitor must diffuse faster than the activator. This is a necessary condition for a Turing instability.

Exactly how much faster the inhibitor must diffuse than the activator depends on the details of the local interactions and we will next study this in a concrete example.

7.4.3 Example

Let’s consider the following concrete reaction-diffusion system: \[\begin{split} \partial_tu&=-u+u^2v+D_1\partial_x^2u=f(u,v)+D_1\partial_x^2u,\\ \partial_tv&=b-u^2v+D_2\partial_x^2v=g(u,v)+D_2\partial_x^2v. \end{split} \tag{7.55}\]

The Jacobian matrix of the local dynamics is \[ A=\begin{pmatrix}\frac{\partial f}{\partial u}&\frac{\partial f}{\partial v}\\ \frac{\partial g}{\partial u}&\frac{\partial g}{\partial v}\end{pmatrix}= \begin{pmatrix}-1+2uv&u^2\\-2uv&-u^2\end{pmatrix}. \tag{7.56}\] We are interested in the coexistence steady state \((u^*,v^*)\) where both species coexist. Imposing that \(f(u^*,v^*)=0\) gives us that \(u^*v^*=1\) and then imposing that \(g(u^*,v^*)=0\) gives us that \(u^*=b\) and hence \(v^*=1/b\). So the coexistence steady state is \[ (u^*,v^*)=\left(b,\frac{1}{b}\right). \tag{7.57}\] The Jacobian matrix evaluated at the coexistence steady state is \[ A=\begin{pmatrix}1&b^2\\-2&-b^2\end{pmatrix}. \tag{7.58}\] We see that the sign structure of this matrix is one of those in Eq. 7.54 that can give rise to a Turing instability. Species \(u\) is the activator and species \(v\) is the inhibitor.

Now let us check each of the conditions a) to d) from Eq. 7.53. Condition a) becomes \(1-b^2<0\) and hence \(b^2>1\). Condition b) becomes \(-b^2+2b^2>0\) which is automatically satisfied. Condition c) becomes \(D_2-b^2D_1>0\) which implies \(D_1b^2<D_2\), which means that the ratio \(d=D_2/D_1>b^2\). Condition d) becomes \((D_2-b^2D_1)^2>4D_1D_2b^2\). Multiplying out the square and bringing everything to one side gives \(b^4D_1^2-6D_1D_2+D_2^2>0\). We can rewrite this in terms of the ratio \(d=D_2/D_1\) as \[ b^4-6db^2+d^2>0. \tag{7.59}\] The graph of the left-hand side against \(d\) is a parabola that opens upwards and has roots at \[ d_\pm=3b^2\pm\sqrt{9b^4-b^4}=b^2(3\pm\sqrt{8}). \tag{7.60}\] The condition in Eq. 7.59 tells us that \(d\) is outside the interval \((d_-,d_+)\). We have also already derived that \(d>b^2>d_-\). So the final condition on the ratio of the diffusion rates is \[ d=\frac{D_2}{D_1}>b^2(3+\sqrt{8}). \tag{7.61}\]

7.4.4 Finite domain

We now restrict space to a finite interval \([0,L]\) and impose no-flux boundary conditions at the boundaries: \[ \partial_xu(0,t)=\partial_xu(L,t)=0,~~~\partial_xv(0,t)=\partial_xv(L,t)=0. \tag{7.62}\] This implies that also the perturbations away from the homogeneous steady state must satisfy these boundary conditions: \[ \partial_x\xi(0,t)=\partial_x\xi(L,t)=0,~~~\partial_x\eta(0,t)=\partial_x\eta(L,t)=0. \tag{7.63}\] Substituting our harmonic wave Ansatz Eq. 7.44 into the boundary condition \(\partial_x\xi(0,t)=0\) gives \(k\cos(\alpha)=0\) and hence \(\alpha=\pi/2\). Substituting it into the boundary condition \(\partial_x\xi(L,t)=0\) gives \(k\cos(kL+\pi/2)=0\). This means that \(kL\) must be a multiple of \(\pi\) and hence \(k\) must be an integer multiple of \(\pi/L\). We can write this as \[ k=k_n=\frac{n\pi}{L},~~~n=1,2,... \tag{7.64}\] The boundary conditions for \(\eta\) are satisfied for the same set of wave numbers \(k_n\).

We had seen that only perturbations with wave numbers in the interval \([k_-,k_+]\) can grow, with \(k_\pm\) given by Eq. 7.50. So the condition for a Turing instability in our bounded domain is that there is at least one \(n\) for which \(k_n\) is in the interval \([k_-,k_+]\). This gives the condition \(k_1=\pi/L<k_+\) or equivalently \(L>\pi/k_+\). Turing instabilities can only occur in domains that are large enough.

7.5 Directed motion

So far we have only discussed random motion. We now want to consider directed motion. There are many reasons for individuals to move in a directed way. For example, they might move towards a food source or away from a predator. This kind of movement is referred to as taxis: the movement towards (or away from) a higher concentration of something (prey, predator, light, nutrient). If the something is a chemical this is called chemotaxis. Chemotaxis is a very important mechanism in biology. For example, immune cells move towards sites of infection, sperm cells move towards the egg, and bacteria move towards nutrients.

To model this, we introduce the density \(n(\mathbf{x},t)\) of individuals and the density \(a(\mathbf{x},t)\) of the chemical, which we will assume is an attractant, i.e., individuals move towards higher concentrations of the chemical. The individuals move in response to the gradient of the chemical, so the chemotactic movement is described by a flux \[ \mathbf{J}_c=\chi n\mathbf{\nabla} a, \tag{7.65}\] where \(\chi>0\) is the chemotactic sensitivity. Compare this with the diffusive flux \[ \mathbf{J}_d=-D\mathbf{\nabla} n. \tag{7.66}\]

The equation for the density of individuals is then \[ \partial_t n=f(n)-\mathbf{\nabla}\cdot\mathbf{J}, \tag{7.67}\] where \(f(n)\) describes the local population dynamics. This is just the higher-dimensional version of Eq. 7.5. We can now plug in the expressions for the fluxes and get \[ \partial_t n=f(n)-\mathbf{\nabla}\cdot(\chi n\mathbf{\nabla} a) +\mathbf{\nabla}\cdot(D_n\mathbf{\nabla} n). \tag{7.68}\] The attractant \(a\) is governed by the equation \[ \partial_t a=g(a,n)+\mathbf{\nabla}\cdot(D_a\mathbf{\nabla} n). \tag{7.69}\] The local dynamics of the attractant is described by \(g(a,n)\) and the diffusion of the attractant has diffusion coefficient \(D_a\).

7.5.1 Slime mould aggregation

A fascinating example of chemotaxis is the aggregation of slime moulds. Slime moulds are single-celled organisms that live in the soil and feed on bacteria. When the food supply runs out, the slime moulds aggregate to form a fruiting body that releases spores. The aggregation is guided by a chemical attractant. The slime moulds move towards higher concentrations of the attractant. The attractant is a chemical that is released by the slime moulds themselves.

When modelling the slime mould aggregation we can neglect the population dynamics which are not relevant for the aggregation process and thus set \(f(n)=0\). For the local dynamics of the attractant we use \[ g(a,n)= h n - q a, \tag{7.70}\] where \(h>0\) is the rate at which the attractant is produced by the slime moulds and \(q>0\) is the rate at which the attractant is degraded. We assume that the chemotactic sensitivity \(\chi\) and the diffusion coefficients \(D_n\) and \(D_a\) are positive constants. We also restrict ourselves again to motion in one dimension. Then the equations for the slime mould aggregation are: \[\begin{split} \partial_t n&=-\chi\partial_x(n\partial_x a)+D_n\partial_x^2 n,\\ \partial_t a&=hn-q a+D_a\partial_x^2 a. \end{split} \tag{7.71}\] The homogeneous steady state where \(n(x,t)=n^*\) and \(a(x,t)=a^*\) has to satisfies \(h n^*=q a^*\). We can now linearise the equations around the homogeneous steady state by setting \(n(x,t)=n^*+\xi(x,t)\) and \(a(x,t)=a^*+\eta(x,t)\) with \(\xi\) and \(\eta\) small. We get \[\begin{split} \partial_t \xi&=-\chi n^*\partial_x^2\eta+D_n\partial_x^2\xi,\\ \partial_t \eta&=h\xi-q\eta+D_a\partial_x^2\eta. \end{split} \tag{7.72}\] We write this linear system in matrix form as \[ \partial_t\begin{pmatrix}\xi\\\eta\end{pmatrix}=\begin{pmatrix}D_n\partial_x^2&-\chi n^*\partial_x^2\\h&-q+D_a\partial_x^2\end{pmatrix}\begin{pmatrix}\xi\\\eta\end{pmatrix}. \tag{7.73}\] To solve this we make the Ansatz \[ \begin{pmatrix}\xi\\\eta\end{pmatrix}=\begin{pmatrix}B_1\\B_2\end{pmatrix}e^{\sigma_k t}\sin(kx+\alpha). \tag{7.74}\] Plugging this into Eq. 7.73 and dividing by \(e^{\sigma_k t}\sin(kx+\alpha)\) gives \[ \sigma\begin{pmatrix}B_1\\B_2\end{pmatrix}=\begin{pmatrix}-D_nk^2&\chi n^*k^2\\h&-q-D_ak^2\end{pmatrix}\begin{pmatrix}B_1\\B_2\end{pmatrix}. \tag{7.75}\] We again denote the matrix above as \(A(k)\). The condition for the homogeneous steady state to be stable is that the real parts of the eigenvalues \(\sigma_k\) are all negative. The condition for the homogeneous steady state to be unstable is that there is at least one \(k\) for which the real part of \(\sigma_k\) is positive.

To determine the stability we again calculate the trace and determinant of the matrix in Eq. 7.75. The trace is \[ \text{tr}(A(k))=-(D_n+D_a)k^2-q. \tag{7.76}\] This is manifestly negative for all \(k\), so this does not yet indicate any instability. The determinant is \[ \det(A(k))=D_nD_ak^4+D_nqk^2-h\chi n^* k^2. \tag{7.77}\] This is negative and hence gives an instability if \[ D_nD_ak^2<h\chi n^* -D_nq \tag{7.78}\] Because the left-hand side in this inequality is positive, the right-hand side must be positive as well. This gives us the condition \[ D_n<\frac{h\chi n^*}{q}. \tag{7.79}\] This makes sense: if the individuals diffuse too fast, the aggregation can not happen. We can also rewrite the condition for the aggregation as \[ h>\frac{qD_n}{\chi n^*}. \tag{7.80}\] So aggregation only happens if the rate of production of the attractant is larger enough.

7.6 Exercises

Fishing model with diffusion

Exercise 7.1  

A population of fish \(F(x,t)\) in a river of width \(L\) with banks at \(x=0\) and \(x=L\) can be modelled by the partial differential equation \[ \frac{\partial F}{\partial t} = rF\left( 1-\frac{F}{K} \right) + D\frac{\partial^2 F}{\partial x^2}, \tag{7.81}\] where \(r\), \(K\) and \(D\) are positive constants. No-flux boundary conditions are applied at \(x=L\). At precisely \(x=0\) some shore-based fishermen catch all of the fish. We wish to find the minimum width of the river to ensure the fish population does not collapse.

  1. First determine the spatially uniform steady states and indicate their stability.

  2. Linearise the system for small \(F\).

  3. By considering a solution of the form \(F(x,t) = e^{\lambda t} \left( A \cos kx + B \sin kx \right)\) show that \[ \lambda = \lambda_n =: r-k_n^2D, \tag{7.82}\] where \[ k = k_n =: \frac{(2n+1)\pi }{2L}, ~~~ n = 0,~1,~2,... \]

  4. Hence, determine the condition on \(L\) for the fish population not to collapse.


Travelling wave in 1-species reaction-diffusion model

Exercise 7.2  

A reaction-diffusion population model has the form \[ \frac{\partial u}{\partial t}=f(u)+D\frac{\partial^2u}{\partial x^2} \tag{7.83}\] where \(D>0\) and where \(f\) satisfies \(f(0)=f(1)=0\) and \(f(u)>0\) for \(u\in(0,1)\).

  1. Convert this equation into travelling wave form by making the Ansatz that \(u(x,t)=U(z)\) and \(v(x,t)=V(z)\) with \(z=x-ct\).

  2. Assuming that a solution exists such that \(U(-\infty)=1\), \(U(\infty)=0\) show, by linearising the equation at the leading edge, that for a biologically realistic solution the wave speed \(c\) satisfies \(c\geq 2\sqrt{Df'(0)}\).

  3. Suppose that, instead, \(f(u)=0\). Show that the equation in travelling wave form becomes \(DU''+c\,U'=0\). Give the general solution of this. Does the solution look realistic for a function that represents a population?


* Travelling wave in 2-species reaction-diffusion model

Exercise 7.3  

Consider the system \[ \begin{split} \frac{\partial u}{\partial t}&=-u^2v \\ \frac{\partial v}{\partial t}&=u^2v-\rho v+\frac{\partial^2v}{\partial x^2} \end{split} \tag{7.84}\] where \(\rho\) is a positive constant. Convert this system into travelling wave form by making the Ansatz that \(u(x,t)=U(z)\) and \(v(x,t)=V(z)\) with \(z=x-ct\), and show in particular that \[ -c\frac{dV}{dz}=c\frac{dU}{dz}-\frac{\rho c}{U^2}\frac{dU}{dz}+\frac{d^2V}{dz^2}. \tag{7.85}\] Consider a solution of the travelling wave equations such that \(U(\infty)=1\), \(U(-\infty)=a\), \(V(\pm\infty)=0\) for some \(a<1\). By integrating the above equation from \(-\infty\) to \(\infty\) and imposing the boundary conditions, determine the value of \(a\). Also determine a lower bound on the wave velocity \(c\) by linearising around the leading edge of the wave.


SIR model with logistic growth

Exercise 7.4  

Assume that in the absence of rabies, the fox population is described by a logistic model with intrinsic growth rate \(r\) and carrying capacity \(K\). Adding this population dynamics to the SIR model gives the equations \[\begin{split} \partial_t S&=rS\left(1-\frac{S}{K}\right)-\beta SI\\ \partial_t I&=\beta SI-\gamma I. \end{split} \tag{7.86}\] We ignore the removed component \(R\), which in this case would correspond to dead foxes.

  1. Show that this model can be written in non-dimensionalised form as \[\begin{split} \partial_{\tilde{t}} u&=bu(1-u)-uv\\ \partial_{\tilde{t}} v&=uv-mv. \end{split} \tag{7.87}\]

  2. Determine the steady state solutions (fixed points). Under what condition on the parameters is there an endemic state?

Now assume that both susceptible and infected foxes move around randomly, but at different rates. This adds diffusion terms to the SIR model, where now \(S=S(x,t)\) and \(I=I(x,t)\) denote densities in space. \[\begin{split} \partial_t S&=rS\left(1-\frac{S}{K}\right)-\beta SI+D_1\partial_x^2S\\\ \partial_t I&=\beta SI-\gamma I+D_2\partial_x^2I. \end{split} \tag{7.88}\]

  1. Give expressions for \(\tilde{x}\) and \(d\) such that this model can be written in non-dimensionalised form as \[\begin{split} \partial_{\tilde{t}} u&=bu(1-u)-uv+\partial_{\tilde{x}}^2u\\ \partial_{\tilde{t}} v&=uv-m v+d\,\partial_{\tilde{x}}^2v \end{split} \tag{7.89}\]

We will now work with these non-dimensional equations but drop the tildes to avoid clutter.

  1. Make the travelling wave Ansatz \[ u(x,t)=A(z),~~~v(x,t)=B(z) \tag{7.90}\] with \(z=x-ct\) for some \(c>0\) and derive the system of ODEs describing the functions \(A\) and \(B\).

  2. If \(A(\infty)=1\), what are \(A(-\infty), B(\infty)\) and \(B(-\infty)\)? Make a sketch of \(A(z)\) and \(B(z)\) and indicate the direction of travel with an arrow.

  3. By linearising about the leading edge of the wave where \(B\) is very small, determine a lower limit on the wave speed \(c\).


Derive Turing instability

Exercise 7.5 Consider the reaction-diffusion model \[\begin{split} \frac{\partial u}{\partial t}&=a-u+u^2v+D_1\frac{\partial^2u}{\partial x^2} \\ \frac{\partial v}{\partial t}&=b-u^2v+D_2\frac{\partial^2v}{\partial x^2} \end{split} \tag{7.91}\] where \(b>0\) and \(a+b>0\).

Show that a spatially uniform steady state solution \((u^*,v^*)\) exists for this model and is given by \((u^*,v^*)=(a+b\,,\;b/(a+b)^2)\). Show that the conditions for this steady state to be driven unstable by diffusion are that the three inequalities \[ b-a<(a+b)^3\] \[ \left[D_2\left(\frac{b-a}{a+b}\right)-D_1(a+b)^2\right]^2>4D_1D_2(a+b)^2 \tag{7.92}\] \[ D_2\left(\frac{b-a}{a+b}\right)-D_1(a+b)^2>0 \tag{7.93}\] should all hold simultaneously.


* Conditions for Turing instability

Exercise 7.6 Consider the reaction-diffusion model \[\begin{split} \frac{\partial c_1}{\partial t}&=\delta -kc_1 - c_1c_2^2 + D_1\frac{\partial^2c_1}{\partial x^2} \\ \frac{\partial c_2}{\partial t}&=kc_1+c_1c_2^2 - c_2 + D_2\frac{\partial^2c_2}{\partial x^2} \end{split} \tag{7.94}\] where \(k>0\) and \(\delta>0\).

Show that a spatially uniform steady state solution \((c_1^*,c_2^*)\) exists for this model and give the conditions for this steady state to be driven unstable by diffusion. You may use any results derived in the lectures.


Slime mould with boundary

Exercise 7.7 In the model for the aggregation of slime mould amoebae suppose that the spatial domain is \(0\leq x\leq L\) rather than \(-\infty<x<\infty\). Show that the conditions for aggregation to occur are \(\chi a^*f>k\mu\) and \[ L>\pi\sqrt{\frac{D\mu}{\chi a^* f-k\mu}}. \tag{7.95}\]