Widespread biochemical reaction networks enable Turing patterns without imposed feedback

Abstract

Understanding self-organized pattern formation is fundamental to biology. In 1952, Alan Turing proposed a pattern-enabling mechanism in reaction-diffusion systems containing chemical species later conceptualized as activators and inhibitors that are involved in feedback loops. However, identifying pattern-enabling regulatory systems with the concept of feedback loops has been a long-standing challenge. To date, very few pattern-enabling circuits have been discovered experimentally. This is in stark contrast to ubiquitous periodic patterns and symmetry in biology. In this work, we systematically study Turing patterns in 23 elementary biochemical networks without assigning any activator or inhibitor. These mass action models describe post-synthesis interactions applicable to most proteins and RNAs in multicellular organisms. Strikingly, we find ten simple reaction networks capable of generating Turing patterns. While these network models are consistent with Turing's theory mathematically, there is no apparent connection between them and commonly used activator-feedback intuition. Instead, we identify a unifying network motif that enables Turing patterns via regulated degradation pathways with flexible diffusion rate constants of individual molecules. Our work reveals widespread biochemical systems for pattern formation, and it provides an alternative approach to tackle the challenge of identifying pattern-enabling biological systems.

Fig. 1. A comparison of previous approaches and this work.

a Examples of Turing patterns. b Previous approaches of using signed directed graphs (gene regulatory networks) to generate intuition and to identify specific molecular systems governing Turing patterns. c Common biochemical reaction networks based on synthesis, degradation and binding that do not have general correspondence to signed directed graphs. The goal of this study is to identify widespread reaction networks capable of generating Turing patterns.

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⁻⁴ 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.

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.

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.

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 μm² s⁻¹ 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.

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.