Mode Isolation of RDSs

# A computational approach for mode isolation for reaction-diffusion systems on arbitrary geometries

## Abstract.

In this article we present a computational framework for isolating spatial patterns arising in the steady states of reaction-diffusion systems. Such systems have been used to model many different phenomena in areas such as developmental and cancer biology, cell motility and material science. Often one is interested in identifying parameters which will lead to a particular pattern. To attempt to answer this, we compute eigenpairs of the Laplacian on a variety of domains and use linear stability analysis to determine parameter values for the system that will lead to spatially inhomogeneous steady states whose patterns correspond to particular eigenfunctions. This method has previously been used on domains and surfaces where the eigenvalues and eigenfunctions are found analytically in closed form. Our contribution to this methodology is that we numerically compute eigenpairs on arbitrary domains and surfaces. Here we present various examples and demonstrate that mode isolation is straightforward especially for low eigenvalues. Additionally we see that if two or more eigenvalues are in a permissible range then the inhomogeneous steady state can be a linear combination of the respective eigenfunctions. Finally we show an example which suggests that pattern formation is robust on similar surfaces in cases that the surface either has or does not have a boundary.

## 1. Introduction

In his seminal work, turing1952 presented an elegant mathematical theory of reaction-diffusion type for pattern formation in developmental biology. He showed that, via a symmetry breaking, a homogeneous state which is linearly stable in the absence of diffusion may be driven unstable in the presence of diffusion to give rise to the emergence of a spatially inhomogeneous pattern. This process is now well known as diffusion-driven instability or Turing instability. Since then, reaction-diffusion systems have been proposed and applied to model many phenomena including cancer invasion and angiogenesis in cancer biology (chap; chaplain; gatenby), pattern formation in developmental biology (hunding; maini), wound healing in biomedicine (dale; sherratt), cell motility (mogilner; moek; uduak) and material science (bozzini; krinsky) among many others. Despite their numerous applications, Turing’s theory of pattern formation has been widely criticised mainly due to the lack of robustness of the model system to changes in the parameters as well as the lack of experimental evidence of the existence of so-called morphogens with varying diffusivities. Only recently has the existence of chemical morphogens been experimentally validated in hair follicle pattern formation by sick2006.

To-date mode selection and parameter identification for reaction-diffusion systems have been mainly carried out on regular planar domains and surfaces where the eigenvalue problem can be analytically solved to yield analytical forms of the wave numbers as well as their corresponding eigenfunctions (ano; madzvamuse2003; uduak). In this work, we will depart from this framework and extend computationally mode selection and parameter identification to include arbitrary domains and stationary surfaces. First, we will solve the eigenvalue problem numerically using finite elements on planar domains or surface finite elements on smooth surfaces, respectively, to obtain the eigenmodes and their corresponding eigenfunctions. Here, we employ the Krylov-Schur algorithm (stewart) for solving the resulting algebraic system arising from the finite element discretisation. Second, we then pick an eigenmode to which we apply the necessary and sufficient conditions for Turing diffusion-driven instability in order to isolate reaction-kinetic model parameter values within a reaction-diffusion system. This process can be loosely thought of as an inverse problem for model parameter identification. Once the parameter values are isolated, the full reaction-diffusion system is then solved with these isolated parameter values to obtain an inhomogeneous spatially varying solution which is then compared to the numerically computed eigenfunction on the domain or surface. Alternatively, one could pose the following problem to which this methodology will provide insightful information which is otherwise out of reach with the current methodology: Given a biological pattern on a domain or surface and a plausible reaction-diffusion system, what are the model parameter values within this reaction-diffusion system that will give rise to the observed pattern? This article provides a theoretical and computational framework to answer such a question.

It must be observed that the eigenvalue problem and the reaction-diffusion system are both solved by a similar numerical method, the finite element method in multi-dimensions (johnson). The finite element method is well known for its capability to deal with complex irregular geometries (barreira; elliot; venkataraman). Alternative numerical methods such as finite differences (beckett), spectral methods (chap; ruuth) and finite volume methods among others could be used but with considerable efforts in dealing with geometrical complexities. As mentioned above one interpretation of our approach is that it provides a means of estimating parameter values such that the pattern predicted by linear stability analysis is close to a desired pattern. It must be noted that in many cases the steady state pattern may not be an eigenfunction (or a linear combination of the eigenfunctions) of the Laplacian on the given domain. This is since the nonlinear terms play a role in the resultant steady state pattern (murray2003). In such a setting our approach may provide parameters which serve as a suitable initial guess for a more advanced parameter identification algorithm (croft; garvie)

