Abstract
Background
Modeling of molecular networks is necessary to understand their dynamical properties. While a wealth of information on molecular connectivity is available, there are still relatively few data regarding the precise stoichiometry and kinetics of the biochemical reactions underlying most molecular networks. This imbalance has limited the development of dynamical models of biological networks to a small number of wellcharacterized systems. To overcome this problem, we wanted to develop a methodology that would systematically create dynamical models of regulatory networks where the flow of information is known but the biochemical reactions are not. There are already diverse methodologies for modeling regulatory networks, but we aimed to create a method that could be completely standardized, i.e. independent of the network under study, so as to use it systematically.
Results
We developed a set of equations that can be used to translate the graph of any regulatory network into a continuous dynamical system. Furthermore, it is also possible to locate its stable steady states. The method is based on the construction of two dynamical systems for a given network, one discrete and one continuous. The stable steady states of the discrete system can be found analytically, so they are used to locate the stable steady states of the continuous system numerically. To provide an example of the applicability of the method, we used it to model the regulatory network controlling T helper cell differentiation.
Conclusion
The proposed equations have a form that permit any regulatory network to be translated into a continuous dynamical system, and also find its steady stable states. We showed that by applying the method to the T helper regulatory network it is possible to find its known states of activation, which correspond the molecular profiles observed in the precursor and effector cell types.
Background
The increasing use of high throughput technologies in different areas of biology has generated vast amounts of molecular data. This has, in turn, fueled the drive to incorporate such data into pathways and networks of interactions, so as to provide a context within which molecules operate. As a result, a wealth of connectivity information is available for multiple biological systems, and this has been used to understand some global properties of biological networks, including connectivity distribution [1], recurring motifs [2] and modularity [3]. Such information, while valuable, provides only a static snapshot of a network. For a better understanding of the functionality of a given network it is important to study its dynamical properties. The consideration of dynamics allows us to answer questions related to the number, nature and stability of the possible patterns of activation, the contribution of individual molecules or interactions to establishing such patterns, and the possibility of simulating the effects of loss or gainoffunction mutations, for example.
Mathematical modeling of metabolic networks requires specification of the biochemical reactions involved. Each reaction has to incorporate the appropriate stoichiometric coefficients to account for the principle of mass conservation. This characteristic simplifies modeling, because it implies that at equilibrium every node of the metabolic network has a total mass flux of zero [4,5]. There are cases, however, where the underlying biochemical reactions are not known for large parts of a pathway, but the direction of the flow of information is known, which is the case for socalled regulatory networks (see for example [6,7]). In these cases, the directionality of signaling is sufficient for developing mathematical models of how the patterns of activation and inhibition determine the state of activation of the network (for a review, see [8]).
When cells receive external stimuli such as hormones, mechanical forces, changes in osmolarity, membrane potential etc., there is an internal response in the form of multiple intracellular signals that may be buffered or may eventually be integrated to trigger a global cellular response, such as growth, cell division, differentiation, apoptosis, secretion etc. Modeling the underlying molecular networks as dynamical systems can capture this channeling of signals into coherent and clearly identifiable stable cellular behaviors, or cellular states. Indeed, qualitative and semiquantitative dynamical models provide valuable information about the global properties of regulatory networks. The stable steady states of a dynamical system can be interpreted as the set of all possible stable patterns of expression that can be attained within the particular biological network that is being modeled. The advantages of focusing the modeling on the stable steady states of the network are twofold. First, it reduces the quantity of experimental data required to construct a model, e.g. kinetic and rate limiting step constants, because there is no need to describe the transitory response of the network under different conditions, only the final states. Second, it is easier to verify the predictions of the model experimentally, since it requires stable cellular states to be identified; that is, longterm patterns of activation and not shortlived transitory states of activation that may be difficult to determine experimentally.
In this paper we propose a method for generating qualitative models of regulatory networks in the form of continuous dynamical systems. The method also permits the stable steady states of the system to be localized. The procedure is based on the parallel construction of two dynamical systems, one discrete and one continuous, for the same network, as summarized in Figure 1. The characteristic that distinguishes our method from others used to model regulatory networks (as summarized in [8]) is that the equations used here, and the method deployed to analyze them, are completely standardized, i.e. they are not networkspecific. This feature permits systematic application and complete automation of the whole process, thus speeding up the analysis of the dynamical properties of regulatory networks. Moreover, in contrast to methodologies for the automatic analysis of biochemical networks (as in [9]; for example), our method can be applied to networks for which there is a lack of stoichiometric information. Indeed, the method requires as sole input the information regarding the nature and directionality of the regulatory interactions. We provide an example of the applicability of our method, using it to create a dynamical model for the regulatory network that controls the differentiation of T helper (Th) cells.
Figure 1. Methodology. Schematic representation of the method for systematically constructing a dynamical model of a regulatory network and finding its stable steady states.
Results and discussion
Equations 1 and 3 (see Methods) provide the means for transforming a static graph representation of a regulatory network into two versions of a dynamical system, a discrete and a continuous description, respectively. As an example, we applied these equations to the Th regulatory network, shown in Figure 2. Briefly, the vertebrate immune system contains diverse cell populations, including antigen presenting cells, natural killer cells, and B and T lymphocytes. T lymphocytes are classified as either T helper cells (Th) or T cytotoxic cells (Tc). T helper cells take part in cell and antibodymediated immune responses by secreting various cytokines, and they are further subdivided into precursor Th0 cells and effector Th1 and Th2 cells, depending on the array of cytokines that they secrete [10]. The network that controls the differentiation from Th0 towards the Th1 or Th2 phenotypes is rather complex, and discrete modeling has been used to understand its dynamical properties [11,12]. In this work we used an updated version of the Th network, the molecular basis of which is included in the Methods. Also, we implement for the first time a continuous model of the Th network.
Figure 2. The Th network. The regulatory network that controls the differentiation process of T helper cells. Positive regulatory interactions are in green and negative interactions in red.
By applying Equation 1 to the network in Figure 2, we obtained Equation 2, which constitutes the discrete version of the dynamical system representing the Th network. Similarly, the continuous version of the Th network was obtained by applying Equation 3 to the network in Figure 2. In this case, however, some of the resulting equations are too large to be presented inside the main text, so we included them as the 1. Moreover, instead of just typing the equations, we decided to present them in a format that might be used directly to run simulations. The continuous dynamical system of the Th network is included as a plain text file that is able to run on the numerical computation software package GNU Octave http://www.octave.org webcite.
Additional File 1. The file contains the set of differential equations describing the continuous version of the Th model. It is a plain text file formatted for running simulations using the GNU Octave package http://www.octave.org webcite
Format: TXT Size: KB Download file
The high nonlinearity of Equation 3 implies that the continuous version of the dynamical model has to be studied numerically. In contrast, the discrete version can be studied analytically by using generalized logical analysis, allowing all its stable steady states to be located (see Methods). In our example, the discrete system described by Equation 2 has three stable steady states (see Table 1). Importantly, these states correspond to the molecular profiles observed in Th0, Th1 and Th2 cells. Indeed, the first stable steady state reflects the pattern of Th0 cells, which are precursor cells that do not produce any of the cytokines included in the model (IFNβ, IFNγ, IL10, IL12, IL18 and IL4). The second steady state represents Th1 cells, which show high levels of activation for IFNγ, IFNγR, SOCS1 and Tbet, and with low (although not zero) levels of JAK1 and STAT1. Finally, the third steady state corresponds to the activation observed in Th2 cells, with high levels of activation for GATA3, IL10, IL10R, IL4, IL4R, STAT3 and STAT6.
Table 1. Stable steady states of the dynamical systems. ^{a}
Equation 3 defines a highly nonlinear continuous dynamical system. In contrast with the discrete system, these continuous equations have to be studied numerically. Numerical methods for solving differential equations require the specification of an initial state, since they proceed via iterations. In our method, we propose to use the stable steady states of the discrete system as the initial states to solve the continuous system that results from application of equation 3 to a given network. We used a standard numerical simulation method to solve the continuous version of the Th model (see Methods). Starting alternatively from each of the three stable steady states found in the discrete model, i.e. the Th0, Th1 and Th2 states, the continuous system was solved numerically until it converged. The continuous system converged to values that could be compared directly with the stable steady states of the discrete system (Table 1). Note that the Th0 and Th2 stable steady states fall in exactly the same position for both the discrete and the continuous dynamical systems, and in close proximity for the Th1 state. This finding highlights the similarity in qualitative behavior of the two models constructed using equations 1 and 3, despite their different mathematical frameworks.
Despite the qualitative similarity between the discrete and continuous systems, there is no guarantee that the continuous dynamical system has only three stable steady states; there might be others without a counterpart in the discrete system. To address this possibility, we carried out a statistical study by finding the stable steady states reached by the continuous system starting from a large number of initial states. The continuous system was run 50,000 times, each time with the nodes in a random initial state within the closed interval between 0 and 1. In all cases, the system converged to one of only three different regions (Table 2), corresponding to the abovementioned Th0, Th1 and Th2 states. These results still do not eliminate the possibility that other stable steady states exist in the continuous system. Nevertheless, they show that if such additional stable steady states exist, their basin of attractions is relatively small or restricted to a small region of the state space.
Table 2. Regions of the state space reached by the continuous version of the Th model, as revealed by a large number of simulations starting from a random initial state. ^{a}
The three steady states of the continuous system are stable, since they can resist small perturbations, which create transitory responses that eventually disappear. Figure 3a shows a simulation where the system starts in its Th0 state and is then perturbed by sudden changes in the values of IFNγ and IL4 consecutively. Note that the system is capable of absorbing the perturbations, returning to the original Th0 state. If a perturbation is large enough, however, it may move the system from one stable steady state to another. If the system is in the Th0 state and IFNγ is transiently changed to it highest possible value, namely 1, the whole system reacts and moves to its Th1 state (Figure 3b). A large second perturbation by IL4, now occurring when the system is in its Th1 state, does not push the system into another stable steady state, showing the stability of the Th1 state. Conversely, if the large perturbation of IL4 occurs when the system is in the Th0 state, it moves the system towards the Th2 state (Figure 3c). In this case, a second perturbation, now in IFNγ, creates a transitory response that is not strong enough to move the system away from the Th2 state, showing the stability of this steady state. These changes from one stable steady state to another reflect the biological capacities of IFNγ and IL4 to act as key signals driving differentiation from Th0 towards Th1 and Th2 cells, respectively[10]. Furthermore, note that the Th1 and Th2 steady states are more resistant to large perturbations than the Th0 state, a characteristic that represents the stability of Th1 and Th2 cells under different experimental conditions.
Figure 3. Stability of the steady states of the continuous model of the Th network. a. The Th0 state is stable under small perturbations. b. A large perturbation on IFNγ is able to move the system from the Th0 to the Th1 steady state. This latter state is stable to perturbations. c. A large perturbation of IL4 moves the system from the Th0 state to the Th2 state, which is stable. For clarity, only the responses of key cytokines and transcription factors are plotted. The time is represented in arbitrary units.
The whole process resulted in the creation of a model with qualitative characteristics fully comparable to those found in the experimental Th system. Notably, the model used default values for all parameters. Indeed, the continuous dynamical system of the Th network has a total of 58 parameters, all of which were set to the default value of 1, and one parameter (the gain of the sigmoids) with a default value of 10. This set of default values sufficed to capture the correct qualitative behavior of the biological system, namely, the existence of three stable steady states that represent Th0, Th1 and Th2 cells. Readers can run simulations on the model by using the equations provided in the "Th_continuous_model.octave.txt" file. The file was written to allow easy modification of the initial states for the simulations, as well as the values of all parameters.
Analysis of previously published regulatory networks related to Th cell differentiation
We wanted to compare the results from our method (Figure 1) as applied to our proposed network (Figure 2) with some other similar networks. The objective of this comparison is to show that our method imposes no restrictions on the number of steady states in the models. Therefore, if the procedure is applied to wrongly reconstructed networks, the results will not reflect the general characteristics of the biological system. While there have been multiple attempts to reconstruct the signaling pathways behind the process of Th cell differentiation, they have all been carried out to describe the molecular components of the process, but not to study the dynamical behavior of the network. As a result, most of the schematic representations of these pathways are not presented as regulatory networks, but as collections of molecules with different degrees of ambiguity to describe their regulatory interactions. To circumvent this problem, we chose four pathways with low numbers of regulatory ambiguities and translated them as signaling networks (Figures 4 through 7).
Figure 4. Alternative Th network. T helper pathway published in [69], reinterpreted as a signaling network.
Figure 5. Alternative Th network. T helper pathway published in [70], reinterpreted as a signaling network.
Figure 6. Alternative Th network. T helper pathway published in [43], reinterpreted as a signaling network.
Figure 7. Alternative Th network. T helper pathway published in [71], reinterpreted as a signaling network.
The methodology introduced in this paper was applied to the four reinterpreted networks for Th cell differentiation. The stable steady states of the resulting discrete and continuous models are presented in Tables 3 through 6. Notice that none of these four alternative networks could generate the three stable steady states representing Th0, Th1 and Th2 cells. Two networks reached only two stable steady states, while two others reached more than three. Notably, all these four networks included one state representing the Th0 state, and at least one representing the Th2 state. The absence of a Th1 state in two of the networks might reflect the lack of a full characterization of the IFNγ signaling pathway at the time of writing the corresponding papers.
Table 3. Stable steady states of the signaling network in Figure 4
Table 4. Stable steady states of the signaling network in Figure 5
Table 5. Stable steady states of the signaling network in Figure 6
Table 6. Stable steady states of the signaling network in Figure 7
It is important to note that the failure of these four alternative networks to capture the three states representing Th cells is not attributable to the use of very simplistic and/or outdated data. Indeed, the network in Figure 6 comes from a relatively recent review, while that in Figure 7 is rather complex and contains five more nodes than our own proposed network (Figure 2). All this stresses the importance of using a correctly reconstructed network to develop dynamical models, either with our approach or any other.
Conclusion
There is a great deal of interest in the reconstruction and analysis of regulatory networks. Unfortunately, kinetic information about the elements that constitute a network or pathway is not easily gathered, and hence the analysis of its dynamical properties (via simulation packages such as [13]) is severely restricted to a small set of wellcharacterized systems. Moreover, the translation from a static to a dynamical representation normally requires the use of a networkspecific set of equations to represent the expression or concentration of every molecule in the system.
We herein propose a method for generating a system of ordinary differential equations to construct a model of a regulatory network. Since the equations can be unambiguously applied to any signaling or regulatory network, the construction and analysis of the model can be carried out systematically. Moreover, the process of finding the stable steady states is based on the application of an analytical method (generalized logical analysis [14,15] on a discrete version of the model), followed by a numerical method (on the continuous version) starting from specific initial states (the results obtained from the logical analysis). This characteristic allows a fully automated implementation of our methodology for modeling. In order to construct the equations of the continuous dynamical system with the exclusive use of the topological information from the network, the equations have to incorporate a set of default values for all the parameters. Therefore, the resulting model is not optimized in any sense. However, the advantage of using Equation 3 is that the user can later modify the parameters so as to refine the performance of the model, approximating it to the known behavior of the biological system under study. In this way, the user has a range of possibilities, from a purely qualitative model to one that is highly quantitative.
There are studies that compare the dynamical behavior of discrete and continuous dynamical systems. Hence, it is known that while the steady state of a Boolean model will correspond qualitatively to an analogous steady state in a continuous approach, the reverse is not necessarily true. Moreover, periodic solutions in one representation may be absent in the other [16]. This discrepancy between the discrete and continuous models is more evident for steady states where at least one of the nodes has an activation state precisely at, or near, its threshold of activation. Because of this characteristic, discrete and continuous models for a given regulatory network differ in the total number of steady states [17]. For this reason, our method focuses on the study of only one type of steady state; namely, the regular stationary points [18]. These points do not have variables near an activation threshold, and they are always stable steady states. Moreover, it has been shown that this type of stable steady state can be found in discrete models, and then used to locate their analogous states in continuous models of a given genetic regulatory network [19].
It is beyond the scope of this paper to present a detailed mathematical analysis of the dynamical system described by Equation 3. Instead, we present a framework that can help to speed up the analysis of the qualitative behavior of signaling networks. Under this perspective, the usefulness of our method will ultimately be determined through building and analyzing concrete models. To show the capabilities of our proposed methodology, we applied it to analysis of the regulatory network that controls differentiation in T helper cells. This biological system was well suited to evaluating our methodology because the network contains several known components, and it has three alternative stable patterns of activation. Moreover, it is of great interest to understand the behavior of this network, given the role of T helper cell subsets in immunity and pathology [20]. Our method applied to the Th network generated a model with the same qualitative behavior as the biological system. Specifically, the model has three stable states of activation, which can be interpreted as the states of activation found in Th0, Th1 and Th2 cells. In addition, the system is capable of being moved from the Th0 state to either the Th1 or Th2 states, given a sufficiently large IFNγ or IL4 signal, respectively. This characteristic reflects the known qualitative properties of IFNγ and IL4 as key cytokines that control the fate of T helper cell differentiation.
Regarding the numerical values returned by the model, it is not possible yet to evaluate their accuracy, given that (to our knowledge) no quantitative experimental data are available for this biological system. The resulting model, then, should be considered as a qualitative representation of the system. However, representing the nodes in the network as normalized continuous variables will eventually permit an easy comparison with quantitative experimental data whenever they become available. Towards this end, the equations in our methodology define a sigmoid function, with values ranging from 0 to 1, regardless of the values of assigned to the parameters in the equations. This characteristic has been used before to represent and model the response of signaling pathways [21,22]. It is important to note, however, that the modification of the parameters allow the model to be fitted against experimental data.
One benefit of a mathematical model of a particular biological network is the possibility of predicting the behavior of complex experimental setups. Therefore, it is important to be aware of its limitations beforehand, to avoid generating experimental data that cannot be handled by the model. The method we present in this paper has been developed to obtain the number and relative position of the stable steady states of a regulatory network. Equations 1 and 3 include a number of parameters that allow the response of the model to be finetuned, but the equations were not designed to describe the transitory responses of molecules with great detail. Therefore, failure to predict a stable steady state with high numerical accuracy should not be interpreted as a failure of the approach presented here. By contrast, failure to describe and/or predict the number and approximate location of stable steady states under a wide range of values for the parameters would call the validity of the reconstruction of a particular network into question. Here, however, it is essential to establish the validity of the network used as input. Indeed, we applied our method to four alternative forms of the network that regulates Th cell differentiation. The alternative networks (Figures 4 through 7) were taken from previously published attempts to discover the molecular basis of this differentiation process. Originally, such networks were not developed with the idea of studying dynamical properties. It is not surprising, then, that these networks do not reflect the existence of three stable steady states, representing the molecular states of Th0, Th1 and Th2 cells, respectively. In these cases, the failure to find the correct stable steady states is not a problem in the modeling methodology, but a problem in the inference of the regulatory network.
In conclusion, we have shown that the creation of a dynamical model of a regulatory network can be considerably simplified with the aid of a standardized set of equations, where the feature that distinguishes one molecule from another is the number of regulatory inputs. Such standardization permits a continuous dynamical system to be systematically and analytically constructed together with a basic analysis of its global properties, based exclusively on the information provided by the connectivity of the network. While the use of a standardized set of functions to model a network may severely restrict the capability to fit specific datasets, we believe that the loss in flexibility is balanced by the possibility of rapidly developing models and gaining knowledge of the dynamical behavior of a network, especially in those cases where few kinetic data are available. Thus, we provide a method for incorporating the dynamical perspective in the analysis of regulatory networks, using the topological information of a network, without the need to collect extensive timeseries or kinetic data.
Methods
Molecular basis of the Th network topology
The following paragraphs detail the evidence used to infer the topology of the Th regulatory network, updating the data summarized in [11]. Th1 cells are producers of IFNγ [10,23], which acts on its target cells by binding to a cellmembrane receptor [2426] to start a signaling cascade, which involves JAK1 and STAT1 [2729]. STAT1 can be activated by a number of ligands besides IFNγ, but importantly, it cannot be activated by IL4 [30], which is a major Th2 signal. In contrast, STAT1 plays a role in modulating IL4, being an intermediate in the negative regulation of IFNγ exerted on IL4 expression [31]. Different signals converge in STAT1, among them that of IFNβ/IFNβR [32]. The IFNγ signaling continues downstream to activate SOCS1 in a STAT1dependent pathway [33,34]. SOCS1, in turn, influences both the IFNγ and IL4 pathways. On the one hand, SOCS1 is a negative regulator of IFNγ signaling, blocking the interaction of IFNγR and STAT1 [35] due to direct inhibition of JAK1 [29,36]. On the other hand, SOCS1 blocks the IL4R/STAT6 pathway [37]. SOCS1 is, therefore, a key element for the inhibition from the IFNγ to the IL4 pathway. Th1 cells express high levels of SOCS1 mRNA, while it is barely detectable in Th0 and Th2 cells [38]. Finally, another key molecule is Tbet, which is a transcription factor detected in Th1 but not Th0 or Th2 cells. Tbet expression is upregulated by IFNγ in a STAT1dependent mechanism [39]. Importantly, Tbet is an inhibitor of GATA3 [40], an activator of IFNγ [40] and activator of Tbet itself [41,42].
Th2 cells express IL4, which is the major known determinant of the Th2 phenotype itself [43]. IL4 binds to its receptor, IL4R, which is preferentially expressed in Th2 cells [23,44]. The IL4R signaling is transduced by STAT6, which in turn activates GATA3 [10]. GATA3, in turn, is capable of inducing IL4 [45], thus establishing a feedback loop. The influence from the IL4 pathway on the IFNγ pathway seems to be mediated by GATA3 via STAT4 [46]. Like Tbet, GATA3 also presents a selfactivation loop [4749].
IL12 and IL18 are two molecules that affect the IFNγ pathway. IL12 is a cytokine produced by monocytes and dendritic cells and promotes the development of Th1 cells [50]. The IL12 receptor is present in its functional form in Th0 and Th1 but not Th2 cells [51]. IL12R signaling is mediated by STAT4 [52], which is able to activate IFNγ [41,46,53]. The IL12 signaling pathway can be blocked by IL4 by the STAT6 dependent downregulation of one subunit of IL12R [54]. IL18 is a cytokine produced by many cell types and promotes IFNγ production in Th cells [55]. It acts upon binding to its receptor, IL18R, which acts through IRAK [56]. IL12 and IL18 act synergistically to increase IFNγ production, but using different pathways [57,58]. Finally, IL4 is able to block IL18 signaling in a STAT6 dependent manner [59].
IL10 is a cytokine actively produced by Th2 cells, and it inhibits cytokine production by Th1 cells. As with the other cytokines mentioned above, IL10 acts upon binding to a cell surface receptor, IL10R, which in turn activates the STAT signaling system [60]. In particular, it has been shown that the functioning of IL10 signaling is dependent upon the presence of STAT3 [61]. As for the signals affecting IL10 expression, it has been shown that IL4 enhances IL10 gene expression in Th2 but not Th1 cells [62]. This requirement implies that the intracellular signaling from IL4 to IL10 should pass through a Th2 specific molecule, which from the molecules considered here can only be GATA3. Finally, IL10 has been shown to be a very powerful inhibitor of IFNγ production [60,63].
Cytokine gene expression in T cells is induced by the activation of the T cell receptor (TCR) by ligand binding. Different signaling pathways are activated by the TCR [64]. Among these is the pathway including the NFAT family of transcription factors, which are implicated in the T cell activationdependent regulation of numerous cytokines. A constitutively active form of one of the NFAT proteins, specifically NFATc1, increases the expression of IFNγ [65]. Importantly, the same experimental procedure does not affect the expression of IL4. All this indicates that the NFAT family members play a central role in the TCRinduced expression of cytokines during Th cell differentiation, especially in the Th1 pathway.
The discrete dynamical system
The discrete system represents the network as a series of interconnected elements that have only two possible states of activation, 0 (or inactive) and 1 (or active). Given this property, the network is completely described by the following set of Boolean equations:
Equation 1.
A node x in the network can have only one of three possible forms depending on whether it has activator and inhibitor input nodes, or only activators, or only inhibitors. In the first case, i.e. form § in Eqn.1, the Boolean function can be read as: x will be active in the next time step if at this time any of its activators and none of its inhibitors are acting upon it. Similarly, form §§ can be translated as: x will be active if any of its activators is acting upon it. And finally, form §§§ reads as: x will be active if none of its inhibitors are acting upon it. Note than in all cases inhibitors are strong enough to change the state of a node from 1 to 0, while activators are strong enough to change the state of a node from 0 to 1 if no inhibitor is acting on the node of reference. The three alternative forms of representing a node in Equation 1 imply two possible default states of activation, i.e. the state of a node when there are neither activators nor inhibitors acting upon it. If the connectivity of the node includes either only positive inputs, or both positive and negative inputs, then the node has an inactive state by default. Alternatively, if the connectivity of a node has only negative inputs, then the node has an active state by default.
The Th network (Figure 2) can be converted into a discrete dynamical system using Equation 1. The resulting system of equations is as follows:
Equation 2.
GATA3(t + 1) = (GATA3(t) ∨ STAT6(t)) ∧ ¬(T  bet(t))
IFN  βR(t + 1) = IFN  β(t)
IFN  γ(t + 1) = (IRAK(t) ∨ NFAT(t) ∨ STAT  4(t) ∨ T  bet(t)) ∧ ¬(STAT3(t))
IFN  γR(t + 1) = IFN  γ(t)
IL  10(t + 1) = GATA3(t)
IL  10R(t + 1) = IL  10(t)
IL  12R(t + 1) = IL  12(t)
IL  18R(t + 1) = IL  18(t) ∧ ¬(STAT6(t))
IL  4(t + 1) = GATA3(t) ∧ ¬(STAT1(t))
IL  4R(t + 1) = IL  4(t) ∧ ¬(SOCS1(t))
IRAK(t + 1) = IL  18R(t)
JAK1(t + 1) = IFN  γR(t) ∧ ¬(SOCS1(t))
NFAT(t + 1) = TCR(t)
SOCS1(t + 1) = STAT1(t) ∨ T  bet(t)
STAT1(t + 1) = IFN  βR(t) ∨ JAK1(t)
STAT3(t + 1) = IL  10R(t)
STAT4(t + 1) = IL  12R(t) ∧ ¬(GATA3(t))
STAT6(t + 1) = IL  4R(t)
T  bet(t + 1) = (STAT1(t) ∨ T  bet(t)) ∧ ¬(GATA3(t))
Notice that there are only 19 equations out of a total of 23 elements in the Th network. The reason is that four elements, namely IFNβ, IL12, IL18 and TCR, do not have inputs. These four elements are thus treated as constants, since there are no interactions that regulate their behavior. Throughout the text, these four elements are considered as having a value of 0.
Stable steady states of the discrete system
The discrete dynamical system defined by Equation 2 can be solved in different ways to find its attractors, depending on how to update the vector state X(t) to its successor, X(t+1). By far the easiest method for solving the equations is the synchronous approach (as in [66,67]). This method, however, can generate spurious results (see [14]). Hence, we use generalized logical analysis to find all the steady states of the system [15]. Generalized logical analysis allows us to find all the steady states of a discrete dynamical system by evaluating the functionality of the feedback loops, also known as circuits, in the system. In this case, the Th network (Figure 2) contains a total of 27 circuits (Table 7), 24 positive and 3 negative. Depending on the set of parameters used, positive feedback loops can generate multistationarity, while negative feedback loops can generate damped or sustained oscillations. Generalized logical analysis is a wellestablished method and the reader may find indepth explanations elsewhere [14,15,18].
Table 7. Circuits of the Th network ^{a}
The continuous dynamical system
To describe the network as a continuous dynamical system, we use the following set of ordinary differential equations:
Equation 3.
The righthand side of the differential equation comprises two parts: an activation function and a term for decay. Activation is a sigmoid function of ω, which represents the total input to the node. The equation of the sigmoid was chosen so as to pass through the two points (0,0) and (1,1), regardless of the value of its gain, h; see Figure 8. The bounding of a node x to the closed interval [0,1] implies that its level of activation should be interpreted as a normalized, not an absolute, value. This characteristic permits direct comparison between the discrete and the continuous dynamical systems, since in both formalisms the minimum and maximum levels of activation are 0 and 1. Subsequently, the second part of the equation is a decay term, which for simplicity is directly proportional to the level of activation of the node.
Figure 8. Activation of a node as a function of its total input, ω. Equation 3 ensures that the activation of a node has the form of a sigmoid, bounded in the interval [0,1] regardless of the values of h.
The total input to a node, represented by ω, is a combination of the multiple activatory and inhibitory interactions acting upon the node of reference. In the general case, different nodes have different connectivities; hence it is necessary to write a function ω so that it can describe different combinations of activatory and inhibitory inputs. For this reason, ω has three possible forms in Equation 3. If a node x_{i }is regulated by both activators and inhibitors, then the first form, §, is used. However, if is regulated exclusively by activators, form §§ is used instead. Finally, the form §§§ is used if x_{i }has only negative regulators. In all cases, the total input is a combination of weighted activators and/or inhibitors, where the weights are represented by the α and β parameters for the activators and inhibitors, respectively. The mathematical form ω was chosen so as to be monotonic and to be bounded in the closed interval [0,1] given that 0≤x≤1, α>0 and β>0. Figure 9 shows the behavior of ω when a node is controlled only by one activator. Notice that regardless of the value of α, the function is monotonically increasing and bounded to [0,1]. The reason for choosing a monotonic bounded function for ω is to preserve the sigmoid form of the total activation acting upon a node x_{i}, irrespective of the number and nature of the regulatory inputs acting upon it. Indeed, Figure 10 shows the total activation of a node x_{i }controlled by one positive regulation with different weights. Notice that the total activation retains a bounded sigmoid form independently of the value of α. This same qualitative behavior for total activation on a node x_{i }is observed if it is regulated only by inhibitors. Figure 11 shows ω as a function of one inhibitor, plotted for different strengths of interaction. In this case, the total input to x_{i }is still a bounded sigmoid regardless of the value of the parameter β (see Figure 12). This general qualitative behavior persists even with a mixture of activatory and inhibitory inputs acting upon a node. Figure 13 presents the total activation of a node x_{i }as a function of two regulatory inputs, one positive and one negative. Notice again that the equation warrants a bounded sigmoid form for the total input to a node.
Figure 9. Total input to a node, ω, as a function of one positive input, x^{a}. The value of ω is a bounded function in the interval [0,1] regardless of the interaction weight of the positive input, α.
Figure 10. Activation of a node as a function of one positive input. The activation of a node in response to one positive input, plotted for various possible interaction weights.
Figure 11. Total input to a node, ω, as a function of one negative input, x^{i}. The value of ω is a bounded function in the interval [0,1] regardless of the interaction weight of the negative input, β.
Figure 12. Activation of a node as a function of one negative input. The activation of a node in response to one negative input, plotted for various possible interaction weights.
Figure 13. Activation of a node as a function two inputs, one positive and one negative. The strength of the interactions are equal for the activation and the inhibition, α = β = 1.
Once a network is translated to a dynamical system using Equation 3, it is necessary to specify values for all parameters. For a system with n nodes and m interactions, there are m+2n parameters. However, there are usually insufficient experimental data to assign realistic values for each and every one of the parameters. Nevertheless, it is possible to use a series of default values for all the parameters in Equation 3. The reason is that, as we showed in the previous paragraph, the equations have the same qualitative shape for any value assigned to the parameters. Hence, for the sake of simplicity, it is possible to assign the same values to most of the parameters, as a first approach. For the present study on the Th model, we use a value of 1 for all αs, βs and γs; and we use h = 10, since we currently lack quantitative data to estimate more realistic values. Moreover, the use of default values ensures the possibility of creating the dynamical system in a fully automated way. Nonetheless, after the initial construction and analysis of the resulting system, the modeler may modify the values of the parameters so as to finetune the dynamical behavior of the equations, whenever more experimental quantitative data become available. The continuous dynamical system of the Th model, constructed with the use of Equation 3, yields a system of 23 equations, which is included in the file "Th_continuous_model.octave.txt".
Stable steady states of the continuous system
Nonlinear systems of ordinary differential equations are studied numerically. Hence the continuous dynamical system defined by Equation 3 poses the problem of how to find all its stable steady states without using very timeconsuming and computingintensive methods. This is where the creation of two dynamical systems of the same network, one discrete and one continuous, bears fruit. Since a Boolean (step) function is a limiting case of a very steep sigmoid curve, networks made of binary elements share many qualitative features with systems modeled using continuous functions [68]. Indeed, it has been shown [19] that the qualitative information resulted from generalized logical analysis can be directly used to find the number, nature and approximate location of the steady states of a system of differential equations representing the same network. We therefore decided to use this characteristic to speed up the process of finding all the stable steady states in the continuous dynamical system. Specifically, the stable steady states of the discrete system are used as initial states to solve the differential equations, running them until the system converges to its own stable steady states. Calculating the convergence of a system of ordinary differential equations from a given initial state is a straightforward procedure using any numerical solver. For our simulations we used the lsode function of the GNU Octave package http://www.octave.org webcite, stopping the numerical integration when all the variables of the system changed by less than 10^{4 }for at least 10 consecutive steps of the procedure. The final values of the variables in the system are considered to be the stable steady states of the continuous model of the network.
Implementation
The methodology was fully implemented in a java program, and it has been tested under a linux environment using java version 1.5.0 (JRE 5.0), as well as octave version 2.1.34. The bytecode version of the program is included as 2.
Additional File 2. The file is a java program that implements the methodology described in this paper; it requires a working installation of GNU Octave http://www.octave.org webcite. The program takes as input a plain text file containing the topology of the network to analyze, with the following format: MoleculeA > MoleculeB MoleculeB  MoleculeA The output of the program is a stream of plain text formatted for GNU Octave.
Format: JAR Size: 33KB Download file
Competing interests
The author(s) declare that they have no competing interests.
Authors' contributions
LM inferred the regulatory network, created the equations, developed the methods and wrote the paper. IX made a substantial contribution to the design and development of the methods, revised the intellectual content, and helped in drafting the manuscript.
Acknowledgements
We want to thank Massimo de Francesco, Mark Ibberson, Caroline JohnsonLeger, Maria Karmirantzou, Lukasz Salwinski, François Talabot and Francisca Zanoguera for their valuable comments and suggestions.
References

