Identifying network topologies that can generate turing pattern

Journal of Theoretical Biology 408 (2016) 88–96

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Identifying network topologies that can generate turing pattern M. Mocarlo Zheng a,1, Bin Shao a,b,1, Qi Ouyang a,b,n a b

The State Key Laboratory for Artificial Microstructures and Mesoscopic Physics, School of Physics, Peking University, Beijing 100871, China The Center for Quantitative Biology and Peking-Tsinghua Center for Life Sciences, Peking University, Beijing 100871, China

H I G H L I G H T S


A framework is developed to systematic identify networks that are capable for Turing pattern formation.
We enumerated all possible 2, 3-node genetic regulatory networks to access their ability to generate Turing instability.
Topologies with fewer parameter constraints are further analyzed to uncover the design principles of robust networks.

art ic l e i nf o

a b s t r a c t

Article history: Received 15 March 2016 Received in revised form 15 July 2016 Accepted 8 August 2016 Available online 9 August 2016

Turing pattern provides a paradigm of non-equilibrium self-organization in reaction-diffusion systems. On the basis of many mathematical studies, it has been proposed that various biological development processes use Turing instability to achieve periodic patterns. In this paper, we introduce a framework to systematic identify network topologies that are capable for Turing pattern formation. All possible 2, 3-node genetic regulatory networks are enumerated and linear stability analysis is applied to access their ability to generate Turing instability. We find that all 3-node networks that can achieve Turing pattern can be mapped to either pure or cross activator-inhibitor mechanisms, and the pure activator-inhibitor system is more robust for Turing pattern formation than the other one. Additional linkages can further increase the performance of the circuit by either introducing another core topology or complementing existing regulations. Moreover, we find that addition of a fixed node enables the formation of Turing pattern even when the diffusion coefficients of two morphogens are fairly close to each other. Our results provide the design principle of robust circuits for Turing pattern generation and can be further applied for systematically exploring other bifurcation phenomena. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Turing pattern Network enumeration Nonlinear dynamic analysis Robustness

1. Introduction A most fascinating aspect in biological development is the emergence of complex spatial pattern that composed of highly organized differential cells. Under the control of molecular and cell-cell interaction networks, differentiation coordinates with growth and gives rise to the development of patterns. Different biological pattern formation processes have been resolved in the molecular level in the past decades. For segment formation in Drosophila, the pattern is controlled by the gradient of the signal called morphogen. Maternal secreted factors such as Bicoid are enriched at one end of embryo and encodes position information along the ante-posterior axis (Driever and Nusslein-Volhard, 1988; Little et al., 2011). In the somitogenesis process of zebrafish, a clock n

Corresponding author. E-mail address: [email protected] (Q. Ouyang). 1 These authors contributed equally to this work.

http://dx.doi.org/10.1016/j.jtbi.2016.08.005 0022-5193/& 2016 Elsevier Ltd. All rights reserved.

and wave front model incorporating gene oscillation and cell growth is used to generate vertebrate segmentation (Holley et al., 2000; Oates et al., 2012). Previous work has shown that scaling of mouse mesoderm patterning is precisely determined by the phase gradient of underlying oscillation of signaling proteins (Lauschke et al., 2013). Despite using local information of morphogen or molecular clock, cells can also use reaction-diffusion (RD) mechanisms to achieve spontaneous pattern formation. Patterns emerged from RD system range from complex pattern of cell colony to sea shell patterns and strip and spots in the animal skin (Mimura et al., 2000; Meinhardt, 2003; Kondo and Miura, 2010). Among all the patterns produced by RD systems, Turing pattern draws special attention since proposed in 1950s (Turing, 1952). A typical Turing pattern system includes two diffusing components, one ‘activator’ and one ‘inhibitor’. Self-activation of activator and lateral inhibition helps generate spontaneous pattern in the system. After Turing's pioneering work, a general mathematical framework has been developed to investigate the Turing instability of regular

M.M. Zheng et al. / Journal of Theoretical Biology 408 (2016) 88–96

