Parameter identification through mode isolation for reaction – diffusion systems on arbitrary geometries

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 natural phenomena in areas such as developmental and cancer biology, cell motility and material science. In many of these applications, often one is interested in identifying parameters which will lead to a particular pattern for a given reaction–diffusion model. 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 examples and demonstrate that mode isolation is straightforward especially for low eigenvalues. Additionally, we show that in some cases the inhomogeneous steady state can be a linear combination of eigenfunctions. Finally,


Introduction
In his seminal work, Turing [59] presented an elegant mathematical theory of reaction-diffusion type for pattern formation in developmental biology.He showed that, via a symmetry breaking process, a homogeneous steady 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 natural phenomena including cancer invasion and angiogenesis in cancer biology [10,11,21], pattern formation in developmental biology [29,42], wound healing in biomedicine [13,52], cell motility [22,43,44] and material science [9,33], among many others.Despite their numerous applications, Turing's theory of pattern formation has been widely criticized 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 Sick et al. [53].
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 [22,36,41].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 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 paper provides a theoretical and computational framework to answer such a question.A recent paper by Dhillon et al. [14] uses a similar approach to model pattern development and presents a multiresolution algorithm for tracing bifurcation branches.
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 [30].The finite element method is well known for its capability to deal with complex irregular geometries [6,18,60].Alternative numerical methods such as finite differences [7], spectral methods [11,50] 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 [47].In such a setting our approach may provide parameters which serve as a suitable initial guess for a more advanced parameter identification algorithm [12,20].
The remainder of this paper is structured as follows.In Sec. 2, we introduce the mathematical model which we study in this work.We summarize the necessary and sufficient conditions for Turing diffusion-driven instability in Sec. 3. We then detail how mode selection and parameter identification are carried out.In Secs. 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 Sec. 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 (where no analytical solutions exist).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.

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 [47].Let Ω ⊂ R m (m = 1, 2, 3) be a simply connected bounded stationary volume for all time t ∈ I = [0, t F ], t F > 0 and ∂Ω be the surface boundary enclosing Ω.Also let u = (u(x, t), v(x, t)) T be a vector of two chemical concentrations at position x ∈ Ω ⊂ R m and time t ∈ I.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 [47,59] to yield the dimensional system where ∆ denotes the usual Cartesian Laplace operator, D u > 0 and D v > 0 are diffusion coefficients.Here, n is the unit outward normal to ∂Ω.Initial conditions are prescribed through non-negative bounded functions u 0 (x) and v 0 (x).In the above, f (u, v) and g(u, v) represent nonlinear reactions.
In the case of surfaces, the Laplace operator is replaced by the Laplace-Beltrami operator ∆ Γ , where Γ is the (smooth) surface.Surface gradients are also employed.This can be described as follows (for further details we refer the interested reader to see [16]).If f : Γ → R is differentiable at x ∈ Γ we can define the tangential gradient of f at x ∈ Γ by Here f is a smooth extension of f : Γ → R to an (n + 1)-dimensional neighborhood U of the surface Γ, so that f| Γ = f .∇ is the gradient in R n+1 and n is the unit normal.The Laplace-Beltrami operator applied to a twice differentiable function f ∈ C 2 (Γ) is then defined by 3) 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 behavior using linear stability theory and numerical methods.Linear stability analysis is one way of determining the behavior 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 [47,59].The mathematical treatment of the derivation of the necessary conditions for diffusion-driven instability requires solving the wellknown eigenvalue problem, with W a solution of where the solution pairs (k (eigenvalues), W k (x) (eigenfunctions) are obtained either analytically on certain spatial domains or numerically for the general case) of this equation can be compared to the spatially inhomogeneous steady state solutions of (2.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 paper "Can one hear the shape of the drum?" by Kac [31].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 et al. [25] 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 [11,36], 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 (2.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.

Conditions for Diffusion-Driven Instability for Reaction-Diffusion Systems
All two-component reaction-diffusion systems of the form (2.1) can be nondimensionalized and scaled to take the form where u = u(x, t), v = v(x, t), d is the ratio of diffusion coefficients, f (u, v) and g(u, v) describe the reaction kinetics.For simplicity, we assume that f and g are continuously differentiable, γ can be described as the relative strength of the reaction terms or alternatively as proportional to the domain size.We have zero flux boundary conditions (homogeneous Neumann) because we want only internal sources of instability, i.e. self-organization of the system.A uniform steady state (u s , v s ) is a fixed point where (u, v) = (u s , v s ), constant in time and space, satisfies (3.1), i.e. (u t , v t )| u=us,v=vs = 0. We can find the steady state by solving The conditions for instability due to diffusion are well known (see, for example, [47]).First, in the absence of diffusion, the steady state (u s , v s ) is linearly stable if and only if the partial derivatives of f and g at (u s , v s ) satisfy Linear stability analysis considering small perturbations from the equilibrium w(x, t) = (û(x, t), v(x, t)) leads us to the system which can be solved by method of separation of variables to yield where W k (x) solves the eigenvalue problem These are modes that will decay with time unless the wavenumber k 2 satisfies this means that instability will occur if where We exploit this range to isolate particular patterns/modes.The unstable modes will correspond to the eigenfunctions of the Laplacian (or Laplace-Beltrami operator) on the chosen domain or surface with the selected boundary conditions and k 2 the associated eigenvalues.The effect of varying d and γ on (3.6) is shown in Fig. 1.In summary the necessary conditions for diffusion driven instability are ) 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 [39,40,47].
1850053-6 Int.J. Biomath.2018.11.Downloaded from www.worldscientific.comby UNIVERSITY OF SUSSEX on 09/26/18.Re-use and distribution is strictly not permitted, except for Open Access articles.) has no roots so no modes can be isolated.As d increases, so does the difference between the two roots hence there is more chance that the value of k we seek will be between k 2 − and k 2 + .Similarly, for a fixed value of d, increasing γ causes both k 2 − and k 2 + to increase.

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