Jeong H, Tombor B, Albert R, Oltvai ZN, Barabási AL: The largescale organization of metabolic networks.
Nature 2000, 407:651654. PubMed Abstract  Publisher Full Text

Milo R, ShenOrr S, Itzkovitz S, Kashtan N, Chklovskii D, Alon U: Network motifs: Simple building blocks of complex networks.
Science 2002, 298:824827. PubMed Abstract  Publisher Full Text

Ravasz E, Somera AL, Mongru DA, Oltvai ZN, Barabási AL: Hierarchical organization of modularity in metabolic networks.
Science 2002, 297:15511555. PubMed Abstract  Publisher Full Text

Covert MW, Schilling CH, Famili I, Edwards JS, Goryanin II, Selkov E, Palsson BO: Metabolic modeling of microbial strains in silico.
Trends Biochem Sci 2001, 26:179186. PubMed Abstract  Publisher Full Text

Herrgård MJ, Covert MW, Palsson BØ: Reconstruction of microbial transcriptional regulatory networks.
Curr Opin Biotechnol 2004, 15:7077. PubMed Abstract  Publisher Full Text

Mendoza L, Thieffry D, AlvarezBuylla ER: Genetic control of flower morphogenesis in Arabidopsis thaliana: a logical analysis.
Bioinformatics 1999, 15:593606. PubMed Abstract  Publisher Full Text

