- Research article
- Open Access
The capacity for multistability in small gene regulatory networks
- Dan Siegal-Gaskins^{1, 2},
- Erich Grotewold^{1, 2} and
- Gregory D Smith^{1, 3}Email author
https://doi.org/10.1186/1752-0509-3-96
© Siegal-Gaskins et al; licensee BioMed Central Ltd. 2009
- Received: 13 April 2009
- Accepted: 21 September 2009
- Published: 21 September 2009
Abstract
Background
Recent years have seen a dramatic increase in the use of mathematical modeling to gain insight into gene regulatory network behavior across many different organisms. In particular, there has been considerable interest in using mathematical tools to understand how multistable regulatory networks may contribute to developmental processes such as cell fate determination. Indeed, such a network may subserve the formation of unicellular leaf hairs (trichomes) in the model plant Arabidopsis thaliana.
Results
In order to investigate the capacity of small gene regulatory networks to generate multiple equilibria, we present a chemical reaction network (CRN)-based modeling formalism and describe a number of methods for CRN analysis in a parameter-free context. These methods are compared and applied to a full set of one-component subnetworks, as well as a large random sample from 40,680 similarly constructed two-component subnetworks. We find that positive feedback and cooperativity mediated by transcription factor (TF) dimerization is a requirement for one-component subnetwork bistability. For subnetworks with two components, the presence of these processes increases the probability that a randomly sampled subnetwork will exhibit multiple equilibria, although we find several examples of bistable two-component subnetworks that do not involve cooperative TF-promoter binding. In the specific case of epidermal differentiation in Arabidopsis, dimerization of the GL3-GL1 complex and cooperative sequential binding of GL3-GL1 to the CPC promoter are each independently sufficient for bistability.
Conclusion
Computational methods utilizing CRN-specific theorems to rule out bistability in small gene regulatory networks are far superior to techniques generally applicable to deterministic ODE systems. Using these methods to conduct an unbiased survey of parameter-free deterministic models of small networks, and the Arabidopsis epidermal cell differentiation subnetwork in particular, we illustrate how future experimental research may be guided by network structure analysis.
Keywords
- Gene Regulatory Network
- Chemical Reaction Network
- Multiple Equilibrium
- Interaction Graph
- Cell Fate Determination
Background
The availability of high-throughput techniques for gene expression analysis and identification of promoter-transcription factor (TF) interactions has led to characterization of the intricate gene regulatory networks that govern organism behavior [1–3]. These networks are composed of a large number of small and topologically distinct subnetworks, including the overrepresented 'network motifs' [4–7]. In recent years, dynamical systems modeling of regulatory and signaling pathways has provided insight into the equilibrium states and transient dynamics of such subnetworks [8, 9]; for example, detailed cellular and subcellular models demonstrate that interconnected positive and negative feedback loops may give rise to the phenomena of oscillations, excitability, and the existence of multiple stable equilibria (e.g., bistability) [10, 11].
Bistability in particular is ubiquitous in biological systems ranging from biochemical networks to ecosystems [12–16]. In bistable systems, graded inputs (e.g., the concentration of a specific hormone) are converted into a discontinuous ON/OFF response [17–20]. Switch-like behavior is also a characteristic of many developmental processes, and it has been suggested that the maintenance of two distinct phenotypic states in the absence of genetic or environmental differences may sometimes be attributed to bistability in an underlying gene network [21].
Previous attempts at modeling this cell fate determination system have aimed at explaining how trichome patterns form out of a field of initially equivalent epidermal cells, but have ignored the question of how the primary identity decision is made [25, 26]. Such models assume an underlying mechanism of either the 'activator-inhibitor' or 'trapping/depletion' type, both of which include positive regulation of GL3 by the GL3-GL1 active complex. Consistent with the activator-inhibitor model [27], it has been shown experimentally that the activators do positively control the diffusible inhibitors [3]. However, although recent work has suggested that a positive feedback loop may be involved in cell fate specification in the root epidermis [28], no such direct positive feedback has been found in the leaf. Indeed, the regulators of trichome patterning in the leaf may in fact be involved in a negative feedback loop, as recently demonstrated for GL3 [29]. This poses the question of whether the regulatory subnetworks known to be involved in trichome initiation have the capacity for bistability, or whether this is unlikely given the absence of experimental evidence for direct positive autoregulation.
Tools for establishing multistability in ODE models of regulatory networks
has N eigenvalues with negative real part when evaluated at the equilibrium concentrations. Since in practice there is poor knowledge of the rate constants and binding affinities that occur in F (c), we are primarily interested in assessing a network's capacity for multiple equilibria or bistability without specifying these parameters, though all rate constants are assumed to be positive. A simple parameter-free graphical representation of the system is given in its interaction graph in which each species is represented as a graph node. A directed positive edge is drawn from species i to j if ∂F_{ j }/∂c_{ i }(c) > 0, a directed negative edge is drawn if ∂F_{ j }/∂c_{ i }(c) < 0. No edge is drawn if ∂F_{ i }/∂c_{ i }(c) = 0.
where S is a stoichiometric matrix indicating the change in the number of molecules of each species in a particular reaction, and ϕ(c) is a vector of the reaction rates (written as polynomial functions of the species concentrations [37]). While both Eq. 1 and Eq. 3 are nonlinear ODEs, the specific form of CRNs enables the application of powerful CRN-specific network analysis techniques. One computational method can rule out multiple equilibria in a CRN by determining if the associated polynomial function is injective, i.e., it maps distinct arguments to distinct values [37, 38]. A related graphical technique involves analysis of the CRN species-reaction (SR) graph [39, 40] (see Additional file 1 for a description of the SR graph and its associated multistationarity theorem). There is also a suite of computational methods utilizing the so-called deficiency theorems for CRNs that are implemented in the publicly-available Chemical Reaction Network Toolbox [41]. These deficiency theorems give conditions for the existence, multiplicity, and stability of CRN steady-states that may often be applied in a parameter-free context; for example, the deficiency zero theorem states that if a CRN has certain topological properties, then within each compatibility class (an invariant manifold in species concentration space in which solutions are bound and determined by initial conditions), there is exactly one steady-state with strictly positive concentrations, and this steady-state is locally asympotically stable [42]. Interested readers may refer to refs. [42–48] for details on the full range of theorems implemented in the CRNT.
This section has briefly introduced five tools that can be used to determine if a given regulatory network has the capacity for multiple equilibria: the Thomas conjecture and Kauffman's multistationarity conditions based on analysis of the network's interaction graph (denoted below as the IG-T and IG-K methods, respectively), a multistationarity theorem based on the structure of the CRN SR graph (denoted by SRG), a computational method that can establish an injective property for the CRN's polynomial function (denoted by INJ), and CRN theory as implemented by the Chemical Reaction Network Toolbox (denoted by CRNT). Below we use these tools to analyze several small networks consisting of one gene and one gene product (the one-component subnetworks), in order to establish both the effectiveness of the various methods and the subnetworks' capacity for multistability. We then investigate a large set of more complex two-component subnetworks, including the well-studied 'double negative' feedback system. These small networks are studied without regard for their frequency of occurrence in the larger regulatory machinery of real biological systems, that is, we do not restrict ourselves to overrepresented network motifs. Lastly, we revisit the epidermal differentiation system in Arabidopsis as a real-world example of the applicability of these techniques.
Results and Discussion
Bistability in two positive autoregulatory subnetworks: a comparison of methods
In Fig. 3A, the motif is specified as a CRN that includes basal production of P (X → X + P), degradation of P (P → ∅), reversible binding of P to the promoter of X (X + P ⇌ XP), and production of P by the TF-gene complex XP (XP → XP + P). In Fig. 3B there is also basal production of P from X and degradation of P, however two TFs now associate to form a homodimer (P + P ⇌ PP), and it is this PP dimer that binds to the promoter of X (X + PP ⇌ XPP) and produces P via the dimer-gene complex XPP (XPP → XPP + P). From these sets of elementary reactions, which following ref. [49] can also be visualized as wiring diagrams (see second row of Fig. 3), one can write an unambiguous system of ODEs that take the form of a CRN (Eq. 3). Although it has been shown that transcription and elements of post-transcriptional control combine to regulate the protein production rate, and that stochasticity exists at all levels of regulation [50], for simplicity we model protein production as a single deterministic process (cf. [51]). We address the role of translation and mRNA degradation in determining the capacity for multiple stable states in the next section.
For both Fig. 3A (the monomer model) and Fig. 3B (the dimer model), the existence of positive circuits in the interaction graph (e.g., X-XP-X and P-XP-P in Fig. 3A and X-XPP-X and P-PP-P in Fig. 3B) means that the IG-T theorem does not preclude multiple equilibria. Neither does the IG-K theorem preclude multiple equilibria; in Fig. 3A, there is a variable nucleus that includes the variable circuit X-P-X, and in Fig. 3B, there are two nuclei of opposite signs (the positive nuclei composed of circuits P-PP-P, X-X, and XPP-XPP, and the negative nuclei composed of circuits P-PP-P and X-XPP-X). (Eliminating one equation using the conserved quantities [X] + [XP] (Fig. 3A) and [X] + [XPP] (Fig. 3B) does not change the results of the IG-T and IG-K analysis; see Additional file 1.) The SRG method is also unable to rule out bistability in these motifs, in Fig. 3A due to the splitting of the X + P complex pair (or c-pair) associated with reaction X → X + P by two even-cycles (one composed of X and X → X + P, and the second composed of X, X → X + P, P, and X + P → XP), and in Fig. 3B because the cycle containing P and P + P ⇌ PP is not a 1-cycle or an odd-cycle. In addition, an injective property cannot be established for either model's polynomial function.
The bistability of the dimer version of the autoregulatory motif, confirmed using CRNT, is consistent with the inability of the IG-T, IG-K, SRG, and INJ methods to rule it out. However, these methods are also unable to rule out bistability in even the simple case of positive autoregulation by a protein monomer, suggesting that they may be of limited use in the analysis of more complex systems.
Analysis of twelve one-component subnetworks
Construction of one-component regulatory subnetworks
Reaction label | Reaction |
---|---|
a | X → X + P |
b | P → ∅ |
c | X + P ⇌ XP |
d | XP → XP + P |
e | P + P ⇌ PP |
f | X + PP ⇌ XPP |
g | XPP → XPP + P |
Survey of one-component regulatory subnetworks
Subnetwork | Dimer formation | Monomer binding | Dimer binding | Multiple equilibria ruled out? | ||||
---|---|---|---|---|---|---|---|---|
IG-T | IG-K | SRG | INJ | CRNT | ||||
ab | no | no | no | yes | yes | yes | yes | yes |
abcd | no | yes ^{+} | no | no | no | no | no | yes |
abc | no | yes ^{-} | no | no | no | no | yes | yes |
abe | yes | no | no | no | no | yes | yes | yes |
abefg | yes | no | yes ^{+} | no | no | no | no | no |
abef | yes | no | yes ^{-} | no | no | no | yes | yes |
abcde | yes | yes ^{+} | no | no | no | no | no | yes |
abcdefg | yes | yes ^{+} | yes ^{+} | no | no | no | no | no |
abcdef | yes | yes ^{+} | yes ^{-} | no | no | no | no | yes |
abce | yes | yes ^{-} | no | no | no | no | yes | yes |
abcef | yes | yes ^{-} | yes ^{-} | no | no | no | yes | yes |
abcefg | yes | yes ^{-} | yes ^{+} | no | no | no | no | no |
In contrast to these graphical methods, the INJ ruled out multistability in six of twelve cases, and the CRNT ruled out multistability in nine of twelve cases. Furthermore, for the three subnetworks where multiple equilibria were not ruled out by the CRNT (abefg, abcdefg, and abcefg), parameter sets leading to bistability were provided. In these two methods we begin to see the strength of the full chemical reaction network theory: with the assumption of mass-action kinetics, the underlying structure of CRNs limits the nonlinearities that appear in the full ODE model [48]. The bifurcation diagram for subnetwork abefg is shown in Fig. 4. Similar results for subnetworks abcdefg and abcefg are shown in Additional file 1, Fig. S5.
As previously mentioned, the twelve subnetworks presented in Table 2 were constructed under the assumption that, for the purpose of analyzing subnetwork equilibria, transcription and translation may be combined into a single operation. We repeated our analysis using one-component subnetworks augmented by elementary reactions that make these processes explicit. The elementary reactions X → X + P, XP → XP + P, XPP → XPP + P were replaced by reactions representing transcription (X → X + R, XP → XP + R, XPP → XPP + R), translation (R → R + P), and mRNA degradation (R → ∅). Application of the CRNT to these augmented one-component systems led to no change in the results of Table 2, and thus we model protein production as a single one-step process for the remainder of this paper.
Survey of two-component regulatory subnetworks
Consideration of small networks consisting of two genes (X_{1} and X_{2}) and two gene products (P_{1} and P_{2}) is a natural extension of the analysis presented above. With a second component included, the number of relevant subnetworks increases to 40,680 (Additional file 1, Table S1 lists the constituent reactions). The subnetworks may be grouped by their number of dimer/promoter association reactions; every two-component subnetwork has between zero and six such reactions (two homodimers and one heterodimer which can each bind to either promoter). For bistability analysis, we applied the INJ method in an automated fashion using Matlab scripts provided by G. Craciun [37]. No attempt was made to apply the IG-T, IG-K, and SRG methods that were found to be uninformative with respect to subnetworks with one-component (Table 2). Approximately 15,000 subnetworks were randomly selected for study, sampled in a manner that ensured that subnetworks with zero to six dimer-binding reactions were represented in the same proportion as in the complete set of subnetworks.
Bootstrap analysis results for two-component regulatory subnetworks
# dimers binding | % of total models | % with multiple equilibria ruled out |
---|---|---|
0 | 0.8 | 26.2 [19.0, 34.4] |
1 | 4.8 | 19.8 [17.1, 22.9] |
2 | 14.5 | 12.6 [11.3, 14.1] |
3 | 25.5 | 8.9 [8.1, 9.9] |
4 | 28.8 | 4.1 [3.5, 4.7] |
5 | 19.1 | 1.4 [1.0, 1.9] |
6 | 6.5 | 0.0 |
In the preceding analysis, we did not distinguish between productive and unproductive TF binding reactions. We define a productive reaction as one in which expression of the corresponding gene is activated. Similarly, an unproductive reaction is one in which gene expression is repressed. When the analysis is applied to subsets of systems lacking either productive or unproductive TF binding, the percentage of sampled subnetworks for which INJ could note rule out bistability increases as the number of binding reactions increases (see Additional file 1, Tables S3 and S4). Note that the INJ method is more successful at ruling out bistability for subnetworks lacking productive TF binding than it is for subnetworks lacking unproductive TF binding, perhaps because bistable networks are enriched in productive TF-promoter interactions.
The CRNT was found to be the most effective method for establishing and ruling out bistability in the one-component subnetworks. However, the current release of the CRNT cannot be automated, and is thus difficult to apply to a large set of two-component systems. For this reason we randomly selected 25 two-component subnetworks and applied the INJ and CRNT methods to compare their effectiveness. The CRNT found 19 subnetworks that do not support multiple equilibria (76% with confidence interval [52%, 88%]), and it generated parameter sets that lead to bistability for 5 of the remaining 6 subnetworks (20% of total with confidence interval [4%, 36%]). This is quite similar to the fraction of one-component subnetworks that were found by the CRNT to be bistable (25%) and well within the confidence interval. Conversely, the INJ method was only able to rule out bistability for 3 of the 19 subnetworks ruled out by the CRNT (Additional file 1, Table S2). Although derived from a small sample, taken together these results suggest that the INJ method is debilitated by network size to a larger degree than the CRNT method.
The canonical reciprocal repression subnetwork
Revisiting the trichome differentiation subnetwork with the CRNT
As described above and shown in Fig. 1, the trichome differentiation subnetwork consists of regulators belonging to the R2R3-MYB (e.g., GL1) and bHLH (e.g., GL3) classes that directly regulate GL2 and CPC among several other targets [3]. As there is no experimental evidence for the direct positive feedback that has been posited in current models of trichome initiation, the question of whether the system can support bistability without it is one of great interest.
CRNs consistent with the core trichome differentiation network
Reactions | Model 1 | Model 2 | Model 3 | Model 4 |
---|---|---|---|---|
X → X + CPC | ✔ | ✔ | ✔ | ✔ |
CPC → ∅ | ✔ | ✔ | ✔ | ✔ |
GL3 + GL1 ⇌ GL3-GL1 | ✔ | ✔ | ✔ | ✔ |
GL3 + CPC ⇌ GL3-CPC | ✔ | ✔ | ✔ | ✔ |
GL3-GL1 + GL3-GL1 ⇌ GL3-GL1-GL3-GL1 | ✔ | ✔ | ||
GL3-GL1 + X ⇌ GL3-GL1-X | ✔ | ✔ | ✔ | |
GL3-GL1 + GL3-GL1-X ⇌ GL3-GL1-GL3-GL1-X | ✔ | |||
GL3-GL1-X → GL3-GL1-X + CPC | ✔ | ✔ | ||
GL3-GL1-GL3-GL1 + X ⇌ GL3-GL1-GL3-GL1-X | ✔ | ✔ | ||
GL3-GL1-GL3-GL1-X → GL3-GL1-GL3-GL1-X + CPC | ✔ | ✔ | ✔ | |
Result of CRNT analysis: | not bistable | bistable | bistable | bistable |
Conclusion
We have used a modeling framework in which small gene regulatory networks are specified as chemical reaction networks (CRNs), with reactions describing transcription factor (TF) production, degradation, dimerization, and association of both monomer and dimer TFs with promoters in a productive or unproductive fashion. We used these methods to survey twelve subnetworks consisting of one TF gene and gene product (the one-component subnetworks) and found that bistability is only exhibited by those that include positive feedback with cooperative TF binding (Table 2). This is not unexpected, as direct positive feedback and non-linearity are often prerequisites for bistability in simple systems [15, 20]. When the analysis tools were applied to a large random sample of two-component subnetworks (from a total of 40,680 consistent with our CRN specification rules; see Additional file 1), we found that the presence of cooperative TF binding makes it more difficult to rule out bistability (Table 3), although it is ruled out more often in subnetworks lacking direct positive feedback than it is for those lacking direct negative feedback (Additional file 1, Tables S3 and S4). Interestingly, we found several examples of bistable two-component subnetworks that do not include promotor-TF dimer interactions at all (Fig. 5 and Additional file 1, Table S2), consistent with prior work demonstrating that sequestration and titration of active TFs into inactive complexes can give rise to non-linearity and bistable behavior [51, 56].
Our application of the various equilibria analysis tools suggests that the interaction graph-based techniques (IG-T and IG-K) that are generally applicable to deterministic ODE systems are rarely informative with respect to small gene regulatory networks. We also found that the deficiency theorems implemented in the Chemical Reaction Network Toolbox (CRNT) are generally more informative than computational (INJ) and graphical (SRG) methods that attempt to rule out multiple equilibria by establishing injectivity of the reaction network polynomial function. We emphasize that the success of the CRNT is partly due to the special form of chemical reaction networks that assume kinetics of the mass-action form. However, this assumption would appear to be a valid one under many working conditions. That the CRNT accurately determined the stability properties of the well-studied positive autoregulatory and reciprocal repression systems (Fig. 3 and Fig. 6, respectively), and that one-component systems found by CRNT to be bistable were confirmed as such with bifurcation diagrams (Fig. 4 and Additional file 1, Fig. S5), is further validation of the CRNT's usefulness.
It is important to note that the analytical methods described herein are not the only ones available for probing the stability properties of a dynamical system. One additional method involves evaluation of the Jacobian matrix (Eq. 2), because the existence of a saddle-node bifurcation leading to stable and unstable equilibria requires the Jacobian to be singular at the critical point c_{ ss }(i.e., det{J(c_{ ss })} = 0) [57]. We do not use continuation methods to directly rule out saddle-node bifurcations in this work, as our attempts at such analysis have not yielded results, and to our knowledge it would be difficult to explore large numbers of regulatory networks in this fashion. Furthermore, deterministic ODEs are only one of several mathematical structures used to model gene regulatory networks; other possibilities include boolean networks, Petri nets, Markov chains, stochastic ODEs, and so on. Still, the parameter-free approach to modeling that we have emphasized here is an attractive one for at least two reasons: (1) there are often multiple assumptions made (such as the validity of Michaelis-Menten or Hill-type expressions) when complex network models are constructed as a generic system of ODEs, and these assumptions may result in inconsistencies [49], and (2) as we have shown, CRN theory facilitates analysis that is usually more fruitful than techniques available for generic ODEs.
The current state of CRN-based modeling and analysis does not include stochastic effects, despite the fact that there is both experimental and theoretical evidence that low concentrations of TFs or binding sites can lead to fluctuations (i.e, noise) that may sometimes be an essential aspect of network function [50, 58–60]. The CRN approach thus represents a first step towards generating mixed models of biological systems that incorporate stochastic effects. Future work could extend CRN-based techniques to include stochastic dynamics consistent with the elementary processes that compose gene regulatory networks [61]. We further anticipate that new releases of the CRNT software package will also be amenable to automation and will allow for CRNs of unlimited size.
Although the modeling formalism used here is a simple one, it is not clear that an increase in complexity would lead to different conclusions; for example, we have already demonstrated that excluding mRNA does not influence the one-component subnetwork results. Two other simplifications made find significant experimental support: (i) that TF dimers are stable to proteolytic degradation (cf. [51, 56]), a reasonable assumption given that protein dimerization typically provides protection against proteolysis by enhancing thermal stability or through blocking access to monomer degradation tags [52], and (ii) the exclusion of higher-order oligomeric TF complexes (such as trimers or tetramers), an assumption appropriate for many eukaryotic TF gene families where dimerization is sufficient for functionality [62]. On the other hand, TF activation can occur by posttranslational modification (e.g., phosphorylation [51]), a situation analogous to cell signaling pathways in which dual phosphorylation/dephosphorylation cycles confer bistability [63]. While it is possible to extend the modeling formalism illustrated here in these and other directions, the additional combinatorial complexity would make surveys of multi-component networks computationally challenging.
Applying the CRNT to the real-world example of Arabidopsis epidermal cell differentiation, we have shown that gene regulatory subnetworks consistent with established components and interactions may have the capacity for bistability. This capacity of course depends on the specific details of the network architecture. In particular, we found that dimerization of the active GL3-GL1 complex and cooperative sequential binding of GL3-GL1 to the CPC promoter are each independently sufficient to generate a bistable subnetwork, even in the absence of direct positive autoregulation. Both of these alternatives have some experimental support. Sequential binding of GL3-GL1 to the CPC promoter is plausible, as multiple TF binding sites in gene promoters are a general characteristic of eukaryotic gene expression [64]. And while dimerization of GL3-GL1 has not been experimentally demonstrated, GL3 can dimerize [65] either through the conserved bHLH domain, as is the case for other bHLH factors [66], or through the conserved C-terminal ACT domain [67]. Bistability may also be achieved if the GL3-GL1 complex is not stable against proteolysis (not shown).
If it can be determined that GL3-GL1 does not in fact dimerize or bind to two separate binding sites in the CPC promoter, and that it is in fact a stable complex, then the current view of the core trichome differentiation network (Fig. 1) must be incomplete. In that case, there are several intriguing possibilities. It may be that one or more of the many additional trichome-related genes that have been recently identified [3] are as central to trichome differentiation as GL3, GL1, and CPC. It is possible that epidermal cell differentiation in Arabidopsis is not in fact subserved by a bistable gene regulatory network. It may also be that some unknown extrinsic mechanism plays a role in Arabidopsis epidermal cell fate determination. A mechanism of this kind may render a bistable network irreversible by causing one equilibrium state to be inaccessible [21], or it may fix a transient activated state before the system can return to a basal equilibrium [68]. There is some evidence for these two latter possibilities, as a differentiated trichome cell is unlikely to revert back to an undifferentiated state. In particular, a recent study has shown that the induction of GL3 function for 4 hrs is sufficient to trigger trichome initiation, with development continuing even after GL3 function is completely removed [29]. Clearly there are a number of possibilities that we have only begun to investigate. Future experiments will be required to distinguish between them and may yet provide insights into this interesting developmental model system.
Methods
Analyzing bistability with the CRNT
The Chemical Reaction Network Toolbox (CRNT) is a stand-alone DOS program that uses CRN theory [42, 44–47] to to determine if a given CRN has the capacity for multiple equilibria; if so, it will often provide example parameter values (see ref. [41]). In Additional file 1 we detail rules designed to specify families of a gene regulatory networks of interest (e.g., the one-component subnetworks shown in Tables 1 and 2, each of which is a CRN). To analyze a subnetwork using CRNT, we list the species, complexes (including the null complex ∅), and reactions, and enter them manually into the CRNT Network Analyst. With the currently available implementation of the CRNT we were restricted to the analysis of subnetworks containing twenty complexes or less. For most subnetworks discussed here, the CRNT Advanced Deficiency Theory option is needed following the initial analysis.
Generation of one- and two-component regulatory subnetworks
One- and two-component subnetworks were generated in the Matlab programming environment (The MathWorks, Inc.). One-component networks are distinguished by whether TF monomer can bind the promoter and, if monomer binding does occur, whether the monomer-promotor complex results in transcription; thus there are three options with respect to monomer binding (no, yes^{-}, yes^{+}; see Table 2). One-component networks are further distinguished by whether or not dimers can form (no, yes) and, if so, whether the dimer-promoter complex is productive (no, yes^{+}, yes^{-}). Because this latter question is only relevant when dimer can form, there are 3 + 3.3 = 12 distinct one-component networks.
The generation of two-component networks involves enumerating all possibilities with regard to homodimer and heterodimer formation, being careful to avoid double counting of networks that are equivalent via symmetry, i.e., when all species (X_{1}, X_{2}, P_{1}, P_{2}, etc.) are relabeled by exchanging the subscripts 1 and 2. The heterodimers P_{1}P_{2} and P_{2}P_{1} were assumed to be equivalent. These distinct possibilities for homodimer and heterodimer formation were subsequently instantiated multiple times to account for each possible pattern of monomer and dimer binding, resulting in 40,680 two-component networks.
Declarations
Acknowledgements
This material is based on work supported by the National Science Foundation under Grants No. DMS-0443843 to GDS and MCB-0418891 and DBI-0701405 to EG, and under Agreement No. 0635561. The author's gratefully acknowledge research leaves during academic year 2007-2008 supported by long-term visitor positions at the Mathematical Biosciences Institute at The Ohio State University (GDS and EG) and by the College of William and Mary (GDS). The authors would also like to thank Joe Pomerening and the anonymous reviewers for their careful reading of the manuscript and their extremely helpful comments, as well as Gheorghe Craciun for instructive conversations and for the Matlab code that tests for injectivity of a CRN polynomial function [37].
Authors’ Affiliations
References
- Lee TI, Rinaldi NJ, Robert F, Odom DT, Bar-Joseph Z, Gerber GK, Hannett NM, Harbison CT, Thompson CM, Simon I, Zeitlinger J, Jennings EG, Murray HL, Gordon DB, Ren B, Wyrick JJ, Tagne JB, Volkert TL, Fraenkel E, Gifford DK, Young RA: Transcriptional regulatory networks in Saccharomyces cerevisiae. Science. 2002, 298 (5594): 799-804. 10.1126/science.1075090View ArticlePubMedGoogle Scholar
- McGrath PT, Lee H, Zhang L, Iniesta AA, Hottes AK, Tan MH, Hillson NJ, Hu P, Shapiro L, McAdams HH: High-throughput identification of transcription start sites, conserved promoter motifs and predicted regulons. Nat Biotech. 2007, 25 (5): 584-592. 10.1038/nbt1294.View ArticleGoogle Scholar
- Morohashi K, Grotewold E: A systems approach reveals regulatory circuitry for Arabidopsis trichome initiation by the GL3 and GL1 selectors. PLoS Genet. 2009, 5 (2): e1000396- 10.1371/journal.pgen.1000396PubMed CentralView ArticlePubMedGoogle Scholar
- Zhang L, King O, Wong S, Goldberg D, Tong A, Lesage G, Andrews B, Bussey H, Boone C, Roth F: Motifs, themes and thematic maps of an integrated Saccharomyces cerevisiae interaction network. J Biol. 2005, 4 (2): 6- 10.1186/jbiol23PubMed CentralView ArticlePubMedGoogle Scholar
- Martinez NJ, Ow MC, Barrasa MI, Hammell M, Sequerra R, Doucette-Stamm L, Roth FP, Ambros VR, Walhout AJM: A C. elegans genome-scale microRNA network contains composite feedback motifs with high flux capacity. Genes Dev. 2008, 22 (18): 2535-2549. 10.1101/gad.1678608PubMed CentralView ArticlePubMedGoogle Scholar
- Milo R, Shen-Orr S, Itzkovitz S, Kashtan N, Chklovskii D, Alon U: Network motifs: Simple building blocks of complex networks. Science. 2002, 298 (5594): 824-827. 10.1126/science.298.5594.824View ArticlePubMedGoogle Scholar
- Shen-Orr S, Milo R, Mangan S, Alon U: Network motifs in the transcriptional regulation network of Escherichia coli. Nat Genet. 2002, 31 (1): 64-68. 10.1038/ng881View ArticlePubMedGoogle Scholar
- Tyson JJ: Bringing cartoons to life. Nature. 2007, 445 (7130): 823-823. 10.1038/445823aView ArticlePubMedGoogle Scholar
- Bennett MR, Volfson D, Tsimring L, Hasty J: Transient dynamics of genetic regulatory networks. Biophysical Journal. 2007, 92 (10): 3501-12. 10.1529/biophysj.106.095638PubMed CentralView ArticlePubMedGoogle Scholar
- Tyson JJ, Chen KC, Novak B: Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell. Curr Opin Cell Biol. 2003, 15 (2): 221-231. 10.1016/S0955-0674(03)00017-6View ArticlePubMedGoogle Scholar
- Rand DA, Shulgin BV, Salazar D, Millar AJ: Design principles underlying circadian clocks. J R Soc Interface. 2004, 1: 119-130. 10.1098/rsif.2004.0014PubMed CentralView ArticlePubMedGoogle Scholar
- Laurent M, Johannin G: Molecular clues to pathogenesis in prion diseases. Histol Histopathol. 1997, 12 (2): 583-594.PubMedGoogle Scholar
- Dubnau D, Losick R: Bistability in bacteria. Mol Microbiol. 2006, 61 (3): 564-572. 10.1111/j.1365-2958.2006.05249.xView ArticlePubMedGoogle Scholar
- Avery S: Cell individuality: the bistability of competence development. Trends Microbiol. 2005, 13 (10): 459-462. 10.1016/j.tim.2005.08.006View ArticlePubMedGoogle Scholar
- Ferrell J: Self-perpetuating states in signal transduction: positive feedback, double-negative feedback and bistability. Curr Opin Cell Biol. 2002, 14 (2): 140-148. 10.1016/S0955-0674(02)00314-9View ArticlePubMedGoogle Scholar
- Rietkerk M, Dekker S, de Ruiter P, Koppel van de J: Self-organized patchiness and catastrophic shifts in ecosystems. Science. 2004, 305 (5692): 1926-1929. 10.1126/science.1101867View ArticlePubMedGoogle Scholar
- Paliwal S, Iglesias PA, Campbell K, Hilioti Z, Groisman A, Levchenko A: MAPK-mediated bimodal gene expression and adaptive gradient sensing in yeast. Nature. 2007, 446 (7131): 46-51. 10.1038/nature05561View ArticlePubMedGoogle Scholar
- Laslo P, Spooner CJ, Warmflash A, Lancki DW, Lee HJ, Sciammas R, Gantner BN, Dinner AR, Singh H: Multilineage transcriptional priming and determination of alternate hematopoietic cell fates. Cell. 2006, 126 (4): 755-766. 10.1016/j.cell.2006.06.052View ArticlePubMedGoogle Scholar
- Ozbudak EM, Thattai M, Lim HN, Shraiman BI, van Oudenaarden A: Multistability in the lactose utilization network of Escherichia coli. Nature. 2004, 427 (6976): 737-740. 10.1038/nature02298View ArticlePubMedGoogle Scholar
- Becskei A, Seraphin B, Serrano L: Positive feedback in eukaryotic gene networks: cell differentiation by graded to binary response conversion. EMBO J. 2001, 20 (10): 2528-2535. 10.1093/emboj/20.10.2528PubMed CentralView ArticlePubMedGoogle Scholar
- Laurent M, Kellershohn N: Multistability: a major means of differentiation and evolution in biological systems. Trends Biochem Sci. 1999, 24 (11): 418-422. 10.1016/S0968-0004(99)01473-5View ArticlePubMedGoogle Scholar
- Schellmann S, Hulskamp M: Epidermal differentiation: trichomes in Arabidopsis as a model system. Int J Dev Biol. 2005, 49 (5-6): 579-584. 10.1387/ijdb.051983ssView ArticlePubMedGoogle Scholar
- Zhao M, Morohashi K, Hatlestad G, Grotewold E, Lloyd A: The TTG1-bHLH-MYB complex controls trichome cell fate and patterning through direct targeting of regulatory loci. Development. 2008, 135 (11): 1991-1999. 10.1242/dev.016873View ArticlePubMedGoogle Scholar
- Ishida T, Kurata T, Okada K, Wada T: A genetic regulatory network in the development of trichomes and root hairs. Annu Rev Plant Biol. 2008, 59: 365-386. 10.1146/annurev.arplant.59.032607.092949View ArticlePubMedGoogle Scholar
- Digiuni S, Schellmann S, Geier F, Greese B, Pesch M, Wester K, Dartan B, Mach V, Srinivas BP, Timmer J, Fleck C, Hulskamp M: A competitive complex formation mechanism underlies trichome patterning on Arabidopsis leaves. Mol Syst Biol. 2008, 4: 217- 10.1038/msb.2008.54PubMed CentralView ArticlePubMedGoogle Scholar
- Bouyer D, Geier F, Kragler F, Schnittger A, Pesch M, Wester K, Balkunde R, Timmer J, Fleck C, Hülskamp M: Two-dimensional patterning by a trapping/depletion mechanism: the role of TTG1 and GL3 in Arabidopsis trichome formation. PLoS Biol. 2008, 6 (6): e141- 10.1371/journal.pbio.0060141PubMed CentralView ArticlePubMedGoogle Scholar
- Meinhardt H, Gierer A: Applications of a theory of biological pattern formation based on lateral inhibition. J Cell Sci. 1974, 15 (2): 321-346.PubMedGoogle Scholar
- Kang YH, Kirik V, Hulskamp M, Nam KH, Hagely K, Lee MM, Schiefelbein J: The MYB23 Gene Provides a Positive Feedback Loop for Cell Fate Specification in the Arabidopsis Root Epidermis. Plant Cell. 2009, 21 (4): 1080-94. 10.1105/tpc.108.063180PubMed CentralView ArticlePubMedGoogle Scholar
- Morohashi K, Zhao MZ, Yang ML, Read B, Lloyd A, Lamb R, Grotewold E: Participation of the Arabidopsis bHLH factor GL3 in trichome initiation regulatory events. Plant Physiol. 2007, 145 (3): 736-746. 10.1104/pp.107.104521PubMed CentralView ArticlePubMedGoogle Scholar
- Ideker T, Lauffenburger D: Building with a scaffold: emerging strategies for high- to low-level cellular modeling. Trends Biotechnol. 2003, 21 (6): 255-262. 10.1016/S0167-7799(03)00115-XView ArticlePubMedGoogle Scholar
- Smolen P, Baxter DA, Byrne JH: Mathematical modeling of gene networks. Neuron. 2000, 26 (3): 567-580. 10.1016/S0896-6273(00)81194-0View ArticlePubMedGoogle Scholar
- Thomas R: On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations. Numerical Methods in the Study of Critical Phenomena. Edited by: Della Dora J, Demongeot J, Lacolle B. 1981, 180-193. Berlin: Springer-VerlagView ArticleGoogle Scholar
- Soulé C: Graphic Requirements for Multistationarity. ComPlexUs. 2003, 1 (3): 123-133. 10.1159/000076100.View ArticleGoogle Scholar
- Kaufman M, Soule C, Thomas R: A new necessary condition on interaction graphs for multistationarity. J Theor Biol. 2007, 248 (4): 675-685. 10.1016/j.jtbi.2007.06.016View ArticlePubMedGoogle Scholar
- Eisenfeld J, DeLisi C: On conditions for qualitative instability of regulatory circuits with application to immunological control loops. Mathematics and Computers in Biomedical Applications. Edited by: Eisenfeld J, DeLisi C. 1985, 39-53. Amsterdam: ElsevierGoogle Scholar
- Feinberg M: Chemical oscillations, multiple equilibria and reaction network structure. Dynamics and Modelling of Reactive Systems. Edited by: Conley C. 1980, 59-130. New York: Academic PressGoogle Scholar
- Craciun G, Feinberg M: Multiple equilibria in complex chemical reaction networks: I. The injectivity property. SIAM J Appl Math. 2005, 65 (5): 1526-1546. 10.1137/S0036139904440278.View ArticleGoogle Scholar
- Craciun G, Feinberg M: Multiple equilibria in complex chemical reaction networks: II. The species-reaction graph. SIAM J Appl Math. 2006, 66 (4): 1321-1338. 10.1137/050634177.View ArticleGoogle Scholar
- Craciun G, Tang Y, Feinberg M: Understanding bistability in complex enzyme-driven reaction networks. Proc Natl Acad Sci USA. 2006, 103 (23): 8697-8702. 10.1073/pnas.0602767103PubMed CentralView ArticlePubMedGoogle Scholar
- Schlosser PM, Feinberg M: A theory of multiple steady-states in isothermal homogeneous CFSTRs with many reactions. Chem Eng Sci. 1994, 49 (11): 1749-1767. 10.1016/0009-2509(94)80061-8.View ArticleGoogle Scholar
- Feinberg M, Ellison P: The Chemical Reaction Network Toolbox, Version 1.1. 2000, http://www.che.eng.ohio-state.edu/~feinberg/crnt/Google Scholar
- Feinberg M: Chemical-reaction network structure and the stability of complex isothermal reactors: I. The deficiency-zero and deficiency-one theorems. Chem Eng Sci. 1987, 42 (10): 2229-2268. 10.1016/0009-2509(87)80099-4.View ArticleGoogle Scholar
- Feinberg M: Lectures on chemical reaction networks. 1979, [Delivered at the Mathematics Research Center, University of Wisconsin-Madison], http://www.che.eng.ohio-state.edu/~FEINBERG/LecturesOnReactionNetworks/Google Scholar
- Feinberg M: Chemical reaction network structure and the stability of complex isothermal reactors: II. Multiple steady states for networks of deficiency one. Chem Eng Sci. 1988, 43 (1): 1-25. 10.1016/0009-2509(88)87122-7.View ArticleGoogle Scholar
- Feinberg M: The existence and uniqueness of steady states for a class of chemical reaction networks. Arch Ration Mech An. 1995, 132 (4): 311-370. 10.1007/BF00375614.View ArticleGoogle Scholar
- Feinberg M: Multiple steady states for chemical reaction networks of deficiency one. Arch Ration Mech An. 1995, 132 (4): 371-406. 10.1007/BF00375615.View ArticleGoogle Scholar
- Ellison P, Feinberg M: How catalytic mechanisms reveal themselves in multiple steady-state data: I. Basic principles. J Mol Catal A-Chem. 2000, 154 (1-2): 155-167. 10.1016/S1381-1169(99)00371-4.View ArticleGoogle Scholar
- Gunawardena J: Chemical reaction network theory for in-silico biologists. 2003, [Lecture Notes.], http://www.jeremy-gunawardena.com/papers/crnt.pdfGoogle Scholar
- Sabouri-Ghomi M, Ciliberto A, Kar S, Novak B, Tyson JJ: Antagonism and bistability in protein interaction networks. J Theor Biol. 2008, 250: 209-218. 10.1016/j.jtbi.2007.09.001View ArticlePubMedGoogle Scholar
- Mata J, Marguerat S, Bahler A: Post-transcriptional control of gene expression: a genome-wide perspective. Trends Biochem Sci. 2005, 30 (9): 506-514. 10.1016/j.tibs.2005.07.005View ArticlePubMedGoogle Scholar
- Francois P, Hakim V: Design of genetic networks with specified functions by evolution in silico. Proc Natl Acad Sci USA. 2004, 101 (2): 580-585. 10.1073/pnas.0304532101PubMed CentralView ArticlePubMedGoogle Scholar
- Buchler N, Gerland U, Hwa T: Nonlinear protein degradation and the function of genetic circuits. Proc Natl Acad Sci USA. 2005, 102 (27): 9559-9564. 10.1073/pnas.0409553102PubMed CentralView ArticlePubMedGoogle Scholar
- Monod J, Jacob F: Teleonomic mechanisms in cellular metabolism, growth, and differentiation. Cold Spring Harb Symp Quant Biol. 1961, 26: 389-401.View ArticlePubMedGoogle Scholar
- Cherry J, Adler F: How to make a Biological Switch. J Theor Biol. 2000, 203 (2): 117-133. 10.1006/jtbi.2000.1068View ArticlePubMedGoogle Scholar
- Gardner T, Cantor C, Collins J: Construction of a genetic toggle switch in Escherichia coli. Nature. 2000, 403 (6767): 339-342. 10.1038/35002131View ArticlePubMedGoogle Scholar
- Buchler NE, Louis M: Molecular Titration and Ultrasensitivity in Regulatory Networks. J Mol Biol. 2008, 384 (5): 1106-1119. 10.1016/j.jmb.2008.09.079View ArticlePubMedGoogle Scholar
- Seydel R: Practical Bifurcation and Stability Analysis: From Equilibrium to Chaos. 1994, New York: Springer-Verlag, 2Google Scholar
- Kaern M, Elston TC, Blake WJ, Collins JJ: Stochasticity in gene expression: from theories to phenotypes. Nat Rev Genet. 2005, 6 (6): 451-64. 10.1038/nrg1615View ArticlePubMedGoogle Scholar
- Kepler TB, Elston TC: Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations. Biophys J. 2001, 81 (6): 3116-36. 10.1016/S0006-3495(01)75949-8PubMed CentralView ArticlePubMedGoogle Scholar
- Adalsteinsson D, McMillen D, Elston TC: Biochemical Network Stochastic Simulator (BioNetS): software for stochastic modeling of biochemical networks. BMC Bioinformatics. 2004, 5: 24- 10.1186/1471-2105-5-24PubMed CentralView ArticlePubMedGoogle Scholar
- Keizer J: Statistical Thermodynamics of Nonequilibrium Processes. 1987, Berlin: Springer VerlagView ArticleGoogle Scholar
- Amoutzias GD, Robertson DL, de Peer YV, Oliver SG: Choose your partners: dimerization in eukaryotic transcription factors. Trends Biochem Sci. 2008, 33 (5): 220-229. 10.1016/j.tibs.2008.02.002View ArticlePubMedGoogle Scholar
- Markevich N, Hoek J, Kholodenko B: Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J Cell Biol. 2004, 164 (3): 353-359. 10.1083/jcb.200308060PubMed CentralView ArticlePubMedGoogle Scholar
- Wray GA, Hahn MW, Abouheif E, Balhoff JP, Pizer M, Rockman MV, Romano LA: The evolution of transcriptional regulation in eukaryotes. Mol Biol Evol. 2003, 20 (9): 1377-1419. 10.1093/molbev/msg140View ArticlePubMedGoogle Scholar
- Zhang F, Gonzalez A, Zhao M, Payne C, Lloyd A: A network of redundant bHLH proteins functions in all TTG1-dependent pathways of Arabidopsis. Development. 2003, 130 (20): 4859-4869. 10.1242/dev.00681View ArticlePubMedGoogle Scholar
- Murre C, Baltimore D: The helix-loop-helix motif: Structure and function. Transcriptional Regulation. Edited by: McKnight S, Yamamoto K. 1992, Cold Spring Harbor, New York: Cold Spring Harbor PressGoogle Scholar
- Feller A, Hernandez J, Grotewold E: An ACT-like domain participates in the dimerization of several plant basic-helix-loop-helix transcription factors. J Biol Chem. 2006, 281 (39): 28964-28974. 10.1074/jbc.M603262200View ArticlePubMedGoogle Scholar
- Aldridge BB, Haller G, Sorger PK, Lauffenburger DA: Direct Lyapunov exponent analysis enables parametric study of transient signalling governing cell behaviour. Systems biology. 2006, 153 (6): 425-32.View ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.