Isolating Patterns in Open Reaction–Diffusion Systems

Realistic examples of reaction–diffusion phenomena governing spatial and spatiotemporal pattern formation are rarely isolated systems, either chemically or thermodynamically. However, even formulations of ‘open’ reaction–diffusion systems often neglect the role of domain boundaries. Most idealizations of closed reaction–diffusion systems employ no-flux boundary conditions, and often patterns will form up to, or along, these boundaries. Motivated by boundaries of patterning fields related to the emergence of spatial form in embryonic development, we propose a set of mixed boundary conditions for a two-species reaction–diffusion system which forms inhomogeneous solutions away from the boundary of the domain for a variety of different reaction kinetics, with a prescribed uniform state near the boundary. We show that these boundary conditions can be derived from a larger heterogeneous field, indicating that these conditions can arise naturally if cell signalling or other properties of the medium vary in space. We explain the basic mechanisms behind this pattern localization and demonstrate that it can capture a large range of localized patterning in one, two, and three dimensions and that this framework can be applied to systems involving more than two species. Furthermore, the boundary conditions proposed lead to more symmetrical patterns on the interior of the domain and plausibly capture more realistic boundaries in developmental systems. 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 selection mechanisms.


Introduction
Reaction-diffusion systems (RDS) are employed to model an increasingly wide-range of phenomena (Ball 2001;Kuramoto 2003;Murray 2004). Following Turing's work (Turing 1952), an enormous literature has developed using these models to capture aspects of patterns seen in biology and chemistry (Gierer and Meinhardt 1972;Murray 1981;De Kepper et al. 1991;Cross and Hohenberg 1993;Kondo and Miura 2010;Maini et al. 2012). However, as noted by Turing, these models are idealizations and cannot account for many observations of real developmental systems (Maini et al. 2012). We further extend present efforts to increase the realism, and hence the explanatory power, of this theory by considering the role of boundary conditions in isolating patterns away from domain boundaries. By 'pattern,' we mean a localized solution (a non-random spatial arrangement) of a morphogen above a threshold value, which in embryonic settings would correspond to a sufficient level to induce biologically meaningful change in cell state. Our main aim here is to demonstrate that a simple choice of boundary conditions can guarantee interior localization of patterns. By deriving these boundary conditions from a heterogeneous problem, we argue that such pattern isolation can be understood as arising from different kinds of boundaries due to differential gene expression 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 chemical chaos, are genuinely non-equilibrium thermodynamic processes. Since the early work of Nicolis and Prigogine (Prigogine and Nicolis 1971;Nicolis and Prigogine 1977), substantial further work has extended this thermodynamic perspective of Turing instabilities as a non-equilibrium phenomenon associated with open systems (Cross and Hohenberg 1993;Ross 2008;Falasco et al. 2018;Esposito 2020). However, very little emphasis in these works is put on the role of boundary conditions. In the classical review of pattern selection (Borckmans et al. 1995), for instance, the authors explicitly state that they will not be concerned with the role of boundary conditions. While many authors describe the impact of boundary conditions on such non-equilibrium phenomena (Cross and Hohenberg 1993;Murray 2004), relatively few consider more exotic conditions beyond the standard Neumann, Dirichlet, or periodic settings (Setayeshgar and Cross 1998;Dillon et al. 1994;Klika et al. 2018). Physically, periodic boundary conditions can be justified for a flat approximation of a curved manifold, which can be appropriate for patterns appearing across the entire surface of an organism. However, the typical Neumann boundary conditions used to study pattern formation are very clearly a mathematical idealization of what are often open systems. Arcuri and Murray (1986) have suggested that such idealized boundaries are likely unrealistic and that inhomogeneous boundary conditions can lead to more robust pattern formation (i.e. less sensitivity on initial conditions).
Many developmental systems exhibit localized patterning within a structurally homogeneous field (i.e. a domain without any spatial heterogeneity in diffusive fluxes or reactions due to variations in cell type or arrangement). Examples include ectodermal structures such as hair or teeth (Tucker and Sharpe 2004;Johansson and Headon 2014) and the role of auxin in plant root initiation (Duckett et al. 1994;Fischer et al. 2006; Avitabile et al. 2018). Periodic patterning of the primary hair follicles 1 in mammalian whiskers and hair more generally, for instance, suggests a regionalization of different patterning fields within which RDS may explain the emergence of patterns. However, as noted in Murray (2004), Neumann boundary conditions can exhibit patterning all the way up to the boundary, often leading to partial patterns at the boundary (e.g. half-spots in the case of spotty patterns). Of course, inhomogeneous solutions of nonlinear RDS may exhibit a variety of different behaviours, depending on the parameters and initial data of the system. In many cases (such as in the limit of large diffusivity ratios), it has been shown that spike solutions can approach, and be pinned to, boundaries, particularly points of extremal curvature in multiple spatial dimensions (Iron and Ward 2000a, b;Kolokolnikov and Ward 2004;Miyamoto 2005;Ei and Ishimoto 2013). In other cases, initial data consisting only of internal spots can be shown to remain in a configuration of internal spots due to boundary repulsion (Kolokolnikov et al. 2009;Chen and Ward 2011). Numerically, as shown in Sect. 3, Neumann conditions will generically allow patterns to form up to the boundary for a variety of systems and parameters starting from small perturbations of a homogeneous steady state.
There have been several different approaches to designing RDS which exhibit patterns isolated away from the boundaries. A simple approach pursued by Varea et al. (1997) was to modify the reaction kinetics at a boundary to push the system outside the Turing regime locally, so that patterning was restricted to an interior region. Similar ideas were explored with a variety of more complex spatial heterogeneities in reaction kinetics (Page et al. 2003(Page et al. , 2005, or diffusive fluxes (Benson et al. 1993). Recently, this kind of localized patterning has been justified in the linear regime either in the case of smoothly varying heterogeneity  or jumps (step functions) in the kinetics (Kozák et al. 2019). Similar work has been considered in some nonlinear regimes, where spike pinning, and hence pattern localization, can be obtained via heterogeneity in kinetics (Ward et al. 2002;Avitabile et al. 2018). Another approach, pursued in the context of the positioning of bacterial protein clusters, is to understand pattern selection from a nonlinear theory and look for nonlinearities which exhibit the desired isolated patterns (Murray and Sourjik 2017;Subramanian and Murray 2021).
More generally, several studies have investigated inhomogeneous or mixed boundary conditions. Dirichlet boundary conditions can have a variety of effects on patterns, including modifying parameters for which patterns are observed (Maini and Myerscough 1997). Robin conditions in pattern-forming RDS have recently been shown to impact the number of interior spots (Tzou et al. 2011), and inhomogeneous fluxes of inhibitor can lead to movement of interior pulse solutions in 1-D towards the source of the flux (in addition to changing amplitudes and spacing) (Tzou and Ward 2018). Closer to what we will propose here, Dillon et al. (1994) investigated a variety of different boundary conditions for each species in a two-species RDS, demonstrating that different pattern selection mechanisms could be influenced by the choice of these boundary conditions.
While we are primarily interested in understanding isolated patterns, our work also falls into the larger context of trying to connect Turing's simple patterning mechanism with the complex reality of biological development. As Turing himself said, pattern formation occurs in stages through subsequent patterning (Turing 1952). There are difficult questions of how to idealize such situations to elucidate the fundamental mechanisms at play in any particular stage. One such question is the robustness of Turing patterns, both in terms of the dependence of steady patterned states on initial data, and the strict requirements on parameters needed for RDS to admit Turing-type patterns (Maini et al. 2012;Woolley et al. 2017;Scholes et al. 2019). Domain growth has been suggested as one way to help improve robustness (Crampin et al. 2002;Krause et al. 2019;Van Gorder et al. 2021), as has stochasticity (Woolley et al. 2011).
Here, we will suggest that RDS exhibiting isolated patterning will also be more robust in both enlarging the parameter space within which we see patterns, as discussed in Maini and Myerscough (1997), 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. In particular, boundaries can arise due to either explicit heterogeneity in the tissue (where cells of different types explicitly demarcate patterning regions, due to previous fate determination) or due to more complex and diffuse mechanisms involving cell state, which may be refined into sharp boundaries via, e.g. bistability. Such boundaries can be thought of as an example of positional information (Meinhardt 1983;Green and Sharpe 2015). We briefly discuss the biology of such boundary formation before presenting our modelling framework.

Biological Interlude: Boundaries in Developmental Patterning
Activator and inhibitor species are taken to represent specific molecules (or more accurately, signalling pathways) in biological systems. These undergo processes of synthesis, decay, and diffusion, depending on the 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. The synthesis 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 enzymes that catalyse their formation. These molecules are secreted from cells, then diffuse, either alone or interacting with other proteins, through subcellular structures, or the meshwork of extracellular matrix. Ultimately, these molecules are capable of attaching to receptors on the surface of, or sometimes within, cells in the vicinity. The effect of receptor binding is to trigger a typically multi-step process of signal transduction, ultimately altering rates of gene expression (molecule synthesis), or cell behaviour (such as cell shape or movement) (Bradshaw and Dennis 2009).
Many extracellular proteins that act as signals are themselves bound in the extracellular space by inhibitory or transport proteins. Thus, inactive complexes can be formed, capable of diffusing but unable to bind receptors and elicit a response. An example is the Bone Morphogenetic Protein (BMP) family, several members of which have been implicated in reaction-diffusion periodic patterning systems of skin (Ho et al. 2019;Glover et al. 2017), limb skeleton (Raspopovic et al. 2014), and gut (Walton et al. 2016), and for which a range of distinct extracellular inhibitor proteins have been characterized (Walsh et al. 2010).
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 inhibitor molecules in many systems do not interact directly, but rather through chains of intermediates in signalling pathways. The specificity of molecular interactions and indirect nature of activator-inhibitor interaction permit selective spatial variation in the properties of the activator and inhibitor that we will explore in this paper. Boundaries relevant to reaction-diffusion patterning will be defined by the geometry of the embryo or organ, and by spatial variation in the expression of genes that encode or modulate activator and inhibitor synthesis, diffusion, biological activity, and decay. Geometric boundaries are unavoidable and often prominent in experimental systems at the cut edges of the tissue, as they have a strong influence on pattern formation (Glover et al. 2017).
Variation in gene expression can come about due to differences in tissue structure that cause entirely distinct cell types to be directly apposed. This is observed in composite organs such as skin and gut, in which distinct and non-intermingling cell types with deeply divergent developmental origins are organized into a tightly-packed epithelial sheet and a looser mesenchyme rich in extracellular matrix. In each of these tissues, the cells express a distinct subset of the genome, have different shapes and mobility, and the extracellular space between the cells has very different properties. These two tissue types are separated by a specialized planar matrix called the basement membrane. The locations of boundaries between these components are marked 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 primarily of cells of the same type) can arise from previous patterning events, such as an external organizing region emitting a graded signal that influences part of the patterning field. Graded boundaries may persist in that form, or can be refined into sharp boundaries by, for example, the action of mutually antagonistic factors that create bistable systems, notable in the formation of discrete segments (Briscoe and Small 2015). If the output of such systems influences the production of a receptor or an extracellular inhibitor for a component of the reaction-diffusion system, then across the boundary this factor would be predicted to diffuse, but to lack any biological activity. Boundaries characterized by differences in gene expression may be present in structurally homogeneous tissues and require application of specialized molecular methods to reveal them. For example, in vertebrate embryonic gut regionalization, initially diffuse boundaries become sharper (Li et al. 2009), and BMP receptor expression becomes distinct in different regions or segments along the tract (Smith et al. 2000). In Fig. 1a, we depict both kinds of boundaries, and in (b) we show an example of a boundary observed in the ability of a tissue to respond to BMP within the otherwise-homogeneous gut tract described above. 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 the dumbbell-shaped extracellular signalling molecules. The red and blue rectangular cells represent two distinct cell states, with the blue capable 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 of receptor expression towards the right boundary (which itself may be a hard no-flux boundary), leading to a more diffuse intracellular concentration of the protein u. b In situ hybridization of a chick gut (developmental stage E4.5) with a riboprobe to BMPR1A from Fig. 1(H) of Smith et al. (2000). The purple/blue stain shows the ability of the tissue to receive BMP signals, whereas the white/translucent areas lack this receptor and do not have the same response to BMP (though BMP may still diffuse throughout these regions) (Color Figure Online) Thus, the spatially varying presence or absence of a receptor, a component of its signal transduction pathway, 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 a selective effect on a single species (activator or inhibitor) in the reaction-diffusion system, which will be free to diffuse across the boundary and persist physically, but without the ability to influence reaction kinetics in this region of the field.

Mixed Boundary Conditions
We now consider modelling such boundaries in a reaction-diffusion system. Throughout the paper, we will consider a generic two-species RDS in the domain x ∈ ⊂ R n : where D > 1 is a ratio of diffusion coefficients. In the classical case of activatorinhibitor kinetics, where v is a self-inhibitor, D > 1 is a necessary condition for inhomogeneous steady solutions (at least in convex domains (Kishimoto and Wein-berger 1985)). Hereafter, we restrict attention to this case in our analysis, and consider u and v as presumptive activator and inhibitor, respectively. We will show that inhomogeneous steady solutions (Turing-type patterns and otherwise) can be isolated away from the domain boundary. This can be achieved by using the boundary conditions where R should be approximately a minimum of u in a patterned state (e.g. with Neumann boundary conditions). Specifically, if we compute a stable steady-state solution given Neumann boundary conditions and initial data perturbed from the homogeneous equilibrium, we will take R = min (u). Below we will show that the precise value of R seems to play only a small role in influencing the pattern, at least for some reaction kinetics. While it can have quantitative effects on the size and spacing of localized solutions, our numerical results suggest that taking R less than the homogeneous steady state suffices to lead to interior patterning. In some sense this is obvious, as we are fixing the value of u at the boundary. We will focus on the case where the activator is fixed at the boundary and the inhibitor is not allowed to diffuse (i.e. satisfies zero flux conditions). In principle, our results on interior localization will hold if we swap the boundary conditions between activator and inhibitor, but the lengthscale of the boundary influence will be larger in the case of the inhibitor being fixed (as D > 1, and often D 1, for a pattern-forming RDS). Finally, we note that other values of R may be suitable for the development of interior localization, particularly in the presence of multistability or multiple species interactions; see Sect. 3.3 for examples.
In developmental settings, the discussion in the preceding subsection allows us to consider u and v as behaving differently outside the domain due to boundaries in receptor distribution or gene expression for example; in turn, this can lead to different boundary conditions for the activator and inhibitor. We will justify the boundary conditions given in (3) asymptotically in Sect. 2.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 a chemical which undergoes exothermic or endothermic reactions (Serna et al. 2017;Van Gorder 2020). Throughout, we will refer to (3) as mixed boundary conditions, though other authors have used the term nonscalar boundary conditions (Dillon et al. 1994). We will contrast these boundary conditions with homogeneous Neumann conditions on both species, and refer to the latter simply as Neumann conditions for brevity throughout the manuscript.
The rest of the paper is organized as follows. In Sect. 2, we analytically show how the mixed boundary conditions (3) can be obtained from a heterogeneous problem on an enlarged domain. We also describe how these boundary conditions force pattern isolation away from the boundary in terms of a local pattern-selection mechanism. In Sect. 3, we demonstrate a broad numerical exploration of this interior localization across a variety of two-species reaction kinetics with qualitatively distinct kinds of patterns. We also show that the same mechanism can operate in three spatial dimensions, as well as in systems with more than two species. Finally, in Sect. 4 we conclude by discussing both phenomenological and mechanistic applications of these boundary conditions, as well as a number of further research directions.

Heterogeneity and Isolated Patterning
We now explore two aspects of the boundary conditions (3). First we explain how they arise from an asymptotic analysis of a heterogeneous problem, motivated by regions of different gene expression as discussed in Sect. 1.1. Viewing these conditions as arising from spatial heterogeneity is consistent both with the biological observations discussed there, as well as with the theoretical literature on localization of patterns in heterogeneous media. We then discuss steady-state selection mechanisms observed from different boundary conditions, explaining some aspects of the isolated patterning we expect to achieve with these boundary conditions.

Derivation of Mixed Boundary Conditions from a Heterogeneous Problem
Working non-dimensionally throughout, we consider a larger domain˜ ⊂ R n , and a strict subset of it ⊂˜ such that˜ \ is connected, with all domain boundaries sufficiently smooth. We now consider the non-dimensional heterogeneous reactiondiffusion system, 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 transduction and gene expression) in the medium. Within the interior region, , the normal activator-inhibitor dynamics are present, whereas in the exterior region,˜ \ , the activator relaxes to an equilibrium concentration on a nondimensional timescale of 1/ρ, while the inhibitor may interact with the activator in a different manner g 2 with a timescale T g 2 . Such dynamics are easily obtained in examples where gene expression is spatially modulated due to previous symmetry breaking and fate specification events, leading to heterogeneity in local signalling dynamics or in tissue/cell type as discussed in Sect. 1.1.
One can show that for sufficiently large ρ, the one-dimensional steady-state version of Eqs. (4)-(5) can, to leading asymptotic order, be reduced to considering the steadystate problem (1)-(2) on the smaller domain with the boundary conditions (3). In higher dimensions, however, one must impose additional geometric constraints. For concreteness, to capture the general case we consider the problem in two spatial dimensions, remarking that extensions to higher dimensions follow the same reasoning. We let = ∂ be the boundary between the internal and external regions of˜ (equiva- Fig. 2 Geometry of the full domain˜ , with a region in˜ \ enlarged to show various asymptotic parameters and coordinates along the boundary discussed in the text. Note that η is, in general, a ratio of the maximum thickness of˜ \ to arclength of , but for illustrative purposes we assume the arclength is order unity. Similarly, we assume sufficient smoothness of both internal and exterior boundaries so that the local mean curvature is always order unity or smaller lently, the internal boundary of˜ \ ) and define η as the (non-dimensional) ratio of the maximum thickness (normal to ) of˜ \ to the arclength of . The one-dimensional case follows exactly the same reasoning as given below for two dimensions, without the technicalities related to the geometry.
We proceed to show that for sufficiently large ρ and for˜ \ sufficiently slender and regular, steady states of this system asymptotically reduce to those of the reactiondiffusion system (1)-(3) defined on the interior domain . Additionally, the dynamics of u and v outside the interior domain are simple and (asymptotically) determined analytically. In particular, with η 1 denoting the assumption of slenderness for \ , we consider the asymptotic regime ρ 1 and η 1. We focus on steady states, as transient solutions to the heterogeneous problem (4)-(5) need not satisfy the no-flux boundary condition for v on ∂ in this asymptotic scaling without further restrictions on kinetic timescales. Our results will then show that the steady states of (1)-(3) will be precisely the same as those of (4)-(5). Similarly, we will need to exploit the geometry of a slender domain to prevent problems with transverse gradients in v; in one spatial dimension the restriction of η 1 can be relaxed completely. See Fig. 2 for a geometric depiction.
Our first objective is to determine how the solution within the interior of the boundary of is related to the solution in the exterior of on approaching the boundary . To proceed, we first mollify the transition between the two regions with a function H δ which is centred on this boundary and transitions from zero to unity in the normal direction of the inner boundary of˜ \ . We also assume has a well-defined and bounded curvature. Furthermore, the transition across this boundary is on a (dimensionless) lengthscale of δ η, and since the mollification scale δ is for mathematical convenience, it is taken to be smaller than any physical scale below.
The mollified steady-state equations on the extended domain are We introduce an orthonormal curvilinear coordinate system given by (ζ 1 , ζ 2 ) along the curve , where the level set ζ 1 = 0 corresponds to with the ζ 1 axis perpendicular to and pointing into the interior of˜ \ . We now follow the derivation of an Eikonal equation along the boundary [see, e.g. Keener and Sneyd (1998)]. Using the Einstein summation convention, the chain rule gives and hence, for example, We have that the vector is normal to and, by the need for orthonormality, we fix the scale of ζ 1 so that α is a unit normal and then κ := ∇ · α is the mean curvature of . In particular, assuming the curvature is of order unity relative to the small parameter δ entails the scale of ζ 1 is order unity. Further noting by orthogonality of the (ζ 1 , ζ 2 ) coordinate system, we have where the coefficients c 1 , c 2 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 ζ 1 below. Derivatives of u will have an identical expansion. Given the scale of the transition of H δ , we now consider a prospective boundary layer of width of the scale δ at the boundary and thus we consider Hence, at leading order we have u ζ 1 ζ 1 = v ζ 1 ζ 1 = 0 in the transition region and so v is and similarly for u in this region. However, boundedness of this inner solution when leaving the transition region (i.e. as ζ 1 → ∞) forces P(ζ 2 ) to be zero so that matching to an outer solution is possible, and thus, The independence of u and v from the normal coordinate, ζ 1 , in the prospective boundary layer, means that on matching into this layer on either side of , we attain continuity of u and v. Within the δ-boundary layer, we have ∂v/∂ζ 1 ∼ 0 at leading order; however, this does not mean physical fluxes are essentially zero, since the mollification lengthscale δ is taken to be smaller than any physical lengthscale and v may vary significantly over the latter, with analogous remarks for u. Nonetheless, the collapse of the prospective δ-boundary layer strongly suggests the continuity of u, v across should be supplemented by continuity of flux.
To explicitly deduce this, let z 2 such that |ζ 2 − z 2 | ≤ δ 1/2 be fixed and consider the integral form of the steady-state conservation relation (given by (6)) on the pillbox region Integrating the conservation equation (6) for v over the pillbox and applying Green's theorem give and similarly for u. In particular, the absence of a boundary layer at entails that u, v are bounded, with bounded normal derivatives with bounds independent of δ. Thus, the O(δ) terms emerge from the other flux boundary integrals and the O(ρδ 3/2 ) term emerges from the integrals of sinks and sources over the pillbox, noting that the area of the pillbox P is 4δ 3/2 and the largest source term scales with ρ 1. However, the asymptotic parameter δ is from the mollification of a Heaviside function, and thus, we can take δ 3/2 ρ ∼ 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 δ 1/2 and taking the limit δ → 0, we have continuity of flux, and thus normal derivative, across for v, and analogously for u, as expected.
Knowing how to match across , we now need to consider our next objective, which is to determine the behaviour of the solutions within the region˜ \ to assign boundary conditions for a reduced model on the boundary . We take advantage of the slenderness of˜ \ , characterized by an aspect ratio of its normal extent to the arclength of via η 1. It is most useful to see if non-trivial behaviour may occur in the normal direction as, if present, such behaviour would preclude a homogeneous Neumann condition for v in the reduced system. We thus rescale ζ * 1 = ζ 1 /η, ζ * 2 = ζ 2 , and at leading order in η we have If the characteristic time of the inhibitor kinetics in the exterior, T g 2 , is large when compared to that of diffusion, η 2 D −1 T g 2 , we finally have, from (6), at least once away from by a scale of more than δ η. We have also assumed that the domain of˜ \ is sufficiently regular to ensure κ ∼ o(1/η), which will occur if the radius of curvature is on the lengthscale of 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 boundary can be specified by ζ * by construction and is also a single-valued function that satisfies η|F (ζ * 2 )| 1, where prime denotes derivative. Hence, we do not consider an external boundary consisting of a curve that, for example, has a cusp or varies too rapidly. From this assumption, we have that the normal derivative operator on the exterior boundary of˜ \ is proportional to Thus, working at leading order in η|F (ζ 2 )| 1, homogeneous Neumann conditions on the exterior boundary of˜ \ give v ζ * 1 = 0 at the external boundary, and hence, throughout˜ \ , except possibly δ-close to . As will emerge below, we ultimately consider the leading order behaviour when the asymptotic parameters satisfy 1/δ ρ 1/2 1/η 1. The restriction 1/δ ρ 1/2 entails the scale of the mollification region, δ, is the smallest scale present as required since the mollification is introduced for analytical expedience, rather than as a fundamental feature of the biophysics. The constraint ρ 1/2 1/η ensures the source strength is sufficiently large to enforce approximately homogeneous solutions in the region˜ \ , as observed below. Finally, η 1 is required to render the problem to be effectively one-dimensional at leading order and enforces˜ \ to be slender. To proceed, we consider u in the region˜ \ , whereby at leading order once away from by a scale of more than δ. Noting the Neumann boundary condition on the exterior boundary, Eq. (8) gives 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 ζ * 1 → 0, this solution is continuous with the solution on as the latter approaches for the same value of ζ * 2 , and analogously for the flux. No useful information is gained by continuity of concentration. However, continuity of flux entails ∂u/∂ζ 1 ∼ O(1) as ρ increases on approaching , assuming ∂u/∂ζ 1 in does not blow up with increasing ρ, which is consistent with numerical simulations (below) and parabolic regularity. We then have that Eq. (9) yields and hence, We further have that ρ is sufficiently large to ensure ρ 1/2 1/η, and we have already taken η 1. Noting that cosh( p) is monotonically increasing for p > 0, and Thus, providing R 1/ρ 1/2 , we have u = R is an estimate of the solution throughout \ with asymptotically small relative error. We typically take R ∼ O(1) so that a small relative error in the approximation of u = R in the region˜ \ is assured by our assumption of ρ 1/2 1/η 1. In addition, if R ∼ O(1/ρ 1/2 ) we have u ∼ O(1/ρ 1/2 ), so that the estimate u = R presents with asymptotically small absolute error, even though the 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), on the domain , which has external boundary with ρ 1/2 1/η 1 and sufficiently regular˜ \ . Continuity of flux immediately gives homogeneous Neumann conditions for v on for the reduced system, while continuity is inappropriate for v, since v is not determined in˜ \ , except ultimately via its coupling to the solution on the interior of . In contrast, we have u = R to leading order in˜ \ and thus continuity gives the Dirichlet condition u = R as a boundary condition on for the reduced system. Note the continuity of flux across is still required if considering the full system on˜ , and this would be enforced by the appropriate choice of the integration degree of freedom A(ζ * 2 ) in Eq. (9) via use of the expression for ∂u/∂ζ 1 at ζ 1 = 0 + in Eq. (10). In particular using a Neumann boundary condition for u with the reduced system, Eq. (12) would admit solutions that do not match continuity of u across the boundary. We will explore this idea further in the next subsection, referring specifically to Fig. 4 as a demonstration of this solution selection mechanism under different boundary conditions.
In summary, we have that, providing˜ \ is sufficiently slender and regular, and with ρ sufficiently large, the behaviour of steady solutions of Eqs. (4)-(5) on˜ is given, to asymptotic accuracy, by solutions of Eqs. (12) within the domain with the boundary conditions on that u = R and that the normal derivative of v, that is ∂v/∂n, is zero. We note that the set of steady-state solutions will match between the heterogeneous and reduced problems, but solutions to the time dependent problems could allow for 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 in Sect. 3. Finally, we also remark that the asymptotic restrictions, particularly regarding the slender geometry, can also be observed numerically to have typically little effect on stationary solution behaviours when comparing these two models, which we demonstrate from full numerical simulations in Sect. 3 for specific nonlinear kinetics. The driving force in deriving the conditions (3) is the heterogeneous switching between interacting and weakly or noninteracting species, which can be understood in the context of the biology described in Sect. 1.1.

