Perturbation analysis of a multi-morphogen Turing reaction-diffusion stripe patterning system reveals key regulatory interactions

ABSTRACT Periodic patterning is widespread in development and can be modelled by reaction-diffusion (RD) processes. However, minimal two-component RD descriptions are vastly simpler than the multi-molecular events that actually occur and are often hard to relate to real interactions measured experimentally. Addressing these issues, we investigated the periodic striped patterning of the rugae (transverse ridges) in the mammalian oral palate, focusing on multiple previously implicated pathways: FGF, Hh, Wnt and BMP. For each, we experimentally identified spatial patterns of activity and distinct responses of the system to inhibition. Through numerical and analytical approaches, we were able to constrain substantially the number of network structures consistent with the data. Determination of the dynamics of pattern appearance further revealed its initiation by ‘activators’ FGF and Wnt, and ‘inhibitor’ Hh, whereas BMP and mesenchyme-specific-FGF signalling were incorporated once stripes were formed. This further limited the number of possible networks. Experimental constraint thus limited the number of possible minimal networks to 154, just 0.004% of the number of possible diffusion-driven instability networks. Together, these studies articulate the principles of multi-morphogen RD patterning and demonstrate the utility of perturbation analysis for constraining RD systems. This article has an associated ‘The people behind the papers’ interview.


Supplementary Figures and Tables for Economou et al. " Perturbation analysis of a multimorphogen Turing Reaction-Diffusion stripe patterning system reveals key regulatory interactions"
This document contains • Supplementary Figures S1-S16 (pages 1-10) • Supplementary Tables 1 and 2 (page 11) • Supplementary Notes 1-4 (pages 12-45)   ), red points relate to rugae determined to have fused with another ruga in the inhibitor treated explants, with their position in plotted against the position of that ruga in the inhibitor treated explant.Yellow and green points relate to rugae in the inhibitor treated and untreated explants respectively, where the equivalent ruga could not be identified in the stage matched explant (respectively the untreated or inhibitor treated).For measures of stripe width c) and intensity d) points in red denote posterior rugae (see Methods).e) Percentage of rugae for each treatment lying within one median displacement from the line of equal values (grey), or above (magenta) and below (cyan), for each of the three measures in b-d).
All plots depict quantifications from four pairs of explants.Drug treatments were carried out at the doses specified in the Methods, plus an additional dose, with points pertaining to lower doses shown as circles and higher dose as triangles in all plots.Doses used are: cyclopamine 20 µM and 80 µM, SU-5402 40 µM and 80 µM, IWP-2 10 µM and 50 µM, dorsomorphin 10 µM and 50 µM.Topologies are group as according the the hierarchical clustering in figure S7. b) Schematics showing the ten minimal 5-interaction topologies networks recovered from the parameter search.On left, the net feedback between Wnt and BMP is marked, with magenta for a positive feedback, or cyan for a negative feedback.(It should be noted that for some topologies (such as iv, vii, iii and ix) the feedback is indirect with both the interaction from Wnt to BMP and BMP to Wnt pass through Hh.For example, in iv, while there is no single positive feedback loop, the net interaction between Wnt and BMP is a mutual inhibition and therefore a net positive feedback.)On the right, the presence of a negative feedback between Hh and the autoactivating component (Wnt or BMP) is marked in cyan.Topologies are group as according the the hierarchical clustering in figure S7.Wnt in blue, BMP in green and Hh in red.Coupling between components (see supplementary note 4 and figure 4) is not shown.The ten topologies for which the predicted responses are consistent with experimental observations are marked in red.These are the same topologies as recovered by the RD analysis (see figure 3).The conditions for a diffusion driven instability (DDI) have been extensively analysed.For completeness, we present general conditions for a DDI, as well as specific criteria for assessing whether specific parameterisations of two-and three-component RD systems will give a DDI.

1.1: General conditions for a DDI
We considered a general reaction-diffusion system of the form where u is a vector of N reactant concentrations and D is a diagonal N×N matrix of diffusion coefficients where  ≥ 2, such that Following the well-established linear stability approach to deriving conditions for DDI, as outlined Murray (2003) (see also White andGilligan 1998, Marcon et al. 2016) (N.B. References below), we derive the stability matrix S, for which where q is the wave number of a spatially-periodic perturbation, and J is the Jacobian matrix of the reaction system, where the elements Jij are the first partial derivatives of the components of f, evaluated at a spatially-uniform steady state, such that For diffusion driven instability, the system must be stable without diffusion and unstable with diffusion.Thus, when q=0, all solutions  of the dispersion relation must have negative real part, and there must exist positive values of q 2 for which there is at least one solution  with positive real part.
We will now consider the specific cases of two-and three-component RD systems.

1.2: Two-component system
For a two-component system of the form above, the dispersion relationship takes the form Following the approach of Murray (2003) gives the conditions for a DDI.Stability in the absence of diffusion requires Development: doi:10.1242/dev.190553:Supplementary information Development • Supplementary information and  ""  !! −  "!  !" > 0, and instability with diffusion requires As is well established (see Murray 2003), if J11>0 and J22<0, the two components will be inphase when J12<0 and J21>0, and out-of-phase when J12>0 and J21<0.

1.3: Three-component system
For a three-component system of the form above, the dispersal relationship takes the form Following the approach of White and Gilligan (1998), we considered the Routh-Hurwitz criteria.Stability without diffusion requires and instability with diffusion requires Development: doi:10.1242/dev.190553: Supplementary information

Development • Supplementary information
As detailed by White and Gilligan (1998), the conditions for instability can only be satisfied if Both of the functions in the inequalities are cubic functions of the form where a, b, c and d are constant coefficients.A DDI will occur if there exists a minimum turning point  12 !, for which  12 !> 0, and ( 12 !) < 0 (see White and Gilligan (1998) for details).
RD systems with more than two components can destabilise either as monotonically growing spatially periodic patterns (as for two component systems) or temporally oscillating patterns.
When the system is destabilized with real positive l, the spatial patterns grow monotonically, while if l has non-zero imaginary parts, they will oscillate temporally.White and Gilligan (1998) provide a detailed analysis of conditions for monotonically growing and oscillating patterns.From their analysis, a parameterization will be produce monotonically growing The function in the final condition is quadratic in q 2 , with a minimum at  34# !, such that For each parameterization, a phase relationship was calculated as the eigenvector of the largest positive eigenvalue of the matrix calculated at the turning point  12 !for which  .( ! ) < 0. While this is not necessarily the phase pattern with which the pattern with the fastest growing wavelength will grow, it does give a phase pattern that is consistent with the network topology defined by the reaction matrix.

2: Perturbation predictions from RD reaction terms
Our analysis of the response of RD systems to perturbations suggested that there is a relationship between the feedback loops within the RD system, and the nature of the response of the system to perturbation.Specifically, whether a component increases or decreases its level in response to the inhibition of another component depends on whether the components are in a positive or negative feedback loop.As the feedback loops are encoded by the reaction terms of the system, we considered what role the reaction terms could play in this behaviour.
Preliminary simulations suggested that the response of an established spatially periodic steady state solution of an RD system to perturbation is predicted by the response of the reaction term equilibrium, and whether it increases or decreases in response to equivalent perturbations.We therefore decided to investigate further the factors that govern the responses of the reaction terms.