cellular networks or artificial random networks in which interacting species can diffuse across the links (Othmer and Scriven, 1971; Othmer and Scriven, 1974; Nakao and Mikhailov, 2010). For simple two node regulatory networks, nonlinear dynamics analysis reveals that two mechanisms can be responsible for Turing pattern formation: the pure activator-inhibitor model and the cross activatorinhibitor model (Dillon et al., 1994; Murray, 2003). It also has been suggested that Turing pattern can be generated by self-enhancement in conjunction with long-range inhibition (Meinhardt and Gierer, 1974; Meinhardt and Gierer, 2000). With the development of genetic and molecular techniques over the last decades, it is possible to identify examples of Turing pattern in living systems. For example, researchers have shown that interaction of pigment cells of zebrafish fulfill the condition of Turing pattern, which is further proved by experiments (Nakamasu et al., 2009; Bullara and De Decker, 2015). Digital pattern of mouse is another example in which the dose of Hox gene determines the wave length of the pattern (Sheth et al., 2012; Miura, 2013; Uzkudun et al., 2015). Although various genetic mutant can produce abnormal number of digits, one critical parameter that can adjust wave length (digit number) is firstly discovered which is in accordance with a Turing-type mechanism. Despite the development of mathematical models and increasing evidence for Turing pattern in biological system, a systematic understanding of the robust regulatory networks that are capable for Turing pattern formation is still not categorized. One of the most useful methods to identify topologies that can achieve desired function is the enumeration of networks. Full range of possible topologies is explored and functional topologies with fewer parameter constraints are selected. Current examples include adaptation (Ma et al., 2009), dose-response alignment (Long et al., 2012), self-organizing cell polarization (Chau et al., 2012), and segment polarity in development (Ma et al., 2006). However, the scheme for network selection only uses dynamics information of the variables. It cannot discover networks that can achieve particular bifurcation behavior, such as Turing bifurcation. Instead of searching for particular dynamics, we develop a framework for searching networks that are capable for certain bifurcation behavior. In this paper, we carry out an exhaustive computational analysis of 2,3-node genetic regulatory networks that can generate Turing pattern. Clustering analysis and classification of the simplest networks (defined as the simplest ‘irreducible’ topology) are further applied to identify the structural and mechanical characteristics of the robust networks. We find that all the robust 3-node network can be alternatively classified into the pure activator-inhibitor model or the cross activator-inhibitor model, respectively. Notably, the classical activator-inhibitor mechanism is much more robust than the other mechanism and underlies most of the robust networks. On the basis of two node network, appropriate involvement of a third node could further increase the robustness for both mechanisms and allows the two interacting morphogens diffuse at similar rates.

2. Results

89

control mesendodermal gene expressions and also contribute to an activator-inhibitor system (Muller et al., 2012). In our simulation, we enumerated all 2-node transcriptional networks and 3-node networks with 2 morphogens (3N2M) and 3 morphogens (3N3M) (Fig. 1A, left). We modeled the networks via partial differential equations (Fig. 1A, right). The steady state of the network under randomly sampled parameter set was obtained by numerical solving the nonlinear differential equations. Using the values of three variables at the steady state, we can calculate the Jacobi matrix of the circuit and get the relationship between the largest real eigenvalue max(Re(λ)) and the wave number k. The criterion of Turing instability is that there exists a real positive eigenvalue under some mode k (Fig. 1B). We used Q value to quantify the robustness of each network, which is defined as the number of parameter sets that can achieve Turing instability. A larger Q value indicates that the network can generate Turing pattern over a wider range of parameters (i.e. more robust). To better understand a circuit's underlying mechanism of Turing instability, we introduced the concept ‘Jacobian topology’ (Fig. 1B, dashed box), which corresponds to the sign of partial derivatives in the Jacobian matrix at the steady state. Through computational enumeration of 2-node networks ( 32 × 2 = 81 in all), we found four networks with Q 40 and they can be classified into two mechanisms according to their Jacobian topologies, where the first one is the activator-inhibitor Jacobian topology (AIJT) and the two morphogens have the same phase in the self-organized pattern (Fig. 2A, left, assume D
o Do). The second mechanism is cross activator-inhibitor model (Dillon et al., 1994), in which one node activates itself and inhibit the other node (Fig. 2B, right) and the two morphogens generate anti-phased pattern. We denote it as cross activator-inhibitor Jacobian Topology (cAIJT) in our paper. In the self-organized patterns, the spatial phase difference between the components is ‘0′ or ‘π’, which makes it possible to map a 3-node network into a 2-node Jacobian topology according to the phase difference (as well as the type of mutual regulations) between its morphogens. A related issue that should also be mentioned is the scale of the parameters which will affect the length scale of the pattern we investigated. The diffusion coefficient of proteins depends on spatial environments and varies across different organisms, thus may contribute to the patterns of different scales (Gregor et al., 2005). In the Nodal–Lefty model, the growth factors diffuse extracellularly and the diffusion coefficient ranges from 1 to 19 mm2/s (Muller et al., 2012). Thus we can set the unit of diffusion coefficient to 10 mm2/s in the simulation. The reaction time scale of the system is about 104 s (estimated from the clearance rate). Then starting from a formula about the length scale of the pattern: λ~ DT (Quyang, 2010), we can estimate that the resulting pattern may have a wave length around 102 μm, which is comparable to the size of the blastula embryo. The intracellular diffusion rate of proteins is smaller as compared to the diffusing signal. For example, the diffusion coefficient of protein Bicoid in Drosophila cell cytosol is estimated to be 0.3 mm2/s (Gregor et al., 2007). The smaller diffusion rate, combined with a similar time scale of transcriptional activation, can lead to patterns with smaller length scale which may underlies the development or cell polarity.

2.1. Screening procedures and two-node Jacobian topologies 2.2. Three node networks with two morphogens We used genetic regulatory networks to model the pattern formation process. Transcriptional regulation is crucial in development. In Drosphila embryo, maternal factor Bicoid is a transcription factor which controls expression of the gap genes in thousands of nuclei; over five hundred transcription factors (TFs) in organ development of Drosphila have spatially restricted expression (Hammonds et al., 2013). In sea urchin and mouse, Nodal and Lefty are two growth factors that diffuse freely in space and

