Numerical simulation of RD patterning inhibition for two-component systems. (A) Violin plots showing percentage change in the mean level of components U and V in illustrated AI and SD RD networks upon inhibition of the response to morphogens U and V in RD simulations (plots for inhibition of production are shown in Fig. S9). (B) Networks showing the two possible configurations of Wnt-Hh AI systems and BMP-Hh SD systems, with components coloured according to the equivalent component in A. Associated with each network is a schematic of the response of Hh upon inhibition of each component in the network, based on the predominant response in the simulations. Solid lines indicate levels after inhibition, with dashed lines representing uninhibited states. Horizontal dotted lines represent an arbitrary detection threshold. The two topologies that have responses to inhibition that replicate the experimental observations (see Fig. 1A) are highlighted with a red box. B, BMP; H, Hh; W, Wnt.

‘Stalactite plot’ and numerical simulation identifying subsets of three-component RD systems and their behaviours under inhibition. (A) Topology atlas for the phase group of the three-component network identified in the parameter search, which is consistent with the spatial pattern of Wnt, BMP and Hh. Topologies that can be recovered with fast-diffusing Hh and slow-diffusing Wnt and BMP are in blue. Out of the 45 strongly connected topologies that are found at the bottom of the stalactites, the 18 that are also consistent with the diffusion constraints are outlined in black. For completeness, not strongly-connected topologies are shown in grey along with their relationship to the strongly connected topologies. (B) Graphs showing the 18 Wnt-BMP-Hh networks divided into four-interaction and five-interaction networks. Within each group, networks are numbered according to their position from left to right in A. (C) Heat map showing the percentage of parameter sets in which the level of Hh increases in response to the inhibition of each component in the network in RD simulations. Topologies are grouped according to hierarchical clustering (see Fig. S10). The ten topologies that have a response that is consistent with the experimental data are highlighted with a red box, as are the network diagrams in B. B, BMP; H, Hh; W, Wnt.

Identification of feedback loops and resulting behaviours under inhibition for three-component RD systems. (A) Illustrative two- and three-component networks showing minimal feedback requirements for RD of one or more components forming a single positive feedback loop (highlighted in magenta), with a negative feedback loop being formed through at least one additional component (highlighted in cyan). (B) Summary constraint table showing the response of components in such a minimal network to the inhibition of a component found in either the positive feedback loop alone (magenta ‘plus’ sign), negative feedback loop alone (cyan ‘minus’ sign) or both (‘plus’ and ‘minus’), depending on which loop they are in and the phase relative to the inhibited component (‘in-phase’ or ‘out-of-phase’). For the response of a component to its own inhibition (‘self’), the inhibition of response (‘res.’) and production (‘pro.’) are shown. Upward pointing arrows indicate an increase in the level of a component, whereas downward arrows indicate a decrease. Two arrows of equal size show when the system is unconstrained. Where opposing large and small arrows are shown, the system behaves according to the large arrows, apart from under certain topologies (see Appendix S1, section 4) in which the reverse is seen for certain components. Although components are not constrained, different components showing the unconstrained response are coupled to one another. (C) Illustrative examples showing how additional components can be integrated into a minimal RD system outside of the ‘core’ RD network by forming external loops. External components and interactions are shown in grey, either forming positive or negative feedback loops with the core positive feedback loop (highlighted in magenta and cyan, respectively). Core RD network outlined with dashes. (D) Summary constraint table showing the response of components in such a network to the inhibition of a component that provides either additional positive feedback (magenta ‘plus’ sign) or additional negative feedback (cyan ‘minus’ sign) to the core positive feedback loop. The table shows how this depends on which loop in the core network they are in [response (Res.); production (Pro.); +, ± or −] and the phase relative to the inhibited component. Symbols as in B. For a subset of topologies in which a component is in both loops see Appendix S1, section 4. B, BMP; H, Hh; W, Wnt.

