Ordered Hexagonal Patterns via Notch-Delta Signaling

TL;DR

This paper presents an exact mathematical analysis of pattern formation via Notch-Delta signaling, explaining the prevalence of hexagonal cellular patterns and proposing mechanisms for defect-free pattern generation.

Contribution

It introduces a detailed differential equation model for juxtacrine signaling and analyzes bifurcations to explain observed biological patterns and their formation mechanisms.

Findings

Hexagonal patterns with high Delta at centers are prevalent.

Low cis-coupling leads to novel high Delta and high Notch patterns.

Biological systems require additional mechanisms for defect-free pattern formation.

Abstract

Many developmental processes in biology utilize Notch-Delta signaling to construct an ordered pattern of cellular differentiation. This signaling modality is based on nearest-neighbor contact, as opposed to the more familiar mechanism driven by the release of diffusible ligands. Here, exploiting this "juxtacrine" property, we present an exact treatment of the pattern formation problem via a system of nine coupled ordinary differential equations. The possible patterns that are realized for realistic parameters can be analyzed by considering a co-dimension 2 pitchfork bifurcation of this system. This analysis explains the observed prevalence of hexagonal patterns with high Delta at their center, as opposed to those with central high Notch levels. We show that outside this range of parameters, in particular for low cis-coupling, a novel kind of pattern is produced, where high Delta cells have high Notch as well. It also suggests that the biological system is only weakly first order, so that an additional mechanism is required to generate the observed defect-free patterns. We construct a simple strategy for producing such defect-free patterns.

Equations

$${\dot{N}}_x$$

$${\dot{D}}_x$$

$${\dot{I}}_x$$

$$H_{+}(I)=1+\frac{k_H{I}^{n_{+}}}{I^{n_{+}}+s_0^{n_{+}}};H_{-}(I)=\frac{s_0^{n_{-}}}{I^{n_{-}}+s_0^{n_{-}}}$$

$$k_t{\lambda }_I{N}_0{D}_0={\gamma }_I{I}_0$$

$$N_0=\frac{{\lambda }_N{H}_{+}(I_0)}{(k_c+k_t)D_0}$$

$$I_0=\mathcal{R}(I_0)\equiv \frac{k_t{\lambda }_N(1+k_H\frac{I_0^2}{1+I_0^2})}{{\gamma }_I{k}_0}.$$

$${\dot{u}}_n=-u_n-\gamma u_n^3+{\Omega }_0+ϵu_{\text{ext}}$$

$$\dot{A}$$

$$\dot{B}$$

$$\dot{C}$$

$$0=-H-\gamma H^3+{\Omega }_0+ϵH$$

$$A$$

$$B$$

$$C$$

$${\dot{C}}_1$$

$${\dot{C}}_2$$

$$0=-\frac{\Delta ϵ}{2}+3\gamma H{C}_1-3\gamma C_1^2$$

$${\lambda }_D^c\approx 1.931+0.262\delta -0.0236{\delta }^2+0.000676{\delta }^3$$

$${\dot{C}}_1$$

$${\dot{C}}_2$$

$$\alpha =0.04440ϵ;\sigma =0.000382\delta ;\mu =-\mathrm{0.0003063.}$$

$${\tilde{\lambda }}_N$$

$${\tilde{\lambda }}_D$$

$$N_0$$

$$D_0$$

$$I_0$$

$$k_0\equiv k_c+k_t.$$

$$N_0$$

$$D_0$$

$$I_0=\frac{k_t{\lambda }_I}{{\gamma }_I}\frac{-k_0\left({\tilde{\lambda }}_D+{\tilde{\lambda }}_N\right)-{\gamma }^2-\sqrt{{\left[k_0\left({\tilde{\lambda }}_D+{\tilde{\lambda }}_N\right)+{\gamma }^2\right]}^2-4k_0^2{\tilde{\lambda }}_D{\tilde{\lambda }}_N}}{2k_0^2}\equiv \mathcal{R}(I_0).$$

$${\cal M}_{U}=\left(\begin{array}[]{ccc}-k_{0}D_{0}-\gamma&(-k_{c}+k_{t}/2)N_{0}&\frac{d\tilde{\lambda}_{N}}{dI_{0}}\\ (-k_{c}+k_{t}/2)D_{0}&-k_{0}N_{0}-\gamma&\frac{d\tilde{\lambda}_{D}}{dI_{0}}\\ k_{t}\lambda_{I}D_{0}&-(k_{t}/2)\lambda_{I}N_{0}&-\gamma_{I}\end{array}\right)$$

$$k_t{\lambda }_I{D}_0\left[4\gamma +3(k_0+k_c)N_0\right]\frac{d{\tilde{\lambda }}_N}{d{I}_0}=4{\gamma }_I\gamma \left[\gamma +k_0\left(D_0+N_0\right)\right]+3{\gamma }_I{N}_0{D}_0{k}_t\left(k_0+3k_c\right)+2\frac{d{\tilde{\lambda }}_D}{d{I}_0}{k}_t{\lambda }_I{N}_0\left(\gamma +3k_c{D}_0\right)$$

$$0=ϵ\frac{\partial f_i}{\partial x_j}{x}_{j,1}+\frac{{ϵ}^2}{2}\frac{{\partial }^2{f}_i}{\partial x_j\partial x_k}{x}_{j,1}{x}_{k,1}+\Delta \mu \frac{\partial f_i}{\partial \mu }+\frac{{ϵ}^3}{6}\frac{{\partial }^3{f}_i}{\partial x_j\partial x_k\partial x_l}{x}_{j,1}{x}_{k,1}{x}_{l,1}+ϵ\Delta \mu \frac{{\partial }^2{f}_i}{\partial x_j\partial \mu }{x}_{j,1}$$

$$\frac{{\partial }^2{f}_i}{\partial x_j\partial x_k}{x}_i^L{x}_j^R{x}_k^R=0.$$

