We developed an automated linear stability analysis (Murray, 2003) 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 (Koch and Meinhardt, 1994; Satnoianu et al., 2000; Murray, 2003; Miura and Maini, 2004), and its application beyond two-reactant models has required dedicated theoretical studies for selected networks (Othmer and Scriven, 1971; White and Gilligan, 1998; Klika et al., 2012; Korvasova et al., 2015). 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.

Figure 1 with 4 supplements see all

Download asset Open asset

The software pipeline comprises six steps to identify patterning networks:

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

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

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

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

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

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

Steps 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'.

We 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:

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

We 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.

To 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 (Murray, 2003), 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 (Klika et al., 2012; Korvasova et al., 2015). 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.

Figure 2

Download asset Open asset

Together, 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 (Kondo and Miura, 2010), we propose that long-range lateral inhibition is not required to limit the expansion of the activator (Appendix 3).

To 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.

It has been proposed that Nodal and Lefty implement an activator-inhibitor system that patterns the germ layers and the left-right axis in vertebrates (Chen and Schier, 2001; Shiratori and Hamada, 2006; Shen, 2007; Meinhardt, 2009; Schier, 2009; Kondo and Miura, 2010; Rogers and Schier, 2011; Korvasova et al., 2015) (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 (Müller et al., 2012). The Nodal/Lefty system has been modeled as a two-component activator-inhibitor system (Nakamura et al., 2006; Müller et al., 2012), but the influence of cell-autonomous factors including receptors and the well-characterized intracellular signal transduction cascade via phosphorylated Smad2/3 (Schier, 2009) 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 (Müller et al., 2012). Finally, we selected only reaction-diffusion networks that produced in-phase patterns of Nodal and Lefty, which recapitulate their overlapping expression domains (Schier, 2009). 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 (Chen and Shen, 2004). However, the role and significance of these two alternative mechanisms for Nodal/Lefty-mediated patterning has remained unclear (Cheng et al., 2004; Middleton et al., 2013). 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 (Müller et al., 2012) (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.

Figure 3 with 1 supplement see all

Download asset Open asset

As a second example, we used RDNets to analyze the BMP/Sox9/Wnt (BSW) self-organizing network that underlies digit patterning (Sheth et al., 2012; Raspopovic et al., 2014). 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 (Raspopovic et al., 2014). 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 (Raspopovic et al., 2014) 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.

Figure 4

Download asset Open asset

Although reaction-diffusion mechanisms have a simple network design, they exhibit unique self-organizing capabilities making them appealing for synthetic engineering (Diambra et al., 2015). 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 (Carvalho et al., 2014), and a general gap between abstract models and real sender-receiver reaction-diffusion circuits (Marcon and Sharpe, 2012; Barcena Menendez et al., 2015).

RDNets 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 (Chen and Weiss, 2005). 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 (Murray, 2003), 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 Chen and Weiss (2005) could be obtained by a) destabilizing the signaling hormone and the receptor to increase their turn-over (c₁, c₂), and b) introducing another hormone that signals through the same receptor and implements a negative feedback loop to its own expression or activity (c₃, Figure 5a - right).

Figure 5 with 1 supplement see all

Download asset Open asset

In 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 Chen and Weiss (2005) 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 Chen and Weiss (2005) 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 (Chen and Weiss, 2005; Carvalho et al., 2014).

We 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.

Our 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 (Murray, 2003), 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

This 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.

An 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

where λ is the eigenvalue associated with the reaction-diffusion system, and the coefficients ${\displaystyle a_1,a_2}$ and $a_3$ 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 ${\displaystyle a_1,a_2}$ and $a_3$ 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 (Murray, 2003), a mathematical theorem that finds the necessary and sufficient condition for all negative roots in terms of the polynomial coefficients ${\displaystyle a_1,a_2...a_n}$. However, as the number of network nodes N increases, finding these parameter intervals becomes challenging and tedious because the coefficients ${\displaystyle a_1,a_2...a_n}$ 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 ${\displaystyle a_1,a_2...a_n}$. Finally, Step 6 requires to evaluate which of the $2^k$ possible topologies that exist for a given network are compatible with the pattern-forming conditions derived in Step 5 (see Appendix 1 for details).

The 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.

We 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

where $k_1\mathrm{...}{k}_{NxN}$ are the reaction parameters and $d_1\mathrm{...}{d}_n$ the diffusion parameters, ${\displaystyle f(k_1\dots k_{NxN},d_1\dots d_n)}$ are the pattern-forming conditions of the networks, and ${\displaystyle l_1,l_2}$ and ${\displaystyle d{l}_1,d{l}_2}$ 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.

To 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.

Our 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).

We first generate a list of possible networks with $k$ edges. This is done by systematically selecting all the possible combinations of $N-k$ elements of $\mathbf{𝐉}$ that can be set to zero. The number of combinations is calculated as the binomial coefficient

For each of these combinations we derive a matrix $C$ with $N-k$ elements set to 0 and the remaining elements set to 1. Appendix 1—figure 1 shows an example of a $C$ matrix and its corresponding network in the case of $N=3$ and $k=6$.

Appendix 1—figure 1

Download asset Open asset

Importantly, we define the minimal size $k$ of a reaction-diffusion network of $N$ nodes as the minimal number of edges $k$ for which at least one of the possible networks can satisfy the requirements for Turing instabilities. In the case of $N=2$, it is well known that all the elements of $\mathbf{𝐉}$, that is $2x2=4$, have to be different than zero, and therefore the minimal network size is $k=4$ (Murray, 2003). In other words, all the regulatory edges of 2-node networks have to be present to be able to form a pattern, and therefore there is just one possible network. By analyzing 2-node, 3-node, 4-node, and 5-node networks, we empirically identified the minimal network sizes shown in Appendix 1—Table 1.

From the matrices $C$, we construct a list of adjacency directed graphs that represent each network. These graphs can be constructed by deriving the adjacency matrices from $C$ as the transpose $C^T$. Finally, we use an in-built function of the computer algebra system Mathematica to select matrices $C$ that correspond to strongly connected graphs. This filter guarantees that networks with isolated nodes or nodes that solely act as read-outs are discarded (Appendix 1—figure 2).

Appendix 1—figure 2

Download asset Open asset