Sánchez L, Thieffry D: Segmenting the fly embryo: a logical analysis of the pairrule crossregulatory module.
J Theor Biol 2003, 224:517537. PubMed Abstract  Publisher Full Text

de Jong H: Modeling and simulation of genetic regulatory systems: a literature review.
J Comp Biol 2002, 9:67103. Publisher Full Text

Lok L, Brent R: Automatic generation of cellular reaction networks with Moleculizer 1.0.
Nature Biotechnol 2005, 23:131136. Publisher Full Text

Murphy KM, Reiner SL: The lineage decisions on helper T cells.
Nat Rev Immunol 2002, 2:933944. PubMed Abstract  Publisher Full Text

Mendoza L: A network model for the control of the differentiation process in Th cells.

Remy E, Ruet P, Mendoza L, Thieffry D, Chaouiya C: From Logical Regulatory Graphs to Standard Petri Nets: Dynamic Roles and Functionality of Feedback Circuits.

Mendes P: Biochemistry by numbers: simulation of biochemical pathways with Gepasi 3.
Trends Biochem Sci 1997, 22:361363. PubMed Abstract  Publisher Full Text

Thomas R: Regulatory networks seen as asynchronous automata: a logical description.

Thomas R, Thieffry D, Kaufman M: Dynamical behaviour of biological regulatory networksI. Biological role of feedback loops and practical use of the concept of the loopcharacteristic state.
Bull Math Biol 1995, 57:247276. PubMed Abstract  Publisher Full Text