$$x^{R}=\left(\begin{array}[]{c}4\gamma_{I}(\gamma+k_{0}N_{0})+2k_{t}\lambda_{I}N_{0}\frac{d\tilde{\lambda}_{D}}{dI_{0}}\\ -4\gamma_{I}\left(k_{c}-\frac{k_{t}}{2}\right)D_{0}+4k_{t}\lambda_{I}D_{0}\frac{d\tilde{\lambda}_{D}}{dI_{0}}\\ k_{t}\lambda_{I}D_{0}\left(4\gamma+6\left(k_{c}+\frac{k_{t}}{2}\right)N_{0}\right)\end{array}\right);\qquad x^{L}=\left(\begin{array}[]{c}k_{t}\lambda_{I}D_{0}\left(4\gamma+6\left(k_{c}+\frac{k_{t}}{2}\right)N_{0}\right)\\ -2k_{t}\lambda_{I}N_{0}(\gamma+3k_{c}D_{0})\\ 4\gamma^{2}+4k_{0}\gamma(D_{0}+N_{0})+3k_{t}(k_{0}+3k_{c})D_{0}N_{0}\end{array}\right)$$

Full text

Notch/Delta signaling;hexagonal patterns;pattern formation

I Introduction

Biological cells can exist in a number of distinct phenotypes, even with a fixed genome. These phenotypes arise via multi-stability of the underlying dynamical network controlling cell behavior and allow cells to take on differentiated roles in overall organism function. It is clear that developmental processes must ensure that these phenotypes arise in the right place and time, i.e., ensure the emergence of functional phenotypic patterns.

A well-studied case of such a system is that of Notch-Delta signaling Andersson et al. (2011). Various cells contain Notch transmembrane receptors Gazave et al. (2009) that couple to Notch ligands such as Delta or Jagged on both the same cell (cis-coupling) and neighboring cells (trans-coupling). Because of the manner by which Notch and Delta inhibit each other (see below), their interaction typically leads to an alternating “salt and pepper” structure. This type of patterning is seen in systems ranging from eyes Artavanis-Tsakonas et al. (1995) and ears Neves et al. (2013) to intestines Chen et al. (2017) and livers Kaylan et al. (2018). As a general rule, the high Delta cells are the most specialized ones (for example, the photoreceptors Fanto and Mlodzik (1999)) and are surrounded by less differentiated high Notch supporting cells. Parenthetically, changes in the transcriptional regulation utilizing the Delta-alternative Jagged ligand may be crucial for the role of Notch in cancer metastasis Boareto et al. (2016); Zhu et al. (2013), but here we focus solely on Delta and its interplay with Notch.

In this paper, we study two aspects of the Notch-Delta system on 2d hexagonal arrays of cells. A motivation for this choice is that epithelial layers consist of polygonal cells that roughly form a hexagonal lattice, albeit with some size dispersion and some defects (see Fig. 1). It is worthwhile to first work out how pattern formation works in the more idealized case of the perfect lattice and then afterwards consider possible effects of the irregularities; we do note in passing that some effort has already been devoted to understanding the role of variations in cell size Shaya et al. (2017). First we focus on the existence and stability of hexagonal patterns in this geometry, which allows an exact re-writing of the (ordered) pattern-forming problem as a nine-dimensional dynamical system. We numerically construct the phase diagram of the system, showing the variety of patterns that emerge. Many features of this system can be understood by expanding about a co-dimension two pitchfork bifurcation. Others can be captured by examining various limiting cases of the parameters. The second aspect concerns mechanisms for ordered patterns to emerge from generic initial conditions. Here we identify a possible role for an initiating wave, similar to what has been seen in at least some biological realizations Corson et al. (2017).

The Notch-Delta interaction is an example of juxtacrine (i.e., contact-dependent) signaling. As sketched in Fig. 2, Ligands such as Delta bind Notch receptors and, when this occurs between neighboring cells, leads to the cleavage of the receptor and release of its intracellular domain (NICD). NICD translocates to the nucleus where it transcriptionally up-regulates Notch and down-regulates Delta. The ligand-receptor interaction between molecules on the same cell leads to mutual annihilation with no NICD release Sprinzak et al. (2010). The combination of cis-annihilation and NICD-mediated trans-repression is responsible for the observed lateral inhibition of Delta Collier et al. (1996). We will use a baseline model Formosa-Jordan and Sprinzak (2014) of this process involving three concentrations, $N$ (receptor), $D$ (ligand) and $I$ (NICD),

[TABLE]

Here, positions $x$ refer to locations on an hexagonal lattice (see Fig. 3) and the superscript “ext” refers the average over the six nearest neighbor sites of $x$. The production terms $H_{\pm}$ corresponding to the aforementioned transcriptional regulation are taken to be Hill functions,

[TABLE]

such that $H_{\pm}(0)=1$ and $H_{+}$ is an increasing function that saturates at $1+k_{H}$, while $H_{-}$ is a decreasing function that decays to 0 with increasing $I$. We define a typical set of parameters taken from the literature Boareto et al. (2015): $\gamma=0.1$. $\gamma_{I}=0.5$, $n_{+}=n_{-}=2$, $k_{c}=0.1$, $k_{t}=0.04$, $\lambda_{I}=1$, $k_{H}=1$, $s_{0}=1$, and primarily focus on the role of $\lambda_{N}$ and $\lambda_{D}$. Nevertheless, we will also have occasion to investigate the effects of changing $k_{H}$ and $k_{t}$.

II Uniform patterns and the Region of Instability

We start by considering the simplest type of solution for this system of equations, namely one which is spatially uniform. In general, solving for such a solution requires the simultaneous solution of three nonlinear equations for $N_{0}$, $D_{0}$ and $I_{0}$. In Appendix A, we show that one can reduce the problem to the solution of one rather complicated equation for $I_{0}$; this is given as Eq. A5. For general values of the parameters, this last equation must be solved numerically.

We can gain some analytic insight into the solution space by considering various limiting cases (for full details see Appendices B and C). One convenient such case is that of large $\lambda_{D}$. In that limit, we expect $D_{0}$ to be large, which then forces $N_{0}$ to be small to satisfy the condition arising from the last of Eqns. 1,

[TABLE]

From the Notch equation, we have in this limit

[TABLE]

which then immediately leads to the condition

[TABLE]

with $k_{0}\equiv k_{c}+k_{t}$ and where we have specialized to the case of Hill functions with $n_{\pm}=2$. From this equation we can determine when there is a unique solution versus when multiple solutions exist. The details of this calculation are presented in Appendix D. For our purposes here, we note that the range of parameters for which there is a unique uniform solution ranges over all reasonable values of the system parameters.

