Reaction - diffusion equations

This last section deals with the existence of equilibria in models that govern reaction-diffusion systems. These equilibria are critical points at which the system is in equilibrium state. A question that is of interest is the following: if the equilibrium is perturbed, does the system act such that it returns to the equilibrium state or does it move on to a new equilibrium? If the system returns to the equilibrium state after being perturbed, we refer to the equilibrium as a stable point while if the system evolves into a different state, we refer to the equilibrium as an unstable point. The notion of equilibria and their stability should be familiar to you from Year 1 mathematics in the context of systems of ordinary differential equations (e.g. predator-prey models). A similar theory of the stability of equilibria can be developed for partial differential equations which will be our focus here.

The typical form of a reaction-diffusion equation (1D) is given as:

where $u=u(x, t)$. The $D u_{x x}$ term describes diffusion and $D$ is a diffusion coefficient. The function $f(u)$ describes a process in which $u$ changes rather than just diffusing in space. So, $f(u)$ represents the reaction term. Note that $f$ can also depend on $u_{x}$ or $x$ as well in a more general form of Eq. (2.211).

Typical reactions

This section gives a quick review of equilibria and stability in equations which model typical reactions present in reaction-diffusion systems. Note that in the absence of 'diffusion' terms which are given by spatial changes, the differential equation reduces to an ordinary differential equation describing 'reaction' phenomena, such as:

Eq. (2.212) has equilibria at $u^{\prime}=0$ (where the prime denotes differentiation with respect to time) and their stability depends on the function $f(u)$. For instance, let us look at a simple ODE modelling logistic growth population:

where $r$ and $k$ represent the intrinsic growth rate and the carrying capacity, respectively. We have 2 equilibria occuring at $u^{\prime}=0: u^{*}=0$ and $u^{*}=k$. These are the points in $u$ for which there is no change in $u$ with time. Mathematically, if we have initial conditions that start at, say, $u(0)=k$, then we stay on $u=k$ for all time $t$ since at $u=k, u^{\prime}=0$. Our next question is whether the states 0 and $k$ are stable or unstable. This is determined easily in single equations like Eq. (2.213). We already know where $u^{\prime}=0$; we now need to find where $u^{\prime}>0$ and $u^{\prime}<0$. The significance of this is that if $u^{\prime}>0$, this means that the rate of change of $u$ is positive and hence $u$ is increasing. If $u^{\prime}<0$, then the solutions $u$ are decreasing. This is shown in Fig. 2.9. To the left of $u^{*}=0, f(u)<0$ hence solutions decrease (this is mathematically sound but physically, solutions exist for $u \geq 0$ since $u$ represents population density) while to the right of $u^{*}=0, f(u)>0$ hence solutions increase. Since near $u^{*}=0$, solutions are moving away from the equilibrium, the point is unstable; the motion away from equilibria is indicated by green arrows in Fig. 2.9. With similar arguments, we observe that the point $u^{*}=k$ is stable.

Figure 2.9: A plot of $f(u)=r u\left(1-\frac{u}{k}\right)$ against $u$ [Eq. (2.213)]. The circle markers indicate equilibria at $u^{*}=0, k$ corresponding to $u^{\prime}=f(u)=0$. The 0 point is unstable while the $k$ point is stable.

This approach is qualitative; we focus on the behaviour of solutions near equilibria but we are not concerned with exact solutions in closed form or numerical approximations. The method finds use in various applications in the life and physical sciences as well as engineering settings.

Diffusion

To model diffusion, a classical approach is through application of mass conservation and Fick's law. The linear, one-dimensional diffusion equation is given by:

We solve the diffusion equation (2.214) on a bounded interval; for instance, we can have no-flux boundary conditions

where $[0, L]$ represents the interval in which we wish to find solutions to the problem. Following Section 2.4 this problem can be solved by the separation of variables method which assumes that $u(x, t)$ may be expressed as a product of a function of space only, $X(x)$ and time, $T(t)$.

Single reaction-diffusion equation

Let us go back to Eq. (2.211) and analyse the partial differential equation using the notion of stability of solutions and the method of separation of variables as applied previously in these notes.

Consider the reaction-diffusion equation given by Eq. (2.211) with no-flux boundary conditions,

Let $u(x, t)=u^{*}$ represent a constant equilibrium solution such that $f\left(u^{*}\right)=0$ and assume that $f^{\prime}\left(u^{*}\right)>0$.

In order to assess the stability of the equilibrium, we introduce a small perturbation, $\eta(x, t)$ which is both spatially and temporally dependent, i.e.:

We then substitute Eq. (2.217) into the PDE (2.211) and the boundary conditions (2.216), which gives:

Perturbations are typically very small i.e. $|\eta|<<1$, but if the equilibria are unstable, even tiny perturbations can result in major changes (often undesirable) in the system behaviour. Since $|\eta|<<1$, we proceed by expanding $f\left(u^{*}+\eta\right)$ using a Taylor series and discarding the nonlinear terms. This gives the linearised equation:

which is subject to the boundary conditions included in Eq. (2.218). Note that $f^{\prime}$ is the derivative of $f$ with respect to $u$ and it is evaluated at $u^{*}$; so $f^{\prime}\left(u^{*}\right)$ is simply a constant whose sign dictates the stability of the solution as we will see later.

Observing that we now have a linear PDE given by (2.219) on a bounded interval, we start with the separable solution,

Substituting in Eq. (2.219) gives:

where $a=f^{\prime}\left(u^{*}\right)$. The first equation for $T(t)$ yields $T(t)=c_{1} e^{\lambda t}$. The boundary-value problem for $X(x)$ is given as follows:

Note that for the roots of the associated characteristic equation to be real and distinct, $a-\lambda<0$ but such exponential solutions cannot satisfy the boundary conditions unless they are trivial (i.e. zero exponential solutions). Since we are not interested in trivial solutions we let $a-\lambda \geq 0$ and allow an oscillatory solution (made up sines and cosines). For $\lambda=a$, we obtain a constant solution, $X(x)=c_{2}$ for any nonzero $c_{2}$. For $a-\lambda>0$, we have the general solution,

Applying $X^{\prime}(0)=0$ gives $c_{4}=0$ and $X^{\prime}(\pi)=0$ gives:

The solutions for $\eta(x, t)$ are therefore given by:

The constants $c_{n}$ are determined by the initial perturbation given at $t=0$ by $\eta(x, 0)=f(x)$ where $f(x)$ is arbitrary at this point. Equation (2.225) tells us that the perturbation $\eta$ dies down if all the modes (i.e. for all $n$ ) decay and it will grow if one of the modes grows. The term given by $e^{\lambda_{n} t}$ therefore determines growth or decay; this is decided by the eigenvalues $\lambda_{n}$. For instance, at $n=0$, Eq. (2.225) gives $\eta=c_{0} e^{a t}$ and since $a>0$ (this was assumed in the problem statement), then the perturbation grows exponentially. For $n \neq 0$, perturbation grow if $a>D n^{2}$; this is most likely going to happen when $n$ is low (i.e. for low frequency modes). We can conclude that low frequency modes are destabilizing (growing perturbations render the system unstable) while high frequency modes are stabilizing. Note that while our analysis was carried out for a constant equilibrium solution, it is possible to follow the same method for nonconstant equilibrium solutions as well.

System of two linear reaction-diffusion equations

Generally, much more complex behaviour is observed in systems of equations rather than a single equation. Based on the seminal work of Alan Turing (a British mathematician who contributed to the development of computer science and was able to crack German codes during WWII), reaction-diffusion models can form the basis of pattern formation such as stripes ans spots observed in animal skin.

We now consider two interacting chemicals, $u_{1}$ and $u_{2}$ in a reaction-diffusion system of two PDEs:

The functions $f, g$ are typically nonlinear functions governing the reaction kinetics while $D_{1}$ and $D_{2}$ are diffusion coefficients for each chemical. For each species, we use the same no-flux boundary conditions as in the previous section:

where $i=1,2$.

Now, we assume that when the mixture is well-mixed (i.e. when $D_{1}=D_{2}=0$ ), there exists a state given by $\left(u_{1}^{*}, u_{2}^{*}\right)$ which makes $u_{1 t}=u_{2 t}=0$. It follows that at $\left(u_{1}^{*}, u_{2}^{*}\right)$, we have $f\left(u_{1}^{*}, u_{2}^{*}\right)=g\left(u_{1}^{*}, u_{2}^{*}\right)=0$. This represents the equilibrium state of the system; note that there may exist more than one set of equilibria for a given system. Here, we assume that this equilibrium is stable for the wellmixed case. Again, what we are interested in is, upon perturbing the system, does the system return to equilibrium or not? The approach discussed for the scalar equation (single equation) above applies here as well: we introduce small space- and time-dependent perturbations to identify whether the equilibrium may be destabilized in the presence of diffusion. We introduce the following perturbations in Eqs. (2.226) and (2.227) and boundary conditions (2.228)

where $\eta_{1}$ and $\eta_{2}$ are small disturbances. Substitution of Eqs. (2.229) and linearising using a Taylor series yields the following equations for the spatiotemporal evolution of the perturbations, after discarding any nonlinear terms:

