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 bulk–surface membrane–cytosol interactions (Halatek et al. 2018; Rätz and Röger 2014; Kretschmer and Schwille 2016; Spill et al. 2016; Cusseddu et al. 2018) to epithelial–mesenchymal 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 (Budrene and 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; Halatek and Frey 2018), 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. 2016, 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. 2016, 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 reaction–diffusion 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 m s. 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.

Fig. 1

Example experimental system under consideration. a Here we consider cells growing in culture on top of a solid reservoir of nutrients, such as agar. b The surface (cellular) region denoted has height and contains both reaction and diffusion terms, whereas the bulk (nutrient) region 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)

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
L^	8×10^-3 m-2.2×10^-2 m	See text
‖D^_S‖_∞,‖D^_B‖_∞	10^-10 m^2 s^-1-10^-9 m^2s^-1	Table S8, Grant et al. (2016)
1/τ	8.4×10^-5 s^-1-1.0×10^-3 s^-1	Tables S8, S9, Grant et al. (2016)
H_ε	10^-5 m-10^-4 m	See text
H	10^-3m -10^-2 m	See text
η^	Unknown	−

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; Halatek 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; Halatek et al. 2018; 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 (Halatek and Frey 2018). 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 (Asllani et al. 2014; 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; Halatek and Frey 2018). 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. 2003, 2005; Haim 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; Gorder et al. 2019) 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.

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 where we refer to as the bulk region, and as the surface region (see Fig. 1). We write for the concentrations of reactants in the bulk, and 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:where 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 haveWe non-dimensionalize the above model via concentration, time and length scales corresponding to the reaction kinetics and a unit lengthscale , respectively. Specifically, we define , , where is a diagonal matrix of concentration scales. Equally, we set , where is the timescale of the fastest reaction in the surface and bulk, and . The dimensional scalings are then chosen such thatWe define new dimensionless groupings , , , , and . The nondimensional system is written asThere 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 , one can consider either scaling or scaling 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 88 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) and Fig 3E Grant et al. (2016)). Thus we consider a range of 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 ms ms 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 . Similarly, in the parameter fitting by Grant et al., a reaction timescale on the scale of the faster reactions is such that s s, 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 m m (Levin and Angert 2015), and hence a surface depth on the scale of m 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 mm; estimates for the interfacial permeability, , 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.

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, , 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 as well as the range of are presented as they will be important in the asymptotic analyses below

Parameter	Typical value/range
ε=H_ε/L^	4.5×10^-4-1.3×10^-2
h=H/L^	0.045-1.3
L~=L/L^	1
ε_∗=ε‖D_S^-1J‖∞1/2	1.1×10^-3-0.87
h_∗=h‖D_S^-1J‖∞1/2	0.11-87
ε∗2/3	4.0×10^-7-0.3
‖D_S‖_∞=τ‖D^_S‖_∞/L^^2	2×10^-4-0.2
‖D_B‖_∞∼‖D_S‖_∞	2×10^-4-0.2
η=τη^/L^	Unknown

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 specifyingso that the surface reactions determine the spatially homogeneous steady state concentration in both regions, and our simple constitutive choice of 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 formwhere , 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 satisfywhere the Jacobian, , 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 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 for the zero mode of the transport operator acting in the surface layer, , one hasfor a general 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 requiresOne possible solution occurs for , which generates a spatially homogeneous mode that has already been considered above. Satisfying this equation for other requires mathematical fine-tuning as 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 iswhere 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 (). Furthermore, we note that even in the absence of completeness, 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 we will have possibly different eigenfunctions and . 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 . As we are looking for modes which grow in time, leading to instability, we will impose 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 findwhere denotes the ordinary derivative with respect to y. These spatial functions are required to satisfy the external boundary conditions and the coupling conditions which read,To find suitable that solve the coupled problem (19)–(21), we will make use of the matrix-valued function defined by , for some matrix , as well as . We recall the differentiation identity , 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 defineWe 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 and , and thus functions of and . This dependence on the squares of these matrices follows as if 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, and 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 and are the square roots with positive (or possibly zero) real parts of the eigenvalues of and .

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 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 is diagonal, and from our assumption that , we have that its eigenvalues have positive real part. Therefore, the elements of are given by the hyperbolic cosine of the diagonal elements of , and since these are all positive definite, is invertible.