In the exhaustive computational analysis of three node networks, node A and node B were assigned to be morphogen (DA, DB 4 0) and node C was assigned to be a fixed regulator (DC ¼ 0). We enumerated all 6162 valid 3-node networks and found 651 networks with a nonzero Q value. The distribution of Q value is shown in Fig. 2B. For convenience of further analysis, we swapped node A and node B to ensure that every network is capable of

90

M.M. Zheng et al. / Journal of Theoretical Biology 408 (2016) 88–96

Fig. 1. Scheme of screening for Turing pattern circuits. (A) All possible links for 3-node networks (left) and general RD equations in our research (right). Note the reaction part in the equations includes a regulatory term relating to the network structure and a dilution term. (B) Analysis procedure for a given network. Random parameters were generated using Latin hypercube sampling and linear stability analysis was applied to determine whether the system is capable for Turing pattern formation. We use edges with arrowhead ends to denote positive regulations and edges with bar ends to denote negative regulations. Note that the total self-regulatory effect of node C in Jacobian topology depends on both the non-linear self-activation term and the dilution term.

Turing instability under the condition of DA oDB. There are 59 robust networks with QZ3. We found three dominant edges among these networks, including node A's self-activation, activation from node A to B, and inhibition from node B to A (Fig. 2C, left). These links exactly represent AIJT, considering that node B's self-dilution contributes to a negative term in the Jacobian matrix. The 4 networks with the highest Q value (¼ 5) are shown in Fig. 2C. These networks all belong to the category of AIJT. Robustness of network i, ii and iii is enhanced by the mutual activation between node A and C; robustness of network iv results from the self-inhibition of node B. Furthermore, network i's robustness is further increased by its compatibility of opposite diffusion-coefficient combinations (i.e. DA o DB and DA 4DB ) (see

Discussion). Our results reveal that robust networks tend to include the classical AIJT, and appropriate coupling of a third node C could notably improve the robustness of AIJT. The 651 robust networks can be roughly categorized into five clusters, and the most common edges for each cluster are presented in Fig. 2D. Cluster 1 can be classified as cAIJT, while cluster 2–5 can be classified as AIJT. In those clusters corresponding to AIJT, node C either mediates the inhibition from B to A (cluster 2), or forms a feedback loop with node A (cluster 4 and cluster 5). We also picked out all the simplest networks and classify them into cAIJT and AIJT (Fig. 2E). Among the simplest networks of cAIJT, the fixed node C mainly plays three possible roles in the generation of anti-phased pattern: enhancing the self-activation of morphogen

M.M. Zheng et al. / Journal of Theoretical Biology 408 (2016) 88–96

91

Fig. 2. Two Jacobian Topologies and results of three node networks with two morphogens. (A) Jacobian topologies of AIJT and cAIJT. (B) Q value distribution of robust networks. (C) Most common edges in networks with Q ¼ 3 (left) and the most robust networks with Q¼ 5 (right). (D) Clustering results of topologies of robust networks. Nodes of the networks have been arranged to ensure that the network is capable for pattern formation in the condition of DC ¼0 and DA o DB before analyzing. In the cluster diagram, red color is for activation, green for inhibition and black for no regulation, respectively. Edges of each column are presented in the lower part. (E) The simplest networks in the 3N2M enumeration. Brown dashed box, node C acts as a buffered inhibition node for B; blue dashed box, node C helps generate a self-activation of the slower morphogen A; yellow dashed box, node C helps transmit regulations between morphogen A and B. Each simplest network either belongs to cAIJT or AIJT.

A, enhancing the self-inhibition of morphogen B, and transmitting interactions between morphogen A and B. These regulations via node C complement to the existing regulation and thus increase the robustness of the corresponding network. Among the simplest networks of AIJT, there are at least two roles for the fixed node C, carrying out A's self-activation and transmitting interactions between node A and B. 2.3. Three node networks with three morphogens We next enumerated all the 1934 valid networks in the 3N3M enumeration. The robustness of the networks is quantitatively accessed with the same 100,000 sets of random parameters. We found 475 networks that can achieve Turing instability, the corresponding Q value distribution is shown in Fig. 3A. Still, we swapped the status of nodes in some networks to ensure that every network is capable for Turing instability under the condition of DC oDA oDB . There are 74 robust networks with Q Z 7, and clustering analysis shows three dominant edges among these

networks, including node C's self-activation, activation from node C to B, and inhibition from node B to C (Fig. 3B, left), also representing AIJT. The four most robust networks (Q ¼ 12) are shown in Fig. 3B, right. Networks i and ii contain two 2-node classical AIJT, where both node C and A and node C and B form the AIJT. Network iii is a AIJT composed of node C and B, while node A helps enhance node C's self-activation and the inhibition from node B to node C. Network iv is capable for Turing instability under different combinations of diffusing components, node C and A or node B and A can act as the classical AIJT. Overall, classical AIJT is still the most robust topology for generating Turing pattern in 3N3M system. Networks with higher Q value tend to include more than one classical AIJT. These 475 Turing networks can also be roughly clustered into 5 clusters, and the most common edges of each cluster is drawn in Fig. 4C. We found cluster 1 can be classified into cAIJT, while cluster 2–5 can be classified into AIJT. In cluster 1, node C acts as the cross activator and node B acts as the cross inhibitor. The most common edges of cluster 2–4 all contribute to classical 2-node