We are interested here in hexagonal patterns on our hexagonal lattice. These patterns arise as the uniform state becomes unstable with respect to spatially varying perturbations. In fact, our numerical experiments on Eq. (1) indicate that the first instability of the uniform state for the set of parameters above is almost always to a hexagonal mode, as that arrangement maximizes the average number of ”satisfied” nearest neighbor inhibitory interactions. This is of course a consequence of studying the model on a hexagonal grid, as motivated by the biological application. Then, the immediate question is where in parameter space the uniform state is unstable to a hexagonal pattern. The condition for this is presented in Eq. (19) in Appendix A. The region of instability in the $k_{t}$, $\lambda_{D}$ plane is presented for various values of $\lambda_{N}$, both for $k_{H}=2/3$ in the left panel of Fig. 4 and for $k_{H}=4/3$ in the right panel. For each value of $\lambda_{N}$ considered, the unstable region is that above the corresponding curve, i.e. $\lambda_{D}>\lambda_{D}^{c}(k_{t})$. For each value of $\lambda_{N}$, the value of $\lambda_{D}^{c}$ diverges at two values of $k_{c}$, with the instability only possible between these two values. For $k_{H}=2/3$, the curve $\lambda_{D}^{c}(k_{t})$ doubles back on itself for low $k_{t}$, but not near the upper limiting value of $k_{c}$, whereas for $k_{H}=4/3$, this doubling back happens on both sides. This qualitative change in behavior is explained in Appendix C, and occurs at $k_{H}=1$. Clearly wherever there is a hexagonal instability, a hexagonal solution exists. However, as we shall see, there is also the possibility of coexistence of a stable uniform solution with a hexagonal solution, so that hexagonal solutions extend outside the region of linear instability of the uniform state.

III Hexagonally Ordered Patterns

We now turn to the construction of hexagonally ordered patterns whose existence (but not stability) is guaranteed by the linear instability discussed above. These patterns are invariant under translation with vectors $\pm 6\hat{x}$, $\pm 3\hat{x}\pm 3\sqrt{3}\hat{y}$ where $\hat{x}$ and $\hat{y}$ are unit vectors along the coordinate axes and the unit of length is 1/2 the length of one of the hexagonal sides. From Fig. 3, it is clear that the fields everywhere are completely determined by their values on three sub-lattices that we have labeled A, B and C. This means that the entire problem of hexagonally ordered pattern structure (and their stability with respect to modes invariant under translation of the A-B-C unit cell) is reduced to nine coupled ODE’s.

This exact mapping is very different than what occurs for more traditional pattern formation problems Cross and Hohenberg (1993) such as for example for convection rolls Busse (1978), where the reduction to a set of ODE’s is valid only as an approximation near the bifurcation point. This ODE reduction has been performed previously by Negrete and Oates Negrete, Jr. and Oates (2019) for their simple one-field model of the Notch-Delta system (see below for details on this model), and here it is essential for our analysis of the full realistic model, especially away from the region of the bifurcation. As already shown, at fixed $\lambda_{N}$, the uniform solution with the fields taking on the same values on all three sublattices becomes unstable for $\lambda_{D}>\lambda_{D}^{U}(\lambda_{N})$ via a transcritical bifurcation, i.e. an intersection with a nonuniform solution. On this new branch, the respective values of the fields on two sublattices (say B and C) are identical, differing from the values on the remaining (in this case, A) sublattice. This hexagonally structured solution has a 6-fold hexagonal symmetry about any site on the different (here, A) sublattice. The bifurcation is transcritical in general, because of the lack of any symmetry between positive and negative deviations of the fields from their uniform values. As already mentioned and unlike convection, the transition to rolls does not take place at the same parameter value as that to hexagons as the roll pattern necessarily has a different wavelength on the lattice; hence rolls do not compete with hexagons, at least near the bifurcation. For completeness, the 6 coupled equations that govern these ordered hexagonal patterns are given explictly in Appendix E

IV The Pitchfork Bifurcation and its Unfolding: The Negrete-Oates model

To understand the nature of the origin of hexagonal patterns in a hexagonal lattice system, Negrete and Oates Negrete, Jr. and Oates (2019) introduced a simple one-field model containing the same type of instability. The model is given by

[TABLE]

where again $u_{\textit{ext}}$ is the average of $u$ on the 6 nearest-neighbor sites. For sufficiently negative $\epsilon$, this model has a lattice analog of a negative diffusion constant and hence become unstable to hexagonal patterns. On the three sites of the unit cell introduced above, this system reduces to

[TABLE]

They noted the presence of a pitchfork bifurcation at a specific value of the parameters, namely $\epsilon_{\textrm{PF}}=-2$, $\Omega_{\textrm{PF}}=0$. The latter is directly dictated by the symmetry $u\leftrightarrow-u$ of the equations, which of course only is present at $\Omega_{0}=0$. It is useful to go through the exercise of working out the bifurcation analysis for this model, as the algebraic structure of our three-field more biologically realistic model is quite similar in structure but more complicated in detail.

The homogeneous stationary state $A=B=C\equiv H$ satisfies

[TABLE]

with solution $H\approx\Omega_{0}/(1-\epsilon)$ for small $\Omega_{0}$. At the critical value $\epsilon_{c}=-2-6\gamma H^{2}$, the stability matrix around the homogenous solution has two zero modes, $(\delta A,\delta B,\delta C)=(-2,1,1)$ and $(0,1,-1)$ and the nonzero mode $(1,1,1)$, with eigenvalue $-3$. We will be carrying out a weakly-nonlinear analysis in the neighborhood of the co-dimension two pitchfork bifurcation.point. To obtain this, we write

[TABLE]

where $C_{1}$ and $C_{2}$ are the ${\cal O}(\Omega_{0})$ amplitudes of the zero modes, and $\delta_{2}$ is an ${\cal O}(\Omega_{0})^{3}$ correction. Writing $\epsilon=\epsilon_{c}+\Delta\epsilon$, where $\Delta\epsilon$ is ${\cal O}(\Omega_{0})^{2}$, and expanding to third order in $\Omega_{0}$, we obtain three equations for the time-dependent amplitudes $C_{1}$ and $C_{2}$ (which vary on the slow time scale $\Omega_{0}^{2}$) and $\delta_{2}$. Eliminating $\delta_{2}$, we find the two amplitude equations

[TABLE]