2.1: Shape of N-component response
We considered the response of a system of N components, using a system of linear ODEs equivalent to that used for the full reaction-diffusion system we used before, where The parameters aij represent the weight of the interaction from component uj to component ui, where a positive term indicates an activation and a negative term an inhibition,  $ is a background production term and ci is the linear degradation rate of component ui.
The position of the equilibrium is given by setting = 0 for all i, and as in (Desoer 1960), rearranging this gives where the matrix A is given by and C[] %$ is the cofactor of the element of the j th row and i th column of A.
We considered the effect of inhibiting a component on the equilibrium levels of the system.
Specifically, we looked at the effect of inhibiting component u1 on both the levels of u1 itself, and the levels of other components by looking at the effect on a component u2.As for our simulations of full RD systems we performed two forms of inhibition, namely inhibiting the response of a component to the levels of u1 (analogous to inhibiting a "u1 receptor") and inhibiting the production of u1.

2.1.1: Inhibition of response to u1
We first consider the effect of perturbing the response to u1 by multiplying the weight of all interactions from component u1 (to itself and all other components) by a factor of 1 −  (i.e. with no inhibition when  = 0 and complete inhibition when  = 1).At equilibrium, this gave a system of the form and It should be noted that for the response of u1 to its own inhibition, only the weight of a11 is reduced, as c1 represents a linear degradation term and would therefore be unaffected directly by an inhibitor blocking a receptor.
Rearranging these equations gave equilibrium levels of ui upon inhibition of strength a, where and for i=2,…,N.Matrix C, which does not include the degradation term c1, is given by and C[] %$ is the cofactor of the element of the j th row and i th column of C.

Response of u1:
The response of u1 to its own inhibition is governed by Equation (2.1), which is a function of a with horizontal and vertical asymptotes at respectively.The intercept of the u1-axis lies at Development : doi:10.1242/dev.190553:Supplementary information

Development • Supplementary information
which is the unperturbed level of u1 (i.e. when a = 0).As the unperturbed level of u1 has a positive value, the function in Equation (2.1) can take two forms, depending on whether the vertical asymptote corresponds to a positive or negative value of a (see Figure SN1 a).
Therefore, u1 will increase with its own inhibition if and will decrease if

Development • Supplementary information
Response of u2: The response of u2 to inhibition of u1 is governed by Equation (2.2), which is also a function of a with horizontal and vertical asymptotes, which lie at and respectively.The intercept of the u2-axis lies at which is the unperturbed state (when a = 0) and has a positive value, while the intercept of the a-axis is at . For the remaining three (Figure SN1 b iv-vi), u2 decreases with the inhibition of u1.From the positions of the horizontal asymptote relative to the unperturbed state at  = 0, for conditions i, ii, iv and v while for conditions iii and vi However, from the position of the vertical asymptote, for i, ii and vi while for iii, iv and v Therefore, in conditions i, ii and iii, where u2 increases on the inhibition of u1, while in conditions iv, v and vi, where u2 decreases on the inhibition of u1,

2.1.2: Inhibition of production of u1
We next considered the effect of inhibiting the production of u1 by reducing the weights of all the coefficients affecting the production of u1 (including the constant term, but excluding the degradation coefficient) by a factor of 1-a.At equilibrium, this gives a system of the form for i=2,…,N.Rearranging these equations gives the equilibrium levels of ui upon inhibition of strength a, where (2.9) and for i=2,…,N.As the equation giving the equilibrium value of ui where  ≠1 is the same as that for the inhibition of the response (Equation (2.2)), the shape of the response of u2 to an inhibition of strength a is the same for inhibiting either the production of, or the response to, u1.

Behaviour of u1:
The behaviour of the equilibrium value of u1 is again governed by a function of a (Equation (2.9)) with horizontal and vertical asymptotes, which now lie at and respectively.The intercept of the u1-axis lies at which is the unperturbed state and has a positive value, while the intercept of the a-axis is at Depending on the relative positions of the asymptotes and intercepts, this function can take three forms (Figure SN1 c), for two of which (Figure SN1 c i, iii), the equilibrium value of u1 decreases with its inhibition, while for the remaining one (Figure SN1 c ii), the equilibrium value of u1 increases with its inhibition.Therefore, the equilibrium value of u1 will increase on its inhibition if While the equilibrium value of u1 will decrease on its inhibition if

2.2: Comparison with a numerically modelled RD system
Using this framework, we could compare the responses to inhibition for the full 2-and threecomponent RD systems with the behaviour of the reaction terms.We took the parameterisations of the reaction terms used in the full RD simulations (Figure 3 and Supplementary Figure S11) and inferred the direction of response by calculating the relative position of the asymptotes and intercepts.Comparison of the direction of the response to an inhibition for the reaction terms alone and for the full RD system (Figure SN2) showed that the pattern of the responses is very similar.This suggests that the response to an inhibition is predominantly dictated by the response of the reaction terms.

Development • Supplementary information
What about the small proportion of parameterisations for which the response of the equilibrium values differs between the full RD system and the analysis of the reaction terms alone?Interestingly, repeating the simulations, but reducing the strength of the inhibition by a factor of 0.25 or 0.25 2 , successively increased the discrepancy (Figure SN2).As decreasing the strength of an inhibition decreases the magnitude of the response for the reaction terms, this suggests that the response of the reaction terms will dictate the response of the full RD system provided that the responses are sufficiently large.As the only difference between the reaction term systems and the full RD systems is diffusion, this suggests that if the response to an inhibition of the reaction terms is not sufficiently large, there can be a significant effect of diffusion, which in some cases can change the direction of response generated by the reaction terms alone.
This is likely to be generally true for RD systems based on the following argument.Large perturbations in reaction terms will produce large shifts in the mean value of the periodic pattern while changes in the diffusion terms will flatten or sharpen the peaks and troughs.For example, an increase in the reaction term will increase the value of ui in both peaks and troughs.

Development • Supplementary information
For equilibrium to be restored, the diffusion term  !" !# " "$ ! has to decrease.Since at the peaks, "$ ! is negative, this means it has to become more negative, i.e. the peaks must become sharper.Conversely, the troughs must become flatter (the positive second derivative becomes less positive).These changes in the curvature of the peaks and troughs do not necessarily have a net effect on the average levels that are the same or opposite to the reaction term change.Our simulations, however, say that their effects do not increase at the same rate as those of the reaction terms.The simulations suggest that their role is more pronounced for small reaction term changes although we have not been able to demonstrate this analytically.Development: doi:10.1242/dev.190553

