A systematic search for Turing pattern enabling reaction networks without imposed feedbackTo address the question of whether common biochemical reactions with no apparent feedback loop can generate Turing patterns, we first enumerated 11 basic types of biochemical “complexes” formed by noncovalent interactions between regulatory molecules such as RNAs and proteins (Fig. 2a, 11 icons with various numbers of circles). These complexes correspond to all topologically distinct configurations with up to four subunits (e.g., protein chains and RNA molecules) (denoted as A, B, C and D or the four types of circles in Fig. 2a. see Methods). We next considered elementary reaction networks that lead to the formation of the complexes via binding. Each complex in Fig. 2a can be viewed as the “final product” of each reaction network and is defined as characteristic complex with respect to the network. In addition to binding and unbinding, each network includes constant synthesis of the unbound molecules, degradation of each molecule in each complex, and diffusions of all chemical species (e.g., Fig. 2a callout). The inclusion of production and degradation is because the timescale of pattern formation process in biology is often comparable to or slower than typical half-lives of RNAs and proteins (see Methods). We described each reaction with mass-action kinetics. The 11 characteristic complexes (Fig. 2a) and their associated reactions (e.g. Fig. 2a callout) are fundamental processes applicable to virtually all biomolecules in cells. Of note, these networks are different from previous graph-based approaches for studying Turing systems in that the networks chosen here describe regulations at post-synthesis stages, which typically correspond to post-transcriptional and post-translational processes widespread in biology, whereas transcriptional networks (i.e. genes turning each other ‘on/off’) and their associated graphs were often used in previous studies21,22,23. Our modeling framework also differ from previous mass-action based models24,26, i.e., we explicitly describe complex formation as an elementary and ubiquitous biochemical reaction, and we do not attempt to produce Turing patterns from minimal systems that can be interpreted directly with the activator-negative-feedback concept. Our goal is to test whether each reaction network leading to the formation of characteristic complexes can produce Turing patterns with biologically plausible parameter values.Fig. 2: A screen for pattern-enabling reaction networks.a 10 complexes with various configurations containing 1–4 subunits (circles). Callout shows an example of reaction networks associated with a characteristic complex. Square boxes show pattern-enabling complexes (see B and C for procedures leading to these conclusions). b Illustration of binding reactions in networks (paths, models) leading to complexes shown in A. Each color of arrows in a network corresponds to one binding reaction. Synthesis, degradation, and diffusion of each subunit are not shown but are included in models. Triangles indicate networks analytically shown to be incapable of producing Turing patterns. Boxes show pattern-enabling networks (see C for procedures leading to these conclusions). Inset boxes display examples of stationary patterns illustrating concentrations of free B at time 500 and a box length of 100, simulated with a temporal step size of 10−4 and a spatial grid size of 1.0. The physical interpretations of time and space are described in Supplementary Information. Color gradients are normalized across the 10 models, but the minimum range (max-min) is 0.26 units of concentration. c A flowchart for screening reaction networks capable of producing Turing patterns (see Methods for details). Source data are provided as a Source Data file.While the 11 complexes allow us to examine biochemical reactions with various complexities, each complex may be formed by multiple sets of binding events (e.g. AA + B → AAB or AB + A → AAB). As a result, 23 distinct networks (i.e. reaction paths) leading to the formation of the 11 complexes were included in our analysis (Fig. 2b. Production and degradation not shown). We built a mathematical model for each network with mass-action kinetics using ordinary differential equations (ODEs). We further extended these models to include Fickian diffusion of all molecules, which gave rise to the partial differential equations (PDEs) (see Methods and Supplementary Information). The models contain 1–11 variables representing concentrations of molecular species.Using analytical approaches such as the Routh-Hurwitz Criterion27,28, we found that five networks cannot generate any Turing pattern with arbitrary positive rate constants (see Methods and Supplementary Information) (Fig. 2b, pink triangles. Supplementary Figs. S1–S3). For each of the remaining 18 reaction paths (models), we used a computational pipeline (Fig. 3c. Supplementary Fig. S4) to search for possible parameter sets that can generate Turing patterns. It has been shown that Turing patterns (stationary periodicity in space) often appear below the instability threshold of Hopf bifurcations at the temporally stable regions, whereas limit cycle oscillations (periodicity in time) are observed above the instability threshold of Hopf bifurcation29,30,31. Albeit Turing patterns can be observed even in the absence of Hopf bifurcation, to enhance the efficiency of the computational pipeline we first considered only the reaction, but not the diffusion, terms in our models, and asked whether the ODE system for each reaction network (path) can generate Hopf bifurcations. We selected 10,000 parameter sets for each ODE system and performed numerical continuation for each set to find Hopf bifurcations. These parameter sets were randomly chosen from biologically plausible ranges covering two orders of magnitude (see Methods. Supplementary Table S1). The Hopf bifurcations served as the “encouraging” bifurcation that allowed us to search for Turing bifurcation efficiently. For reaction paths that produced Hopf bifurcations, we added the diffusion terms back into the models, which yielded systems of PDEs. We randomly sampled diffusion coefficients, and identified Turing pattern-enabling parameter sets by analyzing dispersion relations and simulating the PDEs numerically (see Methods) (Fig. 2c)21,22,32,33,34.Fig. 3: Topological requirements for pattern-enabling reaction networks.a Percentages of parameter sets that produced Hopf bifurcations and those that produced Turing patterns for 10 pattern-enabling networks. Light gray area shows 95% confidence interval of the linear regression line (dark gray). b A unified network motif of 10 pattern-enabling reaction networks. c Topological relationships of all pattern-enabling reaction networks. Numbers in parentheses show multiplicities of molecules with the same configurations. d Distributions of terminal points of the complete set of Turing pattern-enabling dispersion curves (505,125 pattern-enabling sets) for detecting possible existence of two types of dispersion curves (presence of the red points’ positional distribution in the positive range would have implied the fractional occurrence of the Type 2 dispersion curves). Source data are provided as a Source Data file.Among the 18 reaction networks with inconclusive algebraic analysis, we found that 10 of them, corresponding to 6 characteristic complexes, produced Hopf bifurcations with biologically plausible parameters (Fig. 2a, b boxes). While none of these paths contains any imposed negative feedback loop often considered a requirement for oscillation (Hopf bifurcation), our results are consistent with recent studies that showed the abilities for post-transcriptional and post-translational reaction networks to generate oscillation without imposed feedback35,36,37. All 10 PDE models derived from the Hopf bifurcation-enabling reaction networks also produced Turing patterns with some combinations of reaction parameters and diffusion coefficients (Fig. 2b, boxes; Fig. 3a; Supplementary Data 1 and Supplementary Data 2). Overall, for the 10 Turing pattern-enabling models, 2% of the randomly generated ODE parameter sets gave rise to Hopf bifurcations, and 0.13% of parameter sets for the full reaction-diffusion systems produced Turing patterns (Fig. 3a. Supplementary Fig. S6). The latter fraction is comparable to a previous study for quantifying the robustness of activator-inhibitor Turing systems21. Nearly all parameter sets with Hopf bifurcation produced Turing patterns with some combinations of diffusion coefficients (Supplementary Fig. S7), and unsurprisingly, there was a positive correlation between percentages of Hopf bifurcations and Turing patterns across the models (Fig. 3a). Although parameter sets that did not produce Hopf bifurcations in ODE models can also generate Turing patterns in PDE models, the percentage of pattern-enabling sets in this group was very low (<3%). We therefore did not include dispersion analysis for the bifurcation-free ODE set in our main computational pipeline for efficiency considerations (Fig. 2c). Taken together, our results show that a wide range of basic biochemical reactions can produce Turing patterns without any imposed feedback loop (see Discussion for a comparison between explicit and implicit feedback).Network topology required for Turing pattern-enabling reactionsOur high-throughput network and parameter scanning allows us to examine the requirements for the reaction networks to achieve Turing patterns. In this section, we describe requirements in terms of network structures. In the next section, we will discuss the requirements of kinetic rate constants.Each of the 10 pattern-enabling reaction networks contain at least two different molecules (Fig. 2b, black and white circles). Note that in the model, different molecules specifically refer to molecules with different parameters such as degradation rate constants (The parameter values differ by up to 100 folds, e.g. Supplementary Fig. S8. See next section for detailed analysis). Furthermore, the characteristic complex of a pattern-enabling reaction network is either heterotrimers or heterotetramers (Fig. 2b). Most interestingly, we found that a simple network motif is present in each pattern-enabling network and absent in each network that failed to produce a pattern (Fig. 3b). This motif contains a monomer that is sequentially involved in the formation of two heterocomplexes in a ‘feed-forward’ manner (Fig. 3B, Monomer X). Reaction network AAB-2 is the simplest reaction network that contains this motif, while other pattern-enabling networks are connected to this motif via three distinct ways: 1) a smaller pattern-enabling network (AAB-2 or its derivative) is a subnetwork of a larger pattern-enabling network (Fig. 3c, orange arrows); 2) a pattern enabling network is expanded by homodimerization of a molecular species to generate another pattern-enabling network (Fig. 3c, green arrows); and 3) a pattern enabling network is a special case of a general pattern-enabling network that has more types of subunits (Fig. 3c, blue arrows).It should be noted that the networks in our study are hypergraphs (each edge that represents binding involves 3 nodes) rather than regular graphs commonly used for defining feed-forward loops. These two types of graphs have very different mathematical properties. We therefore do not claim that we found pattern-enabling feed-forward loops. Nonetheless, the simple motif shown in Fig. 3b (see later sections for its prevalence in biology) is sufficient for distinguishing the pattern-enabling networks from all other ones used in our search, and these results showed that network topology is important for Turing pattern formation even in the absence of imposed feedback and activator-inhibitor identities.It was shown that there are two types of instability for Turing patterns (Type 1 and Type 2) that can be revealed by dispersion analysis. The systems showing Type 1 variant of Turing instability remain stable at small length scale perturbations and characterized by nonzero diffusion of all the molecules comprising the system38, whereas Type 2 instability is associated to the systems comprising one or more non-diffusive molecules21,22,39. Because all pattern-enabling networks have the same origin in terms of their topologies (Fig. 3b, c), we hypothesized that they are of the same Turing pattern type. We therefore analyzed the dispersion curves of all pattern-enabling parameter sets for all networks. We found that all of them produced Turing pattern via Type 1 instability (Fig. 3d. Supplementary Fig. S5) and the finding reinforces the validity of the current PDE model by assuring the sustenance of continuum hypothesis throughout the parameter space38. This result further supports the uniform fundamental mechanism for pattern formation via the 10 basic reaction networks, and it suggests the importance of diffusion in the patterns that we observed.Kinetic rate constant requirement for Turing pattern-enabling reactionsWith the 10 pattern-enabling reaction networks, we asked whether Turing patterns were observed in a specific range of each parameter. Because the 10 networks have different numbers of parameters, we analyzed three major groups of parameters instead of individual parameters: 1) scaled association constants (denoted by K with subscripts representing complexes) reflecting the binding affinities; 2) scaled degradation rate constants for monomer A, C and D (denoted as χA, χC and χD, which were scaled with respect to monomer B’s degradation rate constant); and 3) relative degradation rate constants (also defined as regulated degradation factors, i.e. RDFs) of each molecules in each complex (α, β, γ, and δ with subscripts representing complexes) with respect to the corresponding monomer. We found that while it is more likely for models with greater binding affinities (K), particularly with dimer formation, to generate patterns, patterns were observed in a wide range of K (Fig. 4a, top). In contrast, lower monomer degradation rates with respect to B markedly enhanced the likelihood of pattern formation (Fig. 4a, middle). Nonetheless, patterns were observed across two orders of magnitudes of these rate constants. In addition, we found that different molecules can have significantly different monomer degradation rate constants in those pattern-enabling sets even if the two molecules (e.g. A and C in the ABC-1 model) have similar reactions in the network (Supplementary Fig. S8). Interestingly, we found more pattern-enabling parameter sets with lower degradation rates in complexes (i.e. RDFs) in dimers, whereas the ability for pattern formation did not depend strongly on other individual RDFs (Fig. 4a, bottom).Fig. 4: Distributions of reaction parameters enabling patterns.a Distributions of three types of parameters in all 10 pattern-enabling models. All sampled parameters (gray) and those enabling patterns (colored) are in scales shown in the left and right axes, respectively. The right diagrams show the meanings of the parameters. b Top chart shows percentages of pattern-enabling parameter sets with randomly chosen values from ranges shown on the left (Basal ranges). Heatmap shows the changes of the percentages when the ranges were perturbed. For each perturbation, the computational pipeline shown in Fig. 2c was rerun to obtain the percentages. c–f distributions of RDFs in two types of complexes in all (c) and two selected models (d–f). g Summary of sensitivities of three types of reaction parameters. Source data are provided as a Source Data file.To interrogate the dependence of pattern formation on the cooperativities of parameters, we perturbed the parameter ranges with various imposed restrictions, and performed the computational searches for Turing patterns with the perturbed ranges (Fig. 4b). We found that binding cooperativity is not required for pattern formation, i.e., all complexes in a model can be formed with the same binding affinity without losing the capacity of pattern formation. There was a simple positive correlation between binding affinity and percentages of pattern-enabling parameter sets in all 10 networks (Fig. 4b). Interestingly, highly positive binding cooperativity (i.e. higher binding affinities in higher order complexes) dramatically decreased the percentages of pattern-enabling parameter sets for most networks (Fig. 4b). In contrast, negative binding affinities generally increased the percentages of pattern-enabling parameter sets even if the overall binding affinities were lower than the unperturbed ones (Fig. 4b). We observed that these perturbations had different extents of impact on the 10 pattern-enabling networks. For example, the AABB-3 network had an exceptionally low sensitivity to the perturbation to positive cooperativity. Similar to most other networks, the AAAB-1 network had significantly higher percentages of pattern-enabling parameter sets upon switching to negative binding cooperativity, whereas for the AAAB-3 network, the pattern-enabling ability was reduced with the same perturbation (Fig. 4b).Although differences among monomer degradation rate constants were not strictly required for pattern formation, equalizing all four monomers’ degradation rate constants resulted in significantly lower fractions of pattern-enabling parameter sets, compared to both unperturbed conditions, as well as a perturbed condition allowing a difference between monomer B’s degradation rate constants and that of other monomers (Fig. 4b). Finally, although mechanisms for complex formation can either protect all molecules from degradation or enhance their degradation for most networks (Fig. 4b), equalizing the RDFs for all molecules (i.e. all subunits degrade with the same rate constant in all complexes) completely abolished the ability of pattern formation for all networks (Fig. 4b).The stark contrast between the perturbation by equalizing RDFs and their broad distributions among pattern-enabling parameter sets (Fig. 4a, b, bottom) suggested that pattern formation could be sensitive to the ratios, rather than individual values, of RDFs within each model. We therefore examined the distributions of trimer-to-dimer ratios of RDFs for A and B (αXXX/αXX and βXXX/βXX) under the unperturbed condition, and we found that the pattern-enabling parameter sets tend to have higher trimer-to-dimer RDF ratio for B (Fig. 4c) when parameter sets of all networks were combined. When we compared two networks containing only A and B molecules, AAB-2 and AAAB-4, we found that the pattern formation is generally correlated with higher trimer-to-dimer RDF ratios for AAB-2 but not AAAB-4 (Fig. 4d, e). However, the pattern-enabling parameter sets for the AAAB-4 network had high tetramer-to-dimer RDF ratios for B compared to the background (Fig. 4f). This shows that even though AAB-2 is a subnetwork of AAAB-4, AAB-2’s requirements of RDFs for pattern formation may be replaced by those of RDFs in higher order complexes. In summary, among all three types of rate constants, pattern formation was most sensitive to ratios of RDFs, and it was significantly sensitive to the individual degradation rate constants of monomers. While network topology of binding reactions is important for pattern formation (Fig. 2b), fraction of pattern-enabling parameter sets is not very sensitive to variations of binding affinities (Fig. 4g).We next asked whether Turing pattern formation requires differences among diffusion coefficients (D) of molecular species in a model. For the AAB-2 model, pattern-enabling parameter sets had wide distributions of diffusion coefficients of A and B monomers. There was only moderate skewness of the distributions of DA and DB (Fig. 5a, top callout), suggesting there was not a strong preference for slow- or fast-diffusing monomers. In contrast, the distributions of DAB and DAAB were highly skewed. Patterns seem to be enabled by high DAB (highly negative skewness) and low DAAB (highly positive skewness) (Fig. 5a, top callout), reflecting a physically plausible scenario of critical slowdown of diffusion upon trimer formation. The preference for slow diffusion of the largest complex, and fast diffusion of the second largest complex was also observed with five other pattern-enabling networks as a general trend (Fig. 5a, top cluster in heatmap). In addition, ordered diffusion coefficients from large, slow complexes to small, fast complexes/monomers were observed in nine out of the ten pattern-enabling networks (Supplementary Fig. S9). Nonetheless, two models (AABB-3 and -4) had the opposite trend: pattern-enabling parameter sets generally had higher diffusion coefficients of the largest complex compared to those of the second largest complex (Fig. 5a, bottom callout and heatmap). In all pattern-enabling models, the monomers had weak skewness of diffusion coefficients’ distributions (Fig. 5a, heatmap), suggesting that these parameters are flexible for generating Turing patterns.Fig. 5: Distributions of diffusion parameters enabling patterns.a Histograms show distributions of diffusion coefficients in two pattern-enabling models. All sampled parameters are shown in gray and those enabling patterns are shown in blue and orange. Heatmap shows the skewness of the distributions of pattern-enabling diffusion coefficients for subunits of 3–4 complexes in all 10 models. b, c The left plots show stationary patterns of free B concentrations from simulations with two representative pattern-enabling parameter sets. Values of dimensionless diffusion coefficients shown in icons (one unit of dimensionless diffusion coefficient corresponds to approximately 4.6 μm2 s−1 with a representative length scale. Detailed description regarding the unit of diffusion coefficients is provided in the Supplementary Information). The right plots show perturbed diffusion coefficients and resultant stationary distributions of free B concentration. Source data are provided as a Source Data file.To test the causal relationship between specific diffusion coefficients and pattern formation, we selected two representative pattern-enabling parameter sets for the AAB-2 model which has four molecular species. In Set I (Fig. 5b, left), the two monomers and the dimer AB have the same diffusion coefficient, and the trimer AAB has a lower diffusion coefficient. In Set II (Fig. 5c, left), Subunit B serves as a diffusion facilitator and the diffusion coefficients of the four molecular species were ordered accordingly. In both cases, when we perturbed the parameter set by equalizing the diffusion coefficients of AAB and AB, we found that the patterns were lost (Fig. 5b, c). We found that although differential diffusivity is required to achieve Turing patterns in our models, the differences among diffusion coefficients were moderate in most pattern-enabling parameter sets (Supplementary Fig. S10, Fig. 5b). In addition, our models do not assume any immobile molecules. Therefore, the requirements on differential diffusion rates in our models are less stringent than those in most conventional models of Turing patterns.Omics-level prediction of potential biochemical complexes enabling Turing patternsWhat are possible biomolecules involved in the 6 pattern-enabling complexes (Fig. 2a boxes)? To address this question, we inferred human protein complexes and mRNA-microRNA complexes as examples of high-order complexes with configurations of subunits specified in Fig. 2a. For protein complexes predicted at the proteomic level, we used data integrated from 15,000 mass spectrometry experiments and other experimental data40. For mRNA-microRNA complexes predicted at the transcriptomic level, we used sequence-based target prediction of microRNAs41 (see Methods). We found that almost 9000 human proteins are involved in complexes that match configurations predicted to enable Turing patterns (Fig. 6a. Blue filled bars and black boxes). Similarly, nearly half of human genes’ mRNA products can be potentially involved in high-order mRNA-microRNA complexes that were predicted to be pattern-enabling (Fig. 6a. Filled red bars). In particular, the pattern-enabling configuration ABC and its associated reaction network ABC-1 has most instances for both proteins and mRNAs (Fig. 6a). The predicted complexes not only include some well-known diffusive molecules crucial for development and regeneration (Fig. 6b)42,43,44,45,46,47,48,49, but also contain a wide range of proteins and RNAs that can enable new hypotheses of mechanisms underlying Turing pattern formation (Supplementary Data 3). Due to the lack of measured rate constants in these predicted systems, we were not able to further constrain the list of gene products that enable Turing patterns with experimental data. This qualitative inference of complexes is therefore not sufficient to identify pattern-enabling molecules definitively, and the pool of complex-associated gene products should be viewed as a preliminary set for subsequent screening in future studies. It should be noted, however, that this limitation also applies to the classical activator-inhibitor paradigm for identifying systems for Turing patterns. Our findings on the broad range of proteins and RNAs whose network topologies permit Turing patterns suggest the need of measurements for key parameters, such as degradation rate constants and diffusion coefficients as shown in the previous section, for predicting Turing pattern-enabling systems more accurately.Fig. 6: Omics-level estimation of instances of pattern-enabling complexes.a Numbers of genes whose products are involved in protein (blue) and mRNA-microRNA (red) complexes with indicated configurations inferred from experimental data of protein interaction and sequence complementary respectively. Each bar has a maximum of 20,000, and the filled portion indicates the estimates for the configuration shown on the left. b Examples of protein and mRNA-microRNA complexes that potentially enable patterns. Connecting lines in the protein complex icon indicated experimentally supported pairwise physical interaction. Known biological functions of the proteins and mRNAs are annotated. c Functional enrichment of secreted and membrane proteins involved in high-order configurations that potentially enable patterns (boxes in a). Source data are provided as a Source Data file.Despite the limitation, we asked whether our prediction provides a meaningful list of proteins and RNAs that are involved in biological processes related to pattern formation. We performed functional enrichment analysis for membrane and secreted proteins which are more likely to be transported between cells50. We found that with a background of these proteins from secretome and membrane proteome51, those involved in high-order, pattern-enabling complexes are enriched with functional terms for a wide range of developmental processes (Fig. 6c. See Methods). The mRNAs involved in high-order mRNA-microRNA complexes are also involved in many pattern-formation processes, although the folds of enrichment are moderate due to the large number of mRNAs in the gene set (Supplementary Data 4). These results suggest that pattern formation in development may be a function of proteins in high-order complexes even in the absence of explicit feedback loops of gene regulation. Nonetheless, it should be noted that there are many biological scenarios underpinned by functionally important Turing patterns besides tissue level patterning in development. The length scales of our mathematical models cover both intracellular and tissue-level patterns, so the potential biological functions arising from the complex-driven Turing patterns can be very broad (see Supplementary Information and Discussion).
