High-throughput mathematical analysis identifies Turing networks for patterning with equally diffusing signals

1How cells self-organize to form ordered structures is a central question in developmental biology, and identifying self-organizing mechanisms promises to provide new tools for synthetic biology and regenerative medicine. More than six decades ago, Alan Turing proposed a theoretical model in which interactions between diffusible substances can break the initial symmetry of cell fields to form periodic patterns. Subsequent work from Gierer and Meinhardt postulated that such self-organizing processes require differential diffusivity between a short-range self-enhancing activator and a feedback-induced long-range inhibitor. Numerous studies have proposed models based on these concepts to explain pattern formation during development, including skin appendage specification, lung branching, tooth development, rugae formation, and digit patterning. However, the evidence in support of specific activator-inhibitor pairs has been limited, and few studies have provided experimental support for the differential diffusivity of activators and inhibitors.

2Pattern formation processes are regulated by the interactions between secreted signaling molecules and their receptors that activate complex cell-autonomous signaling events. However, since most reaction-diffusion models have been reduced to abstract networks of two diffusible reactants, the influence of immobile cell-autonomous factors on reaction-diffusion patterning is largely unknown. Previous theoretical studies on selected network topologies have challenged the differential diffusivity requirements and indicated that in the presence of an immobile substance, patterns can form for a wider range of reaction and diffusion parameters. These and other studies suggest that extending models beyond abstract two-node systems can reveal different pattern formation requirements and may uncover new biologically relevant network designs. However, due to the complex mathematical analysis required to identify and understand such systems, extending reaction-diffusion models to more realistic signaling networks has been challenging, and the main assumption in the field has remained that complex models should reduce to simple systems that require an effective differential diffusivity.

3Here, we developed the freely available and user-friendly software RDNets (available at http://www.RDNets.com) to perform a high-throughput mathematical analysis of complex reaction-diffusion networks with non-diffusible components. In comparison to previous numerical studies, this method guarantees completeness, reproducibility, and detailed mechanistic insights into the principles underlying pattern formation. We used RDNets to build a comprehensive catalog of minimal three-node and four-node reaction-diffusion networks that include interactions between diffusible signals and cell-autonomous factors. Our results show that reaction-diffusion systems have three types of requirements for the diffusible signals depending on the network topology: Type I networks require differential diffusivity, Type II networks allow equal diffusivities, and Type III networks allow for unconstrained diffusivity. Overall, 70% of the networks identified by our analysis are of Type II and Type III and thus do not require differential diffusivity to form a spatial pattern. This reveals that realistic reaction-diffusion systems are based on mechanisms that are fundamentally different from the concepts of short-range activation and long-range inhibition based on differential diffusivity that have been predominant in previous models of pattern formation. Our software can be used to explore these new networks and is a unique tool to understand in vivo reaction-diffusion systems and to engineer synthetic circuits with spatial patterning capabilities.

4Understanding how complex gene regulatory networks control cellular behavior is a challenging problem in biology; even small networks can contain regulatory feedbacks that make systems behaviors difficult to predict. Mathematical biology has helped to identify network motifs that underlie basic behaviors such as oscillations, bi-stability or noise reduction, but this approach has been difficult to scale up to more complex networks and behaviors. Previous studies have overcome this obstacle by using numerical simulations to screen for topologies that implement a certain behavior. However, such simulations demand large computational power, their coverage is incomplete, and they do not have the explanatory power of analytical approaches. The ideal tool to analyze the behavior of gene networks should retain the explanatory power of mathematical approaches and yet be able to comprehensively screen for network topologies and the underlying mechanistic principles.

5We have developed the web-based software RDNets (http://www.RDNets.com) to derive a comprehensive catalog of minimal three-node and four-node reaction-diffusion networks and their pattern-forming conditions. Our analysis reveals that networks have different diffusivity requirements depending on the topology. RDNets can constrain candidate topologies with qualitative and quantitative experimental data, making it a convenient tool for users that aim to study developmental patterning networks or to design synthetic reaction-diffusion circuits.

6We developed an automated linear stability analysis to derive the pattern forming conditions of networks with N nodes (Figure 1a, Materials and methods). Linear stability analysis determines whether a system can form a pattern by testing i) if the concentrations of the reactants are stable at steady state, and ii) if diffusion-driven instabilities arise with small perturbations. Because of its mathematical complexity, this type of analysis has been the exclusive domain of mathematicians and systems biologists, and its application beyond two-reactant models has required dedicated theoretical studies for selected networks. To generalize the analysis to networks with more than two nodes, we utilized a modern computer algebra system and developed the software pipeline RDNets that automates the algebraic calculations. Within this framework, secreted molecules like ligands and extracellular inhibitors are represented by diffusible nodes, and cell-autonomous components such as receptors and kinases are represented by non-diffusible nodes. Our software analyzes networks with k interactions between the nodes; these interactions are represented by first order kinetics rates, where a positive rate corresponds to an activation and a negative rate to an inhibition.