where the coefficients $a_{i j}$ are the following partial derivatives evaluated at the equilibrium point, $\left(u_{1}^{*}, u_{2}^{*}\right)$ :

${ }^{4}$ Note that these are nonlinear terms in $\eta_{1}$ and $\eta_{2}$. The idea is if the perturbations are small the nonlinear terms such as $\eta_{1}^{2}$ for example are super small and therefore negligible. Note that the set of equations (2.230) and (2.231) are the analogue of Eq. (2.219) discussed in Subsec. 2.9.3 Again, the method of separation of variables is applicable here and, based on the analysis of the scalar equation in Subsec. 2.9.3 as well as the material discussed in Subsec. 2.3.2, we can use sines and cosines for the spatial frequencies. Since we have no-flux conditions, sines cannot satisfy the boundary conditions so we are left with a cosine series expansion of various frequencies.

To make progress therefore we introduce the following ansatz for the modal solutions of the perturbations:

Note that the argument of the cosine, $n x$ is a result of no-flux conditions at $x=0$ and $x=\pi$; in the more general case where the conditions are given for the finite interval $[0, L]$, the argument is $n \pi x / L$. Further, $\lambda$ here is an unknown that we need to determine to assess whether perturbations grow or decay. Upon substitution of Eqs. (2.233) in Eqs. (2.230) and (2.231) and after cancelling the common factor $e^{\lambda t} \cos (n x)$, we have the following system of algebraic equations for the coefficients $c_{1}$ and $c_{2}$

The above system is linear and homogeneous which implies that the trivial solution $c_{1}=c_{2}=0$ satisfies the two equations. We look for nontrivial coefficients however since trivial coefficients imply no perturbation [this is obvious from Eqs. (2.233)]. What we have is two equations with three unknowns: $c_{1}$, $c_{2}$ and $\lambda$. We rewrite the system of equations in matrix form as follows,

and from linear algebra we know that nontrivial solutions exist for $c_{1}$ and $c_{2}$ if the determinant of the matrix is zero:

This reduces to a quadratic equation for $\lambda$ :

where the constant $p$ is given by $p=n^{2}\left(D_{1}+D_{2}\right)-a_{11}-a_{22}$ and $q$ is given by $q=D_{1} D_{2} n^{4}-\left(a_{11} D_{2}+a_{22} D_{1}\right) n^{2}+a_{11} a_{22}-a_{12} a_{21}$. From Eqs. (2.233), the sign of $\lambda$ dictates whether perturbations grow or decay. From the definitions of $p$ and $q$, the growth or decay of perturbations therefore depends on the equilibrium state since the $a_{i j}$ coefficients are evaluated at the equilibrium point, the diffusion coefficients $D_{1}$ and $D_{2}$ and the frequency $n$. The roots for $\lambda$ are obtained using:

We now draw the following conclusions:

• If $p>0, q>0$ then $\lambda<0$ and therefore perturbations decay and the system returns to equilibrium.
• Let us return to the well-mixed assumption when $D_{1}=D_{2}=0$. The original nonlinear system of PDEs reduces to a nonlinear system of ODEs which is given below:

where primes denote differentiation with respect to time since spatial effects vanish. Now, we assumed that there exists a stable equilibrium state $\left(u_{1}^{*}, u_{2}^{*}\right)$ when the system is well-mixed. The criteria for stability/instability are given by the eigenvalues of the Jacobian 5 of the system when evaluated at the equilibrium point, i.e.

The system is stable if $\operatorname{trace}\left(J\left(u_{1}^{*}, u_{2}^{*}\right)\right)<0$ and $\operatorname{det}\left(J\left(u_{1}^{*}, u_{2}^{*}\right)\right)>0$; the negative trace of the Jacobian (2.240) implies that $a_{11}+a_{22}<0$. Back to the general case now where $D_{1}$ and $D_{2}$ are not zero. Assuming that the system starts off as stable when well-mixed, then $p>0$ (i.e. $p$ is always positive) since

The importance of this is that if the system is stable when well-mixed, then adding diffusion does not affect the stability because of changes in the sign of $p$. What may affect the stability therefore is the sign of $q$. In particular, if $\left(a_{11} D_{2}+a_{22} D_{1}\right) n^{2}>D_{1} D_{2} n^{4}+a_{11} a_{22}-a_{12} a_{21}$, then one of the eigenvalues has a positive real part which leads to the growth of the perturbations thus rendering the equilibrium unstable.

${ }^{5}$ The presence of the Jacobian of the system is a result of linearising the nonlinear system and the criteria for stability arise from eigenvalues with negative real part such that perturbations of the form $e^{\lambda t} h(x)$ decay if $\operatorname{Re}(\lambda)<0$.