92

M.M. Zheng et al. / Journal of Theoretical Biology 408 (2016) 88–96

Fig. 3. Results of three node networks with three morphogens. (A) Q value distribution of the robust networks. (B) The most common edges in networks with Q¼ 7 (left) and the most robust networks with Q ¼12 (right). (C) Clustering results of robust networks. Nodes of the networks have been arranged to ensure that the network is capable for pattern formation in the condition of DC o DA o DB before analyzing. In the cluster diagram, red color is for activation, green for inhibition and black for no regulation, respectively. (D) The simplest networks in the 3N3M enumeration. Blue box, node C and node B forms a cAIJT. yellow box, node C acts as the activator and node B acts as the inhibitor; pink box, node C acts as the activator and node A acts as the inhibitor; green box, node A acts as the activator and node B acts as the inhibitor.

Fig. 4. Design principles of robustness and multiple underlying mechanisms. (A) Examples of an appropriate implement of additional modules. (B) Addition of an edge can either enhance or reduce robustness of the network, depending on its relationship with core topology.

M.M. Zheng et al. / Journal of Theoretical Biology 408 (2016) 88–96

AIJTs with an activator and an inhibitor. In cluster 5, node A acts as an activator in combination with node C and node B is the inhibitor. Our results reveal that AIJT still dominates the underlying mechanism for Turing pattern generation in 3N3M systems. Fig. 4D provides all the simplest networks and classifies them into cAIJT and AIJT. Among all the simplest networks of cAIJT, node A seems to act as a buffer for the inhibition upon node B or help enhance the self-inhibition of node B. Among all the simplest networks of AIJT, there are three combinations in all, AIJT consists of node C and B, AIJT consists of node C and A, and AIJT consists of node A and B.

3. Discussion 3.1. Design principles of robust networks for turing pattern formation Various works have been proposed to design Turing pattern circuit in different biological context. For example, it has been shown that Turing pattern can be generated by oscillator-driven gene networks (Hsia et al., 2012), which may provide guidance for engineering the functional circuits. A working model using the steps of cell metabolism is proposed to produce Turing pattern with a length scale smaller than the cell size (Strier and Ponce Dawson, 2007). However, these studies usually focus on a small number of networks. Only until recently, a high-throughput framework that yield complete solutions for conditions of pattern formation has been proposed to analysis Turing networks (Marcon et al., 2016). In contrast to their approach based on algebra calculation, we here applied a network enumeration approach and provide a systematic investigation of networks capable for Turing pattern formation. All the 3-node networks that are capable for our target function can be classified into two mechanisms. Based on our analysis of the robust 3-node networks, we propose two principles to be considered when designing robust regulatory networks to generate Turing pattern. Firstly, try using as many AIJT topologies in the network as possible. 3of 4 most robust networks of 3N3M enumeration include more than one AIJT core topology. We analyzed all three node networks containing two AIJT core topologies and found that their mean Q value is 4.94, significantly larger than the networks with one AIJT core topology which has a mean Q value of 2.86. The mean Q value of the networks without any AIJT core topology is smaller than 0.1 (Table 1). This analysis indicates that incorporating multiple AIJT core topologies can greatly enhance the robustness of the Turing network. The calculation is based on the enumeration results of 1934 valid 3N3M networks. The AIJT was identified by matching 3 direct regulations ‘A-4A', ‘A- 4I' and ‘I-|A' among all edges, where A (activator) and I (inhibitor) are the selected two nodes within the network. We try all possible A32 =6 selections for each network and count the matching frequency. Avg, average; Std, standard deviation. Secondly, try to add complementary regulations based on the existing core topology. For the left network in Fig. 4A, implementation of an activation edge from node A to node C helps Table 1 Correlation between Q value and number of AIJT cores. AIJT number

Avg of Q value

Std of Q value

Sum of networks

0 1 2 3

0.09 2.86 4.94 4.50

0.46 2.62 3.55 2.89

1426 424 80 4

93