7The software pipeline comprises six steps to identify patterning networks:

8Construction of a list of possible networks of size k.

9Selection of strongly connected networks without isolated nodes or nodes that solely act as read-outs.

10Deletion of symmetric networks, such that isomorphic networks are considered only once.

11Selection of networks that are stable in the absence of diffusion (i.e. homogeneous steady state stability).

12Selection of networks that are unstable in the presence of diffusion (i.e. instability to spatial perturbations).

13Analysis of the possible reaction-diffusion topologies associated with the networks and derivation of the resulting in-phase and out-of-phase patterns.

14Steps 4 and 5 represent the core part of the automated linear stability analysis and involve the majority of analytical computations. In Step 6, our software screens the possible reaction-diffusion topologies associated with a network. A reaction-diffusion network of size k defines only a set of k regulatory links between nodes but does not make any assumption on whether these are activating or inhibiting interactions. In the following, we refer to the possible combination of activating and inhibiting interactions as 'network topologies'.

15We used our software RDNets to systematically explore the effect of cell-autonomous factors in reaction-diffusion models for the generation of self-organizing patterns. We studied two types of networks: a) 3-node networks with two diffusible nodes and one non-diffusible node representing the interaction between two secreted molecules and one signaling pathway, and b) 4-node networks with two diffusible nodes and two non-diffusible nodes representing the interaction between multiple ligands and signaling pathways. Table 1 shows the number of networks identified at each step of our automated mathematical analysis (see Figure 1—figure supplements 1–4 for the complete catalog of the identified reaction-diffusion networks). Our analysis revealed that in the presence of cell-autonomous factors there are three types of networks with different constraints on the diffusible signals:

16where D is the list of diffusion coefficients that are non-zero.

17We found that 70% of the identified networks with non-diffusible nodes are of Type II and Type III (Figure 1b), showing that in the presence of cell-autonomous factors the differential diffusivity requirement is unexpectedly rare. Type III networks have never been characterized before and surprisingly have patterning conditions that are independent of specific diffusion rates. We found that Type III networks are not only numerous but also extremely robust to changes in parameter values compared to Type I and Type II networks (Figure 1b, Materials and methods). Using numerical simulations, we systematically confirmed our mathematical analysis and determined that a network can form all possible combinations of in-phase or out-of-phase periodic patterns depending on the network topology (Figure 1c, Appendix 1). Together, our results show that realistic reaction-diffusion networks are intrinsically robust, do not require differential diffusivity, and have patterning capabilities identical to classical two-node reaction-diffusion models. Importantly, the novel class of Type III networks that we discovered suggests a new mechanism of pattern formation that is independent of short-range activation and long-range inhibition based on differential diffusivity.

