1.2.1 Self-competition

One improvement would be to demand that we include a notion of self-competition within the population. This could, for example, model competition for food or territory. Mathematically, we need a decay term which is small for small \(x\), allowing the population to grow, but dominates the growth term when \(x\) gets larger, thus restricting its growth.

The simplest example of such a model is the logistic equation. Originally introduced by Pierre François Verhulst in 1838, the equation is nonlinear and takes the form \[\begin{equation} \label{logistic} \mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = a x\left(1-\frac{x}{K}\right), \end{equation}\] where \(a>0\) is the usual growth term and, as we shall see, \(K\) is the limiting population or carrying capacity.

Can we guess what this model looks like? See that there is an equilibrium (\(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = 0\)) at \(x=K\). The term \((1-x/K)\) is negative if \(x>K\) and positive if \(x<K\) so we might expect it to either decay towards \(K\) from above and up towards \(K\) if from below.

This sort of analysis turns out to be quite common for population models because it is not always possible to solve them analytically. Thankfully this time we can, however.

1.2.2 Solutions to the logistic equation

Let’s go ahead and now work out the solutions of our model. We can separate [logistic]: \[\int\frac{1}{x(1 - x/K)}\,\mathrm{d}x= at + C.\] We can integrate the first integral using partial fractions (remember those!), \[\frac{1}{x(1 - x/K)} = \frac{1}{x} + \frac{1}{K\left(1-x/K\right)},\] and so \[\log(x) - \log(1-x/K) =a t+C \implies \frac{x}{1-x/K} = A\mathrm{e}^{at},\] so that finally, \[\begin{equation} \label{logsol} x(t) = \frac{A \mathrm{e}^{a t}}{1+ \frac{A}{K}\mathrm{e}^{at}}. \end{equation}\]

1.2 demonstrates the behaviour of this equation.

What about that limiting population we promised? See that as \(t \to \infty\), \(x(t) \to K\): a fact independent of the initial condition.

1.2.3 Equilibria of the logistic equation

Rather than worry about how the population changes, we might only really care where it will end up, given sufficient time. We have already spotted that the system tends to a state where \(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}=0\) at \(x=K\). Is this the only possibility?

No! If we set \(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}=0\) then we have \[a x\left(1-\frac{x}{K}\right) = 0 \implies x=0, K.\] So there is also an unchanging state of \(x=0\) where there is no population. This makes sense, of course: no population means no reproduction. At this stage we make our first definition:

Equilibrium:

An equilibrium of a system is one in which all time derivatives are zero. For example, consider the system \[\mathchoice{\frac{{\mathrm d}^2 u}{{\mathrm d}t^2}}{{\mathrm d}^2 u/{\mathrm d}t^2}{{\mathrm d}^2 u/{\mathrm d}t^2}{{\mathrm d}^2 u/{\mathrm d}t^2} + \mathchoice{\frac{{\mathrm d}u}{{\mathrm d}t}}{{\mathrm d}u/{\mathrm d}t}{{\mathrm d}u/{\mathrm d}t}{{\mathrm d}u/{\mathrm d}t} = u^2+v^2,\qquad \mathchoice{\frac{{\mathrm d}^{3} v}{{\mathrm d}t^{3}}}{{\mathrm d}^{3} v/{\mathrm d}t^{3}}{{\mathrm d}^{3} v/{\mathrm d}t^{3}}{{\mathrm d}^{3} v/{\mathrm d}t^{3}} = uv.\] The equations of equilibrium are \[0 = u^2+v^2,\qquad 0 = uv\] (the only solution to which is \(u=v=0\)). The definition of equilibrium is often (in a dynamical systems context) referred to as a steady state or fixed point.

Permissible/feasible equilibrium:

We earlier demanded that the population is \(\geq 0\). We thus define a permissible or feasible equilibrium as one which satisfies this criteria. It will be important throughout the course that our models have permissible equilibria to be valid. Another idea we will come to discuss is that a good model should not allow a positive initial population to become negative.

1.2.4 A first look at stability

We note that for any positive \(A\), [logsol] will tend to \(K\). We say that the \(x=K\) equilibrium is stable because any small change from \(x=K\) (say \(x= K-\varepsilon\)) will tend back to \(x=K\) if we go forward in time (convince yourself of this by looking again at 1.2). However, if we are at \(x=0\) and there is a sudden small change to \(x=\varepsilon\), perhaps representing a small population migration, if we go forward in time it will grow inexorably towards \(K\). Thus we say that the \(x=0\) equilibrium is unstable.

But did we need to solve the logistic equation, [logistic], to find this? Actually no! Because our differential equation is of the form \(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = f(x)\), we can use a common technique where we simply plot \(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}\) against \(x\) (known as plotting the phase space) and make some observations.

Look at 1.3. Our two equilibria are marked. If you start at a value of \(x\) where \(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}\) is positive, we know that \(x\) increases in time, so you move to the right over time (indicated by the forward-pointing red arrow). And where \(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}\) is negative, you move to the left (the backward-pointing red arrow). This instinctively tells us that you will, at \(t=\infty\), always end up at \(x=K\) unless you start exactly at \(x=0\). Do you agree?

We will make a more mathematically precise notion of stability/instability of equilibria in 2 and 2.2.

1.4.1 Interactions between predator and prey

Instead of introducing self-competition, a second possibility to avoid unbounded growth is to model a second population, \(y(t)\), which represents a second species. In the case of our greenfly population, we have our antagonists in ladybirds (1.5). Since ladybirds prey on greenflies, the greenfly population, \(x\), will decrease proportionally to \[\text{[the number of ladybirds, $y$]} \times \text{[the number of greenflies, $x$],}\] i.e., the number of interactions of the two species which may lead to a sad little greenfly funeral. This law will therefore be in the form \[\begin{aligned} \mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = a x - bxy, \end{aligned}\] with \(b\) the rate at which fatal interactions occur.

But we must then also model the changing ladybird population, \(y(t)\). We assume in the absence of greenflies it will decrease as its food supply has vanished: \[\mathchoice{\frac{{\mathrm d}y}{{\mathrm d}t}}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t} = -cy.\] (Ask yourself: Why is this proportional to \(y\)?) However, it will also grow proportionally to the availability of food, i.e., interactions of the two species (at some rate \(d\)), so \[\mathchoice{\frac{{\mathrm d}y}{{\mathrm d}t}}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t} = -cy + d x y.\] So we have a coupled set of ordinary differential equations, \[\begin{align} \label{lotvol} \mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} &= ax - bxy,\\ \mathchoice{\frac{{\mathrm d}y}{{\mathrm d}t}}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t} &= -cy + dxy, \end{align}\] where, once again, \(ax\) represents prey reproduction with unlimited food, \(-cy\) is predator natural death or emigration, and \(bxy\) and \(dxy\) are interaction terms. Note that \(b\) and \(d\) are not necessarily the same: many prey may have to be eaten for the predator population to grow by one.

This system represents both the individual growth/decay of the species (self-interactions) as well as their mutual interaction. This is the so-called Lotka–Volterra (predator–prey) system developed by the Italian mathematician Vito Volterra (1860–1940) in 1926 to explain the fluctuation of fish populations in the Adriatic Sea (more on Lotka in 4). In more modern theories there are multiple species, each with their own interactions, but we will limit ourselves to this simpler but highly instructive classical system.

An example solution is shown for the parameters \((a,b,c,d) = (\frac23,\frac43,1,1)\), \(x(0)=1\), \(y(0)=1\) in 1.6(a). See how the peaks in the greenfly population naturally increase the ladybird food supply: the ladybird population then increases. In turn this leads to the greenfly population dropping as they get eaten, and this decrease in food supply leads to the ladybird population dropping as food becomes competitive.

This periodic behaviour is made clear using a phase plot, as shown in 1.6(b). In this case we have a parametrised plot \((x(t),y(t))\): a geometric plot of the variables of the system (2D here because we have two variables). Closed curves in phase space indicate a periodic relationship between the two parameters.

1.4.2 Parameter reduction of Lotka–Volterra

The parameters \((a,b,c,d)\) play a key role in determining the system’s behaviour. However, they are not all independent. We can work out how many we actually need by doing a nondimensionalisation.

Thus far we’ve dealt directly with the dimensional form of the differential equations, [lotvol], meaning that the parameters in the equations have relevant dimensions (or units) associated with them. While this is perhaps convenient for making direct predictions from measured values of the parameters, working with dimensional parameters can at best keep the equations looking ‘untidy’, and at worst obscure the mathematical structure of the equations or relevant approximations that can be made.

The process of nondimensionalisation is to separate our variables into a nondimensional bit and a dimensional bit: \[x = \widehat{x}X, \quad y = \widehat{y}Y, \quad t = \widehat{t}T.\] Here, \(x\), \(X\), \(y\) and \(Y\) have units of ‘individuals’, and \(t\) and \(T\) have units of time. Variables with hats indicate they are nondimensional.

If we use this substitution, our equations are \[\begin{aligned} \mathchoice{\frac{{\mathrm d}\widehat x}{{\mathrm d}\widehat t}}{{\mathrm d}\widehat x/{\mathrm d}\widehat t}{{\mathrm d}\widehat x/{\mathrm d}\widehat t}{{\mathrm d}\widehat x/{\mathrm d}\widehat t}\frac{X}{T} &= a\widehat{x}X - b\widehat{x}\widehat{y}XY,\\ \mathchoice{\frac{{\mathrm d}\widehat y}{{\mathrm d}\widehat t}}{{\mathrm d}\widehat y/{\mathrm d}\widehat t}{{\mathrm d}\widehat y/{\mathrm d}\widehat t}{{\mathrm d}\widehat y/{\mathrm d}\widehat t}\frac{Y}{T} &= -c\widehat{y}Y + d\widehat{x}\widehat{y}XY. \end{aligned}\] Rearranging so that there are no dimensions on the left hand sides, we have \[\begin{aligned} \mathchoice{\frac{{\mathrm d}\widehat x}{{\mathrm d}\widehat t}}{{\mathrm d}\widehat x/{\mathrm d}\widehat t}{{\mathrm d}\widehat x/{\mathrm d}\widehat t}{{\mathrm d}\widehat x/{\mathrm d}\widehat t} &= a\widehat{x}T - b\widehat{x}\widehat{y}YT,\\ \mathchoice{\frac{{\mathrm d}\widehat y}{{\mathrm d}\widehat t}}{{\mathrm d}\widehat y/{\mathrm d}\widehat t}{{\mathrm d}\widehat y/{\mathrm d}\widehat t}{{\mathrm d}\widehat y/{\mathrm d}\widehat t} &= -c\widehat{y}T + d\widehat{x}\widehat{y}XT. \end{aligned}\] So now we pick \(T\), \(X\) and \(Y\) to remove as many of our parameters as possible. If we choose \[\begin{equation} \label{TYX-scaling} T = \frac{1}{a}, \quad Y = \frac{a}{b}, \quad X = \frac{c}{d}, \end{equation}\] then the system can be written as \[\begin{align} \label{reducedlotvol} \mathchoice{\frac{{\mathrm d}\widehat{x}}{{\mathrm d}\widehat{t}}}{{\mathrm d}\widehat{x}/{\mathrm d}\widehat{t}}{{\mathrm d}\widehat{x}/{\mathrm d}\widehat{t}}{{\mathrm d}\widehat{x}/{\mathrm d}\widehat{t}} &= \widehat{x} - \widehat{x}\widehat{y},\\ \mathchoice{\frac{{\mathrm d}\widehat{y}}{{\mathrm d}\widehat{t}}}{{\mathrm d}\widehat{y}/{\mathrm d}\widehat{t}}{{\mathrm d}\widehat{y}/{\mathrm d}\widehat{t}}{{\mathrm d}\widehat{y}/{\mathrm d}\widehat{t}} &= \gamma(-\widehat{y} +\widehat{x}\widehat{y}), \end{align}\] where \[\gamma = c/a.\] We see that nondimensionalisation has ‘tidied up’ our equation – we went from four parameters to one. Why might this matter? Well, if you fix a set of parameters you get a solution. Select another set of parameters and you get another solution. So (for given initial conditions) there are as many solutions as there are parameters. If we have decided that the four parameters are positive real numbers, then there is a four-dimensional space of solutions: quite a lot if we want to map out all the system’s behaviour! In fact, we have shown that the parameters relate to each other and that there is, in fact, only a one-dimensional space of solutions, a much simpler search.

The scalings show that the solutions just relate by constant stretching factors. For example, if I choose \(\gamma=1\), I can then set \(c=2\), thus \(a=2\). In fact, the parameters \(b\) and \(d\) are even redundant and we can just choose them (but this is rare, Lotka–Volterra is a slightly odd system). The ratio \(c/d\) just stretches any \(x\) solution (stretches its range/amplitude); \(a/b\) stretches \(y\); and the scaling along \(t\) changes the period of the solutions. That is to say, for a given \(\gamma\) we have a main solution which can be simply scaled to get other solutions, without solving the system again.

1.4.3 Solutions to Lotka–Volterra

We can solve this system using separation of variables. Dividing the two equations (and dropping hats) we obtain \[\begin{equation} \label{sepvareq} \mathchoice{\frac{{\mathrm d}y}{{\mathrm d}x}}{{\mathrm d}y/{\mathrm d}x}{{\mathrm d}y/{\mathrm d}x}{{\mathrm d}y/{\mathrm d}x}= \frac{\gamma y}{x}\left(\frac{-1+x}{1- y}\right)\implies \int \frac{1-y}{y} \,\mathrm{d}y= -\gamma \int \frac{1-x}{x} \,\mathrm{d}x. \end{equation}\] Integrating both sides of [sepvareq] we obtain \[\begin{equation} \label{dynsol} \log y-y =- \gamma(\log x -x) + C, \end{equation}\] where the constant \(C\) can be set by some initial condition, \((x(0),y(0))\). Unfortunately it is not possible to write this relationship in explicit form. This gives us the phase curves determined by the value of the constant \(C\). Parametrising this curve then gives the solutions \(x(t)\) and \(y(t)\), i.e. we could choose some behaviour for \(x(t)\) and [dynsol] will then determine the behaviour of \(y\). We will find that this kind of solution is common to such systems.

1.4.4 Equilibria of Lotka–Volterra

Looking at [reducedlotvol] we can see there are two possible equilibria where \(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = \mathchoice{\frac{{\mathrm d}y}{{\mathrm d}t}}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t}=0\) \(\forall t\): \[\begin{equation} \label{equlibria} (x,y) = (0,0) \quad \mbox{and} \quad (x,y) = (1,1). \end{equation}\] (Dynamical systems fans would call these fixed points or steady state solutions). The \((0,0)\) solution corresponds to both populations being extinct! The second corresponds to the nonzero population densities at which the population sizes will remain fixed.

In 1.7(a), we see the varying behaviour of the closed phase curves of the system. All curves encircle the equilibrium at \((1,1)\) and as the initial conditions get closer to the equilibrium value, the radius of the curve decreases. In 1.7(b), we see the dramatic variety of morphology the parametric curves can exhibit. When the pair \((x(0),y(0))\) are initially close to the equilibrium, the curves have low-amplitude sinusoidal shape, while if \(y(0)\) is initially small, the curves have extremely sharp gradients and dramatic rates of change at the maxima.

1.4.5 Stability of Lotka–Volterra?

We have our dynamic solutions, [dynsol], and the fixed point equilibria, [equlibria]. A number of questions beg to be asked at this point:

1. Can one or both of the species die out if they are both nonzero at some time \(t\)?

2. Can an oscillating pair of populations relax to nonzero fixed values, i.e., do the populations ever settle?

An immediate observation in regards to (i) is that [dynsol] only allows \(x=0\) when \(y=0\) and vice versa, so they would have to become extinct simultaneously. The existence of periodic solutions as shown in 1.7 seems to suggest neither (i) or (ii) can occur, because the system repeats itself cyclically. A solution which decays into equilibrium would have to have a phase space diagram which spiralled inwards. In fact, we have a precise means of determining the answer to such questions which we will discuss in 2.

2.1.1 Equations of motion

Let’s consider a rigid pendulum: a bead of mass \(m\) attached to a rigid rod of length \(\ell\), depicted in 2.1(a). Three forces act on the bead: its weight, \(mg\); tension in the rod, \(T\); and friction, which opposes the motion of the bead.

The distance travelled by the bead is the arclength, \(s=\ell \theta\). The velocity and acceleration are therefore \[\mathchoice{\frac{{\mathrm d}s}{{\mathrm d}t}}{{\mathrm d}s/{\mathrm d}t}{{\mathrm d}s/{\mathrm d}t}{{\mathrm d}s/{\mathrm d}t} = \ell\mathchoice{\frac{{\mathrm d}\theta}{{\mathrm d}t}}{{\mathrm d}\theta/{\mathrm d}t}{{\mathrm d}\theta/{\mathrm d}t}{{\mathrm d}\theta/{\mathrm d}t}, \quad \mathchoice{\frac{{\mathrm d}^2 s}{{\mathrm d}t^2}}{{\mathrm d}^2 s/{\mathrm d}t^2}{{\mathrm d}^2 s/{\mathrm d}t^2}{{\mathrm d}^2 s/{\mathrm d}t^2} = \ell\mathchoice{\frac{{\mathrm d}^2 \theta}{{\mathrm d}t^2}}{{\mathrm d}^2 \theta/{\mathrm d}t^2}{{\mathrm d}^2 \theta/{\mathrm d}t^2}{{\mathrm d}^2 \theta/{\mathrm d}t^2}.\] We say friction is proportional to velocity, \(\nu\ell \, \mathchoice{\frac{{\mathrm d}\theta}{{\mathrm d}t}}{{\mathrm d}\theta/{\mathrm d}t}{{\mathrm d}\theta/{\mathrm d}t}{{\mathrm d}\theta/{\mathrm d}t}\), and points in the direction opposing motion.

If we resolve parallel and perpendicular to the rod, given the rod makes an angle \(\theta\) to the vertical, then: \[\begin{aligned} \text{[parallel]} & & 0 &= T - mg\cos\theta, \\ \text{[perpendicular]} & & m \ell\mathchoice{\frac{{\mathrm d}^2 \theta}{{\mathrm d}t^2}}{{\mathrm d}^2 \theta/{\mathrm d}t^2}{{\mathrm d}^2 \theta/{\mathrm d}t^2}{{\mathrm d}^2 \theta/{\mathrm d}t^2} &= - \nu\ell \mathchoice{\frac{{\mathrm d}\theta}{{\mathrm d}t}}{{\mathrm d}\theta/{\mathrm d}t}{{\mathrm d}\theta/{\mathrm d}t}{{\mathrm d}\theta/{\mathrm d}t} - mg\sin\theta, \end{aligned}\] the second equation of which can be rearranged to give \[\begin{equation} \label{peneq} \mathchoice{\frac{{\mathrm d}^2 \theta}{{\mathrm d}t^2}}{{\mathrm d}^2 \theta/{\mathrm d}t^2}{{\mathrm d}^2 \theta/{\mathrm d}t^2}{{\mathrm d}^2 \theta/{\mathrm d}t^2} + \frac{\nu}{m} \mathchoice{\frac{{\mathrm d}\theta}{{\mathrm d}t}}{{\mathrm d}\theta/{\mathrm d}t}{{\mathrm d}\theta/{\mathrm d}t}{{\mathrm d}\theta/{\mathrm d}t}+ \frac{g}{\ell}\sin\theta = 0. \end{equation}\] This is a nonlinear ODE; our physical intuition tells us that the pendulum will swing with a decreasing amplitude until it relaxes to \(\theta=0\).

2.1.2 Solutions to the pendulum equation

The damped pendulum system is not an integrable system, i.e., there are no general closed-form solutions. This is generally the case for mechanical systems with friction. On the other hand, numerical solutions are simple to obtain.

An example is shown in 2.2(a): the solution is indeed oscillatory with decreasing amplitude. In (b), we see a pendulum with a very high coefficient of friction (a pendulum in fluid, say) where the damping is so strong that the angle never becomes negative: its motion is killed off on the first swing.

2.1.3 Equilibria of the pendulum equation

Accounting for periodicity, there are two equilibria to [peneq]: \[\frac{g}{\ell}\sin\theta = 0\] gives \(\theta(t) = 0\) and \(\theta(t) = \pi\) for all \(t\). The first \(\theta=0\) solution corresponds to a pendulum starting at the bottom of its cycle and not moving.

The second solution is far more interesting. This is when the pendulum rod is vertically upwards as in 2.1(b). The forces in the system are balanced as the rod tension balances gravity and the rotating moments are equal in either direction. But have you tried this? No matter how hard you try, it will eventually fail and the pendulum will start to rotate back to its \(\theta=0\) equilibrium. Why?

Any small variation in the pendulum – a gentle breeze or vibration in the rod – no matter how small, always grows over time. In practice, no system is perfect and such variations always exist. Mathematically we represent small variations by a linear stability analysis.

2.1.4 Order notation, \(\mathcal{O}\)

In asymptotic analysis, we say a function \(g(t,\varepsilon)\) is \(\mathcal{O}(f(t,\varepsilon))\) if \[\lim_{\varepsilon\to 0} \frac{g(t,\varepsilon)}{f(t,\varepsilon)} = C\] with \(C\) bounded.

