Bifurcation Analysis of Reaction Diffusion Systems on Arbitrary Surfaces
Abstract
In this paper we present computational techniques to investigate the solutions of twocomponent, nonlinear reactiondiffusion (RD) systems on arbitrary surfaces. We build on standard techniques for linear and nonlinear analysis of RD systems, and extend them to operate on largescale meshes for arbitrary surfaces. In particular, we use spectral techniques for a linear stability analysis to characterize and directly compose patterns emerging from homogeneities. We develop an implementation using surface finite element methods and a numerical eigenanalysis of the LaplaceBeltrami operator on surface meshes. In addition, we describe a technique to explore solutions of the nonlinear RD equations using numerical continuation. Here, we present a multiresolution approach that allows us to trace solution branches of the nonlinear equations efficiently even for largescale meshes. Finally, we demonstrate the working of our framework for two RD systems with applications in biological pattern formation: a Brusselator model that has been used to model pattern development on growing plant tips, and a chemotactic model for the formation of skin pigmentation patterns. While these models have been used previously on simple geometries, our framework allows us to study the impact of arbitrary geometries on emerging patterns.
∎ \@input@sec.aux \NAT@set@cites
1 Introduction
Reactiondiffusion (RD) is often used to model the development of biological systems, most prominently in the study of biological pattern formation. The mathematical representation of these models results in systems of nonlinear partial differential equations (PDEs). The analysis of these systems of PDEs aim at answering two fundamental questions: (i) what are the possible solutions that satisfy the given system of PDEs, and how can these solutions be discovered systematically; and (ii) which of these solutions are stable against minor perturbations. Two strategies are commonly used to perform this analysis. First, linear stability analysis can predict the emergence of new patterns near trivial, homogeneous solutions (homogeneities) of the PDEs. These emergent patterns correspond to sudden qualitative changes to the state of an RD system and, hence, constitute bifurcations from the homogeneity. Second, nonlinear analysis provides solutions of the nonlinear RD equations far away from the homogeneous steadystates. To construct these solutions, numerical continuation techniques are used to follow continuous branches of solutions starting from initial bifurcation patterns constructed using the linear analysis. Solutions vary gradually with one of the system parameters, called the continuation parameter, along each branch.
Several existing tools allow to delineate solutions for an RD system through linear and nonlinear bifurcation analyses. However, many of them are constrained to work with only simple surface geometries such as rectangles or hemispheres, and at low resolutions. These are serious limitations given that surface geometry plays an important role in pattern formation (Murray, 2003, page 108), that most of the interesting biological domains for RD systems have rather arbitrary shapes. In addition, corresponding patterns are too complex to be resolved with lowresolution meshes.
In this paper, we develop a framework to perform bifurcation analysis for generic RD systems with two components, with or without cross diffusion, acting on arbitrary surfaces. Unlike approaches that iteratively traverse the trivial branch for a desired continuation parameter to detect bifurcation points stepbystep, our proposed framework uses an analysissynthesis approach to directly determine bifurcation points and construct emerging patterns along the trivial branch. We exploit the Hermitian nature of the LaplaceBeltrami (LB) operator, acting on a given arbitrary surface, which enables the computation of a spectral basis for emerging patterns. This allows us to derive formulae to directly compose emergent bifurcation patterns from eigenfunctions (also called eigenmodes or wavemodes) of the LB operator. Similarly, a bifurcation point for a given bifurcation pattern can be computed directly using its eigenvalues. We discuss several boundary conditions for our framework that are common for biological systems and ensure that the LB operator is Hermitian. Unlike detectionbased approaches, our analysissynthesis approach avoids missing out on bifurcations due to potential failures of a test function. In addition, our approach allows for tracing multiple and mixedmode bifurcations apart from simple bifurcations.
Our framework also supports highresolution meshes for triangulated surface domains. We propose a progressive geometric multigrid approach to support multiple levels of mesh resolution. We first trace a branch at the lowest resolution for the surface mesh. Next, we use a simple twostep approach to upsample a solution pattern to the next level by first improving the mesh resolution for a solution, and then resolving the solution using an improved mesh geometry. This twostep approach progressively upsamples and converges a given pattern across multiple levels. Our multiresolution approach decouples branch tracing complexity from the complexity of dealing with a largescale system. Unlike traditional multigrid approaches, our twostep approach does not go back and forth across multiple levels. This lends a high degree of parallelisability to our framework.
We demonstrate the working of our framework for a Brusselator system with zero Dirichlet boundary conditions to study the emergent patterns of cotyledons on a conifer tip (Nagata et al, 2003). We also demonstrate its working with Murray’s chemotactic model for pattern formation of skin pigmentation (Murray and Myerscough, 1991), subject to zero Neumann boundary conditions. In both cases, we illustrate the influence of surface shape on pattern formation using several example geometries. Finally, we evaluate the computational performance of our multiresolution branch tracing approach. In summary, our contributions include:

A framework for analysing twocomponent RD systems with or without cross diffusion that supports arbitrary triangulated surface domains.

A direct analysissynthesis approach for computing emergent patterns and locating bifurcation points along the trivial branch based on spectral analysis of the LaplaceBeltrami operator on arbitrary triangulated surface domains.

A progressive geometric multigrid approach supporting highresolution FEM discretisation to determine patterns along nonlinear branches.

