Theory of mechano-chemical patterning in biphasic biological tissues

The formation of self-organized patterns is key to the morphogenesis of multicellular organisms, although a comprehensive theory of biological pattern formation is still lacking. Here, we propose a minimal model combining tissue mechanics to morphogen turnover and transport in order to explore new routes to patterning. Our active description couples morphogen reaction-diffusion, which impact on cell differentiation and tissue mechanics, to a two-phase poroelastic rheology, where one tissue phase consists of a poroelastic cell network and the other of a permeating extracellular fluid, which provides a feedback by actively transporting morphogens. While this model encompasses previous theories approximating tissues to inert monophasic media, such as Turing's reaction-diffusion model, it overcomes some of their key limitations permitting pattern formation via any two-species biochemical kinetics thanks to mechanically induced cross-diffusion flows. Moreover, we describe a qualitatively different advection-driven Keller-Segel instability which allows for the formation of patterns with a single morphogen, and whose fundamental mode pattern robustly scales with tissue size. We discuss the potential relevance of these findings for tissue morphogenesis.


A. Derivation of the model.
A.1. Tissue spatial organisation. As sketched in Fig. 1A of main text, we model multicellular tissues as biphasic porous media, with a first phase consisting in a three-dimensional elastic meshwork made of interconnected cells, and a second phase being made up of aqueous extracellular fluid permeating in-between cells. The typical lengthscale of the whole tissue is l. Cells, compartmentalized by their plasma membrane have a typical diameter lc and the average intercellular or interstitial space has a characteristic size li. Typically, these three lengthscales are well-separated such that: We adopt a continuum description of the tissue where all physical quantities are averaged over a representative volume element (RVE) of intermediate size l lr lc at which the system can be coarse-grained (1,2). At that scale, we introduce the space coordinate r ∈ Ω, where Ω is the domain occupied by the whole tissue, and t denotes the time coordinate.
A.2. Derivation of the biphasic model. Although we only consider a biphasic media in the main text of this manuscript, a more complete version of the model would introduce three volume fractions within the tissue: φs( r, t) the volume occupied by the skeletal meshwork of cells divided by the volume of the RVE, φc( r, t) the volume of fluid contained within cell membranes divided by the volume of the RVE and φe( r, t) the volume occupied by interstitial fluid divided again by the RVE. In the following, we provide a rationale and orders of magnitude for considering only two phases in the main text.
A difference with the simple biphasic model is that there are here three velocities associated with each of the phases, vs, vc and ve. Importantly, although it could be thought that fluid is only exchanged between cells and extracellular medium (meaning that vc = vs i.e solid phase and intracellular fluid are advected together), this is not the case, as gap junctions allow for direct cell-cell communication and fluid exchange (which means that in a triphasic model there can be non-zero vc with vanishing vs). This point is important to consider when writing intracellular morphogen conservation equations, as vs would be included in an advective term, whereas vc would not, as most morphogens are too large to be transported through gap junctions. In order to simplify the description of the system, we make two important assumptions relative to the cellular meshwork phase. First, as both cells and extracellular space are mostly composed of fluid, we neglect the volume fraction occupied by the cellular meshwork such that φc def = φ and φe 1 − φ. Then, we assume that the velocity of cells moving in the tissue is very slow compared to the velocity of fluid flows in the extracellular space and from cell to cell. Thus, the present study focuses on cell specification patterns considering that the evolution of the tissue topology is quasi-static. To make this approximation clearer, we provide below (Section A.3) simple physical estimates showing that the characteristic timescales associated with tissue remodelling processes, such as cell proliferation and active cell motility, can be larger by an order of magnitude than fluid circulation inside the tissue.
In Section I, we also show how the theory can be extended to account for large scale tissue remodelling occurring during embryonic development when tissue topology is changing on a timescale similar to that of fluid flows inside the tissue.
A.3. Fluid re-circulation vs tissue remodelling. Large scale advective flows of extracellular fluid must be compensated by some recirculation of fluid ensuring global volume conservation (except if there are exchanges with other tissues, which we do not consider here for the sake of simplicity). We first consider the possibility that this re-circulation occurs through intracellular fluid flows, and calculate the time scales associated with this process. In general, there are three characteristic timescales associated with fluid circulation inside a multicellular tissue: τi for fluid flows inside the interstitial space, τc for fluid flows from one cell to the other through gap junctions, and finally τ for the exchanges of fluid between cells and the extracellular fluid through aquaporins.
As we demonstrate in Section A.7 using a poroelastic model of multicellular tissues, τi and τc can be interpreted as diffusive timescales and take the form, τi = ηl 2 p /(κK) and τc = ηl 2 p /(κcK) where η ≈ 10 −3 Pa.s is the viscosity of the aqueous cytosolic fluid (3,4), K ≈ 10 4 Pa is the drained elastic modulus of the tissue (5), lp ≈ 10 −4 m the characteristic size of multicellular patterns in the tissue. κ and κc (expressed in [m 2 ]) are the effective interstitial and cell to cell water permeabilities. In agreement with the typical timescales involved in cell volume changes following an osmotic perturbation (6), we estimate that τ ≈ 10 2 − 10 3 s. In the main text, based on typical tissue interstitial sizes and packing fractions, we estimate the interstitial permeability κ ≈ 10 −15 − 10 −19 m 2 leading to τi ≈ 1 − 10 4 s. We note that this estimate is one of the coarsest one in the model, as the packing fraction of multicellular tissues can vary widely between tissues, developmental stages and species.
Finally, the cell-cell permeability relies on gap junctions which are highly specialized membrane structures made of clusters of transmembrane protein channels named connexons and which connect together the cytoplasms of neighbouring cells (7). The ensuing permeability can be estimated in applying the Darcy-Weisbach equation to a parallel array of connexons mediating the cell to cell fluid flux: taking, lc ≈ 10 −5 m the diameter of a cell, n ≈ 10 5 as the average number of connexons at a cell-cell interface (8)(9)(10)(11), d ≈ 3 × 10 −9 m (7) as the typical diameter of a connexon and δ ≈ 10 −8 m as its length (7,12). The coefficient 32 is the friction factor for a laminar flow (13). Interestingly, this simple theoretical estimate of κc is in line with direct experimental measurements using fluorescent dyes injected in the cytoplasm of cells in epithelial cell monolayers (14,15). Eventually, this yields to an estimate of τc ≈ 4 × 10 2 s. Next, we compare these timescales for direct cell-cell fluid re-circulation to the timescale associated with tissue remodelling. Indeed, an alternative to enforce global volume conservation given extracellular fluid flows would be to have a counter-flow of the celular meshwork. Typically, this would be very slow, as the viscosity of multicellular tissues is very high. Such a timescale has been estimated from rheological measurements made on various multicellular tissues (16,17) and is typically τr ≈ 10 4 s, associated with slow tissue re-arrangement and remodelling. This slow timescale is consistent with the observed timescales for cell division and cell motility, which have been both noted as major sources of tissue remodelling (18). The typical timescale of cell division in embryonic tissues is τ d ≈ 10 4 − 10 5 s. The characteristic velocity of a cell moving within such a tissue is vm ≈ 10 −9 m.s −1 . Over the characteristic size lp ≈ 10 −4 m of tissue patterning, this leads to a typical timescale for motility-driven remodelling of τm ≈ 10 5 s.
In the end, it appears that all three timescales associated with fluid flows are either slightly or much smaller compared to the timescale associated with tissue remodelling and cellular flows. Therefore, in our model, we shall consider tissue remodelling as a quasi-static process and assume a fixed tissue structure during the faster fluid exchanges in-between cells or between cells and extracellular fluid. A more complete description incorporating these effects and wich does not make such a simplifying assumption is presented in Section I. This is motivated by the fact that in some biological systems, during early embryogenesis, cell division could be faster than our estimate, closer to τ d ≈ 10 3 s. In this situation, plastic re-arrangements of the tissue could be potentially more important and play a role in the patterning instabilities we describe.
A.4. Fluid conservation. As fluid can be exchanged between cells and between cells and extracellular medium, we can write, assuming fluid incompressibility, the following mass balance relations: where vc,e( r, t) are velocities associated with fluid flows respectively inside and outside the cells. It should be noted that because these equations are defined at a coarse-grained level, vc refers to fluid flowing directly from one cell to the other, which occurs in multicellular tissues via gap junctions, as discussed in the previous section. Gap junctions have been well-characterised as key regulators of cell-cell coupling, allowing in addition of fluid transport, the movement of inorganic ions and small water-soluble molecules from the cytoplasm of one cell to the other (19). Nevertheless, most morphogens cannot cross through gap junctions, as they are proteins whose size exceeds by far the diameter of connexons lumen (7,12). Crucially, in our model, this means that although extracellular flows ve can advect morphogens, this is not the case of intercellular flows vc. This allows advection to concentrate morphogens locally, as they do not participate in fluid re-circulation. Summing up the two previous mass balance equations, we obtain ∇.(φvc + (1 − φ)ve) = 0 which is generally satisfied if: As φ is usually rather close to one, the cell to cell fluid velocity vc is thus expected to be much smaller than the velocity of fluid flows developing in the extracellular space in-between cells ve. However, in spite of taking relatively low values, vc should not be neglected, as it is central for the model consistency and in enforcing volume conservation in the system. We note here that although the cell volume fraction φ is close to one in various multicellular tissues, our theory does not depend on that assumption, and can in principle accommodate values of φ very different from one. The source term S, on the right-hand side of Eq. 1, originates both in the genetically controlled transport of water between cells and extracellular medium, which occurs through specific transmembrane proteins named aquaporins (20,21) and on the balance between cell division and apoptosis regulating the number of cells in a RVE. In absence of spatial fluid flows in interstitial spaces or at tissue boundaries, we have: where Nc is the number of cells in a RVE, Vc the volume of a single cell and Vr the volume of an RVE. When submitted to an external osmotic shock or mechanical load, the volume of cells transiently increases or decreases prior to relaxing to an homeostatic volume V h . Considering an exponential relaxation, we can write: where τ is the typical timescale associated with a regulatory response to cell volume changes. This time scale can be estimated to be of the order of ∼ 10 2 − 10 3 s in the case of an osmotic shock (6). In our model, as discussed with the orders of magnitude provided in Section A.3, the tissue topology is taken to be fixed and we thus neglect the variations of Nc due to cell division and cell loss. However, we present below a more general argument and show that these phenomena does not impact the expression of the term S in Eq. (1). This argument will be also used to justify the general model presented in Section I where tissue topology becomes dynamic thanks to cell motility and cell turnover.
The simplest evolution law for Nc is given by ∂tNc = s d − Ncka, where s d accounts for cell division and ka is the cell loss rate. This law can be rewritten as is the typical duration of a cell cycle of the order of ∼ 10 4 − 10 5 s (during chick limb bud patterning this time is about 10 h = 3.6 × 10 4 s) and N h = s d τ d is the equilibrium number of cells when cell division balances cell loss. Defining φ h = N h V h /Vr the homeostatic volume fraction of cells and δV and δN some initial perturbations from the homeostatic state, we obtain: Considering typical timescales such that t τ and as τ d τ , the first term dominates the last expression and we finally obtain In agreement with the fact that cell fate decisions are usually accompanied by modifications of the cell cycle and cell morphology, including cell volume changes, we suppose that φ h depends on the concentrations of morphogens (22)(23)(24). Indeed, most of the known morphogens are in fact growth factors which regulate both cell differentiation but also cell growth and cell division (25)(26)(27)(28). We can thus write the source term of the conservation equation as: Ii) is the cell volume fraction with cells in their homeostatic volume state.
As mentioned above, it should be noted for the sake of completeness that we have assumed here a closed-system, where no fluid is exchanged between the tissue and the rest of the organism. Relaxing this assumption would be equivalent to adding another source term, S b , in the conservation equation for 1 − φ (see Eq. 1), without adding it to the conservation equation for φ. This would be an alternative to intracellular fluid re-circulation to allow for the existence of non-zero steady state fluid flows in the tissue.
A.5. Morphogen conservation outside cells. Conservation equations for morphogens in the extracellular space are given by: where γA,I and λA,I are respectively the import and export rates of morphogens through some specific channels and pumps, and D is the common global diffusion coefficient of both morphogens in the intercellular space. D do not refer to a molecular diffusion coefficient but to an homogenized quantity over a RVE.