form a feedback loop between node A and node C, thus contributing to a backup of the self-activation of node A and making the network more robust. In another example, node C mediates inhibition from node B to A, as shown in Fig. 4A, right. The addition of linkages to an existing network can be divided into two entirely different categories. If the addition of a new linkage introduces a new core topology to the original network or adds coherent regulations to the presenting core topology without contradicting any present core topology, it often leads to a beneficial effect. If the additional linkage plays a contradiction role of regulations among core nodes, such addition may lead to a neutral or even deteriorative effect to robustness. For the left network in Fig. 4B, addition of negative regulation from node C to node A helps form a double-inhibition loop and enhance the self-activation of node C, thus increasing the robustness of the network. For the right network in Fig. 4B, the same edge also contributes to the self-activation of node C, yet deteriorates the activation of node B by node C. Therefore, it is a neutral addition to the original network. It also explains why the mean Q value of networks with three AIJTs is even smaller than the networks with two AIJTs (Table 1) since there may be contradiction of regulations when trying cramming three AIJTs into one 3-node network. Similar concept has been termed ‘regulation entropy’ in former works which is used to quantify the coherence of regulations (Wu et al., 2009). 3.2. Relationship between the diffusion coefficients of the two morphogens In the analysis of 3N2M networks, we have mentioned that the involvement of a third node may relax the requirement of unilateral inequality for the diffusion coefficients in classic activatorinhibitor system (DA o DB), which coincides with results proposed in a recent study (Marcon et al., 2016). Based on our results of 3-node transcriptional networks, we found two mechanisms underlying this phenomenon. In the first case, the networks contain both cAIJT and AIJT and it is possible to use the same network to generate either co-phased or anti-phased Turing pattern by adjusting parameters in the dynamic functions. One example is illustrated in Fig. 5A. If the coupling of A and C is strong while selfactivation of B is weak, chances are that node A gets a positive selfregulation while node B get a negative self-regulation in the Jacobian matrix (taking the dilution term into consideration), and we achieve a cAIJT-kind Turing instability. However, if the coupling of A and C is weak while self-activation of B is strong, we get an AIJT-kind Turing instability. In the second mechanism, the network only contains one cAIJT with node A acts as the activator by forming a feedback loop with node C. Interactions between A and C lower the efficient diffusion coefficient of the activator, thus allowing DA 4 DB. In this case, the system always generates co-phased Turing Pattern. One example is demonstrated in Fig. 5A, in which node A and C activates each other. Cooperated diffusion may exist in the interactions between membrane and cytoplasmic proteins. For example, in the polarity process of budding yeast, membrane protein Cdc42 can recruit the cytoplasmic Bem1p complex to the membrane, which enables amplification of small GTP-Cdc42 clusters (Howell et al., 2012). Such mechanism relaxes the requirement of the difference between diffusion rates of proteins and may help to form Turing pattern in larger parameter space. Many previous simulations of Turing pattern formation in biology are based on the typical activator-inhibitor RD system. In most (if not all) cases, the ratio of diffusion coefficients has to be very large (DB/DA 410) for the instability to occur (Othmer and Scriven, 1971; Asai et al., 1999), which may not be the real case because seldom is there a huge difference between the diffusion

94

M.M. Zheng et al. / Journal of Theoretical Biology 408 (2016) 88–96

Fig. 5. Relationship between the Two Diffusion Coefficients. (A) Two mechanisms allowing the compatibility of different diffusion-coefficient combinations. (B) Up, distributions of diffusion-coefficient ratio in 3N2M enumeration. Ratio of diffusion coefficients is generated by dividing the larger diffusion coefficient by the smaller one; bottom, three simplest networks whose diffusion-coefficient ratio is close to 1.

coefficients of two biological morphogens (proteins). Here we further investigate to what extend the involvement of a fixed regulator can narrow down the difference between the diffusion coefficients of two morphogens. All valid 3N2M circuits are analyzed and the corresponding distributions of diffusion-coefficient ratio are plotted in Fig. 5B. To be specific, the involvement of node C is defined as connecting to node A and B with last least one input and one output linkage. Among the circuits without the fixed regulator, the minimal ratio for Turing instability is around 16; when a third fixed node is coupled, however, we could find the ratio can be close to 1. We further selected out three typical circuits with close diffusion-coefficient for demonstration (Fig. 5B, bottom). The effective diffusion rate of node A in network i is reduced owing to its inter-inhibitions with the fixed node. Similarly, in network ii, node B's diffusion is rescaled and cAIJT underlies the pattern formation in this case. The positive feedback between node A and C enhances the robustness of network iii as AIJT. Our results reveal that when a “fixed” node is included, diffusion coefficients of morphogens may need not to vary enormously from each other and thus makes it a more plausible explanation for Turing pattern formation in a biological context. In conclusion, we have demonstrated an enumeration framework for identifying multiple-node network topologies that can achieve Turing pattern. We find that there are two mechanisms which can lead to Turing instability in transcriptional network, in which the classical activator-inhibitor system is more robust for Turing pattern generation. Cluster analysis and identification of the simplest networks are conducted to better characterize the robust networks. Two Design principles of robust network have been proposed, including incorporating as many 2-node core topologies as possible and adding complementary regulations to the core topology. We also show that incorporation of a fixed node can relax the requirement of difference between the diffusion coefficients and thus may be responsible in biological system. Our results of robust networks provide guidance for deciphering possible networks response for Turing pattern in biology and our framework proposed here can be extended for solving more bifurcation problems systematically.