The third step of our automated analysis removes all symmetric networks. Symmetric networks are defined as networks whose graphs are isomorphic. To find isomorphic networks, we extended the in-built isomorphism test function of Mathematica to take into account all reaction or diffusion constraints that may introduce additional differences. An example of a symmetric network and their isomorphisms is given in Appendix 1—figure 3.

Appendix 1—figure 3

Download asset Open asset

This step is a central part of the linear stability analysis and required us to develop a program that derives the characteristic polynomial of a reaction-diffusion system in a symbolic manner. This was done by calculating the determinant of a matrix $A$ defined as

where $\lambda$ is the eigenvalue and $I$ the identity matrix. The coefficients of the characteristic polynomial $({\text{a}}_1,\mathrm{\cdots },{\text{a}}_N)$ are also polynomials in terms of $\mathbf{𝐉}$,$\mathbf{𝐃}$ and $q$. A diffusion-driven instability requires that the solutions of the characteristic polynomial are all negative when $q=0$ and that there is at least one real positive solution when ${\displaystyle q>0}$:

As the number of nodes $N$ increases, finding the analytical solution of Equation (5) becomes unfeasible. As an alternative, the Routh-Hurwitz stability criterion can be used to derive the necessary and sufficient conditions to satisfy Equation (6) by deriving Routh-Hurwitz terms that are obtained by combining the coefficients $({\text{a}}_1,\mathrm{\cdots },{\text{a}}_N)$. Our software RDNets automatically derives the Routh-Hurwitz terms ${\mathrm{\Delta }}_i$ for the general case of degree $N$. For example, in the case of $N=4$ the Routh-Hurwitz terms are

The Routh-Hurwitz stability criterion states that all eigenvalues have a negative real part if and only if ${\mathrm{\Delta }}_i>0$ for i=1 to $N$. As mentioned above, the coefficients ${\text{a}}_1$ to ${\text{a}}_4$ are polynomials in terms of $\mathbf{𝐉}$,$\mathbf{𝐃}$ and $q$, which quickly become complex as $N$ increases. With $N=4$, the general $\mathbf{𝐉}$,$\mathbf{𝐃}$ matrices, the characteristic polynomial and the corresponding coefficients are

It is evident that even in the case with $N=4$ the coefficients have complex formulae. Moreover, when substituted in Equation (8), they give rise to further complex Routh-Hurwitz terms of higher order in $q$. The stability conditions are derived when $q=0$, which simplifies the coefficients. However, even in this case the manual derivation of the conditions that make all Routh-Hurwitz terms positive is challenging. We therefore automated these tedious algebraic calculations using the computer algebra system in Mathematica.

In the previous section we showed that the conditions for the stability of a reaction network can be derived by imposing that all the Routh-Hurwitz determinants are positive at $q=0$. Conversely, the requirements for the existence of diffusion-driven instabilities are obtained by imposing that at least one Routh-Hurwitz determinant turns negative when diffusion is considered: ${\displaystyle \mathrm{\exists }i,|{\mathrm{\Delta }}_i(q)<0}$ for some ${\displaystyle q>0}$.

In addition, the Routh-Hurwitz theorem allows – by checking which of the determinants turn negative – to derive the conditions that lead to stationary Turing patterns or oscillatory patterns. This can be done by considering the necessary and sufficient conditions for the absence of diffusion-driven instabilities presented in previous studies (Kellog, 1972; Cross, 1978; Othmer, 1982). A violation of any of these conditions is necessary and sufficient for the existence of a single real eigenvalue crossing the right half plane and thus for the formation of a stable Turing pattern. This occurs if and only if ${\displaystyle {\mathrm{\Delta }}_N(q)<0}$ for ${\displaystyle q>0}$, which can be rewritten as the simpler equivalent condition ${\displaystyle a_N(q)<0}$ and ${\displaystyle a_k(q)>0}$ for ${\displaystyle q>0}$ and for all ${\displaystyle k}$ different than $N$, by taking into account the identities between the $a_k$ coefficient and $k_{th}$ elementary symmetric functions of the eigenvalues (Horn and Johnson, 1990). Since in this work we are interested in the analysis of stable Turing patterns, we can use these conditions to select for networks that produce a stable pattern and filter out those that produce oscillatory patterns.

As shown in the previous section, these conditions generally lead to complex formulae. Non-diffusible factors are associated to vanishing entries in the diagonal of $\mathbf{𝐃}$, which reduces the complexity of the equations. Nevertheless, the negativity conditions of Routh-Hurwitz determinants (see for example ${\mathrm{\Delta }}_3$ in Equation [8]) remain tedious to derive. However, the algebraic calculations are automated with the computer algebra system in RDNets.

The last step of the analysis consists of finding when the stability conditions and the instability conditions for the existence of Turing patterns are fulfilled simultaneously. Importantly, the computer algebra system allows to handle the combinatorial explosion of algebraic cases systematically in moderate computational time.

The networks derived in Step 5 represent connectivities between nodes (matrix $C$) associated with coefficients of the characteristic polynomial that can satisfy the pattern-forming conditions. These conditions are written in terms of $J$ and $D$ but do note make any explicit assumption on whether elements of the reaction matrix $J$ must have a positive (representing activation) or a negative (representing inhibition) value. We define a network topology as a set of signs associated with reaction rates of $J$ whose correspondent elements in $C$ are set to one. Given a network of size $k$, there are $2^k$ possible topologies that encode all the possible combinations of positive and negative signs for a given matrix $C$. Reaction-diffusion topologies are topologies that can satisfy the pattern forming conditions. Our high-throughput analysis of minimal 3-node and 4-node signaling networks reveals that every unconstrained network is associated with a set of topologies that exhaustively determine all possible in-phase and out-of-phase periodic patterns between network nodes (Appendix 1—figure 4).

Appendix 1—figure 4

Download asset Open asset

RDNets can analyze the list of networks topologies to select for a specific phase of the periodic patterns. This feature requires to perform both analytical and numerical computations that derive the relative sign of the eigenvectors associated with a network topology. First, a random set of parameters that satisfies the pattern forming conditions is obtained for each network topology using an in-built function of the computer algebra system. The parameters are substituted into the stability matrix (Equation [3]) to leave $q$ as the only free parameter, and the solutions of the characteristic polynomial $\lambda$(Equation [5]) are calculated. The solution ${\lambda }_r$ without a complex part is selected, and its dispersion relation is analyzed to identify $q_{max}$ for which the eigenvalue $\lambda$ is maximum using a gradient ascend numerical method.