18To obtain insight into the organizing principles underlying the three types of networks identified by our high-throughput analysis, we developed a novel graph-theoretical formalism to express the pattern forming conditions in terms of network feedbacks rather than reaction parameters (see Materials and methods and Appendix 2). This analysis determines which feedback cycles contribute to the stability and the instability conditions (Figure 2a,b) and defines the topological features that underlie Type I, Type II, and Type III networks. In agreement with previous studies, our analysis confirmed that two-node networks can only simultaneously satisfy the stability and instability conditions when the diffusion ratio d between the inhibitor and the activator is greater than one (Figure 2, left column). This observation has been linked with the widespread belief that reaction-diffusion systems require differential diffusivity to implement short-range auto-activation and long-range inhibition. Our analysis instead suggests that the differential diffusivity requirement arises from the opposite nature of the stability and instability conditions, which require that the destabilizing feedback must be both higher and lower than the stabilizing feedback. Since the diffusion term only appears in the destabilizing condition, it assumes the role of a unique pivot that can satisfy both conditions simultaneously when d > 1. Our results indicate that the presence of non-diffusible nodes allows feedbacks that do not appear in the instability conditions to act as an additional pivot to satisfy both conditions simultaneously by increasing stability. This is the case for most Type II networks (Figure 2, middle column) that contain additional negative feedbacks that allow for equal diffusivities. Importantly, our analysis also reveals that non-diffusible nodes can implement positive feedbacks that can drive the network unstable independently of stabilizing feedbacks and for any diffusion ratio d. This is the case for Type III networks (Figure 2, right column), where the stability and instability conditions are uncoupled and can be simultaneously satisfied for large parameter sets. This is possible because immobile factors can act as 'capacitors' that retain and amplify perturbations independently of the reactants’ diffusion coefficients (see Appendix 3 for details). Such systems represent a fundamentally new pattern formation mechanism that has not been described previously.

19Together, our results show that models based on 'short-range auto-activation and long-range inhibition' implemented by differential diffusivity are only a special case of a general trade-off between stabilizing and destabilizing feedbacks required for pattern formation. The virtually indistinguishable simulations of Type I networks with differential diffusivity and Type II networks with equal diffusivities reveal that the final aspect of the periodic patterns does not reflect a difference in the range of activators and inhibitors but only a difference in their amplitude (Figure 2c, see Appendix 3 for details). Indeed, in other Type II and Type III networks the relationship between the amplitude of activators and inhibitors can even be inverted, such that the perceived range of the activator appears larger than the perceived range of the inhibitor. Therefore, in contrast to previous studies, we propose that long-range lateral inhibition is not required to limit the expansion of the activator (Appendix 3).

20To demonstrate the functionality and applicability of RDNets, we analyzed two known self-organizing developmental patterning networks, the Nodal/Lefty reaction-diffusion system and the BMP/Sox9/Wnt network. In the following, we show how quantitative and qualitative experimental data from these developmental systems can be used to constrain the high-throughput analysis and to characterize the possible underlying patterning topologies.

21It has been proposed that Nodal and Lefty implement an activator-inhibitor system that patterns the germ layers and the left-right axis in vertebrates (Figure 3a). In agreement with this hypothesis, the self-enhancing activator Nodal has been shown to diffuse 7.5 times slower than the feedback-induced inhibitor Lefty in living zebrafish embryos. The Nodal/Lefty system has been modeled as a two-component activator-inhibitor system, but the influence of cell-autonomous factors including receptors and the well-characterized intracellular signal transduction cascade via phosphorylated Smad2/3 has not been studied. We used our software to screen for networks that extend the two-node Nodal/Lefty system with a non-diffusible node corresponding to active Nodal signaling (Figure 3b). The screen was constrained with known qualitative regulatory interactions: a positive feedback loop between Nodal and its signaling, and a promotion of Lefty by Nodal signaling (Figure 3b). Moreover, we constrained the two negative self-regulations on Nodal and Lefty, which represent their clearance from the diffusible pool, with the previously measured clearance rate constants. Finally, we selected only reaction-diffusion networks that produced in-phase patterns of Nodal and Lefty, which recapitulate their overlapping expression domains. With these constraints, our mathematical analysis identified just two possible minimal networks: In one network Lefty inhibits Nodal signaling indirectly at the receptor level, and in the other network Lefty inhibits Nodal directly (Figure 3b). These predictions are in agreement with the two possible mechanisms by which Lefty has been proposed to inhibit Nodal activity: by binding to the Nodal receptor or by directly sequestering Nodal. However, the role and significance of these two alternative mechanisms for Nodal/Lefty-mediated patterning has remained unclear. Our mathematical analysis predicts that the first mechanism (Lefty blocks the receptor complex) determines a Type II network, whereas the second mechanism (Lefty blocks Nodal directly) determines a Type III network. Importantly, both models suggest that the Nodal/Lefty system may form patterns without differential diffusivity of activator and inhibitor. Using the clearance rate constants of Nodal and Lefty as quantitative constraints, our mathematical analysis predicts a possible minimum diffusion ratio d = 0.55 for the Type II network, whereas the Type III network allows for any combination of diffusion coefficients (Figure 3d). The robustness analysis of the networks shows that for unconstrained valued of d, the Type III network is more robust to parameter changes (Figure 3c). However, when we fix the diffusion ratio to the experimentally quantified value (d = 7.5), the Type II network becomes more robust than the Type III network (Figure 3d). This shows that, while Nodal and Lefty do not necessarily need to have different diffusivities to form a pattern, the combination of differential diffusivity and clearance rate constants increases the robustness of the Type II system.