4. Materials and methods 4.1. Experimental procedures for enumeration of networks 4.1.1. Enumeration of 3-node networks For three node directed networks, there are nine possible links (including three self-regulations). Each link has three possibilities: positive regulation, negative regulation, or no regulation. Thus there are 39 ¼19,683 possible networks in total. But we can in fact eliminate many of them by checking the symmetry and connectivity of the networks.


Symmetrical characteristic. In the 3N2M enumeration, there is




a mirror symmetry between node A and node B, so nearly 1/2 of the networks can be eliminated. In the 3N3M enumeration, there is a mirror symmetry and a rotational symmetry among all 3 nodes, so nearly 5/6 of the possible networks can be eliminated. The symmetric networks were deleted by application of group theory. Taking the 3N3M case for example, we defined 6 operators (including identity operator E). We enumerated all the 3-node networks and for every selected network, we applied 5 operators on it (except E), and deleted the resulting 5 networks from the set of networks. Connectivity of the network. We eliminate those networks in which at least one node has no input or output (except selfregulation) edges. Additionally, we don’t take the connectivity of the fixed node C into consideration in the 3N2M enumeration in order to cover all the 2-node networks at the same time.

Finally, we get 6162 networks out of 19,683 for further screening in the 3N2M enumeration and 1934 in our 3N3M enumeration. 4.1.2. Reaction-diffusion equations Our equations are based on the transcriptional regulatory model. We use Hill functions for the regulation terms, which is used ubiquitously in biological models. The induction terms are added together and the inhibition terms of proteins take a multiplication form, suggesting that the inhibitors are independent of

M.M. Zheng et al. / Journal of Theoretical Biology 408 (2016) 88–96

one another and can independently block the expression of the target protein. Then, the model equation for species X reads:

⎛ positive ⎞ ∂X 1 ⎟ = bX + hX ⋅⎜⎜ ∑ ni ⎟ ∂t ⎝ i ∈ IX 1 + ( KYi X /Yi ) ⎠ negative

⋅

∏ j ∈ JX

1

(

1 + Y j/KY j X

nj

)

a11 − k2Dx − λ k

a12

a13

a21

a22 − k2Dy − λ k

a23

a32

k2Dz

a31

− γX⋅X + DX ⋅∇2X (1)

Here X,Y¼A,B,C; IX and JY indicate the sets of the indexes of the nodes which have positive and negative regulations upon node X, respectively. n is the Hill index, γ refers to self-degradation rate and D refers to diffusion coefficient. Due to the non-dimensionalization, we have KAA, KBB, KCC, DA ¼1. And DC equals to zero for all the 2-morphogen systems. The basal expression of proteins b (equals to 0.1) is also added in our equations. In the Hill functions, we set all Hill indexes to 2. We further demonstrate that our results are insensitive to n as long as it is larger than 1 (see Text 1 in the Supplementary Material). 4.1.3. Sampling the parameter space We use Latin hypercube sampling method (Iman et al., 2006) to evenly sample the high-dimensional parameter space. There are maximumly 14 parameters included in the 3N3M sampling and 13parameters in the 3N2M sampling (DC≡0). Moreover, the distribution for each parameter is uniform in logarithmic coordinates within the range: [0.01, 100]. Here we randomly generate 30,000 sets of parameters for the 3N2M enumeration and 100,000 for the 3N3M enumeration, and then test every valid network using all these parameter-vectors. In total, There are 6162  30,000E1.8  108 combinations left to be further analyzed in the 3N2M enumeration and 1934  100,000E1.9  108 combinations in the 3N3M enumeration. 4.1.4. Criteria of Turing-instability There are two steps to check for the Turing instability, getting the nonzero steady state and finding out an unstable mode.


Get a nonzero stable state. We use Newton–Raphson method to

95

a33 −

=0 − λk

(3)

where aij is the relative partial derivatives at steady state. If there is a k that makes Re ⎡⎣ max ( λk ) ⎤⎦ > 0 while Im ⎡⎣ max ( λk ) ⎤⎦ = 0, then this circuit meets the criteria of Turing instability.

4.2. Simulation of turing pattern A customary DR-simulation program is used to test the validity of our Turing instability criteria and network classification. The code is written in MATLAB. The simulation is done within a 256  256 lattice and diffusion is allowed within 4 closest neighbor lattice points in each step. No flow boundary condition is applied. Simulation results are presented in the Supplementary materials. 4.3. Clustering of networks There are nine possible links for a network. For every network, we assign a value to each of the nine links: ‘1′ for positive regulation, ‘ 1′ for negative regulation, and ‘0′ for no regulation. The Hamming distance, which is the number of links that differ in the two networks, is used to cluster the networks. 4.4. Selecting the simplest 3-node networks Firstly, we pick out all robust networks with no more than 6 edges. Secondly, we screen out those networks that can be derived by adding edges to another robust network with less edges, thus we get all the ’irreducible’ networks. Thirdly, we leave out those networks in which node C has only input or output edge in the 3N2M enumeration since these are not actually 3-node networks.

solve the root for the ODE part of Eq.1: ⎛ positive 1 0 = bX + hX ⋅⎜ ∑ ⎜ ⎝ i ∈ IX 1 + KYi X /Yi