where the time derivatives refer to variation on the aforementioned slow time scale. It is easy to see that all the terms in our equations are ${\cal O}(\Omega_{0})^{3}$.

We start by looking for time-independent solutions of these equations. One stationary state of the system is the original homogeneous state with $C_{1}=C_{2}=0$. There is a pair of solutions with $C_{2}=0$, with $C_{1}$ satisfying the quadratic equation

[TABLE]

These two solutions emerge from a saddle-node bifurcation occurring at $C_{1}=H/2$, $\Delta\epsilon=3\gamma H^{2}/2$. This lies on the stable side of the transition and thus the saddle node bifurcation precedes the instability of the homogenous state as $\epsilon$ is decreased.. As $\Delta\epsilon$ further decreases from its saddle-node value, one of the solution branches has increasing $C_{1}$, while the other approaches $C_{1}=0$, i.e., the homogeneous solution, intersecting it at $\Delta\epsilon=0$, the location of the homogeneous instability. It then crosses over to $C_{1}<0$, so that here it has “polarity” opposite to that of the other branch emerging from the saddle-node bifurcation that had $C_{1}>0$ and increasing.

There are other stationary solutions, having $C_{2}\neq 0$. Solving the $\dot{C}_{2}=0$ equation for $C_{2}$ and substituting in the $\dot{C}_{1}=0$ equation yields $C_{1}=-x/2$, where $x$ is the value obtained above for $C_{1}$ when $C_{2}=0$. Substituting this back into the equation for $C_{2}$ then reveals $C_{2}=\pm 3x/2$. For these solutions, the leading order value of $A-H$ is $-2C_{1}=x$, which was the leading order value of $B-H=C-H$ in our original solution. In addition, either $B-H$ or $C-H$ (depending on the sign of $C_{2}$) equals $x$ as well, with the other equalling $-2x$, which was the value of $A-H$ in the original solution. Thus, all these new solutions are simply the previous solutions translated to be centered on $B$ or $C$, instead of $A$. At the linear level, these new solutions are the linearly combinations $-1/2[(-2,1,1)\pm 3(0,1,-1)]$ of the two zero modes of the homogeneous solution. This crossing of this inhomogeneous solutions with the original homogeneous solution represents a two dimensional transcritical bifurcation of the homogeneous solution; two-dimensional here because of the two zero modes of the homogeneous solution. Precisely at $\Omega_{\textrm{PF}}$, i.e., $\Omega_{0}=0$, the saddle node merges with the transcritical point and the two branches meet symmetrically at $\Delta\epsilon=0$, i.e., $\epsilon_{\textrm{PF}}$, as indicated by the vanishing of the term linear in $C_{1}$ in Eq. (9).

We can also easily calculate the stability spectrum of the non-trivial patterned states of this reduced system. The homogeneous solution has a degenerate pair of modes, with growth rates $\Omega_{1}=\Omega_{2}=-\Delta\epsilon/2$, so that is stable below the transition, $\Delta\epsilon>0$ and unstable above. Focusing on the $C_{2}=0$ inhomogeneous solution, it has two eigenvalues, $\Omega_{1}=3\gamma C_{1}(H-2C_{1})$ and $\Omega_{2}=-9\gamma C_{1}H$. From this we can see that (for $\Omega_{0}>0$) on the positive branch both modes are stable. $\Omega_{1}$ crosses 0 at the saddle node, and on the second branch we have 1 stable and one unstable mode. When the second of the non-uniform branches crosses the homogeneous solution at the transcritical point, both $\Omega$’s cross zero and we remain with one stable and one unstable mode, switched compared to those on the previous side. The story for the other, shifted, solutions is of course the same. This picture is more intricate than for the standard co-dimension one pitchfork, where the stable solution on one side of the bifurcation gives rise to one unstable and two symmetry-related stable solutions on the other side.

V The Pitchfork Bifurcation and its Unfolding: The Biological Notch-Delta model

The shared symmetry structure with the Negrete-Oates model above suggests and numerical exploration of Eq. (1) confirms that our more realistic Notch-Delta system also possesses a co-dimension 2 pitchfork bifurcation. For our standard parameters the pitchfork occurs at at $\lambda_{N}^{\textit{PF}}=2.7582488$, $\lambda_{D}^{\textit{PF}}=1.93112$. Specifically, for $\lambda_{N}^{\textit{PF}}$, as is the case for general $\lambda_{N}$, only the uniform solution exists for $\lambda_{D}<\lambda_{D}^{U}$ and it is stable. At $\lambda_{D}^{U}=1.93112$, two additional solutions are born, one a “hexagon” (by definition, a solution where high $D$ is surrounded by high $N$) and one an “anti-hexagon” (high $N$ surrounded by high $D$). What is interesting is that this point occurs not only for physically possible (i.e. positive) values of the parameters, but within the range of parameter values determined (at least roughly) by experiment Sprinzak et al. (2010). Thus the unfolding of the bifurcation gives us detailed information about the pattern formation possibilities realizable in real physiological settings.

Many of the features in the vicinity of the pitchfork bifurcation of the Negrete-Oates model carry over into our more complicated system. Specifically, it can directly be shown by numerical analysis of the 9-dimensional reduced system that the uniform state has 2 (degenerate) unstable modes above the critical $\lambda_{D}$. The emerging hexagon branch is stable, whereas the anti-hexagon has one unstable mode. The instability is with respect to a mixed mode (defined as a mode with all three sublattices having different values) which converts the anti-hexagon to a shifted hexagon. This overall structure is shown in Fig. 5 A, and indeed recapitulates the structure determined for the Negrete-Oates model via the amplitude equation analysis.

At all other values of $\lambda_{N}$, the pitchfork breaks up into a transcritical bifurcation and a saddle-node (see Fig. 5B), also as in the Negrete-Oates model. Again by direct numerical solution, we find that for $\lambda_{N}>\lambda_{N}^{\textit{PF}}$, the uniform solution undergoes a transcritical bifurcation with a unstable anti-hexagon (with respect to a mixed-mode perturbation) on the high $\lambda_{D}$ side and an unstable (with respect to a pure-mode) hexagon on the low $\lambda_{D}$ side. The unstable hexagon then undergoes a saddle-node bifurcation, rendering the hexagon stable; this stable branch then continues on as $\lambda_{D}$ increases. For $\lambda_{D}$ smaller than the saddle-node value, no patterned solution exists. Hence, there exists a range of parameters for which a stable hexagon coexists with the stable uniform solution, a range which widens as $\lambda_{N}$ increases; we will return to this point below. For example, for $\lambda_{N}=3.5$, the transcritical bifurcation in which the uniform state goes unstable is at $\lambda_{D}=2.11356$, whereas the saddle node bifurcation is at $\lambda_{D}=2.1097$. An example of a stable hexagon solution for $\lambda_{N}=3.5$, $\lambda_{D}=2.3$ is shown in Fig. 6.