Two case studies illustrating the effect of arbitrary geometries on pattern formation.
2 Related Work
Emergent Patterns near Homogeneity.
Linear stability analysis with twocomponent RD systems is often employed to study the emergence of new patterns from near homogeneous patternless initial conditions. Nagata et al (2003, 2013) study the relation between the shape and size of a conifer embryo and the emergent cotyledon patterns. They use a Brusselator RD system acting on a parametric family of spherical caps while imposing zero Dirichlet boundary conditions near homogeneity. With this, they explain how the number of emergent cotyledons is simply the selection of a spherical cap harmonic based on the radius of the conifer and its curvature (Nagata et al, 2013). Winters et al (1990) explain emergence of heterogeneous snake skin color patterns with a chemotactic RD model with cross diffusion (i.e., where the flux of one component is driven by gradients in the concentration of the second component). They perform numerical simulations on flat rectangular domains and illustrate the similarity of emergent patterns to the patterns observed in nature, for different snake species. Winters et al. use zero Neumann boundary conditions for their RD system. Similarly, Gambino et al (2013a) discuss pattern formation due to crossdiffusion for LotkaVolterra kinetics between two components in a 2D rectangular domain, commonly used to model predatorprey populations. Kealy and Wollkind (2012) use linear stability analysis to study onset of various spatial Turing patterns for vegetation in an arid flat land using a twocomponent RD system. They also apply a weak nonlinear stability analysis to predict the longterm behaviour of these emerging vegetation patterns. Murray (2003, Chapter ) demonstrates the effect of both geometry and scale on the emergence of patterns under a twocomponent RD system and discusses the relevance of these parameters for explaining animal coat patterns. He derives an analytical form for the stripe and spot patterns that emerge on a tapering cylinder representing an animal tail. Also, he presents the selection of different stripe or patchy patterns (modes) with changes in the size of a planar 2D shape representing an animal coat. Most studies on the emergence of patterns, such as those discussed above, are limited to simple, welldefined surface geometries with analytically defined emergent patterns. Recently, Tuncer et al (2015) have introduced a projected Finite Elements Method for studying pattern formation by RD systems on surfaces that can be approximated analytically and later mapped with Lipschitz continuity onto a sphere. They note that studying RD systems on arbitrary surfaces is rather a “young and emerging research area”— (Tuncer et al, 2015) and that surface geometry is crucial for such studies. Our framework extends studies of RD systems with or without crossdiffusion to arbitrary surface domains without any geometric constraints by directly (numerically) computing emergent patterns on them. Also, it supports several common boundary conditions that arise in biological problems such as homogeneous Dirichlet, Neumann and Robin boundary conditions. As noted earlier, (arbitrary) surface geometry plays an important role in pattern formation. Thus our framework serves as an important tool for studying emergent patterns on ’real’ geometries.
Marginal Stability Analysis.
Marginal stability analysis is often used to study the interaction and mutualexclusivity of two or more emergent wavemodes. It is also used to demarcate and characterise the parameter space of an RD system. Kealy and Wollkind (2012) use analytically defined marginal stability curves to demarcate regions in the parameter space with subjectively different vegetation patterns on arid flat lands (represented as 2D rectangles). Nagata et al (2013) define marginal stability curves for emergent wavemodes in terms of their corresponding eigenvalues and investigate the influence of the spherical cap surface geometry on the pattern of cotyledons development. In particular, they note that changes in the curvature or size of a plant tip may cause a change in the number of cotyledons that develop despite that the concentrations of chemical precursors are fixed. Our framework generalises such marginal stability analysis to arbitrary domains. It numerically computes the eigenvalues for the wavemodes that constitute emergent patterns. These eigenvalues can then be used to plot corresponding marginal stability curves. In its full potential, our framework supports case studies with dynamically changing arbitrary shapes for marginal stability analysis as demonstrated later.
Bifurcations and Branch Tracing.
Analytical solutions for branch tracing are only possible for simple surface domains. Ma and Hu (2014) express branches and patterns for a twocomponent Brusselator model acting on a straight line domain in analytical forms. They prove that, except for the first branch along the continuation parameter dimension, all other branches are unstable. Méndez and Campos (2008) derive analytical expressions for tracing a branch with a single component RD system to predict the survival of an isolated patch of a population in its surrounding hostile environment. Using stability predictions along the branch, they establish that the survival of a population at a very low or negative growth rate depends on its initial density. Instead, we support branch tracing with numerical methods in our framework. Winters et al (1990) and Maini et al (1991) perform branch tracing numerically for their twocomponent RD system with cross diffusion acting on simple 2D rectangular domains. They simulate a diverse range of complex patterns with significant amplitudes for studying snakeskin pigmentations. Yochelis et al (2008) perform numerical branch tracing for a twocomponent GiererMeinhardt RD system on a simplified periodic domain to study cardiovascular calcification patterns. They use the insights gained from branch tracing to characterise the parameter space for further experiments with domains. Chien and Liao (2001) investigate multiple modes bifurcating at a given bifurcation point for a twocomponent Brusselator system subject to Robin boundary conditions. They demonstrate numerical continuation of multiple branches due to mode interactions for a square domain. Paulau (2014) perform numerical branch tracing for a twocomponent FitzHughNagumo (FHN) RD system on a planar domain to study the properties of its localised solutions, i.e. solitons. With this, he establishes the existence and stability of certain first higher order radially symmetric solitons with nonzero azimuthal quantum number^{*}^{*}*i.e, characterising the zero crossings of a circumferential profile of a radially symmetric soliton. which require the third and fifth order nonlinear reaction terms to produce them. Our framework generalises such studies with twocomponent RD systems to arbitrary 2D surfaces.
Other Cases.
In general, studies with emergent patterns, marginal stability or bifurcation analysis may deal with more than twocomponents (Qian and Murray, 2001), coupled layers (Yang et al, 2002; Vasquez, 2013), quasi equilibrium (Rozada et al, 2014), advection (Vasquez, 2013; Satnoianu et al, 2001; Madzvamuse and Zenas George, 2013), a very large number (ZamoraSillero et al, 2011) or range (Lo et al, 2012) of control parameters, shearinduced instability (Vasquez, 2013), growth induced instability (Madzvamuse, 2008), nonlinear diffusion (Gambino et al, 2013b), fractional RD (Gafiychuk et al, 2009) or even nonsteady state (oscillatory) or travelling wave solutions (Draelants et al, 2013; Banerjee and Banerjee, 2012; Wyller et al, 2007; Qiao et al, 2006; Gambino et al, 2012). While our framework may be used directly or with simple modifications for only a few of these general cases, they serve as directions for future work on our framework.
System Discretisation, Sparsity, largescale and Numerical Methods.
To numerically solve RD systems for emergent patterns and bifurcation branches, we must first discretise the PDEs. Several existing techniques allow discretisation based on either discrete geometric operators (Botsch et al, 2010), finite differences for closestpoint methods (Ruuth and Merriman, 2008; Macdonald et al, 2011) or Finite Element Methods (FEM). We employ FEMs owing to their generality and robustness and use Deal.II (Bangerth et al, 2007), a general purpose software library that approximates arbitrary surface geometries with quadrilateral finite elements. With discretisation, a generic RD system is expressed as a linear system of equations where individual scalar elements may still be computed using nonlinear operations on discrete surface variables.
With discretisation and linear stability analysis near homogeneity, the problem of computing emergent patterns translates into a problem of detecting bifurcations and solving a generalised eigenvalue problem at detected bifurcation points (Winters et al, 1990). Seydel (2010) discusses standard test function based techniques for detecting bifurcation points along any branch. These detection methods may fail, as strategies for evaluation (and construction) of test functions are empirical in nature. In contrast to these existing detection based approaches for bifurcations on a general branch, we exploit the spectral properties of a potential pattern emergent from a trivial branch to directly compute its bifurcation point. Seydel also describes methods such as the method of parallel computation to solve for emergent patterns along any branch for switching (see Seydel, 2010, Chapter ). With our simplifications near homogeneity, we directly compose an emergent pattern and exploit the principles of the method of parallel computation to switch over to a new branch. We then use a pseudo arclength continuation approach for branch tracing as presented by Salinger et al (2002) which is also easily scalable to large systems (Salinger et al, 2001).
Several natural steady state patterns are fairly complex and we need a sufficiently large system of discrete variables to study them. These largescale systems result in sparse system matrices and we need iterative algorithms such as Krylov space methods (Gutknecht, 2007) to solve the aforementioned eigenvalue problems (Arbenz et al, 2012). In our framework, we use the Trilinos Project (Heroux et al, 2005) software libraries that implement various numerical algorithms for largescale systems with sparse matrices on parallel processing architectures. When dealing with largescale systems it is important to use preconditioners to avoid numerical instabilities due to poor conditioning of matrices. In our framework, we use algorithms based on incomplete factorization schemes, as implemented in packages like IFPACK and AztecOO for the Trilinos libraries.
For largescale systems under consideration, we need multigrid methods to make iterative algorithms practicable that otherwise generate smooth but stubborn residuals with poor convergence rates (Strang, 2007, Chapter ). Analysing computational costs for iterative Krylov subspace methods is dependent on data, convergence criteria and numerical stability, and not generalisable (Liesen and Strakos, 2012, page 248). Furthermore, for nonsymmetric matrices that arise with twocomponent RD systems, the convergence of popular methods like biconjugate gradient with stabilisation (BiCGStab) and generalised minimal residual (GMRES) that we use, is poorly understood (Simoncini and Szyld, 2007). However, multigrid methods coarsen a grid to increase the spatial frequency of residual errors which helps to accelerate their removal (Strang, 2007, Chapter ). In practice, more than two grid levels and several passes (Vcycle, Wcycle or fullmultigrid (Strang, 2007, Section )) across multiple levels are required. Also, multigrids may be classified as geometric (Landsberg and Voigt, 2010) or algebraic (Falgout, 2006) (Strang, 2007, Chapter ), and used as preconditioners (Sala and Heroux, 2005) or even directly improve efficiency and robustness of branch tracing (Chien and Jeng, 2006). Baker et al (2011) demonstrate that the scalability and performance of algebraic multigrids is highly dependent on the multicore platform architecture and needs expert or empirical programming efforts for optimization. In addition, tightly coupled multilevels for existing multigrid approaches make distribution and scheduling of computational workload across parallel architectures difficult, owing to interdependencies across grid levels and decomposed domain partitions. In contrast, we propose a simple progressive geometric multigrid approach that decouples branch tracing from scale improvements to allow a highly parallelisable implementation. Our approach is based on a multigrid continuation approach by Bank and Chan (1986) where branch tracing is performed at the coarsest level and resultant solutions are refined for higher resolution meshes. We adapt this approach to perform alternative resolution and geometric improvements iteratively and demonstrate its working for a very high mesh resolution ( FEM nodes) for an arbitrary surface. The inherent benefits of this decoupled progressive approach includes selective and parallel refinements of branch solutions.
3 Direct Linear Analysis with Laplacian Eigenbasis
In this section, we formulate our direct approach to compute bifurcation patterns using eigenfunctions of the LaplacianBeltrami operator. We derive a general form for an emergent pattern near homogeneity for a generic twocomponent RD system acting on arbitrary surfaces, with or without crossdiffusion. Also, we derive the constraints to be satisfied by the continuation parameter for locating a bifurcation point. We first perform a simplification of the generic RD system equations near homogeneity into a linear form, to be satisfied by an emergent pattern (Section 3.1). Next, we present a spectral decomposition of potential patterns and the boundary conditions that allow expressing these patterns with orthogonal basis functions (Section 3.2). Then, we substitute a spectrally decomposed potential pattern into the linearised system equations to obtain an explicit general form for an emergent pattern along with expressions and conditions for its spectral coefficients (Section 3.3). We also define the bifurcation point in terms of spectral eigenvalues and system parameters. We discuss three cases of simple, multiple and mixedmode bifurcations and the constraints that they impose on our general derivations, and we discuss how to directly compose bifurcation patterns for these cases (Section 3.4). Also, we discuss common applications such as exclusive mode selection with isotropic domaingrowth .
3.1 Linearising Generic TwoComponent RD Systems
Let us consider a twocomponent general RD system with cross diffusion, defined over an arbitrary surface. It can be expressed mathematically as
(1) 
Here, and are the concentrations of two components over the surface , diffusion is represented with the LaplaceBeltrami operator , and functions and represent nonlinear reaction terms. The scalar coefficients and are positive diffusion rates, and are nonnegative selfdiffusion factors and and are nonnegative crossdiffusion factors for the system (Lou and Ni, 1996). Throughout this paper, we consider only steady state solutions of RD systems. For Equation 1, this implies that we are interested in solutions with . To linearise the RD system defined in Equation 1, we first define its homogeneous steady state as a solution to the simultaneous equations and . Note that depending on the complexity of and (say polynomial order), there may be multiple choices for the homogeneous steady state . Given a steady state , we now perform a Taylor series expansion for the nonlinear reaction terms and for infinitesimal deviations and ,
(2) 
where and are polynomial functions containing the second and higher order terms in and for their respective Taylor series expansions. In other words, and represent the nonlinear part of the reaction terms for the system defined by Equation 1 near the homogeneous steady state . Now substituting , , , and from Equation 2 into Equation 1 and ignoring the nonlinear terms yields
(3) 
with (4)
for and . Here, new diffusion coefficients and reaction coefficients are defined in terms of old coefficients in Equation 1, , and partial derivatives , , and evaluated at . See the supplemental material (SM01.D1) for the definition of these new coefficients along with necessary derivations. We emphasize that for the above linearisation, all nonlinear terms in , , , , and are considered to be negligible and ignored since and near homogeneity.
3.2 Spectral Decomposition and Boundary Conditions
Our framework performs bifurcation analysis near homogeneity using spectral analysis of the LaplaceBeltrami operator . Any surface function such as and in Equation 3 can be expressed in terms of eigenmodes of . Also, if a secondorder linear operator such as is Hermitian, then its eigenmodes form a set of orthonormal basis functions. Using an orthonormal spectral decomposition, we derive the form and conditions for an emergent pattern directly in terms of the eigenmodes and eigenvalues for the LaplacianBeltrami operator.
Now, in order to ensure the orthonormality of the basis functions, we need to consider the boundary conditions that and are subject to. Most biological problems expressed as RD systems are subject to either periodic boundary conditions or zero Dirichlet, Neumann, or Robin boundary conditions near homogeneity, or they deal with closed surfaces without boundaries. We show in the supplemental material (SM01.D2) that all these cases satisfy the Hermitian property of the LaplaceBeltrami operator.
Spectral decomposition.
To generate an orthonormal basis set using the LaplaceBeltrami operator , we must solve the corresponding eigenvalue problem
(5) 
Here, all the eigenvalues are real and nonnegative since is Hermitian. We thus assume that the basis functions form an ordered set, where the ordering index satisfies . Using the basis functions we can express a smooth surface function as
(6) 
and are called spectral coefficients. Next, we leverage the spectral decomposition in Equation 6 to express a pattern emergent near homogeneity for a generic RD system defined in Equation 1.
3.3 Bifurcation Patterns near Homogeneity
While the spectral decomposition suggests that potential patterns could in general contain any superposition of eigenmodes, the conditions near homogeneity impose additional constraints on actual emergent patterns. We now derive an as–general–as–possible form for emergent steadystate patterns of our generic RD system formulation that respects these constraints.
Consider a steadystate bifurcation pattern , where the superscript denotes that this emergent pattern is a bifurcation from the trivial homogeneous solution. Substituting its spectral decomposition from Equation 6 into Equation 3, and setting the temporal derivatives to zero to obtain a steadystate solution, gives us
(7) 
For simplicity, we dropped the superscript from the spectral coefficients and . Simplifying Equation 7 using substitutions , , and imposing linear independence of orthonormal basis functions yields the following relations between the spectral coefficients , , and eigenvalues , (see supplemental material (SM01.D3) for a detailed derivation),
(8) 
For each nonzero pair of and , we derive the following constraint on the corresponding eigenvalue from the above equation,
(9) 
Equation 9 is quadratic in and it admits at most two real valued roots for , say and . This implies that at most two set of and may be combined linearly to form an emergent pattern. Thus the most general form for a bifurcation pattern emergent near homogeneity for an RD system as in Equation 3 is
(10) 
Now, an important requirement for patterns given by Equation 3.3 to emerge due to diffusion driven instabilities is that the RD systems must be linearly stable in absence of diffusion. Murray expresses this requirement in terms of the differentials of the reaction terms in an RD system equation Murray (see 2003, Equation ). For RD systems given by Equation 1, near homogeneity, the stability constraints are expressed in terms of system parameters as
(11) 
To facilitate direct composition of emergent patterns we classify them based on two criteria. First, we classify patterns as (i) exclusive mode selections or (ii) nonexclusive mode selections based on certain conditions for the diffusion induced instability. Diffusion driven instabilities may activate multiple eigenmodes to grow or emerge simultaneously under a fixed set of system parameters. This is captured by the dispersion relation (Murray, 2003, page ), which indicates the range of eigenvalues that are unstable for the fixed system parameters. The dispersion relation provides a growth rate for each eigenvalue , and only eigenmodes with nonnegative growth rates are unstable and may participate in pattern formation.
Figure 1 plots dispersion relations as curves of growth rates over potential eigenvalues . The shapes of the curves are determined by the given system parameters, and varying a free continuation parameter leads to a family of curves as shown here. Circles along the –axis indicate actual eigenvalues from the spectral analysis. In exclusive mode selection, eigenmodes corresponding to exactly one eigenvalue, say , become unstable (they have a nonnegative growth rate) at the bifurcation point. In nonexclusive mode selection, the system becomes unstable to a new set of eigenmodes with eigenvalue , while it remains unstable to other modes with positive growth rates.
Further, similar to Chien and Liao (2001)^{†}^{†}†While Chien and Liao (2001) label all bifurcations made of two or more eigenmodes as mixedmode, we reserve this term only for the bifurcations that involve eigenmodes with different eigenvalues., we classify bifurcations as: (a) simple if only a single wavemode constitutes the emergent pattern, (b) multiple if more than one wavemode constitutes the emergent pattern but all such wavemodes have the same eigenvalue, and (c) mixedmode for the bifurcations with constituent wavemodes of more than one eigenfrequency, say and . Figure 1 on the right illustrates the dispersion relation for a mixedmode pattern at its bifurcation point. Next, we derive different formulae to compose bifurcation patterns directly for these cases.
3.4 Composing Bifurcation Patterns for Continuation
Building on the derivations from Sections 3.1, 3.2 and 3.3, we now derive equations to specifically compose simple, multiple, and mixedmode bifurcations for the RD system in Equation 1 (near homogeneity), both under exclusive and nonexclusive mode selection.
3.4.1 Simple Bifurcations
Simple bifurcation patterns are composed of a single wave such that the algebraic multiplicity of its corresponding eigenvalue is . Let for some , be the steady state bifurcation pattern of interest. We now investigate the conditions that the system parameters must satisfy to branch out a steady state pattern .
Exclusive Mode Selection.
For simple bifurcations in exclusive mode selection problems, the quadratic Equation 9 must have only one real root for satisfying
(12) 
at the bifurcation point. For RD systems without crossdiffusion, and the above relation becomes . Thus for studying exclusive mode selection for RD systems without crossdiffusion, the continuation parameter must be associated with either or .
An interesting class of problems deals with isotropic growth of the domain as the mode selection criterion (Murray, 2003, page ). For such cases, let us introduce a common scale factor for all the parameters to represent isotropic growth ^{‡}^{‡}‡See Murray’s book (Murray, 2003, page ) for an interpretation of the scale factor ., i.e. , , and . Substituting these terms in Equation 12 gives us at the bifurcation point directly in terms of the system parameters as
(13) 
NonExclusive Mode Selection.
For simple bifurcations that emerge nonexclusively, the for a given wavemode may not be a unique root for Equation 9. Nevertheless, for detecting bifurcations near homogeneity with all but one unknown system parameter, we can use Equation 9 with to directly compute the free (continuation) parameter.
Thus, for all simple bifurcations along the trivial branch , we can directly compute a continuation parameter to locate the bifurcation point for a given wavemode by substituting in either Equation 12 (exclusive mode selection), 13 (exclusive mode selection under domain growth) or 9 (nonexclusive mode selection). Finally, with all system parameters known, we can obtain the spectral coefficients and by solving the simultaneous Equations 8. This yields the desired bifurcation pattern , and we can switch to the new branch for tracing (see Section 4.5)^{§}^{§}§Note that Equation 8 is homogeneous and the terms and can only be solved upto a common scale factor, say ..
3.4.2 Multiple Bifurcations
For an eigenvalue with algebraic multiplicity , we have multiple candidate wavemodes that satisfy the conditions and equations that we derived for simple bifurcations. Thus, for such cases, every linear combination of these wavemodes is a bifurcation pattern emergent at the same bifurcation point located using the results from Section 3.4.1. We discuss branch switching and continuation for such multiple bifurcations in more detail in Section 4.
3.4.3 Mixedmode Bifurcations
Let us now consider an as–general–as–possible emergent steady state bifurcation pattern as given in Equation 3.3. This means there are eigenmodes of two different eigenvalues and , and we assume that . Let us define and . Now, since , , Equation 8 implies that all are equal, say , . Similarly, (say), . Substituting these scale factors with their respective eigenvalues in Equation 8 yields^{¶}^{¶}¶See our supplemental material (SM01.D4) for detailed derivation for these coefficients and Equation 15.
(14) 
Now, one of the linear stability requirements in absence of diffusion is that (Equation 11). Expanding this inequality with substitutions from Equation 14 gives us a constraint on the diffusion parameters,
(15) 
which needs to be satisfied to obtain a mixed mode bifurcation. Next, using the fact that and are solutions for in Equation 9, we can solve for the continuation parameter which is the only unknown in the following equation,
(16) 
Finally, with all system parameters determined, we can compute and as
(17) 
In summary, to study a mixedmode bifurcation we start composing a desired emergent pattern by selecting two eigenvalues, and , and the corresponding sets of eigenmodes and . In addition, we freely choose desired spectral coefficients and . We also set values for all but one system parameter such that they satisfy Equation 15. Then, we compute the unknown continuation parameter by solving the (quadratic) Equation 16. In order to continue with branch tracing, the solved bifurcation parameter must be real valued and satisfy preconditions in Equation 11. Else, our framework reports an error. Finally, we compute and using Equation 17, and the spectral coefficients and using , , and . Now, all terms are determined to compose the bifurcation pattern using Equation 3.3.
4 Numerical Method
In this section we describe the numerical implementation of our framework in more detail. We build on the Trilinos and Deal.II libraries to provide an implementation of our proposed method for bifurcation analysis near homogeneity. Our framework uses FEMbased surface discretisation using quadrilateral finite elements (QFEs), and the order of the FEs can be configured up to three. We support the following operations for bifurcation analysis and branch tracing: approximating Laplacian eigenfunctions (Section 4.1), resolving bifurcations (Section 4.2), and branch tracing (Section 4.5). In addition, we will describe a strategy for resolving patterns at higher resolutions in Section 5. Our framework also includes a generic, indirect reference method for branch detection using a test function (see Seydel, 2010, Section ). We use it for various comparisons during the evaluation of our framework.
4.1 Approximating Laplacian Eigenfunctions
We saw in Section 3.3 that the eigenfunctions of the Laplacian of a surface domain are the building blocks for composing bifurcation patterns, given an RD system satisfying one of the boundary conditions discussed in Section 3.2. Hence computing these eigenfunctions is a core functionality of our approach. Interestingly, the eigenfunctions depend only on the shape of the surface domain and the boundary conditions. They are independent, however, of the RD system formulation and its parameters. We use an FEM discretisation of the domain and apply a Galerkin method to discretise the eigenvalue problem in Equation 5 as
(18) 
Here, is a discrete vector representation of the eigenfunction . We obtain the stiffness matrix and the mass matrix by applying a weak formulation integration to the Laplacian operator and the eigenfunction , respectively. As before, is the eigenvalue corresponding to the eigenfunction . Equation 18 is a generalised eigenvalue problem, and in our case we obtain a largescale system of sparse matrices. We use the Anasazi eigensolver package from the Trilinos library for finding the eigenvectors . For largescale general eigenvalue problems a shiftinvert approach is commonly used to solve for a band of eigenvectors as recommended by Lévy and Zhang (2010). However, for RD systems with zero Neumann boundary conditions, is singular and the shiftinvert method cannot be used to compute the lowest frequency band of eigenvectors. To keep things simple, we apply as a preconditioner to both sides and solve the resulting standard eigenvalue problem. We avoid explicit computation of the possibly nonsparse, large matrix with the use of the AztecOO package from the Trilinos library. AztecOO provides an inner loop implementation for each application of matrix to a vector (say) for computing by solving the linear system instead of matrix inversion.
4.2 Resolving Bifurcations
Given the eigenvectors and their respective eigenvalues , we can now compose bifurcation patterns and locate their point of emergence on the trivial branch. Section 3.4 explains how a bifurcation pattern may be categorised as a simple, multiple or mixedmode bifurcation. For simple bifurcations, each basis vector defines a bifurcation pattern, which we denote as . We first locate the bifurcation point corresponding to as follows. Let be the unknown continuation parameter in . Depending on the problem, we use one of the Equations 9, 12 or 13 to compute by plugging in the known parameters from , and (or ). Then, without loss of generality, we set and solve for using Equation 8 to compute . For multiple bifurcations, to compose a bifurcation pattern , we first select a set of eigenvectors . We then use as in the case of simple bifurcations to compute for the corresponding bifurcation point. Next we select an arbitrary set of spectral coefficients to define a desired linear combination of as an emergent pattern and compute the , using Equation 8. Thus we compute the multiple bifurcation pattern . For a mixed mode bifurcation, we pick two sets and and compose an arbitrary linear combination of eigenvectors to define a desired emergent pattern. We then solve for the bifurcation point by substituting all parameters from , and and in Equation 16. With the known bifurcation point and and , we compute and using Equation 17. Next, with the previously defined arbitrary values of and for the desired emergent pattern, we compute and . Finally, we compose the mixed mode bifurcation pattern .
4.3 Reference Method
For evaluation purposes, we also implemented a standard approach for detecting bifurcation points (Seydel, 2010, Chapter ). We call this the reference method since it uses common techniques for different tasks in detecting bifurcation points and computing bifurcation patterns. To detect a bifurcation point, we use a test function which is evaluated at (,) as, , where vector satisfies equation , with . Here, is a unit vector with all but the element set to zero and is the Jacobian matrix for function as discussed for branch tracing in Section 4.5. A bifurcation is detected each time changes its sign from to in an interval, say . Next we perform a midpoint search to locate the crossing for by iteratively evaluating it at and then setting if at or otherwise setting . Upon convergence . This implies that vector lies in the nullspace of the Jacobian matrix at and it is a good approximation to the bifurcation pattern . We thus output and as the next detected bifurcation pattern which may be used for branch tracing in a manner similar to the previous methods. Note that this method is not capable of resolving multiple bifurcations.
4.4 Advantages and Limitations of Proposed Approach
The main advantage of our proposed approach is that finding bifurcation points and patterns is not dependent on the goodness of the tracing test function. Using a test function poses two issues: first, the potential presence of more than one bifurcation in an interval, and second, the lack of a guarantee for any test function to detect the presence of each bifurcation. While there are several strategies to address these two issues, often incompletely, our proposed direct approach completely avoids them. Furthermore, we do not need to perform repeated evaluations of the test function along the trivial branch. As a further advantage, depending on the complexity of the test function and the interval size, our approach reduces computational overheads. We tabulate performance gains due to our approach in Section 6. A third, important advantage of our approach is that it allows the direct composition of infinite multiple and mixedmode bifurcation patterns. This adds considerable possibilities to bifurcation analysis by supporting the exploration of multiple colocated branches.
At present, the main limitation of our proposed approach is that it can only be used to analyse primary branches emerging from the trivial branch.
4.5 Branch Tracing for Nonlinear Analysis
Our framework also supports nonlinear analysis with branch tracing in the faroff nonlinear region for the PDEs. The first task in branch tracing is to switch over to the new branch of patterns characterised by a given bifurcation pattern . While the linear stability analysis gives us , we are interested in a nontrivial solution that fully satisfies the actual nonlinear PDE given in Equation 1, near the bifurcation point (with ) on the trivial branch. With as the homogeneous solution at the bifurcation point , estimating a good starting point with on the new branch is a nontrivial task. The nonlinear solution must be qualitatively similar to the bifurcation pattern ; yet farenough away from the trivial branch to allow continuation without falling back. To achieve this, we propose two improvements over a standard method of parallel computation approach for branch switching (see Seydel, 2010, Section ).
For the first improvement, we suggest using a bordering algorithm (Salinger et al, 2002) to iteratively compute the jump until a successful switch is made. Our key idea is to select a pattern dependent pivot (a discrete FEM node) for fixing the jump size and the direction of parallel computation. As a second improvement which helps multiple and mixed mode bifurcations we propose to apply a strong guidance to the jump at each intermediate step of our iterative switching algorithm. The bifurcation pattern serves as a good guidance for the jump . The details for our proposed improvements are presented in the supplemental material (SM02.A1).
Once we switch over to a new branch by jumping from the trivial solution to the nonlinear solution , we follow it by means of continuation. Our framework uses the LOCA and NOX packages from the Trilinos library to perform pseudo arclength continuation. We use an adaptive approach that updates the step size after each continuation step and impose a tangent scale factor to manoeuvre the direction of the continuation curve as supported by the Trilinos library. In particular, we propose to establish a initial tangent direction which is strongly orthogonal to the trivial branch by attempting a jump from the nonlinear solution to its antithetic solution , instead of jumping from the trivial solution to the nonlinear solution . Again, we present the details of our proposed improvement for branch tracing with other implementation details in the supplemental material (SM02.A1). Also, we provide further details about the configurability of our framework in the supplemental material (SM02.A2).
5 Multiresolution Adaptation
Unlike simple spot and stripe patterns, most interesting biological surface patterns exhibit high shape contour irregularities. The complexity of these patterns is attributed to the surface geometry and the (mid or high) frequency of the constituent eigenmodes (of the LB operator ). Thus, for accurate computation of such complex patterns, the surface geometry must be represented with a highresolution FEM discretisation. As discussed in Section 2, many existing geometric and algebraic multigrid approaches support highresolution meshes but face challenges in either performance, scalability or parallelisability for branch tracing owing to tightly coupled multiple levels. We thus propose a simple multilevel approach in which branch tracing is performed at the baselevel of discretisation and the resultant patterns are upsampled and resolved at higher levels progressively.
Our framework uses a simplified geometric multigrid approach where the surface domain is organized in multiple levels as , with , where is the lowest resolution representation and . We begin with the highest level mesh and generate each lower level mesh from mesh by applying a quadricbased edge collapse decimation algorithm (Garland and Heckbert, 1997; Cignoni et al, 2008). We perform branch tracing at the lowest resolution and progressively upsample resulting patterns up to the highest resolution . Unlike most of the multigrid approaches where all the mesh levels are used simultaneously in a Vcycle or a Wcycle (for the full multigrid approach) we perform complete upsampling of a solution using only two levels at a time. We thus call our approach a progressive geometric multigrid approach.
We propose a twostep approach for upsampling the results between two levels and . For the first step, we increase the mesh resolution for without changing its geometry. We perform an inplane subdivision of each existing (quad) finite element into four to give a new mesh, say . In the second step, we map the geometry of the mesh to . For each solution vector defined over with continuation parameter , we first interpolate it linearly to the higher resolution mesh and then solve the nonlinear system over . In general, we use a Newton method with backtracking to solve for . We compute the direction vector for the Newton method using a biconjugate gradient method with stabilization as implemented in the AztecOO package. For difficult cases, we use a trust region method for solving the above nonlinear system with a GMRES approach for establishing the search direction. Next we perform a similar interpolation from to and then solve for with .
For linear interpolation of a solution across two meshes (say, from to or from to ) we use a projection based mapping scheme between meshes. We project each node from a target mesh (say ) onto the nearest face of the source mesh () and use the barycentric coordinates of the projected node to compute linear interpolation weights. A naive implementation of this mapping scheme has a computational complexity with and as the number of nodes for the source and target meshes respectively. Computational complexity can be reduced to with the use of a kdtree for the nearest face search. Note that we compute the multilevel meshes and the mappings for linear interpolation of solutions only once as a preprocessing step for a given surface domain .
Our twostep approach separates the complexity of resolution improvements from the complexity of geometry improvements for an upsampling task. The results for our progressive geometric multigrid approach are presented in Section 6.
Parallelisability.
Branch tracing is inherently a sequential task, because the next solution on a branch is computed by numerically analysing the previous one. For largescale systems, Salinger et al (2005); Lehoucq and Salinger (2001) suggest using a multiprocessor environment where a large domain () is partitioned into subdomains and each domain is solved on a separate processor. Their domainpartitioning approach is generic but the parallelism is limited to solving for one solution at a time while tracing the branch sequentially. Similarly, Continillo (Continillo et al, 2012) present and discuss a method for parallelisation of the most repeated operation for their reactive dynamical system, i.e. the computation of the Jacobian Matrix . They parallelise the computation of using a cluster with a message passing interface (MPI). Again their approach limits parallelisation to the computation of the next solution.
Our framework has considerable potential for parallelism. The most important of it is the scope for upsampling several patterns on a branch simultaneously. This is a direct outcome of our progressive geometric multigrid approach which completely delineates the branch tracing problem from the convergence of a solution at a higher resolution. Once a branch is traced at the lowest resolution with domain , our framework can be invoked to resolve each of the patterns on the branch independently. This allows several instances of our framework to be launched on a computer cluster to simultaneously solve for all the solutions on a branch at higher resolutions. We evaluate the overhead of such parallel launches due to file reading and mesh data preparation operations in Section 6 and provide theoretical estimates for potential gains from parallel computation of a given branch.
6 Experiments and Results
In this section, we describe our experimental setup, discuss case studies and present results. We use a bit Ubuntu LTS platform running on an Intel Xeon E52630 CPU with Cores @Ghz with GB RAM for all performance evaluations and most other experiments. We demonstrate our framework with two RD system case studies, a Brusselator system for the cotyledon patterning of conifer embryos (Section 6.1) and Murray’s chemotactic model for snakeskin pattern formation (Section 6.2).
6.1 Case Study I: A Brusselator model
Nagata et al (2013) present a study of emergent cotyledon patterns on a plant tip. They model observed patterns as bifurcations for a twocomponent Brusselator RD system. They represent the plant tip with a simple geometric shape, a spherical cap. This spherical cap domain is parameterised by its size factor and a curvature factor as shown in Figure 2. Nagata et al (2013) perform marginal stability analysis to study the relation between the number of emergent cotyledons and the size factor or the curvature factor for a cap. Mathematically, their model is defined as,
(19) 
Here, and are concentrations of two Turing morphogens, , , and are positive rate constants, and are positive diffusion rates and is an infinitesimal element on the boundary for the domain . Nagata et al. linearise Equation 19 near homogeneity with to derive an analytical expression for their continuation parameter (represented by here) at a simple bifurcation point. The continuation parameter is determined by the other system parameters and the eigenvalue of a spherical cap harmonic , which constitutes the emergent pattern. Note that , in turn, depends on the shape and size of the cap, and it can be computed by evaluating an associated Legendre function. This makes it possible to study the marginal stability of emergent patterns with respect to , and .
In our framework, the spherical cap harmonics are just a special case of eigenfunctions for a specific, simple domain. The benefit of our approach, however, is that we can easily generalise the analysis to arbitrarily shaped domains. This could facilitate the discovery of shapeinduced anomalies in emergent patterns, as we discuss in an example later. Further, we simplify the analysis by incorporating the size factor directly into the RD model. We replace the size of the spherical cap with a factor that represents the relative scale of an arbitrary domain with respect to its canonic unit size. We then numerically compute the emergent patterns and corresponding values.
We first modify Equation 19 to include the scale factor in all the reaction terms, that is, in all coefficients except the diffusion rates and . Scaling the reaction relative to the diffusion coefficients has the same effect as scaling the domain (Murray, 2003, page ). We use the notation , where indicates corresponding rates at unit scale, and similarly for the other coefficients. This yields new linearisation parameters
(20) 
Substituting these parameters along with , , , and in Equation 9 gives us the new continuation parameter as
(21) 
Here denotes the eigenvalue for a pattern at unit scale of the domain . It can be shown that under uniform scaling of a surface domain by a factor , the eigenvalue for a given eigenfunction is inversely proportional to . For spherical cap harmonics, for example, eigenvalues are analytically defined as , where is defined in terms of an associated Legendre function dependent on the wavemode and the curvature factor (see Nagata et al, 2013, Equation ). Thus substituting in our modified Brusselator model allows us to incorporate the size factor for studying marginal stability with arbitrary domains.
We compute each eigenmode and its respective eigenvalue numerically only for a representative arbitrary shape at a unit scale. Our modified Brusselator model then enables us to plot marginal stability curves for different bifurcation patterns on this arbitrarily shaped plant tip against the size factor , similar to the plots for spherical caps as presented by Nagata et al (2013) in Figure of their paper.
Experiments.
First we validate our framework against analytically derived results from Nagata et al (2013). We use the same parameter values as given by Nagata et al. for their Figure (i.e. , , , , ) and a spherical cap with radius and curvature factor . For all quantitative evaluations in this section, we compute seven different emergent patterns as shown in Figure 3 and their corresponding bifurcation points using both our proposed method, and the reference method, as explained in Section 4. We begin by examining errors in computing with an FEM mesh of order one with about nodes.
Table 1 provides a numerical comparison of relative errors for both methods against the analytically derived values. For all emergent patterns, our proposed method has a low relative error on the order of when compared with the expected analytical results, and the errors are quantitatively similar to those for the reference method. This implies that our direct approach works well and most of the errors can be explained as approximation errors due to FEM discretisation.
Figure 4 shows relative errors for locating all seven bifurcation points at five mesh resolutions ranging from  vertices. The bifurcation points are ordered according to increasing eigenvalues along the horizontal axis. We found the relative error to fall quickly as a powerlaw function of the mesh resolution for each pattern. Mode is closest to the exclusive mode selection criterion in Equation 12 and it has the lowest relative error. We observe that the errors increase with the distance from this reference mode in the eigenspectrum and expect larger relative errors in locating bifurcation points for higher frequency patterns. Thus for applications requiring high accuracy for studying emergent patterns with arbitrary domains, our multiresolution approach presented in Section 5 is beneficial.
Next, we examine numerical errors in computing bifurcation patterns for our proposed method in comparison with analytically defined spherical cap harmonics (i.e. the ground truths). Figure 4(a) shows relative root mean square (RMS) errors in emergent patterns with different mesh resolutions. For comparison we normalize all spherical cap harmonics to have unit amplitude. Again we see that the accuracy improves with the discretisation resolution with some powerlaw function. We also include the errors for the reference method for a mesh with nodes in Figure 4(a), which indicates that our proposed method is quantitatively consistent with the reference method.
We also studied the impact of the FEM order on RMS errors with about nodes in each case, see Figure 4(b). We denote surface geometries that are approximated with finite elements of order , and , each with about nodes, by , , and . Since FEM of increasing order has an increasing number of nodes per planar finite element, the spherical surface geometry is approximated with a decreasing number of planar elements for higher orders. Comparing the solutions of order FEM on , order on , and order on , we see that FEM order outperforms higher orders for all emergent patterns. At the given resolution of FEM nodes, the increase in error due to the geometric approximation in and apparently outweighs the error reduction due to higher order polynomials. In addition, we also report the error of order FEM on the lower resolution geometries and . As expected this increases the error due to the use of lower order polynomials, but only marginally.
Next, we illustrate the key advantage of our approach, that is, the ability to study the effects of arbitrary shape distortions on emergent patterns. Figure 6 shows one such distorted shape in its second row, which we generate by deforming the circular boundary of the cap and propagating the distortions smoothly over the entire surface. We visualize the effects of these shape distortions on three emergent patterns using a nonlinear color mapping, as illustrated in the figure. Clearly, new patterns show marked deviations from the respective spherical cap harmonics shown at the top row. Deviations from normal emergent patterns often lead to developmental anomalies, which are studied extensively in mathematical biology, refer to Harrison and Von Aderkas (2004) for an example. While our illustration here does not explain any specific anomaly, our framework can certainly be used to study the role of domain shape deviations on actual observed cases.
Finally, we present marginal stability curves for mode in Figure 7. Here we consider three cases: (a) a spherical cap growing isotropically (greenish colour), (b) a distorted cap growing isotropically (brownish colour), and (c) a spherical cap that first grows isotropically until and then progressively distorts with further growth (bluegray colour). We first computed eigenvalues for the spherical cap with , corresponding to a boundary perimeter , and , and for the distorted cap with boundary perimeter . Using these eigenvalues we plot solid curves in Figure 7 for