22As a second example, we used RDNets to analyze the BMP/Sox9/Wnt (BSW) self-organizing network that underlies digit patterning. The expression patterns and the signaling activity of the network components have been well-characterized showing that Sox9 forms periodic expression peaks that are out-of-phase of BMP expression and Wnt activity (Figure 4a). A three-node reaction-diffusion network with two diffusible nodes for the secreted signals BMP and Wnt and a non-diffusible node for the transcription factor Sox9 has previously been derived based on the known regulatory interactions (Figure 4b). It was shown that this network recapitulates the out-of-phase pattern between BMP/Wnt and Sox9 and forms a pattern with extremely low differential diffusivity requirements (d = 1.25). Our comprehensive mathematical analysis reveals that this three-node system is just another topology of the reaction-diffusion network that we analyzed for the extended Nodal/Lefty system; it is therefore a Type II network that can potentially form a pattern even when BMP and Wnt have equal diffusion coefficients. In previous studies, this observation was missed because the clearance rates of BMP and Wnt had been assumed to be identical. However, as we showed in the previous example for Nodal and Lefty, if BMP is cleared faster than Wnt, the diffusion ratio can be equal to or lower than one, d ≤ 1. The three-node BSW model recapitulates the out-of-phase pattern between BMP/Wnt and Sox9, but due to its high abstraction level it does not explain the opposite BMP expression and BMP activity patterns observed in the experimental data (Figure 4a). We therefore used RDNets to screen for more complex models with five nodes that represent all components of the network: two diffusible nodes for BMP (B) and Wnt (W) and three non-diffusible nodes, one for the canonical BMP pathway through pSmad1/5/8 (Sm), one for the intracellular Wnt signaling cascade (β-catenin, β), and one for Sox9 (S). We selected only networks that formed in-phase and out-of-phase patterns reflecting the experimental data (Figure 4a). Previous studies showed that Sox9 is promoted by BMP signaling through pSmad1/5/8 and is inhibited by Wnt through β-catenin. Similar to the Nodal/Lefty example, we constrained the mathematical screen by incorporating these known regulatory interactions. Unexpectedly, the screen revealed that if β-catenin would directly inhibit Sox9, no network could recapitulate the out-of-phase patterns between BMP expression and BMP signaling. By performing a more general screen that left this interaction unconstrained, we found that the opposite BMP expression and signaling patterns can be obtained when β-catenin indirectly inhibits Sox9 through pSmad/1/5/8. RDNets also predicts that the most robust networks include the following additional interactions: i) a negative feedback from Sox9 to Wnt, ii) a negative feedback from pSmad1/5/8 to BMP, and iii) either a positive feedback from β-catenin to Sox9 or a negative feedback from β-catenin on Wnt (Figure 4b, gray arrows). Interestingly, the majority of networks identified by our screen was of Type III, suggesting that the proportion of Type III networks increases when more non-diffusible nodes are added.

23Although reaction-diffusion mechanisms have a simple network design, they exhibit unique self-organizing capabilities making them appealing for synthetic engineering. So far, the synthetic implementation of reaction-diffusion systems has been impeded by the small pattern-forming parameter space of simple two-node models, their requirement for differential diffusivity, and a general gap between abstract models and real sender-receiver reaction-diffusion circuits.

