Instabilities in Reaction-Diffusion Systems

Posted on May 17, 2020 by Andrew Krause

Here I want to provide a very brief introduction to instabilities in some generic spatial equations, and in particular reaction-diffusion systems on continuous and discrete (network) domains. Turing instabilities in such systems provide a way to generate spatially heterogeneous states of such systems which exhibit a variety of interesting and complex structures. Such patterns are often used to try and explain emergent structures in developmental biology, such as the stripes and spots on animal coat patterns.

I will begin with a theoretical approach to the simplest method for the most general problem, and then show how this applies to particular relevant examples. Throughout the following discussion there will be subtle technical points (particularly regarding analytic details such as regularity of functions, topology of sets, etc) but I will omit almost all nonessential technical points to try and just illustrate the ideas in a useful way.

1. Linear Instability Theory in {\mathbb{R}^m}

Consider an autonomous system of ordinary differential equations given by,

\displaystyle \frac{\textrm{d} \mathbf{u}}{\textrm{d} t} = \mathbf{F}(\mathbf{u}), \quad \mathbf{u}(0) = \mathbf{u_0}, \ \ \ \ \ (1)

where {\mathbf{u}(t) \in \mathbb{R}^m} for each {t\geq 0} and {\mathbf{F}} is a sufficiently `nice’ map to guarantee, e.g., global in time existence of solutions. Let {\mathbf{u^*}\in \mathbb{R}^m} be an equilibrium, so that {\mathbf{F}(\mathbf{u^*})=\mathbf{0}}. If we perturb this equilibrium by a small amount, we can track how perturbations grow (or decay) by setting {\mathbf{u} = \mathbf{u^*}+\epsilon \mathbf{v}} where {\epsilon \ll 1}, and substituting this into the ODE above. We get an equation for {\mathbf{v}} which we can Taylor expand as,

\displaystyle \epsilon \frac{\textrm{d} \mathbf{v}}{\textrm{d} t} = \mathbf{F}(\mathbf{u^*}+\epsilon \mathbf{v}) = \mathbf{F}(\mathbf{u^*})+\epsilon\underbrace{\frac{\partial \mathbf{F}}{\partial \mathbf{u}}\left(\mathbf{u^*}\right)}_{\mathbf{J}}\mathbf{v}+O(\epsilon^2)= \epsilon\mathbf{J}\mathbf{v} +O(\epsilon^2), \ \ \ \ \ (2)

where {\mathbf{J}} is the Jacobian of the function {\mathbf{F}} evaluated at the equilibrium {\mathbf{u^*}}.

After canceling an {\epsilon} above, and neglecting the higher-order terms, we have that {\mathbf{v}} satisfies a simple linear system of equations involving the Jacobian of the original nonlinear function at the equilibrium state. From the theory of such linear systems, we know that its long-time behaviour is determined by the eigenvalues of {\mathbf{J}}, say {\lambda_1,\lambda_2,\dots,\lambda_m}. In particular, if we have that all of their real parts are negative, i.e. {\Re(\lambda_i) < 0} for all {i = 1,\dots,m}, then the equilibrium is stable to small perturbations. If instead {\Re(\lambda_i)>0} for at least one eigenvalue, then perturbations will grow, and the original equilibrium is unstable. At the boundary states of marginal stability ({\Re(\lambda_i)=0} for one or more eigenvalues, with {\Re(\lambda_j)<0} for all others), the behaviour of perturbations is harder to classify. Away from these states, it turns out that the linear system above directly tells how trajectories of the original system evolve near the equilibrium, via the Hartman-Grobman Theorem.

So to find if an equilibrium point is linearly unstable, we can simply look for parameters or equilibria that admit eigenvalues with positive real part. Many aspects of different kinds of instabilities are studied in bifurcation theory (particularly local bifurcation theory, which studies instabilities of equilibrium points), but for now we will move on to instabilities in PDEs.

2. General Stability Theory

Here we will consider something analogous to the above linear instability theory for somewhat generic transport-type equations. Instability theory for PDEs can be much more involved, and involve complicated spectral questions which I will not address here, but instead refer you to these introductory notes. In what follows, I will essentially assume that all operators have only point spectra, or that we can otherwise ignore these technicalities.

Consider a vector-valued function {\mathbf{u}(\mathbf{x},t)\in \mathbb{R}^n}, with the spatial domain {\Omega \subset \mathbb{R^m}} which will be taken as either a set of points in the case of discrete operators, or as a manifold sufficiently smooth to do calculus on (you can think of the interval, {\Omega=[0,L]}, for some length {L}, as the simplest example). Here we will also assume this manifold is compact (i.e. bounded and containing its boundary), though other possibilities exist. We then assume that these functions evolve according to the following reaction-transport system,

\displaystyle \underbrace{\frac{\partial \mathbf{u}}{\partial t}}_{\textrm{Rate of change}} = \underbrace{\mathbf{D}\mathcal{L}\mathbf{u}}_{\textrm{Transport}}+ \underbrace{\mathbf{F}(\mathbf{u})}_{\textrm{Reactions}}, \quad \mathbf{u}(\mathbf{x},0) = \mathbf{u_0}(\mathbf{x}), \ \ \ \ \ (3)

where {\mathbf{D}\in \mathbb{R}^{m\times m}} is a transport coefficient matrix, {\mathcal{L}} is a linear spatial operator, {\mathbf{F}} is a nonlinear function of {\mathbf{u}}, and {\mathbf{u_0}} is an initial spatial distribution. We remark that {\mathcal{L}} should be understood to contain whatever boundary conditions one might have, and can in general depend on the spatial variable {\mathbf{x}}, as well as derivatives with respect to this variable (but not in any way on time {t}, as the non-autonomous setting is radically different).

A typical example which we will consider below would be the Laplacian on an interval, so that {\mathcal{L} = \nabla^2 = \partial_x^2} where {x \in [0, L]}. In this case we would also need to augment this operator with boundary conditions, e.g. the Neumann boundary conditions {\partial_x \mathbf{u}(0,t) = \partial_x \mathbf{u}(L,t) = 0}. Here we will think of this operator as really a scalar operator acting linearly on each component of the function {\mathbf{u}}; when this is not the case, and in particular when the spatial components of the system differ between equations, then one typically has trouble constructing eigenfunctions which are necessary to perform any kind of stability analysis, except in special cases. Finally we will also mention that in addition to restrictions of {\mathbf{F}}, one typically needs to consider elliptic operators {\mathcal{L}} and positive-definite matrices {\mathbf{D}} in order to ensure solutions to this system exist as time evolves. Blow-up of solutions and regularity are common mathematical questions that are studied for these kinds of systems. We make one further restriction on {\mathcal{L}} that it must annihilate constant vectors, so that a spatially homogeneous solution is determined purely by the function {\mathbf{F}}.

Next we will discuss stability of an equilibrium to this evolution equation against small spatial perturbations. As in the case of linear stability analysis for ordinary differential equations, perturbing an equilibrium solution will lead to a linear system, and for this reason we consider fundamental solutions to the linear operator {\mathcal{L}}. In the case of a compact domain and a suitable operator, these will be the set of eigenfunctions satisfying

\displaystyle \mathcal{L}w_k + \rho_k w_k = 0. \ \ \ \ \ (4)

Again we recall that for the Laplacian on an interval of length {L} with Neumann conditions, we have that {w_k(x) = \cos(k\pi x/L)} and {\rho_k = (k \pi/L)^2}, with {k=0, 1, 2, \dots}. In more complicated settings, such as {\mathcal{L}} being the Laplace-Beltrami operator on the sphere, or even the Laplacian on a planar square domain, there will be degeneracies such that many eigenfunctions correspond to given eigenvalues. Such degeneracies mean that analysis involving only the eigenvalues will not distinguish different eigenfunctions with the same eigenvalue. We can then assume that a solution to a linear equation involving {\mathcal{L}} can be expanded in a series of these functions as is typically done for separable linear partial differential equations.

With these assumptions in mind, we can then consider an equilibrium {\mathbf{u^*} \in \mathbb{R}^n} to the evolution equation (3) which is spatially homogeneous so that {\mathbf{F}(\mathbf{u^*}) = \mathbf{0}}. We assume the following ansatz for a perturbation to this equilibrium,

\displaystyle \mathbf{u} = \mathbf{u^*}+\varepsilon\mathbf{\hat{u}}(\mathbf{x},t), \ \ \ \ \ (5)

where {\varepsilon \ll 1}. Substituting (5) into (3) and using Taylor’s Theorem we find,

\displaystyle \varepsilon\frac{\partial \mathbf{\hat{u}}}{\partial t} = \varepsilon\mathbf{D}\mathcal{L}\mathbf{\hat{u}} + \mathbf{F}(\mathbf{u^*} + \varepsilon\mathbf{\hat{u}}) = \varepsilon\mathbf{D}\mathcal{L}\mathbf{\hat{u}} +\varepsilon\mathbf{J}\mathbf{\hat{u}} + O(\varepsilon^2), \ \ \ \ \ (6)

where {\mathbf{J}} is again the Jacobian of the nonlinear function {\mathbf{F}} evaluated at the steady state. We see that the the leading terms are satisfied by the equilibrium solution, so that up to corrections of order {\varepsilon}, the perturbation satisfies a linear equation.

Hence using the arguments regarding separability as before, we can expand our perturbation in the form {\mathbf{\hat{u}} = \exp(\lambda_k t)w_k(\mathbf{x})\mathbf{c}_k} and {\mathbf{c}_k} is an undetermined constant vector. By linearity, while we can consider a general perturbation (which is a sum of infinitely many of the eigenfunctions indexed by {k}), they will evolve separately from one another and so studying them in isolation at this point is sufficient. We also note that the exponential here arises naturally in the setting of linear first-order systems, and hence {\lambda} fundamentally determines the long-time behaviour of the perturbation.

Substituting this particular form of the ansatz into (6), and considering only terms of {O(\varepsilon)}, we find the linear algebraic system,

\displaystyle \mathbf{M}_k \mathbf{c}_k = (\lambda\mathbf{I}+\rho_k\mathbf{D}-\mathbf{J})\mathbf{c}_k=\mathbf{0}, \ \ \ \ \ (7)

where {\mathbf{I}} is the identity matrix. In order for this linear system to admit nontrivial solutions, we require that the dimension of the nullspace of this matrix is nonzero. From linear algebra, this says that we must have {\det(\mathbf{M}_k)=0}, which provides an {m}th degree polynomial relationship between {\lambda} and the eigenvalue {\rho_k}. One can then look for things like conditions such that {\Re(\lambda)>0} as {k} varies to find that specific modes will grow in time, and it is these that lead to instabilities which can become spatially structured states (as long as the modes which grow are spatially structured; one typically assumes that the eigenvalues of {\mathbf{J}} are all negative, and hence in the absence of the spatial operator, the homogeneous equilibrium is stable).

3. Linear Analysis of Reaction-Diffusion Systems: Turing instabilities

We now apply the results to the case of instabilities of the homogeneous equilibrium of a two-component reaction-diffusion system on the interval {\Omega = [0,L]} with Neumann boundary conditions. Writing {\mathbf{u} = (u,v)}, we consider the system,

\displaystyle \frac{\partial \mathbf{u}}{\partial t} = \mathbf{D}\nabla^2\mathbf{u}+ \mathbf{F}(\mathbf{u}), \ \ \ \ \ (8)

where the transport operator is now pure diffusion. We have the following quantities,

\displaystyle \mathbf{L} = \frac{\partial^2}{\partial x^2}, \quad \mathbf{F}(u,v) = \begin{pmatrix} f(u,v)\\ g(u,v) \end{pmatrix}, \quad \mathbf{J} = \begin{pmatrix} f_u &f_v \\ g_u & g_v \end{pmatrix}, \quad \mathbf{D} = \begin{pmatrix} d_u & 0\\ 0 & d_v \end{pmatrix}, \ \ \ \ \ (9)

where the subscripts of {f} and {g} denote partial derivatives, and we assume the diffusion coefficients are both positive. Again the eigenfunctions and eigenvalues are given by {w_k = \cos(k \pi x/L)} and {\rho_k = (k \pi/L)^2}. We can then immediately apply the determinant condition above to find a quadratic equation determining {\lambda} for each eigenfunction {k}, given by

\displaystyle \lambda^2+\lambda(\rho_k(d_u+d_v)-(f_u+g_v))+\rho_k^2d_ud_v+\rho_k(f_ud_v+g_vd_u)+f_ug_v-f_vg_u=0. \ \ \ \ \ (10)

While this expression seems complicated, it is easy enough to deduce requirements for a spatial instability. In particular, we note if the linear and constant term in this equation are both positive for some {k}, this is equivalent to there being no values where {\Re(\lambda)>0} and hence the homogeneous equilibrium is stable to perturbations of these modes. Similarly setting {\rho_k=0} (or equivalently, considering the system without diffusion) we can see that for a stable equilibrium in the homogeneous case, we require {f_u + g_v < 0} and {f_ug_v - f_vg_u > 0}, which are the standard conditions for stability of a negative trace and positive determinant of the Jacobian of a planar system. Hence we then see that for any {\rho_k > 0}, the second term can never become negative, and so instability is possible only when the last term becomes negative, Assuming a continuum of eigenvalues {\rho_k}, which is accurate for a sufficiently large domain, one can find algebraic conditions just in terms of model parameters, but here we will leave the eigenvalue {\rho_k} in the condition, as it typically corresponds to the geometry. Murray’s Mathematical Biology gives a far more detailed treatment of these conditions and much more.

More generally, one can replace the Laplacian above by the Laplace-Beltrami operator on a manifold, and show that the growth rates are exactly determined by the above expression. The major difference is that the eigenvalues and eigenfunctions can be much more difficult to compute, and depend on the geometry in interesting ways. For examples of this, as well as simulations of pattern-forming instabilities with other operators, see the videos I have made for growing manifolds and reaction-advection-diffusion equations.

4. Instabilities in Network Reaction-Diffusion Systems

Given an undirected weighted graph (that is, with {m} nodes or vertices {\mathbf{v}} and {\ell} edges {\mathbf{e}}), we can use the adjacency matrix {A_{ij}} to describe a diffusion process between these nodes. In fact, this will involve something called the graph laplacian, which we will denote by {\mathbf{L}}. We think of having {2} species at each node {v_i} reacting with one another given by {\mathbf{u_i}}, and diffusing to neighboring nodes. An analogue of the above reaction-diffusion system, given by (8), on a graph then looks like,

\displaystyle \frac{\textrm{d} \mathbf{u_i}}{\textrm{d} t} = \sum_{i=1}^{m}A_{ij}\mathbf{D}\left(\mathbf{u_j}-\mathbf{u_i}\right)+ \mathbf{F}(\mathbf{u_i}) = \sum_{i=1}^{m}L_{ij}\mathbf{D}\mathbf{u_i}+ \mathbf{F}(\mathbf{u_i}). \ \ \ \ \ (11)

As this is a system of {2m} equations, one could construct a Jacobian from the linearization about an equilibrium of the kinetics and simply compute the eigenvalues. Alternatively, as before, we can exploit the spectral properties of the matrix {\mathbf{L}} to find that perturbations to a homogeneous equilibrium will again have growth rates which satisfy Equation (10). The derivation of this fact replaces the eigenfunctions of {\nabla^2} above with the eigenvectors of {\mathbf{L}}, but is essentially unchanged.

This entry was posted in Uncategorized. Bookmark the permalink.

Leave a comment