Schnakenberg or activator-depleted substrate kinetics
The Schnakenberg kinetics [51] are a condensed version of the well-documented Brusselator model describing a series of autocatalytic reactions also known as activator-depleted models [23,49], and these are characterized by Using the law of mass action and the non-dimensionalization of f and g, within system (3.1),we obtain that where a and b are positive parameters.

Gierer-Meinhart kinetics
One of the models proposed by Gierer and Meinhardt [23] describes a system whereby an "activator" activates the production of an "inhibitor" which inhibits the production of the activator.Again the non-dimensionalized form can be obtained where a and b are positive parameters (representing constant production rate and linear degradation rate, respectively) and k can be thought of as the saturation concentration of u.

Thomas kinetics
The Thomas model [58] is an immobilized-enzyme substrate-inhibition mechanism which can be written in non-dimensional form as where a, ρ, K, α, β are all non-negative parameters.This can be interpreted as in [46] by saying that u and v: (i) are generated by constant production a and αb, respectively, (ii) decay linearly proportional to u and αv, respectively and (iii) are used up in a substrate inhibition manner ρuv 1+u+Ku 2 .

Overview on Mode Isolation for Reaction-Diffusion Systems
The goal of mode isolation is to choose parameters, in our case (d, γ), 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 in one dimension, squares and triangles in [36].In [22], wavenumbers of a visco-elastic model are isolated on the unit disk.
We use similar ideas in this 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.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 and this is due to the spectral theorem for compact self-adjoint operators [8,32,57].
The aim is to have an algorithm to find the parameter values d and γ for a given eigenpair (k 2 , W ) such that only patterns analogous to W will grow.For this, one needs that the corresponding k is in the range defined in (3.8) where and that no other k lies in this range.In other words, the sign of the polynomial c(k 2 ) for a given k determines if the mode will grow.Figure 1 illustrates how the graph of c(k 2 ) changes as d and γ are varied.We define the critical diffusion ratio d c as the root of We find (k 2 , W ) either analytically or numerically.Then we propose the following algorithm described in pseudo-code: , γ > 0, f, g and the k l,n that we wish to be uniquely isolated.
( Note that we cannot have d < d c (because then c(k 2 ) would have no roots) nor γ < 0 (because k 2 > 0).

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 discretization and an implicit-explicit time-stepping scheme for the temporal approximation [35,37,50].In order to compute a finite element approximation, we write the weak formulation of (3.1) as follows: Find u, v ∈ L 2 (0, T ; H 1 (Ω)) such that for all φ, ψ ∈ H 1 (Ω) we have 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 (3.1) may be shown for example by the method of invariant regions proposed and analyzed by Smöller [54].

Spatial discretization
We define the computational domain Ω h by requiring that Ω h be a polyhedral approximation to Ω. Furthermore, we define T h to be a triangulation of Ω h made up of non-degenerate elements κ i , i.e.T h = i {κ i }.Finally, we define the finite element space where we use the Lagrange interpolant I h of the initial data into V h as initial conditions for the scheme.Letting N h be the total number of degrees of freedom of the nodes for the finite element discretization, we can write In order to illustrate a concrete example of the scheme, we focus on the reactiondiffusion system with Schnakenberg kinetics (3.11).The finite element approximation (5.2) with the Schnakenberg kinetics can be written in matrix-vector form as follows: (5.4a) where α = (α 1 , . . ., α N h ) and β = (β 1 , . . ., β N h ) are the coefficient vectors of the finite element functions U and V respectively and M and A are mass and stiffness matrices and H is a load vector with entries given by 1850053-10 Int.J. Biomath.2018.11.Downloaded from www.worldscientific.comby UNIVERSITY OF SUSSEX on 09/26/18.Re-use and distribution is strictly not permitted, except for Open Access articles.

Temporal discretization
For the temporal discretization we employ an IMEX method [35,37,50] 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, f n = f (t n ), the fully discrete scheme we employ reads, for n = 0, 1, . . ., where we use Lagrange interpolant of the initial data into V h as initial conditions for the scheme.This leads us to the following matrix vector form: Since we are interested in convergence to a spatially inhomogeneous steady state, for the stopping criteria we use the L 2 norm of the approximate time-derivative of the discrete solution, stopping the computation if this decreases below a tolerance, usually 10 −9 (see Fig. 2).

Numerical computations
We take the parameter values as shown in Table 1, and discretisation parameters as in Table 2.The uniform states for Schnakenberg kinetics were obtained analytically while for the Gierer-Meinhardt and Thomas reaction kinetics these were calculated computationally using the Newton-Raphson method [2,36].For the initial data we use small quasi-random perturbations around the uniform steady state values.The linear system (5.7) is solved using the conjugate gradient method [4,24,28].

Convergence to a steady state
Figure 2 plots the L 2 norm of the discrete time-derivative of U and V against the elapsed time.To begin with the L 2 norm is large.This quickly decays due to

(b).
There is an initial decay due to diffusion followed by a growth because of the exponentially growing modes which eventually decays, due to the dominant nonlinear terms.
diffusion followed by 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 bind the exponentially growing solution thereby giving rise to a spatially inhomogeneous steady state.

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 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: To approximate the solution we employ the finite element method for the spatial discretization outlined in Sec. 5. We work with the weak formulation of the eigenvalue problem and look for an approximate eigenpairs (W h , k 2 h ) ∈ V h × R + (where V h contains all continuous piecewise linear functions on a given mesh) such that As in (5.4) this may be written in matrix-vector form, we want to find (w, where A and M are stiffness and mass matrices defined in the same way as in Eq. (5.5).This is a generalized eigenvalue problem.We use the package deal.II [4] for its approximation using SLEPc and the Krylov-Schur algorithm.For completeness we give a description of the algorithm employed in Appendix A.

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 Fig. 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  .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 [4] and the references therein.

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 (2.4) on the sphere, we convert the eigenvalue problem into spherical coordinates.The eigenvalue problem in spherical coordinates is as follows [2,45]: with homogeneous Neumann boundary conditions.The solutions of the above eigenvalue problem are well known and are obtained using separation of variables [2,45].Following [2, pp. 424-428] we find an infinite number of discrete solutions of the form ,n r e imφ P m l (cos θ), (8.1)where l, m, n are all integers such that |m| ≤ l ≤ n, A m l,n are constants, J α (x) is a Bessel function of the first kind, i.e.
2 ) 2j+α with Γ(n) = (n−1)!, P m l (x) are the associated Legendre polynomials and j l+ 1 2 ,n are zeros of the differential of the spherical Bessel function.We can find the eigenvalues k 2 l,n = (j l+ 1 2 ,n ) 2 numerically (using the fact that It follows that for each eigenvalue λ l,n = k 2 l,n there are 2l + 1 possible eigenfunctions.Figure 4 shows the eigenfunctions for some selected values of l, m and n.For example k 1,1 = 2.08158 is the first zero of J The spherical Bessel function is given by J 1 = e imφ P m 1 (cos θ) are spherical harmonics whose real parts can be written in Cartesian coordinates as Since the system we are solving is invariant to polarity we can consider these to be equivalent.Figure 4(a) shows a plot of the eigenfunction where as usual r 2 = x 2 + y 2 + z 2 .The second example k 2,1 = 3.34209 corresponds to the eigenfunctions Choosing m = 0, converting the above to Cartesian coordinates and taking the real part gives The plot of the function w 0 21 is shown in Fig. 4(b).

Mode isolation on the sphere
Using the method described in Sec. 4 with all other parameters fixed as in Table 1 we can isolate the wavenumbers for the reaction-diffusion system with Schnakenberg kinetics and these are shown in Table 3.Similarly, for Thomas and Gierer-Meinhart (Table 4).In all these cases the interval [k − , k + ] is centered on k l,n .

Simulations of the reaction-diffusion systems on the unit sphere
Solving using deal.II we use the mesh shown in Fig. 3(a).The timestep is taken to be τ = 10 −3 .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 α 0 = 0.995 + 0.01 , β 0 = 0.895 + 0.01 , ( where is a uniformly distributed random variable between 0 and 1. 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 Figs. 5 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 [5].

Example 2: Ellipse
Eigenmodes on an ellipse have been investigated in various papers [19,26,48,61].Finding the solution involves numerically solving the Mathieu and modified  .Eigenfunctions corresponding to the labeled eigenvalues on an ellipse.These are solutions of (6.1) approximated using deal.II.
Mathieu equations [1].In particular Wu and Shivakumar [61] analytically find the first eigenvalue of ellipses with Dirichlet boundary conditions, of various sizes of ellipse.Using the eigenvalue solver described in Sec.6, with Dirichlet boundary conditions, we can reproduce their results (results not reported for brevity's sake).
In the following we consider Neumann boundary conditions and choose the semimajor axis to be twice the semi-minor axis.
The eigenvalues and eigenfunctions are shown in Fig. 7. Figure 8 shows the converged solutions of the reaction-diffusion system when the chosen values of d and γ isolate the corresponding wavenumbers k 2 i = λ i .It must be observed that the pattern computed will be a scalar multiple of the eigenfunction.This scalar may be negative which results in a reversed pattern (compare for example Figs.7(a

Example 3: Dumbbell
As a third example we consider the dumbbell shaped domain shown in Fig. 3(e).The solver for the eigenvalue problem on this mesh gives the output of eigenvalues and eigenfunctions shown in Fig. 9.The corresponding steady state solutions with the parameters obtained by mode isolation are shown in Fig. 10.

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 (6.1) and (3.1).To approximate solutions in this case, we employ the surface finite element method [6,[15][16][17][18]38].
The eigenpairs on the surface of the unit sphere can be found analytically and are well known and documented in [11] for example.The eigenfunctions are referred to as spherical harmonics.They are the restrictions of the eigenfunctions (8.1) to the surface.The eigenvalues are of the form k 2 = l(l + 1), where l is an integer, and the eigenfunctions are where m and P m l are as described in Sec.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 using 98306 degrees of freedom we obtain the following output of the first 30 eigenvalues computed to six significant figures:  As expected these are the first five values of the form k 2 = l(l + 1) with l = 1, 2, 3, 4, 5.It must be observed that the finite element method is known to be less effective for higher eigenvalues due to the min-max theorem [56].This means we must use a highly refined mesh in order to obtain values that are closer to the analytical values.The eigenfunctions are analogous to those detailed in Sec. 8.As the size of the sphere (radius R) increases, the eigenvalues decrease (specifically they are multiplied by 1 R 2 ).This effect is demonstrated in [34] where they show that fixing other values and increasing R causes higher mode patterns and hence more complex patterns are obtained.

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 Fig. 3(g), 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 Fig. 12. Various patterns are observed including stripes, spots and concentric rings.

Examples 6 and 7: "Eel " shapes
When computing on surfaces, one has to consider whether or not the surface has a boundary.In papers modeling fish or eel patterns (see, for example, [60]), 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, in Figs. 13 and 14, that the eigenvalues and eigenfunctions are very similar and it is possible to isolate similar patterns using the same parameter values.

Quantitative comparisons
By inspecting the plots it can be observed that the modes qualitatively appear to be isolated.To further expand on this, we normalize both the solutions and eigenfunctions to be in the range [−1, 1] then compute the L 2 norm of their difference and 1850053-24 Int.J. Biomath.2018.11.Downloaded from www.worldscientific.comby UNIVERSITY OF SUSSEX on 09/26/18.Re-use and distribution is strictly not permitted, except for Open Access articles.results of these computations are shown in Table 5. Results on the sphere are not possible due to rotational symmetry.It turns out that these L 2 norm differences are small and are due to a number of factors: First, the chosen numerical parameters: the differences get smaller and smaller with further grid refinement.On the other hand, numerical tests seem to suggest that refining the time-step does not make a significant difference in the decrease of the L 2 norms.Second, due to mode clustering, the L 2 norm differences can be affected by small contributions from nearby modes that are residing in the same excitable region.Last, the treatment of the nonlinear terms plays a significant role in the decrease of these L 2 norm differences.

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 (or more generally Laplace-Beltrami), 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 three 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 k 2 but can become slightly more difficult for higher values of k 2 .This is due to the approximation of the nonlinear terms and clustering of the eigenvalues of the linear problem.• When two or more eigenvalues are clustered close to each other it becomes difficult to isolate them computationally as well as analyltically.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 observation is only for the case of zero-flux boundary conditions.Imposing 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 [6,18,35,37,60].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.

(a) γ = 15 (b) d = 10 Fig. 1 .
Fig.1.Here the dispersal relation(3.6)  is plotted (for Schnakenberg kinetics).For a fixed value of γ, when d is below the critical value dc, c(k 2 ) has no roots so no modes can be isolated.As d increases, so does the difference between the two roots hence there is more chance that the value of k we seek will be between k 2 − and k 2 + .Similarly, for a fixed value of d, increasing γ causes both k 2 − and k 2 + to increase.

( 2 )
Compute the dispersal relation(3.6)  for the chosen reaction kinetics (this is independent of the geometry) and the range of admissible wave numbers as a function of d and γ.(3) Compute d * and γ * such that only one of the eigenvalues (wave numbers) computed in step (1) lies 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.

Fig. 2 .
Fig. 2. Plot of the L 2 norm of the discrete time-derivative over time for the example shown in Fig. 8(b).There is an initial decay due to diffusion followed by a growth because of the exponentially growing modes which eventually decays, due to the dominant nonlinear terms.

Fig. 3 .
Fig. 3. Examples of mesh generation for different volumes and surfaces: (a)-(c) Mesh generation on the unit sphere.(d) The ellipse which is a deformation of a circle mesh.(e)-(f) The dumbbell is a deformation of the bulk of a sphere.(g) The "fish" shape is a deformation of the surface of a sphere.(h) An "eel" is modeled by a cylinder with an open boundary and additionally as the same cylinder with added rounded boundaries.1850053-13 Int.J. Biomath.2018.11.Downloaded from www.worldscientific.comby UNIVERSITY OF SUSSEX on 09/26/18.Re-use and distribution is strictly not permitted, except for Open Access articles.

Fig. 3 .
Fig. 3. (Continued )a triangulation Ω h 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 (Γ h = Ω h | d Ω), 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[4] and the references therein.

1 Fig. 4 .
Fig. 4. Analytical solutions to the eigenvalue problem on the unit sphere i.e. (8.1) for selected values of l, m and n.These are plotted using deal.II.1850053-15 Int.J. Biomath.2018.11.Downloaded from www.worldscientific.comby UNIVERSITY OF SUSSEX on 09/26/18.Re-use and distribution is strictly not permitted, except for Open Access articles.

(a) λ 1 Fig. 7
Fig.7.Eigenfunctions corresponding to the labeled eigenvalues on an ellipse.These are solutions of (6.1) approximated using deal.II.

Fig. 8 .
Fig. 8. Converged solutions of system (3.1) with Schnakenberg kinetics(3.11), on an ellipse for the species u, they all match the associated eigenfunctions shown in Fig.7.It must be noted that the pattern can appear to be reversed (e.g. in (a) and (g)), this is due to the choice of the initial conditions.Choosing appropriate initial conditions results in a pattern similar to that shown in Fig.7(a) or 7(g).

Fig. 10 .
Fig. 10.Converged u solutions of system (3.1) with Schnakenberg kinetics (3.11) on a dumbbell.Eigenvalues λ 1 , λ 2 , λ 5 , λ 6 have been isolated, however since λ ≈ λ 4 in (c) we see a linear combination of their eigenfunctions.(e) Shows a reversed mode due to the nature of the random initial conditions as previously described in Sec.8.2.

Fig. 11 .
Fig. 11.Mode isolation for the reaction-diffusion system with Schnakenberg kinetics on the surface of the sphere.(a) The surface finite element solution with given parameters d = 9 and γ = 35 and (b) numerically computed eigenfunction corresponding to eigenvalue λ 9 = 12.0186.
1.3 restricted to the boundary.This shows that the eigenvalue solver gives the required output.Since the results are shown in Sec.8.1.3,we only show one example of mode isolation in Fig. 11.As mentioned in Sec. 3, γ can be thought of as being proportional to the domain size.Here we only consider a sphere with radius one.

Fig. 12 .
Fig. 12. Surface finite element solutions corresponding to the u species of the reaction-diffusion system with Schnakenberg kinetics with the given parameters on the left and numerically computed eigenfunctions corresponding to the given eigenvalue on the right.Again we observe reversed modes as described in Sec.8.2.1850053-23 Int.J. Biomath.2018.11.Downloaded from www.worldscientific.comby UNIVERSITY OF SUSSEX on 09/26/18.Re-use and distribution is strictly not permitted, except for Open Access articles.

Fig. 13 .Fig. 14 .
Fig. 13.Eigenfunctions of the Laplace-Beltrami operator on the "eel" shape with the corresponding eigenvalue.The left column shows the surface without a boundary and the right has a boundary.Note that, although the eigenfunctions are different, λ 23 ≈ λ 24 .
.4b) 1850053-4 This shifts the curve upward so the difference between k 2 If k l,n is uniquely isolated END.If not go to (3).
− increase γ by 1 (this number is arbitrary but should be small).This moves the curve to higher values of k.(3)If k 2 l,n < k 2 + decreaseγ by 1.This moves the curve to lower values of k. (4) If there exists another k * l,n = k l,n such that k 2 − < k * 2 l,n < k 2 + then decrease by d c /100.Output: The appropriate d, γ.

Table 1 .
Parameters for reaction kinetic models and the corresponding uniform steady states.Int.J. Biomath.2018.11.Downloaded from www.worldscientific.comby UNIVERSITY OF SUSSEX on 09/26/18.Re-use and distribution is strictly not permitted, except for Open Access articles.

Table 3 .
Given particular values of d and γ we obtain the values for k − and k + as well as the corresponding wavenumbers that are isolated on the sphere, for the reaction-diffusion system with Schnakenberg kinetics.Int.J. Biomath.2018.11.Downloaded from www.worldscientific.comby UNIVERSITY OF SUSSEX on 09/26/18.Re-use and distribution is strictly not permitted, except for Open Access articles.

Table 5 .
L 2 norm of difference between converged solution and the selected eigenfunction (U − ω k ) is found for the examples shown.