24RDNets provides a comprehensive catalog of reaction-diffusion networks that do not require differential diffusivity of the signaling molecules, which enables bioengineers to explore new mechanisms to form periodic spatial patterns in a robust manner. We demonstrate the utility of RDNets by proposing an extension to an existing synthetic circuit for cell-cell communication in yeast. The original synthetic circuit introduced a diffusible plant hormone, cytokinin isopentenyladenine (IP), and its receptor AtCRE1 into yeast (Figure 5a). This circuit was used to implement a sender-receiver and a quorum sensing mechanism based on a positive feedback loop between IP-signaling and IP (Figure 5a). We used RDNets to identify possible signaling networks that can extend this positive feedback with additional interactions to form a reaction-diffusion pattern. Since at least two diffusible nodes are required to form a pattern, we screened minimal 4-node networks that include the engineered positive feedback and candidate interactions with another diffusible node. In order to look for realistic and easily implementable signaling circuits, we explored only networks with interactions between diffusible nodes through non-diffusible factors representing intracellular signaling cascades. We also imposed self-regulations on diffusible nodes to be exclusively inhibitory, representing decay. With these constraints, our high-throughput analysis identified 16 minimal reaction-diffusion networks (5 Type I, 3 Type II, 8 Type III), of which the Type II and Type III networks were most robust to parameter changes (Figure 5—figure supplement 1). In the following, we demonstrate how the conditions derived by RDNets can be used to engineer the most simple and robust Type II network (Figure 5a - right). In addition to the positive feedback loop, this network contains three additional negative feedbacks: two are self-regulations that correspond to decay, and one is a negative feedback between the newly introduced diffusible node and the non-diffusible node representing the receptor. This network suggests that a simple extension to the circuit developed in could be obtained by a) destabilizing the signaling hormone and the receptor to increase their turn-over (c1, c2), and b) introducing another hormone that signals through the same receptor and implements a negative feedback loop to its own expression or activity (c3, Figure 5a - right).

25In addition to revealing possible topologies, our automated analysis provides the mathematical formulae of pattern forming conditions. This feature together with the specification of quantitative constraints can be used to calculate pattern forming parameter ranges. To determine the strength of the newly introduced negative feedbacks required for pattern formation, we constrained the positive cycle strength with the first-order kinetic rate quantified in by fitting measurements of the signaling activity of IP (Figure 5c). Moreover, we assumed that the diffusion and decay rates are similar. With these constraints, RDNets determined that the newly introduced negative feedback has to be stronger than the positive feedback and decay rates (Figure 5c - right). This could be implemented using the more responsive IP-signaling promoter (TR-SSRE) developed in to drive the expression of the inhibitor. This specific synthetic network represents only one possibility. We find other Type III networks to be even more robust to parameter changes, but they appear to require the design of more complex synthetic circuits (Appendix 4). Once a synthetic network is designed, RDNets can also be used to automatically derive kinetic models that can simulate the reaction-diffusion network (Figure 5d,e, Appendix 5). Numerical simulations can be used to investigate the qualitative aspect of the pattern and its spatial periodicity. In the long term, all these features open new avenues for designing synthetic reaction-diffusion circuits that could be coupled with gene expression to enable complex applications, such as fabrication of spatially patterned three-dimensional biomaterials and tissue engineering in mammalian cells.

26We developed the web-based software RDNets, which exploits a modern computer algebra system to identify new reaction-diffusion networks that reflect realistic signaling systems with diffusible and cell-autonomous factors. Our approach is a new example of high-throughput mathematical analysis, which has several benefits over previous numerical approaches. First, RDNets can be run from most web browsers and does not demand large computational power. Second, our mathematical analysis yields closed-form solutions and is complete, in contrast to numerical simulations that can necessarily only sample from a small region of the entire parameter space. Third, RDNets derives the conditions for pattern formation and therefore provides mechanistic explanatory power to the users. In addition, it helps to identify reaction-diffusion topologies that are in agreement with qualitative and quantitative experimental constraints, which makes it an unprecedented tool for users that aim to study developmental patterning networks or to design reaction-diffusion synthetic circuits.

27Motivated by theoretical studies that showed that non-diffusible factors can influence pattern forming conditions, we used our software to systematically explore the effect of non-diffusible reactants in reaction-diffusion models. Our analytical approach is both comprehensive and informative and reveals that depending on the network topology, reaction-diffusion systems can belong to three classes: Type I systems that require differential diffusivity, Type II systems that can form patterns with equal diffusivity, and Type III systems that form patterns independent of specific diffusion rates. In particular, the novel class of Type III networks has not been described before and challenges models of short-range activation and long-range inhibition based on differential diffusivity that have dominated the field of developmental and theoretical biology for decades (see Appendix 3 for details).