A similar bifurcation structure appears for $\lambda_{N}<\lambda_{N}^{\textit{PF}}$, where now the stable hexagon lies to the right of the transcritical point and the unstable anti-hexagon lies to the left, and it is the one that undergoes a saddle-node bifurcation (see Fig. 7A). The anti-hexagon is born with 1 unstable mode at the transcritical point and turns stable at the saddle-node bifurcation. However, unlike what happened in the Negrete-Oates model, citeNegrete. the stable anti-hexagon branch loses stability to a mixed-mode perturbation; this instability leads to a new pitchfork bifurcation, which is a result of the $B/C$ symmetry breaking. The hexagon, on the other hand, is born with one unstable mode and subsequently becomes stable, also as a result of a mixed-mode pitchfork bifurcation. An overall diagram of the stable phases as a function of the two parameters $\lambda_{N}$ and $\lambda_{D}$ is presented in Fig. 7B The mixed-mode solution branch arising from the hexagon bifurcation is the same solution which arises from the anti-hexagon bifurcation. For example, at $\lambda_{N}=2$, the hexagon becomes stable at $\lambda_{D}=2.0477$ and the anti-hexagon becomes unstable at $\lambda_{D}=2.056$. Thus, there is a very small coexistence region between the hexagon and anti-hexagon solutions. Again, the only solution that survives stably to higher values of $\lambda_{D}$ is the hexagon. This is in accord with the general biological rule given above that the high Delta cells are surrounded by high Notch cells sufficiently far from $\lambda_{N}^{\textit{PF}}$ and its associated $\lambda_{D}^{U}$. From the physics perspective, the explicit lack of symmetry between Notch and Delta as reflected in this model is not eliminated by working close to the co-dimension two bifurcation since the accidental symmetry at this point affects only the leading order term in the amplitude equation, not any of the higher-order ones. This feature is not captured by the simpler one-field model where the model has an exact symmetry at $\Omega_{0}=0$.

We can again use weakly non-linear bifurcation theory to flesh out these numerical findings, working in the immediate vicinity of the pitchfork bifurcation, The bifurcation analysis for this more complicated system of 9 equations can be performed, and after eliminating the seven fast modes (as opposed to just one previously), we get a set of amplitude equations the two slow modes, parameterized by $\delta_{B,C}=N_{B,C}-N_{0}$, where $N_{0}$ is the homogeneous Notch level. Doing this, we get exactly the same bifurcation equations structure as in Eq. 8 with the role of the symmetry breaking parameter $\Omega_{0}$ played by $\delta\equiv\lambda_{N}-\lambda_{N}^{\textit{PF}}$ and the role of the other bifurcation parameter $\Delta\epsilon$ replaced by $\epsilon\equiv(\lambda_{D}-\lambda_{D}^{c})$ where $\lambda_{D}^{c}$ is the location of the instability of the homogeneous solution. Near $\lambda_{N}^{\textit{PF}}$,

[TABLE]

so that $\lambda_{D}^{c}$ increases with $\lambda_{N}$. We can write the final amplitude equations as

[TABLE]

In term of these bifurcation parameters $\delta$ and $\epsilon$, the amplitude equation coefficients are

[TABLE]

This system is of course identical in structure to that we received in the Negrete-Oates model. The new information is the connection of the perturbations $C_{1}$, $C_{2}$ to the deviations from homogeneity of the physical fields $N$, $D$, and $I$.

The analysis of these equation thus follows directly from our previous analysis. Here, the saddle-node point is at $\epsilon=-0.00807\,\delta^{2}$, again below the transition, with $C_{1}^{\textit{SN}}=0.6244\,\delta$, so that for $\delta>0$ the $B$ and $C$ sites has high Notch, so that the $A$ site has low Notch. From the solution for the fast modes, we have $D_{A}=D_{0}+0.0276\,C_{1}$, so the $A$ site has high Delta. Thus the saddle-node solution is what we entitle a “hexagon” solution, which of course a distinction that is meaningless in the one-field model. On the other side of the transcritical point, $C_{1}$ changes sign, and so the $A$ site there has low Notch and high Delta, i.e. an anti-hexagon. For $\delta<0$, i.e. $\lambda_{N}<\lambda_{N}^{\textit{PF}}$, things are reversed, and the saddle node has $C_{1}<0$, so that the saddle-node solution is an anti-hexagon, and the solution on the other side of the transcritical point is a hexagon. This is of course consistent with our detailed numerical findings.

Checking the stability, both homogeneous modes have growth rate $\alpha$, and so are stable for $\lambda_{D}<\lambda_{D}^{c}$ and unstable above the transition. For the $\delta>0$ hexagon, one mode, with growth rate $\Omega_{1}$ is stable on the upper branch $C_{1}>C_{1,\textit{SN}}$ and unstable for $C_{1}$ below the saddle node. The other mode, with growth rate $\Omega_{2}$ is stable both above and below the saddle node. Thus, the hexagon on the upper branch is stable and on the lower branch is once unstable. Across the transcritical point, the two modes switch signs, and the antihexagon also has one unstable mode. For $\delta<0$ on the other hand, the bottom antihexagon is stable and the top antihexagon and the hexagon are once unstable. Precisely at $\lambda_{N}=\lambda_{N}^{\textit{PF}}$, $\Omega_{1}$ and $\Omega_{2}$ both vanish and one has to go to higher order to see that the antihexagon is the unstable solution. Also, the instability of the anti-hexagon to the mixed mode is not present to this order of the amplitude equation analysis.

VI Exotic Solutions