Glass L, Kaufman M: The logical analysis of continuous, nonlinear biochemical control networks.
J Theor Biol 1973, 39:103129. PubMed Abstract  Publisher Full Text

Mochizuki A: An analytical study of the number of steady states in gene regulatory networks.
J Theor Biol 2005, 236:291310. PubMed Abstract  Publisher Full Text

Snoussi EH, Thomas R: Logical identification of all steady states: the concept of feedback loop characteristic states.
Bull Math Biol 1993, 55:973991. Publisher Full Text

Muraille E, Thieffry D, Leo O, Kaufman M: Toxicity and neuroendocrine regulation of the immune response: a model analysis.
J Theor Biol 1996, 183:285305. PubMed Abstract  Publisher Full Text

Singh VK, Mehrotra S, Agarwal SS: The paradigm of Th1 and Th2 cytokines: its relevance to autoimmunity and allergy.
Immunol Res 1999, 20:147161. PubMed Abstract

Hautaniemi S, Kharait S, Iwabu A, Wells A, Lauffenburger DA: Modeling of signalresponse cascades using decision tree analysis.
Bioinformatics 2005, 21:20272035. PubMed Abstract  Publisher Full Text

Sachs K, Perez O, Pe'er D, Lauffenburger DA, Nolan GP: Causal proteinsignaling networks derived from multiparameter singlecell data.
Science 2005, 308:523529. PubMed Abstract  Publisher Full Text

