Bifurcation Analysis of Reaction Diffusion Systems on Arbitrary Surfaces

Bifurcation Analysis of Reaction Diffusion Systems on Arbitrary Surfaces

Daljit Singh J. Dhillon Daljit Singh J. Dhillon Department of Computing,
Imperial College London, London, United Kingdom.
Note: While conducting this work, the first author was a Ph.D. student at the Institute of Computer Science, University of Bern, Bern, Switzerland. Michel C. Milinkovitch Department of Genetics and Evolution, Laboratory of Artificial and Natural Evolution (LANE),
University of Geneva, Geneva, Switzerland
44email: michel.milinkovitch@unige.chMatthias Zwicker Institute of Computer Science,
University of Bern, Bern, Switzerland.
   Michel C. Milinkovitch Daljit Singh J. Dhillon Department of Computing,
Imperial College London, London, United Kingdom.
Note: While conducting this work, the first author was a Ph.D. student at the Institute of Computer Science, University of Bern, Bern, Switzerland. Michel C. Milinkovitch Department of Genetics and Evolution, Laboratory of Artificial and Natural Evolution (LANE),
University of Geneva, Geneva, Switzerland
44email: michel.milinkovitch@unige.chMatthias Zwicker Institute of Computer Science,
University of Bern, Bern, Switzerland.
   Matthias Zwicker Daljit Singh J. Dhillon Department of Computing,
Imperial College London, London, United Kingdom.
Note: While conducting this work, the first author was a Ph.D. student at the Institute of Computer Science, University of Bern, Bern, Switzerland. Michel C. Milinkovitch Department of Genetics and Evolution, Laboratory of Artificial and Natural Evolution (LANE),
University of Geneva, Geneva, Switzerland
44email: michel.milinkovitch@unige.chMatthias Zwicker Institute of Computer Science,
University of Bern, Bern, Switzerland.
Received: date / Accepted: date

In this paper we present computational techniques to investigate the solutions of two-component, nonlinear reaction-diffusion (RD) systems on arbitrary surfaces. We build on standard techniques for linear and nonlinear analysis of RD systems, and extend them to operate on large-scale 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 Laplace-Beltrami 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 large-scale 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