So, for example, let’s have \(g(t,\varepsilon) = C \varepsilon\) and \(f(t,\varepsilon) = \varepsilon\). Well, \[\lim_{\varepsilon\to 0} \frac{g(t,\varepsilon)}{f(t,\varepsilon)} = C,\] so \(C \varepsilon\) is \(\mathcal{O}(\varepsilon)\).

But if \(g(t,\varepsilon) = C \varepsilon^2\) and \(f(t,\varepsilon) = \varepsilon\), \[\lim_{\varepsilon\to 0} \frac{g(t,\varepsilon)}{f(t,\varepsilon)} = 0.\] So \(C \varepsilon^2\) is also \(\mathcal{O}(\varepsilon)\).

On the other hand, if \(g(t,\varepsilon) = C\) and \(f(t,\varepsilon) = \varepsilon\), then \[\lim_{\varepsilon\to 0} \frac{g(t,\varepsilon)}{f(t,\varepsilon)} = \infty.\] So \(C\) is in some way much bigger than \(\varepsilon\) when \(\varepsilon\to 0\). We would say \(C\) is \(\mathcal{O}(1)\) but not \(\mathcal{O}(\varepsilon)\).

For non-polynomial functions the way in is to use Taylor expansion: what if \(g(t,\varepsilon) = \sin\varepsilon\) and \(f(t,\varepsilon) = \varepsilon\)? \[\lim_{\varepsilon\to 0} \frac{g(t,\varepsilon)}{f(t,\varepsilon)} = \lim_{\varepsilon\to 0}\frac{\sin\varepsilon}{\varepsilon} = \lim_{\varepsilon\to 0}\frac{\varepsilon-\varepsilon^3/3!+\cdots}{\varepsilon} = \lim_{\varepsilon\to 0}\left( 1 - \frac{\varepsilon^2}{3!} + \cdots \right) = 1.\] So \(\sin\varepsilon\) is also \(\mathcal{O}(\varepsilon)\).

Long story short: if something is \(\mathcal{O}(\varepsilon^n)\), it is the same size as, or smaller than, \(\varepsilon^n\) when \(\varepsilon\to 0\).

2.1.5 Four steps of a linear stability analysis

The basic steps of a linear stability analysis are as follows

1. Find the system’s equilibria, \(\theta_0\) (we will commonly use the subscript \(0\) to indicate the equilibrium solution).

2. Assume a value which is changed from this equilibrium value by a very small amount, \(\theta = \theta_0 +\varepsilon\theta_1\), with \(\varepsilon\ll 1\). This is supposed to mimic the small vibration in the system. Substitute this into the equation and ignore terms of order \(\varepsilon^2\) and higher. We are left with the behaviour of the system where only small vibrations matter.

3. Solve this system to find out if our small vibration, \(\theta_1\), grows (like it would for \(\theta_0=\pi\) in the pendulum) or decays (as it would for \(\theta_0=0\)).

4. Conclude that the equilibrium is stable if we have decay (small vibrations would disappear) and unstable if they grow.

Let’s apply this to our pendulum.

2.1.6 Linear stability analysis of the damped pendulum

2.1.7 Linear stability analysis of the frictionless pendulum

If there is no friction, \(\nu =0\), steps – will be the same, except that our linearised solution is \[\theta_1(t) = A\mathrm{e}^{\lambda_1 t} + B \mathrm{e}^{\lambda_2 t},\] where \[\lambda_{1} = \sqrt{\mp \frac{g}{\ell}}, \quad \lambda_{2} = - \sqrt{\mp \frac{g}{\ell}}.\]

2.1.8 Summary

1. To test the feasibility (stability) of equilibrium solutions we linearise the equation about the equilibrium state.

2. We analyse the solutions for exponential decay or growth. If there is only decay, then the solution is stable in that it is resistant to small changes. If there is any growth it is unstable as small changes destroy the solution.

3. Some special systems will have neither growth nor decay. In this case we know, from the nonlinear behaviour of the frictionless \(\theta_0=0\) pendulum equation, that this is because neighbouring solutions are periodic.

We can make these conclusions far more general…

3.1.1 Phase portraits

It’s important to keep in mind that the linear analysis about equilibria only gives a local description of the dynamics. To get a better picture of the global dynamics, we’d like to understand how the equilibria are connected, or if closed orbits exist. In one-dimensional systems, \(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = f(x)\), this was easy as we only needed to examine the sign of \(f\) between the equilibria and the phase portrait was clear.

In two dimensions, the phase portrait will typically include:

• The equilibria and some indication of their stability,

• The trajectories near equilibria and periodic orbits that show the flow pattern.

While uniqueness of the solution dictates that trajectories cannot cross, compiling the phase portrait in 2D can be challenging. We’re going to construct the example phase portrait in 3.2. A good recipe looks like this:

1. Draw on the equilibria, coloured in according to their stability from a linear stability analysis

2. Draw the nullclines. These are the curves where either \[\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = f(x,y) = 0 \quad\text{or}\quad \mathchoice{\frac{{\mathrm d}y}{{\mathrm d}t}}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t} = g(x,y) = 0.\] Since the equilibria satisfy both \(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}=0\) and \(\mathchoice{\frac{{\mathrm d}y}{{\mathrm d}t}}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t} = 0\), the nullclines will intersect at the equilibria.

   Particularly helpful is that when trajectories intersect the nullcline, they must be either vertical (\(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = 0\)) or horizontal (\(\mathchoice{\frac{{\mathrm d}y}{{\mathrm d}t}}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t}{{\mathrm d}y/{\mathrm d}t} = 0\)).

3. The nullclines provide the boundary between regions of phase space where \(f>0\) or \(f<0\), and \(g>0\) or \(g<0\). Thus, the nullclines give us some sense of the direction of the trajectories in different regions of the phase plane. However, there is no way to see which sign \(f,g\) have in each region a priori, so we have to test at least one point in each region.

4. Infer the directions from this information. For example, in a region where \(f<0\) and \(g<0\), we’d expect trajectories heading down and to the left.

See how we’ve done this in 3.2 for a made-up scenario. At this point we can try to plot the full trajectories (although I haven’t done so here).

With this in mind, let’s look at Lotka–Volterra and its variants in more detail.

3.3.1 Reflections on our new model

We saw in the original Lotka–Volterra model, [reducedlotvol], that it was impossible for any population to become extinct. In the competitive Lotka–Volterra model, we introduced the notion of independent growth limitation on the individual populations to our system through the logistic model.

This has reintroduced extinction: we now have a model which allows one species to ‘win’, in addition to cases where both populations settle on fixed values. But in doing so we have found limitations. The main one is that we appear to have lost the periodic/cyclic nature of the original Lotka–Volterra model. (This is in fact true on account of the so-called Dulac criterion, but we won’t prove it in this course.)

4.1.1 Logistic growth with harvesting

We saw in 1.2 that the logistic equation has two equilibria, one which is stable and one which is unstable. Let’s now make a small modification to the logistic equation which accounts for constant-rate harvesting, \[\begin{equation} \label{logistic-harvest} \mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = ax\left(1-\frac{x}{K}\right) - H. \end{equation}\] Let’s look at the phase space in 4.1(a).

As \(H\) increases, we can see graphically we go from two equilibria, to one, to none! This corresponds to the roots of \(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = 0\) in [logistic-harvest] becoming equal and then imaginary: the quadratic equation tells us \[\begin{equation} \label{bifurcation-equilibrium-harvesting} x_0 = \frac{K \pm \sqrt{K^2 - 4HK/a}}{2}. \end{equation}\] So if we increase \(H\), as soon as we go over \(H=aK/4\), \(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}\) becomes always negative and we reach extinction in finite time. An interesting result about the dangerous of overharvesting! This discontinuous jump in the long-time solution is known as a catastrophe.

In general, when changes in parameter values in our models lead to qualitative changes in the equilibria and their stability, we say a bifurcation has occurred. The parameter values at which we see the bifurcation is called the bifurcation point.

To emphasise this point: bifurcation points are about parameters. In this harvesting model, we could ask this question more directly: ‘what are the equilibria, \(x_0\), for any given parameter, \(H\)?’ And we’ve actually answered this already, in [bifurcation-equilibrium-harvesting].

So why don’t we plot \(x_0\) as a function of \(H\)? That’s what we’ve done in 4.1(b): you’ll notice that we’ve plotted both of the ‘\(\pm\)’ branches of the solution. Plotting equilibria against parameters is what’s known as a bifurcation diagram and when people draw these, they note which branches refer to stable equilibria, and which refer to unstable equilibria. The stable equilibrium branch gets a solid line and the unstable equilibrium branch gets a dotted line. They can be a bit hard to get your head around to start with, so I will always draw them in purple to remind you that it’s not a phase space.

Can you spot where the bifurcation point is on the diagram? It’s where the unstable and stable lines meet. We’ll see a few more of these diagrams shortly.

4.1.2 Insect outbreak: the spruce budworm

One cannot take a mathematical biology course without encountering the spruce budworm. It shows how a simple, one-dimensional systems can yield rich dynamics with relevant predictions of a biological phenomenon. The spruce budworm, shown in 4.2, is a moth that infects spruce trees (a type of fir tree) in eastern Canada and the US. The moths produce larvae which feed on the needles of the conifers: a serious infestation can lead to complete defoliation of a forest in about four years, and so the budworm is considered one of the most destructive pests in North America.

The budworm is preyed upon by birds that for low budworm population feed only upon the budworms once the budworm population has reached a certain level. Assuming logistic growth for the budworm population in the absence of birds, the budworm population size, \(x\), is governed by the variation on the logistic equation, \[\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = rx\left(1-\frac{x}{K}\right) - \frac{Bx^2}{A^2 + x^2}.\] What does the new term tell us? As \(x \to \infty\), the predation rate becomes \(B > 0\), similar to our harvesting model. Meanwhile, \(A > 0\) provides a measure of the threshold population size where predation suddenly increases.

Nondimensionalising with \[x = A \widehat{x}, \quad t = \frac{A}{B}\widehat{t},\] and dropping hats, we get \[\begin{equation} \label{budworm-nondim} \mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = Rx\left(1-\frac{x}{k}\right) - \frac{x^2}{1+x^2}, \end{equation}\] where \(R=rA/B\) and \(k=K/A\). Once again nondimensionalisation has reduced the parameter space, from four to two!

Find equilibria

The equilibria are given by \(x=0\) and \[\begin{equation} R\left(1-\frac{x}{k}\right) = \frac{x}{1+x^2}, \label{spruce-rk} \end{equation}\] where the logistic growth is balanced by the predation. Thus, the equilibria occur when the line given by the left-hand side of [spruce-rk] intersects the curve given by the right-hand side. Naturally, we’d like to explore how these equilibria vary with \(R\) and \(k\) and, in doing so, we see the benefit of our choice of nondimensionalisation.

The curve given by \(g(x) = x/(1 + x^2)\) does not depend on either parameter, and as for the line \(h(x) = R(1-x/k)\), \(R\) is the \(y\)-intercept, while \(k\) is the where it crosses the \(x\)-axis. This is plotted in 4.3, and as you can see, this makes graphical evaluation quite easy!

Just like in the harvesting model (where we changed the parameter \(H\)), we see that changing \(R\) and \(k\) can lead to us going from three additional equilibria (we already have \(x=0\)), to two, to one: more bifurcations.

4.1.3 General bifurcations of single-species models

In general, a single-species model (that is, a single ODE), can only do so many different things as parameters are varied. In particular, the only long-time behaviour that can occur for a finite population is that it tends to a steady state value.

But what sorts of different things can happen as parameters vary? Well it turns out that in general you can have a number of different scenarios, depending on how many parameters vary, and how ‘generic’ or commonplace we expect these behaviours to be. Rather than formally and systematically listing all of these bifurcations, let’s just consider the two most common examples, and one example of a ‘complex’ bifurcation diagram formed by combining these basic bifurcations together. Importantly, these are not just examples for single-species models, but they can and do occur commonly in larger models of many species.

Saddle node bifurcation

A saddle node bifurcation (sometimes called a fold bifurcation or annihilation point) occurs when an unstable equilibrium and a stable equilibrium ‘collide’ at the same value, and then both cease to exist afterwards. This is what happened in 4.1.1 as \(H\) increased beyond \(aK/4\), our bifurcation point. Confirm you agree with this by looking at 4.1.

The model there was a little awkward to draw the bifurcation diagram for but actually the prototypical model of the saddle node bifurcation is the much simpler equation \[\begin{equation} \label{prototypical-saddle-node} \mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = p+x^2, \end{equation}\] where \(p\) is a parameter.

Draw the phase space of [prototypical-saddle-node] and check how many equilibria it has, and what their stabilities are, as \(p\) varies. Don’t worry about negative \(x\) not being feasible for representing biological populations: here we are looking at the broader mathematical principle.

By plotting the equilibria against the parameter, we create a bifurcation diagram, and you can see it in 4.5(a).

Hysteresis (again)

We can combine these two simple models to find more complex behaviour, the type of which we have already seen in the spruce budworm model.

Let’s take the pitchfork prototype and add the shift from the saddle node prototype: in other words, let’s consider the two-parameter model \[\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = p + x(q-x^2),\] for two parameters \(p\) and \(q\). For this discussion, we’ll fix \(q=1\) but we could do the same analysis as what follows by fixing \(p\). So we have \[\begin{equation} \label{prototypical-hysteresis} \mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = p + x(1-x^2). \end{equation}\] What does the phase space look like and how many equilibria do we have for any \(p\)?

The phase space is drawn in 4.7 and we can see that equilibria appear and disappear as we change \(p\). We can find the equilibria, \(x_0\), in terms of \(p\) by once again solving this equation when \(\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}=0\). Plotting these \(x_0\) when they are real as a function of \(p\) gives us the bifurcation diagram in 4.5(c).

Just like in the spruce budworm, we see we can vary \(p\) such that returning \(p\) to its initial value doesn’t return the system to its original state: hysteresis. If you don’t see why, read again the discussion under the spruce budworm on page and see that it applies exactly to this model.

These three bifurcation phenomena – saddle nodes, pitchforks and hysteresis – will reappear next term: for now, we’re learning about them because with these bifurcation diagrams you can read off essentially the entire system’s behaviour for any given parameter. Pretty cool!

5.1.1 Discrete random walks

As a concrete example, let us consider an individual on a line initially at a position \(x=0\). The individual moves left or right in integer steps with a probability of moving right \(p\) and hence left of \((1-p)\). After \(n\) steps a path can be encoded as LRLLRRR…. If we repeat this process a large number of times we get a set of paths – a set of random walks. In 5.1, we show five random walks of length \(100\) for \(p=0.5\), \(0.6\) and \(0.9\) respectively. Here you see the individual paths; in 5.2 we see the same for 200 walks! For \(p=0.5\), the end positions are reasonably evenly spread about \(x=0\). For \(p=0.6\), the story is similar except that the spread is about the mean positions, \(x=2np-n\), and the spread is not equal either side of this mean. This nonzero mean, indicated by the dotted line, can be interpreted as a natural drift of the set of individuals all starting at the same point due to the probability bias.

5.1.2 Continuum limit of random walks

In our course, we will focus on the continuum limit of such a model. Rather than making the steps of size \(1\), we specify them to be \({\mathrm d}x\) which will be vanishingly small, and each step is taken in a time \({\mathrm d}t\). We let \(c(x,t)\) be the continuous probability that, at a time \(t\), a particle (population member) reaches a displacement \(x\) at a time \(t\). This implies that at a time \(t-{\mathrm d}t\), the particle must have been at either \(x-{\mathrm d}x\) or \(x+{\mathrm d}x\). If \(p\) is the probability that the particle moves right, we therefore have \[c(x,t)= p c(x-{\mathrm d}x,t-{\mathrm d}t) + (1-p) c(x+{\mathrm d}x,t-{\mathrm d}t).\] If we Taylor expand this to \(\mathcal{O}(\varepsilon^2)\), where \(\varepsilon= \max({\mathrm d}t,{\mathrm d}x)\), we obtain \[\begin{aligned} c(x,t) &\approx p\left[c(x,t) - {\mathrm d}x \mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} - {\mathrm d}t \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} + \frac{1}{2}\left( {\mathrm d}x^2 \mathchoice{\frac{\partial^2 c}{\partial x^2}}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2} + 2 \, {\mathrm d}x \, {\mathrm d}t \mathchoice{\frac{\partial^2 c}{\partial x \partial t}}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t} + {\mathrm d}t^2 \mathchoice{\frac{\partial^2 c}{\partial t^2}}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2} \right)\right] \\ & \quad + (1-p)\left[c(x,t) + {\mathrm d}x \mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} - {\mathrm d}t \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} + \frac{1}{2}\left( {\mathrm d}x^2 \mathchoice{\frac{\partial^2 c}{\partial x^2}}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2} - 2 \, {\mathrm d}x \, {\mathrm d}t \mathchoice{\frac{\partial^2 c}{\partial x \partial t}}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t} + {\mathrm d}t^2 \mathchoice{\frac{\partial^2 c}{\partial t^2}}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2} \right)\right] \nonumber\\ &= c(x,t) + (1-2p){\mathrm d}x \mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} - {\mathrm d}t \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} + \frac{1}{2}{\mathrm d}x^2 \mathchoice{\frac{\partial^2 c}{\partial x^2}}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2} - (1-2p)\, {\mathrm d}x \, {\mathrm d}t \mathchoice{\frac{\partial^2 c}{\partial x \partial t}}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t} + \frac12 {\mathrm d}t^2 \mathchoice{\frac{\partial^2 c}{\partial t^2}}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2}, \end{aligned}\] so that to order \(\mathcal{O}(\varepsilon^2)\), \[\mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} = (1-2 p)\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}\mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} + \frac12 \frac{({\mathrm d}x)^2}{{\mathrm d}{t}} \mathchoice{\frac{\partial^2 c}{\partial x^2}}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2} - (1-2p)\,\mathrm{d}x\mathchoice{\frac{\partial^2 c}{\partial x \partial t}}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t}{\partial^2 c /\partial x \partial t} + \frac12 {\mathrm d}t \mathchoice{\frac{\partial^2 c}{\partial t^2}}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2}{\partial^2 c/\partial t^2}.\] Taking the limit \({\mathrm d}x\to 0\) and \({\mathrm d}t\to 0\), the final two terms will vanish. The coefficient \((1-2p)\,\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}\) has the dimensions of a velocity, so we set \((2p-1)\,\mathchoice{\frac{{\mathrm d}x}{{\mathrm d}t}}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t}{{\mathrm d}x/{\mathrm d}t} = v\) (with the flipped sign so that it points in the positive direction for high \(p\) biases) and we define a constant \(D\), \[D = \frac12 \frac{({\mathrm d}x)^2}{{\mathrm d}t},\] which we call the diffusion constant for reasons which will follow. One should note that the ratio \(v/D\) will diverge as \({\mathrm d}x \to 0\) unless \(p\) is very close to \(1/2\), so as the step size gets smaller, the velocity term (the so-called advective term) dominates. Then we have \[\begin{equation} \label{advectiondiffusion} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} = -v\mathchoice{\frac{\partial c}{\partial x}}{\partial c/\partial x}{\partial c/\partial x}{\partial c/\partial x} + D \mathchoice{\frac{\partial^2 c}{\partial x^2}}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}, \end{equation}\] a partial differential equation determining the behaviour of \(c(x,t)\). Note that if \(p=1/2\), there is no velocity: we shall come to associate this term with a drift. The physical description of the term, \(D\), can be seen by deriving [advectiondiffusion] in a continuum setting, as we’ll do now.

5.4.1 The Fourier transform

The Fourier transform of a function \(f(x)\), denoted \(\mathcal{F}(k)[f(x)]\), is given by \[\mathcal{F}(k)[f(x)] = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\mathrm{e}^{-\mathrm{i}k x}f(x)\,\mathrm{d}x.\] Its inverse, \(\mathcal{F}^{-1}(x)[g(k)]\), is \[\mathcal{F}^{-1}(x)[g(k)] = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\mathrm{e}^{\mathrm{i}k x}g(k)\,{\mathrm d}{k}.\] You might have seen a slightly different version of these where only one of the terms has a \(1/2\pi\) factor, but either way, crucially it can be shown that \[\mathcal{F}^{-1}[\mathcal{F}[f(x))]] = f(x).\] You’re not expected to derive this result, merely to be aware of it. The critical property of the Fourier transform which makes it of use in solving PDEs is its effect on derivatives. The following result can be fairly easily demonstrated: \[\begin{equation} \label{fderiv} \mathcal{F}(k)\left[\mathchoice{\frac{\partial^{n} f}{\partial x^{n}}}{\partial^{n} f/\partial x^{n}}{\partial^{n} f/\partial x^{n}}{\partial^{n} f/\partial x^{n}}\right] = (\mathrm{i}k)^{n}\mathcal{F}(k)[f(x)]. \end{equation}\] Again, you don’t need to derive this result, you just need to know how to apply it.

5.4.2 Solving the diffusion equation

