- Open Access
A method for the generation of standardized qualitative dynamical systems of regulatory networks
Theoretical Biology and Medical Modellingvolume 3, Article number: 13 (2006)
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 well-characterized 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.
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.
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.
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 , recurring motifs  and modularity . 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 gain-of-function 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 so-called 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 ).
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 semi-quantitative 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 two-fold. 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, long-term patterns of activation and not short-lived 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 ) is that the equations used here, and the method deployed to analyze them, are completely standardized, i.e. they are not network-specific. 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 ; 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.
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 antibody-mediated immune responses by secreting various cytokines, and they are further sub-divided into precursor Th0 cells and effector Th1 and Th2 cells, depending on the array of cytokines that they secrete . 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.
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 Additional file 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.
The high non-linearity 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-γ, IL-10, IL-12, IL-18 and IL-4). The second steady state represents Th1 cells, which show high levels of activation for IFN-γ, IFN-γR, SOCS1 and T-bet, 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, IL-10, IL-10R, IL-4, IL-4R, STAT3 and STAT6.
Equation 3 defines a highly non-linear 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 above-mentioned 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.
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 IL-4 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 IL-4, 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 IL-4 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 IL-4 to act as key signals driving differentiation from Th0 towards Th1 and Th2 cells, respectively. 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.
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).
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.
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.
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 ) is severely restricted to a small set of well-characterized systems. Moreover, the translation from a static to a dynamical representation normally requires the use of a network-specific 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 . 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 . For this reason, our method focuses on the study of only one type of steady state; namely, the regular stationary points . 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 .
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 . 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 IL-4 signal, respectively. This characteristic reflects the known qualitative properties of IFN-γ and IL-4 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 fine-tuned, 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 time-series or kinetic data.
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 . Th1 cells are producers of IFN-γ [10, 23], which acts on its target cells by binding to a cell-membrane receptor [24–26] to start a signaling cascade, which involves JAK1 and STAT-1 [27–29]. STAT-1 can be activated by a number of ligands besides IFN-γ, but importantly, it cannot be activated by IL-4 , which is a major Th2 signal. In contrast, STAT-1 plays a role in modulating IL-4, being an intermediate in the negative regulation of IFN-γ exerted on IL-4 expression . Different signals converge in STAT-1, among them that of IFN-β/IFN-βR . The IFN-γ signaling continues downstream to activate SOCS-1 in a STAT-1-dependent pathway [33, 34]. SOCS-1, in turn, influences both the IFN-γ and IL-4 pathways. On the one hand, SOCS-1 is a negative regulator of IFN-γ signaling, blocking the interaction of IFN-γR and STAT-1  due to direct inhibition of JAK1 [29, 36]. On the other hand, SOCS-1 blocks the IL-4R/STAT-6 pathway . SOCS-1 is, therefore, a key element for the inhibition from the IFN-γ to the IL-4 pathway. Th1 cells express high levels of SOCS-1 mRNA, while it is barely detectable in Th0 and Th2 cells . Finally, another key molecule is T-bet, which is a transcription factor detected in Th1 but not Th0 or Th2 cells. T-bet expression is upregulated by IFN-γ in a STAT-1-dependent mechanism . Importantly, T-bet is an inhibitor of GATA-3 , an activator of IFN-γ  and activator of T-bet itself [41, 42].
Th2 cells express IL-4, which is the major known determinant of the Th2 phenotype itself . IL-4 binds to its receptor, IL-4R, which is preferentially expressed in Th2 cells [23, 44]. The IL-4R signaling is transduced by STAT-6, which in turn activates GATA-3 . GATA-3, in turn, is capable of inducing IL-4 , thus establishing a feedback loop. The influence from the IL-4 pathway on the IFN-γ pathway seems to be mediated by GATA-3 via STAT-4 . Like T-bet, GATA-3 also presents a self-activation loop [47–49].
IL-12 and IL-18 are two molecules that affect the IFN-γ pathway. IL-12 is a cytokine produced by monocytes and dendritic cells and promotes the development of Th1 cells . The IL-12 receptor is present in its functional form in Th0 and Th1 but not Th2 cells . IL-12R signaling is mediated by STAT-4 , which is able to activate IFN-γ [41, 46, 53]. The IL-12 signaling pathway can be blocked by IL-4 by the STAT-6 dependent down-regulation of one subunit of IL-12R . IL-18 is a cytokine produced by many cell types and promotes IFN-γ production in Th cells . It acts upon binding to its receptor, IL-18R, which acts through IRAK . IL-12 and IL-18 act synergistically to increase IFN-γ production, but using different pathways [57, 58]. Finally, IL-4 is able to block IL-18 signaling in a STAT-6 dependent manner .
IL-10 is a cytokine actively produced by Th2 cells, and it inhibits cytokine production by Th1 cells. As with the other cytokines mentioned above, IL-10 acts upon binding to a cell surface receptor, IL-10R, which in turn activates the STAT signaling system . In particular, it has been shown that the functioning of IL-10 signaling is dependent upon the presence of STAT-3 . As for the signals affecting IL-10 expression, it has been shown that IL-4 enhances IL-10 gene expression in Th2 but not Th1 cells . This requirement implies that the intracellular signaling from IL-4 to IL-10 should pass through a Th2 specific molecule, which from the molecules considered here can only be GATA-3. Finally, IL-10 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 . Among these is the pathway including the NFAT family of transcription factors, which are implicated in the T cell activation-dependent regulation of numerous cytokines. A constitutively active form of one of the NFAT proteins, specifically NFATc1, increases the expression of IFN-γ . Importantly, the same experimental procedure does not affect the expression of IL-4. All this indicates that the NFAT family members play a central role in the TCR-induced 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:
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.
GATA 3(t + 1) = (GATA 3(t) ∨ STAT 6(t)) ∧ ¬(T - bet(t))
IFN - βR(t + 1) = IFN - β(t)
IFN - γ(t + 1) = (IRAK(t) ∨ NFAT(t) ∨ STAT - 4(t) ∨ T - bet(t)) ∧ ¬(STAT 3(t))
IFN - γR(t + 1) = IFN - γ(t)
IL - 10(t + 1) = GATA 3(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) = GATA 3(t) ∧ ¬(STAT 1(t))
IL - 4R(t + 1) = IL - 4(t) ∧ ¬(SOCS 1(t))
IRAK(t + 1) = IL - 18R(t)
JAK 1(t + 1) = IFN - γR(t) ∧ ¬(SOCS 1(t))
NFAT(t + 1) = TCR(t)
SOCS 1(t + 1) = STAT 1(t) ∨ T - bet(t)
STAT 1(t + 1) = IFN - βR(t) ∨ JAK 1(t)
STAT 3(t + 1) = IL - 10R(t)
STAT 4(t + 1) = IL - 12R(t) ∧ ¬(GATA 3(t))
STAT 6(t + 1) = IL - 4R(t)
T - bet(t + 1) = (STAT 1(t) ∨ T - bet(t)) ∧ ¬(GATA 3(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-β, IL-12, IL-18 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 ). Hence, we use generalized logical analysis to find all the steady states of the system . 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 well-established method and the reader may find in-depth explanations elsewhere [14, 15, 18].
The continuous dynamical system
To describe the network as a continuous dynamical system, we use the following set of ordinary differential equations:
The right-hand 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.
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.
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 fine-tune 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 time-consuming and computing-intensive 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 . Indeed, it has been shown  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, 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.
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 Additional file 2.
Jeong H, Tombor B, Albert R, Oltvai ZN, Barabási AL: The large-scale organization of metabolic networks. Nature. 2000, 407: 651-654. 10.1038/35036627.
Milo R, Shen-Orr S, Itzkovitz S, Kashtan N, Chklovskii D, Alon U: Network motifs: Simple building blocks of complex networks. Science. 2002, 298: 824-827. 10.1126/science.298.5594.824.
Ravasz E, Somera AL, Mongru DA, Oltvai ZN, Barabási AL: Hierarchical organization of modularity in metabolic networks. Science. 2002, 297: 1551-1555. 10.1126/science.1073374.
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: 179-186. 10.1016/S0968-0004(00)01754-0.
Herrgård MJ, Covert MW, Palsson BØ: Reconstruction of microbial transcriptional regulatory networks. Curr Opin Biotechnol. 2004, 15: 70-77. 10.1016/j.copbio.2003.11.002.
Mendoza L, Thieffry D, Alvarez-Buylla ER: Genetic control of flower morphogenesis in Arabidopsis thaliana: a logical analysis. Bioinformatics. 1999, 15: 593-606. 10.1093/bioinformatics/15.7.593.
Sánchez L, Thieffry D: Segmenting the fly embryo: a logical analysis of the pair-rule cross-regulatory module. J Theor Biol. 2003, 224: 517-537. 10.1016/S0022-5193(03)00201-7.
de Jong H: Modeling and simulation of genetic regulatory systems: a literature review. J Comp Biol. 2002, 9: 67-103. 10.1089/10665270252833208.
Lok L, Brent R: Automatic generation of cellular reaction networks with Moleculizer 1.0. Nature Biotechnol. 2005, 23: 131-136. 10.1038/nbt1054.
Murphy KM, Reiner SL: The lineage decisions on helper T cells. Nat Rev Immunol. 2002, 2: 933-944. 10.1038/nri954.
Mendoza L: A network model for the control of the differentiation process in Th cells. BioSystems.
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. Transactions on Computational Systems Biology.
Mendes P: Biochemistry by numbers: simulation of biochemical pathways with Gepasi 3. Trends Biochem Sci. 1997, 22: 361-363. 10.1016/S0968-0004(97)01103-1.
Thomas R: Regulatory networks seen as asynchronous automata: a logical description. J Theor Biol. 1991, 153: 1-23.
Thomas R, Thieffry D, Kaufman M: Dynamical behaviour of biological regulatory networks-I. Biological role of feedback loops and practical use of the concept of the loop-characteristic state. Bull Math Biol. 1995, 57: 247-276. 10.1016/0092-8240(94)00036-C.
Glass L, Kaufman M: The logical analysis of continuous, non-linear biochemical control networks. J Theor Biol. 1973, 39: 103-129. 10.1016/0022-5193(73)90208-7.
Mochizuki A: An analytical study of the number of steady states in gene regulatory networks. J Theor Biol. 2005, 236: 291-310. 10.1016/j.jtbi.2005.03.015.
Snoussi EH, Thomas R: Logical identification of all steady states: the concept of feedback loop characteristic states. Bull Math Biol. 1993, 55: 973-991. 10.1016/S0092-8240(05)80199-5.
Muraille E, Thieffry D, Leo O, Kaufman M: Toxicity and neuroendocrine regulation of the immune response: a model analysis. J Theor Biol. 1996, 183: 285-305. 10.1006/jtbi.1996.0221.
Singh VK, Mehrotra S, Agarwal SS: The paradigm of Th1 and Th2 cytokines: its relevance to autoimmunity and allergy. Immunol Res. 1999, 20: 147-161.
Hautaniemi S, Kharait S, Iwabu A, Wells A, Lauffenburger DA: Modeling of signal-response cascades using decision tree analysis. Bioinformatics. 2005, 21: 2027-2035. 10.1093/bioinformatics/bti278.
Sachs K, Perez O, Pe'er D, Lauffenburger DA, Nolan GP: Causal protein-signaling networks derived from multiparameter single-cell data. Science. 2005, 308: 523-529. 10.1126/science.1105809.
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: 1-0022. 10.1186/gb-2001-2-7-research0022.
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: 5627-5631.
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: 206-213.
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: 201-207.
Kotenko SV, Pestka S: Jak-Stat signal transduction pathway through the eyes of cytokine class II receptor complexes. Oncogene. 2000, 19: 2557-2565. 10.1038/sj.onc.1203524.
Kerr IM, Costa-Pereira AP, Lillemeier BF, Strobl B: Of JAKs, STATs, blind watchmakers, jeeps and trains. FEBS Lett. 2003, 546: 1-5. 10.1016/S0014-5793(03)00411-3.
Krebs DL, Hilton DJ: SOCS proteins: negative regulators of cytokine signaling. Stem Cells. 2001, 19: 378-387. 10.1634/stemcells.19-5-378.
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 IL-12 and IL-4. J Immunol. 1998, 160: 3385-3392.
Elser B, Lohoff M, Kock S, Giaisi M, Kirchhoff S, Krammer PH, Li-Weber M: IFN-γ represses IL-4 expression via IRF-1 and IRF-2. Immunity. 2002, 17: 703-712. 10.1016/S1074-7613(02)00471-5.
Goodbourn S, Didcock L, Randal RE: Interferons: cell signalling, immune modulation, antiviral responses and virus countermeasures. J Gen Virol. 2000, 81: 2341-2364.
Chen XP, Losman JA, Rothman P: SOCS proteins, regulators of intracellular signaling. Immunity. 2000, 13: 287-290. 10.1016/S1074-7613(00)00028-5.
Saito H, Morita Y, Fujimoto M, Narazaki M, Naka T, Kishimoto T: IFN regulatory factor-1-mediated transcriptional activation of mouse STAT-induced STAT inhibitor-1 gene promoter by IFN-γ. J Immunol. 2000, 164: 5833-5843.
Diehl S, Anguita J, Hoffmeyer A, Zapton T, Ihle JN, Fikrig E, Rincón M: Inhibition of Th1 differentiation by IL-6 is mediated by SOCS1. Immunity. 2000, 13: 805-815. 10.1016/S1074-7613(00)00078-9.
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 signaling-1 is important for inhibition of cytokine action in vivo. Proc Natl Acad Sci USA. 2001, 98: 13261-13265. 10.1073/pnas.231486498.
Losman JA, Chen XP, Hilton D, Rothman P: Cutting edge: SOCS-1 is a potent inhibitor of IL-4 signal transduction. J Immunol. 1999, 162: 3770-3774.
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: 3181-3187.
Lighvani A., Frucht DM, Jankovic D, Yamane H, Aliberti J, Hissong BD, Nguyen BV, Gadina M, Sher A, Paul WE, O'Shea JJ: T-bet is rapidly induced by interferon-γ in lymphoid and myeloid cells. Proc Natl Acad Sci USA. 2001, 98: 15137-15142. 10.1073/pnas.261570598.
Szabo SJ, Kim ST, Costa GL, Zhang X, Fathman CG, Glimcher LH: A novel transcription factor, T-bet, directs Th1 lineage commitment. Cell. 2000, 100: 655-669. 10.1016/S0092-8674(00)80702-3.
Mullen AC, High FA, Hutchins AS, Lee HW, Villarino AV, Livingston DM, Kung AL, Cereb N, Yao TP, Yang SY, Reiner SL: Role of T-bet in commitment of TH1 cells before IL-12-dependent selection. Science. 2001, 292: 1907-1910. 10.1126/science.1059835.
Zhang Y, Apilado R, Coleman J, Ben-Sasson S, Tsang S, Hu-Li J, Paul WE, Huang H: Interferon γ stabilizes the T helper cell type 1 phenotype. J Exp Med. 2001, 194: 165-172. 10.1084/jem.194.2.165.
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: 147-161. 10.1023/A:1023381027062.
Nelms K, Keegan AD, Zamorano J, Ryan JJ, Paul WE: The IL-4 receptor: signaling mechanisms and biologic functions. Annu Rev Immunol. 1999, 17: 701-738. 10.1146/annurev.immunol.17.1.701.
Ouyang W, Löhning M, Gao Z, Assenmacher M, Ranganath S, Radbruch A, Murphy KM: Stat6-independent GATA-3 autoactivation directs IL-4-independent Th2 development and commitment. Immunity. 2000, 12: 27-37. 10.1016/S1074-7613(00)80156-9.
Usui T, Nishikomori R, Kitani A, Strober W: GATA-3 suppresses Th1 development by downregulation of Stat4 and not through effects on IL-12Rβ2 chain or T-bet. Immunity. 2003, 18: 415-428. 10.1016/S1074-7613(03)00057-8.
Zhou M, Ouyang W, Gong Q, Katz SG, White JM, Orkin SH, Murphy KM: Friend of GATA-1 Represses GATA-3-dependent activity in CD4+ cells. J Exp Med. 2001, 194: 1461-1471. 10.1084/jem.194.10.1461.
Zhou M, Ouyang W: The function role of GATA-3 in Th1 and Th2 differentiation. Immunol Res. 2003, 28: 25-37. 10.1385/IR:28:1:25.
Höfer T, Nathansen H, Löhning M, Radbruch A, Heinrich R: GATA-3 transcriptional imprinting in Th2 lymphocytes: A mathematical model. Proc Natl Acad Sci USA. 2002, 99: 9364-9368. 10.1073/pnas.142284699.
Trinchieri G: Interleukin-12: a proinflammatory cytokine with immunoregulatory functions that bridge innate resistance and antigen-specific adaptive immunity. Annu Rev Immunol. 1995, 13: 251-276.
Szabo SJ, Jacobson NG, Dighe AS, Gubler U, Murphy KM: Developmental commitment to the Th2 lineage by extinction of IL-12 signaling. Immunity. 1995, 2: 665-675. 10.1016/1074-7613(95)90011-X.
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 interleukin-12-mediated responses of natural killer and T cells. Nature. 1996, 382: 171-174. 10.1038/382171a0.
Kaplan MH, Sun YL, Hoey T, Grusby MJ: Impaired IL-12 responses and enhanced development of Th2 cells in Stat4-deficient mice. Nature. 1996, 382: 174-177. 10.1038/382174a0.
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: 817-824. 10.1084/jem.185.5.817.
Swain SL: Interleukin 18: tipping the balance towards a T helper cell 1 response. J Exp Med. 2001, 194: F11-F14. 10.1084/jem.194.3.F11.
Chang JT, Segal BM, Nakanishi K, Okamura H, Shevach EM: The costimulatory effect of IL-18 on the induction of antigen-specific IFN-gamma production by resting T cells is IL-12 dependent and is mediated by up-regulation of the IL-12 receptor beta2 subunit. Eur J Immunol. 2000, 30: 1113-1119. 10.1002/(SICI)1521-4141(200004)30:4<1113::AID-IMMU1113>3.0.CO;2-P.
Akira S: The role of IL-18 in innate immunity. Curr Opin Immunol. 2000, 12: 59-63. 10.1016/S0952-7915(99)00051-5.
Kanakaraj P, Ngo K, Wu Y, Angulo A, Ghazal P, Harris CA, Siekierka JJ, Peterson PA, Fung-Leung WP: Defective interleukin (IL)-18-mediated natural killer and T helper cell type 1 response in IL-1 receptor-associated kinase (IRAK)-deficient mice. J Exp Med. 1999, 189: 1129-1138. 10.1084/jem.189.7.1129.
Smeltz RB, Chen J, Hu-Li 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 IL-4. J Exp Med. 2001, 194: 143-153. 10.1084/jem.194.2.143.
Moore KW, de Waal Malefyt R, Coffman R, O'Garra A: Interleukin-10 and the interleukin-10 receptor. Annu Rev Immunol. 2001, 19: 683-765. 10.1146/annurev.immunol.19.1.683.
Riley JK, Takeda K, Akira S, Schreiber RD: Interleukin-10 receptor signaling through the JAK-STAT pathway. J Biol Chem. 1999, 274: 16513-16521. 10.1074/jbc.274.23.16513.
Schmidt-Weber CB, Alexander SI, Henault LE, James L, Lichtman AH: IL-4 enhances IL-10 gene expression in murine Th2 cells in the absence of TCR engagement. J Immunol. 1999, 162: 238-244.
Skapenko A, Niedobitek GU, Kalden JR, Lipsky PE, Schulze-Koops H: Generation and regulation of human Th1-biased immune response in vivo: A critical role for IL-4 and IL-10. J Immunol. 2004, 172: 6427-6434.
Huang Y, Wange RL: T cell receptor signaling: beyond complex complexes. J Biol Chem. 2004, 279: 28827-28830. 10.1074/jbc.R400012200.
Porter CM, Clipstone NA: Sustained NFAT signaling promotes a Th1-like pattern of gene expression in primary murine CD4+ T cells. J Immunol. 2002, 168: 4936-4945.
Kauffman SA: Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol. 1969, 22: 437-467. 10.1016/0022-5193(69)90015-0.
Kauffman SA: Antichaos and adaptation. Sci Am. 1991, 265: 78-84.
Thomas R: Laws for the dynamics of regulatory networks. Int J Dev Biol. 1998, 42: 479-485.
Muraille E, Leo O: Revisiting the Th1/Th2 paradigm. Scand J Immunol. 1998, 47: 1-9.
Street NE, Mosmann TM: Functional diversity of T lymphocytes due to secretion of different cytokine patterns. FASEB J. 1991, 5: 171-177.
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: 451-494. 10.1146/annurev.immunol.18.1.451.
We want to thank Massimo de Francesco, Mark Ibberson, Caroline Johnson-Leger, Maria Karmirantzou, Lukasz Salwinski, François Talabot and Francisca Zanoguera for their valuable comments and suggestions.
The author(s) declare that they have no competing interests.
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.
Electronic supplementary material
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.