Turing Patterning in Stratified Domains

Reaction–diffusion processes across layered media arise in several scientific domains such as pattern-forming E. coli on agar substrates, epidermal–mesenchymal coupling in development, and symmetry-breaking in cell polarization. We develop a modeling framework for bilayer reaction–diffusion systems and relate it to a range of existing models. We derive conditions for diffusion-driven instability of a spatially homogeneous equilibrium analogous to the classical conditions for a Turing instability in the simplest nontrivial setting where one domain has a standard reaction–diffusion system, and the other permits only diffusion. Due to the transverse coupling between these two regions, standard techniques for computing eigenfunctions of the Laplacian cannot be applied, and so we propose an alternative method to compute the dispersion relation directly. We compare instability conditions with full numerical simulations to demonstrate impacts of the geometry and coupling parameters on patterning, and explore various experimentally relevant asymptotic regimes. In the regime where the first domain is suitably thin, we recover a simple modulation of the standard Turing conditions, and find that often the broad impact of the diffusion-only domain is to reduce the ability of the system to form patterns. We also demonstrate complex impacts of this coupling on pattern formation. For instance, we exhibit non-monotonicity of pattern-forming instabilities with respect to geometric and coupling parameters, and highlight an instability from a nontrivial interaction between kinetics in one domain and diffusion in the other. These results are valuable for informing design choices in applications such as synthetic engineering of Turing patterns, but also for understanding the role of stratified media in modulating pattern-forming processes in developmental biology and beyond.