(




ni

)

⎞ ⎟⋅ ⎟ ⎠

negative

∏ j ∈ JX

Acknowledgments

1

(

1 + Y j/ KY j X

nj

)

− γX⋅X

(2)

where the result converges quickly after several steps. Solver from GSL (Gough, 2009) is used to accelerate the multidimensional root-finding process, which will stop when the residual 3 value f = ∑i = 1 fi is smaller than 10  8. Using the values of the steady state, all partial derivatives of the Jacobian matrix in Fig. 1B are analytically solved. Then we check whether the fixed point is positive and whether all real parts of its three eigenvalues of Jacobi matrix are less than zero, which means that the fixed point is stable. Any fixed point that doesn't fulfill the above conditions is abolished. In addition, we sample 64 different initial values in the 3-dimensional variable space ranging [0.01, 100] for each network, considering that there may be more than one stable point. A stable point with Euclidean distances more than 10  5 to all other stable points are regarded as a valid one, and all valid points are then tested in the next scanning step for unstable mode. Find out an unstable mode. The criterion of Turing pattern formation is that there exists a real positive eigenvalue of Jacobi matrix under some wave number k. Then we scan k within [0.01, 100] and sample 1001 points with uniform interval on logarithmic coordinates. For specific k, we solve the following Jacobian equation to get all three λk :

We thank Jing Li and Jingyi Xi for their kind help. This work is partly supported by NSFC grant (11434001), a MSTC grant (2012AA02A702), and a NFFTBS grant (J0630311).

Appendix A. Supporting information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.jtbi.2016.08.005.

References Asai, R., Taguchi, E., Kume, Y., Saito, M., Kondo, S., 1999. Zebrafish leopard gene as a component of the putative reaction-diffusion system. Mech. Dev. 89 (1), 87–92. Bullara, D., De Decker, Y., 2015. Pigment cell movement is not required for generation of Turing patterns in zebrafish skin. Nat. Commun., 6. http://dx.doi.org/ 10.1038/ncomms7971. Chau, A.H., Walter, J.M., Gerardin, J., Tang, C., Lim, W.A., 2012. Designing synthetic regulatory networks capable of self-organizing cell polarization. Cell 151 (2), 320–332. Dillon, R., Maini, P.K., Othmer, H.G., 1994. Pattern formation in generalized Turing systems. J. Math. Biol. 32 (4), 345–393. Driever, W., Nusslein-Volhard, C., 1988. The bicoid protein determines position in the Drosophila embryo in a concentration-dependent manner. Cell 54 (1), 95–104, Epub 1988/07/01. doi: 0092-8674(88)90183-3 [pii]. PubMed PMID:

96

M.M. Zheng et al. / Journal of Theoretical Biology 408 (2016) 88–96

3383245. Gough, Brian, 2009. GNU Scientific Library Reference Manual. Network Theory Ltd,. Gregor, T., Bialek, W., van Steveninck, R.R.R., Tank, D.W., Wieschaus, E.F., 2005. Diffusion and scaling during early embryonic pattern formation. Proc. Natl. Acad. Sci. USA 102 (51), 18403–18407. Gregor, T., Wieschaus, E.F., McGregor, A.P., Bialek, W., Tank, D.W., 2007. Stability and nuclear dynamics of the bicoid morphogen gradient. Cell 130 (1), 141–152, Epub 2007/07/17. doi: S0092-8674(07)00663-0 [pii]10.1016/j.cell.2007.05.026. PubMed PMID: 17632061; PubMed Central PMCID: PMC2253672. Hammonds, A.S., Bristow, C.A., Fisher, W.W., Weiszmann, R., Wu, S., Hartenstein, V., et al., 2013. Spatial expression of transcription factors in Drosophila embryonic organ development. Genome Biol. 14 (12), R140. http://dx.doi.org/10.1186/ gb-2013-14-12-r140 gb-2013-14-12-r140 [pii]. PubMed PMID: 24359758; PubMed Central PMCID: PMC4053779. Holley, S.A., Geisler, R., Nusslein-Volhard, C., 2000. Control of her1 expression during zebrafish somitogenesis by a delta-dependent oscillator and an independent wave-front activity. Genes Dev. 14 (13), 1678–1690, Epub 2000/07/ 11. PubMed PMID: 10887161; PubMed Central PMCID: PMC316735. Howell, A.S., Jin, M., Wu, C.F., Zyla, T.R., Elston, T.C., Lew, D.J., 2012. Negative feedback enhances robustness in the yeast polarity establishment circuit. Cell 149 (2), 322–333. http://dx.doi.org/10.1016/j.cell.2012.03.012, S0092–8674(12) 00342-X [pii]. PubMed PMID: 22500799; PubMed Central PMCID: PMC3680131. Hsia, J., Holtz, W.J., Huang, D.C., Arcak, M., Maharbiz, M.M., 2012. A feedback quenched oscillator produces turing patterning with one diffuser. PLoS Comput. Biol. 8 (1), e1002331. http://dx.doi.org/10.1371/journal.pcbi.1002331, PCOMPBIOL-d-11-00932 [pii]. PubMed PMID: 22291582; PubMed Central PMCID: PMC3266880. Iman, R.L., Davenport, J.M., Zeigler, D.K., 2006. Latin Hypercube Sampling Technical Report SAND79-1473. Sandia Laboratories, Albuquerque. Kondo, S., Miura, T., 2010. Reaction-diffusion model as a framework for understanding biological pattern formation. Science 329 (5999), 1616–1620. Lauschke, V.M., Tsiairis, C.D., Francois, P., Aulehla, A., 2013. Scaling of embryonic patterning based on phase-gradient encoding. Nature 493 (7430), 101–105. Little, S.C., Tkacik, G., Kneeland, T.B., Wieschaus, E.F., Gregor, T., 2011. The formation of the bicoid morphogen gradient requires protein movement from anteriorly localized mRNA. Plos Biol. 9 (3). Long, Y., Qi, O., Wang, H.L., 2012. Dose-response aligned circuits in signaling systems. PLoS ONE 7 (4). Ma, W., Lai, L., Ouyang, Q., Tang, C., 2006. Robustness and modular design of the Drosophila segment polarity network. Mol. Syst. Biol. 2 (1). Ma, W., Trusina, A., El-Samad, H., Lim, W.A., Tang, C., 2009. Defining network topologies that can achieve biochemical adaptation. Cell 138 (4), 760–773. Marcon, L., Diego, X., Sharpe, J., Muller, P., 2016. High-throughput mathematical analysis identifies Turing networks for patterning with equally diffusing