Hamalainen H, Zhou H, Chou W, Hashizume H, Heller R, Lahesmaa R: Distinct gene expression profiles of human type 1 and type 2 T helper cells.
Genome Biology 2001, 2:10022. BioMed Central Full Text

Groux H, Sornasse T, Cottrez F, de Vries JE, Coffman RL, Roncarolo MG, Yssel H: Induction of human T helper cell type 1 differentiation results in loss of IFNγ receptor βchain expression.
J Immunol 1997, 158:56275631. PubMed Abstract  Publisher Full Text

Novelli F, D'Elios MM, Bernabei P, Ozmen L, Rigamonti L, Almerigogna F, Forni G, Del Prete G: Expression and role in apoptosis of the α and βchains of the IFNγ receptor in human Th1 and Th2 clones.
J Immunol 1997, 159:206213. PubMed Abstract  Publisher Full Text

Rigamonti L, Ariotti S, Losana G, Gradini R, Russo MA, Jouanguy E, Casanova JL, Forni G, Novelli F: Surface expression of the IFNγR2 chain is regulated by intracellular trafficking in human T lymphocytes.
J Immunol 2000, 164:201207. PubMed Abstract  Publisher Full Text

Kotenko SV, Pestka S: JakStat signal transduction pathway through the eyes of cytokine class II receptor complexes.
Oncogene 2000, 19:25572565. PubMed Abstract  Publisher Full Text