Finally, $q_{max}$ is substituted into $q$ in Equation (3), and the eigenvectors of the matrix are calculated (a list of eigenvectors for each solution $\lambda$). Each component of the eigenvector associated with ${\lambda }_r$ represents a reactant, and the relative sign between them represents the phase of the patterns. Reactants with eigenvector components of the same sign will be in-phase, while reactants with eigenvector components of opposite sign will be out-of-phase.

As mentioned in Appendix 1, a diffusion-driven instability occurs when the characteristic polynomial in Equation (5) has all negative real solutions with $q=0$ and further has at least one positive real solution when ${\displaystyle q>0}$. Given the Routh-Hurwitz conditions, the necessary and sufficient conditions to respect these requirements can be formulated in terms of the coefficients $({\text{a}}_1,\mathrm{\cdots },{\text{a}}_N)$ of the characteristic polynomial. In the following, we derive a new graph-theoretical formalism to express the coefficients $a_k$ in terms of cycles of the graph associated with the reaction matrix $\mathbf{𝐉}$.

Let ${\gamma }_k=\{i_1,\mathrm{\dots },i_k\}$ be a sequence of $k$ distinct integers such as ${\displaystyle 1⩽i_1<i_2....<i_k⩽N}$, and let $S_k^N$ be the set of all the different ${\gamma }_k$ sequences of $k$ elements in $\{1,\mathrm{\dots },N\}$. ${\displaystyle {\mathbf{F}}^{RD}({\gamma }_k)}$ denotes the $k$-by-$k$ principal submatrix of ${\mathbf{𝐅}}^{RD}$ given by the coefficients with row and column indices equal to ${\gamma }_k=\{i_1,\mathrm{\dots },i_k\}$. There are $N!\mathrm{/}(N-k)!k!$ different $k$-by-$k$ principal submatrices ${\mathbf{𝐅}}^{RD}({\gamma }_k)$, which are in a one-to-one correspondence with all the different sequences ${\gamma }_k$ in $S_k^N$. The sum of the determinants of all the different $k$-by-$k$ principal submatrices is denoted by

The following identity, which can be verified using the Laplace expansion of the determinant (Horn and Johnson, 1990), expresses the coefficients of the characteristic polynomial in Equation (5) in terms of $E_k$

The new formulation of the coefficients $a_k$ of the characteristic polynomial can be expanded to partially uncouple the contribution of the diffusion and reaction terms:

Uncoupling of the diffusion and reaction contributions is achieved by expanding the minors of order $k$ in Equation (12) and grouping the resulting terms according to the number of entries from the diffusion matrix. In this way, each minor $det[{\mathbf{𝐅}}^{\mathbf{𝐑𝐃}}({\gamma }_k)]$ is expressed as a summatory of products of all possible minors $det[\mathbf{𝐉}({\gamma }_m)]$ of order $m⩽k$ given by the sequences $\gamma =\{i_1,\mathrm{\dots },i_m\}$ multiplied by the $k-m$ coefficients $(d_{j_1}\cdot d_{j_2}\mathrm{\dots }\cdot d_{j_{k-m}})$ given by the complementary set $\overline{\gamma }=\{j_1,\mathrm{\dots },j_{k-m}\}$ in $\mathbf{𝐃}$. The subsets ${\overline{\gamma }}_m$ and ${\gamma }_m$ are complementary sets in ${\gamma }_k$ in the sense that ${\gamma }_m\cap {\overline{\gamma }}_m=\mathrm{\varnothing }$ and ${\gamma }_m\cup {\overline{\gamma }}_m={\gamma }_k$. Then, after some tedious algebra, the following expression is reached:

The first and third summands in Equation (13) stem exclusively from reaction and diffusion terms. The second is formed by minors of the reaction matrix weighted by the complementary coefficients of the diffusion matrix. Next, we show how determinants can be formulated in terms of cycles of the reaction graph associated with a matrix. This allow us to reformulate Equation (13) in terms of feedbacks of the reaction-diffusion network.

The definition of the reaction and interaction graphs closely follows the definition of the Coates graph of a general square matrix (Brualdi and Cvetkovic, 2008). The reaction graph ${\text{G}}_R[\mathbf{𝐉}]$ is a labeled, weighted, directed graph associated to the linearization of the reaction-diffusion system. In a system with $N$ interacting species, ${\text{G}}_R[\mathbf{𝐉}]$ is a graph with $N$ nodes that has a directed edge from node $j$ to node $i$ if ${\mathbf{𝐉}}_{ij}\ne 0$. The weight assigned to the edge is the coefficient ${\mathbf{𝐉}}_{ij}$. Note that according to this definition, the entries ${\mathbf{𝐉}}_{ii}$ in the diagonal of the reaction matrix have associated an edge with $i$ as the initial and terminal node. These type of edges are called loops and account for decay terms and autocatalysis in the reaction. The interaction graph ${\displaystyle {\text{G}}_I[{\mathbf{F}}^{\text{RD}}]}$ is the equivalent of the reaction graph including the diffusion term in ${\mathbf{𝐅}}^{RD}$. If the diffusion matrix is diagonal, both graphs are topologically identical and the only difference lies in the weight of the loops. The following $4$x$4$ matrix $\mathbf{𝐀}$ is used as an example to illustrate these definitions:

According to the definitions given previously, $G_R[\mathbf{𝐀}]$ is the $4$-node graph shown in Appendix 2—figure 1.

Appendix 2—figure 1

Download asset Open asset

The definitions from graph theory introduced next will be necessary to develop the framework for the analysis of the stability of a reaction-diffusion system. The indegree and the outdegree of a node are the number of edges that have this node as the initial or terminal node, respectively. A loop, defined as an edge that originates and ends at the same node, contributes 1 to both the indegree and the outdegree of that node. As an example, in the graph of Appendix 2—figure 1, node $3$ has indegree equal to $2$, for it has an incoming edge from node 1 and the loop. The outdegree of this node is $3$, because there are two edges going to nodes $1$ and $4$ plus the loop contribution (Appendix 2—figure 2).