The aforementioned exact mapping of the ordered pattern equations to a coupled ODE system allows as well for the analytic understanding of a surprising type of pattern not heretofore investigated. If we solve our system at a low value of $k_{c}=10^{-4}$, with $\lambda_{N}=3.5$, $\lambda_{D}=10$ and all other parameters at standard values, we find a hexagonal structure with $N_{A}=52.7>N_{B}=51.0$, while as usual $D_{A}=0.558>D_{B}=0.215$,. We refer to this type of solution as “high-high”, as cells with high Delta also have high Notch. However, we did not find any evidence for this possibility over the range of cis-inhibition parameters proposed by Sprinzak, et al. for typical developmental processes; this is shown in a numerically computed phase diagram (Fig. 8). To better understand this diagram, we show in Appendix F how one can analytically derive the boundary between regular hexagons and high-high solutions in the large $\lambda_{D}$ limit. Specifically, an accurate approximation for the vertical line in the figure is $k_{c}^{*}=k_{t}/2(1+k_{H})$, which equals $0.01$ with our typical choices. This result shows that hexagonal patterns are not dependent on having high cis-inhibition but that the anti-correlation between Delta and Notch cannot be taken for granted in its absence.

VII Initial Value Problem

The existence of stable ordered hexagons leaves open the question of how these patterns can be generated in the noisy biological system with plausible initial conditions Barad et al. (2010); Glass et al. (2016). In particular, it is easy to check numerically that, starting with no pattern for a set of parameters for which the uniform state is linearly unstable, the presence of noise, either in the initial data or in the time evolution, will lead to disordered states with many domain boundaries between hexagon patterns centered on different sublattices. One way out is based on the fact we have shown above that there could exist a parameter range for which there is a subcritical bifurcation to stable hexagons in which case a local perturbation which nucleates the pattern can spread in an ordered manner; this is a standard scenario in many non-living systems Cross and Hohenberg (1993). Crucially, the bifurcation analysis suggests that for biological systems studied to date, there is no significant range of physiologically relevant parameters where propagation would occur into a metastable state. Intuitively, we believe that most biological systems exhibit insufficient parameter control and too high a level of stochasticity for this to be a robust strategy. In some, coupling to additional components could alter the bistability range. There could also be more complex biological mechanisms that, for example, would provide downstream checks that prevent neighboring cells from both developing the same phenotype even if there is some initial defect in the Notch-Delta structure Barad et al. (2011).

A more physics-based possibility is that the system is not put all at once into the unstable state. Rather, we imagine that the system is initially in a regime of parameter space for which the uniform state is stable. Then some external mechanism induces a propagating wave, behind which the parameters are in the unstable region. To exhibit this possibility, we assume that only $\lambda_{D}$ is affected by this wave, and $\lambda_{D}=2$ ahead of the wave and $\lambda_{D}=3.5$ behind the wave. We do not concern ourselves here with the origins or dynamics of this initiation wave, and rather choose a standard tanh waveform, and vary the wave speed $v$. In this regard, our suggestion differs from that of Ref. Pennington and Lubensky (2010); Lubensky et al. (2011), who start with a bistable system with two uniform states - in our proposal, the bistable dynamics is not intrinsically related to the Notch-Delta dynamics. In Fig. 9, we show a pair of simulations of our model augmented by quenched noise. The wave in $\lambda_{D}$ propagates radially outward from an initial point creating an expanding region inside of which the system exhibits Notch-Delta patterning. At large $v$, the parameter shift is essentially instantaneous over a large spatial region and the noise nucleates incommensurate patterns in different parts of the lattice, leading to obvious defects. If the parameter wave is slowed down, the leading edge of the pattern has sufficient time to align itself with the preceding rows before having to itself act as a template for the next radial row. One can do this as well in a planar geometry, as is suggested by some biological data Artavanis-Tsakonas et al. (1995). In some sense, all we have done is transfer the problem to one of creating a bistable system responsible for the parameter dependence. But, we assert that this is relatively easy to accomplish and that decoupling the patterning aspect from the bistable aspect (i.e., the Notch system is not bistable at the physiological parameters) is a robust approach to the elimination of defects.

VIII Summary

In conclusion, we have studied the problem of ordered pattern formation for the realistic Notch-Delta system by mapping it to a 9 degree of freedom ODE system. This mapping facilitates both numerical and analytical progress. We have shown how the presence of a pitchfork bifurcation value $\lambda_{N}^{\textit{PF}}$ close to the physical relevant parameter range organizes the high-Delta centered hexagon pattern as well as the high-Notch centered antihexagon pattern and guarantees that the former is the generic stable structure. Thus, models built on our current understanding of molecular mechanisms do help explain this recurring feature of tissue development. The importance of this is highlighted by the fact that outside the physical range of parameters, alternative correlations between Notch and Delta are possible. Furthermore, we have seen that creating a perfect pattern is a significant challenge in the vast majority of parameter space where the transition from the uniform state to the patterned state is second order. Lastly, we demonstrated that coupling a parameter to an initiation wave could provide a way to meet this challenge.

Acknowledgements.

ET and DAK acknowledge the support of the United States-Israel Binational Science Foundation, Grant no. 2015/619. HL acknowledges the support of the NSF grant no. PHY-1605817. ET acknowledges useful conversations with David Sprinzak.

Appendix A Uniform State Reduction to a Single Equation

When discussing uniform steady-states, it is convenient to reduce the set of three equations to a single nonlinear equation. For simplicity, we will specialize the following discussion to the case of the Hill coefficients $n_{\pm}$ equal to 2. We first introduce the notation

[TABLE]

In terms of these, we have

[TABLE]

where we have also introduced

[TABLE]

Solving for $N_{0}$ and $D_{0}$ yields

[TABLE]

so that $I_{0}$ satisfies

[TABLE]

We next consider the stability condition for a hexagonal perturbation, which in line with the Negrete-Oates model has the form $N_{A}=N_{0}+\delta_{N},N_{B,C}=N_{0}-\frac{1}{2}\delta_{1},D_{A}=D_{0}+\delta_{D},D_{B,C}=D_{0}-\frac{1}{2}\delta_{D},I_{A}=I_{0}+\delta_{I},I_{B,C}=I_{0}-\frac{1}{2}\delta_{I}$. The stability matrix ${\cal M}_{U}$, is then

[TABLE]

The instability sets in when $\textrm{Det}{\cal M}_{U}=0$, i.e., when

[TABLE]