Kerr IM, CostaPereira AP, Lillemeier BF, Strobl B: Of JAKs, STATs, blind watchmakers, jeeps and trains.
FEBS Lett 2003, 546:15. PubMed Abstract  Publisher Full Text

Krebs DL, Hilton DJ: SOCS proteins: negative regulators of cytokine signaling.
Stem Cells 2001, 19:378387. PubMed Abstract  Publisher Full Text

Moriggl R, Kristofic C, Kinzel B, Volarevic S, Groner B, Brinkmann V: Activation of STAT proteins and cytokine genes in human Th1 and Th2 cells generated in the absence of IL12 and IL4.
J Immunol 1998, 160:33853392. PubMed Abstract  Publisher Full Text

Elser B, Lohoff M, Kock S, Giaisi M, Kirchhoff S, Krammer PH, LiWeber M: IFNγ represses IL4 expression via IRF1 and IRF2.
Immunity 2002, 17:703712. PubMed Abstract  Publisher Full Text

Goodbourn S, Didcock L, Randal RE: Interferons: cell signalling, immune modulation, antiviral responses and virus countermeasures.
J Gen Virol 2000, 81:23412364. PubMed Abstract  Publisher Full Text

Chen XP, Losman JA, Rothman P: SOCS proteins, regulators of intracellular signaling.
Immunity 2000, 13:287290. PubMed Abstract  Publisher Full Text