Applying [fderiv] to the first term in [diffprob] we obtain \[\mathchoice{\frac{\partial\mathcal{F}(k)[u]}{\partial t}}{\partial\mathcal{F}(k)[u]/\partial t}{\partial\mathcal{F}(k)[u]/\partial t}{\partial\mathcal{F}(k)[u]/\partial t} = -Dk^2\mathcal{F}(k)[u].\] Just like that, we have transformed a second order PDE into a first order ODE. Its solution is clearly an exponential, \[\mathcal{F}(k)[u(x,t)] = \mathcal{F}(k)[u(x,0)]\mathrm{e}^{-D k^2 t}.\]

To get the real space solution, we must apply the inverse operator. This is the price we pay for using the transform to simplify the system. This is not always straightforward, but in this case this is not so problematic. First we must find the initial condition \(\mathcal{F}(k)[c(x,0)]\), which, from [deltafunc], is the Fourier transform of the Dirac delta function, \[\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\delta(x)\mathrm{e}^{-\mathrm{i}kx}\,\mathrm{d}x= \frac{1}{\sqrt{2\pi}}.\] Thus our final solution is \[u(x,t) = \mathcal{F}^{-1}[\mathcal{F}(k)[u(x,t)]] = \frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{e}^{\mathrm{i}kx}\mathrm{e}^{-D k^2 t}\,{\mathrm d}{k}.\] In order to perform this integral we note the following known result: \[\mathcal{F}(k)[\mathrm{e}^{-\alpha x^2}] = \frac{1}{\sqrt{2\alpha}}\mathrm{e}^{-k^2/4\alpha}.\] We therefore have \[\begin{equation} \label{fundamentalsol} u(x,t) = \frac{1}{2 \sqrt{\pi Dt}} \mathrm{e}^{-x^2/4Dt}. \end{equation}\]

This tells us the density spreads out over time following a Gaussian distribution, as depicted in 5.6.

As we remarked at the start of the chapter, the limit of the random walk as the step size reduces to zero is a Gaussian distribution. We then showed, using a slightly hand-wavy argument, that a limit of this distribution was the advection–diffusion equation. We have now squared the circle so to speak, in the sense that we see that the diffusion equation gives Gaussian-like spreading of a point source. In fact it was Einstein (and Marian Smoluchowski) who first made this link precise at the turn of the previous century.

This solution, [fundamentalsol], is known as the fundamental solution of the diffusion equation. Why is it so fundamental? Because it allows us to answer the question ‘what about different initial conditions?’.

5.4.3 Fundamental solutions (or Green’s functions)

Why is a delta function which spreads out over time called the fundamental solution? In general we might have some linear differential operator \(\mathcal{L}\). When acting on a function \(u(\mathbfit{x},t)\) we want to solve the following initial value problem on a domain \(V\) with boundary \(\partial V\): \[\begin{equation} \label{fullprob} \mathcal{L}u(\mathbfit{x},t) = 0,\quad u(\mathbfit{x},0) = g(\mathbfit{x}),\quad u(\partial{V},t) = 0, \end{equation}\] that is to say, \(g(\mathbfit{x})\) is our initial condition.

For example, the diffusion equation operator is \(\mathcal{L}= \mathchoice{\frac{\partial}{\partial t}}{\partial/\partial t}{\partial/\partial t}{\partial/\partial t}- D\nabla^2\), so that \[\mathcal{L}u(\mathbfit{x},t) = \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} -D\nabla^2 u = 0.\] \(\mathcal{L}\) is linear so for \(a,b\) scalar constants and \(u,v\) scalar densities, \[\mathcal{L}(au+bv) = a\mathcal{L}u + b\mathcal{L}v.\] The fundamental solution or Green’s function of \(\mathcal{L}\) is defined as the solution, \(u_f\), to \[\begin{equation} \label{fundprob} \mathcal{L}u_f(\mathbfit{x},t) = 0, \quad u_f(\mathbfit{x},0) = \delta(\mathbfit{x}), \quad u_f(\partial{V},t) = 0, \end{equation}\] where specifically the initial condition is the delta function. See how our setup for the diffusion equation, [diffprob], is in this form.

We can relate the two initial conditions by the delta function identity, \[\begin{equation} \label{deltaid} \int_{V}\delta(\mathbfit{x}- \mathbfit{s})g(\mathbfit{s})\,{\mathrm d}\mathbfit{s} = g(\mathbfit{x}). \end{equation}\] This combination of the two functions \(\delta\) and \(g\) on the left-hand is known as a convolution. We now propose that the following relationship – another convolution – holds: \[\begin{equation} \label{gensol} \int_{V}u_f(\mathbfit{x}-\mathbfit{s},t) g(\mathbfit{s})\,{\mathrm d}\mathbfit{s} = u(\mathbfit{x},t). \end{equation}\] To see this, we note that, as the domain of integration is fixed, we can take the operator inside (the Reynolds transport theorem again), \[\begin{equation} \label{lingreen} \int_{V}\mathcal{L}u_f(\mathbfit{x}-\mathbfit{s},t)g(\mathbfit{s})\,{\mathrm d}\mathbfit{s} = \mathcal{L}u(\mathbfit{x},t). \end{equation}\] Remembering that all derivatives in \(\mathcal{L}\) will be with respect to \(\mathbfit{x}\), [fundprob] and [fullprob] ensure that both sides vanish. In addition, as \(\lim_{t \to 0} u_f = \delta(\mathbfit{x})\), we will obtain the right initial condition through [deltaid]. Thus, solving [fullprob] can be reduced to solving the fundamental problem, [fundprob]: generally an easier task! Integrating [gensol] then leaves us with the full solution. This can often be of benefit for analytic solutions and asymptotic approximations.

5.4.4 Solving the diffusion equation with different initial conditions

Let’s say we wish to solve the following problem: \[\begin{equation} \label{rectheat} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} = D \mathchoice{\frac{\partial^2 c}{\partial x^2}}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2}{\partial^2 c/\partial x^2},\quad c(x,0) = \left\{\begin{array}{ll} 0 & x< -a\\ 1 & -a \leq x \leq a\\ 0 & x>a \end{array} \right. , \quad c(\pm \infty ,t) = 0. \end{equation}\] The initial condition is a ‘population’ spread evenly over a domain \(x\in[-a,a]\). Using the fundamental solution in [fundamentalsol], and [gensol], the solution to the problem in [rectheat] is \[c(x,t) = \frac{1}{ 2\sqrt{\pi Dt}}\int_{-a}^{a} \mathrm{e}^{-(x-s)^2/4Dt}\,{\mathrm d}{s} = \frac{1}{2}\left[\operatorname{erf}\left(\frac{a-x}{\sqrt{4Dt}}\right) + \operatorname{erf}\left(\frac{a+x}{\sqrt{4Dt}}\right)\right],\] where \[\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}}\int_{0}^{x}\mathrm{e}^{-s^2}{\mathrm d}{s}.\] is the error function, a function with well-known properties. The solution is shown to vary with \(t\) in 5.7.

5.5.1 The porous medium equation

Rather than carrying on with the standard one-dimensional case we now turn to a radially symmetric two-dimensional case. We will again ignore source terms (\(f\equiv 0\)) and we shall assume no velocity \(\mathbfit{v}\), so only diffusive behaviour is present. We seek a fundamental solution with particle density \(Q\), i.e., we are seeking a solution to the problem \[\begin{equation} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} = D_0\boldsymbol{\nabla}\cdot\left[\left(\frac{c}{c_0}\right)^m \boldsymbol{\nabla}c\right],\quad c(\mathbfit{x},0) = Q\delta(\mathbfit{x}), \quad c(\mathbfit{x} \to \infty,t)=0. \label{radial-diff-eqn} \end{equation}\] By \(c(\mathbfit{x} \to \infty,t)\) we mean to imply that the distribution must eventually go to zero in all directions.

In polar coordinates, \((r,\theta)\), you will remember that \[\boldsymbol{\nabla}c = \mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r}\widehat{\mathbfit{r}} + \frac{1}{r}\mathchoice{\frac{\partial c}{\partial\theta}}{\partial c/\partial\theta}{\partial c/\partial\theta}{\partial c/\partial\theta}\widehat{\mathbfit{\theta}},\] and for some vector \(\mathbfit{A} = A_r\widehat{\mathbfit{r}} + A_{\theta}\widehat{\mathbfit{\theta}}\), \[\boldsymbol{\nabla}\cdot\mathbfit{A} = \frac{1}{r}\left(\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}(rA_r) + \mathchoice{\frac{\partial A_{\theta}}{\partial\theta}}{\partial A_{\theta}/\partial\theta}{\partial A_{\theta}/\partial\theta}{\partial A_{\theta}/\partial\theta}\right).\] If we assume radial symmetry, \(c(r,\theta,t) \equiv c(r,t)\), then the right-hand side of [radial-diff-eqn] yields \[\begin{equation} \label{porusradial} D_0\boldsymbol{\nabla}\cdot\left[\left(\frac{c}{c_0}\right)^m\mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r} \widehat{\mathbfit{r}}\right] = \frac{D_0}{r}\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}\left[r\left(\frac{c}{c_0}\right)^m\mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r} \right], \end{equation}\] and our equation becomes have \[\begin{equation} \label{porous-media-eqn} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} = \frac{D_0}{r}\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}\left[r\left(\frac{c}{c_0}\right)^m\mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r} \right]. \end{equation}\] This is an example of the porous media equation. We will shortly show that the solution is \[\begin{equation} \label{porussol} c(r,t) = \left\{ \begin{array}{ll} \displaystyle\frac{c_0}{\lambda(t)}\left[1 - \left(\frac{r}{r_0 \lambda(t)}\right)^2\right]^{1/m} & \mbox{} r\leq r_0 \lambda(t), \\[6pt] 0 & \mbox{} r> r_0 \lambda(t), \end{array} \right. \end{equation}\] where \[\lambda(t) = \left(\frac{t}{t_0}\right)^{1/(2+m)},\quad t_0 = \frac{r_0^2 m}{2 D_0 (m+2)}, \quad r_0 =\frac{Q\mathit{\Gamma}(1/m + 3/2)}{\pi^{1/2}c_0 \mathit{\Gamma}(1/m + 1)},\] and where \(\mathit{\Gamma}(z)\) is the gamma function. The solution has compact support – it is zero outside of a finite radius – and this is a common property of hyperbolic PDEs. The volume, \(Q\), is conserved and we see in 5.8 the solution spreading out over time.

5.5.2 Similarity solutions

Deriving this solution involves using another technique for solving partial differential equation: a symmetry reduction technique known as looking for similarity solutions. The basic idea is that if we can find some transformation of our solution which changes its shape but still solves the basic equations, then we should be able to use this transformation to eliminate one of the independent variables. This is a bit like with Fourier transforms when this elimination turns a PDE into an ODE. Unfortunately, of course, first you must such a scaling, and this is not always easy! But for a significant class of PDEs with one spatial dimension the generic transformation of \(u(x,t)\) given by \[\begin{equation} \label{dilation} U = \lambda^{-\alpha} u, \quad X = \lambda^\beta x, \quad T = \lambda t %u(x,t) = \lambda^{\alpha} v(\lambda^{\beta}x,\lambda t) \end{equation}\] for some suitable constants \(\alpha\) and \(\beta\) will give the original PDE for \(U(X,T)\) (a dilation).

For example, consider the 1D diffusion equation from the previous chapter, \[\mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} = D\mathchoice{\frac{\partial^2 u}{\partial x^2}}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}.\] Applying [dilation] we get \[\lambda^{\alpha+1}\mathchoice{\frac{\partial U}{\partial T}}{\partial U/\partial T}{\partial U/\partial T}{\partial U/\partial T} = D\lambda^{\alpha+2\beta}\mathchoice{\frac{\partial^2 U}{\partial X^2}}{\partial^2 U/\partial X^2}{\partial^2 U/\partial X^2}{\partial^2 U/\partial X^2}.\] We can spot that if \(\beta=1/2\) then \[\mathchoice{\frac{\partial U}{\partial T}}{\partial U/\partial T}{\partial U/\partial T}{\partial U/\partial T} = D \mathchoice{\frac{\partial^2 U}{\partial X^2}}{\partial^2 U/\partial X^2}{\partial^2 U/\partial X^2}{\partial^2 U/\partial X^2},\] i.e., \(U\) satisfies the diffusion equation! So we see a transformation in the form [dilation] can lead to the same differential equation, just in \(U\) instead of \(u\).

Note that \(\alpha\) turned out to be arbitrary here: as we will see, we can pick it to suit the initial conditions.

But why do we care that under this transformation, \(U\) satisfies the same PDE? Well, note that the groupings \[\begin{equation} \eta = \frac{X}{T^\beta} = \frac{x}{t^\beta} \quad \text{and} \quad \xi = UT^\alpha = ut^{\alpha} \label{sim-sol-sub} \end{equation}\] are invariant (i.e. unchanged) under the transformation. By a theorem sometimes known as Morgan’s theorem,³ we can write one invariant as a function of the other, \[\xi = v(\eta),\] for some function \(v\). Substituting in the definitions of \(\eta\) and \(\xi\), we get \[u(x,t) = \frac{1}{t^\alpha} v \left( \frac{x}{t^\beta} \right),\] and this is the key result! This implies that the solution \(u(x,t)\) has the same shape as the ‘master shape’ \(v(\eta)\), just stretched by a factor \(t^{\beta}\) along the \(x\)-axis and squashed by a factor \(t^{\alpha}\) along the \(y\)-axis. This shape can then be appropriately scaled to give the time-varying behaviour of the system (determined by the parameters \(\alpha\) and \(\beta\), depending on the actual equation).

We have seen this behaviour already. Recall the fundamental solution of the diffusion equation, [fundamentalsol] (with \(u\) taking the \(c\) role), \[u(x,t) =\frac{1}{2 \sqrt{\pi Dt}} \mathrm{e}^{-x^2/4Dt}.\] We have seen above that for the diffusion equation, \(\beta=1/2\), so if we write \(\eta = x/t^{1/2}\) and \(v = ut^{\alpha}\) (from [sim-sol-sub]) then \[v(\eta,t) =\frac{t^{\alpha}}{2\sqrt{\pi D t}} \mathrm{e}^{-\eta^2/4D}.\] We said earlier that \(\alpha\) can be set to fit the initial conditions. Well, the fundamental solution has a very specific initial condition which gives form to this solution, and here if we set \(\alpha=1/2\) then \[v(\eta) =\frac{1}{2 \sqrt{\pi D}}\mathrm{e}^{-\eta^2/4D}.\]

So we see our solution \(u(x,t)\) can be written in the form \[u(x,t) = \frac{1}{t^\alpha}v\left(\frac{x}{t^{\beta}}\right) = \frac{1}{t^{1/2}}v\left(\frac{x}{t^{1/2}}\right),\] and the solution will still satisfy the diffusion equation.

As we saw in 5.6, this confirms that the solution \(u(x,t)\) has the same shape as the master shape \(v(\eta)\), stretched by a factor \(t^{1/2}\) along the \(x\)-axis and squashed by a factor \(t^{1/2}\) along the \(y\)-axis.

5.5.3 Back to the porous media equation

Let’s go back to [porous-media-eqn], \[\begin{equation} \label{radialpde} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} = \frac{D_0}{r}\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}\left[r \left(\frac{c}{c_0}\right)^m \mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r}\right], \end{equation}\] and, using the observed similarity solution, make the substitution \[c(r,t) = \frac{1}{t^{\alpha}}v\left(\frac{r}{t^{\beta}}\right).\] Denoting \(\eta=r/t^{\beta}\), so that \(v=v(\eta)\), we first change the left-hand side: \[\begin{align} \mathchoice{\frac{\partial c}{\partial t}}{\partial c/\partial t}{\partial c/\partial t}{\partial c/\partial t} &= -\alpha\frac{1}{t^{\alpha+1}}v + \frac{1}{t^{\alpha}}\mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}\mathchoice{\frac{{\mathrm d}\eta}{{\mathrm d}t}}{{\mathrm d}\eta/{\mathrm d}t}{{\mathrm d}\eta/{\mathrm d}t}{{\mathrm d}\eta/{\mathrm d}t} \\ &= -\alpha\frac{1}{t^{\alpha+1}}v - \frac{r \beta}{t^{\alpha}t^{\beta+1}}\mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}\\ \label{fincaldcdt}& = -\alpha\frac{1}{t^{\alpha+1}}v - \frac{\eta \beta}{t^{\alpha+1}}\mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}. \end{align}\] On the right-hand side, \[\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r} = \frac{1}{t^\beta}\mathchoice{\frac{\partial}{\partial\eta}}{\partial/\partial\eta}{\partial/\partial\eta}{\partial/\partial\eta} \implies \mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r}= \frac{1}{t^{\alpha+\beta}}\mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta},\] and so we have \[\begin{align} \frac{D_0}{r}\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}\left[r \left(\frac{c}{c_0}\right)^m \mathchoice{\frac{\partial c}{\partial r}}{\partial c/\partial r}{\partial c/\partial r}{\partial c/\partial r}\right]&=\frac{D_0}{\eta t^{2\beta}}\mathchoice{\frac{{\mathrm d}}{{\mathrm d}\eta}}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}\left[\frac{t^\beta \eta}{t^{\alpha m}}\left(\frac{v}{c_0}\right)^m \frac{1}{t^{\alpha+\beta}}\mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}\right]\\ \label{finalrhs}&= \frac{1}{t^{\alpha(m+1)+2\beta}}\frac{D_0}{\eta}\mathchoice{\frac{{\mathrm d}}{{\mathrm d}\eta}}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}\left[\eta\left(\frac{v}{c_0}\right)^m \mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}\right]. \end{align}\] Comparing [fincaldcdt] and [finalrhs] we see we want to set the \(t\)-exponents equal to each other, \[\alpha + 1 = \alpha(m+1) + 2\beta \implies 1 = \alpha m + 2\beta.\] We can multiply through by \(t^{\alpha+1}\) to obtain \[\begin{equation} \alpha v + \beta \mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}\eta +\frac{D_0}{\eta}\mathchoice{\frac{{\mathrm d}}{{\mathrm d}\eta}}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}\left[\eta\left(\frac{v}{c_0}\right)^m\mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}\right] =0. \label{rel-between-alpha-and-beta} \end{equation}\] So now we have a relationship between \(\alpha\) and \(\beta\), but still some freedom in choosing one of them.

We now manipulate [rel-between-alpha-and-beta] to simplify it with an appropriate choice of \(\alpha\). In particular, we want to write our equation in the form \(\mathchoice{\frac{{\mathrm d}[\cdots]}{{\mathrm d}\eta}}{{\mathrm d}[\cdots]/{\mathrm d}\eta}{{\mathrm d}[\cdots]/{\mathrm d}\eta}{{\mathrm d}[\cdots]/{\mathrm d}\eta}\) so that we can then integrate and create a first-order equation.

Multiply our equation through by \(\eta\), \[\begin{equation} \label{newode} \alpha v \eta + \beta \mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta} \eta^2 + D_0\mathchoice{\frac{{\mathrm d}}{{\mathrm d}\eta}}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}\left[\eta\left(\frac{v}{c_0}\right)^m \mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}\right]=0. \end{equation}\] The idea is to turn the first two terms into a derivative so we can integrate. You might notice that for some \(\gamma\), \[\gamma\mathchoice{\frac{{\mathrm d}}{{\mathrm d}\eta}}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}[v\eta^2] = \gamma \mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta} \eta^2 + 2\gamma v\eta,\] which is almost the form we want. If we set \(\beta=\gamma\), then we have two simultaneous equations, \[\begin{aligned} 1 & = \alpha m + 2 \gamma,\\ \alpha & = 2\gamma. \end{aligned}\] This gives us \(\alpha = 1/(m+1)\) and \(\gamma=1/[2(m+1)]\), so [newode] becomes \[\frac{1}{2(m+1)}\mathchoice{\frac{{\mathrm d}}{{\mathrm d}\eta}}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}[v\eta^2] + D_0\mathchoice{\frac{{\mathrm d}}{{\mathrm d}\eta}}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}{{\mathrm d}/{\mathrm d}\eta}\left[\eta\left(\frac{v}{c_0}\right)^m \mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}\right] = 0.\] Integrating with respect to \(\eta\) gives \[\frac{1}{2(m+1)}v\eta^2 + D_0\eta\left(\frac{v}{c_0}\right)^m \mathchoice{\frac{{\mathrm d}v}{{\mathrm d}\eta}}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta}{{\mathrm d}v/{\mathrm d}\eta} = C,\] with \(C\) a constant of integration.

We choose the imposed boundary condition that \(v(\eta) \to 0\) sufficiently fast as \(\eta(r) \to \infty\) which implies the constant of integration is \(C=0\). For any finite \(t\), this is equivalent to saying the population vanishes as \(r\to \infty\).

We solve the differential equation using separation of variables, \[\begin{aligned} \int \frac{\eta}{2(m+1)}{\mathrm d}{\eta} &= - \int \frac{D_0}{c_0^m}v^{m-1}{\mathrm d}{v},\\ \implies \frac{\eta^2}{4(m+1)} &=- \frac{D_0}{m c_0^m}v^{m} + A,\\ \implies v &= \frac{m^{1/m} c_0}{D_0^{1/m}}\left(A - \frac{\eta^2}{4(m+1)}\right)^{1/m}. \end{aligned}\] Finally, using \(c= v/t^\alpha\) and \(\eta= r/t^{\beta}\), we get \[c(r,t)= \frac{m^{1/m} c_0}{t^{\alpha}D_0^{1/m}}\left(A - \frac{r^2}{t^{2\beta}4(m+1)}\right)^{1/m},\] where \(\alpha = 1/(m+1)\) and \(\beta=1/[2(m+1)]\).