28We used RDNets to obtain new mechanistic insight into two developmental patterning systems. By using quantitative data to constrain possible patterning networks, we found a Type II and a Type III topology that extend the Nodal/Lefty activator-inhibitor system with realistic cell-autonomous signaling. In such extended networks, Nodal and Lefty do not necessarily need to have different diffusivities to form a pattern. However, our results suggest that the differential diffusivity can contribute to a more robust patterning system in Type II networks with indirect Nodal signaling inhibition. We propose that the high general robustness of the Type III network might have played a role for the evolution of the Nodal/Lefty reaction-diffusion system in the first place, and that the indirect Nodal signaling inhibition of Type II networks might have been fine-tuned during evolution (Figure 3—figure supplement 1). Similarly, we extended the three-node BSW digit patterning model with additional previously characterized cell-autonomous factors and constrained a five-node model with qualitative data. Our analysis identified realistic network topologies that accurately reflect the previously puzzling opposite pattern of BMP ligands and BMP signaling and predicts novel interactions between network components.

29Finally, we used RDNets to design a novel synthetic patterning circuit based on a previously engineered positive feedback module. Identifying a comprehensive catalog of gene networks that can perform a certain behavior has been shown to be a successful strategy to uncover the design space of stripe-forming networks, which can be directly useful to synthetic biology. In particular, this approach permitted a whole family of network mechanisms to be synthetically constructed in bacteria – all capable of forming a gene expression stripe in a bacterial lawn. Similarly, our software provides a comprehensive catalog of reaction-diffusion networks and enables bioengineers to explore new mechanisms to form periodic spatial patterns in a robust manner. These networks explicitly include non-diffusible factors that mediate signaling and are easy to relate with sender-receiver synthetic toolkits. In addition, we found that the majority of realistic reaction-diffusion networks eliminate the differential diffusivity requirement that is difficult to implement synthetically. The possibility to use qualitative and quantitative constraints to screen for reaction-diffusion networks makes RDNets a unique tool to customize patterning systems from initial starting networks. Moreover, the pattern-forming conditions derived by the software can be used to estimate parameter ranges and network robustness. Particularly promising is our finding that each network is associated with a set of topologies that exhaustively determine all the in-phase and out-of-phase relations between periodic patterns (Figure 5d). It is therefore possible to design network topologies that promote the co-localized expression of any desired combination of factors. This will enable novel applications in tissue engineering, where the co-localized expression of differentiating factors can be used to induce specific tissues. Coupled with the three-dimensional pattern-forming capabilities of reaction-diffusion mechanisms (Figure 5e), this could be used to devise new strategies for engineering scaffolds or tissues with complex architecture.

30In summary, our analysis defines new concepts of reaction-diffusion-mediated patterning that are directly relevant for developmental and synthetic biology. We demonstrate three applications of our software RDNets to understand developmental mechanisms and to design synthetic patterning systems, but this approach can be extended to numerous other patterning processes. We therefore expect that RDNets will contribute to the wide-spread use of mathematical biology and that a similar approach could be applied to other dynamical processes such as oscillations and traveling waves.

31We analyzed reaction-diffusion networks represented by a reaction matrix J and a diffusion matrix D of size NxN, where N is equal to the number of nodes. The matrix J corresponds to the Jacobian of the reaction-diffusion system and contains partial derivatives that describe the relative influence of one node on another. Elements of the reaction matrix represent the first order kinetics rates of the regulatory interactions in the network, where a positive rate corresponds to an activation and a negative rate to an inhibition. The matrix D contains the diffusion rates of the reactants along its principal diagonal and is zero otherwise.

32Our analysis aims to identify minimal reaction-diffusion networks, defined as the networks with the minimal number of edges k that can form a reaction-diffusion pattern. In the case of 2-node networks, it has been described that minimal reaction-diffusion networks must have 2x2=4 edges, and therefore only a completely connected network is allowed. This completely connected 2-node network allows for only two reaction-diffusion topologies: the 'activator-inhibitor system' that forms in-phase periodic patterns, and the 'substrate-depleted model' that forms out-of-phase periodic patterns. Our automated approach takes the following inputs through a graphical user interface: the number of network nodes N, constraints on J and D including reaction or diffusion rates set to zero, and the number of regulatory interactions k. This last parameter defines the number of edges that each network should have with an upper bound of NxN edges representing a completely connected network (Figure 1a). This parameter also defines the number of possible networks that are analyzed by the software, which is calculated according to