3: A graphical interpretation of the responses to perturbation
We wanted to understand how the position of a component in a network affects the response of the system to its inhibition.In particular, we wanted to understand how the response depends on the type of feedback loop the component is found in.In a recent study Marcon et al. (2016) used such a graphical approach to understand the conditions for DDI, in terms of stabilizing and destabilizing contributions of different feedback loops.This was based on an interpretation of the determinant as a combination of the cycles defined by the reaction matrix.
As the conditions determining the response to perturbation (see Section 2) are in part dependent on the determinants of matrices A and C, we took a similar approach.
In brief, a matrix M (such as the reaction matrices A and C, or in the approach of Marcon et al. the Jacobian matrix J) has associated graph Gr[M] which is a labelled, weighted, directed graph of N nodes with edges from node j to node i of weight aij, where  $% ≠ 0 are the elements of the matrix M.This is the graphical representation of the reaction network that we have been using throughout to depict the signs of the interactions as a network topology.Such a graph contains cycles, where a path can be traced from any node, through a series of other nodes, back to itself (or directly onto itself without passing through any other nodes).A cycle represents a positive feedback loop if the product of the weights of the edges is positive, and a negative feedback loop if the product is negative.
For any graph Gr[M], a subgraph can be defined containing all the nodes but only a subset of the edges of Gr[M], such that each node receives an input from only one other node and also outputs to only one other node (i.e. each node has both an in-degree and out-degree of 1).Such a subgraph is called a linear spanning subgraph (L-subgraph, l), and is made up of a set of 1 to N cycles through all nodes.The weight of an L-subgraph w(l) is given by the product of the weights of these cycles, with the weight of each cycle itself being the product of its constituent edges.The determinant of the N×N matrix M is given by the sum of the weights of all the Lsubgraphs l of the graph Gr[M], with signs as below

3.1: Response of u1
As demonstrated in Section 2, the response of u1 to its own inhibition depends on the sign and magnitude of the ratio of det[C] to det [A].The sign is determined by the criteria for a DDI (see Section 1).To understand the response of u1 to its own inhibition in graphical terms, we therefore considered these determinants in terms of the feedback loops (or cycles) through u1 in matrix C. From Equation (3.1), each summand in det [] can be written in the form (−1) # (−1) >(?) (), where l is a linear spanning subgraph of Gr [C].As all nodes in an l-graph have an in-degree and an out-degree of 1, for each summand in Equation (3.1) there is a single cycle containing node u1, which comprises a set of n nodes, where 1 ≤  ≤ , with  E a sequence of integers {i1,i2,…,in} which defines the sequence of nodes ui through which the cycles passes, with each such sequence starting at i1=1.If n=1, the cycle forms an interaction from u1 directly onto itself, while if n=N, the cycle passes through all nodes. F ( is the set of all possible  E , and ( E ) is the weight of the cycle defined by  E (positive for a positive feedback loop, negative for a negative feedback loop).For each such cycle, where n<N, the nodes excluded from the cycle form an L-graph ln consisting of 1 to N-n cycles between the remaining nodes, where ( E ) is the number of remaining cycles (with ( E ) = () − 1 as one cycle has been removed).E ) is the product of the weights of the remaining cycles such that () = ( E )( E ), and  ?( is the set of all possible L-graphs  E for a given  E .