Now we define the matrices , , , . By the above argument, we have that is invertible. We then have that (26) can be written (by exchanging rows and using the Schur complement) as,Noting that , we have thatso Eq. (27) is equivalent to,We note that the Turing instability conditions for the surface system in isolation—neglecting spatial structure in y—are precisely that the growth rates computed from have negative real part for , and positive real part for some , 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 (), the surface domain is a rectangle, and hence, the surface eigenfunctions are also cosines of different spatial eigenvalues, which can vary independently from . 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 . 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, , has elements that are of the same order and thus of the order of , so that is of the same scale as for any matrix considered.

Table 3

Thin-surface limits obtained in different asymptotic regimes given  ord(). Note that moving left to right corresponds to an increasing size of , and moving top to bottom corresponds to increasing scales of . No instabilities: ; Isolated (1-D) surface: ; Quadratic : , ; Reduced instability: , ; Averaged condition: Eq. (41)

	Case I.	Case II.	Case III.
	(h_∗≪ε_∗≪1)	(ε_∗≪1, ε_∗/h_∗=ε/h∼ord(1))	(ε_∗≪h_∗)
η≪ε‖J_S‖_∞	Isolated	Isolated	Isolated
η∼ord(ε‖J_S‖_∞)	Isolated	Quadratic condition, Eq. (40)	Reduced instability
η≫ε‖J_S‖_∞	Isolated	Averaged condition, Eq. (41)	No instabilities

Thin-Surface Limits

Here, we consider an asymptotically thin surface, requiring . First note that in the thin-layer limit below, the surface Jacobian only appears viaIn addition, given patterning (i.e., ), the matrix cannot be dominated by the terms or within 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 patterningand also that has an upper bound ( and in particular the term is in fact bounded). Noting the boundedness of we have that is invertible, as for sufficiently small this matrix has a determinant which is asymptotically . In addition, for sufficiently small , we have the Taylor expansionwhere means the same scale as , or smaller, as the surface thickness tends to zero. Thus by right multiplying (26) by and Taylor expanding, we obtain (to leading order) the relationproviding (ensuring the invertibility of and the validity of the Taylor expansion above). Furthermore, noting that together with the relations (30), which give the maximum scale of and show that we haveusing The latter inequality is immediate in the two species case on writing with as then with a trivial generalization to higher number of species. Hence, for conditions associated with patterning, the relative error in the leading-order thin-surface approximation arising from Eq. (31) is , and thus, we require 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 at least for the parameter estimates of Tables 1 and 2, and also relations (30), (33), we also haveHence, an appropriate scale for the largest component of is , which ranges from small to large in Table 2 and thus we proceed to consider simplifications of the expression within the instability condition (28) for small and large values of . For the small limit, a Taylor series expansion immediately gives , with relative corrections on the scale of , and we also haveFor large simplifications, first note that is diagonal, with diagonal components that have positive real parts since 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 and enforce by wavemode selection, which bounds the real part of away from zero.

For with , we have the limitas 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 by the diagonal matrix with positive semi-definite real part for the square root of the diagonals of , and in fact no such square root has zero real part since the diagonals of have positive real part. Consequently, at leading order we have in the large limit that and thus , withFinally, with this definition of , 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 for . This may be summarized together with Eqs. (35) and (36) viaWe 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 , which is a measure of the non-dimensional depth of the bulk relative to the patterning lengthscale. We also consider the case for relative completeness, even though Tables 1 and 2 highlight that is anticipated for experiments with synthetic pattern formation within bacterial populations.

Thin-Surface Asymptotic Regimes with