To get \(A\) and hence the final solution, [porussol], we enforce the fact that the density should be constrained for all time to be constant (it isn’t necessarily conserved in this more complex system). This is a little tedious and I don’t really want you to spend your time doing this. The main concept we’re interested in here is the method of similarity solutions.

8.1.1 A brief review of bifurcation theory

The classification of equilibria above in terms of (linearly) stable or unstable is especially important when we think of ODE models which depend on parameters, and think of how solutions change as parameters are varied. If an equilibrium becomes stable or unstable as a parameter is varied, we refer to this as a bifurcation, and these can come in many different flavours. One can develop a fairly sophisticated theory of bifurcations of systems of the form, \[\mathchoice{\frac{{\mathrm d}\mathbfit{u}}{{\mathrm d}t}}{{\mathrm d}\mathbfit{u}/{\mathrm d}t}{{\mathrm d}\mathbfit{u}/{\mathrm d}t}{{\mathrm d}\mathbfit{u}/{\mathrm d}t} = \mathbfit{f}(\mathbfit{u};\bm{p}),\] where now \(\bm{p} = (p_1,p_2,\dots,p_m)\) is a vector of parameters. You saw several examples of these bifurcations last term, and we will see more this term (and extend the theory to bifurcations of PDEs). Rather than pursue a general treatment, let’s review one particularly important simple example which will play a role later.

Consider the ODE \[\begin{equation} \label{pitchfork-US-ODE} \mathchoice{\frac{{\mathrm d}u}{{\mathrm d}t}}{{\mathrm d}u/{\mathrm d}t}{{\mathrm d}u/{\mathrm d}t}{{\mathrm d}u/{\mathrm d}t} = u(p-u^2), \end{equation}\] where \(p\) is a real parameter. If \(p \leq 0\), then \(u_0=0\) is the only steady state (as we are only working with real functions), and this solution is stable for \(p < 0\). You should check these linear stability calculations yourself. If instead \(p>0\) then (in addition to the zero steady state), \(u_0=\pm \sqrt{p}\) are two new steady states. You can also show that these steady states are always stable (when they exist), and that \(u_0=0\) is unstable for \(p>0\). Hence this zero solution has lost stability by generating two symmetric branches of solutions. The corresponding bifurcation diagram is shown in 1.1, where the name ‘pitchfork bifurcation’ comes from.

This particular bifurcation will come up again this term, and ideally seeing this diagram will help remind you how to read these kinds of plots. Having drawn something like 1.1, we can read off essentially all of the local behaviour of the system (as this is a 1-D system which are always governed by equilibria). If \(p<0\), then the only possibility is that the system approaches \(u_0=0\). If instead \(p>0\), then the steady state \(u_0=0\) is a separatrix so that if \(u(0)>0\), the solution of the ODE approaches \(u_0=\sqrt{p}\), and if \(u(0)<0\) then as \(t \to \infty\), \(u(t) \to -\sqrt{p}\). [pitchfork-US-ODE] can be solved exactly, but the explicit solution itself is less insightful than the diagram. More importantly, this kind of bifurcation arises in a number of analytically intractable models, but linear stability analysis (and hence sketching the full bifurcation diagram) allows us to get important insights into these kinds of systems.

8.1.2 Generalisations and limitations

While this basic approach is powerful, allowing the analysis of complex nonlinear models which cannot be analytically solved, it also has important limitations. This kind of stability analysis is ‘local’ so that it is strictly only valid when \(u\approx u_0\). For an unstable equilibrium (where the solution moves away from \(u_0\)) or for an initial condition not very close, linear stability often does not tell us what happens. Beyond one dimension of state space (a scalar ODE), there are also periodic solutions which can be stable or unstable, and they can undergo bifurcations as parameters vary, though the analysis is harder than for steady state solutions. In three or more dimensions, chaotic solutions can be found, and this kind of local analysis tells us essentially nothing in these cases.

Still, this is a general tool which is widely applicable. In particular, the presentation above only considers first-order systems, but last term you considered ODEs of the form, \[\begin{equation} \label{higher-US-order-US-ODE} \mathchoice{\frac{{\mathrm d}^{n} u}{{\mathrm d}t^{n}}}{{\mathrm d}^{n} u/{\mathrm d}t^{n}}{{\mathrm d}^{n} u/{\mathrm d}t^{n}}{{\mathrm d}^{n} u/{\mathrm d}t^{n}}+F\left (u, \mathchoice{\frac{{\mathrm d}u}{{\mathrm d}t}}{{\mathrm d}u/{\mathrm d}t}{{\mathrm d}u/{\mathrm d}t}{{\mathrm d}u/{\mathrm d}t}, \dots, \mathchoice{\frac{{\mathrm d}^{n-1} u}{{\mathrm d}t^{n-1}}}{{\mathrm d}^{n-1} u/{\mathrm d}t^{n-1}}{{\mathrm d}^{n-1} u/{\mathrm d}t^{n-1}}{{\mathrm d}^{n-1} u/{\mathrm d}t^{n-1}} \right)=0. \end{equation}\] We can rewrite such equations as first-order equations by creating new variables for each time derivative; that is, \[\begin{equation} \label{higher-US-oder-US-system} \mathchoice{\frac{{\mathrm d}u}{{\mathrm d}t}}{{\mathrm d}u/{\mathrm d}t}{{\mathrm d}u/{\mathrm d}t}{{\mathrm d}u/{\mathrm d}t} = u^1, \quad \mathchoice{\frac{{\mathrm d}^{2} u}{{\mathrm d}t^{2}}}{{\mathrm d}^{2} u/{\mathrm d}t^{2}}{{\mathrm d}^{2} u/{\mathrm d}t^{2}}{{\mathrm d}^{2} u/{\mathrm d}t^{2}}=\mathchoice{\frac{{\mathrm d}u^1}{{\mathrm d}t}}{{\mathrm d}u^1/{\mathrm d}t}{{\mathrm d}u^1/{\mathrm d}t}{{\mathrm d}u^1/{\mathrm d}t} = u^2,\,\, \dots, \,\, \mathchoice{\frac{{\mathrm d}^{n-1} u}{{\mathrm d}t^{n-1}}}{{\mathrm d}^{n-1} u/{\mathrm d}t^{n-1}}{{\mathrm d}^{n-1} u/{\mathrm d}t^{n-1}}{{\mathrm d}^{n-1} u/{\mathrm d}t^{n-1}} = \mathchoice{\frac{{\mathrm d}u^{n-2}}{{\mathrm d}t}}{{\mathrm d}u^{n-2}/{\mathrm d}t}{{\mathrm d}u^{n-2}/{\mathrm d}t}{{\mathrm d}u^{n-2}/{\mathrm d}t} = u^{n-1}, \quad \mathchoice{\frac{{\mathrm d}u^{n-1}}{{\mathrm d}t}}{{\mathrm d}u^{n-1}/{\mathrm d}t}{{\mathrm d}u^{n-1}/{\mathrm d}t}{{\mathrm d}u^{n-1}/{\mathrm d}t} = F\left (u, u^1, \dots, u^{n-1} \right). \end{equation}\] If you recall from last term, we can linearise [higher-US-order-US-ODE] by first finding a steady state where all time derivatives are zero. Hence, this is a constant \(u_0\) satisfying \(F(u_0,0,\dots,0)=0\). Substituting in \(u = u^1_0 + \varepsilon u^1_1(t)\), we get a linear ODE in \(u^1_1\) and \(n\) of its derivatives. Solving this ODE via the ansatz \(u^1 = \exp(\lambda t)\) leads to an \(n\)th degree polynomial for the growth rates \(\lambda\), which determine the behaviour of the perturbation. This characteristic polynomial for \(\lambda\) turns out to be exactly what we would get by considering the eigenvalues of the Jacobian of the system [higher-US-oder-US-system].

8.3.1 Food webs & ecological networks

Last term you studied different kinds of populations and their interactions. These included single-species phenomena such as different kinds of growth and Allee effects, as well as multi-species interactions such as predation and competition. In general, ecological systems can be extremely complex, with hundreds or thousands of different populations interacting. Quite a lot of mathematical modelling is about how to think of all of these complex interactions, and what sensible assumptions we might use to simplify such models to their essential components in order to learn something about how these ecosystems behave. We often organise these interactions into networks (also known as graphs⁹) of pairwise effects, such as competition or predation. Such networks are known as food webs, though it is typical that such webs have a layered structure divided into trophic levels; see 1.2 for an example. This is a way of thinking of where energy flows in an ecosystem, although as can be seen, these levels are not always cleanly separated.

9.1.1 Exponential growth

Let’s consider one of the simplest population models given by \(f(u) = ru\) for some growth rate \(r>0\). This should remind you of the very first model of growth covered in the course. Our population density \(u\) then evolves according to \[\begin{equation} \label{exp-US-diff-US-eqn} u_n = ru_{n-1} \implies u_n = r^{n}u_0, \end{equation}\] where we have iterated the system to write down the dependence of \(u_n\) on \(r\) and the initial population, \(u_0\). We see that this model behaves differently depending on \(r\): for \(r<1\), the population decays towards \(u_n=0\) as \(n \to \infty\), whereas for \(r>1\) the population grows exponentially in time. In the case of \(r=1\), the model gives that every real number \(u_0\geq 0\) is an equilibrium point. We have seen such a situation before in ODE models when the Jacobian matrix had a zero eigenvalue. Biologically, we would not expect a parameter to have an exact value, so we can discount the case \(r=1\) and instead just view it as the bifurcation point between growth and extinction.

Phase portraits and other kinds of graphs were valuable for studying ODE models, so we introduce a similar way of understanding scalar difference equations called the cobweb plot. See 2.1 for an example involving the exponential growth model. One iteration of the map given in [exp-US-diff-US-eqn] corresponds to moving vertically to the green line given by \(f(u)=ru\) from the initial condition on the \(x\)-axis. To iterate again, one moves horizontally to the line of slope one, i.e. to the blue line, before again moving vertically to the function \(f(u)\). These plots can take some time to understand, but they can help immediately visualise the change in the stability of \(u_0=0\) as \(r\) crosses the threshold of \(1\).

Finally, note that this model is linear, and so it is perhaps not surprising that we could solve it directly. This will not be true of generic models, but through linear stability analysis, we can study properties of equilibria by reducing our models to this one.

9.1.2 The logistic map

As we did at the start of Michaelmas, we will now consider a simple population model involving logistic growth, given by \[\begin{equation} \label{logistic-US-map} u_n = ru_{n-1}(1-u_{n-1}), \end{equation}\] where we have already nondimensionalised the carrying capacity to be \(1\). This model is often called the logistic map, due to its similarity with the logistic model of population growth. Despite this similarity to the model studied previously, we will see that it has both biological and mathematical properties which are drastically different from its continuous-time counterpart.

Immediately we see that if \(u_{n-1}>1\) at any iteration \(n\), then \(u_n<0\). This means that we only consider initial conditions \(u_0 \in [0,1]\). For all \(r\) with \(0\leq r \leq 4\), this map takes points in \([0,1]\) to points in \([0,1]\) at each iteration (you should be able to show this), and for \(r\) outside of this range, some iterations will take negative values. Therefore, we will restrict our range of \(r\) values to be in \([0,4]\) interval, as we are looking for biologically-meaningful models.

We can start to get an idea of how this model behaves by looking at equilibria and their stability. These satisfy \[\begin{equation} \label{log-US-map-US-equilibria} u^* = ru^*(1-u^*) \implies u^* = 0 \,\,\textrm{ or }\,\, u^* = \frac{r-1}{r}. \end{equation}\] We see that \(u_0=0\) is always a permissible equilibrium, whereas the nonzero solution is only distinct and permissible for \(r>1\).

We cannot solve [logistic-US-map] in general, but we can perform a linear stability analysis as we did for differential equations, by thinking of small perturbations of initial conditions and asking if they grow or decay. We use the ansatz \(u_n = u^* + \varepsilon u_{1,n}\), noting that \(u^*\) is a constant but \(u_{1,n}\) is a sequence in \(n\) (so depends on time). Substituting this ansatz into [logistic-US-map] we find \[u^*+\varepsilon u_{1,n} = r(u^*+\varepsilon u_{1,n-1})(1-u^*-\varepsilon u_{1,n-1}) = ru^*(1-u^*) + \varepsilon r(1-2u^*)u_{1,n-1} - \varepsilon^2 r(u_{1,n-1})^2.\] We notice that the \(\mathcal{O}(1)\) equation, given by \(u^* = ru^*(1-u^*)\), is exactly the equation of an equilibrium, and so will identically cancel out. We can then divide by \(\varepsilon\), and neglect the highest-order terms in \(\varepsilon\) to arrive at the linearised equation, \[\begin{equation} \label{lin-US-log-US-map} u_{1,n} = r(1-2u^*)u_{1,n-1}. \end{equation}\] This equation is exactly [exp-US-diff-US-eqn], though the constant \(r\) in that model is now \(r(1-2u^*)\), and so we know by solving the linear model our equilibrium is stable if \(|r(1-2u^*)|<1\), and unstable if \(|r(1-2u^*)|>1\). For the extinction equilibrium, \(u^*=0\), we then have stability when \(r<1\), and the nonzero equilibrium exists for \(r>1\).

Plugging in the nonzero equilibrium into our stability condition, we find that it is stable given the condition, \[\left|r\left(1-2\frac{r-1}{r}\right)\right|< 1 \implies -1 < -r+2 < 1,\] where we are just writing \(|-r+2|< 1\) as the two separate inequalities, and simplifying the terms. Recall that there were two different routes for [exp-US-diff-US-eqn] to become unstable, corresponding to one side or the other of these inequalities being violated. Once either is no longer satisfied, then \(u_0\) is unstable. We see that, since we are only considering stability given permissibility, \(r>1\) means that the inequality can only break when \(-1 = -r+2\), i.e. when \(r=3\). We can get some insight into this by drawing some cobweb diagrams for \(r\) in these various ranges, which we do in 2.2. I would strongly encourage you to use the links in the caption to explore this model with different initial conditions and values of \(r\).

As expected from our linear stability analysis, we see convergence to \(u_0=0\) for \(r<1\), and convergence to another equilibrium for \(r \in (1,3)\), corresponding exactly to the intersection of the blue line and the green curve. One qualitative change occurs at \(r=2\), where the convergence to the equilibrium goes from being monotonic to oscillatory. This can be predicted from the linear theory, as \(r=2\) is exactly when the coefficient of \(u_{1,n-1}\) in [lin-US-log-US-map] goes from being positive to negative.

What happens at \(r=3\) and beyond? This is a different kind of bifurcation from what occurred at \(r=1\) for the extinction state, as there are no other equilibria for the system to go to. The cobweb plot in 2.2(d) is also not too informative, though direct evaluations can be helpful for getting an idea of what happens. Plugging in \(u_0=0.5\) for \(r=3.2\) we find the sequence, \[u_1=0.8,\, u_2 = 0.512,\, u_3 =0.7995392,\, u_4 = 0.51288405652,\, u_5 \approx 0.79946880348,\, \dots,\] where by \(u_2\), we see that the sequence has converged rapidly to a periodic solution.

We can even find the values of this periodic oscillation directly. What would such a solution look like? It would be some value \(u\) such that, after two iterations, it was unchanged. That is, we are looking for equilibria of the iteration \(u_n = f^2(u_{n-2}) = f(f(u_{n-2}))\). Any equilibrium of the model, i.e. those given in [log-US-map-US-equilibria], is also an equilibrium of this equation, but there may be new ones. In fact we can find them without too much trouble, as they would satisfy the quartic equation \[u^* = f^2(u^*) = r^2u^*(1-u^*))(1-ru^*(1-u^*))),\] where we can find the roots (e.g. by dividing out the two roots we know) to find \[u^* = \frac{1 + r + \sqrt{r^2 - 2 r - 3}}{2 r}, \, \frac{1 + r - \sqrt{r^2 - 2 r - 3}}{2 r}.\] Plugging in \(r=3.2\) these give \(u^* \approx 0.51304\dots\) and \(u_0=0.79945\dots\), so our numerical intuition above was not wrong in this case.

Rather than dwell too long on this example, let’s go back to questions of biology. Do we expect populations to exhibit complicated period-doubling oscillations, which may be somewhat sensitively dependent on their growth rate? Do we expect large populations to lead to negative populations? Really these are more artefacts of the particular form of [logistic-US-map], which is, unlike its continuous counterpart, often not considered a good model of real populations. In the problem sheets you will explore a few examples of possibly realistic population models of this kind, though that isn’t to say that complex oscillations and chaos are not possible kinds of biological dynamics.

9.2.1 Fibonacci, rabbits, and unit conversions

Leonardo of Pisa (given the more well-known name Fibonacci in the 1800s) wrote a book on arithmetic in 1202 which contained the following exercise: After one year, how many rabbits can a single breeding pair produce in total? The following assumptions about rabbit biology are made: 1. Rabbits take one month to mature to reproductive age. 2. Once they reach reproductive age, a pair of rabbits will produce a new pair (one male and one female) each month. One can then graphically (see 2.3(a)) or mentally depict this process over a 12 month period to deduce the following sequence of pairs of rabbits: \[1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144,\] which are the first twelve terms of the Fibonacci sequence¹², where each subsequent term is the sum of the previous two. In essence this captures the two simplistic biological assumptions about the population.

Now, if we consider a difference-equation formulation of this model¹³, we have: \[\begin{equation} \label{fib-US-eqn} F_n = F_{n-1}+F_{n-2}, \quad F_0=0, \quad F_1=1. \end{equation}\] Importantly, this is a linear difference equation, and so we might expect solutions like those seen for the exponential model, [exp-US-diff-US-eqn], to work. That is, we can look for solutions of the form \(F_n = \lambda^n\), and see what values of \(\lambda\) would solve [fib-US-eqn].

Plugging this ansatz in, we can reduce the resulting equation to a quadratic in \(\lambda\) by dividing through by \(\lambda^{n-2}\) to find: \[\begin{equation} \label{fib-US-lambdas} \lambda^n = \lambda^{n-1}+\lambda^{n-2} \implies \lambda^2 = \lambda+1 \implies \lambda = \frac{1\pm \sqrt{5}}{2}. \end{equation}\] As [fib-US-eqn] is a linear equation, we expect that any superposition of these two values of lambda would be a solution. That is, the general solution takes the form, \[F_n = C_1\left(\frac{1+ \sqrt{5}}{2} \right)^n+C_2\left(\frac{1- \sqrt{5}}{2} \right)^n.\] Finally, to determine the values of \(C_1\) and \(C_2\) we use the initial data that \(F^0=0\), \(=F^1=1\) to find, \[C_1+C_2=0, \, \, C_1\left(\frac{1+ \sqrt{5}}{2} \right)+C_2\left(\frac{1- \sqrt{5}}{2} \right) =1 \implies C_1=\frac{1}{\sqrt{5}}, \,\, C_2 = \frac{-1}{\sqrt{5}}.\] This gives the following formula for the \(n\)th Fibonacci number as, \[F_n = \frac{1}{\sqrt{5}}\left[\left(\frac{1+\sqrt{5}}{2}\right)^n-\left(\frac{1-\sqrt{5}}{2}\right)^n\right] = \frac{1}{2^n\sqrt{5}}\left[(1+\sqrt{5})^n-(1-\sqrt{5})^n\right].\]

The numbers \(F_n\), as well as their ratio \(F_n/F_{n-1}\), which rapidly approaches the golden ratio \((1+\sqrt{5})/2\), appear in a huge variety of natural and cultural phenomena. These range from the use of the golden ratio in art, to the appearance of Fibonacci numbers in the number of seeds on the top of many flowers (see 2.3(b) or search for ‘Fibonacci phyllotaxis’), and even in converting miles to kilometres (see 2.1). However, from a biological modelling point of view, it should be clear that linear models will only be very rough caricatures of real populations. Still, they can be useful in making rough predictions, such as estimating the growth rate of populations in the absence of predation or resource limitations.

9.2.2 General linear stability for systems