𝑤(𝑙
Therefore, for any term in the determinant, the loop through u1 can be factored out to give Gathering together determinant terms containing the same cycle through u1, the determinant can be rewritten Development: doi:10.1242/dev.190553 where (Ẽ) is the submatrix of C excluding the columns and rows defined in  E (for n=N, det [(Ẽ)] = 1, as it is the determinant of a 0×0 matrix).

3.1.1: Inhibiting the response to u1
As demonstrated in Section 2, the behaviour of u1 on inhibiting the response to u1 is dependent on whether the ratio HIJ [𝑪] is positive or negative (see Equations (2.3) and (2.4)).Therefore, u1 will increase in response to its own inhibition if and will decrease if where each summand contains a different cycle through u1, each multiplied by the determinant of the submatrix made of components excluded from the cycle.The response of u1 to the inhibition of its response therefore depends on the natures of the feedback loops through u1, and the stability of the components excluded from the loops.

3.1.2: Inhibiting the production of u1
For inhibiting the production of u1, we demonstrated in Section 2 that the response of u1 again depends on det [] and det [] (see inequalities 2.10 -2.12).First considering when N is even, As matrices A and C only differ by the presence or absence of the -c1 in the first matrix term, then where (") is the submatrix of C, excluding the first row and column.It therefore follows that u1 will increase following its own inhibition if and u1 will decrease following its own inhibition if Given that  " > 0, it follows that u1 will increase when its production is inhibited if and will decrease when By the same reasoning, these conditions also hold if N is odd.The response of u1 to the inhibition of its production therefore depends on the stability of the submatrix made from the remaining components.

3.2: Response of u2:
As demonstrated in Section 2 (Equations (2.5) and (2.6)), the response of component u2 in an N-component system to the inhibition of u1 is in part determined by the relative signs and sizes of the determinants det[A] and det [C].The response to inhibition also depends on the terms which are of the same form as a determinant, but with the second column of matrix A and C respectively replaced by the vector of background production terms (b1,b2,…,bN).Therefore, each term in these sums can similarly be interpreted graphically, as a path from an external node ub to component u2, with all components not involved in the path making up a set of cycles with each node having an in-degree and an out-degree of 1.
Therefore, for each summand, every node has a single input and a single output, except ub, which only has a single output and no input, and u2, which only has a single input and no output.
Conditions (2.5) and (2.6) can be rearranged to give and respectively, where C[(̃")] (%-")" is the cofactor of the term aj2 from matrix C in the submatrix (̃") (as this submatrix excludes the first row and column of C, this will be the cofactor of the first column and the (j-1)th row).
The only term in conditions (3.7) and (3.8) whose sign is not constrained is . From Equations (2.1) and (2.9), which describe the levels of u1 under inhibition, for the unperturbed state of u1 to be positive, the term where pn is a sequence of integers {i1,i2,…,in}, where 2 ≤  ≤ , starting with i1=1 and ending with in=2, which defines the series of nodes ui through which the path from u1 to u2 passes,  M ( is the set of all possible such paths, and ( E ) the weight of the path pn.The response of u2 to the inhibition of u1 therefore depends on the weight of the paths between the two nodes, and the stability of the components excluded from the paths.
Development : doi:10.1242/dev.190553, (4.1) which will always have a positive sign.Therefore, unsurprisingly, the condition for instability when there is diffusion cannot be satisfied in this simple system, since  # ( ! ) < 0.
In order to satisfy the conditions for DDI there must be at least one term in  # ( ! ) with a negative sign.From Equation (4.1) this can be achieved by simply adding a single autoactivation term for a component.For example, if u1 directly activates itself, providing that this auto-activation is strong enough to outweigh the degradation term for u1 (in terms of matrix A,  "" −  " > 0), a single term on the diagonal of the Jacobian matrix J11 will go from being negative to positive, resulting in the sign of any term det[( 9 )] for which  9 contains u1 also changing sign.Alternatively, this can also be achieved by the addition of a single positive feedback loop through the set of r nodes defined by the sequence  A of weight ( A ).For any where  _ A is a complementary sequence of integers to  A such that  A ∩  _ A = 0 and  A ∪  _ A =  9 .
Again, provided that the weight of the positive feedback loop is greater than for the degradation terms, det[( 9 )] will change sign.
However, as the sign of any submatrix containing the positive feedback loop will change sign, it follows that with the addition of a single destabilising positive feedback loop, (−1) # det[] < 0.
This will no longer satisfy the condition for stability in the absence of diffusion  # (0) > 0. An additional positive summand in  # (0) is therefore needed to stabilize the system.By the above reasoning, this can most simply be achieved by the addition of a single negative feedback between a set of components

Development • Supplementary information
those for which  A is also a different subset).If the weight of this loop is sufficiently large, then it can overcome the effect of the positive feedback loop, and satisfy  # (0) > 0.
Moreover, as there will still exist  9 for which  A is a subset, but not  O , there will be terms  !(5-9) (−1) 9 det[( 9 )] det[( _ 9 )] which will remain negative, even after the addition of a negative feedback loop to the system.Therefore, with the appropriate weighting of diffusion coefficients, there will exist conditions where  5 ( ! ) < 0 when k=N, but not for other k, while maintaining  # (0) > 0. As such  _ 9 will by definition exclude components of the positive feedback loop through the nodes in  A , to achieve such a weighing, diffusion coefficients among components in the negative will be larger than those in the positive feedback loop.This gives an intuitive understanding for then general idea of local auto-activation and long-range inhibition in RD systems.In summary, a system containing single positive feedback loop coupled to a negative feedback loop is, therefore, sufficient to satisfy the conditions for a DDI for any size system, given appropriate strengths of interactions and diffusion coefficients.
Considering the minimal requirements for DDI in terms of the number of interactions between components in the interaction matrix A (including auto-activation terms aii, but not degradation terms ci), in any such system a positive feedback loop through r components will add r interactions to the system, while the negative feedback loop though the remaining N-r components will add a further N-r+1 terms, resulting in a total of N+1 interactions.For any strongly-connected network, it is clear than fewer than N+1 interactions will not allow sufficient feedback loops to be made to satisfy the criteria for DDI, while other networks made of N+1 interactions (e.g. two positive or two negative feedback loops) cannot satisfy the conditions for DDI.While there are other ways of satisfying these conditions, they will require additional feedback loops, and therefore additional interaction.Therefore, in a system of N components, the combination of a single positive feedback loop through one to N-1 components, and a negative feedback loop through the remaining components is the minimal requirement for a DDI.

4.1.2: The phase of minimal DDI systems
Having established the architecture of minimal topologies, we next sought to establish the relative phase of the components for each of these topologies.This is dictated by the signs of Development: doi:10.1242/dev.190553:Supplementary information

Development • Supplementary information
the components of the eigenvector corresponding to the largest positive eigenvalue of the Jacobian matrix: components with the same sign will be in the same phase.From the derivation of conditions for DDI in Section 1, we have the relationship For all minimal architectures, there is only one component that receives more than one input negative self-interactions).For all remaining components, the rows in the matrix will consist of only two non-zero terms: the degradation term on the diagonal,  $$ which is therefore negative, and an additional input term from another component  $% , whose sign depends on the nature of the input.Therefore, each of these rows will take the form By definition, for a DDI l>0, meaning that  $$ −  $  !−  < 0. Consequently, ui and uj will have the same sign, and therefore be in-phase with each other if Jij is positive, while they will have opposing signs and be out-of-phase for negative Jij.Therefore, the relative phase of successive components in the positive and negative feedback loops can be inferred from the signs of interaction between them, with the exception of the interaction from a component in the negative feedback loop feeding into the positive feedback loop.Moreover, two components will be in-phase if the weight of the path between them is positive, and out-of-phase if negative, unless the path passes from the negative feedback loop into the positive feedback loop, in which case the relationship will be reversed.This is in keeping with recent finding on the relationship between topology and phase by Diego et al.

4.2: Responses to perturbation of minimal topologies
Having established minimal conditions for DDI, we next asked how these minimal systems respond to the inhibition of the different components.In particular, as these systems consist only of a single positive and a single negative feedback loop, we considered how the response to inhibition depends on the placement of a component relative to these two loops.Development: doi:10.1242/dev.190553:Supplementary information Development • Supplementary information

4.2.1: Inhibition of response
For the inhibition of response to u1, in Section 3.1.1,conditions (3.3) and (3.4), it was shown that the response is dependent on the sign of where each summand corresponds to a different feedback loop through u1.For minimal topologies, the sign of each summand depends on two features: whether the feedback loop supported by u1 is a positive or negative feedback loop (given by ( E )), and whether the nodes excluded from this loop support a positive feedback loop (given by det[(Ẽ)]).For each summand, ( E ) will be negative if it represents a negative feedback loop and positive for a positive feedback loop.However, the relationship between the sign of det [(̃E)], and the feedback loops contained depends on the number of nodes in the system described by submatrix (Ẽ) (see Section 1).Considering all possible combinations of odd and even numbers of components in matrix A and submatrix (Ẽ), it can be seen that (−1) E det [(̃E)] /det [] will always be negative if submatrix (̃E) contains a single positive feedback loop, and will be positive if it only contains negative feedback loops.
As minimal systems only contain a single positive feedback loop and a single negative feedback loop, submatrix (Ẽ) will always contain only negative feedback loops as the two loops must share at least one component.Therefore, each summand in Equation (4.2) will be positive if it represents a positive feedback loop and negative if it represents a negative feedback loop.For certain placements in a minimal network, u1 will only fall in one feedback loop, and there will be a single non-zero summand in Equation (4.2).Therefore, if u1 is only found in the positive feedback loop it will decrease on its inhibition, and if it is found only in the negative feedback loop it will increase on its inhibition.However, for some placements in a minimal network u1 will be in both positive and negative feedback loops.In these instances, there will be two summands in Equation (4.2), one positive and the other negative.The response of the system will depend on the relative magnitude of each term, with both an increase and a decrease in u1 possible.These constraints are summarised in Figure 4b.

Development • Supplementary information
As shown in Section 3.1.2,conditions (3.5) and (3.6), the response of u1 to the inhibition of its own production is determined by the sign of det [(̃")] relative to det [𝑨].This in turn, depends on whether the system described by the submatrix excluding u1 contains the positive feedback loop of the minimal network or not.If u1 is excluded from the positive feedback loop, the submatrix determinant det [(")] will therefore contain the single positive feedback loop in the system.As matrix (̃") by definition has one fewer component than matrix A, as detailed in equation 3.1, det [] and det [(")] will have the same sign.However, if u1 is included from the positive feedback loop, regardless of whether it is also part of the negative feedback loop, det [(")] will only represents negative feedback loops, and det [] and det [(")] will have opposite signs.Therefore, if u1 is excluded from the positive feedback loop, it will increase upon inhibition its production, while if it is part of the positive feedback loop (whether or not it is also in the negative feedback loop) it will decrease if its production is inhibited.These constraints are summarized in Figure 4b.