Appendix 2—figure 2

Download asset Open asset

A cycle of length $m$ is a subset of $m$ distinct nodes and $m$ distinct edges that join $i_k$ to $i_{k+1}$ for $k=1,\mathrm{\dots },m$ and an edge from $i_m$ back to $i_1$. By this definition, loops are also cycles of length one. The weight of a cycle $w(c)$ is the product of weights of the edges that form the cycle. Cycles are classified as positive or negative according to the sign of their weight. The graph $G_R[\mathbf{𝐀}]$ used as an example has, aside from the four loops associated to the diagonal terms in $\mathbf{𝐀}$, a negative cycle of length 2 and a negative cycle of length 3. $C_2$ is negative and its weight is ${\displaystyle w(C_2)=-e\cdot c}$, whereas $C_3$ has weight $w(C_3)=-e\cdot f\cdot d$ and is also a negative cycle (Appendix 2—figure 2).

A subgraph of the reaction graph is a directed graph formed by a subset of edges and whose set of nodes are a subset ${\gamma }_k=\left\{i_1\mathrm{\dots }{i}_k\right\}$ of those in $G_R$, with ${\gamma }_k\subset \left\{\mathrm{1\hspace{0.17em}...}N\right\}$. The induced subgraph of ${\gamma }_k$, referred as the I-subgraph $I_{{\gamma }_k}$, is the subgraph of ${\displaystyle G_R[\mathbf{A}]}$ formed by the subset of nodes ${\gamma }_k$ and all the edges that join nodes within this set. The induced subgraph $I_{{\gamma }_k}$ is identical to the graph ${\displaystyle {\text{G}}_{\text{R}}[\mathbf{A}({\gamma }_k)]}$ obtained by applying the definition of the reaction graph to the principal submatrix $\mathbf{𝐀}({\gamma }_k)$, so that all the definitions and properties of the reaction graph carry over its I-subgraphs. As an example, consider the $3$-by-$3$ principal submatrix matrix $\mathbf{𝐀}({\gamma }_3)$ induced by the sequence ${\gamma }_3=\{1,2,4\}$:

The I-subgraph associated with ${\gamma }_k$ (Appendix 2—figure 3) is obtained by applying the definition of the reaction graph to the matrix $\mathbf{𝐀}({\gamma }_3)$ or, equivalently, by erasing from the complete graph $G_R[\mathbf{𝐀}]$ the nodes that do not belong to ${\gamma }_3$ and the edges that do not start and finish in the nodes of ${\gamma }_3$. Note that the induced subgraphs of the interaction graph $G_I[{\mathbf{𝐅}}^{\mathbf{𝐑𝐃}}({\gamma }_k)]$ obtained by considering all possible sequences ${\gamma }_k$ are in a one-to-one correspondence with ${\mathbf{𝐅}}^{\text{RD}}({\gamma }_k)$, the $k$ x $k$ principal submatrices appearing in the expansion of the $k$-th coefficient of the characteristic polynomial in Equations (10) and (11). Likewise, all the terms ${\mathbf{𝐉}}^{\mathbf{𝐑}}({\gamma }_k)$ in Equation (13) of the coefficient $a_k$ correspond to one and only one I-subgraph of the reaction graph $G_I[{\mathbf{𝐉}}^{\mathbf{𝐑}}]$. Hence, the graph definitions introduced previously provide a method to associate a graph to each of the terms in the algebraic stability conditions.

Appendix 2—figure 3

Download asset Open asset

A spanning subgraph is a subgraph that includes all the nodes in $G_R$, but not necessarily all edges. A linear spanning subgraph $\mathrm{\ell }$, also referred to as an L-subgraph, is a spanning subgraph of $G_R$, in which each node has indegree 1 and outdegree 1. This definition implies that an L-subgraph is composed by a set of disjoint cycles and isolated loops, where the cycles are disjoint in the sense that each node belongs to one and only one cycle. The three different L-subgraphs contained $G_R[\mathbf{𝐀}]$ are depicted in Appendix 2—figure 4.

Appendix 2—figure 4

Download asset Open asset

The number of cycles in a L-subgraph is denoted by $s(\mathrm{\ell })$. The weight of a L-subgraph is simply defined as the product of weights of the cycles contained in it:

For the linear spanning subgraphs depicted in Appendix 2—figure 4, the number of cycles and weights are:

The notion of linear spanning subgraphs can be naturally extended to the I-subgraphs of $G_R[\mathbf{𝐀}]$. The L-subgraphs contained in $I_{{\gamma }_k}$ are all the different subgraphs of order $k$ formed by a set of disjoint cycles spanning the $k$ nodes of the induced subgraph.

The L-subgraphs contained in the reaction graph of a matrix are the factors that determine the stability of the system. It has been shown that the graph methodology provides a way to assign an I-subgraph to each of the terms appearing in the algebraic equations that determine the stability of the system. Particularly, these equations are expressed in terms of principal minors $det[{\mathbf{𝐉}}^{\mathbf{𝐑}}({\gamma }_k)]$ of the reaction matrix. Next, it will be shown that the value of each of these minors is determined only by the weights of the L-subgraphs contained in the associated I-subgraph $I_{{\gamma }_k}$. The intimate relationship between the dynamics of a reaction-diffusion system and the cyclical structure of the reaction network is explained by this fact.

The expression of the determinant of an $N$ x $N$ matrix $\mathbf{𝐀}$ as a linear combination of the weights of the L-subgraphs in $G_R[\mathbf{𝐀}]$ is known as the Coates formula:

where the sum goes through all the L- subgraphs in $G_R[\mathbf{𝐀}]$. A sketch of the proof of the Coates formula is given following the work of Chen, 1997 (p. 143), and a more formal proof can be found in Brualdi and Cvetkovic, 2008 (p. 94). The first part of the proof shows that there is a one-to-one correspondence between the non-vanishing terms in the determinant of a matrix and the linear spanning subgraphs in its associated graph. The second part of the proof shows that the sign of the contribution of the non-vanishing terms is also dictated by the structure the linear spanning subgraphs. The classical definition of the determinant of a $N$ x $N$ matrix is:

where the sum is over all the $N!$ permutations $p=\{1,\mathrm{\dots },N\}\to \{i_1,\mathrm{\dots },i_N\}$. The signature of the permutation is given by ${\epsilon }_{i_1,\mathrm{\dots },i_N}$ and is equal to ${\displaystyle +1}$ if $p$ is an even permutation and equal to $-1$ if it is odd. All non-vanishing terms in Equation (19) are a product $a_{1i_1}\cdot \mathrm{\dots }\cdot a_{{}_{N{i}_N}}$ of $N$ coefficients. Each index appears twice, one as a row index and one as a column index, so that each row and column contribute to the product with exactly one coefficient. Hence, the subgraph in $G_I[\mathbf{𝐀}]$ defined by the entries $a_{1i_1}\cdot \mathrm{\dots }\cdot a_{{}_{N{i}_N}}$ has $N$ edges with exactly one edge coming into every node and one edge coming out of every node. Therefore, the subgraph associated to every term in Equation (19) is by definition a linear spanning subgraph in $G_I[\mathbf{𝐀}]$.

Conversely, every linear spanning subgraph $\mathrm{\ell }$ in $G_I[\mathbf{𝐀}]$ has $N$ nodes with indegree and outdegree equal to one. The $N$ edges in $\mathrm{\ell }$ are associated to $N$ coefficients in $\mathbf{𝐀}$, the edge directed to node $j$ being the only one in the $j$-th row, and the edge coming out of node $j$ being the only one in the $j$-th column. Arranging the indexes by increasing row order, the weight of $\mathrm{\ell }$ becomes $w(\mathrm{\ell })=a_{1i_1},\mathrm{\dots },a_{{}_{N{i}_N}}$, showing the correspondence between each linear spanning subgraph in $G_I[\mathbf{𝐀}]$ with one and only one of the permutations $p$ in the definition of the determinant. In this way, a one-to-one correspondence has been established between the L-subgraphs in $G_I[\mathbf{𝐀}]$ and the non-vanishing permutations terms in the $det[\mathbf{𝐀}]$. Explicit calculation of the determinant of the example matrix illustrates the first part of the proof, as det$[\mathbf{𝐀}]$ is proven to be a linear combination of the weights of the L-subgraphs ${\mathrm{\ell }}_1$, ${\mathrm{\ell }}_2$, ${\mathrm{\ell }}_3$ represented in Appendix 2—figure 4:

The one-to-one correspondence provides a convenient way to label a particular L-subgraph by the associated permutation. The permutation $p=\{i_1,\mathrm{\dots },i_N\}$ defines univocally the L-subgraph $\mathrm{\ell }(p)$ as the subgraph of $G_I[\mathbf{𝐀}]$ obtained by selecting the edge from node $i_1$ to node 1, from node $i_2$ to node 2 and generally from node $i_k$ to node $k$ for $k=1,\mathrm{\dots },N$. In this way, the sign of the contribution of an L-subgraph to the determinant can be derived considering the signature of its associated permutation. The L-subgraphs ${\mathrm{\ell }}_1$ and ${\mathrm{\ell }}_3$ of the example correspond to the even permutations $p_1=\{1,2,3,4\}$ and $p_3=\{4,2,1,3\}$. Thus, the corresponding terms in the determinant must have positive sign, as it is confirmed examining the explicit expression in Equation (20). Conversely, ${\mathrm{\ell }}_2$ is associated to the odd permutation $p_2=\{3,2,1,4\}$, and consequently the sign of the corresponding term is negative. In the same way that permutations define univocally a L-subgraph, a cyclic permutation of $k$ integers defines univocally a cycle passing through $k$ nodes in $G_I[\mathbf{𝐀}]$.

The second part of the proof shows how the signature of a permutation ${\epsilon }_p$ is related to the structure of the associated L-subgraph; more precisely, the number of cycles $s(\mathrm{\ell })$ contained in it. The parity of a permutation is the number of transpositions in which it can be decomposed. The decomposition is generally not unique, but the parity is invariant, so that permutations are classified as even or odd according to this number. The signature of a transposition is defined as $-1$, and by extension the signature of a permutation is given by the product of the signatures of its factors. Hence, the signature of even permutations is +1, whereas the sign of odd permutations is -1. The theory of symmetric groups of finite degree establishes that for any permutation $p$ there is a unique decomposition of $p$ as a product of $s(p)$ cyclic permutations (Clarke, 1974):

In turn, any cyclic permutation of $i$ objects can be written as the product of $i-1$ transpositions. Hence, any permutation of $N$ objects can be factorized as $s(p)$ cyclic permutations of $i,j,k,\mathrm{\dots },$ objects, with ${\displaystyle i+j+k+...=N}$. The signature of the permutation, given by the product of the signatures of the cyclic factors is then ${\epsilon }_p={(-1)}^{i-1}{(-1)}^{j-1}{(-1)}^{k-1}\mathrm{\dots }={(-1)}^{N-S(p)}$. Rearranging, the following identity is obtained:

It has been shown that a permutation $p$ corresponds to an L-subgraph $\mathrm{\ell }$, and that the cyclic permutations in $p$ correspond to the cycles in $\mathrm{\ell }$. Thus, replacing the permutation $p$ for $\mathrm{\ell }$ and the number of cyclic permutations in $p$ for the number of cycles in $s(\mathrm{\ell })$ completes the proof of the Coates formula.

The expression of $det[\mathbf{𝐀}]$ in Equation (20) can now be derived strictly from the graphical structure of the associated graph. Indeed, the sign of the contributions of ${\mathrm{\ell }}_1$ and ${\mathrm{\ell }}_3$ are positive because they contain an even number of cycles (four loops in the former case, one loop and one cycle of length 3 in the latter case), whereas ${\mathrm{\ell }}_2$ contains an odd number of cycles (two loops and a cycle of length 2) and accordingly, its contribution is negative.

The Coates formula leads naturally to the definition of the weight of an induced subgraph as the signed sum of the weights of the L-subgraphs contained in it. According to this definition, the weight of the I-subgraph $I_{{\gamma }_k}$ is equal to the determinant of the principal submatrix $A({\gamma }_k)$:

This definition is the last element required to reformulate the stability conditions for a reaction-diffusion system from a graph-theoretical point of view. The coefficient of order $k$ in the characteristic polynomial was expressed in Equations (10–11) as a sum over all the principal minors of order $k$ in ${\mathbf{𝐅}}^{\mathbf{𝐑𝐃}}$. A method to associate a graph $G_I[{\mathbf{𝐅}}^{\mathbf{𝐑𝐃}}]$ to ${\mathbf{𝐅}}^{\mathbf{𝐑𝐃}}$ has been established. Particularly, each $k$ x $k$ principal submatrix ${\mathbf{𝐅}}^{\mathbf{𝐑𝐃}}({\gamma }_{\mathbf{𝐤}})$ corresponds to an induced subgraph $I_{{\gamma }_k}$ of order $k$ in the interaction graph. Furthermore, the associated principal minor $[det{\mathbf{𝐅}}^{\mathbf{𝐑𝐃}}({\gamma }_{\mathbf{𝐤}})]$ is given by the weight $w(I_{{\gamma }_k})$ of the associated I-subgraph in $G_I$. Substitution of this identity restates the algebraic expression of $a_k$ given in Equation (12) as sum of weights of the induced subgraphs as:

where the summation goes over all the I-subgraphs of order $k$ in the interaction graph. Expanding the weight of the I-subgraphs in terms of the L-subgraphs according to Equation (23) leads to

Likewise, substitution of the weights of the I-subgraphs in the reaction graph in Equation (13) leads to the graphical counterpart of the uncoupled expressions of $a_k$:

Examination of these results provides an important insight on the relationship between the feedback structure of the reaction scheme and the stability of the associated reaction-diffusion system. More precisely, the graph-based expressions reveal that every cycle in the network has a defined role in the dynamics of the system. This allows to break down the complete network into smaller functional motifs, thus providing a powerful tool to analyze general networks independently of their complexity or specific parameter values. In the following, we show how this analysis can be applied to the three examples shown in Figure 2.

In the following, we show how the graph-theoretical formalism can be used to analyze the three networks presented in Figure 2. For each network (Appendix 2—figures 5, 6, 7), we present the coefficients of the characteristic polynomial, stability conditions, and instability conditions in terms of cycles. Finally, we highlight the feedback that provides instability and the trade-off between stabilizing and destabilizing feedbacks that underlies the minimum diffusion ratio $d$.

Appendix 2—figure 5

Download asset Open asset

Appendix 2—figure 6

Download asset Open asset

Appendix 2—figure 7

Download asset Open asset

Numerical simulations of the networks identified by our analysis reveal that when a destabilizing cycle comprises only non-diffusible reactants, the reaction-diffusion system does not form a periodic pattern but rather amplifies the noise in the initial conditions. This behavior depends on the specific dispersion relation associated with these systems that shows an asymptotic behavior as the wave-number $q$ increases and leads to an amplification of any of the fluctuations in the initial conditions. Therefore, although these types of reaction-diffusion systems fulfill the Turing instability conditions, they do not amplify a preferential wavelength. This type of behavior has been previously described in reaction-diffusion systems composed of two diffusible reactants and one immobile reactant that activates itself (White and Gilligan, 1998; Klika et al., 2012). Similar behaviors can be observed for indirect non-diffusible auto-activation implemented for example by two immobile nodes that mutually activate or repress each other. In RDNets, these networks are filtered out by default, but they can be re-included in the analysis by selecting the option 'Noise Amplifying Nets'.

The type II network presented in Figure 2 can form a self-organizing pattern even when all mobile reactants have equal diffusivities. Our analysis based on a graph-theoretical formalism demonstrates that this can be achieved by reaction terms that increase the stability of the homogeneous steady state. In particular, the combination of the conditions for homogeneous steady state stability and instability to spatial perturbations shows that the minimum diffusion ratio $d$ required to form a pattern is defined by the ratio between two stabilizing feedbacks ${\text{c}}_1$ and ${\text{c}}_2$:

If the stabilizing feedback ${\text{c}}_1$ is stronger than ${\text{c}}_2$, the system can form self-organizing patterns even with equal diffusivity ($d=1$). These two stabilizing feedbacks correspond to the negative self-regulatory loops of the diffusible nodes v and w. Such negative self-regulatory loops have often been interpreted as decay terms. From the LALI perspective, this suggests the possibility that although the diffusion coefficients of v and w can be equal, their range - i.e. the ratio between diffusion and decay - must be different, such that the decay of the destabilizing node v is greater than the decay of the stabilizing node w:

Within the LALI framework, the network topology suggests that u and v, which mutually activate each other, could behave like the short-range auto-activator, while w could implement the long-range inhibitor. Although v and w have the same diffusivities, w may still behave like a classical long-range inhibitor due to its longer half-life to limit the expansion of the activator. This idea also appears to be consistent with the simulated periodic patterns of v and w that are qualitatively similar to the peaks of an activator and an inhibitor, respectively (Figure 2c and Appendix 6—Table 4). We note, however, that opposite periodic patterns of v and w can be obtained with the same diffusivities and half-lives but using different kinetic parameters to implement the cycle ${\text{c}}_3$ that connects u and w (Appendix 3—figure 2). The LALI concept is therefore not an appropriate framework to describe pattern formation for these networks. Indeed, the LALI concept requires the subjective definition of a subset of nodes that behaves as the activator and a subset of nodes that behave as the inhibitor with the final aim of identifying different effective ranges (Miura, 2007; Korvasova et al., 2015). Such effective ranges need to be defined ad hoc for each network and in cases with equally diffusing reactants appear to be just reformulations of the reaction kinetics that do not contribute to the identification of the general principles that underlie pattern formation.

Appendix 3—figure 2

Download asset Open asset

Other Type II networks are also difficult to relate to classical LALI activator-inhibitor models. For example, the network shown in Appendix 3—figure 3, which has no self-regulatory negative loop on v, allows for a minimum diffusion ratio $d$ that depends on a trade-off between the destabilizing cycle ${\text{c}}_4$ and the stabilizing cycles ${\text{c}}_1$, ${\text{c}}_2$ and ${\text{c}}_3$:

In this case, the two stabilizing terms that define the minimum diffusion ratio $d$ are ${\text{c}}_2^2$ and ${\text{c}}_3$, which correspond to the squared strength of the self-regulatory negative loop of w and to the strength of the negative feedback between v and u, respectively. These terms cannot be related in a straightforward manner to the range of the two diffusible nodes. Within the LALI framework, the network topology suggests that u and v could behave like the short-range auto-activator, while w could behave like the long-range inhibitor. However, numerical simulations show that the system can form periodic patterns that are opposite to those of classical activator-inhibitor models (Appendix 3—figure 3). The LALI framework therefore has little explanatory power to describe the pattern formation processes for Type II networks.

Appendix 3—figure 3

Download asset Open asset

The LALI concept of local auto-activation and long-range lateral inhibition is also difficult to relate to the mechanism that underlies Type III networks. These networks satisfy the instability to spatial perturbations due to their topology and have pattern-forming conditions that only depend on the requirements for homogeneous steady state stability. For example, the Type III network shown in Figure 2 can form a pattern for any diffusion ratio $d$:

These stabilizing terms cannot be related in a straightforward manner to a short activation range and a long inhibition range. In addition, numerical simulations show that depending on the diffusion ratio $d$, which can vary freely, periodic patterns qualitatively opposite to the patterns expected in an activator-inhibitor system can be formed (Figure 2c right and Appendix 3—figure 4a). Moreover, as in the case of Type II networks, the qualitative aspect of the periodic patterns not only changes due to the diffusion ratio $d$ but also depending on the reaction kinetics, e.g. in the different Type III topology shown in Appendix 3—figure 4b.

Appendix 3—figure 4

Download asset Open asset

These observations suggest that the final aspect of self-organizing periodic patterns is not necessarily related with the effective ranges of activators and inhibitors and does not reflect a mechanism that prevents the expansion of activator peaks. We propose instead that the periodic patterns are simply associated with the final growth at which each reactant reaches a dynamic equilibrium, which is determined both by reaction and diffusion terms.

In summary, applying the LALI concept based on differential diffusivity to investigate more complex networks provides little insights into the mechanism that underlies pattern formation. As an alternative, the graph-theoretical formalism presented in this study can be systematically applied to all networks and helps to break down Turing systems into different parts by identifying destabilizing and stabilizing feedbacks. Our analysis highlights that the main role of diffusion is to destabilize the system to spatial perturbations by helping the destabilizing feedbacks to be stronger than the stabilizing feedbacks. In the next section, we show that these findings are consistent with the original reaction-diffusion example described by Turing, where the differential diffusivity is a necessary condition to destabilize the system. This contrasts with interpretations of reaction-diffusion patterning based on LALI models with local auto-activation and lateral inhibition, where the differential diffusivity requirement has been described as a necessary condition to stabilize activator peaks that would be otherwise unstable (Appendix 3—figure 1).

The simultaneous destabilization in opposite directions by differential diffusivity described in the section above is also supported by numerical simulations on more than two cells started with small random deviations from the equilibrium steady state as initial conditions. These simulations do not show the chronological series of events depicted by classical interpretations of the LALI framework based on local auto-activation and lateral inhibition (Appendix 3—figure 1) but rather show the simultaneous formation of activation and inhibition peaks across the whole spatial domain (Appendix 3—figure 8).

Appendix 3—figure 8

Download asset Open asset

The only situation when patterning dynamics in agreement with the LALI interpretation can be observed is when large localized perturbations of the activator are used as initial conditions (Appendix 3—figure 9a). However, these simulations do not only reflect the dynamics of the Turing network but also the dilution and propagation of a localized perturbation, which stimulates the formation of a dissipative soliton (Purwins et al., 2005). This is not the best scenario to investigate the mechanisms that drive pattern formation because the localized perturbation hides the underlying mechanisms that break the symmetry of an initially homogeneous state. The original simulation proposed by Turing (Appendix 3—figure 6) and the one-dimensional simulation shown in Appendix 3—figure 8 suggest that breaking symmetry is not the result of LALI but instead the result of a simultaneous deviation from the equilibrium state in opposite directions due to the differential diffusivity. In agreement with this interpretation, if the magnitude of the localized perturbation is reduced, the simultaneous appearance of inhibitor and activator peaks can be recovered (Appendix 3—figure 9b). The LALI model therefore does not accurately describe the underlying Turing mechanism that can drive the self-organization of periodic patterns in the absence of initial asymmetries.

Appendix 3—figure 9

Download asset Open asset

In the following, we discuss how reaction-diffusion models that contain immobile factors can become unstable to spatial perturbations even when all diffusible reactants have the same diffusivity.

Instability with equally diffusing signals due to immobile reactants

The original example presented by Turing can be modified into a Type III network by introducing an additional reactant Z that participates in a mutual activation with the reactant X and is repressed by Y:

(32) ${\displaystyle \begin{array}{ll}\frac{\mathrm{\partial }\text{X}}{\mathrm{\partial }t}& =5\text{Z}-6\text{Y}+1+d_X{\mathrm{\nabla }}^2\text{X}\\ \frac{\mathrm{\partial }\text{Y}}{\mathrm{\partial }t}& =6\text{X}-7\text{Y}+1+d_Y{\mathrm{\nabla }}^2\text{Y}\\ \frac{\mathrm{\partial }\text{Z}}{\mathrm{\partial }t}& =5\text{X}-6\text{Y}+1+d_Z{\mathrm{\nabla }}^2\text{Z}\end{array}}$

In this network, the auto-catalysis of X is no longer direct as the one in Appendix 3—figure 6a, but it is implemented by the mutual activation with Z (Appendix 3—figure 12a). Similar to the original example of Turing, the system is in equilibrium when $\text{X}=1,\text{Y}=1$ and $\text{Z}=1$. The LALI model does not provide a satisfactory explanation for the pattern formation mechanism in this system, since the inhibitor does not spread faster into the lateral domains of the activator peak due to equal diffusivities (Appendix 3—figure 10 and Appendix 3—figure 11). In the following, we therefore analyze in detail how this modified system is destabilized by amplifying small fluctuations.

Appendix 3—figure 10

Download asset Open asset

Appendix 3—figure 11

Download asset Open asset