Let’s now state the general procedure for linear stability analysis of systems of the form, \[\begin{equation} \label{gen-US-diff-US-sys} \mathbfit{u}_n = \mathbfit{f}(\mathbfit{u}_{n-1}), \end{equation}\] where we are thinking of \(\mathbfit{u}_n\) having \(N\) components. The first step is to compute any equilibria. These are constant vectors \(\mathbfit{u}^*\) which satisfy \(\mathbfit{u}^* = \mathbfit{f}(\mathbfit{u}^*)\). Once such an equilibrium is found, we can perturb it in the usual way by writing \(\mathbfit{u}_n = \mathbfit{u}^*+\varepsilon \mathbfit{u}_{1,n}\) to find, \[\mathbfit{u}^*+\varepsilon \mathbfit{u}_{1,n} = \mathbfit{f}(\mathbfit{u}^*+\varepsilon \mathbfit{u}_{1,n-1}) \approx \mathbfit{f}(\mathbfit{u}^*)+\varepsilon \mathsfbfit{J}\mathbfit{u}_{1,n-1} + O(\varepsilon^2),\] where \(\mathsfbfit{J}\) is the usual Jacobian matrix of \(\mathbfit{f}\) evaluated at the steady state \(\mathbfit{u}^*\). Cancelling the \(O(1)\) terms (that is, using the steady state equation \(\mathbfit{u}^* = \mathbfit{f}(\mathbfit{u}^*)\), dividing through by \(\varepsilon\), and then setting \(\varepsilon=0\), we arrive at the linear system, \[\mathbfit{u}_{1,n} = \mathsfbfit{J}\mathbfit{u}_{1,n-1}.\] This system is in many ways similar to the linearised ODE system [gen-US-lin-US-ODE], in that its solutions are essentially determined by the eigenvalues of \(\mathsfbfit{J}\).

For simplicity, let’s consider the case when \(\mathsfbfit{J}\) is diagonalisable as \(\mathsfbfit{J}=\mathsfbfit{P}\mathsfbfit{D}\mathsfbfit{P}^{-1}\) where \(\mathsfbfit{P}\) is the matrix of eigenvectors and \(\mathsfbfit{D}\) is a diagonal matrix of the eigenvalues. We then have, \[\begin{equation} \label{lin-US-stab-US-diff-US-sys} \mathbfit{u}_{1,n} = \mathsfbfit{J}\mathbfit{u}_{1,n-1} = \mathsfbfit{J}^n\mathbfit{u}_{1,0} = \mathsfbfit{P}\mathsfbfit{D}^n\mathsfbfit{P}^{-1}\mathbfit{u}_{1,0}, \end{equation}\] where we have used the matrix factorisation to cancel out each term in the product, so that only the diagonal matrix \(\mathsfbfit{D}\) gets raised to the power \(n\). Finally, we make a change of basis of the form \(\mathbfit{v}_n = \mathsfbfit{P}^{-1}\mathbfit{u}_{1,n}\) and multiply [lin-US-stab-US-diff-US-sys] by \(\mathsfbfit{P}^{-1}\) to find \[\mathbfit{v}_{n} = \mathsfbfit{D}\mathbfit{v}_{n-1}=\mathsfbfit{D}^n\mathbfit{v}_{0}.\] This is essentially the discrete analogue of the ODE system [diagonalizable-US-lin-US-ode], and the manipulations used about the diagonalisability of \(\mathsfbfit{J}\) are the same. So if we write the components of this vector as \(\mathbfit{v}_n = (v_{1,n},v_{2,n},\dots,v_{N,n})\), these equations are of the form, \[v_{i,n} = \lambda_iv_{i,n-1} = \lambda_i^nv_{i,0},\] where \(\lambda_i\) are the eigenvalues of \(\mathsfbfit{J}\).

The perturbations \(\mathbfit{u}_{1,n}\) are just linear combinations of these \(\mathbfit{v}_n\), so they will grow or decay depending on the size of the eigenvalues \(\lambda_i\). In particular we have the following result:

9.2.3 Fibonacci, rabbits, and unit conversions

Let’s consider a two-species version of [logistic-US-map] corresponding to two competing species which have the same growth rate \(r\) and carrying capacity of \(1\) in the absence of competition. This takes the form, \[\begin{align} \label{LV-US-Diff} u_n = u_{n-1}(r - u_{n-1} - \gamma_1v_{n-1}), \\ v_n = v_{n-1}(r - v_{n-1} - \gamma_2u_{n-1}), \end{align}\] where all parameters are positive. There are four equilibria to this model, which you should be able to find as: \[\begin{equation} \label{diff-US-sys-US-equil} (u^*,v^*) = (0,0),\, (r-1,0),\, (0,r-1),\, \textrm{ or } \,\left(\frac{(\gamma_1-1)(r-1)}{\gamma_1 \gamma_2-1},\frac{(\gamma_2-1)(r-1)}{\gamma_1 \gamma_2-1} \right). \end{equation}\] We compute the Jacobian as, \[\mathsfbfit{J}(u^*,v^*) = \begin{pmatrix} r-2u^*-\gamma_1v^* & -\gamma_1u^*\\ -\gamma_2v^* & r-2v^*-\gamma_2u^* \end{pmatrix}.\]

We consider the extinction steady state, \((0,0)\) first, which is always permissible. The Jacobian is \[\mathsfbfit{J}(0,0) = \begin{pmatrix} r & 0\\ 0 & r \end{pmatrix}.\] which has repeated eigenvalues \(\lambda = r\), so this equilibrium is only stable if \(r<1\).

Next we consider the case of a single species surviving. Without loss of generality, we take \((u^*,v^*)=(r-1,0)\), which is only permissible if \(r>1\). The Jacobian is \[\mathsfbfit{J}(r-1,0) = \begin{pmatrix} 2-r & \gamma_1(1-r)\\ 0 & r-\gamma_2(r-1) \end{pmatrix},\] which has eigenvalues \(\lambda=2-r\) and \(r-\gamma_2(r-1)\). Analysing just the first eigenvalue, we see that this state is unstable for \(0\leq r <1\), and for \(r>3\). So if this state is stable, it will be between \(1 < r < 3\). Note that this is exactly the instability we identified in [logistic-US-map], corresponding to logistic growth leading to period-doubling, and hence setting \(r>3\) would only lead to time-dependent solutions for the single species \(u_n\), rather than invasion or coexistence with the species \(v_n\), at least via this particular bifurcation.

We can manipulate conditions for the second eigenvalue as follows. If \(\lambda > -1\), this says that \[r+\gamma_2(1-r) > -1 \implies \frac{1+r}{r-1} > \gamma_2,\] where all terms are positive for permissible values of \(r>1\). If instead we have \(\lambda < 1\), then \[r+\gamma_2(1-r) < 1 \implies \gamma_2>1,\] so that for this state to be stable we must have, \[\frac{1+r}{r-1} > \gamma_2 > 1,\] giving rise to two possible instabilities as \(\gamma_2\) varies. The bifurcation at \(\gamma_2=1\) has an immediate interpretation that a single species can prevent invasion by a new species as long as it competes sufficiently strongly against the second species. That is, if the population \(u\) is at its equilibrium value of \(r-1\), the introduction of some population \(v\) will just lead to the extinction of \(v\) as long as \(\gamma_2>1\). Something qualitatively similar happens in the ODE model, and is related to competitive exclusion. The other route to instability, when \(\gamma_2>(r+1)/(r-1)\) is different, and corresponds to an oscillatory instability (as here \(\lambda < -1\)) in the perturbation of \(v_n\). However, this is problematic as \(v^*=0\) at this steady state, so this would lead to dynamics which are not permissible.

In principle we can consider the coexistence state given by the fourth value of [diff-US-sys-US-equil]. We then have the Jacobian (after simplifying some of the terms): \[\mathsfbfit{J}(u^*,v^*) = \begin{pmatrix} \displaystyle r-\frac{(r-1)(\gamma_1(3-\gamma_2)-2)}{\gamma_1 \gamma_2-1}& \displaystyle -\gamma_1\frac{(\gamma_1-1)(r-1)}{\gamma_1 \gamma_2-1}\\[5mm] \displaystyle -\gamma_2\frac{(\gamma_2-1)(r-1)}{\gamma_1 \gamma_2-1} & \displaystyle r-2\frac{(r-1)(\gamma_2(3-\gamma_1)-2)}{\gamma_1 \gamma_2-1} \end{pmatrix}.\] One can make some progress studying this system, but it is extremely algebraically tedious, so we will not pursue this here (or expect you to be able to do a calculation this messy).

One advantage to difference equations compared to differential equations is how easy they are to simulate. In 2.4, we consider solutions near \(r=3\), validating the previous comment that this instability really is the same single-species period-doubling bifurcation as before. We can essentially ignore the behaviour of \(v_n\) in this case. On the other hand, 2.5 shows different values of \(\gamma_1\) and \(\gamma_2\) for which coexistence, competitive exclusion, and (unphysical) observations occur, corresponding to the two bifurcations identified above. Although we have not analysed the coexistence equilibrium, we observe that it appears to be stable for \(\gamma_1<1\) and \(\gamma_2<1\), as in the ODE case. Note that in panel (d), the equations for \(u_n\) and \(v_n\) are the same, and hence we have that both \((r-1,0)\) and \((0,r-1)\) are stable. Finally, we note that panel (e) is still within the stable parameter regime of the equilibrium \((r-1,0)\), it exhibits unphysical negative oscillations due to the negative eigenvalue, despite the fact that these oscillations eventually tend towards zero.

There are a number of differences in the analysis of this system, compared to the ODE version of a competitive Lotka–Volterra system. Here it is possible for the coexistence equilibrium to be stable, for instance, if there are insufficient resources for either population to persist (\(r<1\)). The ability of one population at equilibrium to suppress the other population is similar to the ODE model, with only minor algebraic differences in the conditions. However the coexistence equilibrium is substantially more involved to analyse, and in fact some of its instabilities will lead to complex periodic and chaotic solutions which are beyond our scope here.

10.1.1 Spatially homogeneous equilibria – survival in a closed environment

For no-flux boundary conditions (\(a_0=a_L=0\)), any constant solution to \(f(u_0)=0\) satisfies the PDE [rd-US-eqn]. In this case, we can solve [rd-US-eqn-US-linear] by assuming a solution in terms of eigenfunctions of the Laplacian. In this case, these will be the familiar cosine solutions you should be familiar with from last term.

To explore this in detail, let’s consider a model of logistic growth in a spatially-extended environment: \[\begin{equation} \label{rd-US-eqn-US-logistic} \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} = D\mathchoice{\frac{\partial^2 u}{\partial x^2}}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2} + ru\left (1-\frac{u}{K} \right), \quad x \in [0,L], \,\, t > 0, \end{equation}\] again with the no-flux boundary conditions. We could nondimensionalise \(x, t\), and \(u\) to set three of these constants to \(1\) (or to set the domain to be \(x \in [0,1]\)), but for now will study the dimensional model. The two spatially homogeneous equilibria are \(u_0=0,K\). The linearisation about either equilibrium reads, \[\begin{equation} \label{rd-US-eqn-US-linear-US-logistic} \mathchoice{\frac{\partial u_1}{\partial t}}{\partial u_1/\partial t}{\partial u_1/\partial t}{\partial u_1/\partial t} = D\mathchoice{\frac{\partial^2 u_1}{\partial x^2}}{\partial^2 u_1/\partial x^2}{\partial^2 u_1/\partial x^2}{\partial^2 u_1/\partial x^2} + r\left (1-\frac{2u_0}{K} \right )u_1, \quad x \in [0,L], \,\, t > 0. \end{equation}\] Now the key insight is that [rd-US-eqn-US-linear-US-logistic] is a constant-coefficient linear PDE, and it can be solved exactly as we did for the diffusion equation last term, with some additional bookkeeping around the extra constant-coefficient term at the end. Alternatively, we can use the ansatz \[\begin{equation} \label{ansatz} u_1 = A\exp(\lambda_k t)\cos\left (\frac{xk\pi}{L}\right ),\quad k=0,1,2,\dots , \end{equation}\] with \(A\) an (unimportant but nonzero) constant. This ansatz is exactly like expanding solutions of the ODEs [gen-US-lin-US-ODE] in terms of \(\mathbfit{u} = \exp(\lambda_k t)\mathbfit{w}_k\).

Substituting this ansatz into [rd-US-eqn-US-linear-US-logistic] and cancelling the function out, we can directly read off the growth rates: \[\begin{equation} \label{rd-US-logistic-US-linstab} \lambda_k = -D\left (\frac{k\pi}{L} \right)^2 + r \left (1-\frac{2u_0}{K} \right ). \end{equation}\] The spatially homogeneous equation is easy enough to analyse with respect to linear stability analysis of the same two equilibria (and in fact we get these growth rates for free by setting \(D=0\) in the above). For this reason, in the case of no-flux boundary conditions, we see that \(\lambda_0\) is the growth rate associated to the spatially homogeneous model. We also have that \(\lambda_0 > \lambda_1 > \lambda_2 > \dots\), so that the effect of diffusion is to stabilise both equilibria to perturbations. The steady state \(u_0=K\) is therefore stable to all perturbations, whereas \(u_0=0\) is unstable for \(k=0\), but for sufficiently large \(k\), it becomes stable. That is, for \(k>k_0\) where \(k_0 = \sqrt{\frac{r}{D}}\frac{L}{\pi}\), these perturbations decay. In practice, stability to such perturbations does not change the fact that the steady state is stable to larger wavelength ones – a small population in a spatially extended environment will still grow away from the extinction steady state.

Of course, our ansatz above is a single term of the solution (which we will call a mode). As in the ODE case, we can justify considering a single mode at a time via expanding \(u_1(x,t) = \sum_{k=0}^\infty A_k(t)\exp(\lambda_k t)\cos(xk\pi/L)\). The key is that these cosine functions are orthogonal to one another for different \(k\), so that the manipulations done in [lin-US-ODE-US-expansion] effectively follow through (subject to details about defining an appropriate inner-product, and technical details about convergence of infinite sums). The take-home message is that considering a single mode at a time is enough to determine the growth or decay of perturbations to a steady state, just as in the ODE case considering a single eigenvector at a time is sufficient. Note that we do not pose initial conditions for this linear equation, but just assume that all modes have some nonzero perturbation – this will be true in general for an arbitrary small perturbation of a spatially constant equilibrium, but is worth spending some time thinking about.

10.1.2 Spatially homogeneous equilibria – survival in a harsh environment

What if we consider the same PDE [rd-US-eqn-US-logistic] but with homogeneous Dirichlet conditions ([rd-US-bcs] with \(b_0=b_L=0\))? Now the steady state \(u_0=K\) of the ODE model does not satisfy these boundary conditions, and so is not a steady state of the PDE, though the steady state \(u_0=0\) is a solution to the PDE (and a steady state, as it will not evolve in time). What changes with the linear stability of \(u_0=0\)?

The linear equation [rd-US-eqn-US-linear-US-logistic] will be the same, although we cannot use the ansatz involving cosines (they will not satisfy the boundary conditions). We can either use separation of variables to solve this linear PDE, or we can guess what will happen to the spatial part of our solutions if we have solved the diffusion equation with Dirichlet conditions before. Essentially we change our cosines to sines to get, \[u_1 = C\exp(\lambda_k t)\sin\left (\frac{xk\pi}{L}\right ),\quad k=1,2,\dots,\] and our growth rates are hence the same as in [rd-US-logistic-US-linstab], except that we no longer have a mode for \(k=0\) (as in this case our ansatz is trivially zero, and cannot grow or decay).

What does this tell us? As before, we have that the growth rates decrease with increasing \(k\), that is \(\lambda_1 > \lambda_2 > \dots\). In the spatially homogeneous model, the growth rate of a perturbation to the extinction steady state was \(\lambda_0=r>0\) (remember this can be seen from [rd-US-logistic-US-linstab] with \(D=0\) and \(u_0=0\)). But now our largest growth rate, \(\lambda_1\), need not be positive. In particular \[\begin{equation} \label{harsh-US-environment-US-growth}\lambda_1 = -D\left (\frac{\pi}{L} \right)^2 + r \implies \lambda_1<0 \quad \textrm{if} \quad L < L_c = \pi\sqrt{\frac{D}{r}}. \end{equation}\] Hence, on a sufficiently small domain, diffusion actually stabilises the extinction steady state, and a small population will die off. Effectively the Dirichlet boundary conditions here are modelling a harsh environment, and random movement of individuals is sufficient to effectively reduce the carrying capacity of the population to zero. If the movement is faster (larger \(D\)), then this occurs on a larger domain \(L\).

What happens when the parameters are such that extinction is no longer possible? Unlike in the ODE model of logistic growth, there is no spatially homogeneous steady state, and so the population density must go to some other state. In principle this could be a time-dependent state, but it turns out this is not possible for this kind of scalar reaction–diffusion equation. Instead, the population will approach another steady state (which can only be found approximately, e.g. numerically), corresponding to a function \(u_0(x)\) solving [spatial-US-steady-US-state]. We will not go into the details of how to approximate this steady state, but give a sketch of how it might depend on the domain length \(L\) in 3.1.

10.1.3 Spectral analysis of the Laplacian

A substantial amount of mathematical theory has been developed around properties of the Laplacian operator. One key piece of this theory which we will make use of is the study of the eigenvalues and eigenfunctions of this operator. These satisfy \[-\nabla^2 w_\mathbfit{k}(\mathbfit{x}) = \rho_\mathbfit{k} w_\mathbfit{k}(\mathbfit{x}),\] on a ‘nice’ domain (typically smooth boundary and simply connected) \(\mathbfit{x} \in \Omega \in \mathbb{R}^n\), \(n=1,2\), or \(3\). The minus sign is conventional so that, with suitable boundary conditions, the eigenvalues \(\rho_\mathbfit{k}\) can be ordered in a strictly increasing sequence tending to \(+\infty\). Though different authors may use different conventions about this, as well as the notation of the eigenvalues \(\rho_\mathbfit{k}\) as opposed to the indexing variable \(\mathbfit{k}\).

For \(n=1\), in the case of no-flux boundary conditions on a domain of length \(L\), we have \[w_k = \cos \left( \frac{x k \pi}{L} \right), \quad \rho_k = \left( \frac{k \pi}{L} \right)^2, \quad k=0,1,2,\dots,\] and for homogeneous Dirichlet conditions we have, \[w_k = \sin \left( \frac{x k \pi}{L} \right), \quad \rho_k = \left( \frac{k \pi}{L} \right)^2, \quad k=1,2,\dots,\] as discussed in the last subsection.

If we consider now \(n=2\), say for homogeneous Dirichlet conditions on the box domain \(\mathbfit{x} = (x,y) \in [0,L_1]\times[0,L_2]\), we can find via separation of variables, \[w_\mathbfit{k} = \sin \left( \frac{n x \pi}{L_1} \right)\sin \left( \frac{m y \pi}{L_2} \right), \quad \rho_\mathbfit{k} = \left( \frac{n \pi}{L_1} \right)^2+\left( \frac{m \pi}{L_2} \right)^2, \quad n=1,2,\dots, \,\, m=1,2,\dots,\] where \(\mathbfit{k}=(n,m)\) is our indexing set. Similarly we can find solutions when our box has Neumann conditions, or even when the vertical and horizontal boundaries have different conditions. We will review other cases of these kinds of functions as needed, but hopefully you will have already seen examples such as Bessel functions.

We will make use of these solutions noting that we can immediately construct ansatzes like [ansatz] to solve linear PDE involving the Laplacian. In some sense these eigenfunctions immediately diagonalise this operator, in the same way that eigenfunctions of a symmetric matrix allow us to solve one component at a time (hence justifying the exponential in \(\lambda\) time dependence).

10.2.1 Linear patterns (Fourier series)

Before jumping into Turing’s theory in detail, we first look at some aspects of linear PDE which will be relevant. We can generalise linear stability theory above for scalar PDE by exploiting the ansatz given in [ansatz]. In some sense, using the eigenfunctions of the Laplacian (cosines for no-flux boundaries, sines for homogeneous Dirichlet conditions) diagonalises the partial differential operator, in the same way that using eigenvectors of a symmetric Jacobian ‘trivially’ diagonalises it. That means that partial differential equations of the form \[\begin{align} \label{linsys} \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} = D_1\nabla^2u + B_1 u + B_2v,\\ \mathchoice{\frac{\partial v}{\partial t}}{\partial v/\partial t}{\partial v/\partial t}{\partial v/\partial t} = D_2\nabla^2v + B_3 u + B_4 v, \end{align}\] i.e. linear equations, have solutions in the form \[u = A\mathrm{e}^{\mathrm{i}\mathbfit{k} \cdot \mathbfit{x} + \lambda t}, \quad v = B\mathrm{e}^{\mathrm{i}\mathbfit{k} \cdot \mathbfit{x} + \lambda t},\] where again we are only expanding a single mode and implicitly using orthogonality to derive an equation for one mode at a time. We have also replaced the cosine by a complex exponential, to account for other possible boundary conditions¹⁸.

Substituting these solutions in we have the algebraic equations: \[\lambda\left(\begin{array}{c}u\\v\end{array}\right)= \begin{pmatrix}-D_1\mathbfit{k}\cdot\mathbfit{k}+ B_1 & B_2 \\ B_3 & -D_2\mathbfit{k}\cdot\mathbfit{k}+ B_4 \end{pmatrix}\left(\begin{array}{c}u\\v\end{array}\right) = \mathsfbfit{M}\left(\begin{array}{c}u\\v\end{array}\right).\] The fundamental theorem of linear algebra states we can only have non-zero \((u,v)\) if \(\det(\mathsfbfit{M}-\lambda \mathsfbfit{I})=0\). Another way to see this is that \(\lambda\) is an eigenvalue of \(\mathsfbfit{M}\). This determinant condition will gives us a quadratic equation in \(\lambda\), which we solve to obtain \(\lambda\) as a function of \(k\) (the so-called dispersion relation). We will derive this again in the next chapter, and spend a lot of time understanding it. As in the scalar case, these growth rates depend eigenvalues of the Laplacian, which can depend heavily on the boundary conditions and geometry of the spatial domain.