4.2.3: Response of u2
Conditions for the response of component u2 to inhibition of component u1 were defined in supplemental note 3, conditions 3.7 and 3.8.Writing these in terms of paths between u1 and u2 (see equation 3.9) the response of u2 depends on the sign of where each summand contains a different path from u1 to u2.Written in this form, it is apparent that the response of u2 to the inhibition of u1 is again dependent on two topological features of the minimal system: the weight of the path from component u1 to u2 (given by the term ( E )) and whether the nodes excluded from this path support a positive feedback loop (given by det[(Ẽ)]).
If there only exists a single path from u1 to u2, there will only be a single non-zero term in Equation (4.3).As shown in the above analysis of the phase relationship (Section 4.1.2),the sign of the path ( 5 ) will be positive if u1 and u2 are in-phase and negative if out-of-phase.The only exception is when u1 is excluded from the positive feedback loop but the path to u2 passes into )(−1) E det [(Ẽ)] will be positive if det [(Ẽ)] contains the destabilising positive feedback loop, and negative otherwise.
Therefore, considering all possible placements of u1 and u2 relative to the two loops of the minimal system it can be seen that if u1 is in the positive feedback loop (and therefore the positive feedback loop is not found in submatrix (̃E)), u2 will decrease if it is in-phase with u1, and increase if out-of-phase.However, if u1 is excluded from the positive feedback loop, u2 will increase if it is in-phase with u1, and decrease if out-of-phase, as either the path to u2 will pass from the negative feedback loop into the positive feedback loop, or this path will exclude components from the positive feedback loop.These constraints are summarized in Figure 4b.
It should be noted that in certain instances when u1 and u2 lie in both the positive and negative feedback loop, two paths can be traced from u1 to u2: one consisting of interactions only found in the positive feedback loop, and the other comprising ones only found in the negative feedback loop.There will therefore be two summands in Equation ( 4.3), one associated with each path, and they will be of opposing signs.Therefore, the response of u2 to the inhibition of u1 will depend on the relative weight of each path, and it can either increase or decrease.
It should be noted that while the response of u2 to the inhibition of u1 is not constrained by topology in these instances, it is coupled with the responses of other components in the system.
Factoring out all coefficients shared between the two summands leaves ‰−1) P-Q(P -) ( P -Š det[(̃P-)] + ‰−1) P-Q(P . )( P .Š det[(P.)],(4.4) where matrix W is the submatrix ( P ) comprising the w components supporting the two paths defined by  P .The terms ( P -) and ( P .) give the weights of the paths through the positive and negative loops, respectively, ( R ) and ( -) the numbers of components in the respective paths, and (̃P-) and (̃P-) the submatrices of W comprising the components excluded from the respective paths.This same sum is also left when considering the response ) ( P .Š det[(P.)]•, then u1 will decrease on its own inhibition, as will any in-phase components, with out-of-phase components increasing.However, when the inequality is reversed and u1 increases on its own inhibition, the response of any component to which the two paths can be traced will be coupled, increasing if it is in-phase, and decreasing if out-of-phase.

4.3: Adding nodes to minimal topologies
In the above analysis of the criteria for a minimal topology (Section 4.1.1)we only considered strongly-connected networks, where the two feedback loops comprise all N components.For systems where the two loops only pass through a subset of components, if the remaining components respond passively (i.e.do not feed back into the system) each component will only add a single interaction.Like the strongly-connected networks these networks also have two loops and N+1 interactions.
These passive components can be wired back into the minimal RD 'core' by adding interactions from the passive components back to the core.If the weights of any additional interactions, or the weights of additional loops generated by them, are not excessively large, the criteria for DDI will still be met, and the phase relationship of the related not strongly-connected minimal topology will also be maintained.While the resulting topologies are not strictly minimal in terms of the number of loops or interactions, they comprise additional strongly-connected networks which satisfy the conditions for DDI.Moreover, as the removal of an interaction will either result in the system either no longer satisfying the criteria for DDI or no longer being strongly-connected, we will also consider such topologies as an additional form of minimal topology.
We will now consider the effect of inhibition on networks with additional components wired in outside of the RD 'core'.We first will consider the case where there is a single component u1, external to the core, and the effect of inhibiting this component on its own levels, before

Development • Supplementary information
considering the effect of inhibiting u1 on core components.We then consider the effect of external interactions on the core.

4.3.1: Effect on self
The addition of a single component external to the RD core will generally result in the addition of a single positive feedback loop to the system (we will discuss below the subset of cases where the addition of a single component results in the addition of two loops).

Inhibition of response
As for the minimal topologies (Section 4.2.1) when the response to u1 is inhibited, the behaviour of u1 depends on the sign of the sum in Equation (4.2).This sum will consist of a single term, and as before, the sign of this term depends in part on whether or not the submatrix (Ẽ) contains the destabilising core positive feedback loop, and also whether this loop supports positive or negative feedback.We will ask how the wiring of this loop relative to the core positive feedback loop, along with the weighting of the loop, determines the response to inhibition.
We will first consider when the feedback loop through u1 includes components from the core positive feedback loop.In this case submatrix (̃E) contains only negative feedback loops and (−1) E det [(Ẽ)] /det [] will be positive.If u1 itself supports a positive feedback loop, condition (3.4) will be satisfied and u1 will decrease on its own inhibition.Conversely, if u1 supports a negative feedback loop, condition (3.3) will be satisfied and u1 will increase on its own inhibition.These relationships are reversed if the loop through u1 does not pass through any components in the core positive feedback loop.If this is the case, submatrix (̃E) will contain the core positive feedback loops and (−1) E det [(Ẽ)] /det [] will be negative (otherwise conditions for DDI cannot be met -see Section 1).In this instance, if u1 supports a positive feedback loop, condition (3.3) will now be satisfied and u1 will increase on its own inhibition, while if u1 supports a negative feedback loop, condition (3.4) will be satisfied and u1 will decrease on its own inhibition.
Interestingly, the response of u1 to its inhibition does not correlate directly with the feedback loop u1 supports, but rather with the net feedback provided to the core positive feedback loop while if it provides a net negative feedback, it will increase on its own inhibition.These constraints are summarized in Figure 4d.

Inhibition of production
As demonstrated in Section 3, the response to the inhibition of production is dependent on the sign of det [(")] relative to det [] (condition (3.5) or (3.6)).As det[] and det[] only differ by the presence of the degradation coefficient c1 in the first matrix term it follows that In the above analysis of the inhibition of response, we demonstrated how network topology (i.e. the sign of the loop containing u1, and how this feeds into the core RD network) determines the sign of (−1) E ( E )det [(̃E)] /det [𝑨].From these arguments, it follows that if u1 supports a net positive feedback, − " det[(")]/ det [] must be negative.As − " is negative, condition (3.6) will be satisfied and u1 will decrease.However, if u1 supports a net negative feedback, the sign of − " det[(")]/ det [] is not constrained and therefore both conditions (3.5) and (3.6) can be satisfied.Therefore, u1 will decrease on the inhibition of its production if it provides a net positive feedback to the core positive feedback loop, but depending on the strength of the different feedback loops, u1 will either increase or decrease when its production is inhibited if it provides a net positive feedback.These constraints are summarized in Figure 4d.
It should be noted that in some instances where u1 feeds into and out of the core RD network at components which are found in both core positive and negative feedback loops, it is possible that a second loop supported by u1, through components only found in the core negative feedback loop, will also exist.In this case the two terms describing alternate paths through the positive feedback loop and negative loop (see Section 4.4) will again remain after shared components are factored out.As discussed above, in the cases where the magnitude of the path through the negative feedback loop is greater than the path through the positive feedback loop, the relationships between the response of u1 to its own inhibition and the feedbacks it supports described here will also be inverted.We will now consider the how inhibiting the external component u1 affects a core component u2.We will consider the effect of whether u1 feeds directly into the positive or negative feedback loop of the core RD network, and on whether u2 itself lies in the core positive or negative feedback loop.As for minimal systems (Section 4.2.2),u2 will decrease on the inhibition of u1 if the sign of the path from u1 to u2 is positive, unless the submatrix determinant det [(̃E)] contains the positive feedback loop of the core RD network.Likewise, u2 will increase if the path from u1 to u2 is negative, unless the path excludes the core RD network positive feedback loop.
If u1 feeds into the core RD network through the positive feedback loop, det [(̃E)] will not contain the core positive feedback loop.The response will therefore depend on the sign of ( 5 ), such that if u1 itself supports a positive feedback loop, the weight of the path from u1 to u2 will be positive if u1 and u2 are in-phase, and negative if they are out-of-phase.Therefore, u2 will decrease if it is in-phase with u1 and increase if it is out of phase.The opposite will occur if u1 is in a negative feedback loop, with the level of u2 increasing if it is in-phase with u1 and decreasing if it is out-of-phase.
We next consider the case when u1 feeds into the core through the negative feedback loop.If u2 is also found in this loop, and the path from u1 does not pass through the core positive feedback, the sign of the path will be positive if the components are in-phase and negative if they are outof-phase, provided u1 supports a positive feedback loop (and therefore provides a net negative feedback to the core positive feedback loop).The reverse is true if u1 supports a negative feedback loop (and thus a net positive feedback).As the core positive feedback components are found in det [(̃E)], if u2 is in-phase with u1 it will decrease on inhibition if u1 provides a net positive feedback, and out-of-phase components will increase, with the reverse true if u1 provides a net negative feedback.The same is also true when u1 feeds into the core negative feedback loop, but the path to u2 passes into or through the core positive feedback loop.When the core positive feedback components are found no longer in submatrix (Ẽ) and therefore the sign of det [(Ẽ)] is reversed, the relationship between the signs of the path from u1 to u2 and the phase relationship are now inverted by passing back into the core positive feedback loop, maintaining the above relationship.

Development • Supplementary information
Therefore, regardless of where the path from u1 feeds into the core RD network, if u1 provides a net negative feedback loop, u2 will increase if it is in-phase with u1, and decrease if it is outof-phase.The opposite will be true if u2 supports a net positive feedback.These constraints are summarized in Figure 4d.Again, there can be instances where two paths of opposing sign can be traced from u1 to u2 if u1 feeds into a component found in both core loops and u2 is also found in both loops.The response is again dependent on the sign of Equation (4.4), with the above responses inverted if the magnitude of the term relating to the path through the negative feedback loop is greater than those relating to the positive feedback loop, with all components receiving input from the two paths responding in a coupled manner.
4.3.3:Effect of adding external components to an RD system Finally, we considered the effect adding external components to existing constraints on the system.We asked whether the addition of an external component could change the behaviours within the minimal RD core.As shown above, a major feature determining the response of a component in a minimal RD network either to its own inhibition, or the inhibition of another component, is whether or not various submatrices contain the destabilizing positive feedback loop of the core RD network.In some instances, the addition of loops can change the stability of these submatrices, either making stable submatrices unstable by adding a positive feedback loop, or making unstable submatrices stable by the addition of further negative feedback.
However, it should be noted that as the various loops and paths that determine the behaviours within the core will still be present, only certain placements of feedback loops have such effects, and that they only affect certain components under particular parameterisations.There will still exists parameterisations where the responses described above (summarized in Figure 4b) are maintained.
The constraints applying to the addition of one external component will, therefore, apply to any number of external components, provided that the magnitudes of additional interactions are not sufficiently large to change the stability of any submatrices.This also opens up the possibility of considering potential interactions between external components u1 and u2 (for simplicity we will only consider the cases where the path between u1 and u2 passes through the core).The path will contain the path from u1 to a core component, plus an additional interaction to u2.Therefore, the response of u2 to the inhibition of u1 will be the same as that of the preceding core component in the path if the final interaction is positive, and the opposite if it is Development: doi:10.1242/dev.190553:Supplementary information Development • Supplementary information negative.The same is true of the response of external u2 to a core u1.Consequently, if u1 is involved in positive feedback (either as the core loop or as a source of additional positive feedback), then u2 will decrease if it is in-phase with u1, and increase if it is out-of-phase, while if u1 provides negative feedback the reverse is true.As has been discussed above, when two paths can be traced from u2 to u1, if the path through the negative feedback loop has a greater weighting these relationships will be inverted coupled to other such components.
Taken together, the effects of inhibiting u1 on itself, and on u2 primarily depend on whether u1 falls in the core positive feedback loop, or provides an additional net positive feedback, or whether it provides negative feedback.While for some parameterisations of some topologies, certain components will behave differently, these constraints on the behavior of minimal systems provide a framework to consider the behavior of RD networks in response to inhibition experiments.

Figure S1 :
Figure S1: Efficacy controls for inhibition of Wnt and BMP signalling In situ hybridisations on E13.5 palatal shelf explants for Wnt target gene Lef1 and BMP target gene Id1, cultured for 24 hrs in the small molecule inhibitor IWP-2 and dorsomorphin respectively.Contralateral shelves shown as vehicle controls.Well known domains of target gene expression in palatal region are lost upon inhibitor treatment; see loss of strong Lef1 expression in tooth upon IWP-2 treatment, and loss of strong Id1 staining at medial palate edge (see arrows).Darker regions in inhibitor treated specimens represent background staining in areas of thickened tissue.Anterior to the right, medial up.Scalebar = 200 µm.

Figure S2 :
Figure S2: Identifying rugae from Shh in situ hybridisation of whole mount explant cultures a) Shh in situ hybridisations on E13.5 palatal shelf explants cultured for 24 hrs.Red rectangle shows 150 µm band along which Shh staining intensity was measured.Anterior to right, medial up.Scalebar = 200 µm.b) Plot of mean Shh staining intensity for each AP position along the red rectangle in a), against the AP position relative to the posterior of the image.Blue dots and coarse dashed lines denote positions of ruga 8 and ruga 1, identified manually on specimen in a).Green and red dots denote half heights and one third heights respectively for each peak, as illustrated by fine dashed lines on anterior of ruga 2. Arrows denote the next position in the array with the same intensity as marked by the red dots, with red arrows identifying potentially fused rugae.c) Plot of mean Shh staining intensity as in b) with blue bars denoting the positions of the rugae, with anterior and posterior boundaries determined by the blue dashed lines.Black dashed lines denote position of rugae, taken as the manually identified positions for ruga 8 and ruga 1, and the ruga midpoints for the remaining rugae.