signals. Elife, 5. Meinhardt, H., 2003. The Algorithmic Beauty of Sea Shells, 3rd ed Springer-Verlag, New York. Meinhardt, H., Gierer, A., 1974. Applications of a theory of biological pattern formation based on lateral inhibition. J. Cell Sci. 15 (2), 321–346. Meinhardt, H., Gierer, A., 2000. Pattern formation by local self-activation and lateral inhibition. Bioessays 22 (8), 753–760. Mimura, M., Sakaguchi, H., Matsushita, M., 2000. Reaction-diffusion modelling of bacterial colony patterns. Physica A: Stat. Mech. its Appl. 282 (1–2), 283–303. Miura, T., 2013. Turing and Wolpert work together during limb development. Sci. Signal 6 (270), pe14. http://dx.doi.org/10.1126/scisignal.2004038, scisignal.2004038 [pii]. PubMed PMID: 23572146. Muller, P., Rogers, K.W., Jordan, B.M., Lee, J.S., Robson, D., Ramanathan, S., et al., 2012. Differential diffusivity of nodal and lefty underlies a reaction-diffusion patterning system. Science 336 (6082), 721–724. Murray, J.D., 2003. Mathematical Biology. Springer,. Nakamasu, A., Takahashi, G., Kanbe, A., Kondo, S., 2009. Interactions between zebrafish pigment cells responsible for the generation of Turing patterns. Proc. Natl. Acad. Sci. 106 (21), 8429–8434. Nakao, H., Mikhailov, A.S., 2010. Turing patterns in network-organized activatorinhibitor systems. Nat. Phys. 6 (7), 544–550. Oates, A.C., Morelli, L.G., Ares, S., 2012. Patterning embryos with oscillations: structure, function and dynamics of the vertebrate segmentation clock. Development 139 (4), 625–639. Othmer, H.G., Scriven, L.E., 1971. Instability and dynamic pattern in cellular networks. J. Theor. Biol. 32 (3), 507–537. Othmer, H.G., Scriven, L.E., 1974. Non-linear aspects of dynamic pattern in cellular networks. J. Theor. Biol. 43 (1), 83–112. Quyang, Q., 2010. Introduction to Nonlinear Science and Pattern Formation. Peking University Press, ISBN 13: 9787301159316. Sheth, R., Marcon, L., Bastida, M.F., Junco, M., Quintana, L., Dahn, R., et al., 2012. Hox genes regulate digit patterning by controlling the wavelength of a turing-type mechanism. Science 338 (6113), 1476–1480. Strier, D.E., Ponce Dawson, S., 2007. Turing patterns inside cells. PLoS ONE 2 (10), e1053. http://dx.doi.org/10.1371/journal.pone.0001053, PubMed PMID: 17940616; PubMed Central PMCID: PMC2013944. Turing, A., 1952. The chemical theory of morphogenesis. Philos. Trans. Roy. Soc. 13, 1. Uzkudun, M., Marcon, L., Sharpe, J., 2015. Data-driven modelling of a gene regulatory network for cell fate decisions in the growing limb bud. Mol. Syst. Biol. 11 (7), 815. http://dx.doi.org/10.15252/msb.20145882, PubMed PMID: 26174932. Wu, Y., Zhang, X., Yu, J., Ouyang, Q., 2009. Identification of a topological characteristic responsible for the biological robustness of regulatory networks. PLoS Comput. Biol. 5 (7), e1000442.