The remainder of this article is structured as follows. In Section 2 we introduce the mathematical model which we study in this work. We summarise the necessary and sufficient conditions for Turing diffusion-driven instability in Section 3. We then detail how mode selection and parameter identification are carried out. In Sections 4 and 6 we outline the new theoretical and computational framework for mode selection and parameter identification. The use of the finite element method is described in Section 5. We then give specific examples in 2- and 3-dimensions for regular (by which we mean domains on which analytic expressions for the eigenfunctions are available) as well as general domains and surfaces. We discuss the implications of our framework in the context of current methodologies and conclude that given a biological pattern and a reaction-diffusion system, our approach provides a useful tool for estimating parameter values which may give rise to the observed pattern.

## 2. Mathematical model framework

In order to illustrate with clarity the novelty of our approach, we first introduce the standard theoretical framework for reaction-diffusion systems in multi-dimensions (murray2003). Let be a simply connected bounded stationary volume for all time , and be the surface boundary enclosing . Also let be a vector of two chemical concentrations at position and time . The evolution equations for reaction-diffusion systems in the absence of cross-diffusion can be obtained from the application of the law of mass conservation and the extended Fick’s first law (murray2003; turing1952) to yield the dimensional system

 (1) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩{ut=DuΔu+f(u,v),\parvt=DvΔv+g(u,v),x∈Ω,t>0,\boldmathn⋅∇u=\boldmathn⋅∇v=0,xon∂Ω,t≥0,u(x,0)=u0(x),andv(x,0)=v0(x),xonΩ,t=0,

where denotes the usual cartesian Laplace operator, and are diffusion coefficients. Here, is the unit outward normal to . Initial conditions are prescribed through non-negative bounded functions and . In the above, and represent nonlinear reactions.

In the case of surfaces, the Laplace operator is replaced by the Laplace-Beltrami operator , where is the (smooth) surface. This can be described as follows (For more details we refer the interested reader to see dziukelliott). If is differentiable at we can define the tangential gradient of at by

 (2) ∇Γf=∇¯f−∇¯f⋅\boldmathn% \boldmathn.

Here is a smooth extension of to an -dimensional neighbourhood of the surface , so that . is the gradient in and is the unit normal. The Laplace-Beltrami operator applied to a twice differentiable function is given by

 (3) ΔΓf=∇Γ⋅∇Γf.

It must be observed that if the surface does not have a boundary, no boundary conditions are needed. If the surface has a boundary, we assume homogeneous Neumann boundary conditions.

Since the reaction terms are nonlinear, analytical solutions cannot normally be obtained. Therefore we investigate solution behaviour using linear stability theory and numerical methods. Linear stability analysis is one way of determining the behaviour of a nonlinear system near a given stationary point, normally a uniform steady state, of the given system. The idea is to find under what conditions on the nonlinear reaction kinetics is the uniform steady state linearly asymptotically stable in the absence of diffusion. When diffusion is introduced, the uniform steady state is driven unstable in what is now known as the process of diffusion-driven instability with the system converging to a spatially inhomogeneous steady state, thereby giving rise to patterning (murray2003; turing1952). The mathematical treatment of the derivation of the necessary conditions for diffusion-driven instability requires solving the well known eigenvalue problem, with a solution of

 (4a) ΔW+k2W=0,\boldmathx∈Ω, (4b) (\boldmathn⋅∇)W=0,\boldmathx∈∂Ω,

where the solution pairs ( (eigenvalues), (eigenfunctions) obtained either analytically on certain spatial domains or numerically for the general case) of this vector equation can be compared to the spatially inhomogeneous steady state solutions of (1), with good agreement expected near primary bifurcation points.

This approach is generally called mode isolation. The most famous exploration of this problem is the celebrated article ”Can one hear the shape of the drum?” by Mark Kac (kac). The question being asked is if one knows all the eigenvalues of the eigenvalue problem is it possible to determine the domain? It was later proven by Gordon, Webb and Wolpert (gww) that the answer is no and they gave examples of distinct regions with identical eigenvalues.

Other work concerned with mode isolation and linear stability theory for reaction-diffusion systems can be found in chap and ano, here the validation has been mainly restricted to special domains and volumes where the eigenvalue problem can be solved analytically. In this work we will depart from this framework, instead we will compute approximations of the eigenpairs on arbitrary, simply connected domains, volumes and surfaces. We then use these eigenvalues to calculate, by use of the Turing-parameter space restrictions, appropriate model parameter values. This approach can be thought to be analogous to an inverse parameter identification approach whereby, given the eigenvalues and eigenfunctions solving the eigenvalue problem (4), find model parameter values that would give rise to an inhomogeneous spatially varying solution similar to that exhibited by the eigenfunction. To confirm numerical predictions, we use the computed model parameter values to solve the full nonlinear reaction-diffusion systems and compare approximated eigenfunctions on these arbitrary domains, volumes and surfaces to the spatially inhomogeneous solutions obtained numerically.

To proceed, next we show the two-component form which we will work with and state the conditions for diffusion-driven instability. These will help us to isolate particular modes.

## 3. Conditions for diffusion driven instability for reaction-diffusion systems

All two component reaction-diffusion systems of the form (1) can be non-dimensionalised and scaled to take the form

 (5a) ut=γf(u,v)+Δu,vt=γg(u,v)+dΔv,\boldmathx∈Ω⊂Rn,t∈[0,∞], (5b) (\boldmathn⋅∇)(uv)=0\boldmathx∈∂Ωt∈[0,∞], (5c) u(\boldmathx,0),v(\boldmathx,0) % given,

where , is the ratio of diffusion coefficients, and describe the reaction kinetics. For simplicity, we assume that and are continuously differentiable, can be described as the relative strength of the reaction terms or alternatively as the domain size. We have zero flux boundary conditions (homogeneous Neumann) because we want only internal sources of instability, ie. self-organisation of the system. A uniform steady state is a fixed point where , constant in time and space, satisfies (5), i.e. . We can find the steady state by solving .
The conditions for instability due to diffusion are well known (see, for example murray2003). Firstly, in the absence of diffusion, the steady state is linearly stable if and only if the partial derivatives of and at satisfy

 (6) fu+gv<0 and fugv−fvgu>0.

Linear stability analysis considering small perturbations from the equilibrium leads us to the system

 (7) wt=γ(fufvgugv)w+(100d)Δw,

which can be solved by method of separation of variables to yield

 (8) w(\boldmathx,t)=∑kckeλtWk(% \boldmathx),

where solve the eigenvalue problem

 (9a) ΔW+k2W=0\boldmathx∈Ω, (9b) (\boldmathn⋅∇)W=0\boldmathx∈∂Ω.

These are modes that will decay with time unless the wavenumber satisfies

 (10) c(k2)=d(k2)2−γ(dfu+gv)k2+γ2(fugv−fvgu)<0,

this means that instability will occur if

 (11) dfu+gv>0,(dfu+gv)2−4d(fugv−fvgu)>0

and lies in the range where

 (12) k2±=γ(dfu+gv)±√(dfu+gv)2−4d(fugv−fvgu)2d.

We exploit this range to isolate particular patterns/modes. The unstable modes will correspond to the eigenfunctions of the Laplacian (or Laplace-Beltrami) on the chosen domain or surface with the selected boundary conditions and the associated eigenvalues. The effect of varying and on (10) is shown in Figure 1.
In summary the necessary conditions for diffusion driven instability are

 (13a) fu+gv <0, fugv−fvgu>0, (13b) dfu+gv >0, (dfu+gv)2−4d(fugv−fvgu)>0.

Additionally, the sufficient conditions for patterning formation are that one must be able to isolate distinct real wave numbers and that the domain must be large enough (madzvamuse2010; madzvamuse2015; murray2003).

### 3.1. Examples of reaction kinetics

For illustrative purposes, we consider three classical reaction kinetics as summarised below. The work presented in this article holds true for other similar reaction kinetics capable of generating Turing patterns.

#### Schnakenberg or activator-depleted substrate Kinetics

The Schnakenberg kinetics (schnak) are a condensed version of the well documented Brusselator model describing a series of autocatalytic reactions also known as activator-depleted models (giermein; prig), and can characterised by

 (14) A⇌XB+X→Y+D2X+Y→3X.

Using the Law of Mass Action and the non-dimensionalisation of and , within system (5), we obtain that

 (15) f(u,v)=a−u+u2v and g(u,v)=b−u2v,

where and are positive parameters.

#### Gierer-Meinhart Kinetics

One of the models proposed by giermein describes an system whereby an ”activator” activates the production of an ”inhibitor” which inhibits the production of the activator. Again the non-dimensionalised form can be obtained

 (16) f(u,v)=a−bu+u2v(1+ku2),andg(u,v)=u2−v,

where and are positive parameters (representing constant production rate and linear degradation respectively) and can be thought of as the saturation concentration of .

#### Thomas Kinetics

The Thomas model (thomas) is an immobilized-enzyme substrate-inhibition mechanism which can be written in non-dimensional form as

 (17) f(u,v)=a−u−ρuv1+u+Ku2,g(u,v)=αb−αv−ρuv1+u+Ku2,

where , , , , are all non-negative parameters. This can be interpreted as in murray1982 by saying that and

• are generated by constant production and respectively,

• decay linearly proportional to and respectively and

• are used up in a substrate inhibition manner .

## 4. Overview on mode isolation for reaction-diffusion systems

The goal of mode isolation is to choose parameters, in our case (), so that a trajectory starting from a small random perturbation from the steady state will evolve into a spatial pattern generated by one that corresponds, or at least is close to, a chosen eigenfunction of the Laplacian on that domain. Wavenumber isolation of reaction-diffusion systems is described by ano in one dimension, squares and triangles. In uduak wavenumbers of a visco-elastic model are isolated on the unit disk. We use similar ideas in the present work. The basic steps are as follows.

1. Determine a subset of eigenpairs of the Laplacian with suitable boundary conditions on the domain. For special domains this can be done analytically but in general must be done numerically.

2. Compute the dispersal relation (10) for the chosen reaction kinetics (this is independent of the geometry) and the range of admissible wave numbers as a function of and .

3. Compute and such that only one of the eigenvalues (wave numbers) computed in step 1 is in the range.

4. In order to compare with the patterned state, solve the reaction-diffusion system numerically with computed parameter values and compare with the numerically computed eigenfunctions.

It is possible to implement the above procedure simply because if a domain is bounded and the boundary is sufficiently regular, the Neumann Laplacian has a discrete spectrum of infinitely many non-negative eigenvalues with no finite accumulation point

 (18) 0<λ1≤λ2≤⋯,λn→∞

and this is due to the spectral theorem for compact self-adjoint operators (benguria, 2016; kreyszig; taylor).

The aim is to have an algorithm to find the parameter values and for a given eigenpair such that only patterns analogous to will grow. For this, one needs that the corresponding is in the range defined in (12)

 (19) γL=k2−

where

 (20a) L=(dfu+gv)−√(dfu+gv)2−4d(fugv−fvgu)2d, (20b) R=(dfu+gv)+√(dfu+gv)2−4d(fugv−fvgu)2d,

and that no other is in this range.

In other words, the sign of the polynomial for a given determines if the mode will grow. Figure 1 illustrates how the graph of changes as and are varied. We define the critical diffusion ratio as the root of

 (21) d2cf2u+2(2fvgu−fugv)dc+g2v=0.

We find either analytically or numerically. Then we propose the following algorithm described in pseudo-code:

Input: , , , and the that we wish to be uniquely isolated.

1. Compute and from (19).

2. If increase by 1 (this number is arbitrary but should be small). This moves the curve to higher values of .

3. If decrease by 1. This moves the curve to lower values of .

4. If there exists another such that then decrease by . This shifts the curve upwards so the difference between and is smaller.

5. If is uniquely isolated END. If not go to 3.

Output: The appropriate .
Note that we cannot have (because then would have no roots) nor (because ).

## 5. Finite element method for reaction diffusion systems

In order to validate that our mode isolation algorithm does indeed isolate the desired unstable mode, we will simulate the reaction-diffusion systems under consideration with the computed parameter values. To do this we employ a finite element method for the space discretisation and an implicit-explicit time-stepping scheme for the temporal approximation (lakkis2013; madzvamuse2006; ruuth).

In order to compute a finite element approximation, we write the weak formulation of (5) as follows: Find such that for all we have

 (22) {∫Ωutϕ+∫Ω∇u⋅∇ϕ=γ∫Ωf(u,v)ϕ,∫Ωvtϕ+d∫Ω∇v⋅∇ϕ=γ∫Ωg(u,v)ϕ,\boldmathx∈Ω,t>0.

In this work we shall assume the well posedness of the weak formulation above. We note that for suitable parameter values existence and uniqueness of a classical solution, and hence a weak solution, to (5) may be shown for example by the method of invariant regions proposed and analysed by Smöller (smoller).

### 5.1. Spatial discretisation

We define the computational domain by requiring that is a polyhedral approximation to . We define to be a triangulation of made up of non-degenerate elements , i.e., . We define the finite element space . The semidiscrete (space discrete) finite element approximation to (22) seeks a pair such that

 (23) {∫ΩhUtϕ+∫Ωh∇U⋅∇ϕ=γ∫ΩhIh[f(U,V)]ϕ,∫ΩhVtϕ+d∫Ωh∇V⋅∇ϕ=γ∫ΩhIh[g(U,V)]ϕ,∀ϕ∈Vh,

where we use the Lagrange interpolant of the initial data into as initial conditions for the scheme. In order to illustrate a concrete example of the scheme, we focus on the reaction-diffusion system with Schnakenberg kinetics (15). The finite element approximation (23) with the Schnakenberg kinetics can be written in matrix-vector form as follows

 (24a) Mαt+Aα=γ[aH−Mα+M(α)2β], (24b) Mβt+dAβ=γ[bH−M(α)2β],

where and are the coefficient vectors of the finite element functions and respectively and

 Mi,j=∫Ωhϕiϕj,Ai,j=∫Ωh∇ϕi⋅∇ϕjandHj=∫Ωhϕj.

### 5.2. Temporal discretisation

For the temporal discretisation we employ an IMEX method (lakkis2013; madzvamuse2006; ruuth) in which the diffusive term is treated implicitly and the reaction terms are treated explicitly, for simplicity we employ a uniform timestep . Introducing the shorthand for a time discrete sequence of functions, , the fully discrete scheme we employ reads, for , given find such that, ,

 (25) {∫Ωh1τ(Un+1−Un)ϕ+∫Ωh∇Un+1⋅∇ϕ=γ∫ΩhIh[f(Un,Vn)]ϕ,∫Ωh1τ(Vn+1−Vn)ϕ+d∫Ωh∇Vn+1⋅∇ϕ=γ∫ΩhIh[g(Un,Vn)]ϕ,

where we use Lagrange interpolant of the initial data into as initial conditions for the scheme. This leads us to the following matrix vector form

 (26a) (1τM+A)αm+1=γ[aH−Mαm+M(αm)2βm]+1τMαm, (26b) (1τM+dA)βm+1=γ[bH−M(αm)2βm]+1τMβm.

Since we are interested in convergence to a spatially inhomogeneous steady state, for the stopping criteria we use the norm of the approximate time derivative of the discrete solution, stopping the computation if this decreases below some tolerance (see Figure 2).

### 5.3. Numerical computations

We take the parameter values as shown in Table 1, for the initial data we use small quasi-random perturbations around the uniform steady state values. The linear system (26) is solved using the conjugate gradient method (dealii; golub; CG).

### 5.4. Convergence to a steady state

Figure 2 plots the norm of the discrete time derivative of and against the elapsed time. To begin with the difference is large. This quickly decays due to diffusion then there is a rapid growth, because of the exponentially growing modes. The time derivative eventually starts to decay due to the effects of the nonlinear terms that act to bound the exponentially growing solution thereby giving rise to a spatially inhomogeneous steady state.

## 6. Isolating modes on general domains

On arbitrary domains, analytical solutions for the eigenvalue problem are not typically available but approximate eigenpairs can be computed numerically. Numerically approximating these pairs is a significant challenge. In general, as we are only typically interested in a small number of eigenpairs, it is not necessary to find all solution pairs, however for our approach to mode isolation to remain applicable, it is important that we obtain consecutive pairs.
As previously stated, the eigenvalue problem we wish to solve is as follows,

 (27) {ΔW+k2W=0,\boldmathx∈Ω,(\boldmathn⋅∇)W=0,\boldmathx∈∂Ω.

To approximate the solution we employ the finite element method for the spatial discretisation outlined in Section 5. We work with the weak formulation of the eigenvalue problem and look for an approximate eigenpairs (where contains all continuous piecewise linear functions on a given mesh) such that

 (28) ∫Ω∇Wh⋅∇ϕ=k2∫ΩWh⋅ϕ,∀ϕ∈Vh.

As in (24) this may be written in matrix-vector form, we want to find , where is the dimension of such that

 (29) Aα=k2Mα,

where and are stiffness and mass matrices defined respectively, by

 (30) Ai,j=∫Ωh∇ϕi⋅∇ϕjandMi,j=∫Ωhϕiϕj.

This is a generalised eigenvalue problem. We use the package deal.II (dealii) for its approximation using SLEPc and the Krylov-Schur algorithm. For completeness we give a description of the algorithm employed in Appendix A.

## 7. Mesh generation

All the mesh generation is carried out using the deal.II library. We use hexahedral meshes for the volumes and quadrilaterals for the ellipse and surfaces. In Figure 3 we exhibit different meshes generated by this package on which we will carry out computations. We also consider smooth surfaces; these meshes are generated by creating a triangulation of the bulk of the domain then the surface triangulation is defined by collecting the faces of the elements of the bulk triangulation that lie on the surface (), i.e., the surface mesh is the trace of the volume mesh (in the example of the cylinder with open ends we use only the elements on the curved surface). For this reason the equations are not being approximated on the actual surface but on an approximation of it. For more details on surface mesh generation the reader is referred to dealii and the references therein.

## 8. Comparisons of eigenfunctions and spatially inhomogeneous steady states

### 8.1. Example 1: Sphere

We start by considering the unit sphere, a domain for which the eigenvalue problem can be solved analytically.

#### Eigenvalues and eigenfunctions of the Laplacian in the bulk of the unit sphere

In order to solve (4) on the sphere, we convert the eigenvalue problem into spherical coordinates. The eigenvalue problem in spherical coordinates is as follows (arfken; morimoto),

 Δw+k2w=1r2∂∂r(r2∂w∂r)+1r2sinθ∂∂θ(sinθ∂w∂θ)+1r2sinθ∂2w∂ϕ2+k2w=0,

with homogeneous Neumann boundary conditions. The solutions of the above eigenvalue problem are well known and are obtained using separation of variables (arfken; morimoto). Following arfken (p. 424-428) we find an infinite number of solutions of the form

 wml,n(r,θ,ϕ)=Aml,nJl+12(j′l+12,nr)eimϕPml(cosθ), where ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩l,m,n all integers such that |m|≤l≤n,Aml,n are constants,Jα(x)=∑∞j=0(−1)jj!Γ(1+j+α)(x2)2j+α with Γ(n)=(n−1)! (i.e. a Bessel function of the first kind% ),Pml(x) are associated Legendre polynomials,j′l+12,n are zeros of the differential of the % spherical Bessel function.

We can find the eigenvalues numerically (using the fact that ). It follows that for each eigenvalue there are possible eigenfunctions. Figure 4 shows the eigenfunctions for some selected values of , and . For example is the first zero of and corresponds to the eigenfunctions

 wm1,1(r,θ,ϕ)=J32(k1,1r)eimϕPm1(cosθ), with m=−1,0,1.

The spherical Bessel function is given by . Meanwhile are spherical harmonics whose real parts can be written in cartesian coordinates as , and . Since the system we are solving is not sensitive to polarity we can consider these to be equivalent. Figure 4LABEL:sub@fig:sphef1 shows a plot of the eigenfunction

 w11,1=(sin(k1,1r)(k1,1r)2−sin(k1,1r)(k1,1r)2)⋅xr,

where as usual . The second example, corresponds to the eigenfunctions

 wm2,1(r,θ,ϕ)=J52(k2,1r)eimϕPm2(cosθ), with −l≤m≤l.

Choosing , converting the above to cartesian coordinates and taking the real part gives

 w02,1(x,y,z)=((3k22,1r2−1)sin(k2,1r)k2,1r−3cos(k2,1r)k22,1r2)(14√5π⋅−x2−y2+2z2r2).

The function is plotted in Figure 4LABEL:sub@fig:sphef2.

#### Mode isolation on the sphere

Using the method described in Section 4 and the values given in Table 1 we can isolate the wavenumbers for the reaction-diffusion system with Schnakenberg kinetics and these are shown in Table 2. We can do the same for Thomas and Gierer-Meinhart (Table 3). In these cases the interval is centered on .

#### Simulations of the reaction-diffusion systems on the unit sphere

Solving using deal.II we use the mesh shown in Figure 3LABEL:sub@fig:sphmesh. The timestep is taken to be . We take the initial conditions to be a small random perturbation from the previously computed homogeneous steady state. So for the reaction-diffusion system with Schnakenberg kinetics, at each point in the grid we set the initial conditions to be:

 (32) α0=0.995+0.01ϵ,β0=0.895+0.01ϵ,

where is a uniformly distributed random variable between and .
For each eigenvalue there are a number of different eigenfunctions. Computing using the values obtained with mode isolation, the solution converges to either one of the eigenfunctions or a linear combination. These converged solutions are shown in Figures 6 and 6. It is possible to force the solution to converge to an eigenfunction (which it does not appear to with random initial perturbation) by making a suitable choice of initial condition, for example a perturbation of the desired eigenfunction, suitably scaled. Hence, in the case where multiple wave numbers are excited, pattern selection is heavily influenced by the choice of initial conditions which act as the basin of attraction, one of the major criticisms of Turing’s theory for pattern formation (BardandLauder1974).

### 8.2. Example 2: Ellipse

Eigenmodes on an ellipse have been investigated in various articles (fox; greben; neves; wu). Finding the solution involves numerically solving the Mathieu and modified Mathieu equations (hbk). In particular wu analytically find the first eigenvalue of ellipses with Dirichlet boundary conditions, of various sizes of ellipse. Using the eigenvalue solver described in Section 6, with Dirichlet boundary conditions, we can reproduce their results (results not reported in the interests of brevity). In the following we consider Neumann conditions and choose the semimajor axis to be twice the semiminor axis. The eigenvalues and eigenfunctions are shown in Figure 8. Figure 8 shows the converged solutions of the reaction diffusion system when the chosen values of and isolate the corresponding wavenumbers .

### 8.3. Example 3: Dumbbell

As a third example we consider the dumbbell shaped domain shown in Figure 3LABEL:sub@fig:dumbbellgrid. The solver for the eigenvalue problem on this mesh gives the output of eigenvalues and eigenfunctions shown in Figure 10. The corresponding steady state solution with the parameters obtained by mode isolation are shown in Figure 10.

### 8.4. Example 4: Surface of a sphere

In all the previous examples we considered bulk, volumetric domains. In this example we have a curved surface as the domain. This means using the Laplace Beltrami operator instead of the Laplacian in (27) and (5). To approximate solutions in this case, we employ the surface finite element method (barreira; dziuk; dziukelliott; elliot2015; elliot; chung).
The eigenpairs on the surface of the unit sphere can be found analytically and are well known and documented in chap for example. The eigenfunctions are referred to as spherical harmonics. They are the restriction of the eigenfunctions (31) to the surface. The eigenvalues are of the form , where is an integer, and the eigenfunctions are

 (33) wml(θ,ϕ)=AmleimϕPml(cosθ),

where and are as in Section 8.1.3. Therefore we can test the performance of the eigenvalue problem solver with this example. Using the eigenvalue solver on an approximated mesh of the surface of the sphere we obtain the following output of the first 30 eigenvalues computed to 4 decimal places

As expected these are the first 5 values of the form with . The values are not exact because the mesh is an approximation of the actual surface of the sphere. The eigenfunctions are analogous to those detailed in Section 8.1.3 restricted to the boundary. This shows that the eigenvalue solver gives the required output. Since the results are shown in Section 8.1.3 we only show one example of mode isolation in Figure 11.

### 8.5. Example 5: ”fish” surface

We now consider a smooth surface on which no analytical expression for the eigenpairs is available, the surface is taken to be diffeomorphic to the sphere and is shown in Figure 3LABEL:sub@fig:smmesh, it is meant to (very loosely) mimic the shape of a fish. We found the first 100 eigenpairs then chose several to isolate. These are shown in Figure 12. Various patterns are observed including stripes, spots and concentric rings.

### 8.6. Example 6 and 7 ”eel” shapes

When computing on surfaces, one has to consider whether or not the surface has a boundary. In papers modelling fish or eel patterns (see for example venkataraman), a surface with a boundary is often used. To investigate whether having a boundary is significant in this example we consider a surface with and without boundary. We see that the eigenvalues and eigenfunctions are very similar and it is possible to isolate similar patterns using the same parameter values.

## 9. Conclusion and further challenges

In this paper we have considered reaction-diffusion systems and have presented a framework for isolating particular spatially inhomogeneous patterns. The method involves finding eigenpairs of the Laplacian and computing parameters such that when the reaction-diffusion system is solved numerically, only patterns analogous to a particular eigenfunction will grow. In previous works the eigenvalue problem is solved analytically whereas in this paper both the eigenvalue problem and the reaction-diffusion system are solved using the finite element method. Advances in numerical software mean that we can find 100 eigenpairs in a few minutes and we have demonstrated that these eigenpairs match analytical results. The approach is shown to work for 3 different examples of nonlinear reaction kinetics and on a variety of domains and surfaces. In summary, the main observations are:

• Mode isolation is straightforward for low values of but can become slightly more difficult for higher values of . This is due to the approximation of the nonlinear terms and clustering of the eigenvalues of a linear problem.

• When two or more eigenvalues are clustered close to each other it becomes difficult to isolate them computationally. If two or more eigenvalues are in the permissible range then the inhomogeneous steady state could be a linear combination of the corresponding eigenfunctions.

• We display an example of two surfaces where pattern formation appears to be robust despite the fact one has a boundary while the other does not. An interesting investigation would be to see if this can be true for other geometries. Note that this is only the case for zero-flux boundary conditions. Imposing Dirichlet or Robin-type boundary conditions would result in substantially different patterns.

In this paper we have only considered stationary domains/volumes and surfaces. However the domains of biological processes generally evolve with time (barreira; elliot; lakkis2013; madzvamuse2006; venkataraman). This adds more complexity to solving the reaction-diffusion systems. An interesting and natural extension of this work would be to introduce domain growth and surface evolution. For this extension, studies on the effects of initial conditions would also be worthwhile.

## Data management

All the computational data output is included in the present manuscript.

## Acknowledgements

This work (LM) was supported by an EPSRC Doctoral Training Centre Studentship through the University of Sussex. CV and AM acknowledge support from the Leverhulme Trust Research Project Grant (RPG-2014-149) and the EPSRC grant (EP/J016780/1). This research was partly undertaken whilst LM, CV and AM were participants in the Isaac Newton Institute Program, Coupling Geometric PDEs with Physics for Cell Morphology, Motility and Pattern Formation. This work (AM) has received funding from the European Union Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement (No 642866). AM was partially supported by a grant from the Simons Foundation. LM acknowledges the support from the University of Sussex ITS for computational purposes.

## Appendix A The Krylov-Schur algorithm

The Krylov-Schur algorithm was introduced by stewart and is an alternative to the method of arnoldi. The aim of the algorithm is to compute a number of eigenpairs of a given square matrix A.
The basic Arnoldi algorithm has input matrix and initial vector of norm 1 ( will make up the columns of an matrix ) and output such that

 (34) AVm=VmHm+fe∗m,β=||f||2.

A Krylov decompostion is a generalised version of this and is given by

 (35) AVm=VmBm+vm+1b∗m+1,

where is not necessarily upper Hessenberg and is an arbitrary vector. The Krylov-Schur method is described in the SLEPc Technical Report [str7] as follows \blockquote Input: Matrix , initial vector , and dimension of the subspace
Output: A partial Schur decompostion

• Normalize

• Initialize

• Restart loop

• Perform steps of Arnoldi with deflation

• Reduce (part of the output of the Arnoldi algorithm) to (quasi-)triangular form,

• Sort the or diagonal blocks:

• Compute eigenpairs of ,

• Compute residual norm estimates,

• Exit if enough converged eigenpairs, otherwise lock newly converged vectors

• Choose ( (the number of currently converged eigenpairs) ) and set

• Compute (the leading subvector of ) and insert in the appropriate positions of

• end

If the eigenpairs of (ie solutions of ) are then the approximate eigenvalues of are and eigenvectors are . In our problem we have (29) (the generalised eigenvalue problem) instead of , and here one works with a spectral transformation or instead of .

121664