For example we might have considered a 2D Cartesian domain with coordinates (\(x_1,x_2\)) and domain lengths \(L_1\) and \(L_2\) and homogeneous Dirichlet boundary conditions \(u(0,x_2)=u(L_1,x_2)= u(x_1,0)=u(x_1,L_2)=0\) and similar for \(v\), so that we would require \[k_1L_1 = n\pi,\quad k_2 L_2 = m\pi,\] so only the part of the solution \(u = A\mathrm{e}^{\mathrm{i}\mathbfit{k} \cdot\mathbfit{x} + \lambda t}\) which can satisfy the boundary conditions is \[u = A_c\mathrm{e}^{\lambda t}\sin\left(\frac{n\pi}{L_1}x_1\right)\sin\left(\frac{m\pi}{L_2} x_2\right).\] Note that the notation here is different, partly to show you alternative ways of viewing these spatial eigenfunctions (and because, e.g., Murray uses this notation with \(\mathbfit{k}\)). We can note that \(\mathbfit{k}\cdot \mathbfit{k} = |\mathbfit{k}|^2=\rho_k\) is how the eigenvalues in this notation relate to what we had before.

This allows for a countable infinity of possible solutions based on choices \(n\) and \(m\). A general solution to this system for the morphogen \(u\) would then take the form \[u(\mathbfit{x},t) = \sum_{n=0}^{\infty}\sum_{m=0}^{\infty}A_{nm}\mathrm{e}^{\lambda(n,m)t}\sin\left(\frac{n\pi }{L_1}x_1\right)\sin\left(\frac{m\pi }{L_2}x_2\right),\] where we recognise \(\lambda\) would now be a function of \(n\) and \(m\). The initial condition \(u(x_1,x_2,0) = u_0\) can be used to define the coefficients \(A_{ij}\) through¹⁹ \[A_{nm} = \frac{\displaystyle \int_{0}^{L_1}\int_{0}^{L_2}u(\mathbfit{x},0)\sin\left(\frac{n\pi }{L_1}x_1\right)\sin\left(\frac{m\pi }{L_2}x_2\right){\mathrm d}{x_1}{\mathrm d}{x_2}} {\displaystyle \int_{0}^{L_1}\int_{0}^{L_2}\sin\left(\frac{n\pi }{L_1}x_1\right)^2\sin\left(\frac{m\pi }{L_2}x_2\right)^2{\mathrm d}{x_1}{\mathrm d}{x_2}}.\] As before, we assume ‘generic’ perturbations such that \(A_{mn}\) are all nonzero, though some may be very small for any given initial perturbation. We also know that if \(\lambda(n,m)<0\) the mode \((n,m)\) will decay to zero exponentially in time. So pretty quickly for some \(t>0\) only terms for which \(\lambda(n,m)\geq0\) will remain visible, and these will grow rapidly as time evolves.

Let’s say for argument’s sake that only the \(n=10\) and \(m=10\) mode has \(\lambda(n,m)\geq0\). Then at some time \(t\gg 0\) the solution can be approximated by \[u(\mathbfit{x},t) \approx A_{10\,10}\mathrm{e}^{\lambda(10,10)t}\sin\left(\frac{10\pi }{L_1}x_1\right)\sin\left(\frac{10\pi }{L_2}x_2 \right).\] Let’s say further that if \(u>0\) the morphogen triggers cells which radiate white to grow (if \(u<0\) the cells remain dark). The pattern produced by this function is a checkerboard pattern, shown in 3.3(a). If we are a little more demanding and require \(u>0.5\) then this becomes a regular spotted pattern in 3.3(b) (assuming \(t\) is not too large). If instead we consider only the mode \(n=1\), \(m=10\), we get stripes for \(u>0\) – see 3.3(c) – and if \(u >0.5\) these stripes become restricted ‘fin’ type patterns. Of course many of these details depend on exactly when we decide to ‘threshold’ the morphogen, and at what level. The main point is that since solutions to linear PDE grow exponentially in time, we can make decent approximations by considering only those modes which grow fastest.

In general a Fourier series can produce any reasonable scalar function \(f(x_1,x_2)\) but these spotted and striped patterns seem so common in nature. Turing asked why? As we shall see it is typical of reaction–diffusion equations that only a small number of nodes, \((n,m)\), values are unstable and grow.

10.2.2 Nonlinearity

Of course if \(\lambda(n,m)>0\) then in theory \(u\) can grow without bound. However, most realistic models are nonlinear and permit spatially varying equilibria with bounded values, as we shall see this includes spotted and striped type patterns. Of course we can always linearise nonlinear systems and look for which modes might grow. This is the essence of Turing’s analysis which tells us something about the set of patterns which might form in reaction–diffusion systems. We expect a restricted space of unstable modes \((n,m)\) which can form spotted/striped patterns; the nonlinear effects then tend to take over and stabilise the growth of the pattern until it relaxes. The question of how much the pattern changes is one which we will come back to throughout the term.

10.2.3 Domain shape

The above discussion assumed we were in a nice rectangular or square domain. I have yet to come across a square animal! The shape of the domain, as well as the boundary conditions, can play a vital role in the allowed patterns the system forms. For example, one issue we cover this term is pattern formation in a circular domain. We shall see later that the solutions to a system like [linsys], when our domain is circular (in a polar coordinate system \((r,\theta)\)), take the form \[u(r, \theta,t) = \sum_{n=1}^{\infty}\sum_{k_m^n}A_{mn}J_n(k_m^n r)\cos(n\theta) + B_{mn}J_n(k_m^n r)\sin(n\theta),\] where \(J_n(k_m^n r)\) is the \(n^\text{th}\) Bessel function. Bessel functions vary sinusoidally, but unlike sine and cosine the amplitude of variation decays away with \(r\). We see in 3.4(c) this tends to produce patterns which form on rings surrounding the centre of the domain (the Bessel function amplitudes decide where). With stripes appearing as circles and spots lying on radii around a fixed central spot. In (d) we show an example of a pattern formed on an annular domain. It is an almost periodic set of ‘stripes’. We will see later that this pattern can plausibly explain branching type behaviour; in fact reaction diffusion mechanisms have been used to explain the splitting of branches we see in trees, or the bifurcations of blood vessels throughout our bodies.

11.2.1 Interlude: boundary conditions and eigenfunctions

Next we consider, separately, the cases where \(\operatorname{Re}(\lambda_0)<0\) (so that the spatially homogeneous equilibrium is stable in the absence of diffusion), and \(\operatorname{Re}(\lambda_\mathbfit{k})>0\) for some wave number \(\mathbfit{k}\). This will lead to a set of four conditions on the parameters and the geometry in order to ensure this kind of instability occurs.

11.2.2 Growing inhomogeneous modes

Next we seek for inhomogeneous unstable modes \(\mathbfit{k} = \left(n_1\pi/L_1,n_2\pi/L_2,\dots n_n\pi/L_n\right)\) which are unstable, so they don’t vanish and instead grow (at least one \(\lambda_\mathbfit{k}\) has positive real part). Since \((F_u+G_v)<0\) the trace \((F_u +G_v)-\rho_k(D+1)\) is less than zero. Thus we can only have positive \(\lambda_\mathbfit{k}\) if \(\det(\mathsfbfit{M})<0\) so that one of the roots must be positive, that is to say \[\begin{equation} \label{h-US-eq} h(\rho_k) = D\rho_k^2 - \rho_k(G_v+ D F_u)+ (F_uG_v -G_uF_v) < 0. \end{equation}\] Since \(D\rho_k^2>0\) and \((F_uG_v -G_uF_v)>0\) (so the zeroth order mode vanishes), this can only occur if \[G_v+ D F_u>0.\]

This condition is necessary but not sufficient to guarantee that \(\det(\mathsfbfit{M})<0\) for some \((\rho_k\neq 0)\). To help understand this, look at panel (a) of 4.3, ignoring for the moment the parameters \(d\), \(d_c\) etc (and noting that this author has used \(k^2\) in place of our \(\rho_k\)). We need the minimum value of \(h\) to be negative to ensure an instability. To compute this condition, we will make an important approximation – namely that we can treat \(h(\rho_k)\) as a function of a continuous variable, \(\rho_k\), despite the fact that really \(\rho_k\) can only take a discrete set of values. For domains which are sufficiently large, this approximation is OK, as the modes will scale like \(1/L^2\) (see, for example [mode-US-numbers]).

With this approximation in mind, and since we know that \(h\) has a single global minimum, we compute the wave number at this minimum as: \[h'(\rho_k^*) = 2D\rho_k^* - (G_v+ D F_u)=0 \implies \rho_k^ = \frac{G_v+ D F_u}{2D},\] and so we have that the minimum value of \(h\) is: \[h(\rho_k^*) = \frac{(G_v+ D F_u)^2}{4D} - \frac{(G_v+ D F_u)^2}{2D}+ (F_uG_v -G_uF_v) = -\frac{(G_v+ D F_u)^2}{4D} + (F_uG_v -G_uF_v).\] Rearranging this, our condition for instability (that \(h(\rho_k^*)<0\)) is equivalent to: \[(G_v+ D F_u)^2 - 4D(F_uG_v -G_uF_v) > 0.\]

11.3.1 The Turing instability requires short-range activation, long-range inhibition

The first and third condition imply that \(F_u\) and \(G_v\) must be of opposing sign. To see this we note that if they had the same sign then \(F_u + G_v<0\) tells us they would have to be negative, but then the third condition would say \[D \vert F_u \vert + \vert G_v\vert < 0 \implies D <-\frac{\vert G_v\vert}{\vert F_u\vert}.\] But physically we need \(D>0\) (negative diffusion does not make sense) so this cannot be permissible. Therefore, if \(F_u>0\) we must have \(G_v<0\) and vice versa.

One can see that this condition (of opposing \(F_u\) and \(G_v\)) make physical sense. The ‘reaction’ \(F\) promotes \(u\) and \(G\) promotes \(v\). Consider for example the case \(F_u>0\) and \(G_v<0\). To induce a separation in concentration peaks, we want it to be the case that where a small increase in \(u\) leads to an acceleration in the production of \(u\) to coincide with small increases in \(v\) leading to a drop in \(v\) so that (in this case) \(u\) can become dominant over \(v\). Decreases in \(u\) and \(v\) should lead to \(v\) being favoured. If both \(F_u\) and \(G_v\) are the same sign they will both tend to grow and decay simultaneously, this will not promote separation and hence pattern formation.

The biochemical terminology here is that the \(u\) species is a self-activator (assuming without loss of generality that \(F_u>0\)), and that \(v\) is a self-inhibitor. Often this is shortened to ‘activator and inhibitor’ dynamics, though it is important to note that we have only determined how these species interact with themselves (and only at the steady state \((u_0,v_0)\)).

What does this requirement say about the permissible values of \(D\)? Let’s fix \(F_u>0\) so that \(u\) is the activator, so we must have \(G_v<0\). We then have from the first and the third inequalities in [turingconditions]: \[F_u + G_v<0 \iff F_u < -G_v \implies \vert F_u \vert < \vert G_v \vert \implies \frac{\vert G_v\vert }{\vert F_u \vert} >1,\] and \[DF_u+G_v > 0 \implies D\vert F_u\vert +\vert G_v\vert >0 \implies D>\frac{\vert G_v\vert }{\vert F_u \vert} >1.\] A similar argument shows that if \(G_v>0\), we must have \(F_u<0\) and \(D<1\). Note that these inequalities are strict, i.e. \(D\) cannot equal \(1\). This is why the Turing instability is often referred to as a diffusion-driven instability, as it requires that the diffusion rates of the two species are different. We can also see that the inhibitor in either case must always diffuse more quickly than the activator. This is referred to as short-range activation, long-range inhibition, or sometimes as LALI: local activation, lateral inhibition.

11.3.2 Restrictions on sign structures of the Jacobian

What about the other terms in the Jacobian, \(F_v\) and \(G_u\)? Well since \(F_u\) and \(G_v\) have opposite signs, we must have \(F_uG_v<0\). But the second condition means that we must have \(F_uG_v -G_uF_v> 0\), so this inequality can never be satisfied unless \(F_v\) and \(G_u\) have opposite signs. Therefore, the Jacobian matrix must take one of the following forms: \[\mathsfbfit{J} = \begin{pmatrix} + & -\\ + & - \end{pmatrix},\quad \begin{pmatrix} + & +\\ - & - \end{pmatrix},\quad \begin{pmatrix} - & -\\ + & + \end{pmatrix},\quad \begin{pmatrix} - & +\\ - & + \end{pmatrix}.\]

Of course, one can always just consider the first two as the latter two are just a relabelling of \(u\) and \(v\) in a sense, whereas the first two have qualitatively different behaviours. In particular, you can show that the coefficient of the linear term of \(u_1\) and \(v_1\) will be in-phase or out-of-phase, depending on which of these two forms it takes. While \(u_1\) and \(v_1\) will have different forms (and especially the functions satisfying the nonlinear reaction–diffusion equations, \(u\) and \(v\), will generally appear very different from one another), the maxima and minima of their concentrations will coincide or be exactly out-of-phase with one another.

To see this in and out of phase behaviour, consider the second equation in [turinglin]: \[\lambda_k v_1 = G_u u_1 + G_v v_1 - \rho_k D v_1.\] Let’s also assume we have \(F_u>0\) and \(G_v<0\) (so that \(D>1\) for a pattern-forming instability). Moving the \(v_1\) terms to the left side, we have: \[(\lambda_\mathbfit{k} - G_v + D \rho_k)v_1 = G_u u_1.\] For an instability, we must have \(\lambda_\mathbfit{k}>0\), and by assumption the other two terms on the left are positive, so \(v_1\) and \(u_1\) are related in sign by \(G_u\). That is, if \(G_u>0\), then \(u_1\) and \(v_1\) must be of the same sign, and hence their spatial oscillations must be in-phase. If \(G_u<0\), they must be out-of-phase, as there is now a sign difference between them, so anywhere \(u_1\) is positive, \(v_1\) must be negative and vice-versa. See 4.4.

This has useful biological applications, as often an experimentalist has identified a candidate activator and inhibitor, and knows if they have peaks (that is, regions of high concentration) in the same place with one another, or in different places. This can help determine the correct form for a model of their system, as often we do not understand the molecular biology well enough to write down what forms \(F\) and \(G\) should take. Similar remarks can be made about ecological applications, particularly predator–prey interactions.

See chapters 2 and 3 of the second volume of Murray’s Mathematical Biology for a much more thorough discussion of the biological interpretation of these Turing conditions.

11.5.1 The peak unstable mode

For a fixed wavenumber, Equation [turinglambda] gives the growth rate \(\lambda(\rho_k)\). If we assume that the domain is sufficiently large, we can compute a value of \(\rho_k\) for which \(\lambda(\rho_k)\) is maximal. This will then be the mode which grows the fastest, at least as long as the linear approximation is valid. To find this we differentiate \(\lambda(\rho_k)\) with respect to \(\rho_k\) and solve for \(\mathchoice{\frac{{\mathrm d}\lambda}{{\mathrm d}\rho_k}}{{\mathrm d}\lambda/{\mathrm d}\rho_k}{{\mathrm d}\lambda/{\mathrm d}\rho_k}{{\mathrm d}\lambda/{\mathrm d}\rho_k}=0\). Since \(\lambda(\rho_k)\) takes the form \[\begin{aligned} &\frac{1}{2}\left[(F_u +G_v) -\rho_k(D+1) + \left\{\left((F_u +G_v) -\rho_k(D+1) \right)^2- 4\det(\mathsfbfit{M})\right\}^{1/2}\right],\\ &\det(\mathsfbfit{M}) = D\rho_k^2 - \rho_k(DF_u +G_v) + (F_uG_v -F_vG_u). \end{aligned}\] This is not a trivial task! With no little effort solving \(\mathchoice{\frac{{\mathrm d}\lambda}{{\mathrm d}\rho_k}}{{\mathrm d}\lambda/{\mathrm d}\rho_k}{{\mathrm d}\lambda/{\mathrm d}\rho_k}{{\mathrm d}\lambda/{\mathrm d}\rho_k}=0\) gives \[\begin{equation} \label{kmaximal} \rho_k^* = \frac{1}{D-1}\left[(D+1)\left(\frac{-F_v G_u}{D}\right)^{1/2}-F_u+G_v\right], \end{equation}\] as the fastest growing mode. Alongside bounding the range of unstable modes, this computation often matches observed patterns from full numerical simulations reasonably well, and so despite the algebraic difficulty, it can be a useful heuristic for predicting the wavelength between pattern elements (e.g. the spacing between spots) given model parameters.

11.6.1 Summary of the conditions

We can summarise the conditions by saying that for the homogeneous steady state to be stable, we simply require, \[\begin{equation} \label{hom-US-ss-US-stab-US-schnack} (a+b)^3>b-a, \end{equation}\] and that for a diffusion-driven instability, we also need, \[(a+b)^3<D(b-a), \quad \textrm{and} \quad (D(b-a)-(a+b)^3)^2>4D(a+b)^4.\] As noted, the requirement that \(b>a\) for \(F_u\) to be positive can be seen in the first of these two conditions. We also see that, as long as \(b-a>0\), both of these conditions will be automatically satisfied in the limit of \(D \to \infty\). We can confirm this by plotting the conditions numerically.

In 4.5, we plot a region in \((a,b)\) parameter space for different values of \(D\). The black region corresponds to values of \(a\) and \(b\) where the Turing conditions are satisfied, and so for sufficiently large domains we suspect those parameters to give rise to patterned solutions. The grey regions are where the homogeneous equilibrium remains stable, and the white region is where the homogeneous equilibrium is unstable to homogeneous (\(\rho_0=0\)) perturbations, or equivalently where condition [hom-US-ss-US-stab-US-schnack] is violated. The blue dashed line is where \(a=b\). With some work, you can actually draw these boundary curves by hand, though this is much simpler to do using modern numerical software as I have done here.

As predicted, as \(D\) increases the entire region with a stable homogeneous equilibrium in the absence of diffusion (the grey region) becomes Turing unstable as long as \(b>a\), as given by the blue line. On the other hand, for \(D=5\) we see a relatively small Turing unstable region which shrinks rapidly and disappears for \(D<2.5\) or so. Again this has important implications in terms of requiring a massive difference in the diffusion rates of our two morphogens in order to realise a ‘large’ Turing space. This is a major obstacle of robustness of this theory, as we would really want a large and robust parameter space where we see these instabilities, which is not very sensitive to small random fluctuations.

11.6.2 Mode selection via parameter controls (non-examinable)

We can also consider a means by which we exert control on the patterns formed by this system, as long as we can control the kinetic parameters (and thereby the Jacobian). Such control approaches have been used in chemical applications of Turing systems, as well as in more modern synthetic biology applications.

The Turing conditions tell us that when \[(DF_u +G_v)^2 = 4D(F_uG_v-F_vG_u),\] the system is on the border of allowing turning instabilities. This occurs when \[\rho_k= \frac{D F_u + G_v}{2D}.\] So only one mode can potentially be unstable if we are just passed this threshold. One could for example demand this occur for an \(n_1=5,n_2=5\) mode in a 2-D system (creating a checkerboard pattern like in 3.3(a)). In the example we are looking at \[\rho_k= \frac{D(b-a)-(a+b)^3}{2D(a+b)} = 25\pi^2\left(\frac{1}{L_1^2} +\frac{1}{L_2^2} \right),\] which would, in principle, allow us to pick parameters such that the system is unstable exactly at that mode. In fact the amplitude of the mode can, to some degree, also be controlled, though this is beyond our scope here. In fact recent efforts have allowed for the design of elaborate Turing patterns controlled via spatial heterogeneity and nonlinear spot and stripe selection mechanisms, though again you will not be expected to be familiar with these. See this paper for further discussion. An example is 4.6.

12.2.1 Homogeneous equilibrium & linear stability analysis

Despite the geometry being different, this problem is otherwise the same as that in 3.1.2. While there are two spatially homogeneous equilibria in the absence of the diffusion term, given by \(u_0=0\) or \(1\), the latter does not satisfy the boundary condition at \(r=L_r\). So we can only consider instabilities around \(u_0=0\). In the absence of diffusion this steady state is unstable, though this may not be the case when diffusion is considered.

The linearisation of the equation is the same as we’ve seen before. Writing \(u = u_0 + \varepsilon u_1 = \varepsilon u_1\) we find, \[\begin{equation} \label{Petri-US-Linear} \mathchoice{\frac{\partial u_1}{\partial t}}{\partial u_1/\partial t}{\partial u_1/\partial t}{\partial u_1/\partial t} = D\nabla^2 u_1 + u_1(1-2u_0) = D\nabla^2 u_1 + u_1. \end{equation}\] As before, we substitute in \(u_1 = \exp(\lambda_k t)w_k(r,\theta)\) to find that the perturbation’s growth rates are given by, \[\begin{equation} \label{Disp-US-Bessel-US-Scalar} \lambda_k u_1 = -\rho_k D u_1 + u_1\implies \lambda_k = -\rho_k D + 1, \end{equation}\] where the spatial eigenvalue \(\rho_k\) is determined from the equation, \[\begin{equation} \label{eig-US-petri} \nabla^2w_k + \rho_k w_k = \frac{1}{r}\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}\left(r \mathchoice{\frac{\partial w_k}{\partial r}}{\partial w_k/\partial r}{\partial w_k/\partial r}{\partial w_k/\partial r}\right) + \frac{1}{r^2}\mathchoice{\frac{\partial^2 w_k}{\partial\theta^2}}{\partial^2 w_k/\partial\theta^2}{\partial^2 w_k/\partial\theta^2}{\partial^2 w_k/\partial\theta^2} + \rho_k w_k = 0. \end{equation}\]