An example of patterning when is given in the next subsection, but here we consider thin-surface asymptotics with on fixing  ord where ord is defined to mean both and not 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 towhere the scale of the non-dimensional permeability is unknown, and the possible values of are wide-ranging. 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 , in order to understand the impact of the bulk on the classical single-domain situation.

Case I : Noting that is equivalent to by definition, in this limit Eq. (38) reduces to,However the determinant with the hyperbolic tangent term cannot generate a root with and thus patterning. In particular, in Appendix A.2, following Eq. (53), it is shown that when one also has With for patterning, let , where is a bulk diffusion coefficient. Thus is an eigenvalue of , and all eigenvalues of this matrix are of this form. Furthermore, we have where z is an eigenvalue of and satisfies However, for the hyperbolic tangent term in Eq. (39) to generate a root, at least one eigenvalue of , that is one such z, must satisfy a contradiction, thus showing there are no roots from the determinant involving the hyperbolic tangent. Hence, noting is positive definite, the only roots are those of the isolated Turing modes, independent of , and determined purely from .

Case II , : This limit corresponds to the entire domain being thin with respect to the lengthscale in the x direction (). In this case, we have,Equation (40) is a slight modification of the isolated surface Turing conditions in 1D given by , and can similarly be written as an nth-order polynomial in . Further, if , then the conditions for instability are precisely those for an isolated surface. Similarly, if , 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 , 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 , 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 , for a natural number q.

Another plausibly relevant case of Eq. (41) is if , i.e., the surface and bulk diffusivities are the same. Here, the dispersion relation is that of the classic case except and are both scaled by . Hence the allowed values of are those of the classic case divided by , 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.

Case III : This case proceeds similarly regardless of whether or . Noting , from Eqs. (29) and (30), by square rooting the first of relations (34), and Eq. (37), that is , we haveNoting as deduced just below Eq. (33), and we thus haveHence the final term from relation (32), that is , may always be dropped relative to the first term, that is . This reveals that the instability condition simplifies toThe hyperbolic tangent does not contribute to instability, by analogous reasoning to Case I. Thus, there is no instability unless is concomitantly small alongside . Writing out , 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 is increased, we observe a trend of moving from the isolated surface system to a reduced, or average system, and eventually, for and , 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.

Further Thin-Surface Asymptotic Regimes with

There exist nontrivial asymptotic limits which are not described by  ord, which can lead to instabilities not captured in Table 3, as we now show. In particular, with and we then have and finally . Hence from Eq. (32), after expanding for small z, we have at leading orderwith relative corrections of which is required to be much less than unity. Further, noting we have already assumed , we thus additionally require for Eq. (44) to hold. For the range of parameters detailed in Table 2, we have , and thus we have Eq. (44) is typically valid for parameters associated with synthetic patterning in bacterial colonies, but not always.