Figure S3 :
Figure S3: Effects of pathway inhibition on the pattern of the rugae a) Mean number of rugae as marked by clear peaks of Shh expression in E13.5 palatal shelves cultured for 24 hrs with the specified inhibitor (+) or vehicle (-).b-d) Plots of Shh stripe position b), width c) and stripe intensity d) for all ruga from explants for specified inhibitor treatment, plotted against values for the equivalent ruga in stage matched untreated (vehicle) control (see Methods for details of quantifications).Plots from paired untreated explants are also shown to show levels of variability in the absence of inhibitor.Solid grey lines denote equal values, dashed lines represent one median displacement around the line of equal values for untreated explants (see Methods for details).Magenta regions represent larger, and cyan regions smaller, values in inhibitor treated explants.For stripe position b), red points relate to rugae determined to have fused with another ruga in the inhibitor treated explants, with their position in plotted against the position of that ruga in the inhibitor treated explant.Yellow and green points relate to rugae in the inhibitor treated and untreated explants respectively, where the equivalent ruga could not be identified in the stage matched explant (respectively the untreated or inhibitor treated).For measures of stripe width c) and intensity d) points in red denote posterior rugae (see Methods).e) Percentage of rugae for each treatment lying within one median displacement from the line of equal values (grey), or above (magenta) and below (cyan), for each of the three measures in b-d).All plots depict quantifications from four pairs of explants.Drug treatments were carried out at the doses specified in the Methods, plus an additional dose, with points pertaining to lower doses shown as circles and higher dose as triangles in all plots.Doses used are: cyclopamine 20 µM and 80 µM, SU-5402 40 µM and 80 µM, IWP-2 10 µM and 50 µM, dorsomorphin 10 µM and 50 µM.