Reaction-diffusion (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 steady-states. 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 low-resolution 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 step-by-step, our proposed framework uses an analysis-synthesis approach to directly determine bifurcation points and construct emerging patterns along the trivial branch. We exploit the Hermitian nature of the Laplace-Beltrami (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 detection-based approaches, our analysis-synthesis approach avoids missing out on bifurcations due to potential failures of a test function. In addition, our approach allows for tracing multiple and mixed-mode bifurcations apart from simple bifurcations.

Our framework also supports high-resolution 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 two-step 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 two-step 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 large-scale system. Unlike traditional multigrid approaches, our two-step 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 two-component RD systems with or without cross diffusion that supports arbitrary triangulated surface domains.

  • A direct analysis-synthesis approach for computing emergent patterns and locating bifurcation points along the trivial branch based on spectral analysis of the Laplace-Beltrami operator on arbitrary triangulated surface domains.

  • A progressive geometric multigrid approach supporting high-resolution 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 two-component 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 cross-diffusion for Lotka-Volterra kinetics between two components in a 2D rectangular domain, commonly used to model predator-prey 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 two-component RD system. They also apply a weak nonlinear stability analysis to predict the long-term behaviour of these emerging vegetation patterns.  Murray (2003, Chapter ) demonstrates the effect of both geometry and scale on the emergence of patterns under a two-component 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, well-defined 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 cross-diffusion 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 mutual-exclusivity 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 two-component 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 two-component 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 two-component Gierer-Meinhardt 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 two-component 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 two-component FitzHugh-Nagumo (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 non-zero 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 two-component RD systems to arbitrary 2D surfaces.

Other Cases.

In general, studies with emergent patterns, marginal stability or bifurcation analysis may deal with more than two-components (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 (Zamora-Sillero et al, 2011) or range (Lo et al, 2012) of control parameters, shear-induced instability (Vasquez, 2013), growth induced instability (Madzvamuse, 2008), nonlinear diffusion (Gambino et al, 2013b), fractional RD (Gafiychuk et al, 2009) or even non-steady 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, large-scale 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 closest-point 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 large-scale 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 large-scale systems with sparse matrices on parallel processing architectures. When dealing with large-scale 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 large-scale 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 non-symmetric matrices that arise with two-component RD systems, the convergence of popular methods like biconjugate gradient with stabilisation (Bi-CGStab) 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 (V-cycle, W-cycle or full-multigrid (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 Laplacian-Beltrami operator. We derive a general form for an emergent pattern near homogeneity for a generic two-component RD system acting on arbitrary surfaces, with or without cross-diffusion. 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 mixed-mode 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 domain-growth .

3.1 Linearising Generic Two-Component RD Systems

Let us consider a two-component general RD system with cross diffusion, defined over an arbitrary surface. It can be expressed mathematically as


Here, and are the concentrations of two components over the surface , diffusion is represented with the Laplace-Beltrami operator , and functions and represent nonlinear reaction terms. The scalar coefficients and are positive diffusion rates, and are non-negative self-diffusion factors and and are non-negative cross-diffusion 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 ,


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


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 non-linear 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 Laplace-Beltrami operator . Any surface function such as and in Equation 3 can be expressed in terms of eigenmodes of . Also, if a second-order 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 Laplacian-Beltrami 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 Laplace-Beltrami operator.

Spectral decomposition.

To generate an orthonormal basis set using the Laplace-Beltrami operator , we must solve the corresponding eigenvalue problem


Here, all the eigenvalues are real and non-negative 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


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 steady-state patterns of our generic RD system formulation that respects these constraints.

Consider a steady-state 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 steady-state solution, gives us


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),


For each non-zero pair of and , we derive the following constraint on the corresponding eigenvalue from the above equation,


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


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


Figure 1: Dispersion relation curves for the first three bifurcations for Murray’s chemotactic RD system (Winters et al, 1990) with as a bifurcation parameter. Potential eigenvalues are represented along the -axis and the corresponding temporal growth rates are represented along the -axis. Actual eigenvalues of the LB operator are shown as circles. As we move along the trivial branch by increasing the value of the bifurcation parameter , different eigenmodes (represented by bigger brown circles) become unstable to branch out new bifurcations.

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) non-exclusive 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 non-negative 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 non-negative growth rate) at the bifurcation point. In non-exclusive 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 mixed-mode, 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) mixed-mode for the bifurcations with constituent wavemodes of more than one eigenfrequency, say and . Figure 1 on the right illustrates the dispersion relation for a mixed-mode 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.13.2 and 3.3, we now derive equations to specifically compose simple, multiple, and mixed-mode bifurcations for the RD system in Equation 1 (near homogeneity), both under exclusive and non-exclusive 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


at the bifurcation point. For RD systems without cross-diffusion, and the above relation becomes . Thus for studying exclusive mode selection for RD systems without cross-diffusion, 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

Non-Exclusive Mode Selection.

For simple bifurcations that emerge non-exclusively, 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 (non-exclusive 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 Mixed-mode 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 yieldsSee our supplemental material (SM01.D4) for detailed derivation for these coefficients and Equation 15.


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,


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,


Finally, with all system parameters determined, we can compute and as


In summary, to study a mixed-mode 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 FEM-based surface discretisation using quadrilateral finite elements (Q-FEs), 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


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 large-scale system of sparse matrices. We use the Anasazi eigensolver package from the Trilinos library for finding the eigenvectors . For large-scale general eigenvalue problems a shift-invert 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 shift-invert 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 non-sparse, 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 mixed-mode 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 912 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 mid-point 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 null-space 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 mixed-mode bifurcation patterns. This adds considerable possibilities to bifurcation analysis by supporting the exploration of multiple co-located 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 far-off 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 non-trivial 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 non-trivial task. The nonlinear solution must be qualitatively similar to the bifurcation pattern ; yet far-enough 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 arc-length 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 Multi-resolution 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 high-resolution FEM discretisation. As discussed in Section 2, many existing geometric and algebraic multigrid approaches support high-resolution meshes but face challenges in either performance, scalability or parallelisability for branch tracing owing to tightly coupled multiple levels. We thus propose a simple multi-level approach in which branch tracing is performed at the base-level 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 quadric-based 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 V-cycle or a W-cycle (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 two-step 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 in-plane 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 kd-tree 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 two-step 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.


Branch tracing is inherently a sequential task, because the next solution on a branch is computed by numerically analysing the previous one. For large-scale systems, Salinger et al (2005); Lehoucq and Salinger (2001) suggest using a multi-processor environment where a large domain () is partitioned into sub-domains and each domain is solved on a separate processor. Their domain-partitioning 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 E5-2630 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 snake-skin 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 two-component 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,


Figure 2: A spherical cap.

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 shape-induced 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


Substituting these parameters along with , , , and in Equation 9 gives us the new continuation parameter as


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.


Figure 3: Spherical cap harmonics in order of their emergence along the trivial branch with for a spherical cap.

Mode Bifurcation Point Location () Relative Error : (m, n) Analytical Reference Proposed RE0 RE1 A0 A A1 (A - A0) / A0 (A1 - A0) / A0 (5, 1) 0.76520 0.76528 0.76528 1.07E-04 1.07E-04 (0, 3) 0.76382 0.76403 0.76403 2.77E-04 2.76E-04 (2, 2) 0.76171 0.76200 0.76199 3.84E-04 3.74E-04 (3, 2) 0.76049 0.75997 0.75997 -6.88E-04 -6.91E-04 (6, 1) 0.75552 0.75475 0.75475 -1.02E-03 -1.02E-03 (1, 3) 0.75406 0.75343 0.75343 -8.41E-04 -8.37E-04 (4, 1) 0.74744 0.74792 0.74793 6.41E-04 6.49E-04

Table 1: Errors in locating bifurcation points for the Brusselator model acting on a spherical cap domain. We compare the proposed and reference methods vis-a-vis analytically defined results with lookup tables from Bauer (1986), for several emergent modes.

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.

Figure 4: Error in computing bifurcation points for seven emergent wavemodes using different mesh resolutions. The plot shows relative errors in comparison with analytically derived values for different wavemodes in order of their corresponding eigenvalues.

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 power-law 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 multi-resolution 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 power-law 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.

(a) Different resolutions

(b) Different FEM order with mesh
Figure 5: Relative mean error in computing the wavemode pattern .

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.

Figure 6: Effects of shape distortion on emergent patterns.

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 (blue-gray 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