We proceed by noting that the first and third terms of Eq. (44) are and the second is . WritingWe then have 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 , 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 .

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 , 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 given bywith , . The spatially homogeneous steady state is given by . Unless otherwise stated, we will assume equal diffusion coefficients between the surface and the bulk given by the diagonal matrices . 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 , 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 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 exists, and tested this against full 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 , 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 as a continuous function of the spectral parameter in the x direction, , as is commonly done (Murray 2003). For , 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 (which can be computed in the standard way (Murray 2003)). We can then consider the maximum value of (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 size such that a particular mode 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 and each influence the shape of these dispersion curves greatly. We now compare these curves to the predictions in Table 3. For Case I (), we see that is almost unchanged to the standard case up to small corrections not captured by the asymptotics. In Case II (), 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 () is consistent with the asymptotics of Table 3 whenever except for and at relatively small values of . 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 with ord. As a consistency test of this suggested mechanism, in Fig. 2d we replace 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 and sufficiently large, fall below the axis as expected. This is true for different scalar multiples of the identity, such as 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 .

Fig. 2

Dispersion relations in the x coordinate computed via (28) for a continuous variable , using the parameters , , with surface diffusion parameters , and . In a–c we take , though in d we set , corresponding to equal bulk diffusion between species. The solid lines correspond to for different values of for the bulk–surface condition, whereas the dashed line corresponds to the single-domain classical case. For c, we anticipate there is an instability for relatively low and large due to surface–bulk interaction instabilities, as exemplified in Sect. 4.4 for and (Color figure online)

As in the classical case, we expect that for sufficiently large domains, any region where should admit a patterned state. We confirmed this using for each of the dispersion curves, finding that they admitted patterned solutions for long time simulations if and only if for some region in -space. Similar to the classical case, the layered model is always observed to stable at 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 aswhere is simply a positive definite (dimensional) scaling parameter. Note that and for , if and only if is spatially homogeneous. While we do not anticipate this metric to be quantitatively comparable to , 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 near the onset of instability. The value c is taken so that when 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 and when is large, though the functional will still indicate whether or not predicts pattern formation, as well as the scaling of pattern heterogeneity as a function of near the onset of instability.

To use this heterogeneity functional, the full system (8)–(12) was solved until a final time of to ensure a good representation of the steady state pattern. The initial data were taken to be and with and random fields such that at each value of (x, y), they are independently and identically distributed normal random variables with zero mean and variance . The equations were simulated using the COMSOL Multiphysics software (http://www.comsol.com) with at least second-order triangular finite elements. A non-uniform mesh was constructed such that the surface region 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 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 . Additionally, we see a jump in the value of the heterogeneity between and . We plot values of in Fig. 4 across this jump to demonstrate that this discontinuity in the value of the spatial heterogeneity 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.

Fig. 3

Non-trivial dependency of Turing instabilities on geometric parameters. Plots of given by Eq. (28) in blue computed across 250 values of for different parameter combinations, and plots of given by (46) in red asterisks for 100 values of . The other parameters were taken as , , , , . The constant c in was fixed per set of parameters/panel to match the maxima of and across to qualitatively compare these metrics. The parameter sets corresponding to and gave qualitatively the same results as in panel (c) with for all (Color figure online)

Fig. 4

One-dimensional plots of corresponding to parameters in Fig. 3a for two values of in the top two panels, and plots of the corresponding below (with in all cases). The surface concentration 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 and ) (Color figure online)

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 occurs approximately for within the range . We show examples of the mode selection process from Fig. 3d in Fig. 5. For small , we see stable multiple-spike solutions that are essentially confined to the surface. As increases further to , 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 . 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 ) 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).

Fig. 5

Plots of and corresponding to parameters in Fig. 3d for three values of , and . Here, and (Color figure online)

In all of Fig. 3 we observe that appears to asymptotically approach a fixed value for either (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 change in the number of unstable modes was apparent for variation in , though the structure of the solution does change.

Fig. 6

Plots of and corresponding to parameters in Fig. 3 with and for three values of , and (Color figure online)

Fig. 7

Plots of and corresponding to parameters in Fig. 3 except that and for three values of , and (Color figure online)

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 (), 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 () 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 . 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 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 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 non-dimensional parameter groupings, and . 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 . Noting that is chosen such that , one can deduce more generally that is the dimensional timescale of surface reaction. Hence, the final parameter grouping is the ratio of the interface permeability to the surface velocity scale, , 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 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 surface–bulk 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 () that these do not enlarge the Turing space. Note that in systems where diffusion varies significantly between domains (e.g., non-diffusible proteins in the surface) the parameter space that admits pattern formation can increase with increasing bulk size (see, for instance, Halatek and Frey (2018); 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 , 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 , where is the growth rate of a classical Turing mode, one can tune the geometry to select different spatial eigenvalues 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 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 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 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 piecewise-constant reaction–diffusion problems, providing another perspective on heterogeneous reaction–diffusion systems.

Another related generalization would be to study discrete or hybrid discrete-continuum 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 (Asllani et al. 2014; 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.