Figure S4 :
Figure S4: Titration of FGF inhibitor SU-5402 a) Shh expression in palatal shelves explanted at E13.5 and cultured with the FGF receptor inhibitor SU-5402 for 24 hrs at the indicated doses, alongside the contralateral shelves which were cultured in DMSO.Scalebar = 200 µm.b) Plot showing positive correlation between concentration of SU-5402 and Shh stripe width.c) Plot showing positive correlation between concentration of SU-5402 and Shh stripe greyscale intensity.Note that as more intense staining has a lower greyscale value, greyscale axis is plotted in reverse to show a decrease in staining intensity with increased concentration of inhibitor.Points in b) and c) depict ruga 4 to 2 from two embryos per concentration.Other rugae are not plotted as measurements were not calculated for ruga 1 (see methods), while more posterior rugae appear very variable even in control treatments (see Results).Pearson's product-moment correlation coefficient (r) is given for each plot in b) and c).

Figure S5 :
Figure S5: Quantification of striped rugal expression a) In situ hybridisation showing Shh expression in sagittal sections taken at 35 µm increments across a E13.5 palatal shelf, ordered from medial most (red outline) to lateral most (cyan outline).The positions of ruga 3 and ruga 8 are indicated.Dotted lines illustrate the extent of the palatal epithelium used for quantifications.Anterior to right.Scalebar = 200 µm.b) Intensity profiles for each section aligned by the position of ruga 3. Peaks in Shh intensity at ruga 3 and ruga 8 are indicated.Profiles coloured as in a).c) Mean intensity profile from the four sections in a) and b).Shaded area represents 1 +/-sd.

Figure S6 :
Figure S6: Additional pathway marker quantification In situ hybridisation of sagittal section through E13.5 palatal shelf for additional markers of pathway activity.Dotted lines illustrate the extent of the palatal epithelium and the underlying mesenchyme used for quantifications.Anterior to right.The intensity profile averaged across the palatal shelf shown for each specimen from which illustrated in situ is taken for gene of interest (coloured trace) and Shh (grey trace) for the epithelium and mesenchyme.Shaded areas represent 1 sd around gene of interest (for clarity of presentation, variation around Shh trace is not shown).For each marker, the number of specimens showing the observed pattern are: Ptc1 41, Erm 38, Spry2 36 and Axin2 43.Scalebars = 200 µm.

Figure S7 :
Figure S7: Adding a component as a mediator of an interactiona) Network topology comprising three components U, V and W.This topology can be interpreted as a two-component activator-inhibitor system between 'master morphogens' U and V, where the positive interaction from U to V is 'mediated' by a double inhibition (ie a positive interaction) through W (see green arrow in U-V AI system).Alternatively, this same three-component topology can be interpreted in a similar manner as a substrate -depletion system with 'master morphogens' U and W, where V 'mediates' the positive interaction from W to U (see red arrow in U-W SD system).b) All possible ways that an interaction in a two-component Wnt-Hh AI network (marked in green) can be 'mediated' by BMP in the adjacent three-component Wnt-BMP-Hh RD networks.Interaction in the three-component RD networks maintain the net sign of the marker interactions in the two-component networks; the phase of BMP relative to Wnt and Hh is not considered.c) Same as in b) for Wnt mediating interaction in a two-component BMP-Hh SD network (with interactions 'mediated' by Wnt in blue).

Figure S9 :
Figure S9: Inhibition of production produces the same outcomes as inhibition of response Violin plots showing percentage change in the mean level of components U and V in illustrated activator-inhibitor (AI) and substrate-depletion (SD) RD networks, on inhibition of the production of component U and V in RD simulations.

Figure S10 :
Figure S10: Three-component topologies and exclusion by successive constraints a) Three-component network topology showing the nine possible interactions captured in the reaction matrix.b) Enumeration of topologies recovered from parameter search under different constraints as detailed in text.

Figure S11 :
Figure S11: Responses of three-component networks to perturbation Full heat map showing the percentage of parameter sets where the level of each component increases in response to the inhibition of each component in the network in RD simulations.Topology names are as in figure 3. Hierarchical clustering shows that the behaviours of the topologies fall into the main categories for which all three components show similar behaviour.

Figure S12 :
Figure S12: Feedback loops in three-component networks a) Schematics showing the eight minimal 4-interaction RD networks recovered from the parameter search, with either the single positive feedback loop marked in magenta (left) or the single negative feedback loop marked in cyan (right).Topologies are group as according the the hierarchical clustering in figure S7. b) Schematics showing the ten minimal 5-interaction topologies networks recovered from the parameter search.On left, the net feedback between Wnt and BMP is marked, with magenta for a positive feedback, or cyan for a negative feedback.(It should be noted that for some topologies (such as iv, vii, iii and ix) the feedback is indirect with both the interaction from Wnt to BMP and BMP to Wnt pass through Hh.For example, in iv, while there is no single positive feedback loop, the net interaction between Wnt and BMP is a mutual inhibition and therefore a net positive feedback.)On the right, the presence of a negative feedback between Hh and the autoactivating component (Wnt or BMP) is marked in cyan.Topologies are group as according the the hierarchical clustering in figureS7.Wnt in blue, BMP in green and Hh in red.