Saito H, Morita Y, Fujimoto M, Narazaki M, Naka T, Kishimoto T: IFN regulatory factor1mediated transcriptional activation of mouse STATinduced STAT inhibitor1 gene promoter by IFNγ.
J Immunol 2000, 164:58335843. PubMed Abstract  Publisher Full Text

Diehl S, Anguita J, Hoffmeyer A, Zapton T, Ihle JN, Fikrig E, Rincón M: Inhibition of Th1 differentiation by IL6 is mediated by SOCS1.
Immunity 2000, 13:805815. PubMed Abstract  Publisher Full Text

Zhang JG, Metcalf D, Rakar S, Asimakis M, Greenhalgh CJ, Willson TA, Starr R, Nicholson SE, Carter W, Alexander WS, Hilton J: The SOCS box of suppressor of cytokine signaling1 is important for inhibition of cytokine action in vivo.
Proc Natl Acad Sci USA 2001, 98:1326113265. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Losman JA, Chen XP, Hilton D, Rothman P: Cutting edge: SOCS1 is a potent inhibitor of IL4 signal transduction.
J Immunol 1999, 162:37703774. PubMed Abstract  Publisher Full Text

Egwuagu CE, Yu CR, Zhang M, Mahdi RM, Kim SJ, Gery I: Suppressors of cytokine signaling proteins are differentially expresses in Th1 and Th2 cells: implications for the Th cell lineage commitment and maintenance.
J Immunol 2002, 168:31813187. PubMed Abstract  Publisher Full Text

Lighvani A., Frucht DM, Jankovic D, Yamane H, Aliberti J, Hissong BD, Nguyen BV, Gadina M, Sher A, Paul WE, O'Shea JJ: Tbet is rapidly induced by interferonγ in lymphoid and myeloid cells.
Proc Natl Acad Sci USA 2001, 98:1513715142. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Szabo SJ, Kim ST, Costa GL, Zhang X, Fathman CG, Glimcher LH: A novel transcription factor, Tbet, directs Th1 lineage commitment.
Cell 2000, 100:655669. PubMed Abstract  Publisher Full Text

Mullen AC, High FA, Hutchins AS, Lee HW, Villarino AV, Livingston DM, Kung AL, Cereb N, Yao TP, Yang SY, Reiner SL: Role of Tbet in commitment of TH1 cells before IL12dependent selection.
Science 2001, 292:19071910. PubMed Abstract  Publisher Full Text

Zhang Y, Apilado R, Coleman J, BenSasson S, Tsang S, HuLi J, Paul WE, Huang H: Interferon γ stabilizes the T helper cell type 1 phenotype.
J Exp Med 2001, 194:165172. PubMed Abstract  Publisher Full Text

Agnello D, Lankford CSR, Bream J, Morinobu A, Gadina M, O'Shea J, Frucht DM: Cytokines and transcription factors that regulate T helper cell differentiation: new players and new insights.
J Clin Immunol 2003, 23:147161. PubMed Abstract  Publisher Full Text

Nelms K, Keegan AD, Zamorano J, Ryan JJ, Paul WE: The IL4 receptor: signaling mechanisms and biologic functions.
Annu Rev Immunol 1999, 17:701738. PubMed Abstract  Publisher Full Text