To derive the condition for the bifurcation point to be a pitchfork, we consider a general nonlinear system of steady-state equations, $f_{i}({x_{j}},\mu)=0$, where the functions $f_{i}$ depend on the variables ${x_{j}}$ and the bifurcation parameter $\mu$. The steady-state at the bifurcation point, $\mu_{0}$ is denoted $x_{j,0}$. At a slightly shifted value of $\mu=\mu_{0}+\Delta\mu$, the solution is given by $x_{j}=x_{j,0}+\epsilon x_{j,1}$. Expanding the steady-state equation to first order in $\mu$ and third order in $\\ epsilon$ yields

[TABLE]

where repeated indices are summed over and all derivatives are evaluated at the unperturbed solution. At the critical point, $x_{i,1}$ is given by the right zero-mode eigenvector of the linear stability operator $\frac{\partial f_{i}}{\partial x_{j}}$, $x_{j,1}=x^{R}_{j}$. At a pitchfork bifurcation, the second derivative terms have no projection on the zero-mode, so that

[TABLE]

The right and left eigenvectors are given by

[TABLE]

Substituting this, the condition for a pitchfork bifurcation reads

[TABLE]

For completeness, we note that the condition for a saddle-node bifurcation is in general

[TABLE]

Using the expressions for the eigenvectors above, one can use this to derive a fairly lengthy expression for the saddle-node condition.

Appendix B The Large $\lambda_{D}$ Limit and the Region of Hexagonal Instability

We showed in the main text the region of hexagonal instability in the $k_{c}$, $\lambda_{D}$ plane for various $\lambda_{N}$. One striking feature of this diagram is that the phase boundaries become vertical at two critical values of $k_{c}$, so that $\lambda_{D}^{c}(k_{c})$ diverges at these values. We can derive this analytically by solve the steady-state equations in the large $\lambda_{D}$ limit.

In this limit, $D_{0}\ll 1$ and $N_{0}\ll 1$, we have

[TABLE]

finite, which provides a closed equation for $I_{0}$. Then, the stability condition Eq. (19) reduces to

[TABLE]

Doing the algebra, we find that the critical $I_{0}$ satisfies the equation

[TABLE]

This equation typically has two positive roots for $I_{0}$, which when substituted back into Eq. (17) yields the critical value of $k_{0}/k_{c}=1+k_{t}/k_{c}$. We also see that the critical $k_{t}/k_{c}$ only depends on the lumped parameter $\xi\equiv\lambda_{I}\lambda_{N}/\gamma_{I}$ and $k_{H}$. For large $\xi$, things simplify further and we find the solutions

[TABLE]

corresponding to

[TABLE]

In Fig. B1, we plot the two solutions for $k_{t}/k_{c}$ as a function of $\lambda_{N}$ for the standard parameters $k_{H}=1$, $\lambda_{I}=1$, $\gamma_{I}=1/2$. We see that for $\xi\gtrsim 1$, there are two critical values of $k_{t}/k_{c}$, one near 0 and the other near 2, between which the uniform solution is unstable to a hexagonal pattern. This behavior is manifest in the numerical unstable region plotted in Fig. 4 in the main text.

Appendix C The Large $\lambda_{N}$ Limit and the Region of Hexagonal Instability

Similarly, we can compute the region of hexagonal instability for large $\lambda_{N}$. Here, $N_{0}$ is large and $D_{0}$ is small, with the product being of order 1:

[TABLE]

Then, the stability criterion reads

[TABLE]

which reduces to

[TABLE]

Eliminating $k_{t}/k_{c}$ yields an equation for $I_{0}$ in terms of $k_{H}$ and the lumped parameter $\zeta=\lambda_{D}\lambda_{I}/\gamma_{I}$

[TABLE]

yielding for $I_{0}$,

[TABLE]

Substituting this into Eq. (30) and solving for $\zeta$ yields

[TABLE]

This solution, which specifies the critical value of $\zeta$ (or $\lambda_{D}$ for given $\lambda_{I}$ and $\gamma_{I}$) is physical only if $k_{c}/k_{t}>(2-k_{H})/2(2k_{H}-1)$. This relation, plotted in Fig. 4 for $k_{H}=2/3$ and $k_{H}=4/3$ shows that $\zeta$ decreases from infinty as $k_{t}$ increases from 0, reaches a minimum and then increases, diverging as $k_{t}/k_{c}$ approaches $(2-k_{H})/2(2k_{H}-1)$. The graph is similar for all $1/2\leq k_{H}<2$. For $k_{H}<1/2$ there is no critical $\lambda_{D}$ for large $\lambda_{N}$, whereas for $k_{H}>2$, there is a critical $\lambda_{D}$ for all $k_{t}/k_{c}$.

It should be noted that our results for large $\lambda_{D}$ indicated that that the point of divergence of $\lambda_{D}$ is close to $k_{t}=2k_{c}$ for large $\lambda_{N}$. However, our current large $\lambda_{N}$ calculation indicates that $\lambda_{D}$ diverges at $k_{t}=[(2-k_{H})/2(2k_{H}-1)]k_{c}$. This indicates that the large $\lambda_{D}$ and large $\lambda_{N}$ limits do not commute. This is reflected especially in the right panel of Fig. 4 in the main text, where for $k_{H}=4/3$, we see that for large finite $\lambda_{N}$, $\lambda_{D}$ appears to be diverging at the location predicted by the large $\lambda_{N}$ calculation, but once $\lambda_{D}$ is sufficiently large, it turns back at the curve eventually approaches the $k_{c}=2k_{t}$ line.

Appendix D Multiple Uniform States

We have found through numerical experiments parameter regions that show the existence of multiple uniform states. Here we want to clarify where in parameter space such solutions exist. We start by again assuming large $\lambda_{D}$. We get, following the same scaling as above,

[TABLE]

There are either 1 or 3 solutions of this nonlinear equation. The bifurcation point is where these three solutions collapse to 1. At this point,

[TABLE]

Denoting $Q\equiv\frac{k_{t}\lambda_{N}}{\gamma_{I}k_{0}}$, this reads

[TABLE]

Thus, from the second of these, at the bifurcation, $I_{0}^{*}=\sqrt{3}/3$. The first then implies $k_{H}^{*}=\frac{8\sqrt{3}}{9Q*}$. Plugging these into Eq. 36 gives

[TABLE]

so that

[TABLE]