Integration of two out-of-phase FGF morphogens explains FGF-inhibition effects and allows the prediction of behaviours under inhibition of five-component RD systems. (A) Example of a network where the inhibition of epithelial FGF (eF) and mesenchymal FGF (mF) would be predicted to have opposing effects on the levels of Hh (H) according to analysis of reaction terms, alongside an illustration of how this is determined through the constraints imposed by the different feedback loops in the system [positive feedback loop alone (magenta ‘plus’ sign), negative feedback loop alone (cyan ‘minus’ sign) or both (‘plus’ and ‘minus’), depending on which loop they are in and the phase relative to the inhibited component (‘in-phase’ or ‘out-of-phase’)]. Positive and negative feedback loops are shown in magenta and cyan, respectively. (B) Simulation of single and combined inhibition of the two FGF components (dashed lines indicate uninhibited state). Simulations carried out as detailed in Materials and Methods. u₁ is mF, u₂ is eF and u₃ is H, with a₁₂=−0.019, a₁₃=−0.034, a₂₁=−0.019, a₃₂=−0.022, b₁=0.064, b₂=0.037, b₃=0.068, c₁=0.004, c₂=0.011, c₃=0.039, fmax₁=0.008, fmax₂=0.022, fmax₃=0.072, D₁=1.03, D₂=1.71 and D₃=7.63. Where not specified a_ij=0. Initial conditions drawn from a random distribution, as described in Materials and Methods. (C) Map showing all possible responses of Hh to inhibition of each of the five components for all 39,755 predicted minimal topologies (red, increase in Hh; blue, decrease in Hh). Topologies arranged into 31 different groups of responses, outlined in black. Two sets of responses that constitute 3945 topologies showing the observed responses (Wnt down, BMP up, Hh up, and mFGF and eFGF opposing responses) highlighted in the green box. (D) Map showing the interactions making up the 3945 topologies identified by perturbation analysis (positive interactions in magenta, negative in cyan, no interaction in white). Horizontal black line separates topologies giving different patterns of responses of Hh in response to mFGF and eFGF inhibition (upper group of topologies correspond to the upper group of highlighted topologies in C). B, BMP; H, Hh; W, Wnt.

Periodic gene expression kymographs reveal an early Wnt-eFGF-Hh initiating ‘core’ system with mFGF and BMP integrated later. (A) Kymograph showing the pattern of expression of Shh through time. White arrowheads indicate the approximate onset of Shh expression for each ruga. (B) Plot of mean normalised intensity of Shh staining for each AP position relative to ruga 8. The magenta bar denotes the period of onset of Shh expression, bounded by the minimum in staining intensity and the position at which the staining intensity plateaus (vertical dashed lines). Horizontal dashed line shows level at which Shh intensity plateaus (as determined by the mean intensity of the anterior third of the palate). Shaded area represents 1 s.d. (C) Kymographs showing the pattern of expression of indicated target genes through time (magenta) and their expression relative to Shh (green) for rugae 3, 4 and 5. White arrowheads as in A, and orange arrowheads indicate approximate positions of the change in the expression pattern of each target gene associated with each ruga. Mesenchymal Gli1 and Id1 expression resemble the epithelial patterns (see Fig. S16) (horizontal dark bands in the red channel are stages for which too few specimens were obtained to allow interpolation). (D) Plot of the Spearman’s rank correlation coefficient for the intensity of Shh staining and the marked target gene across all time points for each position relative to ruga 8, indicating when, relative to the onset of Shh expression, the spatial pattern for each target gene emerges. Horizontal dashed lines represent maximal correlation coefficient, calculated over the anterior third of the palate (Fig. S15D) and half this value. Vertical dashed lines represent AP position where this level is first reached. Lower dashed lines for Gli1 show where the half maximal level is reached for the opposite correlation (i.e. out-of-phase rather than in-phase). Shaded area represents 95% confidence interval from bootstrapping (see Materials and Methods). The magenta bar represents the period of onset of Shh expression, as determined in B. (E) Sequence whereby different targets come into (upward pointing arrowheads) or out-of (downward pointing arrowheads) phase with Shh relative to the distance from ruga 8, based on the correlation analysis in D. The magenta bar represents the period of onset of Shh expression, as determined in B. a.u., arbitrary units.