33This number represents the possible subsets of size k that can be taken from J and corresponds to the number of possible networks of size k.

34An important part of the automated high-throughput mathematical analysis is the derivation of the characteristic polynomial, a mathematical expression that determines the stability of the reaction-diffusion system, which is calculated from the determinant of a matrix that combines J and D, the 'wave number' q, and the eigenvalue λ. For 3-node networks, the characteristic polynomial has the form

35where λ is the eigenvalue associated with the reaction-diffusion system, and the coefficients and are polynomials formulated in terms of the elements of J, D and q (see Appendix 1). The eigenvalue λ determines the stability of the network: negative real solutions of λ represent a system that is stable around its steady state, while a positive real solution of λ represents an unstable system. The variable q that appears in the coefficients and is the wave number that is introduced by the linear stability analysis and is multiplied for D. For values q>0, this parameter defines the periodicity of the reaction-diffusion pattern. Step 4 of our pipeline entails finding the ranges of the reaction parameters in J and diffusion parameters in D for each network, for which the solutions λ are all negative when q=0. Similarly, Step 5 requires finding parameter intervals, for which at least one solution λ has a positive real part when q>0. For characteristic polynomials of degree higher than 2, this is usually done by using the Routh-Hurwitz stability criterion, a mathematical theorem that finds the necessary and sufficient condition for all negative roots in terms of the polynomial coefficients. However, as the number of network nodes N increases, finding these parameter intervals becomes challenging and tedious because the coefficients are also complex polynomials of high degree in q. We used a computer algebra system to automatically derive and analyze the Routh-Hurwitz criterion in terms of the coefficients. Finally, Step 6 requires to evaluate which of the possible topologies that exist for a given network are compatible with the pattern-forming conditions derived in Step 5 (see Appendix 1 for details).

36The complete analysis of minimal networks is limited by the existence of analytical solutions. According to the Abel-Ruffini theorem, there is no general algebraic solution for systems with more than four nodes. However, in practice many five-node networks can be solved if the constraints specified in the input of RDNets lead to a simplification of the coefficients of the characteristic polynomial, as is the case for the five-node digit patterning network (Figure 4). Analytical approaches become also challenging when further diffusible nodes are added and when minimal models are extended with additional interactions.

37We analyzed the robustness of the networks by calculating the volume of the parameter space that respects the pattern-forming condition in relation to the unit length multidimensional space of all the possible parameter values. This robustness measure corresponds to the probability of randomly picking pattern-forming parameters. The pattern-forming parameter volume is calculated with a multiple integral of the pattern-forming conditions over all the parameters of the reaction-diffusion networks, in the form

38where are the reaction parameters and the diffusion parameters, are the pattern-forming conditions of the networks, and and are the limits of reaction and diffusion variables that are set respectively to (-0.5, 0.5) and (0, 1) representing a multidimensional parameter space of unit side length.

39To investigate the topological basis of Type I, Type II, and Type III networks, we developed a new theoretical framework based on graph theory that can be used to rewrite the pattern-forming conditions in terms of network feedbacks rather than their reaction rates. Further details of this theory are provided in Appendix 2.

40Our web-based software RDNets was written in Mathematica (Wolfram Research Inc., Champaign, Illinois) and is available at http://www.RDNets.com. RDNets requires only the installation of the freely available Wolfram CDF player plugin (http://www.wolfram.com/cdf-player/). A simple graphical interface can be used to specify inputs and constraints and to run the high-throughput mathematical analysis. Constraints can be specified by clicking on the nodes or edges of the networks, or by providing specific values for the corresponding parameters (see User Guide available at http://www.RDNets.com). These constraints are automatically translated into mathematical formulae that are coupled with the symbolic linear stability analysis performed by the computer algebra system. The graphical user interface can be used to explore and simulate the list of reaction-diffusion topologies given as output of the linear stability analysis. Additional constraints can be progressively added to the analysis to refresh the output and to narrow down the list of candidate topologies (see User Guide available at http://www.RDNets.com).