A.6. Morphogen conservation inside cells. Conservation equations for morphogens inside cells read:
where f and g are the morphogens turnover rates, non-linear functions of morphogen concentrations describing their production and degradation by cells. Note here, that these reaction terms could also be a function of extracellular morphogen concentrations, describing a biological situation where signaling at the membrane also directs internal morphogen production and degradation.
While such an addition would be biologically relevant, it nevertheless does not alter the mathematical formulation of our final model Eq. (10) as explained in Section C.
The vector field [f, g] is assumed to have a single stable equilibrium point ( det t > 0 and tr t < 0 such that the eigenvalues of the jacobian matrix Here, we do not consider the complex case where morghogens can undergo transitions between various stable chemical states. Contrary to Eq. (2), advective terms are absent in Eq. (3), as known morphogens are too large to be transported through gap junctions (7,12). Indeed, because of this impossibility of direct cell to cell transport of morphogens, advective terms would equate the transport terms ∇(φAivs) and ∇(φIivs), where vs( r, t) is the cellular meshwork velocity, which we neglect here in ignoring cellular motion. Such an assumption is critically discussed in Section I. A.7. Effective mechanical behaviour of the tissue. To complete the model described by Eqs. (1,2,3), we need to relate the interstitial fluid velocity, ve, with the cell volume fraction φ. Such a link is modelled using the classical poroelastic theory (29,30).
When biochemical reactions involving morphogens are at equilibrium, Eq. (1) has an homogeneous solution given by: This state of the tissue constitutes a "reference configuration" where the deformation of the multicellular network is considered to be null. Supposing small variations from this reference state, we can write an effective quadratic mechanical free energy for the cellular network, F , coupling the tensorial infinitesimal strain i with porosity variations δφ = φ − φ * . The free energy of a poroelastic tissue reads: where G is the shear modulus of the tissue, K is its drained modulus, Ku its undrained modulus and α its Biot coefficient.
We derive the non-dissipative mechanical stress in the tissue, σ, and the interstitial fluid pressure p by differentiation of the free energy: where s denotes the identity matrix. In this classical theory, mechanical energy dissipation is assumed to happen mainly through viscous friction of the extra-cellular fluid permeating in between the cells and is described using Darcy's law: as, again, we neglect vs in this framework, but this time as compared to the extracellular fluid velocity ve. The tissue permeability is denoted κ ([m 2 ]) and η is the viscosity of the permeating aqueous fluid. Finally, force balance within the tissue implies that ∇.σ = 0, which using the first constitutive relation Eq. (4) leads to: Coupling that last equation with the second equation of Eq. (4) and assuming, in line with the typical orders of magnitude for biological tissues (31, 32) that Ku K G and α 1, we obtain the simple relation: relating ve and φ where we can introduce the hydrodynamic diffusion coefficient of the extracellular fluid through the tissue Dm = Kκ/η, which is, as explained in the main text, a driving force of the instability, counter-acting morphogen Fickian diffusion which acts as a resisting force.
B. Model of an active biphasic tissue. Coupling Eqs. (1,2,3,5), the full set of equations describing the chemo-mechanical behaviour of an active biphasic multicellular tissue can be written as: This system of coupled equations can also be written in a more compact matrix form: where X i,e = (Ai,e, Ii,e) are the morphogen concentration vectors, F (X i ) = (f (Ai, Ii), g(Ai, Ii)) is the morphogen turnover rates vector, and £ = λA 0 0 λI and = γA 0 0 γI are the diagonal matrices of transmembrane import and export rates of morphogens. The model has as a single homogeneous steady state solution given by:

P. Recho, A. Hallou and E. Hannezo
where u = £ −1 is the matrix of the transmembrane transport equilibrium constants.
To investigate the stability of this solution and fully specify the problem, we need, in addition of Eq. (6), to define initial and boundary conditions. The initial condition is typically chosen as a random perturbation of the homogeneous solution Eq. (7). For the boundary conditions on the frontier ∂Ω of the domain Ω, we impose a no-flux of fluid and morphogens: where n is the outward normal of Ω. Such boundary conditions are natural as they correspond to an isolated system and we wish to consider endogenous pattern formation that are not externally driven.
Eventually, we also need to specify the non-linear function f, g and φ h as a function of the internal morphogens concentrations. For f and g, we will essentially consider two cases in this paper: For Fig. 2 of the main text, we use the classical activator-inhibitor Gierer-Meinhardt scheme (33-35) (see section "Turing-Keller-Segel instabilities" in the main text): where ρ is the rate of activation and inhibition and τA,I the timescales of degradation of A and I. Within this scheme, Fig. 4 of the main text, to illustrate the instability induced by the hydrodynamic cross-diffusion (see section "Crossdiffusion Turing instabilities" in the main text), we construct a very simple inhibitor-inhibitor scheme where: is a non-linear function chosen to ensure that Ai and Ii remain bounded and positive (H denotes the Heaviside function). Parameters used for the figure numerical simulation are βKAτA = 1, Amin/A * i = Imin/I * i = 1 and lA/l = 0.1. Generically, φ h (Ai, Ii) could be a rather complex function. However, for a linear stability analysis of patterning instabilities, we can restrict its functional form to a first-order linear expansion of the cell volume fraction around the reference state: The term φ * accounts for mechanisms controlling cell volume and proliferation that are independent of morphogens levels, such as mechanical stress caused by active cell contractility (36), osmotic effects due to passive water flows trough aquaporins (20,21) and ions efflux/influx by transmembrane ion channels and pumps (6). The χA,I terms account for the sensitivity of homeostatic cell volume fraction to intracellular morphogens concentrations and thus depend on how morphogens regulate both the volume of individual cells but also their proliferation. We also note that although this linear expansion is valid to study the onset of the instability, this does not constrain φ to be bound between 0 and 1, which is required from a physical perspective, and arises from non-linear terms. This becomes important for numerical simulations, where non-linear terms are crucial to stabilize the amplitude of the instability, so we chose for all numerical integrations: , which is compatible with the linear approximation mentioned above. C. Simplified model. To achieve mathematical tractability and discuss the solutions of the model in a biophysically relevent context, it is possible to reduce Eq. (6) by assuming that the import rates £ t are much larger than the export rates, such that u 1 is small and that the import/export dynamics is faster than morphogens turnover. From a biological standpoint (see also Section G.3 for further details) , £ ≈ 10 −2 s −1 t ≈ 10 −4 s −1 , (37) and we explain in Section F why we additionally consider the limit u 1. Note however that we demonstrate in Section G.1 that the Keller-Segel instability that we disclose in the single morphogen case does not depend on the simplifying double limit considered in this Section. When £ t, morphognen import and export almost balance each other such that X e uX i . As a result the external concentrations of morphogens are small compared to the internal ones. However, exchange terms in Eq. (6) remain finite because −1 (X e − uX i ) is an indeterminate product between a large quantity ( −1 ) and a small one (X e − uX i ). To resolve this indeterminate form, we sum up the first and second equations of Eq. (6) such that the exchange terms exactly cancel out each other and, as X e uX i , we obtain a closed system for the internal morphogens concentration: For such a system to be non-trivial, we assume that both Dmu and Du remain finite. This implies that the diffusion Dm is large and for lm = √ Dmτ to remain finite as well, we also need to assume that τ ∼ τAKA 1 is small. As a result, the cell volume fraction relaxes infinitely fast and we finally obtain the simplified model: [10] which couples a reaction-advection-diffusion equation for the morphogens levels with an elliptic equation controlling the drift as a function of the morphogen levels. Interestingly, the mathematical structure of this model is reminiscent of some previous models which have been proposed in the completely different context of cytoskeleton self-organisation (38). At this stage, we note that the robustness of the model with respect to two hypothesis can be demonstrated. In Eq. (6), we can add to the second equation a degradation term for the external morphogen concentrations − extX e and we can also suppose that, in the third equation φ(X i , X e ) also depends on X e . As X e = uX i , this will not change the form of Eq. (10) but only effectively rescale the expressions of F and φ h . Importantly, when χA,I = 0 (no sensitivity of the volume fraction on morphogen concentrations), the second equation uncouples from the first one such that φ ≡ φ * and Eq. (10) reduces to the classical Turing case: which is fully studied in (39) with the only difference that here the global diffusion of morphogens is now rescaled by the transmembrane transport equilibrium constants (see section "Orders of magnitude on morphogen transport" in the main text).
D. Stability analysis of the simplified model. We are now interested in determining if we can obtain a stationary spatial pattern for morphogens concentrations in perturbing the system Eq. (10) in the vicinity of its homogeneous equilibrium state. To do so, we write: 1, and we obtain at the first order in : We now introduce the sets of eigenvalues {λ k } k≥1 and eigenvectors {U k ( r)} k≥1 of the minus Laplace operator defined by the problem: By classical theorems (39), {λ k } is a countable set of positive and increasing values, and {U k } forms a basis on which we can project the solutions of the linearised system : where ω k is the growth rate associated to the k th growth mode. For example, when Ω is a one dimensional segment of size l we have λ k = (πk/l) 2 and U k (x) = cos(πkx/l). From this ansatz, we obtain that: Thus, unstable spatial modes are those for which a ω k with a positive real part solves the dispersion relation: [12] To be in line with the classical Turing's model, we impose that the homogeneous solution in the absence of spatial terms (λ k = 0) is stable. Thus the two solutions ω 0 of det ω 0 (g + φ * s) − t = 0, have negative real parts. Such condition is satisfied as soon as φ * + tr g > 0, det t > 0 and φ * tr t + det ttr (t −1 g) < 0. [13] Otherwise, the homogeneous solution can lose stability even without spatial terms and the steady state would then oscillate in time. These conditions generalize the classical ones (39), det t > 0 and tr t < 0 obtained for χA,I = 0.
In the main text (cf. Fig. 2A), we study the dispersion relation Eq. (12) in the case where t is computed based on the classical activator-inhibitor scheme Eq. (8) and χI = 0 since conditions Eq. (13) are automatically satisfied in this case. With such choice, the instability of the homogeneous state is a pitchfork bifurcation, meaning that at the onset of instability, the imaginary part of the unstable(s) mode(s) ω k is zero.
Among the unstable λ k (i.e. those associated with a positive ω k ), we can select the most unstable one (i.e. the one associated with the largest value of ω k ) in the dispersion relation Eq. (12). This wavelength λu is the one which grows with the largest rate and is therefore likely to be the observed one in a numerical simulation of Eq. (10). However, this result should not be considered as rigorous given that there may be other non-linear effects selecting the wavelength which escape this analysis. To complement Fig. 2A of the main text, we show on Fig. S1 a contour plot of the values of λu in the (KI /KA, Dm/D) phase space. In the Turing limit (captured by Dm/D 1 in Fig. 2A of the main text), λu depends very weakly on the ratio KI /KA given that this ratio has to be large for the instability to occur in the first place. Therefore λu is mostly controlled by the turnover rates (fixed in Fig. S1) in this case. In the Keller-Segel limit (captured by KI /KA 1 on Fig 2(a) of the main text), the most unstable wavelength depends on both the turnover rates (fixed) and the ratio Dm/D which explains its variation along this axis. Note that the number of patterns shown in direct numerical simulations shown in Fig. 2B of the main text satisfy this linear prediction of the total number of patterns because they are close to the onset of instability.
To complement this analysis, we also show that considering non-vanishing χI , but remaining within the admissible bounds of Eq. (13), can lead to the existence of a third phase representing a Hopf bifurcation where the solution becomes oscillatory. Although interesting, this goes beyond the scope of our analysis, which concentrates on steady-state spatial patterns, and live-imaging would be needed to verify whether complex time-oscillations can occur during pattern formation in embryogenesis. We illustrate however this dynamic behavior in Fig. S2, in which we also show direct numerical simulations of Eq. (10) to justify this claim. We did not find other type of instabilities than the Turing, Keller-Segel and oscillatory ones regardless of the sign of χA and χI in the case of a classical activator-inhibitor scheme.

E. Cross-diffusion Turing patterns.
Another case of special interest is when the cell volume fraction sensitivities to morphogen levels are negative (χA,I < 0), eliminating the possibility of up-hill morphogen diffusion at the origin of the Keller-Segel instability in Eq. (12).
To simply explain the physical nature of the instability developing here, we first set τ = 0 (leading to lm = 0) such that the tissue volume fraction always locally remains at its homeostatic value. We also assume that the sensitivities of the volume fraction are small χA, I 1, but that the effective fickian and hydrodynamic diffusivities Du ∼ Dmug remain finite and of the same order of magnitude. Here the aim is not to quantatively justify these limits from a biological standpoint, given the lack of quantitative measurements of χA,I , but our goal is to give a qualitative explanation of the instability in a mathematically tractable setting (although the core mechanism we describe does not depend on this assumption).
In this limit, the conditions for linear stability of the homogeneous system Eq. (13) reduces to the classical ones while the dispersion relation Eq. (12) reads: Interestingly, this relation of dispersion has been studied for monophasic reaction-diffusion systems with ad hoc cross-diffusion terms (40). We also refer to (41) for further exposure of the potential physical origin of cross-diffusion terms in reaction-diffusion systems. As shown in (40), such cross-diffusion terms result in a dramatic broadening of the patterns phase space. In particular, any two-morphogen reaction scheme (inhibitor-inhibitor, activator-activator, etc.) has the potential to generate spatial patterns and not just the classical activator-inhibitor scheme.
F. Generalized Turing instability case. We now get back to the general case of system Eq. (6) but assume that χA,I = 0 such that cell volume fraction is not sensitive to morphogen concentration. In this case, the last equation uncouples and relaxes to φ ≡ φ * such that Eq. (6) reduces to, [15] The linear stability of the homogeneous solution X i ≡ X * i , X e ≡ uX * i of Eq. (15) reduces to the dispersion relation [ 16] In the main text, we consider the double limit: While it is clear from a biological standpoint that , £ ≈ 10 −2 s −1 t ≈ 10 −4 s −1 , (37), it is not clear that £.
However, we note that the Turing instability typically holds in this limiting case. As a matter of fact, the opposite limit where £ does not lead to patterns. Indeed, in the case £ t, Eq. (16) reduces to: which cannot lead to a patterning instability in general as it corresponds to a case where the diffusion coefficients for the activator and the inhibitor are the same. This result is interestesting as it might have particular biological implications, leading for example to the prediction that cells might generate stable spatial patterns of morphogens mostly in up-or down-regulating their endocytosis machinery in order to modulate the transmembrane transport equilibrium constants and consequently the values of morphogen effective diffusion coefficients.
G. Single morphogen case. We now consider the general case where a single morphogen exists which we chose by convention to be A. In that case, Eq. (6) simplifies to: [17] where f (Ai) = (A * i − Ai)/τA is the linearized form of f .

G.1. Robustness with respect to the main text limit.
In the main text, system Eq. (17) is studied in the limit where γA λA 1/τA and τ /τA ∼ KA 1 and thus reduces to: [18] However, an explicit criterion for the linear stability of Eq. (17) can be obtained where instead of such complex limit, we simply postulate (only for sake of simplicity) that φ * 1. While γA, λA 1/τA is verified in biological situations, this is particularly helpful to demonstrate that the instability we find in the single morphogen case does not rely on the assumption γA λA. The condition for an instability to occur in Eq. (17) (when the characteristic domain size l √ DτA) reads, which can be further reduced to the condition formulated in the main paper: [20] in the biologically meaningful case (see Section G.3) where KAγAτA 1. To further demonstrate this claim we present on

G.2. Scaling properties of Keller-Segel patterns.
In this section, we consider only the special case where τA → ∞ and KA → 0 while KADm, KAD and τ Dm remain finite. This limit corresponds to a case where the morphogen turnover rate is very slow but where the effective fickian and hydrodynamic diffusivities remain nonetheless finite. In such case, Eq. (17) reduces to a classical Keller-Segel system (42): but with a slightly different dynamics (φAi instead of Ai in the time derivative) and a non-linear dependence of φ on Ai through φ h . A very similar equation was derived and studied in the context of active gels (38,43).
Integrating the first equation of Eq. (21), steady states of Eq. (21) are given by where h = |Ω| −1 Ω h( r)d r denotes the spatial average over the tissue domain Ω. For simplicity, we now consider Eq. (22) on the finite segment [0, l] for which λ k = (πk/l) 2 , and investigate the possible bifurcations arising from the homogeneous state when the volume fraction sensitive parameter χA increases. Using the continuation software AUTO (44), we numerically compute the first three branches of such bifurcation diagram on Fig. S4 (a). It features an infinite number of critical points χ k A = D/Dm(1 + l 2 m λ k ) from where inhomogeneous solutions bifurcate, see Fig. S4 (c). Similar to (38,43), we numerically show that starting from an arbitrary distribution of Ai and if an instability exists, the system Eq. (21) will always coarsen to a single "two-zones" pattern corresponding to the fundamental mode of the instability λ1. The reason for this is that morphogen concentration creates a non-local convergent (since χA positive) flow recruiting even more morphogens over the hydrodynamic length lm and further increasing the flow. This positive feedback loop would lead to infinite localization of the morphogen concentration if morphogen diffusion was not preventing such a blow up by penalizing the creation of sharp gradients in concentration. If χA is small, the diffusion dominates and the concentration of morphogen is homogeneous but when χA reaches the critical threshold χ 1 A , the morphogen self-organize into a fundamental pattern where active advective transport balances diffusion. This is known as the Keller-Segel instability. Because the positive feedback loop is non-local, the morphogens concentration always coarse-grains into the fundamental mode. Interestingly, we found using a classical normal mode analysis, and confirmed numerically, that the nature of this first bifurcation branch is of pitchfork type but can be supercritical (second order phase transition) or subcritical (first order phase transition) cf. Fig. S4 (b). In practice, this means that two steady states corresponding to a non-patterned and a patterned state may co-exist for a range of parameters in the subcritical case. The transition from one state to the other could be then controlled by additional regulatory pathways.
Next, we define the fundamental pattern size as lpat where the Heaviside function selects the part of the morphogen profile larger than its average value. Therefore lpat = 0 for the homogeneous solution and becomes non-zero when a pattern forms. We show on Fig. S5 that this pattern scales as a function of the domain size for various values of χA. These results demonstrate that when an instability exists, a single "two zones" pattern always forms and its size scales with an increase in tissue domain size. This scaling property is at odds with what is observed for a classical Turing instability, where additional patterns form as the domain size increases. Thus, in this limit, our model constitutes a suitable theory to explain the scaling properties of various "two-zones" developmental patterns such as the dorso-ventral pattern of Xenopus (45), the antero-posterior pattern of planarian (46) or the left-right pattern of zebrafish (47).

G.3. Biological relevance of the Keller-Segel instability.
Here, we briefly review the orders of magnitude from the main text, to justify in particular that the Keller-Segel-Turing phase diagram proposed in Fig. 2A of the main text is relevant given the reported orders of magnitude for the model parameters. We focus here on the Keller-Segel instability and Eq. (19). The hypothesis KAγAτA 1 is rather straightforward to fulfil, as typical values measured in the literature are τA ≈ 10 4 − 10 5 s (47), while γA ≈ 10 −1 − 10 −2 s −1 (37). Therefore, even for small values of KA, the approximation is expected to hold (as values at which it would break down, KA < 10 −3 , would imply that almost no morphogen is present in the extracellular space anymore). This shows that the instability condition derived in the main paper for the Keller-Segel instability is in fact robust with respect to the assumption that KA 1, as opposed to the Turing instability with re-normalized diffusion coefficients which require KA 1.
Thus, one can adopt the simplified instability criteria Eq. (20), where D Dm , τ K A τ A and √ χA must be compared. Because χA is the least experimentally known of the model parameters, we start by focusing on the first two quantities. Orders of magnitude from the main text allow us to estimate τ /(KAτA) ≈ 0.01 − 0.1. Moreover, Dm can adopt a wide range of values given the wide range of length scales of extracellular gaps in different tissues, but its value is consistently larger to much larger than D, with D/Dm ≈ 10 −4 − 1, which tends to say that that τ /(KAτA) is the main factor resisting the Keller-Segel instability in the limit of large hydrodynamic diffusion Dm (i.e. when the intercellular spaces are ≈ 0.1 µm). In this case, the cell volume sensitivity to morphogen concentration, χA, must be larger than 0.01 − 0.1 for the Keller-Segel instability to develop, which, although relatively little data exists on variations of φ in embryonic tissues, is a rather lax criterion, as it represents typical variations of φ of only 1 − 10% with morphogen concentration. One should note that for very compact tissues, where values of D/Dm could be close to unity, the instability could require unrealistically large values of χA and thus would be screened.

H. Linking FRAP and FCS measurements of the diffusion coefficient.
Fluorescence recovery after photobleaching (FRAP) is a method which has been long used to experimentally determine diffusion coefficients (48). More particularly, it has been employed recently in several in vivo systems such as the early Zebrafish (47) and Xenopus (45) embryos or the Drosophila wing disk (49) in order to measure tissue-wide effective diffusion coefficients of various morphogens. Nevertheless, establishing a correspondence between these FRAP measured diffusion coefficients and the effective diffusion coefficient defined in our modeling framework requires some care. Indeed, in such FRAP experiments, the endogeneous molecule of interest is usually replaced by a genetically engineered copy which has been fused to a fluorescent protein such as the GFP. Thus, if one wants to study simultaneously two interacting morphogens, one needs to genetically engineer two constructs with two different fluorescent protein tags. In practice this is rarely the case and the diffusive properties of each morphogen are studied idependently (45,47). Moreover, such an approach does not directly discriminate the internal and external fractions of a given morphogen. Therefore, the intensities of fluorescence are proportional to the total concentration of a given morphogen : At = φAi + (1 − φ)Ae for A and It = φIi + (1 − φ)Ie for I. Combining the conservation equations Eq. (2) and Eq. (3), we obtain: For a small enough bleached region, the contribution of morphogen turnover rates can be neglected and the fluorescence recovery is mostly driven by transport, as shown in the case of the Zebrafish embryo (47). For the sake of simplicity, we shall thus proceed under this assumption, although some other in vivo systems may require to account for morphogen turnover rates (50).
The classical Turing case corresponds to a fixed volume fraction φ = φ * implying ve = 0. As equations for A and I uncouple, we can focus on a single morphogen only (say A)and the fluorescence At satisfies: [24] From Eq. (24), we can deduce the following dispersion relation: where λ k are the eigenvalues of the Laplace operator over the bleached region of typical size L b and ω the rates associated with fluorescence recovery. Non dimensionalizing times with γ −1 A and introducing the non-dimensional quantity ∆ = D/(γAL 2 b ) we can rewrite the above equation as: whereω are non-dimensional rates and ρ k = L 2 b λ k is a positive non-dimensional geometric factor which increases quadratically with the order of the eigenvalue considered. For instance if the bleached region is spherical, ρ k , is the discrete set of zeros of the Bessel functions of the first kind (48). If the bleached region is ∼ 100 µm (47), we have that L 2 b γA ∼ 10 −9 m 2 s −1 suggesting that ∆ 1 in Eq. (25). Then at first order we havē ω2 is an homogeneous clearance rate of fluorescence that is important only at very short timescales whileω2 corresponds to the measured diffusion driven recovery with effective diffusion coefficient, To relate the effective diffusion coefficient D FRAP to the local diffusion coefficient D Fick measured by Fluorescence Correlation Spectroscopy (FCS), one should recall that D is already a homogenized over a RVE and as such, depends on the tissue packing fraction. A classical way to relate D and φ * is through the Maxwell formula (51): Maxwell's result formally only holds in the dilute approximation where φ * is close to zero but in fact corresponds to the Hashin-Shtrikman (52) upper-bound which is known to hold outside the dilute limit: Such bound can be computed irrespective of the microscopic details of tissue geometry meaning that the actual effect of tortuosity is in reality even more pronounced. One can then directly compare the diffusion coefficient measured by FRAP and the one measured by FCS. For instance, in our first case where ∆ 1, we obtain: Interestingly, Lefty diffusion coefficient in the zebrafish embryo varies only by a factor 2 between FCS and FRAP measurements (47), D FRAP /D Fick 0.5. This can hold in the limit of K Lefty 1, since in this case we obtain D FRAP ≤ D Fick 1+φ * /2 , which is consistent with the statement of Muller et al. in (47) that there is "an approximately two-to three-fold decrease in diffusivity in the presence of cells compared to unhindered diffusion in aqueous solution". Note that this last bound (K Lefty 1 limit) also corresponds to the one expected for a passive extracellular molecule which does not chemically interact with cells such as an injected dye or a secreted recombinant fluorescent protein (see (50) for experiments of this kind in early Zebrafish embryos). However, if K Lefty 1 (concentrations outside of the same order of magnitude or smaller than concentration inside or immobilized), given the large cell volume fraction inferred in this system (φ * ≈ 0.8 − 0.9)(47), we expect the decrease in effective diffusion due to tortuosity to become much larger: D FRAP /D Fick ≤ 0.1. This argues for the need for precise measurements of K Lefty concomitant to φ. Another possibility is to consider the possibility of active contributions to morphogen transport (ve = 0 in Eq. (23)) since hydrodynamical diffusion would then, to leading order, rescale passive diffusion in an arbitrary manner as we show below.
If φ is no longer assumed to be homogeneous, variations of φ may contribute to the fluorescence recovery of the bleached region and one needs to fit the intensities At and It in Eq. (23) (with f = g = 0). A quantitative re-interpretation of the measurements of (50) in this context requires the spatial knowledge of φ and is thus left to future work. As a proof of principle, a simple case can be investigated where χI = 0 and τ → 0. In this situation Eq. (24) becomes which at linear order reads, The generalized dispersion relation then reads, Next, assuming χA 1 while χADm ∼ D are of the same order of magnitude, we can generalize the previous result obtained for χA = 0 to This shows that the active hydrodynamic transport could contribute to rationalize the diffusion coefficients measured by FRAP.
In particular, the inequality between D FRAP and D Fick becomes which could explain, among other possibilities of active morphogens transport (53), why in some cases even with φ * 0.85 and KA 1, D FRAP is only slightly smaller than D Fick .

I. Generalisation of the model in presence of tissue rearrangements.
In model Eq. (6), we have systematically assumed that changes of the tissue topology are much slower than the fluid circulation (vs 0) based on the estimates presented in Section A.3. However, these estimates may substantially vary depending on the biological system and cellular flows have been shown to take place at timescales comparable to fluid flows for various morphogentic processes (54). We therefore present below an extension of the theory accounting for such flows.
A first consequence of a non-negligible vs was already discussed in Section A.6 as Eq. (3) now contains a transport term The next step is to obtain the expression of ve as it enters in Eq. (1) and Eq. (2). For this, we have to generalize the mechanical framework presented in Section A.7 to obtain the analogue of Eq. (5) when the tissue can rearrange. We can do so by introducing a new variable is accounting for the non elastic strain present in the tissue due to the cells rearagements (55).
Thus, this strain is connected to the cell flow through the kinetic relation: dis dt = ∂tis + vs.∇is = ∇vs. [28] In this case, the constitutive behavior Eq. (4) becomes [29] If we continue to assume that the interstitial fluid flow in the tissue is described by Darcy's law and that the physiologically relevant limit Ku K G and α 1 still holds, we obtain using force balance that: Eliminating ve and denoting s = tr(is) we can rewrite the full model as [30] Such model generalizes Eq. (6) in the presence of a given cellular flow vs. Large scale cellular flows are ubiquitous during developmental morphogenesis (56) but their genesis is far from being understood and a general assessment of their influence on our system is not presently within reach. Thus, for the sake of clarity and to extract the key biophysical ingredients generating the instabilities described in the main text, we have decided to work in the limit where vs = 0.
However, we first note that if the cell flow is not associated with volumetric changes in the tissue (i.e. ∇.vs = 0) as it has been verified for example during the ventral furrow formation of Drosophilia (57) or in limb bud morphogenesis (58), then s = 0. In this case, the homogeneous solution Eq. (7) is still a solution of Eq. (30) which can be rewritten, . [31] where d/dt = ∂t + vs.∇ is the material derivative with respect to the cellular flow. As such, the only effect of the cellular flow is to add a transport term to Eq. (6). The impact of vs on the stability properties of Eq. (6) presented in the main paper can then be investigated case by case either by modeling vs or by at least checking its order of magnitude. For example for the limb bud morphogenesis vs ≈ 0.01µm/s (58,59), while from our estimates the typical velocity of the interstitial flow driven by spatial changes in the porosity Dm/τ ≈ 0.1 − 10µm/s vs. We therefore do not expect that cell rearrangements will substantially modify the thresholds of patterning instabilities given in the main text, as they are happening at a slower rate than our morphogen-driven pattern formation mechanism. A word of caution for this specific system though, as we note that our model requires the existence of at least a transient elastic tissue structure i.e the existence of some transient level of physical connections between the cells, either directly via some adhesion or indirectly via the Extracellular Matrix (ECM). Such an hypothesis and the measurements of the ensuing elastic moduli would need to be carefully scrutinized for validating the application of our model to this specific example.
A second simple modeling assumption that can be used to generalize our results is to suppose that morphogens are driving the cellular flow according to vs = ϑ.∇X e [32] where ϑ = (ϑA, ϑI ) is a vector of chemotactic sensitivities of the cell motion with respect to the morphogens concentrations. Such a chemotactic role of morphogens has been demonstrated in particular examples such as Shh in Drosophila (60). In this case, Eq. (30) is a closed system of equations which stability can be investigated. Clearly, all the instabilities presented in the paper will remain valid for small enough ϑ. In the more general case, new instabilities associated with chemotaxis will appear. In particular, we can illustrate this in the simple case where only one morphogen is present as in Section G.1 and in the limit considered in the main paper where γA λA 1/τA and τ /τA ∼ KA 1. If that case, if χA = 0, we have shown that no instability occurs in Eq. (18) as Eq. (20) is not satisfied. This is no longer true if we consider that ϑA > 0 as we find the onset of an oscillatory (Hopf) instability when: While the first inequality ensures that chemotaxis is not too strong, such that it would make an infinite number of high modes unstable by surpassing Fickian diffusion, the second inequality shows that chemotaxis alone can be enough to destabilize the regularizing effects of both Fickian and hydrodynamic diffusion as well as the morphogen chemical turnover. In general, that type of instability will combine with the ones presented in the main text but a systematic investigation of such cross-effects is left to future works.

Historical background
In 1948, Turing's long time interest for biology took a new twist, largely influenced by the works of D'Arcy Thompson and Conrad Waddington (61). In this context, he decided to study pattern formation and morphogenesis during embryonic development, developing computer programs to numerically solve the complex partial differential equations resulting of his abstraction of the biological question. In a letter dated of 12 February 1951, he wrote to a former colleague and friend, Mike Woodger (62): "Our new machine is to start arriving on Monday. I am hoping as one of the first jobs to do something about 'chemical embryology'." These "jobs" would be the foundations of Turing's famous paper "The Chemical Basis of Morphogenesis" (63). In this article, he proposed an elegant mathematical model of self-organisation, where two reacting chemicals can spontaneously form periodic spatial patterns through an instability driven by their difference in diffusivity. This idea was not only very influencial amongst biologists such as J.M. Smith (64) and L. Wolpert (65), or physicists such as I. Prigogine (66), F. Crick (67) and H. Meinhardt (35), but was at the stem of a new class of PDE models which became central in the development of modern mathematical biology (39,68). However, given the lack of molecular and cellular evidences available at that time, in a lot of studies, the emphasis was often placed on the striking similarity of the patterns obtained by the model with those observed in biological systems such as animal coat and skin pigmentation patterns, rather than on discussing the biological validity of the modeling assumptions. Quite astonishingly, in a letter to Alan Turing dated of 11 September 1952, Waddington already raised similar concerns about the applicability of Turing's reaction-diffusion model to biological developmental systems, questioning its limitation to reproduce some observed behaviours in embryonic development such as pattern scaling with tissue size or the generation of a spatial pattern of discrete cell types (69): "I was extremely interested to read your recent paper on the chemical basis of morphogenesis. [. . . ] I rather doubt, however, whether the kind of processes with which you were concerned play a very role in the fundamental morphogenesis which occurs in early stages of development. [. . . ] Even in a case like the regeneration of tentacles in Hydra the final result seems to me more regular that one would expect from your type of mechanism. For instance, if one splits two hydras longitudinally and places them together so that the split edges heal and forms a single tube of double the normal diameter, this enlarged tube still regenerates the normal number of tentacles, although it would, of course, contain twice the normal 'chemical wave-lengths' around its circumference." Turing was not naive about the limitations and underlying assumptions of his model and in fact already hints, in his article, at possible extensions to his theory considering the structure and mechanical properties of biological tissues (63): "[...] one proceeds as with a physical theory and defines an entity called 'the state of the system'[. . . ] The description of the state consists of two parts, the mechanical and the chemical [. . . ] In determining the changes of state one should take into account: (i) The changes of position and velocity as given by Newton's laws of motion. (ii)The stresses as given by the elasticities and motions, also taking into account the osmotic pressures as given from the chemical data (iii) The chemical reactions (iv) The diffusion of the chemical substances. The region in which this diffusion is possible is given from the mechanical data. [. . . ] The interdependence of the chemical and mechanical data adds enormously to the difficulty, and attention will therefore be confined, so far as is possible, to cases where these can be separated." The original work of Hans Othmer on the transition between discrete and continuous models of pattern formation in tissues in the 1970's (70) is maybe the closest example of this new research program tangentially set by Turing in his paper. In introduction of a review article entitled "Current problems in pattern formation" Othmer wrote (71): "In Turing's work and virtually all subsequent work, the only mode of transport considered is diffusion. Furthermore, the model systems all deal with structureless, tightly coupled cells or their continuum analogs. Some future work should be directed toward: (a) the analysis of other modes of transport (b) more realistic models of cell and tissue structure and, (c) networks of cells that communicate only indirectly via the external medium." As early attempts to fulfil this program, we would like to highlight the work of L.E. Scriven, who used the coupling of mechanical forces with chemical reactions to describe spontaneous fluid transport (72), paving the way to modern theories of active gels, and the pioneering work of G.F. Oster, one of the fathers of mechanobiology, on limb growth and patterning (73).
A. Methods. Linear stability analysis is performed numerically using conventional algebra packages in Mathematica. Numerical integration of models Eq. (10), Eq. (17) and Eq. (21) are performed in 1D with a custom made code written in Matlab and available on request to the authors. In brief, the code uses a mass conserving finite volume stencil in space and a first order implicit scheme in time.