To see how things look in the vicinity of the bifurcation point, we write $Q=Q^{*}+\delta_{Q}$, $k_{H}=k_{H}^{*}+\delta_{k}$, $I_{0}=I_{0}^{*}+\delta_{I}$. Expanding to third order in $\delta_{I}$, Eq. (36) translates to

[TABLE]

In the last term, we can drop all the $\delta_{k}$, $\delta_{Q}$ pieces as being higher order. Then, for the second term to be the same order as the third, we have to have that $\delta_{k}$, $\delta_{Q}$ are both of order $\delta_{I}^{2}$, and we can drop the last piece of the second term. To make the first term also of the same order, its leading contribution must vanish, so that $\delta_{Q}=-\sqrt{3}/108\delta_{k}+\alpha\delta_{k}^{3/2}$. Then, the equation reduces to

[TABLE]

To have three solutions, we must have $\delta_{k}>0$, so that the derivative of this equation has two solutions. Then, then are solutions in a band around $\delta_{Q}=-\sqrt{3}/108\delta_{k}$, which is negative. Performing the numerical solution for large but finite $\lambda_{D}$ yields the same conclusions, with $k_{H}^{*}$ rising as $\lambda_{D}$ falls. As we take $\lambda_{D}$ smaller and smaller, $k_{H}^{*}$ continues to rise, and eventually diverges at a finite value of $\lambda_{D}$. This is shown in Fig. D1. We will calculate this value momentarily, but the implication is clear: there are only multiple uniform solutions for $k_{H}>8$, which is much larger than physical. This is all the more so since $\lambda_{D}$ is finite, and the multiple uniform solution lower bound on $k_{H}$ is thus in fact much higher.

Our last task is to determine the minimal value, $\lambda_{D}^{*}$ below which multiple uniform solutions are impossible. The numerics show that as $\lambda_{D}$ approaches its critical value, $I_{0}^{*}$ approaches 0 as $k_{H}^{*}$ diverges, with $k_{H}^{*}(I_{0}^{*})^{2}$ remaining finite. In this limit, denoting $x=k_{H}I_{0}^{2}$, we have

[TABLE]

We introduce the notation $\sigma=\frac{\gamma_{D}\gamma_{N}}{k_{0}\lambda_{N}}$, $A=\frac{\lambda_{D}}{\lambda_{N}}+\sigma+1$, $B=4\frac{\lambda_{D}}{\lambda_{N}}$ and rewrite

[TABLE]

At the bifurcation point, we have $R_{0}(I_{0})=0$, $R_{0}^{\prime}(I_{0})=0$, and $R_{0}^{\prime\prime}(I_{0})=0$, where the derivatives are not w.r.t. $x$. Proceeding, we find the solution

[TABLE]

Given $B^{*}$, since

[TABLE]

we have that

[TABLE]

This corresponds to a critical value of $\lambda_{D}/\lambda_{N}$ at which $k_{H}^{*}$ diverges of

[TABLE]

Examining this, we see that there is a maximum of $\lambda_{D}/\lambda_{N}$ for $\sigma=1$. Also, there is a change of behavior for $\sigma>8$. Above this value, $k_{H}^{*}$ remains finite down to $\lambda_{D}=0$, and $\lambda_{D}=0$ it is a decreasing function of $\sigma$, approaching $10+4\sqrt{6}$ from above as $\sigma\to\infty$, In general, as $\sigma$ gets large, $k_{H}^{*}$ is a monotonically decreasing function of $\lambda_{D}/(\lambda_{N}\sigma)$, going between $10+4\sqrt{6}$ at $\lambda_{D}/\sigma\to 0$ and 8 as $\lambda_{D}/\sigma\to\infty$:

[TABLE]

Appendix E Hexagonal Symmetry Steady-State Equation

As discussed in the main text, one can reduce the full system of equations to a coupled ODE system. If we specialize to the case of steady-state patterns with hexagonal symmetry, so that the $B$ and $C$ sublattices are equivalent, we obtain the 6-dimensional system

[TABLE]

Appendix F Hi-Hi Solutions

In the main text, we discussed the fact that the model as presented can support anomalous solutions in which the levels of Notch and Delta are positively correlated in the different sub-lattice sites instead of being anti-correlated. To study this in more detail, we look analytically at the large $\lambda_{D}$ limit and look for the critical curve along which Notch on the $A$ and $B$ sublattice sites are equal. This curve should demarcate the boundary between regular hexagonal solutions and what we are calling Hi-Hi solutions. Studying the numerics (data not shown), we see that $D_{A}$ and $D_{B}$ are large, of order $\lambda_{D}$ and $N_{A}=N_{B}\equiv N$ is small, of order $1/\lambda_{D}$. Writing $N_{A}=N_{B}=N_{0}/\lambda_{D}$, $D_{A}=\lambda_{D}D_{A0}$, $D_{B}=\lambda_{D}D_{B0}$ and specializing to the case of Hill coefficient 2 in the transcriptional/regulation terms, we have to solve the system

[TABLE]

Here we have chosen units such that $s_{0}=1$. This system admits two exact solutions, in terms of $\eta\equiv\frac{\lambda_{N}\lambda_{I}(1+k_{H})}{\gamma_{I}}$:

[TABLE]

The two solutions merge at $\eta=\sqrt{12}$ and do not exist for $\eta<\sqrt{12}$. One key result is that the critical $k_{c}$ is proportional to $k_{t}$. For large $\eta$, which is typical for our parameters where $\gamma=0.1$, things simplify tremendously and we have for the “$-$” branch:

[TABLE]

In fact, for $\eta=\sqrt{12}$ which is its lower limit, we get $k_{c}/k_{t}=(14+5k_{H})/28/(1+k_{H})$, which is only slightly higher than the infinite $\eta$ result.

Going off this border line numerically, we find that the solution to the left of the line is a Hi-Hi solution, whereas to the right is a regular hexagon. The branch starting from the “$-$” border is stable whereas starting from the “$+$” border is unstable. Thus, the two solutions are completely distinct, even in the region to the right of both borders. The region of stable Hi-Hi solutions extends all the way down to $k_{c}=0$. The relevant phase diagram is shown in Fig. 8. We show the line above which stable Hi-Hi solutions exist, and below which stable Hi-Lo hexagons exist, both of which coexist with a stable uniform solution. We also show the saddle-node bifurcation line, below which only the uniform solution is stable.