12.2.2 Solving the spatial eigenvalue problem

To determine the growth rates, we have to solve the Laplacian eigenvalue problem to find the spatial eigenvalues \(\rho_k\). We assume a separable solution of the form \(w_k(r,\theta) = R(r)P(\theta)\). Substituting this into [eig-US-petri] we find, \[\frac{1}{r}\mathchoice{\frac{\partial}{\partial r}}{\partial/\partial r}{\partial/\partial r}{\partial/\partial r}\left(r \mathchoice{\frac{\partial(RP)}{\partial r}}{\partial(RP)/\partial r}{\partial(RP)/\partial r}{\partial(RP)/\partial r}\right) + \frac{1}{r^2}\mathchoice{\frac{\partial^2 (RP)}{\partial\theta^2}}{\partial^2 (RP)/\partial\theta^2}{\partial^2 (RP)/\partial\theta^2}{\partial^2 (RP)/\partial\theta^2} + \rho_k RP = P\frac{1}{r}\mathchoice{\frac{{\mathrm d}}{{\mathrm d}r}}{{\mathrm d}/{\mathrm d}r}{{\mathrm d}/{\mathrm d}r}{{\mathrm d}/{\mathrm d}r}\left(r \mathchoice{\frac{{\mathrm d}R}{{\mathrm d}r}}{{\mathrm d}R/{\mathrm d}r}{{\mathrm d}R/{\mathrm d}r}{{\mathrm d}R/{\mathrm d}r}\right) + R\frac{1}{r^2}\mathchoice{\frac{{\mathrm d}^2 P}{{\mathrm d}\theta^2}}{{\mathrm d}^2 P/{\mathrm d}\theta^2}{{\mathrm d}^2 P/{\mathrm d}\theta^2}{{\mathrm d}^2 P/{\mathrm d}\theta^2} + \rho_k RP= 0.\] Multiplying both sides of this equation by \(r^2/RP\) we have, \[\frac{r}{R}\mathchoice{\frac{{\mathrm d}}{{\mathrm d}r}}{{\mathrm d}/{\mathrm d}r}{{\mathrm d}/{\mathrm d}r}{{\mathrm d}/{\mathrm d}r}\left(r \mathchoice{\frac{{\mathrm d}R}{{\mathrm d}r}}{{\mathrm d}R/{\mathrm d}r}{{\mathrm d}R/{\mathrm d}r}{{\mathrm d}R/{\mathrm d}r}\right) + \frac{1}{P}\mathchoice{\frac{{\mathrm d}^2 P}{{\mathrm d}\theta^2}}{{\mathrm d}^2 P/{\mathrm d}\theta^2}{{\mathrm d}^2 P/{\mathrm d}\theta^2}{{\mathrm d}^2 P/{\mathrm d}\theta^2} + \rho_kr^2 = 0 \implies \frac{r}{R}\mathchoice{\frac{{\mathrm d}}{{\mathrm d}r}}{{\mathrm d}/{\mathrm d}r}{{\mathrm d}/{\mathrm d}r}{{\mathrm d}/{\mathrm d}r}\left(r \mathchoice{\frac{{\mathrm d}R}{{\mathrm d}r}}{{\mathrm d}R/{\mathrm d}r}{{\mathrm d}R/{\mathrm d}r}{{\mathrm d}R/{\mathrm d}r}\right) + \rho_kr^2 = - \frac{1}{P}\mathchoice{\frac{{\mathrm d}^2 P}{{\mathrm d}\theta^2}}{{\mathrm d}^2 P/{\mathrm d}\theta^2}{{\mathrm d}^2 P/{\mathrm d}\theta^2}{{\mathrm d}^2 P/{\mathrm d}\theta^2} = C_\theta,\] where we have separated the \(r\) and \(\theta\) functions, and hence must have that both sides of this equation are constant, which we set to be equal to \(C_\theta\).

We then have the equations, \[\begin{align} \label{R-US-eq} -\frac{1}{Rr}\mathchoice{\frac{{\mathrm d}}{{\mathrm d}r}}{{\mathrm d}/{\mathrm d}r}{{\mathrm d}/{\mathrm d}r}{{\mathrm d}/{\mathrm d}r}\left(r \mathchoice{\frac{{\mathrm d}R}{{\mathrm d}r}}{{\mathrm d}R/{\mathrm d}r}{{\mathrm d}R/{\mathrm d}r}{{\mathrm d}R/{\mathrm d}r}\right) & = \rho_k-\frac{C_\theta}{r^2},\\ - \frac{1}{P}\mathchoice{\frac{{\mathrm d}^2 P}{{\mathrm d}\theta^2}}{{\mathrm d}^2 P/{\mathrm d}\theta^2}{{\mathrm d}^2 P/{\mathrm d}\theta^2}{{\mathrm d}^2 P/{\mathrm d}\theta^2} = C_\theta, \end{align}\] and the boundary conditions, \[R(L_r)=0, \quad P(0) = P(2\pi), \quad \mathchoice{\frac{{\mathrm d}P}{{\mathrm d}\theta}}{{\mathrm d}P/{\mathrm d}\theta}{{\mathrm d}P/{\mathrm d}\theta}{{\mathrm d}P/{\mathrm d}\theta}(0) = \mathchoice{\frac{{\mathrm d}P}{{\mathrm d}\theta}}{{\mathrm d}P/{\mathrm d}\theta}{{\mathrm d}P/{\mathrm d}\theta}{{\mathrm d}P/{\mathrm d}\theta}(2\pi).\] We now proceed to solve these two ODEs.

12.2.3 Finding \(P\): periodic boundary conditions

We write the equation for \(P\) as, \[\mathchoice{\frac{{\mathrm d}^2 P}{{\mathrm d}\theta^2}}{{\mathrm d}^2 P/{\mathrm d}\theta^2}{{\mathrm d}^2 P/{\mathrm d}\theta^2}{{\mathrm d}^2 P/{\mathrm d}\theta^2} +C_\theta P= 0.\] We can go through the cases of \(C_\theta\) taking different signs, but again should note that exponentials and linear functions will fail to satisfy the periodic boundary conditions. This can be shown in detail by applying the boundary conditions as in the Neumann case in 5.1. So we must have \(C_\theta \geq 0\), with the solution for \(C_\theta=0\) being a constant.

We then write our solutions as, \[P(\theta) = A\cos\left(\sqrt{C_\theta} \theta\right) + B \sin\left(\sqrt{C_\theta} \theta\right).\] Using the first boundary condition we have: \[A\cos(0) + B \sin(0) = A = A\cos\left(2\sqrt{C_\theta} \pi\right) + B \sin\left(2\sqrt{C_\theta} \pi\right),\] where we need the \(\cos\) term to equal \(1\) and the \(B\sin\) term to go to \(0\). For the former we require \(2\sqrt{C_\theta} \pi=2 n \pi,\) so \(\sqrt{C_\theta}\) must be an integer, that is \(\sqrt{C_\theta} = n\), \(n=0,1,2,3,\dots\).

The second boundary condition reads, \[-A\sqrt{C_\theta}\sin(0) + B \sqrt{C_\theta}\cos(0) = B\sqrt{C_\theta} = A\sqrt{C_\theta}\sin\left(2\sqrt{C_\theta} \pi\right) + B \sqrt{C_\theta}\cos\left(2\sqrt{C_\theta} \pi\right).\] As before, this requires that \(\sqrt{C_\theta} = n\), an integer. So the solution for \(P\) is, \[P(\theta) = A\cos(n \theta) + B \sin(n \theta),\quad n=0,1,2,\dots,\] where we note that \(C_\theta = n^2\). The inclusion of \(n=0\) covers the case that \(C_\theta=0\), which has the constant solution \(P =A\).

12.2.4 Finding \(R\): the Bessel equation

The equation for \(R\) is given by, after expanding the derivative terms and multiplying by \(r^2\), \[\begin{equation} \label{R-US-eq-US-Bessel} r^2\mathchoice{\frac{{\mathrm d}^2 R}{{\mathrm d}r^2}}{{\mathrm d}^2 R/{\mathrm d}r^2}{{\mathrm d}^2 R/{\mathrm d}r^2}{{\mathrm d}^2 R/{\mathrm d}r^2} + r \mathchoice{\frac{{\mathrm d}R}{{\mathrm d}r}}{{\mathrm d}R/{\mathrm d}r}{{\mathrm d}R/{\mathrm d}r}{{\mathrm d}R/{\mathrm d}r} +( \rho_kr^2-n^2)R = 0, \end{equation}\] where the \(n^2\) comes from \(C_\theta\) in [R-US-eq] and hence represents a coupling of the two boundary-value-problems. It turns out that this equation does not admit closed-form solutions involving the usual algebraic and trigonometric functions, and instead has solutions in terms of Bessel functions. These functions are complicated, but we give a quick overview here – do not overly worry about the details as you will not be expected to be an expert on these special functions, but just have some sense of how they look.

If we scale the spatial variable as \(x = \sqrt{\rho_k}r\), we find that [R-US-eq-US-Bessel] transforms as, \[x^2\mathchoice{\frac{{\mathrm d}^2 R}{{\mathrm d}x^2}}{{\mathrm d}^2 R/{\mathrm d}x^2}{{\mathrm d}^2 R/{\mathrm d}x^2}{{\mathrm d}^2 R/{\mathrm d}x^2} + x \mathchoice{\frac{{\mathrm d}R}{{\mathrm d}x}}{{\mathrm d}R/{\mathrm d}x}{{\mathrm d}R/{\mathrm d}x}{{\mathrm d}R/{\mathrm d}x} +( x^2-n^2)R = 0,\] which is Bessel’s equation. It has two linearly independent sets of solutions given by Bessel functions of the first kind \(J_n(x)\) and Bessel functions of the second kind \(Y_n(x)\). For integer \(n\), we can define these functions via integrals, \[\begin{align} J_n(x) &= \frac{1}{\pi}\int_0^\pi \cos(n\tau-x\sin(\tau))\,{\mathrm d}\tau = \frac{1}{2\pi}\int_{-\pi}^\pi \mathrm{e}^{\mathrm{i}(x\sin(\tau)-n\tau)} {\mathrm d}\tau,\\ Y_n(x) &= \frac{1}{\pi}\int_0^\pi\sin(x\sin(\tau)-n\tau)\,{\mathrm d}\tau-\frac{1}{\pi}\int_0^\infty(\mathrm{e}^{n\tau}+(-1)^n \mathrm{e}^{-n\tau})\mathrm{e}^{-x\sinh(\tau)}{\mathrm d}\tau\label{Y-US-eq-US-1} \end{align}\] There are many other ways to write these functions. For example, they can also be written as, \[\begin{align} \label{J-US-Series} J_n(x) &= \sum_{m=0}^\infty \frac{(-1)^m}{m!(m+n)!}\left(\frac{x}{2} \right)^{2m+n}\\ Y_n(x) &= \lim_{\alpha \to n} \frac{J_\alpha(x)\cos(\alpha \pi)-J_{-\alpha}(x)}{\sin(\alpha \pi)}. \end{align}\] We will not go through these functions in any more detail except to note that \(Y_n(x)\) becomes infinite at \(x=0\). You can see this by looking at [Y-US-eq-US-1] for \(x \ll 1\): \[Y_n(x) \approx \frac{1}{\pi}\int_0^\pi\sin(-n\tau)\,{\mathrm d}\tau-\frac{1}{\pi}\int_0^\infty(\mathrm{e}^{n\tau}+(-1)^n \mathrm{e}^{-n\tau})\,{\mathrm d}\tau \approx -\frac{1}{\pi}\int_0^\infty \mathrm{e}^{n\tau}{\mathrm d}\tau \to -\infty,\] where the \(\exp(n\tau)\) term in the second integral diverges. Even for \(n=0\), this term will diverge as the integral is itself unbounded. We show plots of these two functions in 5.3.

Returning to \(R\), and remembering that we have scaled \(x = \sqrt{\rho_k}r\), we have that for each \(n\), solutions to [R-US-eq-US-Bessel] take the form, \[R_n(r) = A_nJ_n(\sqrt{\rho_k}r) + B_nY_n(\sqrt{\rho_k}r).\] As we have mentioned, \(Y_n\) becomes unbounded for \(r=0\). As this point is inside our domain, we must therefore take \(B_n=0\) to have bounded solutions. Finally, we do not yet have the values for \(\rho_k\). These can be obtained by using the \(R\) boundary condition, \[\begin{equation} \label{R-US-Bessel-US-BC} R_n(L_r)=0 \implies J_n(\sqrt{\rho_k}L_r) = 0. \end{equation}\] We see that, for each \(n\) coming from the \(\theta\) part of the problem, the Bessel function will oscillate around zero. In fact it will have infinitely many zeroes, so that for each fixed \(n\), there will be infinitely many values of \(\rho_k\) that satisfy [R-US-Bessel-US-BC]. As these values may be different for each \(n\), we modify the subscript of our spatial eigenvalue to include this dependence and write \(\rho_{k,n}\). Hence the solution to [R-US-eq-US-Bessel] in full generality can be written as, \[R_n(r) = \sum_{k=1}^\infty A_{k,n}J_n(\sqrt{\rho_{k,n}}r).\] We have chosen to index these zeroes starting at \(k=1\), but this choice is somewhat arbitrary, and just indicates that this first zero of \(J_k(x)\) is not in general \(x=0\) itself (though it will be for \(n>0\)). These values of \(\rho_{k,n}\) do not have closed-form expressions, and must be solved for numerically or, in special cases such as \(k \to \infty\), using analytic approximations. We will not pursue this for this course.

12.2.5 Putting it all together

Finally we can write our solution to the linear problem governing perturbations [Petri-US-Linear] in general as, \[u_1(t,r,\theta) = \sum_{n=0}^\infty \sum_{k=1}^\infty\exp\left(\lambda_{k,n} t \right)J_n\left(\sqrt{\rho_{k,n}}r \right)\left(A_{k,n}\cos(n\theta) + B_{k,n}\sin(n\theta) \right),\] where we can update the growth rate relationship [Disp-US-Bessel-US-Scalar] to incorporate these two indices as, \[\lambda_{k,n} = -\rho_{k,n} D + 1.\] So what can we say about stability? Recalling the results from 3.1.2, we found that there was a domain-size dependent switch between growth and extinction of the population with these harsh (Dirichlet) boundary conditions. Do we observe the same thing here?

Yes, in fact we do. Consider the smallest zero of the Bessel function \(J_n(x)\) shown in 5.3(a). This zero occurs for \(n=0\). If we let \(x_0\) be the first such zero, that is \(J_0(x_0)=0\), then we have that \(\sqrt{\rho_{1,0}}L_r = x_0\) determines the value of the smallest spatial eigenvalue. So \(\rho_{1,0} = (x_0/L_r)^2\). We then have by the growth rate given above that if, \[L_r < \sqrt{D}x_0 \implies \lambda_{1,0}<0,\] then the extinction state will be stable, as all other spatial eigenvalues are larger (and hence all other growth rates are smaller). While the notation is messier, this is exactly what happened in 3.1.2, and qualitatively the result is identical – compare this to [harsh-US-environment-US-growth] noting that here we have taken \(r=1\). The main difference is that the critical domain size no longer has a nice analytic expression, and instead corresponds to the first zero of the Bessel function \(J_0(x)\).

12.3.1 Comparison with regeneration rates

Experimental data shown in 5.7(a) indicate that there is a finite range of Calcium concentrations (modelled by the parameter \(b\) in our model) for which full hair growth (regeneration) can occur. We also see that, within this range the hair spacing decreases as the calcium concentration increases, quickly at first then more gradually. Further still the amplitude of the pattern decreases as the concentration approaches both ends of this range.

In order for patterns to form we must enforce all the Turing conditions. We have already established for this system that, for fixed \(D\), there is a range of values of \(b\) which satisfy the Turing conditions. We can, with some effort, parametrise the boundaries given in 4.5, which correspond to values of \((a,b)\) precisely where some growth rate \(\lambda_k=0\), and hence on one side we observe Turing instability, and on the other we observe stability of the spatially homogeneous state (or homogeneous instability, for the white region). We show an example of how this prediction compares for fixed \(a\) and \(D\) values in 5.7.

So in summary the model predicts that, for the right parameters

• There is some finite range \(b\in [b_\text{min},b_\text{max}]\) for which growth can occur.²⁵

• The spacing of individual pattern elements, which become the hairs, scales with the size of the domain.

This model is an early example of a reaction–diffusion pattern mechanism actually being able to explain measurable phenomena, a first step into the goal of giving biologists confidence in the mechanism.

13.1.1 Snake skin pattern formation via chemotaxis

Relevant reading: Murray book I, chapter 11.4; and book II, chapter 4.11.

Snake coats can exhibit a wide variety of spotted and striped patterns, both along the length of the snake and those around the its cross-section. While some species exhibit almost identical markings, snakes within the same family can also exhibit both length wise and width wise stripes, or mixtures of skin patterns; see 6.1 for some examples.

Murray talks at some length about the biology of alligator stripe patterns in chapter 4, noting that there has been much experimental work on this mechanism. As Murray relates in section 4.11, there is some evidence that a similar mechanism could be at work in snake skin patterning. The basic assumptions of such a mechanism are:

1. The pattern is fixed in the dermis, below the outer epidermis which is what we see. This is important as epidermal cells do not move much, whereas those in the dermis are much more motile.

2. Pigment producing cells in this region exhibit phenotypes leading to different coloured regions, and this pigment production is dependent on some chemical (the ‘morphogen’) reaching a critical concentration.

3. The motion and concentrations of these cells proceeds by a chemotactic mechanism.

Given these assumptions, we now write down and analyse a model of chemotaxis-driven pattern formation.

13.1.2 General cross-diffusion systems

Briefly, we note that the above analysis generalises substantially to other kinds of directed transport. We consider the nonlinear cross-diffusion system, \[\begin\{align\} \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} &= \bm{\nabla}\cdot \left (d_1(u,v)\bm{\nabla}u + d_2(u,v) \bm{\nabla}v \right ) + F(u,v), \\ \mathchoice{\frac{\partial v}{\partial t}}{\partial v/\partial t}{\partial v/\partial t}{\partial v/\partial t} &= \bm{\nabla}\cdot \left (d_3(u,v)\bm{\nabla}u + d_4(u,v) \bm{\nabla}v \right ) + G(u,v), \end{align} where the $d_i$ and $F$, $G$ are given sufficiently smooth functions.\footnote{For systems like these, PDE-theory questions of existence and regularity of solutions become quite a bit harder, and blow-up phenomena are much more common. Here we will assume the functions are such that we can neglect these issues, but it is worth noting that one has to choose these nonlinear functions with some care.} We now that \cref{scaledsnake} is a special case of this model with $d_1=D_u$, $d_3=0$, $d_4=D_v$ and $d_2 = -\alpha u$, and $u$ replaced by $n$, $v$ replace by $c$. As in the case of chemotaxis, if we impose Neumann conditions on $u$ and $v$, then this automatically gives no-flux conditions. This can be seen because the flux of $u$, for example, is $\mathbfit{J}_u =d_1(u,v)\bm{\nabla}u + d_2(u,v) \bm{\nabla}v$. The converse, that no-flux implies Neumann conditions, is not in general true. However, it turns out that for all permissible values of $u$ and $v$, we need the matrix, \begin{equation} \mathsfbfit{D} = \begin{pmatrix} d_1(u,v) & d_2(u,v) \\ d_3(u,v) & d_4(u,v) \end{pmatrix},\] to be positive-definite in order to guarantee existence of sufficiently smooth solutions²⁷. One can then show that imposing no-flux conditions must imply Neumann conditions.

Homogeneous equilibria of this system satisfy \(F(u_0,v_0)=G(u_0,v_0)=0\). We can perform the usual perturbation analysis for linear stability by setting \(u = u_0 + \varepsilon u_1\), \(v = v_0 + \varepsilon v_1\). After some careful algebra, we find that the linear perturbations satisfy, \[\begin\{align\} \mathchoice{\frac{\partial u_1}{\partial t}}{\partial u_1/\partial t}{\partial u_1/\partial t}{\partial u_1/\partial t} &= d_1(u_0,v_0)\nabla^2 u_1 + d_2(u_0,v_0) \nabla^2 v_1 + F_u(u_0,v_0)u_1+F_v(u_0,v_0)v_1, \\ \mathchoice{\frac{\partial v_1}{\partial t}}{\partial v_1/\partial t}{\partial v_1/\partial t}{\partial v_1/\partial t} &= d_3(u_0,v_0)\nabla^2 u_1 + d_4(u_0,v_0) \nabla^2 v_1 + G_u(u_0,v_0)u_1+G_v(u_0,v_0)v_1. \end{align} As before, we can write this as a system of equations by putting $\mathbfit{u}_1 = (u_1,v_1)$. We then have, \begin{equation} \mathchoice{\frac{\partial\mathbfit{u}_1}{\partial t}}{\partial\mathbfit{u}_1/\partial t}{\partial\mathbfit{u}_1/\partial t}{\partial\mathbfit{u}_1/\partial t} = \mathsfbfit{D}\mathbfit{u}_1 + \mathsfbfit{J}\mathbfit{u}_1, \quad \mathsfbfit{J} = \begin{pmatrix} F_u & F_v \\ G_u & G_v \end{pmatrix},\] where the functions in the matrices are evaluated at the homogeneous steady state, \((u_0,v_0)\). For simplicity we will omit this dependence in all of the partial derivatives of \(F\) and \(G\), and in the \(d_i\).