We performed numerical simulations of a pair of cells starting from a slight perturbation of the equilibrium state: $\text{X}=1.04,\text{Y}=1.02,\text{Z}=1.06$ in the first cell and $\text{X}=0.96,\text{Y}=0.98,\text{Z}=0.94$ in the second cell. When no diffusion is considered ($d_X=d_Y=d_Z=0$), the system returns to the equilibrium owing to its intrinsic stability (Appendix 3—figure 12b). As in the original example of Turing (see above), the equilibrium is reached by transiently increasing or decreasing the absolute concentrations of X, Y, and Z depending on the initial perturbations, but a reduction of the relative difference between X and Y concentrations (${\mathrm{\Delta }}_c$) brings the system back to equilibrium (Appendix 3—figure 12b). Similarly, if all reactants diffuse equally ($d_X=d_Y=d_Z=1$), the system quickly returns to its equilibrium since diffusion helps to redistribute X, Y, and Z towards homogeneity (Appendix 3—figure 6c).

Appendix 3—figure 12

Download asset Open asset

However, if X and Y diffuse equally ($d_X=d_Y=1$) and Z is immobile ($d_Z=0$), the relative difference between X and Y (${\mathrm{\Delta }}_c$) is progressively increased, and both cells keep on deviating from equilibrium (Appendix 3—figure 6d). In the original example presented by Turing (see above), the deviation from equilibrium was guaranteed by differential diffusivity that implemented a larger flow of Y from the first to the second cell. In contrast, when X and Y have the same diffusivity, the flow of X instead is larger due to the higher concentration gradient of X between the two cells. This leads to a progressive decrease in ${\mathrm{\Delta }}_c$ and therefore to a return to equilibrium (Appendix 3—figure 6c). In the simulation of the three-node network shown in Appendix 3—figure 12d, X and Y have the same diffusivity and consequently a larger flow of X is also observed. However, in this case the immobile reactant Z together with the mobile reactant Y can compensate for the larger flow of X. This is possible because Z is not redistributed by diffusion and together with a small flow of Y can promote and inhibit X in the first and second cell respectively (see first equation in Equation (32) and dashed blue line in Appendix 3—figure 12d). In addition, the small flow of Y form the first to the second cell allows Z to further deviate from equilibrium because of its positive feedback with X (see third equation in Equation [32]), which re-enforces the effect of Z over time and maintains the divergence of the entire system from equilibrium. Immobile reactants therefore fulfill a role as 'capacitors' that can integrate the effect of diffusing reactants to destabilize the reaction-diffusion system and to drive self-organizing pattern formation.

Inspecting the example shown in Appendix 3—figure 12 from a LALI perspective, it could be speculated that the effect of the immobile reactant in Type III networks is just equivalent to a reduction of the effective range of the auto-activation to satisfy the differential diffusivity independently of the activator diffusion rate. This simplification, however, neglects the important role of the specific reaction terms associated with the immobile reactant. Indeed, in the example shown in Appendix 3—figure 12d, the immobile reactant Z does not just reduce the range of activator peaks; rather its dynamics - and in particular its ability to diverge fast from equilibrium - are an integral part of the mechanism that allows the peaks to form in the first place. These dynamics do not only depend on the mutual activation between X and Z but also on the the inhibition of Z by Y, whose diffusion stimulates Z to diverge from equilibrium. Interpreting these dynamics as an effective change in spatial range of the activator appears to provide little insights into the general mechanism that underlie pattern formation. Indeed, in other Type III networks, e.g. more complex networks with more than one immobile reactant, the identification of similar effective ranges requires the definition of alternative ad hoc approximations, which ultimately correspond just to a reformulation of the reaction kinetics of the network.

In conclusion, Type III networks can form a pattern even when the diffusing signals have the same range due to the capacitor effect of immobile reactants and associated reaction terms. Remarkably, this can happen even when the activator diffuses more than the inhibitor (Appendix 3—figure 13). This in an intrinsic property of Type III topologies, where the immobile reactant acts as a buffer to amplify any small advantages or disadvantages over the inhibitor diverging quickly from the equilibrium state. In agreement with this finding, if the fast divergence of the immobile reactant is limited by a negative self-regulation, most networks are of Type II, meaning that their ability to compensate for the range of the mobile reactants depends on the relative speed at which each reactant grows. For example, in the Type II network shown in Appendix 3—figure 3, this can happen only when the growth of the activator is slowed down sufficiently by the inhibitor through the negative cycle ${\text{c}}_3$.

Appendix 3—figure 13

Download asset Open asset

The examples presented in this Appendix highlight that in classical two-reactant Turing models the differential diffusivity destabilizes the equilibrium state by maintaining an imbalance between reactants, which drives a further deviation from equilibrium. Importantly, the reaction terms of the Turing system guarantee that the deviation happens simultaneously above and below the equilibrium state. In the activator-inhibitor model, for example, the differential diffusivity not only gives an advantage to the self-enhancement of the activator but simultaneously to the auto-inhibition of the inhibitor. Therefore, in accordance with a recent proposal (Klika et al., 2012), we suggest that the negative self-regulation of the inhibitor has a more important role than previously assumed.

In agreement with these observations, one-dimensional simulations like the one shown in Appendix 3—figure 8 reveal that the periodic patterns of Turing systems are formed with a simultaneous appearance of activation and inhibition peaks and that the patterning dynamics do not follow the sequence of events described by the LALI mechanism based on local auto-activation and lateral-inhibition. The periodic patterns therefore do not reflect a longer range of the inhibitor to limit the auto-activation. Instead, we propose that the periodic patterns of both the activator and the inhibitor reflect only one range, usually referred to as the wavelength, which is determined by the differential diffusion but also by reaction terms. According to this view, the periodic patterns formed in the activator-inhibitor model reflect different amplitudes of activator and inhibitor levels - rather than a difference in the ranges - that depend both on differential diffusivity and reaction terms.

We therefore propose that the role of immobile reactants in Type II and Type III networks is not to implement an effective difference in the ranges of local auto-activation and long-range inhibition, but rather to help the system to diverge from equilibrium, which is normally achieved by differential diffusivity in classical two-component Turing systems. We find that immobile factors can help to destabilize the system, since they are not subjected to the equilibrating effect of diffusion and therefore fulfill a role as 'capacitors' that can integrate the effect of diffusing reactants to destabilize the reaction-diffusion system by quickly amplifying perturbations.