Ouyang W, Löhning M, Gao Z, Assenmacher M, Ranganath S, Radbruch A, Murphy KM: Stat6independent GATA3 autoactivation directs IL4independent Th2 development and commitment.
Immunity 2000, 12:2737. PubMed Abstract  Publisher Full Text

Usui T, Nishikomori R, Kitani A, Strober W: GATA3 suppresses Th1 development by downregulation of Stat4 and not through effects on IL12Rβ2 chain or Tbet.
Immunity 2003, 18:415428. PubMed Abstract  Publisher Full Text

Zhou M, Ouyang W, Gong Q, Katz SG, White JM, Orkin SH, Murphy KM: Friend of GATA1 Represses GATA3dependent activity in CD4+ cells.
J Exp Med 2001, 194:14611471. PubMed Abstract  Publisher Full Text

Zhou M, Ouyang W: The function role of GATA3 in Th1 and Th2 differentiation.
Immunol Res 2003, 28:2537. PubMed Abstract  Publisher Full Text

Höfer T, Nathansen H, Löhning M, Radbruch A, Heinrich R: GATA3 transcriptional imprinting in Th2 lymphocytes: A mathematical model.
Proc Natl Acad Sci USA 2002, 99:93649368. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Trinchieri G: Interleukin12: a proinflammatory cytokine with immunoregulatory functions that bridge innate resistance and antigenspecific adaptive immunity.
Annu Rev Immunol 1995, 13:251276. PubMed Abstract  Publisher Full Text

Szabo SJ, Jacobson NG, Dighe AS, Gubler U, Murphy KM: Developmental commitment to the Th2 lineage by extinction of IL12 signaling.
Immunity 1995, 2:665675. PubMed Abstract  Publisher Full Text

Thierfelder WE, van Deursen JM, Yamamoto K, Tripp RA, Sarawar SR, Carson RT, Sangster MY, Vignali DA, Doherty PC, Grosveld GC, Ihle JN: Requirement for Stat4 in interleukin12mediated responses of natural killer and T cells.
Nature 1996, 382:171174. PubMed Abstract  Publisher Full Text

Kaplan MH, Sun YL, Hoey T, Grusby MJ: Impaired IL12 responses and enhanced development of Th2 cells in Stat4deficient mice.
Nature 1996, 382:174177. PubMed Abstract  Publisher Full Text

Szabo SJ, Dighe AS, Gubler U, Murphy KM: Regulation of the interleukin (IL)12R β2 subunit expression in developing T helper 1 (Th1) and Th2 cells.
J Exp Med 1997, 185:817824. PubMed Abstract  Publisher Full Text

Swain SL: Interleukin 18: tipping the balance towards a T helper cell 1 response.
J Exp Med 2001, 194:F11F14. PubMed Abstract  Publisher Full Text

Chang JT, Segal BM, Nakanishi K, Okamura H, Shevach EM: The costimulatory effect of IL18 on the induction of antigenspecific IFNgamma production by resting T cells is IL12 dependent and is mediated by upregulation of the IL12 receptor beta2 subunit.
Eur J Immunol 2000, 30:11131119. PubMed Abstract  Publisher Full Text

Akira S: The role of IL18 in innate immunity.
Curr Opin Immunol 2000, 12:5963. PubMed Abstract  Publisher Full Text

Kanakaraj P, Ngo K, Wu Y, Angulo A, Ghazal P, Harris CA, Siekierka JJ, Peterson PA, FungLeung WP: Defective interleukin (IL)18mediated natural killer and T helper cell type 1 response in IL1 receptorassociated kinase (IRAK)deficient mice.
J Exp Med 1999, 189:11291138. PubMed Abstract  Publisher Full Text

Smeltz RB, Chen J, HuLi J, Shevach EM: Regulation of interleukin (IL)18 receptor α chain expression on CD4+ T cells during T helper (Th)1/Th2 differentiation: critical downregulatory role of IL4.
J Exp Med 2001, 194:143153. PubMed Abstract  Publisher Full Text

Moore KW, de Waal Malefyt R, Coffman R, O'Garra A: Interleukin10 and the interleukin10 receptor.
Annu Rev Immunol 2001, 19:683765. PubMed Abstract  Publisher Full Text

Riley JK, Takeda K, Akira S, Schreiber RD: Interleukin10 receptor signaling through the JAKSTAT pathway.
J Biol Chem 1999, 274:1651316521. PubMed Abstract  Publisher Full Text

SchmidtWeber CB, Alexander SI, Henault LE, James L, Lichtman AH: IL4 enhances IL10 gene expression in murine Th2 cells in the absence of TCR engagement.
J Immunol 1999, 162:238244. PubMed Abstract  Publisher Full Text

Skapenko A, Niedobitek GU, Kalden JR, Lipsky PE, SchulzeKoops H: Generation and regulation of human Th1biased immune response in vivo: A critical role for IL4 and IL10.
J Immunol 2004, 172:64276434. PubMed Abstract  Publisher Full Text

Huang Y, Wange RL: T cell receptor signaling: beyond complex complexes.
J Biol Chem 2004, 279:2882728830. PubMed Abstract  Publisher Full Text

Porter CM, Clipstone NA: Sustained NFAT signaling promotes a Th1like pattern of gene expression in primary murine CD4^{+ }T cells.
J Immunol 2002, 168:49364945. PubMed Abstract  Publisher Full Text

Kauffman SA: Metabolic stability and epigenesis in randomly constructed genetic nets.
J Theor Biol 1969, 22:437467. PubMed Abstract  Publisher Full Text

Kauffman SA: Antichaos and adaptation.
Sci Am 1991, 265:7884. PubMed Abstract

Thomas R: Laws for the dynamics of regulatory networks.
Int J Dev Biol 1998, 42:479485. PubMed Abstract  Publisher Full Text

Muraille E, Leo O: Revisiting the Th1/Th2 paradigm.
Scand J Immunol 1998, 47:19. PubMed Abstract  Publisher Full Text

Street NE, Mosmann TM: Functional diversity of T lymphocytes due to secretion of different cytokine patterns.
FASEB J 1991, 5:171177. PubMed Abstract  Publisher Full Text

Murphy KM, Ouyang W, Farrar JD, Yang J, Ranganath S, Asnagli H, Afkarian M, Murphy TL: Signaling and transcription in T helper development.
Annu Rev Immunol 2000, 18:451494. PubMed Abstract  Publisher Full Text