Figure S13 :
Figure S13: Using topology to constrain the response to perturbation of minimal three-component networks a) Illustrative examples of minimal 4-interaction three-component RD networks showing the three different ways a single positive and negative feedback loop can be combined.The positive feedback loop can pass through one (topology 3) or two (topologies 1 and 2) components, and the negative feedback can pass through two (topology 1) or three (topologies 2 and 3) components.Positive feedback loops in magenta, negative feedback loops in cyan.b) Illustrative examples of 5-interaction three-component RD networks showing the five different ways that an external component can be wired into a two-component RD core, such that the new topology is not an elaboration of one of the minimal architectures in a).Core components and interactions in black, with external in grey.The nature of the net feedback between the external component (BMP in the illustrative examples) and the core positive feedback component (Wnt in the illustrative examples) is shown as a magenta loop for positive feedback and a cyan loop for negative feedback.c) Examples of three-component RD networks showing the six different orientations that a set of 3-component RD loops can take, using the 4-interaction minimal architecture in a) 1 as an illustrative example.The feedback loops can be rotated into three different orientations (1, 3 and 5), as shown by the positions of the two feedback loops relative to the dashed black line.In each orientation, the architectures can be reflected around the dashed black line to give a further three orientations (2,4 and 6).d) Response of Hh to inhibition of Wnt, BMP and Hh for all 45 possible minimal three-component Wnt-BMP-Hh topologies based on reaction term analysis.Topologies are specified by the combination of a topology (1 to 8 from a) and b)) and an orientation (1 to 6 from c)).Topology architecture 8 is symmetrical upon reflection and therefore two the topology under two possible orientations is the same.The direction of response of Hh shown by arrowhead, with up arrow indicating increase and down arrow decrease.Where the response is unconstrained a white bar is shown.Coupling between components (see supplementary note 4 and figure4) is not shown.The ten topologies for which the predicted responses are consistent with experimental observations are marked in red.These are the same topologies as recovered by the RD analysis (see figure3).

Figure S14 :
Figure S14: Calibration curves for converting palatal length into time a) Plot of time litters of mice were collected vs. the mean embryo mass for that time (n = 3 litters per time point).b) Plot of embryonic mass (taken as the midpoint of the bin) against average ruga 3 to ruga 8 length the the weight bin.For a) and b), equations of best fit quadratics are given which allowing the conversion of palatal lengths to time in embryonic days.

Figure S15 :
Figure S15: Illustrative example of how the appearance of periodic target gene expression relative to the onset of rugal Shh expression is identified a) Kymographs of epithelial Shh (green) and epithelial Id1 (magenta) as in figure 6a.b) Plots of normalised staining intensity for Shh (green) and Id1 (magenta) against time in embryonic days at four different distances from ruga 8 (i-iv), as marked by arrows and white dashed lines in a).c) Plots of normalised staining intensity for Shh against Id1 for 600 equal time steps from 12.5 to 14.0 embryonic days at four different distances from ruga 8 (i-iv), as in b).The Spearman rank correlation coefficient (r) is given for each distance.d) Plot of Spearman rank correlation coefficient (r) against distance from ruga 8 for Id1.Shaded area represents 95 % confidence interval from bootstrapping (see Methods).Solid line represents mean correlation coefficient calculate over anterior third of palatal tissue.Dashed lines demonstrate how the position relative to ruga 8 where half this value is reached is used to determine the onset of out-of-phase Id1 expression, marked by orange dashed lines in a).

Figure S16 :
Figure S16: Kymographs showing patterns of mesenchymal expression through time for Gli1 and Id1 Kymographs showing the pattern of expression of mesenchymal Gli1 and Id1 through time (magenta) and their expression relative to Shh (green) for rugae 3, 4 and 5. White arrowheads indicate the approximate onset of Shh expression for each ruga, and orange arrowheads indicate approximate positions of the change in the expression pattern of each target gene associated with each ruga.

.
Development: doi:10.1242/dev.190553:Supplementary information Development • Supplementary information Solving the system det [ − ] = 0 gives the dispersion relation for the eigenvalues  of S: Figure SN1 a) Graphs showing the two possible responses of u1 to the inhibition of u1 by strength a.The vertical asymptote (dashed line) can lie at a positive (i) or negative (ii) value of a.The horizontal asymptote lies along the a-axis.b) Graphs showing the six possible responses of u2 to the inhibition of u1 by strength a.The vertical asymptote (vertical dashed line) can lie at a positive (i, ii and vi) or negative (iii, iv and v) value of a.The horizontal asymptote (horizontal dashed line) can lie at a positive (ii, iii, v and vi) or negative (i and iv) value of u2.The intercept of the a-axis can lie at a positive (ii, iv and vi) or negative (i, iii and v) value of a. c) Graphs showing the three possible responses of u1to the inhibition of u1 by strength a.The vertical asymptote (vertical dashed line) can lie at a positive (ii and iii) or negative (i) value of a.The horizontal asymptote (horizontal dashed line) can lie at a positive (ii and iii) or negative (i) value of u1.The intercept of the a-axis can lies at a=1. α :10.1242/dev.190553:Supplementary information on the relative values of the asymptotes and intercepts, this function can take six forms (Figure SN1 b), for three of which (Figure SN1 b i-iii), u2 increases with inhibition of u1.
Figure SN2Full maps showing the percentage of parameter sets where the level of each component increases in response to the inhibition of each component in the network in reaction term analysis and in RD simulations.The left hand heat map shows the response for reaction term analysis for the parameter sets used in three-component RD analysis, evaluated using criteria from Section 2. The second heat map, showing the Reaction-Diffusion system with inhibition strength amax, is a reproduction of the heat map from Figure3, and the heat maps with inhibition strength 0.25.amax and 0.25 2 .amaxare from reruns of the three-component RD simulations, replacing the inhibition parameter value a with values of 0.25.a and 0.25 2 .arespectively.Pathway topology names are as in Figure3.
:10.1242/dev.190553:Supplementary information    Development • Supplementary informationwhere s(l) is the number of cycles in the L-spanning subgraph l.Equation (3.1) is known as the Coates formula; for more details on the derivation, seeMarcon et al. (2016) and references therein.

Table S2 : Behaviour of simulations for three-component RD topologies.
Table showing for each topology (as defined in figure3b) the behavior of RD simulation of the three components Wnt (W), BMP (B) and Hh (H), as detailed in supplementary table S1.As stated in Methods, for Wnt inhibition is at the level of production, while for BMP and Hh inhibition is at the level of response.

information Appendix S1. Supplementary Notes for Economou et al. " Perturbation analysis of a multi-morphogen Turing Reaction-Diffusion stripe patterning system reveals key regulatory interactions" 1: Review of Conditions for DDI in 2-3-and N-component systems
To satisfy the condition (2.10), det [] < 0, and as a consequence, det [] > det[].Otherwise, det [] > 0, and in this case, to also satisfy condition (2.11), it follows that again det [] > det[].However, when condition (2.12) holds, it follows that det [] < det[].Therefore, when det [] < det [], u1 increases when its production is inhibited, while when det[] > det [], u1 decreases when its production is inhibited.
or through the positive feedback loop, in which case the relationship between the nodes is reversed.As detailed in Equation (3.1), the relationship between the sign of det [(Ẽ)] and the feedback loops contained depends on the number of nodes in the system described by submatrix (Ẽ) (see Section 1).Again, considering all possible combinations of odd and even numbers of components in submatrices (Ẽ) and path ( 5 ), and given that ∑ C[] %"  % of u1 to its own inhibition if it falls in both feedback loops, and shared terms are factored out of Equation (4.2).In both instances, if the weight through the positive feedback loop is greater than that through the negative feedback loop, that is (i.e. the weight the path from u1 to a core positive feedback component multiplied by the weight of the return path).If u1 supports a net positive feedback, it will decrease on its own inhibition,