Introduction
Since Turing's initial insights into reaction-diffusion-driven morphogenesis (Turing 1952), a substantial research effort has elucidated various mathematical and biophysical aspects of such symmetry-breaking instabilities leading from homogeneity to patterned states (De Kepper et al. 1991;Cross and Hohenberg 1993;Maini et al. 2012;Kondo and Miura 2010;Green and Sharpe 2015;Woolley 2014). An important and well-studied aspect of these instabilities is the underlying geometry, which can influence both the stability of a homogeneous state, as well as the subsequent mode selection of emergent patterns (Murray 2003). However, one less well-studied aspect of geometry is the coupling between layered spatial domains, which can arise in a variety of settings and is the primary object of interest in this paper.
Reaction diffusion processes arise in a diversity of layered settings, from bulksurface membrane-cytosol interactions Rätz and Röger 2014;Kretschmer and Schwille 2016;Spill et al. 2016;Cusseddu et al. 2018) to epithelialmesenchymal couplings in developing skin (Cruywagen and Murray 1992;Shaw and Murray 1990). Synthetic experiments involving pattern formation in monolayers also exhibit clearly stratified regions of cells and culture medium (Sekine et al. 2018). Additionally, many experiments involving bacterial pattern formation are performed using colonies on the surface of a substrate, such as agar Berg 1991, 1995). Such systems either use natural chemotaxis mechanisms to initiate spatial pattern formation of the bacterial density itself Tyson et al. (1999), or instead use synthetic bacteria re-engineered to express additional quorum-sensing pathways that spatially coordinate patterns in gene expression (Basu et al. 2005;Tabor et al. 2009;Grant et al. 2016). Other examples are synthetically reconstituted protein interaction systems with bulk-membrane coupling such as the Min system (Loose et al. 2008;Kretschmer and Schwille 2016;Frey et al. 2018), where molecular interactions (Denk et al. 2018;Glock et al. 2019), or in vitro system geometries (Wu et al. 2016;Brauns et al. 2020;, are modified to stimulate changes in the observed protein patterns. Examples of particular contemporary interest include the use of bacterial colonies as exemplars of synthetic multicellular communication and self-organization (Balagaddé et al. 2008;Dalchau et al. 2012;Payne et al. 2013;Grant et al. 2016;Karig et al. 2018), for example using modified E. coli with engineered quorum-sensing signalling on the surface of an agar plate (Grant et al. 2016;Payne et al. 2013;Cao et al. 2016Cao et al. , 2017. Some of these systems take advantage of the geometry of colony growth and nutrient diffusion to influence pattern formation (Payne et al. 2013;Cao et al. 2016Cao et al. , 2017 while in other systems the bacteria are confined (Grant et al. 2016;Boehm et al. 2018), but the signalling molecules can diffuse into the inert agar layer below the chemically active colonies. The impact of this leaching on the prospects of a Turing instability in experimentally relevant geometries has not been fully characterized and is a key motivation for our study.
Our first objective is to develop a two-domain model of reaction-diffusion processes coupled in a stratified bilayer and to determine conditions for the Turing instability, on the assumption that the upper region is sufficiently substantive in the transverse direction to merit continuum modeling. Such a model is also applicable to a variety of other settings beyond multilayered bacterial pattern formation, such as developing skin. Our second objective is to focus on the Turing instability for multilayered bacterial systems, where signalling molecules only diffuse in the lower (agar) layer and especially where the upper layer is asymptotically thin relative to the scale of the pattern and the depth of the lower layer. The main biological motivation is to determine to what extent the diffusive bulk helps, or hinders, the ability of an engineered system to exhibit Turing-type patterning. In terms of model development, domain-coupled reaction-diffusion systems broadly fall into three major types: instantaneously coupled, bulk-surface models and bulk-bulk models. The first type are models where the components occupy the same physical space (or the reactions occur in thin regions where a homogenization approximation is sensible) (Yang et al. 2002;Epstein et al. 2007;Yang and Epstein 2003;Fujita and Kawaguchi 2013). Such models are essentially just larger reactiondiffusion systems with linear coupling between subsystems, and amenable to block matrix analysis in the study of Turing instabilities (Catllá et al. 2012), but do not capture the spatial separation of the domains. When applied to layered media, these models effectively assume vertical transport between distinct layers (such as that of Fig. 1) is instantaneous. However, considering physical scales representative of synthetic pattern formation experiments using E. coli (Grant et al. 2016;Boehm et al. 2018), and summarized below in Table 1, one has an agar block with a depth of a few millimeters, say three, and a diffusion rate on the scale of 4 × 10 −10 m 2 s −1 . Thus the timescale for vertical transport is in the region of 375 minutes, which is short compared to the timescales on which experimental measurements of the equilibrated system are recorded (1500-3000 min, Grant et al. (2016); Boehm et al. (2018)) but far from instantaneous. Hence, such models are inappropriate for the motivating examples here.
A second class of model considers bulk-surface coupling, where one component is confined to the boundary of the main bulk domain, and reactants flow between the two regions, such as in the case of proteins diffusing in the cytoplasm and binding on the cell membrane (Rätz and Röger 2014;Madzvamuse et al. 2015;Spill et al. 2016;Cusseddu et al. 2018;Paquin-Lefebvre et al. 2018;Frey et al. 2018). There is substantial recent interest in such models, from very theoretical results on existence and fast reaction-limiting behavior (Rätz 2015;Anguige and Röger 2017;Hausberg and Röger 2018), to spike dynamics (Gomez et al. 2018) and a myriad of applications to understanding cell polarity (Thalmeier et al. 2016;Kretschmer and Schwille 2016;Geßele et al. 2020). One particularly well studied example is the pole-to-pole Min protein oscillation in E. coli, which has the biological function of guiding the cell division machinery to midcell (Kretschmer and Schwille 2016). Such intracellular protein patterning systems have been studied experimentally and theoretically in a wide range of system geometries, such as spherical (Klünder et al. 2013;Levine and Rappel 2005), elliptical (Halatek and Frey 2012;Wu et al. 2016;Geßele et al. 2020), and planar membrane geometries . A striking feature of these examples is that the geometry itself has a major impact on pattern formation and pattern selection, which has been confirmed experimentally (Wu et al. 2016;Brauns et al. 2020). More generally, the Turing instability has also been studied in the context of membrane cytosol models (Rätz and Röger 2014;Madzvamuse et al. 2015). Overall, linear stability analysis (as used by Turing) is highly applicable to such membrane cytosol systems because the nonlinear interactions are typically restricted to the lower-dimensional membrane surface. The dynamics in the extended bulk are typically linear such that a general solution (or a good approximation) can be obtained analytically and used to satisfy the linearized reactive boundary condition. This is justified because the membrane can be considered as a surface with no transverse extent, and so transverse gradients only play a role in the cytosolic layer close to the membrane surface. However, in multilayered cellular systems, the transverse lengthscales are at least that of many cells, and hence transverse gradients cannot be neglected a priori and thus should be accommodated in the modeling. Models accounting for this represent the final class, with two separated spatial domains with an interface and suitable coupling boundary conditions. From the perspective of pattern formation, this kind of model has only been subject to recent numerical exploration (Vilaca et al. 2019), though it is used in the derivation of the second class-bulk surface models-given appropriate distinguished limits and scaling assumptions (for example, Chapman et al. (2016), Fussell et al. (2019).
Hence, we will develop models of the latter type, with an exploration of the conditions for the Turing instability, and their detailed study in the context of a stratified model with an inert underlying agar layer. Turing instabilities of reaction-diffusion systems have been studied on a variety of complex spatial domains such as compact manifolds (Varea et al. 1999;Chaplain et al. 2001), networks Ide et al. 2016;Nakao and Mikhailov 2010), and many of the aforementioned complex system geometries (Halatek and Frey 2012;Klünder et al. 2013;. The primary difficulty in such cases, compared to the textbook example of a continuous line, is determining the corresponding set of eigenfunctions and eigenvalues of the spatial transport operators, which for some system geometries do not need to coincide between domains (e.g., in the surface bulk elliptical case (Halatek and Frey 2012)). In such cases, approximate solutions for the system's eigenfunctions need to be derived that are orthogonal in the patterning layer. Examples that deviate even further from the classical case are growing domains (Crampin et al. 1999;Plaza et al. 2004;Krause et al. 2019;Sánchez-Garduño et al. 2019) and spatially heterogeneous reaction-diffusion processes (Benson et al. 1998;Page et al. 2003Page et al. , 2005Haim et al. 2015;Kolokolnikov and Wei 2018), for which the canonical approach does not work. In such cases, novel approaches to pattern-forming instabilities have recently been developed for growth (Madzvamuse et al. 2010;) and heterogeneity (Krause et al. 2020) under certain simplifications, but such analyses are quite different to the classical case. In a similar direction, as part of our objective in exploring the Turing instability for layered reaction-diffusion systems, we will aim to demonstrate a much richer diversity of structure in the resulting dispersion relations (and hence instability conditions), compared to classical counterparts. The surface (cellular) region denoted S has height H ε and contains both reaction and diffusion terms, whereas the bulk (nutrient) region B is of height H and is assumed to have no reactions, but permits diffusion. Both have lateral extent L, with no-flux conditions on all boundaries except for the interface between the two regions, where a coupling condition is applied (Color figure online) As an outline, in Sect. 2, we present a two-domain layered model, where each domain consists of closed two-dimensional rectangular regions, coupled through a single shared boundary (See Fig. 1) and briefly discuss how it can be reduced to a variety of other models. We focus on a special case of a two-domain model where we assume linear coupling and no reactions in the second (bulk) region, but note that our analysis can be applied with relatively simple modifications to more general cases. In Sect. 3, we develop an approach to linear stability analysis of homogeneous states. In Sect. 4, we derive a variety of asymptotic results regarding our dispersion relation, especially considering limits that are of particular relevance for synthetic pattern formation in E. coli colonies. We further explore these results and other parameter regimes numerically in Sect. 5. Finally, we discuss our results in Sect. 6.

Two-Region Model
We consider a layered two-domain model where each domain is governed by a different reaction-diffusion system. We consider several interacting species in these two domains, which we write as = S B where we refer to B = [0, L] × [0, H ] as the bulk region, and S = [0, L] × [H , H + H ε ] as the surface region (see Fig. 1). We writeû B ∈ R n for the concentrations of reactants in the bulk, andû S ∈ R n for the concentrations of reactants in the surface region. For simplicity, we consider a simple one-dimensional lateral geometry (orthogonal to the direction of the coupling condition), but note that the geometric details in the lateral direction(s) can be easily extended to much more complicated geometries, as long as eigenfunctions of the Laplacian in these directions are separable from the transverse coordinate y. We only consider reactions on the surface layer and assume the bulk only permits diffusion.
We have the following equations for the species concentrations in the bulk and surface regions: whereD S ,D B are positive definite diagonal matrices. We further specify Neumann (no-flux) boundary conditions on the outer boundaries as, and lastly coupling conditions on the interior boundary which conserve fluxes and take the form, whereĝ is a given function determining the transport between the surface and the bulk region, andη is a rate of transport across the boundary. Essentially, all of the forthcoming analysis can be carried out with a generalĝ, as linearization will also linearize this function. For brevity and concreteness, we will henceforth assume a linear transport law, so that we havê We non-dimensionalize the above model via concentration, time and length scales corresponding to the reaction kinetics and a unit lengthscaleL, respectively. Specifically, we defineû S = U u S ,û B = U u B , where U is a diagonal matrix of concentration scales. Equally, we sett = τ t, where τ is the timescale of the fastest reaction in the surface and bulk, andx =L x. The dimensional scalings are then chosen such that We define new dimensionless groupings h = H /L, ε = H ε /L, D S = τD S /(L 2 ), D B = τD B /(L 2 ),L = L/L and η = τη/L. The nondimensional system is written as There are several distinguished limits of the nondimensional system (8)-(12) that reduce the model to different cases already present in the literature. In the limit ε → 0, one can consider either scaling η ∼ O(ε) or scaling f S ∼ O(ε −1 ) in order to reduce the system to a bulk-surface model, which is well-studied in the literature (though primarily in radial geometries) (Levine and Rappel 2005;Rätz and Röger 2014;Rätz 2015;Madzvamuse et al. 2015;Cusseddu et al. 2018;Gomez et al. 2018;Paquin-Lefebvre et al. 2018). The second scaling, indicating that the surface timescale is rapid, can be related to assumptions regarding rapid surface reactions used to justify reactive boundary conditions from the microscopic viewpoint (Chapman et al. 2016). Finally, another limit is the case of infinite permeability, η → ∞, wherein the concentrations and fluxes are continuous across the interface. In this case, the system can be seen as a single-domain model with a step function heterogeneity, which has been studied extensively as an example of spatially heterogeneous reaction-diffusion systems (Benson et al. 1998;Page et al. 2003;Kozák et al. 2019). Nonetheless, pattern formation in the system above, as well as several other distinguished limits, has not been analyzed yet in the literature.
In Table 1, we give the dimensional parameter scales to be considered in our framework, taken from the key motivating example of synthetic patterning in E. coli bacterial colonies on an agar substrate (Grant et al. 2016;Boehm et al. 2018). While such experiments can be conducted with a variety of settings, an overall restriction on the variation of these parameters is motivated by the range of the physical scales in these studies. Here, bacteria are plated in squares of about 1 mm (Methods, Boehm et al. 2018) with patterning cells considered in an 8×8 grid in one study (Supplementary Information, Grant et al. 2016) and more generally the patterning fields are observed across about 22 such squares (Fig 5B, Boehm et al. (2018) andFig 3E Grant et al. (2016)). Thus we consider a range ofL ∼ 8 − −22 mm. For the diffusion matrices, the infinity (max) norm · ∞ is presented, i.e., the maximum value of the matrix's components. From Grant et al. (Supplementary Material, Tables S8, S9, Grant et al. (2016)), diffusion coefficients have been estimated in the range D S ∞ ∼ 10 −10 m 2 s −1 − 10 −9 m 2 s −1 by model fitting to exemplar results. As this is also the scale of diffusion (or slightly more than the scale) for the signalling molecule EGF in water (Nauman et al. 2007), the same scale is used for D B ∞ . Similarly, in the parameter fitting by Grant et al., a reaction timescale on the scale of the faster reactions is such that 1/τ ∼ 8.4 × 10 −5 s −1 − 10 −3 s −1 , with the range arising from the use of different model kinetics in parameter fitting. We further assume 10-50 layers of bacteria, with an E. coli bacterium size scale of about 10 −6 m−2 × 10 −6 m (Levin and Angert 2015), and hence a surface depth on the scale of H ε ∼ 10 −5 m−10 −4 m. Finally, the depth of the bulk is highly variable and easily changed upwards from the millimeter scale, and so H is taken with the range of 1 − 10 mm; estimates for the interfacial permeabil- Table 1 Numerical scales of various dimensional parameters and parameter groupings in SI units, based on patterning in synthetic pattern formation with E. coli bacterial colonies, using physical scales motivated by the studies of Grant et al. (2016) and Boehm et al. (2018) Parameter Range Justification 10 −10 m 2 s −1 − 10 −9 m 2 s −1 Table S8 Table 2 Numerical scales of various non-dimensional parameters and parameter groupings, motivated by the physical scales of synthetic pattern formation with E. coli bacterial colonies, in the studies of Grant et al. (2016) and Boehm et al. (2018). For matrices, the infinity (max) norm · ∞ is used, which is the modulus of the matrix component with largest magnitude. For the non-dimensional matrix Jacobian, this norm is taken to be of order unity as the timescale is non-dimensionalized relative to τ , a representative timescale associated with a fast reaction in the system. The non-dimensional lengthscale,L, is retained symbolically throughout the presentation to facilitate determining the impact of this scale, though it is unity for these scalings. The parameter scales ε * and h * as well as the range of ε 2 * /3 are presented as they will be important in the asymptotic analyses below

Parameter
Typical value/range ity,η, are currently unavailable. These dimensional parameter estimates generate the non-dimensional scales of Table 2, which will guide the asymptotic and numerical investigations presented below.

Linear Stability Analysis
For a linear stability analysis of homogeneous equilibria of (8)-(12), we require the steady states to this system, which arise from specifying so that the surface reactions determine the spatially homogeneous steady state concentration in both regions, and our simple constitutive choice of g implies that these concentrations must be equal. We will focus exclusively on the case of absence of reactions in the bulk, as motivated by the underlying inert agar layer in synthetic pattern formation with E. coli bacterial experiments (Grant et al. 2016;Boehm et al. 2018) and which requires only a root of the surface kinetics for there to be a steady state. We proceed by considering perturbations to this steady state of the form where |σ | 1, and in general the bulk and surface perturbations are n-dimensional functions, where n is the number of species. We substitute these perturbations into Eqs. (8)-(12) to find, from Eqs. (8)-(9), that the perturbations will satisfy where the Jacobian, J S = ∂ f S /∂ u S ∈ R n×n , is evaluated at the steady state concentrations. We also have the coupling condition from Eq. (12) given by,

Spatially Homogeneous Perturbations
We now consider the appropriate generalization of stability in the absence of transport, as is typically assumed in a Turing-type analysis. However, unless the reaction kinetics are the same in both domains, which is not true in our setting, then spatially homogeneous perturbations are not consistent with Eqs. (13)-(15). Such perturbations will not remain homogeneous under time evolution due to the coupling condition (15), except in the mathematically fine-tuned case where the homogeneous surface perturbation is along an eigenvector of J S with eigenvalue 0. Previous studies of more complex systems (beyond those considered in textbook Turing models) also highlight that, when generalizing the conditions that arise from the stability of the homogeneous steady state with respect to spatially homogeneous perturbations, one must also consider a perturbation with respect to the zero mode(s) of the transport operator (Klika et al. 2018). However, given the assumption of completeness, i.e., that separable solutions in x, y and t span the space of possible solutions, as generally used in linear stability theory, the existence of zero modes of the transport operator also requires mathematical fine-tuning. In particular, with ∇ 2 u 0 S = 0 for the zero mode of the transport operator acting in the surface layer, u 0 S , one has for a general A and q a natural number on enforcing the zero flux boundary conditions. There is a directly analogous expression for the zero mode of the transport operator within the bulk region. However, after rearrangement, the interfacial condition at y = h requires One possible solution occurs for k q = 0, which generates a spatially homogeneous mode that has already been considered above. Satisfying this equation for other k q requires mathematical fine-tuning as k q is already constrained to a set of zero measure and all other parameters are either geometrical or biophysical in origin.
Hence, to summarize, in contrast to the textbook Turing case, unless the surface and bulk kinetics are the same, constraints on the parameters do not arise from the constraint of stability to homogeneous perturbations; instead, this stability always holds, at least in the absence of a mathematical fine-tuning of parameters and the possibility of such fine-tuning is neglected below.

Spatially Inhomogeneous Perturbations
To proceed, we assume a separable solution in x, y and t for the linearized system (13)-(15). With the usual assumption of a uniform temporal growth rate λ, the perturbation ansatz for a single mode of a separable solution is where k q = qπ/L for q a natural number (including 0). Assuming completeness of such a set of modes, then linear superposition entails that an arbitrary function may be expanded via a weighted linear sum of individual modes. Hence, the question of stability of a linear perturbation reduces to the same question for single modes, as in the standard textbook analysis (Murray 2003), without the need for the modes to be orthogonal. Noting that homogeneous modes are not feasible, as shown above, we proceed to consider whether any heterogeneous modes exhibit instability (R(λ) > 0). Furthermore, we note that even in the absence of completeness, R(λ) > 0 still provides a sufficient condition for instability, though it is not strictly necessary. We will see later that our conditions are not refuted by comparisons with numerics, and so we anticipate that the set of modes we construct is at least generic if not a complete basis. We note that for each q, there may be many distinct λ, and corresponding to each distinct pair of (q, λ) we will have possibly different eigenfunctions s, and b. We will suppress this dependence in the following, but it is important to keep in mind that the following analysis applies for a given pair (q, λ). As we are looking for modes which grow in time, leading to instability, we will impose R(λ) > 0 in the following. Substituting these expansions into (13)- (14), we find that a given mode satisfies, After multiplying these equations by the inverse of the diffusion matrices and rearranging, we find where denotes the ordinary derivative with respect to y. These spatial functions are required to satisfy the external boundary conditions b (0) = s (h + ε) = 0 and the coupling conditions which read, To find suitable (λ, q) that solve the coupled problem (19)-(21), we will make use of the matrix-valued function defined by cosh(M) = (exp(M) + exp(−M))/2, for some matrix M, as well as sinh(M) = (exp(M)−exp(−M))/2. We recall the differentiation identity cosh(y M) = M sinh(y M), which follows from this definition. We now seek to take the square root of the matrices on the right-hand side of Eqs. (19) and (20) and thus define We next consider solutions to Eqs. (19) and (20) via hyperbolic matrix functions. As we will observe (e.g., Eq. (26) and the resulting dispersion relation), these matrices will always be in terms of functions that can be expressed in terms of even powers of M B and M S , and thus functions of M 2 B and M 2 S . This dependence on the squares of these matrices follows as if f : C → C is analytic, then a matrix-valued function can be defined via a power series in the matrix argument (Higham 2008).
The hyperbolic functions we will use are meromorphic with poles away from 0, and hence the ambiguity in defining the square root matrices, M S and M B does not play a role. Without loss of generality, we will consider the principal square roots of the matrices for definiteness, so that eigenvalues of M B and M S are the square roots with positive (or possibly zero) real parts of the eigenvalues of M 2 B and M 2 S . Proceeding, we then have the following solutions to Eqs. (19) and (20) given by the hyperbolic matrix functions: for some nonzero constant vectors α, β. We note these functions satisfy the no-flux conditions at the top and bottom boundaries by construction. We now use the coupling conditions (15) to determine a condition for nontrivial α and β. These read, We then have, writing Eqs. (24) and (25) as a 2n × 2n block matrix, the following condition for nontrivial solutions to this system: As this condition involves transcendental functions of λ, we note that in general for a fixed spatial mode q, there will be infinitely many values of λ for which Eq. (26) is satisfied. Equivalently, q only differentiates between eigenmodes in the x direction, but cannot do so in y, and so these eigenmodes must be captured via multiplicity in λ. While the condition given by Eq. (26) is in principle computable, it is difficult to use to gain insight into Turing-like instabilities. Even simplifying the determinant condition is nontrivial, as the four blocks will not in general commute, so we now exploit the assumption of no reactions in the bulk to simplify this condition. We have that M 2 B is diagonal, and from our assumption that R(λ) > 0, we have that its eigenvalues have positive real part. Therefore, the elements of cosh(M B ) are given by the hyperbolic cosine of the diagonal elements of M B , and since these are all positive By the above argument, we have that B is invertible. We then have that (26) can be written (by exchanging rows and using the Schur complement) as, Noting that (28) We note that the Turing instability conditions for the surface system in isolationneglecting spatial structure in y-are precisely that the growth rates λ computed from det(M S ) = 0 have negative real part for k 0 = 0, and positive real part for some k q > 0, and so this matrix encodes directly the classical case in this way. Furthermore, for a fixed q, and with fixed model parameters, we expect that condition (28) admits infinitely many distinct values of λ. The intuition for this is that in the uncoupled case (η = 0), the surface domain is a rectangle, and hence, the surface eigenfunctions s(y) are also cosines of different spatial eigenvalues, which can vary independently from k q . However, we know of no method to compute analytical expressions for such spatial eigenvalues in the coupled case, and so instead use condition (28) to compute λ directly, remaining aware of the inherent multiplicity. To further understand the dispersion relation given by (28), and how it relates to classical conditions for Turing instabilities, we now pursue several asymptotic reductions.

Instability Conditions in Thin-Surface Regimes
In this section we compute instability conditions from Eq. (28) for a variety of distinguished limits modeling a thin-surface region, as motivated by synthetic patterning in bacterial populations. First, we mention even simpler reductions of the system, as a consistency check of our dispersion relation. We show that patterning is equivalent in the limit of decoupling the interaction of the surface and bulk regions, that is for sufficiently small η 1. This is pursued in Appendix A, where the classical Turing conditions are recovered as the surface system becomes isolated, as required. In addition, in Appendix A, we also demonstrate that no patterning can occur for classical Turing kinetics once all diffusion coefficients are equal in each of the regions, a direct analogue of the well-known result that the classical Turing instability requires differential transport.
Noting that the full system is too rich to investigate in generality and that the non-dimensional surface depth parameters, ε and ε * are small in Table 2 for the motivating example of synthetic pattern formation in E. coli colonies, we proceed below to studying pattern formation instabilities with thin-surface asymptotics. In the experimental setting of Grant et al. (2016), the bacterial layer is always relatively thin, owing to transport constraints in the bacteria, though the agar layer can take different bulk heights. For this reason, after first introducing a thin-surface limit of the dispersion relation (26) in Sect. 4.1, we consider subsequent limits of large or small bulk thickness, h, in Sect. 4.2. We anticipate that the permeability of the interface, η, is large in these experiments but do not have quantitative estimates, and so also consider our asymptotics across varying values of this parameter. In Sect. 4.3, we derive asymptotic results under a regular asymptotic assumption on k and λ (i.e., that they remain comparable with non-asymptotic terms in the dispersion relation), and collect these results in Table 3. Finally, in Sect. 4.4, we give an example of distinguished limits where this asymptotic assumption breaks down. Throughout the following, we implicitly assume that the surface Jacobian, J S , has elements that are of the same order and thus of the order of J S ∞ , so that J S A ∞ is of the same scale as J S ∞ A ∞ for any matrix A considered.

Thin-Surface Limits
Here, we consider an asymptotically thin surface, requiring ε M S ∞ = H ε M S ∞ /L 1. First note that in the thin-layer limit below, the surface Jacobian J S only appears via In addition, given patterning (i.e., R(λ) > 0), the matrix J S cannot be dominated by the terms λI n or k 2 q D S within D S M 2 S , since then the reaction kinetics are subleading in the requirements for patterning, which thus contain only terms associated with pure diffusion at leading order. However, pure diffusion cannot induce patterning, as demonstrated in Appendix B. Thus we conclude that, given patterning and also that ||M S || ∞ has an upper bound ( and in particular the k 2 q term is in fact bounded). Noting the boundedness of M S we have that cosh(ε M S ) is invertible, as for ε sufficiently small this matrix has a determinant which is asymptotically 1 + ε 2 trace(M 2 S )/2 > 0. In addition, for sufficiently small ε, we have the Taylor expansion where O(ε 2 M S 2 ∞ /3) means the same scale as ε 2 M S 2 ∞ /3, or smaller, as the surface thickness tends to zero. Thus by right multiplying (26) by cosh(ε M S ) −1 and Taylor expanding, we obtain (to leading order) the relation providing ε 2 M S 2 ∞ /3 1 (ensuring the invertibility of cosh(ε M S ) and the validity of the Taylor expansion above). Furthermore, noting that ε * = ε D −1 S J S 1/2 ∞ together with the relations (30), which give the maximum scale of k 2 q and show that |λ| The latter inequality is immediate in the two species case on writing D S = diag(a, aξ) with ξ ≤ 1, as then D S −1 ∞ = 1/(aξ) ≥ 1/a = 1/ D S ∞ , with a trivial generalization to higher number of species. Hence, for conditions associated with patterning, the relative error in the leading-order thinsurface approximation arising from Eq. (31) is ε 2 * /3, and thus, we require ε 2 * /3 1. Despite the very large range of potential parameters in Table 2, the scales for synthetic patterning in bacterial colonies are consistent with this bound.

Consideration of Bulk Depth h
Noting D S ≈ D B at least for the parameter estimates of Tables 1 and 2, and also relations (30), (33), we also have Hence, an appropriate scale for the largest ∞ , which ranges from small to large in Table 2 For large h * simplifications, first note that M 2 B is diagonal, with diagonal components that have positive real parts since R(λ) > 0 as we require instability. Furthermore, similar to the synthetic patterning explored in experimental studies (Grant et al. 2016;Boehm et al. 2018), we are interested in lateral patterning (in the x direction of Fig. 1), thus, we take k 2 q > 0 and enforce k 2 q ≥ π/L by wavemode selection, which bounds the real part of M 2 B away from zero. For z ∈ C with R(z) = 0, we have the limit z tanh(z) → Sign(R(z))z as |z| → ∞, as may be deduced by writing z in terms of its real and imaginary parts, with subsequent use of the properties of trigonometric and hyperbolic functions. In addition, we have, without loss of generality, defined M B by the diagonal matrix with positive semidefinite real part for the square root of the diagonals of M 2 B , and in fact no such square root has zero real part since the diagonals of M 2 B have positive real part. Consequently, at leading order we have in the large h Finally, with this definition of M B , which is diagonal with terms whose real parts are bound away from zero, we also have that the diagonal elements, and hence the matrix norm, do not blow up on taking the hyperbolic tangent (all of its singularities lie on the imaginary axis) and thus tanh(h M B ) ∞ ∼ O(1) for h * ∼ O(1). This may be summarized together with Eqs. (35) and (36) via We are now in a position to consider the small ε * , thin surface, limit of the instability condition given by Eq. (28), considering the full range of values of h * , which is a measure of the non-dimensional depth of the bulk relative to the patterning lengthscale.
We also consider the case h * ∼ O(ε * ) for relative completeness, even though Tables 1 and 2 highlight that h * ε * is anticipated for experiments with synthetic pattern formation within bacterial populations.

Thin-Surface Asymptotic Regimes with k
An example of patterning when k 2 q D S ∞ , |λ| ord( J S ∞ ) is given in the next subsection, but here we consider thin-surface asymptotics with ε 2 . Hence, we are considering pattern formation that occurs on the timescales of the kinetics with a lengthscale associated with the timescale of the kinetics and the (largest) diffusion scale and, as previously noted, this simplifies the instability condition (28) at leading order to where the scale of the non-dimensional permeability η is unknown, and the possible values Hence, there are several nontrivial distinguished limits, which we proceed to document. Where possible, we will also relate these limits to the isolated surface case, where λ is determined by the dispersion relation det(M 2 S ) = det(λI n + k 2 q D S − J S ) = 0, in order to understand the impact of the bulk on the classical single-domain situation.
Case I h * ε * 1: Noting that h * ε * is equivalent to h ε by definition, in this limit Eq. (38) reduces to, However the determinant with the hyperbolic tangent term cannot generate a root with R(λ) > 0 and thus patterning. In particular, in Appendix A.2, following Eq. (53), it is shown that when R(z 2 ) > 0 one also has R(z tanh(z)) > 0. With R(λ) > 0 for patterning, let z 2 = h 2 (k 2 q +λ/d B ), where d B is a bulk diffusion coefficient. Thus z 2 is an eigenvalue of h 2 M 2 B , and all eigenvalues of this matrix are of this form. Furthermore, we have R(z 2 ) > 0 where z is an eigenvalue of h M B and satisfies R(z tanh(z)) > 0. However, for the hyperbolic tangent term in Eq. (39) to generate a root, at least one eigenvalue of h M B , that is one such z, must satisfy R(z tanh(z)) < 0, a contradiction, thus showing there are no roots from the determinant involving the hyperbolic tangent. Hence, noting D S is positive definite, the only roots are those of the isolated Turing modes, independent of η, and determined purely from det M 2 S = 0.
Case II ε * 1, h * /ε * = h/ε =ĥ ∼ ord(1): This limit corresponds to the entire domain being thin with respect to the lengthscale in the x direction (L). In this case, we have, Equation (40) is a slight modification of the isolated surface Turing conditions in 1D given by det(M 2 S ) = 0, and can similarly be written as an nth-order polynomial in λ. Further, if η ε D B M 2 B ∞ ∼ ord(ε J S ∞ ), then the conditions for instability are precisely those for an isolated surface. Similarly, if η = ord(ε D B M 2 B ∞ ) ∼ ord(ε J S ∞ ), then we are left with a "quadratic" dispersion relation, which does not simplify from the form given in (40) ("quadratic" meaning this dispersion relation will give a polynomial of order 2n for λ, compared to the standard nth order polynomial). In general such a relation could lead to quite different values of λ from the isolated case, though we will not analyze it further here. If η ε D B M 2 B ∞ ∼ ord(ε J S ∞ ), we then have the instability condition, which can be seen as a homogenization, or averaging, of the bulk and surface layers. Such an averaged dispersion relation has the potential to increase the ability of the system to pattern compared to the isolated case by, e.g., introducing, or increasing, the differential diffusion between species.
In some other (experimentally relevant) cases, this averaged system will decrease the ability of the system to pattern compared to the isolated case. For instance, the necessary differential diffusion for Turing patterning may be due to, for example, substrate binding (Korvasova et al. 2015) that is only present in the surface system. In an inert bulk region, there are fewer physical scenarios where differential diffusion is likely as most biological proteins are roughly the same size. In such a case, we have that D B = c B I, so that (41) can be rearranged to given, which we can see as a shrinking and shifting to the left a root λ coming from the isolated case. Effectively then, such a scenario leads to a smaller instability region in parameter space, subject to the wavemode selection constraint that k q = qπ/L, for a natural number q.
Another plausibly relevant case of Eq. (41) is if D S = D B , i.e., the surface and bulk diffusivities are the same. Here, the dispersion relation is that of the classic case except λ and k 2 q are both scaled by (1 +ĥ). Hence the allowed values of (λ, k 2 q ) are those of the classic case divided by (1 +ĥ), which shrinks the range of the allowed patterning wavenumbers relative to the classic case and thus leads to a smaller Turing space compared to the isolated surface system, though again subject to the wavemode selection constraint.
Noting D S −1 ∞ D S ∞ ≥ 1, as deduced just below Eq. (33), and ε * = Hence the final term from relation (32) The hyperbolic tangent does not contribute to instability, by analogous reasoning to Case I. Thus, there is no instability unless η is concomitantly small alongside ε D S M 2 S ∞ ∼ ε J S ∞ . Writing out M 2 S , we see that the impact of the bulk on the surface system is simply to shift the eigenvalues to the left in the complex plane by the quantity η/ε, and hence the Turing space for this system is strictly smaller than the Turing space for an isolated one-dimensional system with surface kinetics.
Collecting all of these various limits together in Table 3, we can see a pattern emerging. As η or h * is increased, we observe a trend of moving from the isolated surface system to a reduced, or average system, and eventually, for h * /ε * ∼ h/ε 1 and η ε J S ∞ , to no patterning being permitted. While it is not true in general that the averaged case or the quadratic case correspond to a reduced ability for a system to pattern, we anticipate that this is the case for most standard Turing systems, and hence there is a broadly monotonic decrease on the ability of a system to pattern as the bulk becomes larger or the boundary more permeable. This has concomitant implications for prospective multilayered Turing systems, such as the experimental studies involving bacterial patterning which motivate this study. Table 3 Thin-surface limits obtained in different asymptotic regimes given k 2 q D S ∞ , |λ| ∼ ord( J S ∞ ). Note that moving left to right corresponds to an increasing size of h * , and moving top to bottom corresponds to increasing scales of η.
We proceed by noting that the first and third terms of Eq. (44) are ord(ε J S ∞ ) and the second is ord ε 2 J S 2 ∞ /η . Writing We then have h D B M 2 B = εĥ(K 2 q D B + μI n ), and can factor an ε from Eq. (44) to obtain, which, in general, can admit nontrivial instabilities due to the coupling of the surface and the bulk. In particular, when η ε J S ∞ , so that the second term is no longer retained in the leading order, we find that the growth rates μ are given as the eigenvalues of J S /ĥ − K 2 q D B . This matrix resembles the classical isolated-surface case except with a scaling of the kinetics byĥ and the appearance of the bulk diffusion parameters, rather than those in the surface. Hence, we can use usual methods (e.g., the Routh-Hurwitz criterion) to determine parameters that lead to instability in this case, noting that any values of λ associated with instability will be of modulus ord(ε 1/2 ), and hence will be associated with slow growing modes. Additionally, we anticipate that such modes will also exhibit small amplitude patterns, as is typical due to center-manifold reduction near Turing-type bifurcations (Cross and Hohenberg 1993), and hence may not be visible against experimental noise. While other distinguished limits may exist which do not fall into the classifications given in Table 3, for brevity we do not pursue a systematic classification of these here. In the next section, we will show that almost all numerically computed dispersion relations given by condition (28) fall within the asymptotics given in Table 3, with the exception of the case given in Eq. (45) which was found numerically first, and subsequently motivated the above scaling.

Numerical Exploration of Example Systems
As an example of these dynamics we consider the Schnakenberg kinetics for surface reactants u S = (u S , v S ) given by The spatially homogeneous steady state is given by u * S = u * B = a + b, b/(a + b) 2 . Unless otherwise stated, we will assume equal diffusion coefficients between the surface and the bulk given by the diagonal matrices D S = D B = diag(d u , d v ). Without bulk reactions, and given linear interfacial conditions as summarized by Eq. (12) with the relations (6) and (7), we can immediately apply condition (28) to determine whether, or not, we expect a solution to pattern, and then compare these predictions with numerical simulations of the full nonlinear system.
Numerically computing λ from condition (28) is substantially more involved than typical Turing-type analyses (e.g., for polynomial dispersion relations (Murray 2003) due to the transcendental nature of this determinant condition. In particular, we expect that for any given wavemode in the x direction given by k q = qπ/L, for a natural number q, we have infinitely many distinct values of λ. These essentially correspond to the wavemodes in the y direction which we have found only implicitly in our construction of the dispersion relation. So to determine if, for a given set of parameters, condition (28) admits a value of λ with R(λ) > 0 we resort to numerical heuristics. While fast general-purpose methods exist for root finding of polynomials over the complex numbers (Verschelde 1999), we are unaware of similar methods for more complicated functions. In lieu of this, we developed a set of numerical heuristics to accurately determine whether or not a value of λ with R(λ) > 0 exists, and tested this against full For c, we anticipate there is an instability for relatively low k q /π and large η due to surface-bulk interaction instabilities, as exemplified in Sect. 4.4 for h ∼ ord(ε 1/2 ) and ε 2 * /(3ε 1/2 ) 1 (Color figure online) numerical simulations. We make use of the MATLAB function PatternSearch as well as a deflation algorithm based on Muller's method to find many candidate roots with positive real part (Muller 1956;Conte and De Boor 2017), and then discard any which are spurious. Throughout this section, we denote the largest such root by max(R(λ)), noting that even in the classical case this maximum is needed as there are generically n distinct values of λ. We first consider numerical constructions of dispersion relations for the small asymptotic limits described in the previous section. Here we consider R(λ) as a continuous function of the spectral parameter in the x direction, k q , as is commonly done (Murray 2003). For ε, ε * 1, we have that the isolated reaction-diffusion system can admit growth rates λ comparable to a classical one dimensional reaction-diffusion system, which we will denote by λ C (which can be computed in the standard way (Murray 2003)). We can then consider the maximum value of R(λ) (across all values of λ found from condition (28)), and compare this to the isolated case. We have confirmed these dispersion relations against full numerical simulations by simulating on a domain of lateral sizeL such that a particular mode k q = qπ/L is admissible, and observing a patterned solution.
We plot these dispersion curves in Fig. 2 for a variety of the geometric and coupling parameters. As anticipated, the coupling strength η and geometric parameters h, h * and ε, ε * each influence the shape of these dispersion curves greatly. We now compare these curves to the predictions in Table 3. For Case I (h ε), we see that max(R(λ)) is almost unchanged to the standard case up to small corrections not captured by the asymptotics. In Case II (h ∼ ε), we observe approximate equivalence of the dispersion curve to the isolated case for small η, and an apparent change in the dispersion relation for increasing η. The Case III behavior (h ε) is consistent with the asymptotics of Table 3 whenever R(λ) > 0 except for η ∼ ε J S ∞ and η ε J S ∞ at relatively small values of k q /π . Given these constraints, this mismatch is anticipated to be due to the interaction between the surface kinetics and the bulk diffusion, as described in Sect. 4.4 given the thin surface approximation ε 2 * /(3 1/2 ) 1 with k 2 q || D S || ∞ , |λ| ∼ord(ε 1/2 || J S || ∞ ). As a consistency test of this suggested mechanism, in Fig. 2d we replace D B by a scaled identity matrix so that differential diffusion in the bulk is no longer present, and we see that all of the dispersion curves, for smaller values of k q /π and η sufficiently large, fall below the axis as expected. This is true for different scalar multiples of the identity, such as D B = d v I 2 where the dispersion curves were even more stable. We remark that considering other parameters demonstrates that this nontrivial bulk-surface interaction can lead to a non-monotonic behavior of the dispersion relation with respect to η.
As in the classical case, we expect that for sufficiently large domains, any region where R(λ) > 0 should admit a patterned state. We confirmed this usingL = 100 for each of the dispersion curves, finding that they admitted patterned solutions for long time simulations if and only if R(λ) > 0 for some region in k q -space. Similar to the classical case, the layered model is always observed to stable at k q = 0 though with a local maximum at this point, in contrast to the behavior of the classical Turing instability dispersion relation.
To compare these dispersion relations against numerical simulations of the full nonlinear system, we compute a heterogeneity functional determining how far a solution is from a homogeneous state (Berding 1987). For simplicity, and because the surface layer is of primary interest in synthetic pattern formation within bacterial colonies, we only consider the heterogeneity of the activator in the surface. We define the heterogeneity functional as where c > 0 is simply a positive definite (dimensional) scaling parameter. Note that F h (u S ) ≥ 0 and for u S ∈ C 1 , F(u S ) = 0 if and only if u S is spatially homogeneous. While we do not anticipate this metric to be quantitatively comparable to max(R(λ)), we note that near the boundary of a Turing instability, the amplitude of patterns and their growth rates in time both scale with the distance from the bifurcation point, typically as a square root of the growth rate (Cross and Hohenberg 1993). Hence this functional should at least qualitatively scale with the growth of max(R(λ)) near the onset of instability. The value c is taken so that F h (u S ) = max(R(λ)) when η = 0 for scaling purposes. We note that these plots are intended to demonstrate qualitative, rather than quantitative, behavior near the onset of instability. In particular, we anticipate quantitative disagreement between F h (u S ) and max(R(λ)) when η is large, though the functional will still indicate whether or not max(R(λ)) predicts pattern formation, as well as the scaling of pattern heterogeneity as a function of max(R(λ)) near the onset of instability.
To use this heterogeneity functional, the full system (8)-(12) was solved until a final time of t = 10 5 to ensure a good representation of the steady state pattern. The initial data were taken to be u 0 = u * (1 + ξ u (x, y)) and v 0 = v * (1 + ξ v (x, y)) with ξ u and ξ v random fields such that at each value of (x, y), they are independently and identically distributed normal random variables with zero mean and variance 10 −4 . The equations were simulated using the COMSOL Multiphysics ® software (http:// www.comsol.com) with at least 2 × 10 4 second-order triangular finite elements. A non-uniform mesh was constructed such that the surface region S was resolved with at least 10 distinct triangular elements in any vertical cross section. Convergence was checked in spatial and temporal discretizations, and a relative tolerance of 10 −5 was given to the adaptive timestepping algorithm.
In Fig. 3 we give examples of this heterogeneity functional across the ranges of the geometric parameters ε, h, and η, alongside predictions from the instability condition (28). As anticipated by the asymptotics, for very small ε (Fig. 3a), we see the system fails to support spatial patterns for η ≥ 3.9 × 10 −4 . Additionally, we see a jump in the value of the heterogeneity between η = 8 × 10 −5 and η = 10 −4 . We plot values of u B in Fig. 4 across this jump to demonstrate that this discontinuity in the value of the spatial heterogeneity F h (u S ) for these parameters is due to different nonlinear modes emerging as parameters are varied, and so it is sensible that it is not captured in the linear analysis. Other discontinuities in the plots of the heterogeneity functional in Fig. 3 are similarly due to different patterned states being selected, and we do not further explore pattern multistability or dependence on initial data here.
In Fig. 3b, d, we see a region of intermediate values of η for which no patterning occurs, and more broadly across all of Fig. 3 we see that a minimal value of max(R(λ)) occurs approximately for η within the range (10 −3 , 1). We show examples of the mode selection process from Fig. 3d in Fig. 5. For small η = 10 −4 , we see stable multiplespike solutions that are essentially confined to the surface. As η increases further to 10 −1 , a single-spike solution is observed, at a smaller amplitude as the dispersion relation has just crossed the instability threshold given in Fig. 3d. Further increases to large η lead to stable spike solutions that remain essentially vertically homogeneous in the surface, but have small transverse variations in the bulk due to the change in reaction kinetics across the interface, as illustrated for η = 10 5 . Further increasing η sharpens these spike solutions across the domain, but does not impact the number of modes. Besides the discontinuities in the heterogeneity due to nonlinear mode selection, there is often a good match between the linear analysis (e.g. value of max(R(λ)) and the heterogeneity, which can be expected near to the Turing bifurcation points in simpler settings due to the existence of normal forms of the pattern amplitude (Cross and Hohenberg 1993). In all of Fig. 3 we observe that max(R(λ)) appears to asymptotically approach a fixed value for either η → 0 (which corresponds to the static Turing conditions) or η → ∞, with the latter always being smaller than the former, though this may just be a feature of the parameters explored here. However, in 3c (and the other cases noted in the caption), we observe that an instability occurs for all values of η, which is confirmed by numerical simulations of the full system.
As a further example which helps visualize the impact of varying the geometric parameters and coupling constant η, we observe patterns primarily confined to the surface but with some interaction with the bulk in Fig. 6. Again some mode selection effects are present (two vs three spot solutions for small and larger values of η respectively), though due to generic aspects of multistability in two spatial dimensional systems (Dewel et al. 1995), we suspect these depend somewhat on initial data, rather than just parameter values. Finally, in Fig. 7, we give an example where no Fig. 4 One-dimensional plots of u S corresponding to parameters in Fig. 3a for two values of η in the top two panels, and plots of the corresponding u B below (withL = 1 in all cases). The surface concentration u S is effectively homogeneous in the y direction, and so is essentially a one-dimensional pattern, shown above. Note that the bulk concentrations are almost homogeneous, whereas the surface concentrations are not (compare the scales of u S and u B ) (Color figure online) Fig. 3d for three values of η, andL = 1. Here, ε = 10 −2 and h = 10 −1 (Color figure online) change in the number of unstable modes was apparent for variation in η, though the structure of the solution does change.

Fig. 5 Plots of u S and u B corresponding to parameters in
Within Turing-unstable regimes, the surface largely drives the structure of the modes and hence the patterns can be thought of as quasi-one-dimensional (Figs. 6, 7). The permeability η does control how much structure there is, both in the bulk in general and in the surface modes' variation in the y direction, though in all cases the largest spatial variation is along the lateral coordinate x. For the largest permeability we explored (η = 10 5 ), we see that the sizes of the surface and bulk can have a significant impact on the relative shape of the solutions in the bulk region (cf Figs. 5c, 6c, and 7c). In particular, we see that the deepest part of the bulk (y = 0) in Figs. 5c and 7c maintain a fairly distinct periodic pattern between high and low activator concentrations, whereas the larger bulk in Fig. 6c is substantially more homogeneous at y = 0. We also note that for intermediate values of η, Figs. 6b and 7b have the largest visual gradients in the activator in the surface layer, consistent with the intermediate-η values having significant impacts on the predicted values of max(Rλ) in Fig. 2. This further demonstrates nontrivial impacts of the bulk geometry on the structure of emergent patterns, and such leeching into the bulk may be useful to help quantify its impact in synthetic systems.

Discussion
Motivated by recent interest in a range of biological contexts, we have developed and analyzed a general class of reaction-diffusion models of pattern formation in stratified media, though with an absence of reactions in the bulk and a linear coupling between the layers. We have derived a criterion for pattern-forming instability in such media, given by Eq. (26). In Appendix A we showed that the absence of differential transport within each layer entails no patterning for these systems, in direct analogy to the classical Turing instability. We have also demonstrated a range of interesting behaviors via asymptotic reductions in thin domains, and numerical simulations. In particular, this setting of a linearly coupled system with no reactions in the bulk with a thin-surface layer is also of significant biological interest, as several groups are using bacterial colonies on inert substrates as a medium for engineered pattern formation via synthetic biology (Grant et al. 2016;Boehm et al. 2018;Karig et al. 2018). However, as far as we are aware, there is little theoretical understanding of how the inert substrate impacts the surface reaction-diffusion systems in these kinds of geometries. Additionally, to accurately model the real complexity of these experimental systems we would need to account for intracellular (i.e., non-diffusible) proteins which play a role in the reactions, as our reaction-diffusion framework only captures the dynamics of diffusible signalling molecules.
Nevertheless, even in the simplified setting of an inert bulk and a thin surface, the computed instability criteria are much richer than in the classical case. For instance, the nine distinguished limits for |λ|, k 2 q D S ∞ ∼ ord( J S ∞ ) given in Table 3 demonstrate a variety of behaviors not predicted by analyzing the surface reaction-diffusion system alone, as is typical in applications. In addition, these distinguished limits, though emergent from a complex multi-parameter system, depend on only three nondimensional parameter groupings, ε * , h * and η/(ε J S ∞ ). The first two of these, respectively, are the surface and bulk depth relative to the lateral lengthscale, i.e. the basic geometry. The final grouping is η/(ε J S ∞ ) = τη/(H ε J S ∞ ). Noting that τ is chosen such that J S ∞ ∼ ord(1), one can deduce more generally that τ/ J S ∞ is the dimensional timescale of surface reaction. Hence, the final parameter grouping is the ratio of the interface permeability to the surface velocity scale, J S ∞ H ε /τ , with the latter in turn given by the ratio of the surface depth and reaction timescale.
We further note that our instability condition (28), recovers the usual features of Turing instabilities, such as requiring differential diffusion for their onset, and reducing to the polynomial dispersion relation when the bulk becomes uncoupled from the surface. The explicit coupling between bulk diffusion and surface reactions given by (45) when |λ|, k 2 q D S ∞ ∼ ord(ε 1/2 J S ∞ ) suggests additional distinguished limits from those in Table 3; the associated instabilities possess slower growth rates, but nonetheless highlight substantial and non-trivial impacts of the bulk on the system's ability to pattern. We anticipate that there are other examples of nontrivial surfacebulk coupling-driven instabilities, as suggested in the discussion of the Averaged and Quadratic cases in Table 3, but leave investigation of these to further work.
Broadly, our asymptotic and numerical results on thin surfaces suggest that the presence of the inert bulk generally decreases the ability of the surface system to undergo a Turing instability compared to an isolated system. The exceptional cases, such as the homogenized limit (41) and the explicit coupling in Eq. (45), can in principle lead to larger Turing spaces, though we have shown in some realistic cases such as equal bulk diffusions ( D B = I n ) that these do not enlarge the Turing space. Note that in systems where diffusion varies significantly between domains (e.g., nondiffusible proteins in the surface) the parameter space that admits pattern formation can increase with increasing bulk size (see, for instance, Brauns et al. (2020)). Exploring such interplays will be the focus of future work.
Our results suggest that experiments should aim to design large and robust parameter regimes using classical criteria for pattern-formation (e.g., using design approaches such as in Dalchau et al. (2012)), as diffusion into the bulk region will likely decrease the size of such Turing spaces. We have shown that even in cases where the broad influence of the bulk is to decrease the ability of the system to pattern, such a decrease will be non-monotonic in the geometric and transport parameters of the bulk region in general, as illustrated with the non-dimensional bulk depth, h and permeability, η. Many of the parameters may not be controllable, though one can often choose an agar height h above a certain minimal threshold. The results in Table 3 broadly suggest that the agar layer should be made as thin as possible to limit the impact on a system's ability to pattern. There may also be opportunities to decrease the permeability into the bulk, η, by using thicker filter paper or modifying the pore size or density, which would also reduce the negative impact of the bulk on pattern formation, though due to metabolic constraints (as the agar is primarily a nutrient) this too may be somewhat limited. We do note that there are important experimental controls in the genetic circuits encoded in the nonlinear reaction kinetics, which we have only caricatured in this study by considering the two-species case with only diffusible morphogens. Finally, we have shown instances of instability such that the bulk domain is a necessary component to drive an otherwise stable surface system to a patterned state (e.g., Eq. (45) and the following discussion), though we leave systematic analysis of such instabilities for future work. This route to instability does not contradict the preceding suggestions about reducing η and h, as it is likely inadmissible for bacterial pattern formation on agar. Such an experimental setting entails that it is reasonable to assume D B ∝ I, and hence by (45) we see that bulk diffusion will not drive an instability in this case.
We remark that the mode selection phenomena we have illustrated (e.g., in Fig. 5) can be understood in the context of finite-size effects, which are well-studied in the classical case (Murray 2003). Namely, given a dispersion relation for R(λ C (k q )), where λ C is the growth rate of a classical Turing mode, one can tune the geometry to select different spatial eigenvalues k q to give rise to non-monotonic effects as, for instance, the domain size is increased. However, here the effects are more subtle as we cannot explicitly compute the relationship between λ and eigenvalues of the full spatial operator, and so can only implicitly observe these effects. Nevertheless, these mode selection effects appear to be more prevalent compared to classical cases as they require very small domains and other fine-tuning (Murray 2003). Additionally, in our setting mode selection effects appear to be more prevalent across a wide range of geometric parameters, whereas the classical cases have been studied almost entirely in terms of a scalar length, and are generally restricted in parameter regimes where they occur. In particular, we conjecture that the non-monotonic dependence of max(R(λ)) on η seen in Fig. 3 is due to these effects, as we see different modes being excited on either side of this region in Fig. 5.
There are numerous extensions of these results that are worth pursuing. In the example setting of bacterial colony formation on an agar substrate, one might need to augment the bulk evolution with a degradation reaction. We remark that such a simple addition leads to substantial complexity as, if the surface equilibrium is nonzero, then there does not exist a homogeneous equilibrium across the whole coupled system (a degradation reaction in the bulk by itself will always lead to a homogeneous zero equilibrium concentration). We anticipate that the mathematical structure in this case will be even more intricate. A simpler addition, also of relevance to bacterial patterning on agar, would be the inclusion of non-diffusible reactants in the surface region. This approach would also pave the way to account for all gene regulatory dynamics in a quantitative model based on mass action kinetics. In such a case, we can apply techniques to incorporate the impact of such reactants on the surface reaction kinetics directly (in the linearized system) (Klika et al. 2012). Along similar lines, more complicated transport functions g across the membrane can be studied, again leading to new possibilities of differential transport, which can easily be added to the analysis implemented here. We have also assumed that the same number of species diffuse throughout both domains, but in principle one can generalize this by introducing different coupling functions g for the surface and bulk boundary conditions, presently given in Eq. (12). Such an analysis is broadly similar though there are several key details to account for, so we leave this for further work.
There are many biological examples of physical layered media with reactions in multiple different spatial domains, such as in the epithelial-mesenchymal coupling during the development of the skin in mammals (Vilaca et al. 2019). For example, in the study of hair follicle morphogenesis, a substantial amount of biochemical research has implicated Turing-type instabilities in the formation of follicle primordia (Mou et al. 2006). More recently, it has been suggested that a simple activator-inhibitor system is insufficient to capture the dynamical complexity in hair follicle patterning, and so suggestions have been made that such patterns arise due to many coupled processes, which will undoubtedly occur across the different domains of the epithelium and the developing mesenchyme (Glover et al. 2017). Similar remarks can be made about many kinds of skin and other organ patterning events across a range of species, suggesting that general methodologies for stratified reaction-diffusion systems would be useful to elucidate underlying physicochemical mechanisms. A related layered system is the synthetic pattern formation studied in a monolayer of HEK293 cells grown beneath a culture medium (Sekine et al. 2018), where presumably bulk diffusion plays a significant role in transporting signalling molecules.
While we have explored exemplar reaction-diffusion systems in such coupled domains, there are more general transport mechanisms that could be studied. Both chemotaxis and a range of mechanical taxis, as well as mechanical forces, could be included in such a model. We note that a numerical study (Vilaca et al. 2019) has made some progress towards such a model. The linear stability analysis for such problems is involved, but the approach presented here generalizes to these settings. Of course, in the absence of a homogeneous steady state, one must develop new methods for the analysis of pattern-forming instabilities. This has been done recently for heterogeneous steady states (Krause et al. 2020), but extending such an analysis to these coupled geometries is nontrivial. Mathematically, the limit of η → ∞ can be thought of as a step function heterogeneity, as explored in Kozák et al. (2019), so that the systems studied here are also in some sense a generalization of piecewiseconstant reaction-diffusion problems, providing another perspective on heterogeneous reaction-diffusion systems.
Another related generalization would be to study discrete or hybrid discretecontinuum formulations of these kinds of layered media, such as the recent hybrid Turing-type model proposed in Macfarlane et al. (2020). Turing's original paper contained a study of discrete cells (Turing 1952), which was later extended in Othmer and Scriven (1971) and more recently in Nakao and Mikhailov (2010) to reaction-diffusion systems on discrete networks. Such a formulation has been extended to consider multiplex networks, themselves a model of discrete layered media (Gomez et al. 2013), within which Turing pattern formation has also been studied Kouvaris et al. 2015). Such systems deserve exploration on their own, in addition to relating them to spatially continuous analogues of Turing systems in stratified media.
Finally, we mention that one could generalize from our setting of two planar domains to many more coupled domains, or to more complicated geometric settings, including those relevant for more realistic models of development, such as in the blastula stage or later stages of epithelial-mesenchymal development on complicated morphologies. While our approach may be generalizable to very different geometric settings, the dispersion relation we have found in this simple case is already somewhat difficult to analyze, and full numerical simulations may be more expedient. Nevertheless, analytically tractable results for this family of problems are valuable in understanding the role of coupled domain structures in pattern formation, as such scenarios are ubiquitous in biological settings. coefficients in each region for classical Turing patterning systems. In particular, (i) provides a useful consistency check of the modeling framework, while (ii) confirms that in the absence of differential transport within at least one layer or between them, patterning cannot occur for classical Turing systems, in direct analogy to the behavior of the single layer classical Turing instability.

A.1 Patterning in the Limit of an Isolated Surface System, via Sufficiently Small Non-dimensionalized Permeability
We consider the patterning conditions for the isolated surface system in the case that η → 0. Here, we denote the eigenvalues of M B as μ Bp , the eigenvalues of M S as μ Sp , and the eigenvalues of J S as ν p , where p = 1, 2, . . . , n in all cases, where n is again the number of reactants. By continuity of the determinant, Eq. (28) reduces for sufficiently small η to the condition so that the bulk and surface components decouple. As these hyperbolic trigonometric functions are at worst meromorphic, the zero of the determinants occur for eigenvalues of the matrices M B , M S that, respectively, are roots of z tanh(hz) = 0, z sinh(εz) = 0 for z ∈ C. Hence, Eq. (47) is satisfied whenever μ Bp = jiπ/h, or μ Sp = jiπ/ε, with j ∈ {0, 1, 2, . . .}, natural. For the bulk component, this implies that det(M 2 B + ( jπ/h) 2 I n ) = 0 which, as this is a diagonal matrix, implies that the allowed growth rates are given by λ p = −D Bp (k 2 q + ( jπ/h) 2 ) = −D Bp π 2 (q 2 + ( j/h) 2 ) ≤ 0, for q, j ∈ {0, 1, . . .}, p ∈ {1, 2 . . . , n}, with D Bp denoting the pth component of the diffusion matrix D B , and recalling that k q = qπ/L. Hence these solutions do not drive an instability, noting that although these eigenvalues formally break the assumption we made that R(λ) > 0 used to deduce Eq. (28), we have Eq. (47) is precisely (26) in the limit η → 0, and, hence, we have not used this assumption to determine these eigenvalues. Similarly, by the preceding discussion of the eigenvalues μ Sp of M S , we have a condition for the existence of nontrivial eigenvectors for these eigenvalues given by det(M 2 S − μ 2 Sp I n ) = det((q 2 + j/ε)π 2 I n + D −1 S (λI n − J S )) = 0, which we recognize as exactly the condition one would find to compute λ in a traditional n species reaction-diffusion system posed on a rectangle with Neumann boundary data (Klika et al. 2012). In practice finding all such solutions for a given eigenpair q, j is a straightforward numerical problem. While these computations are more easily implemented by noting that the η = 0 limit does separate into two uncoupled regions that can be analyzed via standard methods, these results serve as a useful consistency check of Eq. (28).

A.2 Identical Diffusion Coefficients Within Regions
Below we assume that the surface Jacobian, J S , can be diagonalized, noting that diagonalizable matrices are a dense subset of complex-valued matrices and thus the results derived below therefore hold in general due to continuity. In particular, we now show that if there is identical diffusion for different species within each region, then the spatially homogeneous steady state is linearly stable to the perturbations given by Eq. (16) for kinetics in the surface layer that allow the classical Turing instability.
To proceed, we let D S = c S I n and D B = c B I n , so that M 2 B = (k 2 q + λ/c B )I n . As J S can be diagonalized, we can also diagonalize M 2 S = (k 2 q + λ/c S )I n − (1/c S ) J S . Making this additional assumption, we rewrite these matrices using the n + 1 scalar values m 2 B = k 2 q +λ/c B and m 2 p = k 2 q +λ/c S −ν p /c S , where ν p are the eigenvalues of J s and p = 1, 2, . . . , n. We note that these scalars will never be zero as k 2 q is real and non-negative, we require R(λ) > 0 and, restricting ourselves to kinetics that exhibit the classical Turing instability, we have stable kinetics in the absence of diffusion, and thus R(ν p ) < 0.
With the assumption of identical diffusion coefficients, all of the matrices in Eq. (28) are diagonal and proportional to the identity except for M 2 S , which can be diagonalized using the eigenvectors of J S . Noting that all functions involving M S in (28) are both even functions and analytic, we can simultaneously diagonalize the entire condition using these eigenvectors. Doing so, we arrive at the following n scalar conditions from Eq. (28), c B m B tanh(hm B ) cosh(εm p ) + c S η m p sinh(εm p ) + c S m p sinh(εm p ) = 0, (50) where we note that only one of these conditions must be satisfied, and hence they lead to independent roots for λ. We observe that cosh(εm p ) = 0 cannot occur once m p does not have zero real part, which is enforced by assumptions. Hence we can divide Eq. (50) by this factor to find, c S c B m B m p η tanh(hm B ) tanh(εm p ) + c S m p tanh(εm p ) + c B m B tanh(hm B ) = 0.
The first inequality holds as the real function x sinh(x) is even and, for x > 0, monotonic increasing, while the second inequality follows using the fact sinh(|y|) ≥ |y| ≥ sin(|y|) for all real y and the odd parity of y, sinh(y), sin(y). The case for x < 0 is then inherited from the case x > 0 by parity. Hence, the left hand side of Eq. (52) must have a strictly positive real part, and so this equation can never be satisfied. Therefore, in the case of identical diffusion coefficients within each region, there are no values of λ with positive real part.

B Pure Diffusion Does Not Induce Patterning
We, now, show that diffusion alone, in the absence of reaction terms, cannot induce patterning in the modeling framework. In particular, in the absence of reaction kinetics, and given our linear interfacial conditions, each species decouples and can be considered in isolation. Without loss of generality we consider surface and bulk concentrations of the first species, denoted c s and c B below with respective diffusion coefficients d B , d S in each region. Then, with B and S denoting the bulk and surface regions, as in Fig. 1, and the subscript t denoting a time derivative, we have where the surface integrals in the fourth line vanish, courtesy of the zero flux boundary conditions (10) and (11), and the interfacial conditions (12), on noting relations (6) and (7). This can also be recognized as a free energy inequality, or equivalently an entropy inequality, corresponding to the second law of thermodynamics given a Fickian diffusive flux as a constitutive relation (Gurtin et al. 2013). Hence a standard measure of heterogeneity cannot increase and thus initial conditions that are close to homogeneous (in the sense of a suitable Sobolov norm) cannot induce patterns.