Interior Patterning as Steady-State Selection
Here we further elucidate why these boundary conditions will lead to confinement of patterned states on the interior of the domain, using a particular example of reaction kinetics in a one-dimensional domain. We consider patterned steady states as subsets of all admissible ones in the periodic case and show how specific subsets of this set are selected by different boundary conditions. The motivation here is the discussion of equivariant bifurcation theory and the impact of symmetry on solution branches in Dillon et al. (1994). Specifically, we numerically 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 specific minimum value of u from the Neumann solutions. In particular, solutions with these mixed boundary conditions are forced to pattern away from the boundaries as they are exactly chosen to match a Neumann solution with patterning only on the interior. We numerically observe that small changes in R from this value 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 different 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 MAT-LAB function 'fsolve' to compute steady-state solutions to the discretized system on n = 2000 points. Writing (u * , v * ) as the homogeneous steady state, we consider an initial condition of the form (u * (1 + 10 −2 F i ), v * (1 + 10 −2 G i )), where F i , G i are independent normally distributed random variables with unit variance at each node i. Starting with periodic boundary conditions, we evolve this initial condition for t = 10 4  Table 1 under four choices of boundary conditions. Ten simulations are shown where in a, the initial simulation uses random initial data, but each subsequent simulation is initialized by cyclically shifting this solution 5% along the length L, and a stable steady state is found (by cyclic symmetry, these are automatically stable 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 = 100 units in time, and then used to find a nearby stable steady state. Note that the Neumann solutions are a subset of the periodic ones. Finally the first set of mixed steady states c are, up to numerical accuracy, a subset of the Neumann steady states b, whereas the second set a are close to, but not the same as, these. Note the colours used in c, d correspond to the last plotted solution, as all ten solution curves fall onto the same points. Parameters used were a = 0.1, b = 1.7, c = 1, D = 20, and domain length L = 25 (Color Figure Online) units of time, and then use its value to find a nearby stable steady state (which is always almost identical to the last transient state, suggesting it is reachable from random initial data). From this steady-state solution, we generate nine others, shown in Fig. 3a, by shifting this steady-state pattern cyclically. By the translational invariance of the Laplacian with periodic boundary conditions, any stable shifted steady-state solution with any shift is 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 three other sets of boundary conditions and then again find stable steady states after evolution in time. We use Neumann conditions in Fig. 3b, 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 = 0 in (d). We see precisely the pattern selection as described,  Table 1 under four choices of boundary conditions, exactly as in Fig. 3 except taking a domain length of L = 20. Note that in b, the solutions in red and blue have different amplitudes than those in light yellow and green (Color Figure  Online) with the Neumann solutions corresponding to a shifted version of those in the periodic case (hence only admitting two shifted solutions, rather than infinitely many). The mixed conditions with R set as the minimum value of u with Neumann conditions select precisely one solution from the two Neumann solutions (and no others). Finally, using R = 0 in (d) locally changes the solution from (c) near the boundary, but qualitatively does not change the solution (nor are any new solutions found in this case). Simulations with random initial conditions, as well as repeating this procedure with 100 different periodic solutions as initial data, yield precisely the same curves shown in (b)-(d) as these are the only attracting steady states for all of these different initial data. More generally, solutions satisfying Neumann boundary conditions need not give rise to periodic patterns, particularly if the domain permits half-integer modes (Murray 2004). In Fig. 4, we give an example of this following the same procedure as in Fig. 3 with a reduced domain length. Two of the four solutions given in Fig. 4b do not correspond to solutions with periodic boundary conditions shown in (a), possessing a different amplitude. Setting R as the minimal value of u from these simulations with Neumann conditions (considering only the solutions which are also periodic, shown in green in (b)), 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 influence on the solution near the boundary, as can be seen in (d). The mixed boundary conditions further exclude half-integer mode solutions.
The reduction in the admissible set of patterns has only been shown numerically for this particular example, and one would need to employ more sophisticated mathematical techniques to demonstrate this for all kinetics and parameters, and beyond one-dimensional examples. Nevertheless, the PDEs in question are local, and so any solution where v is extremal and u is minimal at the same points which satisfy both Neumann 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 that reaction-diffusion patterns will have both species sharing extremal points, but this can be shown near a Turing bifurcation via linearization, and is observed (at least approximately) in many systems numerically beyond the bifurcation (Dillon et al. 1994). We note that in Figs. 3d and 4d, the value of R chosen selects a steady solution with a nonzero boundary flux, suggesting that u is being depleted at the boundaries due to reactions outside of the domain. Hence, this is why we view these as a type of 'open' system, which are common in biological systems though less well-studied in general compared to closed systems. The pattern selection mechanism described here also suggests that peaks of the activator u should be approximately half of a wavelength away from the boundary, and we will directly test this prediction across a range of parameters and kinetics in Sect. 3.1.

Demonstrations of Robust Isolated Patterning
We now give example simulations of RDS with the boundary conditions (3), as well as those with Neumann boundary conditions. Additionally, we compare these with simulations of the heterogeneous problem (4)-(5) with Neumann conditions on the outer boundary, to demonstrate that our asymptotic reduction to the mixed boundary conditions is valid for large values of the relaxation rate ρ and suitably thin extended domains. Finally we will also give example simulations in more complicated geometries and higher-dimensional domains, as well as systems with more than two species, showing that the isolated patterning emergent from our mixed boundary conditions extends beyond these initial examples.
For all of these simulations, we used the commercial finite-element software COM-SOL. In the two-dimensional simulations, we used a minimum of 5 × 10 4 triangular second-order finite elements, and in the three-dimensional simulations we used at least 10 5 tetrahedral elements. We always chose parameters such that there was a unique homogeneous steady state for the kinetics, given by f and g for homogeneous Neumann boundary conditions, 2 and we perturbed this state by multiplying it by a normally distributed random variable of standard deviation 10 −2 and unit mean, sampled identically and independently at each finite element (as in Sect. 2.2 in the onedimensional case). All simulations shown throughout this section are given at t = 10 5 units of time, by which point they were within numerical tolerances of a stationary solution. The timestepping was implemented using a standard adaptive backwardsdifferentiation formula of orders 1 through 5. Refinements in space were used to check convergence for particular simulations, and finite-difference simulations were carried out in MATLAB (using the five-point Laplacian stencil) to check convergence to the same steady-state pattern from identical initial data.

Interior Patterning Across Two-Species Kinetics
We now describe results from simulations of each of the reaction kinetics and parameter regimes shown in Table 1. These are given in Table 2, where the first column indicates the parameter set, the second column simulations with Neumann boundary conditions, the third column simulations with conditions (3), and finally the fourth column simulations of the enlarged heterogeneous system (4)-(5). For each of these simulations, we used the same realization of the random initial condition, as all of these systems can admit multistability of different inhomogeneous steady states depending on the initial data (Borckmans et al. 1995;Jensen et al. 1993), and results are displayed 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, we can clearly see that the mixed boundary conditions lead to the inhomogeneous patterns being confined to the interior of the domain. Additionally, it is visually apparent in each case that the distance from the peak of an outermost structure (spot or stripe) to the boundary is approximately half of the wavelength between interior structures, despite the fact that R was not set to be exactly the minimum pattern value as in Sect. 2.2. Finally, as the same initial perturbations were used across all simulations, comparing the mixed boundary conditions to the heterogeneous ones, we see essentially identical steady-state solutions in all cases except for the FitzHugh-Nagumo kinetics, which admit small differences between the steady states selected. We also considered how the asymptotic reduction holds for different domain geometries and for smaller values of ρ, but for brevity omit this analysis (largely as it depends heavily on the kinetics and type of pattern studied). Essentially these results confirm our analysis in Sect. 2.1 that the mixed conditions approximate such heterogeneous problems, at least once any transients have relaxed. We note that these boundaries have discontinuous derivatives at the corners, but simulations on circular domains have identical properties. Similarly, simulations with much larger external regions (e.g. when˜ is [−L, 2L] × [−L, 2L], so that the domain length is three times larger) give steady states with interior patterning, as in those on the 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 further here.
These boundary conditions have some influence on the structure of the resulting patterns beyond the boundaries, especially in the case of the labyrinthine patterns shown in Table 2. The stripes in the Schnakenberg example show a clear confine-   (Schnakenberg 1979); GM: Gierer-Meinhardt (Gierer and Meinhardt 1972); TH: Thomas (Kernevez et al. 1979); FHN: FitzHugh-Nagumo (FitzHugh 1955(FitzHugh , 1961Nagumo et al. 1962). Note that the R values chosen here are roughly close to minimal values of typical patterns, but not exactly these values. Note that the FitzHugh-Nagumo kinetics model voltage potential relative to a ground state via the variable u and so need not admit positive solutions of this variable Table 2 Concentrations of u from simulations of the kinetics and parameters given in Table 1 with Neumann conditions, the mixed boundary conditions (3), and in a heterogeneous system given by (4) ment due to the square geometry. In the labyrinthine case for Thomas kinetics, we see that the (roughly) hexagonal arrangement of inverted spots in the Neumann case is forced into a square arrangement in the centre of the domain with the mixed conditions, with a stripe-like structure now surrounding this region. The FitzHugh-Nagumo structures are qualitatively similar on the interior of the domain, with random-looking labyrinthine stripes, but the boundary is clearly marked with a striped region having some local effects on the meandering of the stripes inside. As these labyrinthine solutions are more global in structure than localized spots, we expect less isotropic patterns and more of an influence of the confining geometry and initial conditions. This is true even in the Neumann case, as we can see by the orientation of the Schnakenberg stripes aligning in a large part of the domain with one of the axes. We also remark that the alignment of labyrinthine patterns with Neumann boundary conditions can be sensitive to initial conditions, and one can obtain quantitatively large differences for different realizations of random initial conditions. In contrast, the Schnakenberg and Thomas kinetics with our mixed boundary conditions gave consistently the same final steady state seemingly independent of a variety of initial conditions (though the positioning of the interior structures of the FitzHugh-Nagumo steady states did depend on initial conditions). Depending on the application, robust alignment may or may not be desirable. We will further explore interior patterning and pattern selection effects in the next section by considering more complicated domains. We remark that the solutions obtained are reasonably robust to choices of R, which were taken to be close to the minimal value of u in the Neumann solutions. Specifically, variations of 0.1 in any direction had no effect except in locally changing the value of u near the boundary, and larger variations in R led to local changes with some impact on pattern selection, but no real qualitative influence on the structure of solutions away from the boundary. The lengthscales L were chosen to be sufficiently large to avoid finite-size effects of a few localized structures. Finally, we remark that the lengthscale of the outer domain had no influence on the resulting patterns when it was increased to a width of 0.5L, suggesting a robustness to the qualitative effect extending beyond the thin-domain assumption made in Sect. 2.1. As this assumption is not needed in the 1-D case, we suspect it is a technicality not affecting the generic situation and could in principle be relaxed, though we do not pursue this here.

Higher-Dimensions and Complex Geometries
We now demonstrate the effects of more complicated domain geometries on patterns with the mixed boundary conditions (3). As seen in the simulations in the previous section, both spot and labyrinthine patterns were more ordered with these conditions, compared to the Neumann case. In particular, the labyrinthine 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 more isotropic character. We now show how these labyrinthine patterns are influenced by more complicated domain boundaries with Neumann and our mixed conditions. We then demonstrate that interior confinement, 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 on parameter regimes which give rise to labyrinthine solutions, comparing qualitative features such as stripe orientation and defects. We remark that for periodic boundary conditions and small-amplitude patterns near the onset of a Turing instability, there is a dichotomy between stripe and spot solutions (Ermentrout 1991), determined by signs and relative magnitudes of quadratic and cubic interactions. More recent work has shown that such a simple classification does not hold beyond the weakly nonlinear regime, and qualitatively different kinds of patterns can be simultaneously stable or interact (Borckmans et al. 1995;Jensen et al. 1993;Bozzini et al. 2015). We observe here that the largest qualitative differences are determined by the geometry (and the initial data in the Neumann setting) in the labyrinthine cases. We focus on defects (regions where the stripes are no longer contiguous at the highest value of u, which appear like dappled spots), in addition to stripe orientation.
First we consider circular and elliptical domains in Fig. 5. As anticipated, solutions with Neumann boundary conditions in Fig. 5a-c exhibit orientations which are seemingly independent of the domain geometry, though there is some correlation of stripe orientation. Stripe defects appear in all cases with Neumann 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 boundary conditions in Fig. 5d-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 stripe curvature 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, depending on what orientation most stripes take. In contrast, if one considers a contiguous region of higher-thanaverage values of u as a stripe, then the confined patterns generated by mixed boundary conditions all have exactly four of these. This is despite the fact that the largest ellipse occupies three times the area of the circle, and these observations appear to be robust to different random perturbations of the initial conditions. Finally, these results are relatively insensitive to the value of R, as simulations with R = 0.2 and R = 0.7 (not shown) gave rise to qualitatively identical patterns, with the only observable difference being the minimal value taken. These values were chosen as the minimum of u in the Neumann cases was u ≈ 0.62.
Next we consider more complicated domains by looking at those formed by perturbing the circle with periodic polar functions. We consider the following parametric domain boundary: x(s) = L cos(s) 1 + γ sin(6s), y(s) = L sin(s) 1 + γ sin(6s), for s ∈ [0, 2π ], (13) where L is now a measure of the domain length (the radius of the circle for γ = 0), and γ ∈ [0, 1) signifies the deviation from this circle. The sin(6s) term makes the perturbations resemble a 6-lobed polar rose or rhodonea curve.
In Fig. 6, we demonstrate patterns on such domains for three different values of γ . As in the elliptical case, the solutions with only Neumann conditions in (a)-(c) have stripes  Table 1 with a = 0.1, b = 1.8, c = 1, D = 10 3 , and, in the mixed case, R = 0. The maximum value of u (light yellow) is 23.56, and the minimum (dark blue) is 0 (Color Figure  Online) with similar alignments, and defects when they meet stripes with different orientations. Very little direct effect of the boundaries can be 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 Fig. 6 clearly conform to the boundary, with defects again appearing at regions where the stripes (but not necessarily the boundaries) might be expected to have large curvature. As before, while the structure of the Neumann solutions can vary tremendously depending on the initial data, the mixed boundary conditions lead to consistent patterns even when the initial data are different. In the centre of each of these domains, independent of the value of γ , is a pattern of seven inverted spots (one surrounded by six others). These simulations demonstrate an intriguing possibility of selecting for hybrid stripe and spot solutions in the domain through the discrete symmetry of the geometry. As in the elliptical case, qualitatively similar patterns are seen for R = 0.2 and R = 0.7. We now consider the same comparisons for three-dimensional patterns. Classification of patterns in higher dimensions is even more difficult than in the two-dimensional case, though there has been work on specific model systems (De Wit et al. 1992Shoji and Yamada 2007;Leda et al. 2009). As the number of Turing-type structures that can exist in three dimensions is vast, we only give two examples to demonstrate the influence of these boundary conditions. To visualize these patterns, we will only plot the solution on the tetrahedral finite elements above some threshold, so that the observed structures correspond to regions of high concentration of u.
In Fig. 7, we give examples of localized sphere-like structures which emerge again in the Schnakenberg kinetics in a rectangular domain. As expected from the twodimensional examples, spheres in the Neumann case form hemispheres and quadrants (quarter-spheres), along the boundary. In contrast, the mixed 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-dimensional case. As in the case of spots in the two-dimensional examples given in Table 2, it is clear that the patterns with these boundary conditions are substantially more ordered  Table 1 with a = 0.1, b = 1.7, c = 1, D = 30, and, in the mixed case, R = 0. The maximum value of u (light yellow) is 3.58, and the minimum (dark blue) is 0 (Color Figure Online) and symmetrical. Additionally, we find a robustness of pattern structure across different random initial conditions which contrasts with the case of Neumann boundary conditions.
In Fig. 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 included to help picture them beyond the regions of high activator concentration. The Neumann case admits complicated tubelike structures inside the domain, as well as forming partial tubes along the boundary. In contrast, the confined patterns resemble a face-centred cubic lattice with tubular connections between sphere-like regions of high activator concentration. This example 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 removed from the left boundary. In a, c, a rectangular region of height 0.7L and width 0.05L has been removed from the domain, and in b, d an elliptical region of semi-major axis of 0.35L and semi-minor axis 0.2L 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, c, and the elliptical cut boundaries in b, d. Simulations shown are on a domain of size L = 60 using the Thomas kinetics and parameters from Table 1 as labelled, with maximal and minimal values of u as in Table 2 clearly shows a difference between the seemingly more random appearance of patterns in the Neumann case, compared to the extremely structured confined patterns arising from the mixed boundary conditions. Lastly in this section, we consider domains which are composed of both Neumann and the mixed boundary conditions (3). Such a setting may arise from cutting a piece of tissue from a patterning field and then inserting an impermeable material, with local no-flux boundaries along the cut. In Fig. 9, we give examples of spots and stripes in the Thomas kinetics from Table 2, where we have shown rectangular and elliptical 'cuts' in the domain along the left boundary, so that these boundaries have Neumann conditions but all others have the mixed conditions. We see half-spots forming only along these Neumann boundaries in panels (a) and (b). In panels (c) and (d) we observe a more pronounced selection effect of the Neumann boundary conditions, where the domain with the small rectangle removed has an internal structure consistent with the corresponding simulation from Table 2 (though with patterns forming up to the Neumann boundary), but the elliptical cut induces a selection of horizontal stripe patterns even away from this boundary. Such a selection mechanism could be a useful approach to validating these boundary conditions in experimental systems.

Extensions to Many-Species
Interior 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 when restricting to patterns arising from Turing instabilities (Pearson and Bruno 1992;Satnoianu et al. 2000), and such multispecies kinetics have been studied extensively in recent years in developmental contexts (Klika et al. 2012;Diego et al. 2018;Scholes et al. 2019). 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 (Kishimoto and Weinberger 1985;Kurowski et al. 2017), though three-species models do (Taylor et al. 2020). We consider the system where D i > 0 are the diffusion coefficients of each species and a i j ∈ R are the interaction coefficients. This is a particular non-dimensionalization of a standard generalized Lotka-Volterra model with diffusion modelling 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 the prey; see Taylor et al. (2020) for more details about the ecological interpretation. We will use the following extension of the two-species mixed boundary conditions (3), u 1 (x, t) = R, n · ∇u 2 = n · ∇u 3 = 0, for all x ∈ ∂ .
We will also consider homogeneous Neumann conditions (on all species) for comparison, again simply referring to these as Neumann conditions. We show simulations of this system with parameters giving spots and labyrinthine patterns in Fig. 10 for both Neumann boundary conditions, and the mixed boundary conditions (15). As in the two species simulations above, we see the mixed boundary conditions forcing the patterned regions of u 1 into the interior of the domain, and some pattern modulation due to nonlinear interactions between patterned regions. For instance, in (b) we see a few smaller spots due to crowding, and in (d) we see a clear mode selection from the boundaries as shown before in the labyrinthine examples. Of course both of these kinds of patterns can also occur with Neumann conditions, but may be more prevalent for mixed boundary conditions as they confine patterns to the interior. Besides these influences, the patterns are qualitatively similar on the interior of the domain, especially between (a) and (b). Simulations in an extended domain, analogous to the system (4)-(5), gave the same results as using (15) (not shown).
We have chosen the species with a fixed Dirichlet boundary condition here to be u 1 , which has the smallest diffusion coefficient. We now explore applying the Dirichlet condition to the other two species. As the system is multistable in the absence of diffusion, we can exploit this to choose a value of R = 1/a ii , where i = 2 or 3 indexes Fig. 10 Values of u 1 from simulations of (14) with either Neumann or the mixed boundary conditions (15). We used a 11 = 1.1, a 12 = a 13 = a 21 = 0.8, a 23 = 0.8, a 31 = 1.7, a 32 = −1, a 33 = 1.1, D 1 = 0.1, D 2 = 100, and D 3 = 1, with a 22 = 1 in a, b and a 22 = 1.3 in c, d. was taken to be a square domain of side length L = 300. u 2 is in phase with u 1 , and u 3 is out of phase with these. The colours correspond to a maximal value of 0.7 (light yellow) and a minimal value of 0 (dark blue). Simulations are shown at time t = 10 5 , and were initialized using normal perturbations of the homogeneous steady state as described in the beginning of the section. See Figure 1 of Taylor et al. (2020) for an ecological description and interpretation of the parameters. We took R = 0 in the case of mixed boundary conditions. (Color Figure Online).

Fig. 11
Values of u 1 from simulations of (14) with mixed boundary conditions u 2 = 1/a 22 , n · ∇u 1 = n · ∇u 3 = 0 in a, b, and u 3 = 1/a 33 , n · ∇u 1 = n · ∇u 2 = 0 in c, d. We used a 11 = 1.1, a 12 = a 13 = a 21 = 0.8, a 23 = 0.8, a 31 = 1.7, a 32 = −1, a 33 = 1.1, D 1 = 0.1, D 2 = 100, and D 3 = 1, with a 22 = 1 in a, c and a 22 = 1.3 in b, d. was taken to be a square domain of side length L = 300. u 2 is in phase with u 1 , and u 3 is out of phase with these. The colours correspond to a maximal value of 0.7 (light yellow) and a minimal value of 0 (dark blue). Simulations are shown at time t = 10 5 , and were initialize using normal perturbations of the homogeneous steady state as described in the beginning of the section. (Color Figure  Online). the species with the Dirichlet condition, which forces u 1 ≈ 0 in the neighbourhood of the boundary (as both (0, 1/a 22 , 0) and (0, 0, 1/a 33 ) are steady states of the kinetic part of (14)). We show these simulations in Fig. 11. In the case of u 3 fixed at the boundary, with the intermediate diffusion value, we see qualitatively similar solutions between Figs. 10b and 11c, as well as 10d and 11d, with only small defects or pattern selection differences between them. In the case of u 2 being fixed at the boundary, however, we see much larger regions of homogeneity extending from the boundary in Fig. 11a, b, as the fixed value of u 2 can diffuse much more readily throughout the domain due to the size of its diffusive flux. This is as anticipated in the discussion at the end of Sect. 1.2. We remark that taking R = 1 has no influence on the qualitative solutions observed (not shown).
While this is only a single example of a multispecies model, we anticipate that the generalization to n species will have qualitatively similar effects. Without loss of generality, the arguments given in Sect. 2.1 can be extended to the n species case as long as all constraints are satisfied. In particular, we anticipate that placing the Dirichlet condition on any single species will confine the patterns to the interior via the same mechanism suggested in Sect. 2.2. We do note, however, that interactions with multistability can lead to more complex behaviour, and may require a choice of R different from what is stated in Sect. 1.2. For instance, if we took R = 0 in the examples shown in Fig. 11, we would see a similar localization of the interior pattern, but the region around it would have u 1 ≈ 0.7 (i.e. a maximal value) rather than u 1 ≈ 0. Lastly, we remark that while this system exhibits multistability, Turing instabilities only make sense around the coexistence equilibrium of all three species (as they cannot exist around two-species coexistence equilibria, and the extinction steady state would lead to negative values of at least one species). We leave exploration of more general multistable systems to future work.

Discussion
Isolation of patterns in reaction-diffusion systems (RDS) on the interior of a domain can be motivated both phenomenologically, by observing patterning fields which are clearly subsets of a given region, or mechanistically in terms of boundaries in gene expression. Such gene expression boundaries themselves can either emerge due to previous cell fate specification or local signalling dynamics. Here we have developed a novel method of inducing such isolated patterns with a simple change in the boundary conditions governing the RDS. We have justified these boundary conditions by considering heterogeneous RDS modelling a change in reaction dynamics explicitly and shown that such systems can reduce asymptotically to RDS satisfying these mixed boundary conditions at steady state, given weak geometric and kinetic constraints. This is consistent with the theoretical literature showing pattern localization and modulation in heterogeneous environments (Varea et al. 1997;Page et al. 2003Page et al. , 2005Kozák et al. 2019), as well as the regionalization of patterning seen experimentally (Johansson and Headon 2014). These conditions give a simple way of modelling these localization phenomena, without having to resort to modelling the heterogeneity that may have given rise to isolated patterning regions. We have further shown that isolated patterning can be obtained quite generically using these mixed boundary conditions across a range of reaction-diffusion systems, geometries, spatial dimensions, and even in systems of more than two interacting species.
As described in the introduction, there are still major gaps in our understanding of RDS, particularly as models in developmental biology. While the original formulation of Turing's theory is a powerful and simple way of obtaining periodically patterned states, it suffers from a number of robustness problems, in part due to this simplicity (Maini et al. 2012;Woolley et al. 2017;Scholes et al. 2019). In addition to demonstrating isolated patterning, our results indicate some level of robustness in pattern selection, partly explained in Sect. 2.2. Specifically, the numerically observed steady states in the case of homogeneous Neumann conditions can generally depend on initial data, as multistability of steady states is common, especially in two or more spatial dimensions (Borckmans et al. 1995). In contrast, there were far fewer such cases of multistability of patterns observed from simulations with the mixed boundary condi-tions (3). Of course, more work is needed to precisely explain this phenomenon, but it is consistent with the discussion of how bifurcation branches of solutions change under different boundary conditions in Dillon et al. (1994). Additionally, while we did not explore this in great detail here, the Turing space for these mixed boundary conditions can be enlarged, as discussed in Maini and Myerscough (1997), plausibly helping to overcome the constraints 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 simulation investigated, we observed that these boundary conditions led to a substantially faster convergence towards the steady-state pattern. In the case of labyrinthine patterns, these conditions often led to orders of magnitude faster convergence. Specifically, for the Schnakenberg kinetics with labyrinthine parameters in Table 1, the solution with conditions (3) was within plotting accuracy of the final steady state within O(10 2 ) units of time, whereas the solution with homogeneous Neumann conditions only attained the final stripe orientations shown in Table 2 after O(10 4 ) units of time. We conjecture that this is not, at least primarily, due to a change in linear stability of the final pattern, but is because of fewer steady states satisfying these boundary conditions, and hence less shadowing of unstable equilibria. These results are consistent with earlier results on inhomogeneous Robin conditions (Arcuri and Murray 1986). As with the above discussion of counting equilibria, we leave further investigation to future work.
Finally, we mention that our results also demonstrate a novel application of domain geometry in robustly giving rise to patterns of qualitatively different types. Specifically, in Fig. 6 we see that inverted spots can be selected on the interior of these domains, with stripes closer to the boundaries, due to the symmetry of the domain geometry. As far as we are aware, this is a novel local pattern selection mechanism quite distinct from those exploiting a reduction in dimension to transition from spots to stripes (Murray 1981(Murray , 2004. In general, the sensitivity of patterns with respect to geometry, and the use of domain geometry to influence patterning, has not been studied as extensively as other aspects of RDS, especially for more general boundary conditions. Many periodic patterns in development are spatially localized within some region of an otherwise homogeneous tissue. There are at least three distinct ways that such regionalization can occur, as described in Sect. 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 of a 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 one species, and whether it promotes the formation of spots nearby (and so the boundary lies close to the edge of 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 boundary is assumed to be repelling, without having precise knowledge of what the field is doing near or beyond the boundaries where periodic patterns are observed. The examples given in Sect. 3.2, and in particular the 'cut' domains shown in Fig. 9, provide ways of validating these isolating boundary conditions via experiments which alter domain geometry, as we have demonstrated strong effects of such geometry on the emergent patterns.
As mentioned in the introduction, Neumann boundary conditions can also exhibit isolated patterning if the initial data consist of internal spot solutions (Kolokolnikov et al. 2009;Chen and Ward 2011). It would be interesting to quantify if such internal states are generic in some sense, by comparing and contrasting their basins of attraction from states with boundary spots. Such a quantification of 'generic' behaviour would be valuable to inform questions of biological robustness of these mechanisms. There are also important open questions raised about pattern formation multistability with these mixed boundary conditions, and how one might choose the value of R in more general systems, as discussed in Sect. 3.3.
There are numerous other extensions of the basic ideas here, but we restrict attention to five possible areas of future work. As above, there is work to be done in more rigorously justifying our conclusions regarding equilibria, pattern selection, and transient dynamics.
As the examples in Table 1 all have large diffusivity ratios D, the shadowlimit approaches of Ward et al. (2002) (and many other references) could be used to provide further insight on interior localization. On the biological side, exploring patterning field boundaries with these kinds of models in mind could help us understand the key aspects of the detailed biological mechanisms at play. There are also important questions, raised in Esposito (2020) and elsewhere, about the underlying thermodynamics of RDS, and the role that boundaries play. We also remark that there are examples of more complicated multi-domain and bulk-surface models, where the role of boundary conditions and geometry have major impacts, and there are important connections between those models and the spatially heterogeneous ones studied in Sect. 2.1 . Finally, we have restricted our attention to studying stationary solutions, but RDS can exhibit a wide variety of spatiotemporal behaviours, such as spatiotemporal chaos and spiral waves observed in excitable systems (Sánchez-Garduno et al. 2019). Preliminary simulations of such dynamics suggest that our mixed boundary conditions can also lead to localized behaviour in these systems, with possible applications in electrophysiology and other areas [e.g. little is known regarding how boundaries influence re-entrant wave dynamics (Clayton et al. 2011)]. While there has been great progress in extending Turing's insights into pattern formation and RDS, there is substantial work left, especially in terms of elucidating the connections with real embryonic development.