Isolating Patterns in Open Reaction-Diffusion Systems
Andrew L. Krause, Václav Klika, Philip K. Maini, Denis Headon, Eamonn A. Gaffney
NNoname manuscript No. (will be inserted by the editor)
Isolating Patterns in Open Reaction-Diffusion Systems
Andrew L. Krause · V´aclav Klika · Philip K. Maini · DenisHeadon · Eamonn A. Gaffney
Received: date / Accepted: date
Abstract
Realistic examples of reaction-diffusion phenomena governing spatial and spatiotemporal pat-tern formation are rarely isolated systems, either chemically or thermodynamically. However, even formu-lations of ‘open’ reaction-diffusion systems often neglect the role of domain boundaries. Most idealizationsof closed reaction-diffusion systems employ no-flux boundary conditions, and often patterns will form upto or along these boundaries. Motivated by boundaries of patterning fields related to emergence of spa-tial form in embryonic development, we propose a set of mixed boundary conditions for a two-speciesreaction-diffusion system which forms inhomogeneous solutions away from the boundary of the domainfor a variety of different reaction kinetics, with a prescribed uniform state near the boundary. We show thatthese boundary conditions can be derived from a larger heterogeneous field, indicating that these condi-tions can arise naturally if cell signalling or other properties of the medium vary in space. We explain thebasic mechanisms behind this pattern localization, and demonstrate that it can capture a large range oflocalized patterning in one, two, and three dimensions, and that this framework can be applied to systemsinvolving more than two species. Furthermore, the boundary conditions proposed lead to more symmetri-cal patterns on the interior of the domain, and plausibly capture more realistic boundaries in developmentalsystems. Finally, we show that these isolated patterns are more robust to fluctuations in initial conditions,and that they allow intriguing possibilities of pattern selection via geometry, distinct from known selectionmechanisms.
Keywords
Pattern formation · mixed boundary conditions · open reaction-diffusion systems Reaction-diffusion systems (RDS) are employed to model an increasingly wide-range of phenomena [2, 41,50]. Following Turing’s work [71], an enormous literature has developed using these models to captureaspects of patterns seen in biology and chemistry [24, 49, 11, 10, 36, 46]. However, as noted by Turing, thesemodels are idealizations, and cannot account for many observations of real developmental systems [46]. Wefurther extend present efforts to increase the realism, and hence the explanatory power, of this theory byconsidering the role of boundary conditions in isolating patterns away from domain boundaries. Our mainaim here is to demonstrate that a simple choice of boundary conditions can guarantee interior localizationof patterns. By deriving these boundary conditions from a heterogeneous problem, we argue that such
A. L. Krause · P. K. Maini · E. A. GaffneyWolfson Centre for Mathematical Biology, Mathematical Institute, University of Oxford, Andrew Wiles Building, Radcliffe Observa-tory Quarter, Woodstock Road, Oxford, OX2 6GG, United KingdomV. KlikaDepartment of Mathematics, FNSPE, Czech Technical University in Prague, Trojanova 13, 120 00 Praha, Czech RepublicD. HeadonThe Roslin Institute and Royal (Dick) School of Veterinary Studies, University of Edinburgh, Easter Bush Campus, Midlothian EH259RG, United Kingdom a r X i v : . [ n li n . PS ] S e p A. L. Krause et al. pattern isolation can be understood as arising from different kinds of boundaries due to differential geneexpression across a field, with patterns forming only in some sub-region of the domain.Stationary Turing patterns, as well as other complex behaviours of RDS such as oscillations and chem-ical chaos, are genuinely non-equilibrium thermodynamic processes. Since the early work of Nicolis andPrigogine [57, 53], substantial further work has extended this thermodynamic perspective of Turing insta-bilities as a non-equilibrium phenomenon associated with open systems [10, 59, 20, 19]. However, very littleemphasis in these works is put on the role of boundary conditions. In the classical review of pattern selec-tion [4], for instance, the authors explicitly state that they will not be concerned with the role of boundaryconditions. While many authors describe the impact of boundary conditions on such non-equilibrium phe-nomena [10, 50], relatively few consider more exotic conditions beyond the standard Neumann, Dirichlet,or periodic settings [65, 15, 35]. Physically, periodic boundary conditions can be justified for a flat approx-imation of a curved manifold, which can be appropriate for patterns appearing across the entire surfaceof an organism. However, the typical Neumann boundary conditions used to study pattern formation arevery clearly a mathematical idealization of what are often open systems.Many developmental systems exhibit localized patterning within a structurally homogeneous field, forexample ectodermal structures such as hair or teeth [70, 30], or in the role of auxin in pinning plant rootinitiation [16, 21, 1]. Periodic patterning of mammalian whiskers and hair more generally, for instance,suggests a regionalisation of different patterning fields within which RDS may explain the emergence ofpatterns. However, as noted in [50], Neumann boundary conditions can exhibit patterning all the way up tothe boundary, often leading to partial patterns at the boundary (e.g. half-spots in the case of spotty patterns).This can be understood from the linear theory as the unstable eigenmodes will achieve maximal or minimalvalues at the boundary (by the Maximum Principle). Of course, inhomogeneous solutions of nonlinear RDSmay exhibit different behaviour, though in many cases it has been shown that spike solutions will approachand be pinned to boundaries, and particularly points of extremal curvature in multiple spatial dimensions[28, 48, 17]. Numerically, as we will show in Section 3, Neumann conditions will generically allow patternsto form up to the boundary for a variety of systems and parameters.There have been several different approaches to designing RDS which exhibit patterns isolated awayfrom the boundaries. A simple approach pursued by [74] was to modify the reaction kinetics at a bound-ary to push the system outside the Turing regime locally, so that patterning was restricted to an interiorregion. Similar ideas were explored with a variety of more complex spatial heterogeneities in reaction ki-netics [54, 55], or diffusive fluxes [3]. Recently, this kind of localized patterning has been justified in the linearregime either in the case of smoothly varying heterogeneity [40] or jumps (step functions) in the kinetics[37]. Similar work has been considered in some nonlinear regimes, where spike pinning, and hence patternlocalization, can be obtained via heterogeneity [77, 1]. Another approach, pursued in the context of the posi-tioning of bacterial protein clusters, is to understand pattern selection from a nonlinear theory, and look fornonlinearities which exhibit the desired isolated patterns [51, 68]. Finally, some authors have investigatedinhomogeneous or mixed boundary conditions. Dirichlet boundary conditions can have a variety of effectson patterns, including modifying parameters for which patterns are observed [45]. Closer to what we willpropose here, [15] investigated a variety of different boundary conditions for each species in a two-speciesRDS, demonstrating that different pattern selection mechanisms could be influenced by the choice of theseboundary conditions.While we are primarily interested in understanding isolated patterns, our work also falls into the largercontext of trying to connect Turing’s simple patterning mechanism with the complex reality of biological de-velopment. As Turing himself said, pattern formation occurs in stages through subsequent patterning [71].There are difficult questions of how to idealize such situations to elucidate the fundamental mechanismsat play in any particular stage. One such question is the robustness of Turing patterns, both in terms of thedependence of steady patterned states on initial data, and the strict requirements on parameters neededfor RDS to admit Turing-type patterns [46, 79, 63]. Domain growth has been suggested as one way to helpimprove robustness [9, 38, 73], as has stochasticity [78]. Here, we will suggest that RDS exhibiting isolatedpatterning will also be more robust in both enlarging the parameter space within which we see patterns, asdiscussed in [45], and in admitting fewer possible steady state solutions.In developmental settings, the source of boundaries between regions which exhibit periodic patterning,and those which do not, can vary between tissue type and the specific morphogen signalling dynamics. Inparticular, boundaries can arise due to either explicit heterogeneity in the tissue (where cells of different solating Patterns in Open Reaction-Diffusion Systems 3 types explicitly demarcate patterning regions, due to previous fate determination) or due to more complexand diffuse mechanisms involving cell state, which may be refined into sharp boundaries via, e.g., bistabil-ity. Such boundaries can be thought of as an example of positional information [47, 26]. We briefly discussthe biology of such boundary formation before presenting our modelling framework.1.1 Biological Interlude: Boundaries in Developmental PatterningActivator and inhibitor species are taken to represent specific molecules (or more accurately, signalling net-works) in biological systems. These undergo processes of synthesis, decay, and diffusion, depending onthe type of molecule employed as a signal. Most intercellular communication is achieved through gene-encoded polypeptides (proteins), with some smaller molecules also playing roles as diffusible signals. Thesynthesis of these signals is largely controlled by rates of gene expression (transcription and translation)within cells, directly in the case of proteins, and for small molecules indirectly through production of en-zymes that catalyse their formation. These molecules are secreted from cells, then diffuse, either alone orinteracting with other proteins, through subcellular structures, or the meshwork of extracellular matrix. Ul-timately, these molecules are capable of attaching to receptors on the surface of, or sometimes within, cellsin the vicinity. The effect of receptor binding is to trigger a typically multi-step process of signal transduc-tion, ultimately altering rates of gene expression (molecule synthesis), or cell behaviour (such as cell shapeor movement) [6].Many extracellular proteins that act as signals are themselves bound in the extracellular space by in-hibitory or transport proteins. Thus inactive complexes can be formed, capable of diffusing but unable tobind receptors and elicit a response. An example is the Bone Morphogenetic Protein (BMP) family, severalmembers of which have been implicated in reaction-diffusion periodic patterning systems of skin [27, 25],limb skeleton [58], and gut [76], and for which a range of distinct extracellular inhibitor proteins have beencharacterised [75].Two points relevant to this work that arise from these general features of cell-cell communication are i)that the specificity of molecular interactions is very high, and ii) that the candidate activator and inhibitormolecules in many systems do not interact directly, but rather through chains of intermediates in signallingpathways. The specificity of molecular interactions and indirect nature of activator-inhibitor interactionpermit selective spatial variation in the properties of the activator and inhibitor that we will explore inthis paper. Boundaries relevant to reaction-diffusion patterning will be defined by the geometry of theembryo or organ, and by spatial variation in the expression of genes that encode or modulate activator andinhibitor synthesis, diffusion, biological activity, and decay. Geometric boundaries are unavoidable andoften prominent in experimental systems at the cut edges of the tissue, as they have a strong influence onpattern formation [25].Variation in gene expression can come about due to differences in tissue structure that cause entirelydistinct cell types to be directly apposed. This is marked in composite organs such as skin and gut, in whichdistinct and non-intermingling cell types with deeply divergent developmental origins are organised into atightly-packed epithelial sheet and a looser mesenchyme rich in extracellular matrix. In each of these tissuesthe cells express a distinct subset of the genome, have different shapes and mobility, and the extracellularspace between the cells has very different properties. These two tissue types are separated by a specialisedplanar matrix called the basement membrane. The locations of boundaries between these components aremarked by abrupt changes in tissue structure that are readily visible at the anatomical or histological levels.On a finer scale, differences in gene expression within a particular tissue (that is, composed primarilyof cells of the same type), can arise from previous patterning events, such as an external organising regionemitting a graded signal that influences part of the patterning field. Graded boundaries may persist in thatform, or can be refined into sharp boundaries by, for example, the action of mutually antagonistic factorsthat create bistable systems, notably formation of discrete segments [7]. If the output of such systems in-fluences the production of a receptor or an extracellular inhibitor for a component of the reaction-diffusionsystem, then across the boundary this factor would be predicted to diffuse, but to lack any biological activ-ity. Boundaries characterised by differences in gene expression may be present in structurally homogeneoustissues and require application of specialised molecular methods to reveal them. For example, in vertebrateembryonic gut regionalisation, initially diffuse boundaries become sharper [44], and BMP receptor expres-
A. L. Krause et al.(a) (b)
Fig. 1 (a) A graphical representation of three kinds of boundaries, leading to the concentration profile of the intracellular protein u shown at the top. The production of this protein is stimulated as a result of the activation of black receptors after binding to thedumbbell-shaped extracellular signalling molecules. The red and blue rectangular cells represent two distinct cell states, with the bluecapable of expressing the receptor. The green oblate cells instead represent a different cell type, separated by a basement membrane(which, in our depiction, is assumed to not influence the intracellular concentration of u ). Finally the blue cells also have a gradient ofreceptor expression towards the right boundary (which itself may be a hard no-flux boundary), leading to a more diffuse intracellularconcentration of the protein u . (b) In situ hybridization of a chick gut (developmental stage E4.5) with a riboprobe to BMPR1A fromFig. 1(H) of [67]. The purple/blue stain shows the ability of the tissue to receive BMP signals, whereas the white/translucent areaslack this receptor and do not have the same response to BMP (though BMP may still diffuse throughout these regions). sion becomes distinct in different regions or segments along the tract [67]. In Figure 1(a) we depict bothkinds of boundaries, and in (b) we show an example of a boundary observed in the ability of a tissue torespond to BMP within the otherwise-homogeneous gut tract described above.Thus, the spatially varying presence or absence of a receptor, a component of its signal transductionpathway, or an extracellular inhibitor, can itself constitute a boundary relevant to the behaviour of a reaction-diffusion patterning system, even in a structurally homogeneous field of cells. Such a boundary can have aselective effect on a single species (activator or inhibitor) in the reaction-diffusion system, which will be freeto diffuse across the boundary and persist physically, but without the ability to influence reaction kineticsin this region of the field.1.2 Mixed Boundary ConditionsWe now consider modelling such boundaries in a reaction-diffusion system. Throughout the paper we willconsider a generic two-species RDS in the domain x ∈ Ω ⊂ R n : ∂ u ∂ t = ∇ u + f ( u , v ) , (1) ∂ v ∂ t = D ∇ v + g ( u , v ) , (2)where D > v isa self-inhibitor, D > u and v as presumptiveactivator and inhibitor, respectively.We will show that inhomogeneous steady solutions (Turing-type patterns and otherwise) can be isolatedaway from the domain boundary. This can be achieved by using the boundary conditions, u = R , n · ∇ v =
0, for all x ∈ ∂ Ω , (3) solating Patterns in Open Reaction-Diffusion Systems 5 where R should be approximately a minimum of u in a patterned state (e.g. with Neumann boundaryconditions). Below we will show that the precise value of R seems to play only a small role in influencingthe pattern, at least for some reaction kinetics. In developmental settings, the discussion in the precedingsubsection allows us to consider u and v as behaving differently outside the domain Ω due to boundaries inreceptor distribution or gene expression for example; in turn this can lead to different boundary conditionsfor the activator and inhibitor. We will justify the boundary conditions given in (3) asymptotically in Section2.1 by explicitly considering a spatially heterogeneous extension of the system defined on a larger domain.We also remark that these conditions arise naturally in the case that u represents temperature, and v achemical which undergoes exothermic or endothermic reactions [64, 72]. Throughout we will refer to (3) asmixed boundary conditions, though other authors have used the term nonscalar boundary conditions [15].We will contrast these boundary conditions with homogeneous Neumann conditions on both species, andrefer to the latter simply as Neumann conditions for brevity throughout the manuscript.The rest of the paper is organized as follows. In Section 2 we analytically show how the mixed bound-ary conditions (3) can be obtained from a heterogeneous problem on an enlarged domain. We also describehow these boundary conditions force pattern isolation away from the boundary in terms of a local pattern-selection mechanism. In Section 3, we demonstrate a broad numerical exploration of this interior localiza-tion across a variety of two-species reaction kinetics with qualitatively distinct kinds of patterns. We alsoshow that the same mechanism can operate in three spatial dimensions, as well as in systems with morethan two species. Finally, in Section 4 we conclude by discussing both phenomenological and mechanisticapplications of these boundary conditions, as well as a number of further research directions. We now explore two aspects of the boundary conditions (3). First we explain how they arise from an asymp-totic analysis of a heterogeneous problem, motivated by regions of different gene expression discussed asin Section 1.1. Viewing these conditions as arising from spatial heterogeneity is consistent both with the bi-ological observations discussed there, as well as with the theoretical literature on localization of patterns inheterogeneous media. We then discuss steady state selection mechanisms observed from different bound-ary conditions, explaining some aspects of the isolated patterning we expect to achieve with these boundaryconditions.2.1 Derivation of Mixed Boundary Conditions from a Heterogeneous ProblemWe consider a larger domain ˜ Ω ⊂ R n , and a strict subset of it Ω ⊂ ˜ Ω such that ˜ Ω \ Ω is connected, withall domain boundaries sufficiently smooth. We now consider the non-dimensional heterogeneous reaction-diffusion system, ∂ u ∂ t = ∇ u + (cid:40) f ( u , v ) , for x ∈ Ω , ρ ( R − u ) , for x ∈ ˜ Ω \ Ω , (4) ∂ v ∂ t = D ∇ v + (cid:40) g ( u , v ) , for x ∈ Ω , g ( u , v ) , for x ∈ ˜ Ω \ Ω . (5)For concreteness we consider Neumann conditions for both species on ∂ ˜ Ω .This heterogeneous system can be understood as two regions of different activity (e.g. signal trans-duction and gene expression) in the medium. Within the interior region, Ω , the normal activator-inhibitordynamics are present, whereas in the exterior region, ˜ Ω \ Ω , the activator relaxes to an equilibrium con-centration on a timescale of 1/ ρ , while the inhibitor may interact with the activator in a different manner g with a timescale T g . Such dynamics are easily obtained in examples where gene expression is spatiallymodulated due to previous symmetry breaking and fate specification events, leading to heterogeneity inlocal signalling dynamics or in tissue/cell type as discussed in Section 1.1.One can show that for sufficiently large ρ , the one-dimensional steady state version of Equations (4)-(5)can, to leading asymptotic order, be reduced to considering the steady state problem (1)-(2) on the smaller A. L. Krause et al.
Fig. 2
The geometry of the full domain ˜ Ω , with a region in ˜ Ω \ Ω enlarged to show the various asymptotic parameters and coordinatesalong the boundary discussed in the text. Note that η is in general a ratio of the maximum thickness of ˜ Ω \ Ω to arclength of Γ , but forillustrative purposes we assume the arclength is order unity. Similarly we assume sufficient smoothness of both internal and exteriorboundaries so that the local mean curvature is always order unity or smaller. domain Ω with the boundary conditions (3). In higher dimensions, however, one must impose additionalgeometric constraints. For concreteness, to capture the general case we consider the problem in two spatialdimensions, remarking that extensions to higher dimensions follow the same reasoning. We let Γ = ∂ Ω bethe boundary between the internal and external regions of ˜ Ω (equivalently, the internal boundary of ˜ Ω \ Ω ),and define η as the ratio of the maximum thickness (normal to Γ ) of ˜ Ω \ Ω to the arclength of Γ .We proceed to show that for sufficiently large ρ and for ˜ Ω \ Ω sufficiently slender and regular, steadystates of this system asymptotically reduce to those of the reaction-diffusion system (1)-(3) defined on the in-terior domain Ω . Additionally, the dynamics of u and v outside the interior domain are simple and (asymp-totically) determined analytically. In particular, with η (cid:28) Ω \ Ω , we consider the asymptotic regime ρ (cid:29) η (cid:28)
1. We focus on steady states, as transient solu-tions to the heterogeneous problem (4)-(5) need not satisfy the no-flux boundary condition for v on ∂ Ω inthis asymptotic scaling without further restrictions on kinetic timescales. Our results will then show thatthe steady states of (1)-(3) will be precisely the same as those of (4)-(5). Similarly we will need to exploit thegeometry of a slender domain to prevent problems with transverse gradients in v ; in one spatial dimensionthe restriction of η (cid:28) Γ is related tothe solution in the exterior of Γ on approaching the boundary Γ . To proceed, we first mollify the transitionbetween the two regions with a function H δ which is centred on this boundary and transitions from zeroto unity in the normal direction of the inner boundary of ˜ Ω \ Ω . We also assume Γ has a well-defined andbounded curvature. Furthermore the transition across this boundary is on a lengthscale of δ (cid:28) η .The mollified steady state equations on the extended domain are0 = ∇ u + H δ f ( u , v ) + ( − H δ ) ρ ( R − u ) , 0 = D ∇ v + H δ g ( u , v ) . (6)We introduce an orthonormal curvilinear coordinate system given by ( ζ , ζ ) along the curve Γ , where thelevel set ζ = Γ with the ζ axis perpendicular to Γ and pointing into the interior of ˜ Ω \ Ω .We now follow the derivation of an Eikonal equation along the boundary (see, e.g., [31]). Using the Einsteinsummation convention, the chain rule gives ∂∂ x i = ∂ζ j ∂ x i ∂∂ζ j = : α ij ∂∂ζ j and hence, for example, ∇ v = α ip α iq ∂ v ∂ζ p ∂ζ q + ∂α ip ∂ x i ∂ v ∂ζ p . solating Patterns in Open Reaction-Diffusion Systems 7 We have that the vector ∇ x ζ = (cid:18) ∂ζ ∂ x , ∂ζ ∂ x (cid:19) = ( α , α ) = α ,is normal to Γ and, by the need for orthonormality, we fix the scale of ζ so that α is a unit normal and then κ : = ∇ · α is the mean curvature of Γ . In particular assuming the curvature is of order unity relative to thesmall parameter δ entails the scale of ζ is order unity. Further noting α i α i = ∇ x ζ · ∇ x ζ = ( ζ , ζ ) coordinate system, we have ∇ v = ∂ v ∂ ζ + κ ∂ v ∂ζ + c ∂ v ∂ζ + c ∂ v ∂ζ ,where the coefficients c , c come from the Jacobian of the transformation and are of the same order as κ or smaller. Hence, these terms will be asymptotically small once we rescale ζ below. Derivatives of u willhave an identical expansion.Given the scale of the transition of H δ , we now consider a prospective boundary layer of width of thescale δ at the boundary Γ and thus we consider ζ (cid:48) = ζ / δ , ζ (cid:48) = ζ . Hence at leading order we have u ζ (cid:48) ζ (cid:48) = v ζ (cid:48) ζ (cid:48) = v = A ( ζ (cid:48) ) ζ (cid:48) + B ( ζ (cid:48) ) ,and similarly for u in this region. However boundedness of leaving the transition region forces A to be zero,and thus u = u ( ζ (cid:48) ) = u ( ζ ) , v = v ( ζ (cid:48) ) = v ( ζ ) . The independence of u and v from the normal coordinate, ζ (cid:48) , in the prospective boundary layer, means that on matching into this layer on either side of Γ , we attaincontinuity of u and v . While one cannot match derivatives in the same manner in general, the collapseof the prospective boundary layer entails the boundary condition of continuity should be supplementedby continuity of flux. To explicitly deduce this, let z such that | ζ − z | ≤ δ be fixed and consider theintegral form of the steady state conservation relation (given by (6)) on the pillbox region P = { ( ζ , ζ ) | ζ ∈ [ − δ , δ ] , ζ − z ∈ [ − δ , δ } ] .Integrating the conservation equation (6) for v over the pillbox and applying Green’s Theorem gives0 = D (cid:90) z + δ z − δ − ∂ v ∂ζ ( − δ , ¯ ζ ) d ¯ ζ + D (cid:90) z + δ z − δ ∂ v ∂ζ ( δ , ¯ ζ ) d ¯ ζ + O ( δ ) + O ( ρδ ) , (7)and similarly for u . In particular, the absence of a boundary layer at Γ entails that u , v are bounded, withbounded normal derivatives with bounds independent of δ . Thus the O ( δ ) terms emerge from the other fluxboundary integrals and the O ( ρδ ) term emerges from the integrals of sinks and sources over the pillbox,noting that the area of the pillbox P is 4 δ and the largest source term scales with ρ (cid:29)
1. However, theasymptotic parameter δ is from the mollification of a Heaviside function and thus we can take δ ρ ∼ O ( δ ) ,since δ has been introduced for mathematical convenience and can be taken as small as required, whereas ρ is constrained by the underlying biophysical limits of source production in the interpretation of the model.Hence, using the integral mean value theorem on the above integrals, dividing by δ and taking the limit δ →
0, we have continuity of flux, and thus normal derivative, across Γ for v , and analogously for u , asexpected.Knowing how to match across Γ we now need to consider our next objective, which is to determinethe behaviour of the solutions within the region ˜ Ω \ Ω to assign boundary conditions for a reduced modelwithin the boundary Γ . We take advantage of the slenderness of ˜ Ω \ Ω , characterised by an aspect ratio ofits normal extent to the arclength of Γ via η (cid:28)
1. It is most useful to see if non-trivial behaviour may occurin the normal direction as, if present, such behaviour would preclude a homogeneous Neumann conditionfor v in the reduced system. We thus rescale ζ ∗ = ζ / η , ζ ∗ = ζ , and at leading order in η we have D (cid:16) η − v ζ ∗ ζ ∗ + O (cid:16) η − (cid:17)(cid:17) = g ( u , v ) . A. L. Krause et al.
If the characteristic time of the inhibitor kinetics in the exterior, T g , is large when compared to that ofdiffusion, η D − (cid:28) T g , we finally have, from (6), v ζ ∗ ζ ∗ = Γ by a scale of more than δ (cid:28) η . We have also assumed that the domain of ˜ Ω \ Ω is sufficiently regular to ensure κ ∼ o ( η ) , which will occur if the radius of curvature is on the lengthscaleof the perimeter of ˜ Ω \ Ω , rather than the smaller lengthscale of its thickness.We now need to make a further regularity assumption about ˜ Ω \ Ω , namely that its external boundarycan be specified by ζ ∗ = F ( ζ ∗ ) where F ∼ O ( ) by construction, and is also a single-valued function thatsatisfies η | F (cid:48) ( ζ ∗ ) | (cid:28)
1, where prime denotes derivative. Hence we do not consider an external boundaryconsisting of a curve that, for example, has a cusp or varies too rapidly. From this assumption, we have thatthe normal derivative operator on the exterior boundary of ˜ Ω \ Ω is proportional to ∂∂ζ ∗ − η F (cid:48) ( ζ ∗ ) ∂∂ζ ∗ ≈ ∂∂ζ ∗ .Thus, working at leading order in η | F (cid:48) ( ζ ) | (cid:28)
1, homogeneous Neumann conditions on the exterior bound-ary of ˜ Ω \ Ω give v ζ ∗ = v = v ( ζ ∗ ) = v ( ζ ) throughout ˜ Ω \ Ω ,except possibly δ -close to Γ .As will emerge below, we ultimately consider the leading order behaviour when the asymptotic pa-rameters satisfy 1/ δ (cid:29) ρ (cid:29) η (cid:29)
1. The restriction 1/ δ (cid:29) ρ entails the scale of the mollificationregion, δ , is the smallest scale present as required since the mollification is introduced for analytical ex-pedience, rather than as a fundamental feature of the biophysics. The constraint ρ (cid:29) η ensures thesource strength is sufficiently large to enforce approximately homogeneous solutions in the region ˜ Ω \ Ω ,as observed below. Finally η (cid:28) Ω \ Ω to be slender. To proceed, we consider u in the region ˜ Ω \ Ω , whereby atleading order 0 = u ζ ∗ ζ ∗ + η ρ ( R − u ) , (8)once away from Γ by a scale of more than δ . Noting the Neumann boundary condition on the exteriorboundary, Equation (8) gives u ( ζ ∗ , ζ ∗ ) = R + A ( ζ ∗ ) cosh (cid:16) ρ η ( F ( ζ ∗ ) − ζ ∗ ) (cid:17) . (9)For δ sufficiently small, the above solution will only vary by an asymptotically small amount across the δ -scale mollification region, as there is no boundary layer. Hence, as ζ ∗ →
0, this solution is continuous withthe solution on Ω as the latter approaches Γ for the same value of ζ ∗ , and analogously for the flux. No usefulinformation is gained by continuity of concentration. However, continuity of flux entails ∂ u / ∂ζ ∼ O ( ) as ρ increases on approaching Γ , assuming ∂ u / ∂ζ in Ω does not blow up with increasing ρ , which is consistentwith numerical simulations (below) and parabolic regularity. We then have that Equation (9) yields ∂ u ∂ζ (cid:12)(cid:12)(cid:12)(cid:12) ζ = + = ∂ζ ∗ ∂ζ ∂ u ∂ζ ∗ (cid:12)(cid:12)(cid:12)(cid:12) ζ ∗ = + = − ρ A ( ζ ∗ ) sinh ( ρ η F ( ζ ∗ )) ∼ O ( ) , (10)and hence | A ( ζ ∗ ) | ∼ O (cid:18) ρ ( ρ η F ( ζ ∗ )) (cid:19) . (11)We further have that ρ is sufficiently large to ensure ρ (cid:29) η , and we have already taken η (cid:28) ( p ) is monotonically increasing for p >
0, and F ( ζ ∗ ) ≥ ζ ∗ ≥ F ( ζ ∗ ) > Ω \ Ω , wethen have for solution (9) that | A ( ζ ∗ ) cosh (cid:16) ρ η ( F ( ζ ∗ ) − ζ ∗ ) (cid:17) | ≤ | A ( ζ ∗ ) | cosh (cid:16) ρ η F ( ζ ∗ ) (cid:17) ∼ O (cid:18) ρ coth ( ρ η F ( ζ ∗ )) (cid:19) ∼ O (cid:18) ρ (cid:19) . (12) solating Patterns in Open Reaction-Diffusion Systems 9 Thus, providing R (cid:29) ρ , we have u = R is an estimate of the solution throughout ˜ Ω \ Ω with asymp-totically small relative error. We typically take R ∼ O ( ) so that small relative error in the approximation of u = R in the region ˜ Ω \ Ω is assured by our assumption of ρ (cid:29) η (cid:29)
1. In addition, if R ∼ O ( ρ ) we have u ∼ O ( ρ ) , so that the estimate u = R presents with asymptotically small absolute error, eventhough relative error is of order unity or possibly larger.Finally, we consider the effective boundary conditions for the steady states of the reduced system (1)-(2),0 = ∇ u + f ( u , v ) , 0 = D ∇ v + g ( u , v ) , (13)on the domain Ω , which has external boundary Γ with ρ (cid:29) η (cid:29) Ω \ Ω .Continuity of flux immediately gives homogeneous Neumann conditions for v , while continuity is inap-propriate for v , since v is not determined in ˜ Ω \ Ω , except ultimately via its coupling to the solution onthe interior of Γ . Furthermore, continuity at Γ gives a boundary condition of u = R on Γ for Equation(13). In contrast, continuity of flux for u on approaching Γ is inappropriate as it will generate a homoge-neous Neumann boundary condition which can lead to behaviours that differ from the original system,Equations (4)-(5) in general, as no-flux boundary conditions lead to a superset of solutions containing thoseapproaching constant values as a special case. We will explore this idea further in the next subsection, refer-ring specifically to Fig. 4 as a demonstration of this solution selection mechanism under different boundaryconditions.In summary, we have that, providing ˜ Ω \ Ω is sufficiently slender and regular, and with ρ sufficientlylarge, the behaviour of steady solutions of Equations (4)-(5) on ˜ Ω is given, to asymptotic accuracy, by so-lutions of Equations (13) within the domain Ω with the boundary conditions on Γ that u = R and that thenormal derivative of v , that is ∂ v / ∂ n , is zero. We note that the set of steady state solutions will match be-tween the heterogeneous and reduced problems, but solutions to the time dependent problems could allowfor solution selection based on transients, leading to a disagreement between the two models. Nevertheless,we expect this to be rare or only have small effects in typical cases, which we demonstrate numerically inSection 3. Finally we also remark that the asymptotic restrictions, particularly regarding the slender geom-etry, can also be observed numerically to have typically little effect on stationary solution behaviours whencomparing these two models, as we will discuss from full numerical simulations in Section 3 for specificnonlinear kinetics. The driving force in deriving the conditions (3) is the heterogeneous switching betweeninteracting and weakly or non-interacting species, which can be understood in the context of the biologydescribed in Section 1.1.2.2 Interior Patterning as Steady State SelectionHere we further elucidate why these boundary conditions will lead to confinement of patterned states onthe interior of the domain, using a particular example of reaction kinetics in a one-dimensional domain. Weconsider patterned steady states as subsets of all admissible ones in the periodic case, and show how specificsubsets of this set are selected by different boundary conditions. The motivation here is the discussion ofequivariant bifurcation theory and the impact of symmetry on solution branches in [15]. Specifically, wenumerically show that there are typically many solutions that satisfy Neumann or periodic conditions for u and v , but only a subset of these satisfy the mixed conditions (3) when R is chosen to coincide with a specificminimum value of u from the Neumann solutions. In particular, solutions with these mixed boundaryconditions are forced to pattern away from the boundaries as they are exactly chosen to match a Neumannsolution with patterning only on the interior. We numerically observe that small changes in R from thisvalue lead to small local changes at the boundary, but broadly similar solution structures overall.To demonstrate this concretely, we simulate the Schnakenberg kinetics from Table 1 (though with differ-ent domain lengths, L ) using periodic, Neumann, and mixed boundary conditions. We employ a method-of-lines approach to simulate the RDS, using the standard three-point stencil to discretize the Laplacian.We use the Matlab function ‘ode15s’ to evolve these ordinary differential equations in time, and the Matlabfunction ‘fsolve’ to compute steady state solutions to the discretized system on n = t = units of time, andthen use its value to find a nearby stable steady state (which is always almost identical to the last transientstate, suggesting it is reachable from random initial data). From this steady state solution, we generate nine (a) Periodic (b) Neumann (c) Mixed, R = min ( u Neumann ) (d) Mixed, R = Fig. 3
Steady state solutions from simulations of the Schnakenberg kinetics in Table 1 under four choices of boundary conditions. Tensimulations are shown where in (a), the initial simulation uses random initial data, but each subsequent simulation is initialized bycyclically shifting this solution 5% along the length L , and a stable steady state is found (by cyclic symmetry, these are automaticallystable steady states as the first one was). For (b)-(d), the initial data are taken as each of the steady states in (a), evolved for t = a = b = c = D =
20, and domain length L = others, shown in Figure 3(a), by shifting this steady state pattern cyclically. By the translational invarianceof the Laplacian with periodic boundary conditions, any stable shifted steady state solution with any shiftis still a stable steady state solution, and we confirm this numerically.We use these ten shifted periodic solutions as initial conditions for simulating the RDS using threeother sets of boundary conditions, and then again find stable steady states after evolution in time. We useNeumann conditions in Figure 3(b), mixed conditions given in (3) with R equal to the minimal value of u from the Neumann case in (c), and conditions (3) with R = R setas the minimum value of u with Neumann conditions select precisely one solution from the two Neumannsolutions (and no others). Finally, using R = solating Patterns in Open Reaction-Diffusion Systems 11 (a) Periodic (b) Neumann (c) Mixed, R = min ( u Neumann ∩ periodic ) (d) Mixed, R = Fig. 4
Steady state solutions from simulations of the Schnakenberg kinetics in Table 1 under four choices of boundary conditions,exactly as in Figure 3 except taking a domain length of L =
20. Note that in (b), the solutions in red and blue have different amplitudesthan those in yellow and green. in Figure 4(b) do not correspond to solutions with periodic boundary conditions shown in (a), possessinga different amplitude. Setting R as minimal value of u from these Neumann conditions, and restricting tothe set which satisfy periodic solutions, we obtain a single solution with mixed conditions in (c). Setting R equal to either the other minimum of u from the Neumann conditions, or to 0, has only a small influenceon the solution near the boundary, as can be seen in (d). The mixed boundary conditions further excludehalf-integer mode solutions.The reduction of the admissible set of patterns has only been shown numerically for this particularexample, and one would need to employ more sophisticated mathematical techniques to demonstrate thisfor all kinetics and parameters, and beyond one-dimensional examples. Nevertheless, the PDEs in questionare local, and so any solution where v is extremal and u is minimal at the same points which satisfy bothNeumann and periodic boundary conditions will also satisfy the conditions (3) when R is suitably chosen(i.e. exactly or approximately equal to a minimum of a Neumann solution). It is not always the case thatreaction-diffusion patterns will have both species sharing extremal points, but this can be shown near aTuring bifurcation via linearization, and is observed (at least approximately) in many systems numericallybeyond the bifurcation [15]. The pattern selection mechanism described here also suggests that peaks of theactivator u should be approximately half of a wavelength away from the boundary, and we will directlytest this prediction across a range of parameters and kinetics in Section 3.1. f g R L D a b c d e SCH (Spots) a − cu + u v b − u v a − cu + u v b − u v a + u v − bu u − cv a − u − buv + u + cu d ( e − v ) − buv + u + cu a − u − buv + u + cu d ( e − v ) − buv + u + cu u − au + v − b c − u − dv -1.6 250 50 0.33 0.6 0.6 0.99 1 Table 1
A list of reaction kinetics and non-dimensional parameter regimes simulated, where L is a domain length. SCH: Schnakenberg[62]; GM: Gierer-Meinhardt [24]; TH: Thomas [32]; FHN: FitzHugh-Nagumo [22, 23, 52]. Note that the R values chosen here areroughly close to minimal values of typical patterns, but not exactly these values. Note that the FitzHugh-Nagumo kinetics modelvoltage potential relative to a ground state via the variable u , and so need not admit positive solutions of this variable. We now give example simulations of RDS with the boundary conditions (3), as well as those with Neu-mann boundary conditions. Additionally, we compare these with simulations of the heterogeneous prob-lem (4)-(5) with Neumann conditions on the outer boundary, to demonstrate that our asymptotic reduc-tion to the mixed boundary conditions is valid for large values of the relaxation rate ρ and suitably thinextended domains. Finally we will also give example simulations in more complicated geometries andhigher-dimensional domains, as well as systems with more than two species, showing that the isolatedpatterning emergent from our mixed boundary conditions extends beyond these initial examples.For all of these simulations, we used the commercial finite-element software COMSOL. In the two-dimensional simulations, we used a minimum of 5 × triangular second-order finite elements, and inthe three-dimensional simulations we used at least 10 tetrahedral elements. We always chose parameterssuch that there was a unique homogeneous steady state for the kinetics, given by f and g for homogeneousNeumann boundary conditions , and we perturbed this state by multiplying it by a normally distributedrandom variable of standard deviation 10 − and unit mean, sampled identically and independently at eachfinite element. All simulations shown throughout this section are given at t = units of time, by whichpoint they were within numerical tolerances of a stationary solution. The timestepping was implementedusing a standard adaptive backwards-differentiation formula of orders 1 through 5. Refinements in spacewere used to check convergence for particular simulations, and finite-difference simulations were carriedout in Matlab (using the five-point Laplacian stencil) to check convergence to the same steady state patternfrom identical initial data.3.1 Interior Patterning across Two-Species KineticsWe now describe results from simulations of each of the reaction kinetics and parameter regimes shown inTable 1. These are given in Table 2, where the first column indicates the parameter set, the second columnsimulations with Neumann boundary conditions, the third column simulations with conditions (3), andfinally the fourth column shows simulations of the enlarged heterogeneous system (4)-(5). For each of thesesimulations we used the same realization of the random initial condition, as all of these systems can admitmultistability of different inhomogeneous steady states depending on the initial data [4, 29], and results aredisplayed for sufficiently long times to ensure transient behaviour has relaxed onto steady states.Comparing the Neumann simulations with those using our mixed boundary conditions in Table 2, wecan clearly see that the mixed boundary conditions lead to the inhomogeneous patterns being confined tothe interior of the domain. Additionally, it is visually apparent in each case that the distance from the peakof an outermost structure (spot or stripe) to the boundary is approximately half of the wavelength betweeninterior structures, despite the fact that R was not set to be exactly the minimum pattern value as in Sec-tion 2.2. Finally, as the same initial perturbations were used across all simulations, comparing the mixed Note that such an initial condition need not satisfy the boundary conditions (3), but will effectively instantly relax to solutionssatisfying these conditions at the boundary under time evolution, as the system is parabolic.solating Patterns in Open Reaction-Diffusion Systems 13 boundary conditions to the heterogeneous ones, we see essentially identical steady state solutions in allcases except for the FitzHugh-Nagumo kinetics, which admit small differences between the steady statesselected. We also considered how the asymptotic reduction holds for different domain geometries and forsmaller values of ρ , but for brevity omit this analysis (largely as it depends heavily on the kinetics and typeof pattern studied). Essentially these results confirm our analysis in Section 2.1 that the mixed conditionsapproximate such heterogeneous problems at least once any transients have relaxed. We note that theseboundaries have discontinuous derivatives at the corners, but simulations on circular domains have identi-cal properties. Similarly, simulations with much larger external regions (e.g. when ˜ Ω is [ − L , 2 L ] × [ − L , 2 L ] ,so that the domain length is three times larger) give steady states with interior patterning, as in those onthe interior domain with mixed boundary conditions in Table 2. This suggests that, at least in some cases,the geometric assumptions on the exterior domain can be relaxed, though we will not pursue this furtherhere.These boundary conditions have some influence on the structure of the resulting patterns beyond theboundaries, especially in the case of the labyrinthine patterns shown in Table 2. The stripes in the Schnaken-berg example show a clear confinement due to the square geometry. In the labyrinthine case for Thomaskinetics, we see that the (roughly) hexagonal arrangement of inverted spots in the Neumann case is forcedinto a square arrangement in the centre of the domain with the mixed conditions, with a stripe-like structurenow surrounding this region. The FitzHugh-Nagumo structures are qualitatively similar on the interior ofthe domain, with random-looking labyrinthine stripes, but the boundary is clearly marked with a stripedregion having some local effects on the meandering of the stripes inside. As these labyrinthine solutionsare more global in structure than localized spots, we expect less isotropic patterns and more of an influenceof the confining geometry and initial conditions. This is true even in the Neumann case, as we can see bythe orientation of the Schnakenberg stripes aligning in a large part of the domain with one of the axes. Wealso remark that the alignment of labyrinthine patterns with Neumann boundary conditions can be sen-sitive to initial conditions, and one can obtain quantitatively large differences for different realizations ofrandom initial conditions. In contrast, the Schnakenberg and Thomas kinetics with our mixed boundaryconditions gave consistently the same final steady state seemingly independent of a variety of initial condi-tions (though the positioning of the interior structures of the FitzHugh-Nagumo steady states did dependon initial conditions). Depending on the application, robust alignment may or may not be desirable. Wewill further explore interior patterning and pattern selection effects in the next section by considering morecomplicated domains.3.2 Higher-Dimensions & Complex GeometriesWe now demonstrate the effects of more complicated domain geometries on patterns with the mixedboundary conditions (3). As seen in the simulations in the previous section, both spot and labyrinthinepatterns were more ordered with these conditions, compared to the Neumann case. In particular, thelabyrinthine patterns shown in Table 2 conform to the geometry in the case of the mixed conditions (3),whereas the Neumann patterns were seemingly less affected by the boundary, having a broadly moreisotropic character. We now show how these labyrinthine patterns are influenced by more complicateddomain boundaries with Neumann and our mixed conditions. We then demonstrate that interior confine-ment, as well as these pattern selection effects, extend to higher spatial dimensions.Throughout these examples, we will use the Schnakenberg kinetics from Table 1, and primarily focus onparameter regimes which give rise to labyrinthine solutions, comparing qualitative features such as stripeorientation and defects. We remark that for periodic boundary conditions and small-amplitude patternsnear the onset of a Turing instability, there is a dichotomy between stripe and spot solutions [18], deter-mined by signs and relative magnitudes of quadratic and cubic interactions. More recent work has shownthat such a simple classification does not hold beyond the weakly nonlinear regime, and qualitatively dif-ferent kinds of patterns can be simultaneously stable or interact [4, 29, 5]. We observe here that the largestqualitative differences are determined by the geometry (and the initial data in the Neumann setting) in thelabyrinthine cases. We focus on defects (regions where the stripes are no longer contiguous at the highestvalue of u , which appear like dappled spots), in addition to stripe orientation. S C H ( S p o t s ) S C H ( L a b y r i n t h i n e ) G M T H ( S p o t s ) T H ( L a b r y n t h i n e ) F H N Table 2
Concentrations of u from simulations of the kinetics and parameters given in Table 1 with Neumann conditions, the mixedboundary conditions (3), and in a heterogeneous system given by (4)-(5) with ρ = and Neumann conditions on the exteriorboundary. Maximal values of u are in the light yellow regions, and minimal values in the dark blue, though precise numerical valuesvary between kinetics and parameters, which are specified in full in Table 1. The domain in columns two and three was the square Ω = [ L ] , and the heterogeneous domain was the enlarged square ˜ Ω = [ − L , 1.05 L ] , using the same interior domain Ω asdescribed near Equations (4)-(5). This final column has an interior domain Ω corresponding to the black square, where the boundariesin the second and third columns are located. For each set of kinetics, the colour scale is fixed across the row, and the approximateminimum/maximum (blue/yellow) values are as follows: SCH (Spots) 0/23.56; SCH (Labyrinthine) 0/3.58; GM 0/163.9; TH (Spots)1.56/37.4; TH (Stripes) 2/16.2; FHN -1.65/1.65solating Patterns in Open Reaction-Diffusion Systems 15(a) Neumann (b) Neumann (c) Neumann(d) Mixed (e) Mixed (f) Mixed Fig. 5
Values of u computed with Neumann and the mixed boundary conditions (3) on the interior of circles and ellipses. Simulationsshown are on a domain with a semi-minor axis of 40, and semi-major axes of 40 in (a) and (d), 80 in (b) and (e), and 120 in (c) and (f).We used the Schnakenberg kinetics in Table 1 with a = b = c = D =
20, and R = First we consider circular and elliptical domains in Figure 5. As anticipated, solutions with Neumannboundary conditions in Figure 5(a)-(c) exhibit orientations which are seemingly independent of the domaingeometry, though there is some correlation of stripe orientation. Stripe defects appear in all cases withNeumann boundary conditions, mostly near the boundaries where stripes of different orientation intersect.In comparison, there is clearly more symmetry in the confined patterns generated by mixed boundaryconditions in Figure 5(d)-(f), which broadly maintain orientations consistent with the domain geometry.Defects with mixed boundary conditions only occur in the elliptical cases, and seem to be where the stripecurvature is greatest.In each of the cases with Neumann boundary conditions, there are slightly different numbers of stripes(depending on what one classifies as a stripe), with (c) having the largest number for these simulations.Different random initial conditions can lead to substantial variation in the number of such structures, de-pending on what orientation most stripes take. In contrast, if one considers a contiguous region of higher-than-average values of u as a stripe, then the confined patterns generated by mixed boundary conditionsall have exactly four of these. This is despite the fact that the largest ellipse occupies three times the areaof the circle, and these observations appear to be robust to different random perturbations of the initialconditions. Finally, these results are relatively insensitive to the value of R , as simulations with R = R = u in the Neumann cases was u ≈ x ( s ) = L cos ( s ) (cid:113) + γ sin ( s ) , y ( s ) = L sin ( s ) (cid:113) + γ sin ( s ) , for s ∈ [
0, 2 π ] , (14)where L is now a measure of the domain length (the radius of the circle for γ = γ ∈ [
0, 1 ) signifiesthe deviation from this circle. The sin ( s ) term makes the perturbations resemble a 6-lobed polar rose orrhodonea curve.In Figure 6 we demonstrate patterns on such domains for three different values of γ . As in the ellipti-cal case, the solutions with only Neumann conditions in (a)-(c) have stripes with similar alignments, anddefects when they meet stripes with different orientations. Very little direct effect of the boundaries canbe seen in these simulations, with the patterns broadly consisting of seemingly randomly oriented stripes(with orientations depending on initial data). In contrast, the confined patterns in Figure 6 clearly conform γ = γ = γ = γ = γ = γ = Fig. 6
Values of u computed with Neumann boundary conditions, as well as the mixed boundary conditions (3) on the interior of thedomain given parametrically in (14). Simulations shown are on a domain with size characterized by L =
70 using the Schnakenbergkinetics in Table 1 with a = b = c = D =
20, and, in the mixed case, R =
0. Colour scales were fixed between a minimum(blue) of 0 and a maximum (yellow) of 2.91. to the boundary, with defects again appearing at regions where the stripes (but not necessarily the bound-aries) might be expected to have large curvature. As before, while the structure of the Neumann solutionscan vary tremendously depending on the initial data, the mixed boundary conditions lead to consistentpatterns even when the initial data are different. In the centre of each of these domains, independent of thevalue of γ , is a pattern of seven inverted spots (one surrounded by six others). These simulations demon-strate an intriguing possibility of selecting for hybrid stripe and spot solutions in the domain throughthe discrete symmetry of the geometry. As in the elliptical case, qualitatively similar patterns are seen for R = R = u .In Figure 7, we give examples of localized sphere-like structures which emerge again in the Schnaken-berg kinetics in a rectangular domain. As expected from the two-dimensional examples, spheres in theNeumann case form hemispheres and quadrants (quarter-spheres), along the boundary. In contrast, themixed boundary conditions lead to symmetrical spheres of the solution u along the centre of the domain,approximately half of a wavelength away from the boundary in any direction as in the two-dimensionalcase. As in the case of spots in the two-dimensional examples given in Table 2, it is clear that the patternswith these boundary conditions are substantially more ordered and symmetrical. Additionally, we find a solating Patterns in Open Reaction-Diffusion Systems 17(a) Neumann (b) Mixed Fig. 7
Thresholded values of u from simulations of (1)-(2) using purely Neumann boundary conditions in (a), and the mixed conditions(3) in (b). Triangular elements were only shown for values of u > × ×
60 using the Schnakenberg kinetics in Table 1 with a = b = c = D = , and, inthe mixed case, R =
0. The maximum value of u (yellow) is 23.56, and the minimum (blue) is 0. robustness of pattern structure across different random initial conditions which contrasts with the case ofNeumann boundary conditions.In Figure 8, we give an example where more lamella-like structures (analogous to stripes in two-dimensions)appear. These patterns are more complicated and harder to visualize, so cross-sections have been includedto help picture them beyond the regions of high activator concentration. The Neumann case admits com-plicated tube-like structures inside the domain, as well as forming partial tubes along the boundary. In con-trast, the confined patterns resemble a face-centred cubic lattice with tubular connections between sphere-like regions of high activator concentration. This example clearly shows a difference between the seeminglymore random appearance of patterns in the Neumann case, compared to the extremely structured confinedpatterns arising from the mixed boundary conditions.Lastly in this section, we consider domains which are composed of both Neumann and the mixedboundary conditions (3). Such a setting may arise from cutting a piece of tissue from a patterning fieldand then inserting an impermeable material, with local no-flux boundaries along the cut. In Figure 9, wegive examples of spots and stripes in the Thomas kinetics from Table 2, where we have shown rectangularand elliptical ‘cuts’ in the domain along the left boundary, so that these boundaries have Neumann con-ditions but all others have the mixed conditions. We see half-spots forming only along these Neumannboundaries in panels (a) and (b). In panels (c) and (d) we observe a more pronounced selection effect ofthe Neumann boundary conditions, where the domain with the small rectangle removed has an internalstructure consistent with the corresponding simulation from Table 2 (though with patterns forming up tothe Neumann boundary), but the elliptical cut induces a selection of horizontal stripe patterns even awayfrom this boundary. Such a selection mechanism could be a useful approach to validating these boundaryconditions in experimental systems.3.3 Extensions to Many-SpeciesInterior pattern localization from mixed boundary conditions can be extended beyond two-species systems.Pattern formation in three and more species systems is far richer than in the two-species case, even whenrestricting to patterns arising from Turing instabilities [56, 61], and such multispecies kinetics have beenstudied extensively in recent years in developmental contexts [34, 14, 63]. We will give a brief example,motivated by pattern formation in three-species Lotka-Volterra systems. We note that two-species Lotka-Volterra systems do not admit patterns except in the case of convex domains and bistability [33, 42], thoughthree-species models do [69].We consider the system, ∂ u i ∂ t = D i ∇ u i + u i (cid:32) − ∑ j = a ij u j (cid:33) , for i =
1, 2, 3, (15)
Fig. 8
Thresholded values of u from simulations of (1)-(2) using Neumann boundary conditions in (a) (with slices shown in (b)), andthe mixed conditions (3) in (c) (with slices shown in (d)). For the plots in (a) and (c), triangular elements were only shown for valuesof u > × ×
25 using theSchnakenberg kinetics in Table 1 with a = b = c = D =
30, and, in the mixed case, R =
0. The maximum value of u (yellow)is 3.58, and the minimum (blue) is 0. where D i > a ij ∈ R are the interaction coefficients.This is a particular nondimensionalisation of a standard generalised Lotka-Volterra model with diffusionmodelling random dispersal. The form of interactions between species permits competition for resources,as well as predation from a generalist or intraguild predator who is also competing for resources with theprey; see [69] for more details about the ecological interpretation. We will use the following extension of thetwo-species mixed boundary conditions (3), u ( x , t ) = R , n · ∇ u = n · ∇ u =
0, for all x ∈ ∂ Ω . (16)We will also consider homogeneous Neumann conditions (on all species) for comparison, again simplyreferring to these as Neumann conditions.We show simulations of this system with parameters giving spots and labyrinthine patterns in Figure10 for both Neumann boundary conditions, and the mixed boundary conditions (16). As in the two speciessimulations above, we see the mixed boundary conditions forcing the patterned regions of u into the inte-rior of the domain, and some pattern modulation due to nonlinear interactions between patterned regions. solating Patterns in Open Reaction-Diffusion Systems 19(a) TH (Spots) (b) TH (spots)(c) TH (stripes) (d) TH (Stripes) Fig. 9
Values of u from simulations of (1)-(2). The domain is the same as in Table 2 except that a region of the domain has been removedfrom the left boundary. In (a) and (c), a rectangular region of height 0.7 L and width 0.05 L has been removed from the domain, and in(b) and (d) an elliptical region of semi-major axis of 0.35 L and semi-minor axis 0.2 L has been removed. We use the mixed conditions(3) along all of the square boundaries, and Neumann conditions along the three rectangular cut boundaries in (a) and (c), and theelliptical cut boundaries in (b) and (d). Simulations shown are on a domain of size L =
60 using the Thomas kinetics and parametersfrom Table 1 as labelled, with maximal and minimal values of u as in Table 2. For instance, in (b) we see a few smaller spots due to crowding, and in (d) we see a clear mode selectionfrom the boundaries as shown before in the labyrinthine examples. Of course both of these kinds of pat-terns can also occur with Neumann conditions, but may be more prevalent for mixed boundary conditionsas they confine patterns to the interior. Besides these influences, the patterns are qualitatively similar onthe interior of the domain, especially between (a) and (b). Finally, we note that simulations in an extendeddomain, analogous to the system (4)-(5), gave the same results as using (16) (not shown).While this is only a single example of a multispecies model, we anticipate that the generalisation to n species will have qualitatively similar effects. Without loss of generality, the arguments given in Section 2.1can be extended to the n species case as long as all constraints are satisfied. In particular, we anticipate thatplacing the Dirichlet condition on any single species will confine the patterns to the interior via the samemechanism suggested in Section (2.2). Isolation of patterns in reaction-diffusion systems (RDS) on the interior of a domain can be motivated bothphenomenologically, by observing patterning fields which are clearly subsets of a given region, or mech-anistically in terms of boundaries in gene expression. Such gene expression boundaries themselves caneither emerge due to previous cell fate specification or local signalling dynamics. Here we have developeda novel method of inducing such isolated patterns with a simple change in the boundary conditions gov-erning the RDS. We have justified these boundary conditions by considering heterogeneous RDS modellinga change in reaction dynamics explicitly, and shown that such systems can reduce asymptotically to RDSsatisfying these mixed boundary conditions at steady state, given weak geometric and kinetic constraints.
Fig. 10
Values of u from simulations of (15) with either Neumann or the mixed boundary conditions (16). We used a = a = a = a = a = a = a = − a = D = D = D =
1, with a = a = Ω was taken to be a square domain of side length L = u is in phase with u , and u is out of phase with these. The colourscorrespond to a maximal value of 0.7 (light yellow) and a minimal value of 0 (dark blue). Simulations are shown at time t = , andwere initialized using normal perturbations of the homogeneous steady state as described in the beginning of the section. See Figure1 of [69] for an ecological description and interpretation of the parameters. We took R = This is consistent with the theoretical literature showing pattern localization and modulation in heteroge-neous environments [74, 54, 55, 40, 37], as well as the regionalisation of patterning seen experimentally [30].These conditions give a simple way of modelling these localization phenomena, without having to resortto modelling the heterogeneity that may have given rise to isolated patterning regions. We have furthershown that isolated patterning can be obtained quite generically using these mixed boundary conditionsacross a range of reaction-diffusion systems, geometries, spatial dimensions, and even in systems of morethan two interacting species.As described in the introduction, there are still major gaps in our understanding of RDS, particularlyas models in developmental biology. While the original formulation of Turing’s theory is a powerful andsimple way of obtaining periodically patterned states, it suffers from a number of robustness problems, inpart due to this simplicity [46, 79, 63]. In addition to demonstrating isolated patterning, our results indicatesome level of robustness in pattern selection, partly explained in Section 2.2. Specifically, the numericallyobserved steady states in the case of homogeneous Neumann conditions can generally depend on initialdata, as multistability of steady states is common, especially in two or more spatial dimensions [4]. Incontrast, there were far fewer such cases of multistability of patterns observed from simulations with themixed boundary conditions (3). Of course more work is needed to precisely explain this phenomenon, but itis consistent with the discussion of how bifurcation branches of solutions change under different boundaryconditions in [15]. Additionally, while we did not explore this in great detail here, the Turing space forthese mixed boundary conditions can be enlarged, as discussed in [45], plausibly helping to overcome theconstraints in finding parameters which admit patterned solutions.While we focused on analyzing terminal steady state patterns, we remark that the boundary conditions(3) often had transient impacts which may be relevant in applications. In particular, in every simulationinvestigated, we observed that these boundary conditions led to a substantially faster convergence towardsthe steady state pattern. In the case of labyrinthine patterns, these conditions often led to orders of magni-tude faster convergence. Specifically, for the Schnakenberg kinetics with labyrinthine parameters in Table1, the solution with conditions (3) was within plotting accuracy of the final steady state within O ( ) unitsof time, whereas the solution with homogeneous Neumann conditions only attained the final stripe orien-tations shown in Table 2 after O ( ) units of time. We conjecture that this is not, at least primarily, dueto a change in linear stability of the final pattern, but is due to fewer steady states satisfying these bound-ary conditions, and hence less shadowing of unstable equilibria. As with the above discussion of countingequilibria, we leave further investigation to future work.Finally, we mention that our results also demonstrate a novel application of domain geometry in ro-bustly giving rise to patterns of qualitatively different types. Specifically, in Figure 6 we see that invertedspots can be selected on the interior of these domains, with stripes closer to the boundaries, due to thesymmetry of the domain geometry. As far as we are aware this is a novel local pattern selection mechanism solating Patterns in Open Reaction-Diffusion Systems 21 quite distinct from those exploiting a reduction of dimension to transition from spots to stripes [49, 50].In general the sensitivity of patterns with respect to geometry, and the use of domain geometry to influ-ence patterning, have not been studied as extensively as other aspects of RDS, especially for more generalboundary conditions.Many periodic patterns in development, such as mammalian hair, are spatially localized within someregion of an otherwise homogeneous tissue. There are at least three distinct ways that such regionalisationcan occur, described in Section 1.1. Determining precisely what is happening at such boundaries is difficult,and often their exact geometry is unknown. If, for example, a periodic pattern of spots forms in one part ofa field but not another, it does not immediately reveal where the functionally relevant boundary might be,whether it is graded or sharp, whether it is relevant for activator and inhibitor equally or selective for onespecies, and whether it promotes the formation of spots nearby (and so the boundary lies close to the edgeof the spotted area) or repels them (and so it lies at some distance from the margin of the spotted area).The boundary conditions presented here give a simple way to model such phenomena, when the boundaryis assumed to be repelling, without having precise knowledge of what the field is doing near or beyondthe boundaries where periodic patterns are observed. The examples given in Section 3.2, and in particularthe ‘cut’ domains shown in Figure 9, provide ways of validating these isolating boundary conditions viaexperiments which alter domain geometry, as we have demonstrated strong effects of such geometry onthe emergent patterns.There are numerous extensions of the basic ideas here, but we restrict attention to five possible areasof future work. As above, there is work to be done in more rigorously justifying our conclusions regard-ing equilibria, pattern selection, and transient dynamics. On the biological side, exploring patterning fieldboundaries with these kinds of models in mind could help us understand the key aspects of the detailedbiological mechanisms at play. There are also important questions, raised in [19] and elsewhere, about theunderlying thermodynamics of RDS, and the role that boundaries play. We also remark that there are exam-ples of more complicated multi-domain and bulk-surface models, where the role of boundary conditionsand geometry have major impacts, and there are important connections between those models and the spa-tially heterogeneous ones studied in Section 2.1 [39]. Finally we have restricted our attention to studyingstationary solutions, but RDS can exhibit a wide variety of spatiotemporal behaviours, such as spatiotem-poral chaos and spiral waves observed in excitable systems [60]. Preliminary simulations of such dynamicssuggests that our mixed boundary conditions can also lead to localized behaviour in these systems, withpossible applications in electrophysiology and other areas (e.g. little is known regarding how boundariesinfluence re-entrant wave dynamics [8]). While there has been great progress in extending Turing’s insightsinto pattern formation and RDS, there is substantial work left, especially in terms of elucidating its connec-tions with real development. References
1. D. Avitabile, V. F. Bre ˜na Medina, and M. J. Ward. Spot dynamics in a reaction-diffusion model of plantroot hair initiation.
SIAM Journal on Applied Mathematics , 78(1):291–319, 2018.2. P. Ball.
The Self-made Tapestry: Pattern Formation in Nature . Oxford University Press, 2001.3. D. L. Benson, J. A. Sherratt, and P. K. Maini. Diffusion driven instability in an inhomogeneous domain.
Bulletin of Mathematical Biology , 55(2):365–384, 1993.4. P. Borckmans, G. Dewel, A. De Wit, and D. Walgraef. Turing bifurcations and pattern selection. In
Chemical waves and patterns , pages 323–363. Springer, 1995.5. B. Bozzini, G. Gambino, D. Lacitignola, S. Lupo, M. Sammartino, and I. Sgura. Weakly nonlinear anal-ysis of Turing patterns in a morphochemical model for metal growth.
Comput. Math. Appl. , 70(8):1948–1969, 2015.6. R. A. Bradshaw and E. A. Dennis.
Handbook of Cell Signaling . Academic press, 2009.7. J. Briscoe and S. Small. Morphogen rules: design principles of gradient-mediated embryo patterning.
Development , 142(23):3996–4009, 2015.8. R. Clayton, O. Bernus, E. Cherry, H. Dierckx, F. H. Fenton, L. Mirabella, A. V. Panfilov, F. B. Sachse,G. Seemann, and H. Zhang. Models of cardiac tissue electrophysiology: progress, challenges and openquestions.
Progress in biophysics and molecular biology , 104(1-3):22–48, 2011.
9. E. J. Crampin, W. W. Hackborn, and P. K. Maini. Pattern formation in reaction-diffusion models withnonuniform domain growth.
Bulletin of Mathematical Biology , 64(4):747–769, 2002.10. M. C. Cross and P. C. Hohenberg. Pattern formation outside of equilibrium.
Reviews of modern physics ,65(3):851, 1993.11. P. De Kepper, V. Castets, E. Dulos, and J. Boissonade. Turing-type chemical patterns in the chlorite-iodide-malonic acid reaction.
Physica D: Nonlinear Phenomena , 49(1-2):161–169, 1991.12. A. De Wit, P. Borckmans, and G. Dewel. Twist grain boundaries in three-dimensional lamellar Turingstructures.
Proceedings of the National Academy of Sciences , 94(24):12765–12768, 1997.13. A. De Wit, G. Dewel, P. Borckmans, and D. Walgraef. Three-dimensional dissipative structures inreaction-diffusion systems.
Physica D: Nonlinear Phenomena , 61(1-4):289–296, 1992.14. X. Diego, L. Marcon, P. M ¨uller, and J. Sharpe. Key features of Turing systems are determined purely bynetwork topology.
Physical Review X , 8(2):021071, 2018.15. R. Dillon, P. Maini, and H. Othmer. Pattern formation in generalized Turing systems.
Journal of Mathe-matical Biology , 32(4):345–393, 1994.16. C. Duckett, C. Grierson, P. Linstead, K. Schneider, E. Lawson, C. Dean, S. Poethig, and K. Roberts. Clonalrelationships and cell patterning in the root epidermis of arabidopsis.
Development , 120(9):2465–2474,1994.17. S.-I. Ei and T. Ishimoto. Dynamics and interactions of spikes on smoothly curved boundaries forreaction–diffusion systems in 2d.
Japan Journal of Industrial and Applied Mathematics , 30(1):69–90, 2013.18. B. Ermentrout. Stripes or spots? Nonlinear effects in bifurcation of reaction-diffusion equations on thesquare.
Proc. Math. Phys. Sci. , 434(1891):413–417, 1991.19. M. Esposito. Open questions on nonequilibrium thermodynamics of chemical reaction networks. arXivpreprint arXiv:2007.00719 , 2020.20. G. Falasco, R. Rao, and M. Esposito. Information thermodynamics of Turing patterns.
Physical ReviewLetters , 121(10):108301, 2018.21. U. Fischer, Y. Ikeda, K. Ljung, O. Serralbo, M. Singh, R. Heidstra, K. Palme, B. Scheres, and M. Grebe.Vectorial information for arabidopsis planar polarity is mediated by combined aux1, ein2, and gnomactivity.
Current Biology , 16(21):2143–2149, 2006.22. R. FitzHugh. Mathematical models of threshold phenomena in the nerve membrane.
The Bulletin ofMathematical Biophysics , 17(4):257–278, 1955.23. R. FitzHugh. Impulses and physiological states in theoretical models of nerve membrane.
Biophysicaljournal , 1(6):445, 1961.24. A. Gierer and H. Meinhardt. A theory of biological pattern formation.
Kybernetik , 12(1):30–39, 1972.25. J. D. Glover, K. L. Wells, F. Matth¨aus, K. J. Painter, W. Ho, J. Riddell, J. A. Johansson, M. J. Ford, C. A.Jahoda, V. Klika, et al. Hierarchical patterning modes orchestrate hair follicle morphogenesis.
PLoSbiology , 15(7):e2002117, 2017.26. J. B. A. Green and J. Sharpe. Positional information and reaction-diffusion: two big ideas in develop-mental biology combine.
Development , 142(7):1203–1211, 2015.27. W. K. Ho, L. Freem, D. Zhao, K. J. Painter, T. E. Woolley, E. A. Gaffney, M. J. McGrew, A. Tzika, M. C.Milinkovitch, P. Schneider, et al. Feather arrays are patterned by interacting signalling and cell densitywaves.
PLoS biology , 17(2):e3000132, 2019.28. D. Iron and M. J. Ward. A metastable spike solution for a nonlocal reaction-diffusion model.
SIAMJournal on Applied Mathematics , 60(3):778–802, 2000.29. O. Jensen, V. O. Pannbacker, G. Dewel, and P. Borckmans. Subcritical transitions to Turing structures.
Physics Letters A , 179(2):91–96, 1993.30. J. A. Johansson and D. J. Headon. Regionalisation of the skin. In
Seminars in cell & developmental biology ,volume 25, pages 3–10. Elsevier, 2014.31. J. P. Keener and J. Sneyd.
Mathematical Physiology , volume 1. Springer.32. J. Kernevez, G. Joly, M. Duban, B. Bunow, and D. Thomas. Hysteresis, oscillations, and pattern forma-tion in realistic immobilized enzyme systems.
Journal of Mathematical Biology , 7(1):41–56, 1979.33. K. Kishimoto and H. F. Weinberger. The spatial homogeneity of stable equilibria of some reaction-diffusion systems on convex domains.
Journal of Differential Equations , 58(1):15–21, 1985.34. V. Klika, R. E. Baker, D. Headon, and E. A. Gaffney. The influence of receptor-mediated interactions onreaction-diffusion mechanisms of cellular self-organisation.
Bulletin of Mathematical Biology , 74(4):935– solating Patterns in Open Reaction-Diffusion Systems 23
SIAM Journal on Applied Mathematics , 78(5):2298–2322, 2018.36. S. Kondo and T. Miura. Reaction-diffusion model as a framework for understanding biological patternformation. science , 329(5999):1616–1620, 2010.37. M. Koz´ak, E. A. Gaffney, and V. Klika. Pattern formation in reaction-diffusion systems with piecewisekinetic modulation: An example study of heterogeneous kinetics.
Physical Review E , 100(4):042220, 2019.38. A. L. Krause, M. A. Ellis, and R. A. Van Gorder. Influence of curvature, growth, and anisotropy on theevolution of Turing patterns on growing manifolds.
Bulletin of Mathematical Biology , 81(3):759–799, 2019.39. A. L. Krause, V. Klika, J. Halatek, P. K. Grant, T. E. Woolley, N. Dalchau, and E. A. Gaffney. Turingpatterning in stratified domains.
Bulletin of Mathematical Biology , 2020. In press.40. A. L. Krause, V. Klika, T. E. Woolley, and E. A. Gaffney. From one pattern into another: Analysis ofTuring patterns in heterogeneous domains via WKBJ.
Journal of the Royal Society Interface , 17:20190621,2020.41. Y. Kuramoto.
Chemical Oscillations, Waves, and Turbulence . Dover books on chemistry. Dover Publica-tions, 2003.42. L. Kurowski, A. L. Krause, H. Mizuguchi, P. Grindrod, and R. A. Van Gorder. Two-species migrationand clustering in two-dimensional domains.
Bulletin of Mathematical Biology , 79(10):2302–2333, 2017.43. M. Leda, V. K. Vanag, and I. R. Epstein. Instabilities of a three-dimensional localized spot.
PhysicalReview E , 80(6):066204, 2009.44. X. Li, A. M. Udager, C. Hu, X. T. Qiao, N. Richards, and D. L. Gumucio. Dynamic patterning at thepylorus: Formation of an epithelial intestine–stomach boundary in late fetal life.
Developmental dynamics ,238(12):3205–3217, 2009.45. P. Maini and M. Myerscough. Boundary-driven instability.
Applied Mathematics Letters , 10(1):1–4, 1997.46. P. K. Maini, T. E. Woolley, R. E. Baker, E. A. Gaffney, and S. S. Lee. Turing’s model for biological patternformation and the robustness problem.
Interface focus , 2(4):487–496, 2012.47. H. Meinhardt. A boundary model for pattern formation in vertebrate limbs.
Development , 76(1):115–137,1983.48. Y. Miyamoto. Stability of a boundary spike layer for the gierer–meinhardt system.
European Journal ofApplied Mathematics , 16(4):467–491, 2005.49. J. D. Murray. A pre-pattern formation mechanism for animal coat markings.
Journal of Theoretical Biology ,88(1):161–199, 1981.50. J. D. Murray.
Mathematical Biology. II. Spatial models and biomedical applications . Interdisciplinary appliedmathematics. Springer, New York, 2004.51. S. M. Murray and V. Sourjik. Self-organization and positioning of bacterial protein clusters.
NaturePhysics , 13(10):1006–1013, 2017.52. J. Nagumo, S. Arimoto, and S. Yoshizawa. An active pulse transmission line simulating nerve axon.
Proceedings of the IRE , 50(10):2061–2070, 1962.53. G. Nicolis and I. Prigogine.
Self-Organization in Nonequilibrium Systems: From Dissipative Structures toOrder Through Fluctuations . A Wiley-Interscience publication. Wiley, 1977.54. K. Page, P. K. Maini, and N. A. M. Monk. Pattern formation in spatially heterogeneous Turing reaction–diffusion models.
Physica D: Nonlinear Phenomena , 181(1-2):80–101, 2003.55. K. M. Page, P. K. Maini, and N. A. M. Monk. Complex pattern formation in reaction–diffusion systemswith spatially varying parameters.
Physica D: Nonlinear Phenomena , 202(1-2):95–115, 2005.56. J. E. Pearson and W. J. Bruno. Pattern formation in an n+ q component reaction–diffusion system.
Chaos:An Interdisciplinary Journal of Nonlinear Science , 2(4):513–524, 1992.57. I. Prigogine and G. Nicolis. Biological order, structure and instabilities.
Quarterly Reviews of Biophysics ,4(2-3):107–148, 1971.58. J. Raspopovic, L. Marcon, L. Russo, and J. Sharpe. Digit patterning is controlled by a bmp-sox9-wntTuring network modulated by morphogen gradients.
Science , 345(6196):566–570, 2014.59. J. Ross.
Thermodynamic and Stochastic Theory of Reaction–Diffusion Systems , pages 41–58. Springer BerlinHeidelberg, Berlin, Heidelberg, 2008.60. F. S´anchez-Garduno, A. L. Krause, J. A. Castillo, and P. Padilla. Turing–hopf patterns on growing do-mains: the torus and the sphere.
Journal of theoretical biology , 481:136–150, 2019.
61. R. A. Satnoianu, M. Menzinger, and P. K. Maini. Turing instabilities in general systems.
Journal ofMathematical Biology , 41(6):493–512, 2000.62. J. Schnakenberg. Simple chemical reaction systems with limit cycle behaviour.
Journal of theoreticalbiology , 81(3):389–400, 1979.63. N. S. Scholes, D. Schnoerr, M. Isalan, and M. P. Stumpf. A comprehensive network atlas reveals thatTuring patterns are common but not robust.
Cell systems , 9(3):243–257, 2019.64. H. Serna, A. P. Mu ˜nuzuri, and D. Barrag´an. Thermodynamic and morphological characterizationof Turing patterns in non-isothermal reaction–diffusion systems.
Physical Chemistry Chemical Physics ,19(22):14401–14411, 2017.65. S. Setayeshgar and M. Cross. Turing instability in a boundary-fed system.
Physical Review E , 58(4):4485,1998.66. H. Shoji and K. Yamada. Most stable patterns among three-dimensional Turing patterns.
Japan journalof industrial and applied mathematics , 24(1):67, 2007.67. D. M. Smith, C. Nielsen, C. J. Tabin, and D. J. Roberts. Roles of bmp signaling and nkx2. 5 in patterningat the chick midgut-foregut boundary.
Development , 127(17):3671–3681, 2000.68. S. Subramanian and S. M. Murray. Pattern selection in reaction diffusion systems. arXiv preprintarXiv:2005.07940 , 2020.69. N. P. Taylor, H. Kim, A. L. Krause, and R. A. Van Gorder. A non-local cross-diffusion model of popula-tion dynamics I: Emergent spatial and spatiotemporal patterns.
Bulletin of Mathematical Biology , 82(112),2020.70. A. Tucker and P. Sharpe. The cutting-edge of mammalian development; how the embryo makes teeth.
Nature Reviews Genetics , 5(7):499–508, 2004.71. A. M. Turing. The chemical basis of morphogenesis.
Philosophical Transactions of the Royal Society ofLondon. Series B, Biological Sciences , 237(641):37–72, 1952.72. R. A. Van Gorder. Influence of temperature on Turing pattern formation.
Proc. R. Soc. A. , 476:20200356,2020.73. R. A. Van Gorder, V. Klika, and A. L. Krause. Turing conditions for pattern forming systems on evolvingmanifolds. 2020. arXiv:1904.09683 [nlin.PS].74. C. Varea, J. Arag ´on, and R. Barrio. Confined Turing patterns in growing systems.
Physical Review E ,56(1):1250, 1997.75. D. W. Walsh, C. Godson, D. P. Brazil, and F. Martin. Extracellular bmp-antagonist regulation in devel-opment and disease: tied up in knots.
Trends in cell biology , 20(5):244–256, 2010.76. K. D. Walton, M. Whidden, ˚A. Kolterud, S. K. Shoffner, M. J. Czerwinski, J. Kushwaha, N. Parmar,D. Chandhrasekhar, A. M. Freddo, S. Schnell, et al. Villification in the mouse: Bmp signals controlintestinal villus patterning.
Development , 143(3):427–436, 2016.77. M. J. Ward, D. McInerney, P. Houston, D. Gavaghan, and P. Maini. The dynamics and pinning of a spikefor a reaction-diffusion system.
SIAM Journal on Applied Mathematics , 62(4):1297–1328, 2002.78. T. E. Woolley, R. E. Baker, E. A. Gaffney, and P. K. Maini. Stochastic reaction and diffusion on growingdomains: understanding the breakdown of robust pattern formation.
Physical Review E , 84(4):046216,2011.79. T. E. Woolley, R. E. Baker, and P. K. Maini. Turing’s theory of morphogenesis: where we started, wherewe are and where we want to go. In