We can expand these perturbations as \(\mathbfit{u}_1 \propto \mathrm{e}^{\lambda_k t}w_k(\mathbfit{x})\), with the \(w_k\) again satisfying \(\nabla^2 w_k + \rho_k w_k = 0\). Substituting this in, we find that the growth rates \(\lambda_k\) are eigenvalues of the matrix, \[\mathsfbfit{M} = -\rho_k \mathsfbfit{D} + \mathsfbfit{J} = \begin{pmatrix} F_u-\rho_k d_1 & F_v-\rho_k d_2 \\ G_u-\rho_k d_3 & G_v-\rho_k d_4 \end{pmatrix}.\] We then have that the growth rates satisfy, \[\lambda_k^2 - \operatorname{tr}(\mathsfbfit{M})\lambda_k + \det(\mathsfbfit{M}) = 0, \quad \operatorname{tr}(\mathsfbfit{M}) = F_u+G_v-\rho_k(d_1+d_4),\] and we have, \[h(\rho_k) = \det(\mathsfbfit{M}) = \rho_k^2(d_1d_4-d_2d_3) + \rho_k(d_2G_u+d_3F_v - d_1G_v - d_4F_u) + F_uG_v - F_vG_u.\]

We will now derive conditions for pattern formation. We first assume that \(\mathsfbfit{J}\) has eigenvalues with negative real part, and hence we have, \[\operatorname{tr}(\mathsfbfit{J}) = F_u + G_v < 0,\quad \det(\mathsfbfit{J}) = F_uG_v - F_vG_u > 0.\] By assumption \(\mathsfbfit{D}\) is positive-definite, and hence it must have a positive trace and determinant. So we have, \[d_1+d_4 > 0, \quad d_1d_4-d_2d_3 > 0.\] From this last condition we have that the coefficient of \(\rho_k^2\) in \(h(\rho_k)\) is positive, and the constant term in this function is also equivalent to \(\det(\mathsfbfit{J})\), which is positive. We also have that \(\operatorname{tr}(\mathsfbfit{M}) < \operatorname{tr}(\mathsfbfit{J}) < 0\), and hence must again focus our attention on values of \(\rho_k\) where \(h(\rho_k)<0\).

To find the minimum value of \(h\) We compute, \[\mathchoice{\frac{{\mathrm d}h}{{\mathrm d}\rho_k}}{{\mathrm d}h/{\mathrm d}\rho_k}{{\mathrm d}h/{\mathrm d}\rho_k}{{\mathrm d}h/{\mathrm d}\rho_k} = 2(d_1d_4-d_2d_3)\rho_k^* + d_2G_u+d_3F_v - d_1G_v - d_4F_u = 0,\] and hence, \[\rho_k^* = \frac{d_1G_v + d_4F_u-d_2G_u-d_3F_v}{2(d_1d_4-d_2d_3)}.\] For this to be positive we must have, \[d_1G_v + d_4F_u > d_2G_u+d_3F_v,\] which is a generalisation of the third Turing condition given in [turingconditions] (in this notation, set \(d_1=1\), \(d_4=D\), and \(d_2=d_3=0\) to exactly recover the result in the reaction–diffusion case).

Substituting this critical value in and simplifying, we find \[h(\rho_k^*) = -\frac{(d_1G_v + d_4F_u-d_2G_u-d_3F_v)^2}{4(d_1d_4-d_2d_3)} + F_uG_v - F_vG_u.\] Hence for \(h(\rho_k^*)<0\) we must have, \[(d_1G_v + d_4F_u-d_2G_u-d_3F_v)^2-4(d_1d_4-d_2d_3)(F_uG_v - F_vG_u) > 0.\] This is precisely a generalisation of the fourth Turing condition given in [turingconditions]; to see this, again set \(d_1=1\), \(d_4=D\), and \(d_2=d_3=0\) to exactly recover the result in the reaction–diffusion case.

On Problem Sheet 8 there is an example where you will use these conditions to find conditions for pattern formation in a predator–prey model with prey moving towards the predator. Examples of these kinds of ambush predators are given in 6.3.

13.2.1 Relation to PDEs

We can relate these nonlocal models to PDEs, like the reaction–diffusion equations we have seen, by choosing a particular form of the kernel \(K\) and interaction term \(F\) which models a kind of ‘local’ interaction. We set \(K(x,y) = K(x-y)=K(y-x)=K(|x-y|)\), so that the kernel is only a function of the distance between two points. We will also assume that this influence decays rapidly as the distance increases; that is \(K(x-y) \to 0\) for \(|x-y| \gg 1\). An example of such a kernel is, \[\begin{equation} \label{example-US-kernel} K(|x-y|) = \mathrm{e}^{-s(x-y)^2}. \end{equation}\] In this case, especially for \(s \gg 1\), this function approaches a very sharp and narrow limit around the point \(x=y\) (which will increasingly approximate a \(\delta\)-function for large \(s\)) so that only nearby regions interact. For the interaction, we simply take the linear function \(F(u) = u\).

We will now use our assumptions on this kernel to perform a Taylor series of the integral term in [nonlocal-US-PDE]. Before we do this, we note some properties of this kind of kernel. By the assumed symmetry of the argument, we see that \(K(y)\) is an even function of \(y\). Hence we have, \[\begin{equation} \label{K-US-symmetry} \int_{-\infty}^{\infty}y^{2m+1}K(y)\,{\mathrm d}y = 0, \quad \textrm{for } m=0,1,2,3,\dots \end{equation}\] We also define the moments of the kernel (which are just numbers) by, \[K_{2m} = \frac{1}{(2m)!}\int_{-\infty}^{\infty}y^{2m}K(y)\, {\mathrm d}y, \quad \textrm{for } m=0,1,2,3,\dots\]

Now we are ready to expand the integral term in [nonlocal-US-PDE], under the assumption that the kernel decays sufficiently quickly (for example, \(s \gg 1\) in [example-US-kernel]). In this case, we can think of the quantity \(z = x-y\) as being ‘small’ and use a Taylor expansion. We have that the integral term looks like, \[\begin{align} \label{horrific-US-expansion} &\int_{-\infty}^{\infty} K(x-y) u(y,t) \, {\mathrm d}y = \int_{-\infty}^{\infty} K(z) u(x-z,t) \, {\mathrm d}z \\ &=\int_{-\infty}^{\infty} K(z) \left [u(x,t)-z\mathchoice{\frac{\partial u}{\partial x}}{\partial u/\partial x}{\partial u/\partial x}{\partial u/\partial x}(x,t)+\frac{z^2}{2}\mathchoice{\frac{\partial^2 u}{\partial x^2}}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}(x,t)-\frac{z^3}{6}\mathchoice{\frac{\partial^{3} u}{\partial x^{3}}}{\partial^{3} u/\partial x^{3}}{\partial^{3} u/\partial x^{3}}{\partial^{3} u/\partial x^{3}}(x,t) + \frac{z^4}{24}\mathchoice{\frac{\partial^{4} u}{\partial x^{4}}}{\partial^{4} u/\partial x^{4}}{\partial^{4} u/\partial x^{4}}{\partial^{4} u/\partial x^{4}}(x,t) +\cdots \right] {\mathrm d}z\\ &=u(x,t)\int_{-\infty}^{\infty} K(z) \, {\mathrm d}z-\mathchoice{\frac{\partial u}{\partial x}}{\partial u/\partial x}{\partial u/\partial x}{\partial u/\partial x}(x,t)\int_{-\infty}^{\infty} zK(z) \, {\mathrm d}z +\mathchoice{\frac{\partial^2 u}{\partial x^2}}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}(x,t)\int_{-\infty}^{\infty} \frac{z^2}{2}K(z) \, {\mathrm d}z+\cdots\\ &=u(x,t)K_0 +\mathchoice{\frac{\partial^2 u}{\partial x^2}}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}(x,t)K_2+\mathchoice{\frac{\partial^{4} u}{\partial x^{4}}}{\partial^{4} u/\partial x^{4}}{\partial^{4} u/\partial x^{4}}{\partial^{4} u/\partial x^{4}}(x,t)K_4+\mathchoice{\frac{\partial^{6} u}{\partial x^{6}}}{\partial^{6} u/\partial x^{6}}{\partial^{6} u/\partial x^{6}}{\partial^{6} u/\partial x^{6}}(x,t)K_6+\cdots, \end{align}\] where all of the odd terms vanish due to the symmetry of \(K\) (see [K-US-symmetry]). This expansion looks scary, but the point is just that if the integral kernel \(K\) is sufficiently smooth and decays quickly enough, then you can approximate the integral term by spatial derivatives.

We note this expansion only makes sense if the \(K_{2m}\) get smaller sufficiently quickly. For [example-US-kernel], we can compute (not trivial but not immensely hard): \[K_{2m} = \frac{\sqrt{\pi}}{\sqrt{s}(4s)^m m!}.\] These terms tend towards \(0\), and generally quite rapidly, though there are some technicalities in ensuring a convergent expansion for [horrific-US-expansion].²⁸ We will not dwell on these, and instead study what happens when we truncate the series. To see how quickly these terms converge we consider \(s=1\), and compute that \(K_2 \approx 0.443\), \(K_4 \approx 0.055\), \(K_6 \approx 0.005\), \(K_8 \approx 0.0003\), etc. We note that the signs of these moments can be different, especially when negative interactions are modelled.

Given the numerical evidence above, and for simplicity, we will truncate this expansion at the fourth order to demonstrate the impact of this kind of nonlocality on instabilities. Our equation then looks like, \[\begin{equation} \label{biharmonic-US-PDE-US-pre} \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} = uK_0 +K_2\mathchoice{\frac{\partial^2 u}{\partial x^2}}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}+K_4\mathchoice{\frac{\partial^{4} u}{\partial x^{4}}}{\partial^{4} u/\partial x^{4}}{\partial^{4} u/\partial x^{4}}{\partial^{4} u/\partial x^{4}} + G(u). \end{equation}\] We can simplify this by writing \(H(u) = uK_0+G(u)\) to find, \[\begin{equation} \label{biharmonic-US-PDE} \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} = K_2\mathchoice{\frac{\partial^2 u}{\partial x^2}}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}+K_4\mathchoice{\frac{\partial^{4} u}{\partial x^{4}}}{\partial^{4} u/\partial x^{4}}{\partial^{4} u/\partial x^{4}}{\partial^{4} u/\partial x^{4}} + H(u). \end{equation}\] We will study this equation assuming that it is posed on a finite spatial domain, \(x \in [0,L]\), in order to avoid some technicalities about eigenfunctions of the Laplacian defined on \(\mathbb{R}\). This can be made sense of if the above derivation was done on a bounded domain, and the limit taken more carefully near the boundaries. It also helps explain why we think of higher-order derivatives as being less ‘local,’ than lower-order ones.

We can impose any boundary conditions that make biological sense, but as the equation is now higher-order in space, we must have two boundary conditions at each endpoint. We will discuss different choices of these boundary conditions in the next subsection.

You have in fact seen an equation with higher-order spatial derivatives before, namely, the beam equation from Michaelmas. This equation had the form, \[\mathchoice{\frac{\partial^{4} w}{\partial s^{4}}}{\partial^{4} w/\partial s^{4}}{\partial^{4} w/\partial s^{4}}{\partial^{4} w/\partial s^{4}}+\frac{N}{EI}\mathchoice{\frac{\partial^2 w}{\partial s^2}}{\partial^2 w/\partial s^2}{\partial^2 w/\partial s^2}{\partial^2 w/\partial s^2}+\rho A\mathchoice{\frac{\partial^2 w}{\partial t^2}}{\partial^2 w/\partial t^2}{\partial^2 w/\partial t^2}{\partial^2 w/\partial t^2} = 0,\] (where \(w\) was the unknown displacement and \(s\) was the spatial variable). This is not quite the same form as [biharmonic-US-PDE], as the beam equation has a second order derivative in time, but you can in fact study it using the methods of this section by reducing it to a first-order system.

13.2.2 Linearisation and spectral analysis

We can again study homogeneous equilibria of [biharmonic-US-PDE] by finding numbers \(u_0\) such that \(H(u_0)=0\), as these will identically satisfy the right-hand side of the equation being zero. Following the usual procedure of writing \(u = u_0 + \varepsilon u_1\) and substituting this into [biharmonic-US-PDE], we find that perturbations satisfy, \[\mathchoice{\frac{\partial u_1}{\partial t}}{\partial u_1/\partial t}{\partial u_1/\partial t}{\partial u_1/\partial t} = K_2\mathchoice{\frac{\partial^2 u_1}{\partial x^2}}{\partial^2 u_1/\partial x^2}{\partial^2 u_1/\partial x^2}{\partial^2 u_1/\partial x^2}+K_4\mathchoice{\frac{\partial^{4} u_1}{\partial x^{4}}}{\partial^{4} u_1/\partial x^{4}}{\partial^{4} u_1/\partial x^{4}}{\partial^{4} u_1/\partial x^{4}} + H'(u_0)u_1.\] If we perform a separation-of-variables solution, we can use an ansatz of the form \(u_1 = e^{\lambda_k t}f(x)\) to find an eigenvalue equation for the function \(f\) of the form, \[\begin{equation} \label{eigen-US-biharmonic} \lambda_k f = K_2\mathchoice{\frac{{\mathrm d}^2 f}{{\mathrm d}x^2}}{{\mathrm d}^2 f/{\mathrm d}x^2}{{\mathrm d}^2 f/{\mathrm d}x^2}{{\mathrm d}^2 f/{\mathrm d}x^2}+K_4\mathchoice{\frac{{\mathrm d}^{4} f}{{\mathrm d}x^{4}}}{{\mathrm d}^{4} f/{\mathrm d}x^{4}}{{\mathrm d}^{4} f/{\mathrm d}x^{4}}{{\mathrm d}^{4} f/{\mathrm d}x^{4}} + H'(u_0)f. \end{equation}\] This is a fourth-order ODE which can be solved to find four linearly independent solutions. Using the boundary conditions, as we did in the reaction–diffusion case, we can compute the value of the constant \(\lambda_k\). But this quickly becomes very tedious!

Alternatively, the following boundary conditions allow the use of a clever trick: \[\begin{equation} \label{biharmonic-US-BCs} \mathchoice{\frac{\partial u}{\partial x}}{\partial u/\partial x}{\partial u/\partial x}{\partial u/\partial x}(0,t) = 0, \quad \mathchoice{\frac{\partial u}{\partial x}}{\partial u/\partial x}{\partial u/\partial x}{\partial u/\partial x}(L,t) = 0, \quad \mathchoice{\frac{\partial^{3} u}{\partial x^{3}}}{\partial^{3} u/\partial x^{3}}{\partial^{3} u/\partial x^{3}}{\partial^{3} u/\partial x^{3}}(0,t) = 0, \quad \mathchoice{\frac{\partial^{3} u}{\partial x^{3}}}{\partial^{3} u/\partial x^{3}}{\partial^{3} u/\partial x^{3}}{\partial^{3} u/\partial x^{3}}(L,t) = 0. \end{equation}\] We let \(w_k(x) = \cos(k \pi x/L)\) be the usual Laplacian eigenfunctions satisfying Neumann boundary conditions, and recall \[\mathchoice{\frac{\partial^2 w_k}{\partial x^2}}{\partial^2 w_k/\partial x^2}{\partial^2 w_k/\partial x^2}{\partial^2 w_k/\partial x^2} + \rho_k w_k = 0,\] with \(\rho_k = (k\pi/L)^2\). We then note that these functions \(w_k(x)\) satisfy all of the boundary conditions given in [biharmonic-US-BCs], and hence we can use the ansatz²⁹ \(u_1 = e^{\lambda_k t}w_k(x)\) to find, \[\begin{equation} \label{biharmonic-US-growth-US-rates} \lambda_k = -\rho_kK_2+\rho_k^2K_4 + H'(u_0). \end{equation}\]

If we assume, as in the Turing analysis from preceding chapters, that the system is stable in the absence of transport, then we must have \(\lambda_0 = H'(u_0)<0\) at the steady state. We also require that \(\operatorname{Re}(\lambda_k)<0\) as \(k \to \infty\), as otherwise we obtain instability of arbitrarily fine-scale patterns, which effectively break our assumptions of a spatial continuum. This implies that \(K_4<0\). Hence the only possibility for some finite range of \(k\) to lead to \(\operatorname{Re}(\lambda_k)>0\) is if \(K_2<0\).

13.2.3 Example: Non-local spatial patterning in a mechanical strain model

Let’s consider an example interaction scheme where we observe spatial pattern formation in a relatively simple system. From the above analysis, we need \(K_2<0\) and \(K_4<0\). One kind of model which can, for the right parameters, admit these signs is a kernel with oscillations between activation and inhibition in space. That is, we choose a \(K\) such that nearby populations benefit one another, but ones farther away inhibit one another. See 6.4 for a sketch of such an interaction kernel. An example is given by the function, \[\begin{equation} %\label{example-US-kernel} K(x-y) = \mathrm{e}^{-(x-y)^2}(\cos(x-y)-\alpha). \end{equation}\]

This kernel has the moments, \[K_2 = \frac{\sqrt{\pi}}{8}\left(\mathrm{e}^{-\frac{1}{4}}-2\alpha\right), \quad K_4 = \frac{\sqrt{\pi}}{384}\left(\mathrm{e}^{-\frac{1}{4}}-12\alpha\right),\] which are both negative for \(2\alpha>\mathrm{e}^{-\frac{1}{4}}\).

An important interpretation of such a structure is that the population \(u\) acts as its own activator nearby, but its own inhibitor far away. Hence it can exhibit the same kind of structures we observed in chapter 3 for reaction–diffusion systems, but in a single species. Such interactions exist in mechanical and neural interactions, where cells interact over long distances via filopodia or mechanical signals, as well as in flocking and aggregation phenomena.

We consider a simple form of a mechanical interaction as an example of this kind of instability analysis. If we let \(u\) denote the deformation of a sheet of cells, such as the epithelium or outermost layer of your skin, we can consider such an interaction kernel with \(K_2 = -a <0\) and \(K_4 = -b<0\) (with \(a,b\) positive) to represent the influence of curvature on the cell layer itself. In the absence of curvature, we assume the cells want to relax to an undeformed state (that is, \(u_0=0\) is the only homogeneous equilibrium). A simple model of this would take \(H(u) = -u-u^3\). Our equation would then be, \[\begin{equation} \label{mechanical} \mathchoice{\frac{\partial u}{\partial t}}{\partial u/\partial t}{\partial u/\partial t}{\partial u/\partial t} = -a\mathchoice{\frac{\partial^2 u}{\partial x^2}}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}{\partial^2 u/\partial x^2}-b\mathchoice{\frac{\partial^{4} u}{\partial x^{4}}}{\partial^{4} u/\partial x^{4}}{\partial^{4} u/\partial x^{4}}{\partial^{4} u/\partial x^{4}} -u-u^3, \end{equation}\] and we see that linearising it is equivalent to just dropping the \(-u^3\) term. Hence by [biharmonic-US-growth-US-rates], we have \[\lambda_k = \rho_ka-\rho_k^2b -1.\]

This dispersion relation is quite simple, and it is easy to calculate, e.g., the range of unstable wavenumbers, and the fastest growing mode. We compute the range of unstable wavenumbers by solving \(\lambda_k=0\) as a function of \(\rho_k\) to find, \[\rho_k^\pm = \frac{a\pm \sqrt{a^2-4b}}{2b},\] from which we see that \(a^2 > 4b\) is necessary for instability (as otherwise the quadratic given by \(\lambda_k\) would never become positive). If this condition is satisfied, then the values of \(\rho_k\) for which \(\lambda_k=0\) are real and positive.

If we maximise \(\lambda_k\) as a function of \(\rho_k\) we find the fastest growing mode as, \[\lambda_k(\rho_k^*) =\lambda_k\left(\frac{a}{2b}\right) = \left(\frac{a^2}{2b}\right)-\left(\frac{a^2}{4b}\right) -1 = \frac{a^2}{4b}-1,\] which is positive as long as \(a^2 > 4b\). You can see simulations of this model in 6.5 for two different domain lengths, in order to see how the instability regimes scale with size. Importantly, as with reaction–diffusion systems, the linear analysis can approximately tell you what wavelength the pattern will be, but it cannot say whether spots or stripes form.