Transforming Boolean models to continuous models: methodology and application to T-cell receptor signaling

doi:10.1186/1752-0509-3-98

BMC Syst Biol. 2009; 3: 98.

Published online 2009 September 28. doi: 10.1186/1752-0509-3-98

PMCID: PMC2764636

Transforming Boolean models to continuous models: methodology and application to T-cell receptor signaling

Dominik M Wittmann,^1,² Jan Krumsiek,¹ Julio Saez-Rodriguez,^3,⁴ Douglas A Lauffenburger,³ Steffen Klamt,⁵ and Fabian J Theis^1,^2,⁶

¹Institute for Bioinformatics and Systems Biology, Helmholtz Zentrum München - German Research Center for Environmental Health, 85764 Neuherberg, Germany

²Department of Mathematical Science, Technische Universität München, 85747 Garching, Germany

³Biological Engineering Department, M.I.T., Cambridge MA 02139, USA

⁴Department of Systems Biology, Harvard Medical School, Boston MA 02115, USA

⁵Max Planck Institute for Dynamics of Complex Technical Systems, 39106 Magdeburg, Germany

⁶Max Planck Institute for Dynamics and Self-Organisation, 37073 Göttingen, Germany

Corresponding author.

Dominik M Wittmann: dominik.wittmann/at/helmholtz-muenchen.de; Jan Krumsiek: jan.krumsiek/at/helmholtz-muenchen.de; Julio Saez-Rodriguez: julio/at/hms.harvard.edu; Douglas A Lauffenburger: lauffen/at/mit.edu; Steffen Klamt: klamt/at/mpi-magdeburg.mpg.de; Fabian J Theis: fabian.theis/at/helmholtz-muenchen.de

Received January 19, 2009; Accepted September 28, 2009.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article has been cited by other articles in PMC.

Abstract

Background

The understanding of regulatory and signaling networks has long been a core objective in Systems Biology. Knowledge about these networks is mainly of qualitative nature, which allows the construction of Boolean models, where the state of a component is either 'off' or 'on'. While often able to capture the essential behavior of a network, these models can never reproduce detailed time courses of concentration levels.

Nowadays however, experiments yield more and more quantitative data. An obvious question therefore is how qualitative models can be used to explain and predict the outcome of these experiments.

Results

In this contribution we present a canonical way of transforming Boolean into continuous models, where the use of multivariate polynomial interpolation allows transformation of logic operations into a system of ordinary differential equations (ODE). The method is standardized and can readily be applied to large networks. Other, more limited approaches to this task are briefly reviewed and compared. Moreover, we discuss and generalize existing theoretical results on the relation between Boolean and continuous models. As a test case a logical model is transformed into an extensive continuous ODE model describing the activation of T-cells. We discuss how parameters for this model can be determined such that quantitative experimental results are explained and predicted, including time-courses for multiple ligand concentrations and binding affinities of different ligands. This shows that from the continuous model we may obtain biological insights not evident from the discrete one.

Conclusion

The presented approach will facilitate the interaction between modeling and experiments. Moreover, it provides a straightforward way to apply quantitative analysis methods to qualitatively described systems.

Background

Close interaction between experiments and mathematical models has proven to be a powerful research approach in Systems Biology. Especially the modeling of regulatory and signaling networks, however, is typically hampered by a lack of information about mechanistic details, as often one can only determine the interactions of the involved species in a qualitative way. The current shift of focus in Systems Biology from single signal transduction pathways to networks of pathways exacerbates this lack of information even more. Therefore, the creation of mass action based models that accurately describe the underlying biochemistry is typically restricted to small well-studied subsystems.

Large-scale models of regulatory or signaling networks are often so-called Boolean models [1]. In fact, these models can be seen as the mathematically rigorous representation of qualitative biological knowledge. Their components, henceforth called species, can have only discrete states, typically two; these may be referred to as 0 and 1, 'off' and 'on', 'deactivated' and 'activated', etc. Time is discretized and the state of a species at time t + 1 is a function of the states of the species at time t. Although being a crude simplification of biological reality, Boolean models are often able to reproduce the qualitative behavior of a system [2-8]. Naturally, Boolean models can neither describe continuous concentration levels nor realistic time scales. For this reason, they cannot be used to explain and predict the outcome of biological experiments that yield quantitative data. However, with increasing emphasis on these quantitative experiments the need for precisely this kind of model arises. In this contribution, we present and exemplify a practicable solution to this problem: a standardized method for accurately converting any Boolean model into a continuous model. This transformation fills the gap in the modeling process shown in Figure Figure1.1. It allows construction of a continuous model from qualitative knowledge by representing this knowledge as a Boolean model and then transforming this discrete model into a continuous one. The continuous model can now be used to explain experimental results and to design and optimize further experiments. The results of these experiments, in turn, help to refine the model.

Figure 1

Overview of the modeling process. The typically qualitative biological knowledge is mathematically rigorously represented as a Boolean model, that is then converted into a continuous model. This continuous model can be used to explain quantitative experimental (more ...)

Boolean models are a very coarse description of biochemical processes. They phenomenologically describe observed dependencies often leaving out still unknown players or intermediate steps. As our transformation requires no additional information the resulting continuous models are, of course, still phenomenological models. We can automatically create these continuous phenomenological models out of a Boolean model, but we cannot create a mass action law without additional knowledge on the biochemistry. Our method is a top-down approach for (large) networks with incomplete mechanistic knowledge — derived e.g. from pathway databases — where predictive kinetic modeling is infeasible. The main point that we want to make in this contribution is that also these phenomenological models can be used in Systems Biology to explain and predict quantitative experimental results.

To this end, we focus on a large Boolean model of T-cell activation proposed by Klamt et al. [6]. These cells play a pivotal role in the immune system. When foreign antigens bind to their receptors, signaling cascades are triggered within the T-cell leading to an activation of several transcription factors. The logical structure of the model is shown in Figure Figure2A.2A. There are three inputs: the T-cell receptor TCR, the coreceptor CD4 and an input for CD45; as well as four outputs: the transcription factors CRE, AP1, NFkB and NFAT. The rephosphorylation of PAG-Csk by Fyn and cCbl mediated degradation are known to be slow processes compared to the other interactions. This is modeled by activating the feedback loops Fyn → PAG-Csk and ZAP-70 → cCbl only at a later stage. Therefore three scenarios are defined:

Figure 2

T-cell model. (A) Structure of the Boolean model as shown in CellNetAnalyzer [24]. (B) Hill functions with parameter n = 3 and different thresholds k = 0.3, k_fast= 0.1 and k_slow= 0.8. (C) Subnet of the T-cell model for scenario 2 (only activation) and (more ...)

scenario 1, resting state: All inputs are set to 0, feedback loops are deactivated.

scenario 2, early events: Inputs are set to 1, feedback loops are still deactivated.

scenario 3, mid-time events: Inputs are still 1, feedback loops are active.

A qualitative analysis of an expanded version of this Boolean model by Saez-Rodriguez et al. [8] yielded new and non-obvious signaling pathways. There are also quantitative models covering aspects of T-cell signaling in mechanistic detail, e.g. [9]. In contrast thereto, the continuous model we obtain in the following describes the T-cell signaling cascade on a larger scale yet in less mechanistic detail.

We apply the transformation method to the T-cell model. The resulting continuous model is able to fully reproduce the behavior of the Boolean model. Moreover, we show that it can easily include and deal with the different time scales of interactions in the three scenarios. Hence the continuous model is indeed a generalization of the Boolean model with richer dynamic properties. This is in line with previous findings [10,11], indicating that the qualitative behavior of a discrete model is reproduced by its continuous homologue. We can further corroborate this hypothesis by generalizing existing theoretical results on the steady-states of discrete and continuous models.

The crucial question is, of course, whether the continuous model derived from a Boolean network is indeed competent to explain an aspect of biological reality in a precise quantitative fashion. We answer this question in the affirmative by showing that our T-cell model reproduces time courses of concentration levels measured for three different ligand concentrations [12]. Moreover, it is able to predict binding affinities of different ligands from their induced signaling profiles. The fact that the model can differentiate between more than two different concentrations shows that we have definitely left the Boolean (binary) world behind.

Results

Representation of Boolean functions and models

A Boolean model consists of

• N species X₁, X₂, ...., X_N, e.g. genes, proteins, etc., each represented by a variable x_itaking values in {0, 1},

• for each species X_ia set of species An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i1.gif

that influence x_iand

• for each species X_ian update function

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i2.gif

giving the value of x_iat the next time step for every possible combination of An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i3.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i3.gif

In the following we think of B_ias a function on the vertices of the unit cube. Since vertex An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i4.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i4.gif

corresponds to the AND gate

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i5.gif

we can write B_iin the form

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i6.gif

This is a so-called sum-of-product representation of B_i. These representations are especially convenient, as they allow to graphically represent our models using interaction hypergraphs [6]. The idea is to represent each product (AND gate) in the sum-of-product form of B_iby a directed hyperedge between a set of start nodes An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i7.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i7.gif

and the end node X_i. Each pair (s, X_i), s [set membership]

S, carries a sign — '+' or '-' — depending on whether there is a factor s or ¬s in the product. All incoming hyperedges at node X_iare then a graphical representation of B_i. This is further illustrated in Additional data file 1.

The general approach to making discrete models continuous

The first step for obtaining a continuous model from a Boolean one is to replace the discrete variables x_iby continuous variables An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i8.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i8.gif

taking values in [0, 1], i.e. we normalize concentrations to the unit interval. Consequently, we have to 'extend' the functions B_ito functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

: [0, 1]^N→ [0, 1]. We call the functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

continuous homologues of the Boolean functions B_i. The crucial point is, how the transformation of a Boolean function into a continuous homologue is performed. We will address this issue in the next section and continue with the outline of the general approach.

In the second step we have to specify how to build the actual continuous model, for which there are two possibilities. The most straightforward is probably to proceed analogously to the Boolean model, i.e. to use discrete time steps and to compute the value of An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i8.gif

at time t + 1 by

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i10.gif

(1)

In numeric simulations the discretization of time is obviously irrelevant. It complicates, however, the detection of small-scale continuous effects and is a serious drawback in the further investigation of the model by analytical methods.

Another way to build the continuous model is to try to mimic biological reality more closely: mRNAs, proteins, etc. are produced at a certain rate and are at the same time degraded. We assume the production of X_ito be given by An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

, and the degradation to be proportional to An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i8.gif

. Then the development of An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i8.gif

over time is governed by the ordinary differential equation (ODE)

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i11.gif

(2)

where τ_ican be interpreted as the life-time of species X_i. Note that due to the normalization of concentrations to the unit interval we have only one parameter for production as well as decay. To clarify this, assume that An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i12.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i12.gif

denotes the non-normalized concentrations and An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i13.gif

the corresponding production functions. Then a general system of ODEs is of the form An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i14.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i14.gif

, where α_iis the production rate and γ_ithe decay rate of species X_i. Since the maximal concentration of X_iis α_i/γ_iwe have the relation An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i15.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i15.gif

. It now follows that

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i16.gif

and by setting τ_i= 1/γ_iand An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i17.gif

we obtain the ODEs (2).

Herein we focus on model (2), as this model can be further analyzed using the rich and mathematically rigorous theory of ODEs. Note furthermore that model (1) can be considered a special case of model (2) after numeric integration.

Continuous homologues of Boolean functions

As already mentioned, the key point is how a continuous homologue An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

can be obtained from a Boolean function B_iin a computationally efficient manner. A suitable transformation has to satisfy three conditions:

Accuracy: It has to be accurate, which means that An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

and B_imust agree on the vertices of the unit cube, i.e. for values in An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i18.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i18.gif

. Functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

satisfying this condition are called perfect continuous homologues. As will be shown, for perfect continuous homologues discrete and continuous models exhibit a similar steady-state behavior.

Good analytical properties: The functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

should have good analytical properties such as smoothness in order to allow and facilitate a mathematical analysis of the system of ODEs.

Minimality and uniqueness: The functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

should be the unique minimal solution in their interpolation class.

The three transformations we propose in the following are all based on multivariate polynomial interpolation [13,14] (see Methods). Here An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

is defined as a polynomial in the variables An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i19.gif

that agrees with B_ion the vertices of the unit cube. As will be shown in the Methods section, this technique satisfies all three of the above requirements. There are other approaches which we shortly review and compare in the Discussion section.

BooleCube

In a first step, we define the functions

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i20.gif

(3)

by linearly interpolating the functions B_iusing the technique of multivariate polynomial interpolation as explained in the Methods section. These functions are called BooleCubes. By substituting An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i21.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i21.gif

for B_iin equation (2), we can then define a system of ODEs that describes the temporal development of the An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i8.gif

HillCube

The functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i21.gif

are affine multilinear, i.e. for each 1 ≤ j ≤ N_iand fixed An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i22.gif

, k ≠ j, there exist constants a, b [set membership]

such that

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i23.gif

. Molecular interactions, however, are known to show a switch-like behavior, which can be modeled using sigmoid shaped Hill functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i24.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i24.gif

[15] (see Additional data file 2). The two parameters n and k have a clear biological meaning. The Hill coefficient n determines the slope of the curve and is a measure of the cooperativity of the interaction. The parameter k corresponds to the threshold in the Boolean model, above which one defines the state of a species as 'on'. Mathematically speaking, it is the value at which the activation is half maximal.

We now define a Hill function f_ijwith parameters n_ijand k_ijfor every interaction and define new functions

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i25.gif

(4)

which we call HillCubes. Plots of the HillCubes of all 16 two-variable Boolean gates can be found in Additional data file 3. Now a new system of ODEs can be defined by replacing the An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

in equation (2) by the HillCubes An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i26.gif

normalized HillCube

Note that Hill functions never assume the value 1, but approach it asymptotically. Hence, the HillCubes are not perfect homologues of the Boolean update functions B_i. If this is desired a simple solution is to normalize the Hill functions to the unit interval. This yields another (perfect) continuous homologue of the

Boolean functions B_i

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i27.gif

(5)

which we call normalized HillCubes, and thus defines also a new continuous ODE model.

Theoretical results about the relation between discrete and continuous models

A natural and interesting question is how similar the discrete and the continuous model are. It can easily be shown (see Methods and [11]) that, whenever the An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

in equation (2) are perfect homologues of the Boolean update functions B_i, a steady-state of the Boolean model will also be a steady-state of the continuous system. This result can be applied to the BooleCube model as well as the normalized HillCube model, but not to the non-normalized HillCube model. Therefore, Boolean steady-states will in general not be steady-states of the non-normalized HillCube system. The question is if Boolean steady-states still correspond to 'similar' steady-states of this continuous model, at least for certain parameters. Using the implicit function theorem we were able to show that this is indeed the case (see Methods). The reverse statement is, of course, not true, as in the continuous model additional (stable) steady-states may arise. Besides the steady-state behavior, monotony properties are a further important characteristic of a dynamical system. These properties determine the effect of a down- or up-regulation of a certain species on the other species. Due to its accuracy the presented transformation method preserves the monotony properties of a Boolean network. In the Discussion section we illustrate this further using the T-cell model as an example.

We show these results as they justify our transformation approach by demonstrating that it preserves essential properties of the Boolean model. For a deeper mathematical investigation of the relation between discrete and continuous models we refer the interested reader to the rich literature on this subject, e.g. [11,16,17]

Simulation of the Boolean T-cell model

Figure Figure3A3A shows a simulation of the Boolean T-cell model using synchronous updates. In the beginning the three inputs were 0 and the system was in its resting state (scenario 1). At t = 10 the inputs were manually switched on and the signaling cascade was triggered off, which at t = 24 led to the activation of all four transcription factors CRE, AP1, NFkB and NFAT (scenario 2). At t = 57 we activated the feedback loops Fyn → PAG-Csk and ZAP-70 → cCbl. Consequently the cascade was blocked and at t = 67 all transcription factors were again deactivated. So, essentially the simulation used three different Boolean models at times [0, 10), [10, 57) and [57, 100]. Note that in the Boolean model the time point for the activation of the feedback loops could be chosen arbitrarily. As is explained below, in the continuous model this time point was determined by our choice of the kinetic parameters. We chose t = 57 since then deactivation of all transcription factors occurred at around the same time point in both models.

Figure 3

Results from the Boolean and continuous simulations of the T-cell model using the manually determined parameter set. The vertical dashed lines in (A, B, F) mark (from left to right): the switching on of the inputs, the total activation of the transcription (more ...)

After deactivation oscillations occurred, artefacts of the synchronous updating. When the species were updated asynchronously according to some permutation, these oscillations could be observed for 3078 out of 10000 randomly sampled permutations; in the other cases a steady-state was reached.

Continuous T-cell model

We applied the transformation technique to the Boolean T-cell model. Using non-normalized HillCubes we obtained a large quantitative model of T-cell activation (see Methods). In a first step, we manually determined approximate parameters for which the continuous model reproduces the behavior of the Boolean model during all three scenarios. The information about fast and slow interactions was thereby encoded in the values of specific parameters. Consequently, the continuous model was able to explain both, the activation as well as the deactivation of T-cells, without any alteration of the network topology between different scenarios. Note that this was necessary in the Boolean case. We set all Hill coefficients to n = 3, the thresholds to k = 0.3 and the life-times to τ = 1, with the exception of Fyn and ZAP-70. Here we knew that the interactions Fyn → PAG-Csk and ZAP-70 → cCbl operate on a slower time scale. Therefore the thresholds for these interactions were set to k_slow= 0.8 whereas the thresholds for the interactions Fyn → TCR-phos, ZAP-70 → LAT-phosp, ZAP-70 → PLCg (act) and ZAP-70 → Itk were set to k_fast= 0.1. The Hill functions for the three thresholds k, k_slowand k_fastare displayed in Figure Figure2B.2B. Finally, we set the life-times of Fyn and ZAP-70 to τ_Fyn= τ_ZAP-70= 10 to enlarge the time gap between the two switching points.

The numeric simulation of the continuous T-cell model with the manually determined parameters is shown in Figure Figure3B.3B. At first, the cell was again in its resting state with all three inputs turned off. Then at t = 10 we manually switched on all inputs and the signaling cascade was triggered off showing an expression profile very similar to the Boolean simulation (Figure (Figure3A).3A). We observed a total activation of the four transcription factors at around t = 21 just like in the Boolean simulation. One can clearly see that Fyn and ZAP-70 were activated more slowly. Nonetheless TCR-phos, LAT-phosp, PLCg (act) and Itk were instantly turned on due to the low threshold k_fast. Only when at around t = 25 the concentrations of Fyn and ZAP-70 were high enough, the feedback loops Fyn → PAG-Csk and ZAP-70 → cCbl became active and began to switch on PAG-Csk and cCbl. While PAG-Csk reached a constant medium expression level, cCbl was only weakly and transiently expressed. This, however, sufficed to switch off the cascade and ultimately at around t = 67 all transcription factors were again deactivated. In contrast to the Boolean model, the continuous model did not exhibit an oscillatory behavior but reached a steady-state after deactivation of the signaling cascade. The species with long expression periods in the Boolean oscillation (Tcr bind, PAG Csk, Fyn, TCR phos and IkB) were expressed at high or medium levels whereas the species with short expression periods were not expressed at all.

Comparison of the discrete and the continuous T-cell model

Although one can argue that there is no real time scale in Boolean networks, we compared both models by substracting the discrete from the continuous time course (Figure (Figure3C).3C). The activation of the feedback loops in the Boolean model was conveniently chosen such that we have a total deactivation of the transcription factors at around t = 67 in both models. While we observed an almost perfect agreement of the time courses of the three inputs, there was a huge difference in the time courses of the species involved in the two feedback loops. This was not surprising, considering that these loops are regulated differently in both models. The species downstream of the regulatory loops (from LAT phosp to CRE in Figure Figure3C)3C) showed again a similar expression pattern.

We then analyzed this part of the cascade more deeply. Figure Figure3C3C also shows the correlation between the discrete and the continuous expression pattern of the species downstream of the feedback loops at each time point 0 <t < 100. We observed a high correlation in the more stationary phases (resting state, activated state and deactivated state) and a significant drop of correlation during the transitions between these phases. This met our expectation that the two models show the same qualitative but a different dynamic behavior.

Ratio of activation and deactivation time

When looking at Figures Figures3A3A and and3B,3B, a striking difference between the dynamics of both models is that in the discrete model activation and deactivation took approximately the same time, whereas in the continuous model activation was a much faster process than deactivation. This can also be seen from the red and blue 'steps' in upper Figure Figure3C3C during the activation and deactivation phases. To confirm that this was not merely an artefact of our choice of parameters, we calculated and analyzed the ratio between activation and deactivation time for different parameters. We defined

• the beginning of activation t_ActBegas the time when the first species (of the lower 29) reaches 5% of its maximal value,

• the end of activation t_ActEndas the time when the last species reaches 95% of its maximal value,

• the beginning of deactivation t_DeactBegas the time after t_ActEndwhen the first species drops under 95% of its maximal value,

• the end of deactivation t_DeactEndas the time when the last species drops under 5% of its maximal value,

• and, finally, the ratio of interest

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i28.gif

In the Boolean model we could easily compute ρ_b= 1 implying equally fast activation and deactivation. In the continuous model the crucial parameters were the life-times τ_Fynand τ_ZAP-70on the one hand and the concentration thresholds k_fastand k_slowon the other hand. The remaining parameters were set to n = 3, k = 0.3 and τ = 1.

First, we computed ρ for fixed k_fastand k_slowand different τ_Fyn= τ_ZAP-70. The result is shown in Figure Figure3D.3D. For τ_Fyn= τ_ZAP-70= 1 the cascade was not activated properly. For larger values we observed a decrease in ρ implying that an increase of τ_Fynand τ_ZAP-70prolonged the deactivation phase. This was to be expected — longer life-times resulted in a lessened increase of the decisive elements Fyn and ZAP-70 in the regulatory loop.

Second, we analyzed the effect of k_fastand k_slowfor fixed τ_Fyn= τ_ZAP-70= 10. The result is shown in Figure Figure3E.3E. Only for parameters k_fast [double less-than sign]

k_slow< 1 the cascade was activated properly. This agrees well with the fact that the difference between these two parameters is responsible for the delayed activation of the feedback loops. If it was not big enough the cascade was being deactivated before it had been fully activated. The greater this difference was, i.e. the farther we go away from the diagonal in Figure Figure3E,3E, the smaller ρ got, implying that a later activation of the feedback loops prolonged the deactivation phase. However, despite all these influences of the parameters, we observed much smaller ratios ρ in the continuous model than in the Boolean model. The average ρ's in Figures Figures3D3D and and3E3E were An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i29.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i29.gif

= 0.24 and An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i30.gif

= 0.27, respectively, which were both significantly smaller then ρ_b= 1. This suggests that ρ [double less-than sign]

1 may be an invariant of the dynamical system of biological importance.

Explanation of experimental data using the continuous T-cell model

The crucial question is if our model can reproduce real data. To test this we used the data set presented in [12] (see Methods). It describes the dynamics of the activation of key signaling elements upon activation of the TCR by three different ligands with varying affinity, Q144, Y144 and L144. We considered, in particular, the ligand L144 for which experiments were performed with three different ligand concentrations. Using global minimization of the model fit with respect to the data, we were able to determine a set of parameters for which our model reasonably reproduced the experimental data (see Methods and Table Table1).1). Figures 4A-D show the corresponding simulated time courses. We see that the model was able to approximate the time courses of ERK and IKK well. The fit of JNK was also acceptable, although the high measured concentrations at time t = 0 constituted a problem, as the model was naturally unable to reproduce them due to the delaying effect of the signaling cascade. In the case of NFAT, the model was unable to reproduce the non-monotone dependence on the ligand concentration. This suggests that the network structure cannot be reconciled with the non-monotone dependence of the two key signaling molecules JNK and especially NFAT with respect to the ligand concentration. A non-monotone response of T-cell signals has been reported in other contexts [18], and is consistent with the role of T-cells: their response has to be exquisitely regulated so as to reply only to a particular stimulus; an uncoupled response between the JNK (and p38) and ERK MAP Kinases has also been observed [19]. The regulation of NFAT and JNK, and more generally of TCR-induced signaling, is complex and not yet fully understood. To mention just two examples, processes such as Ras localization [20] have been shown to play a key role, but are not included in the model, and the regulation of calcium, which governs NFAT behavior, is more complex than described in the model. It is out of the scope of this paper to investigate this intriguing behavior in detail, but this result illustrates the power of our approach to gain new biological insight: taking exactly the same knowledge encoded in the discrete model and fitting it to quantitative data, we were able to identify the incompleteness of our model in an aspect that we could not have explored with a discrete model.

Table 1

Best fit parameter set for ligand L144.

Figure 4

Results of the parameter fit. (A-D) Simulation of the continuous model (solid lines) for high (red), medium (blue) and low (green) concentrations of ligand L144 and experimentally measured concentrations, cf. markers '+', '×' and '○'. (more ...)

Best fit parameter set

Subsequently, we analyzed the distribution of the Hill parameters within the best fit parameter set for ligand L144 (Table (Table1).1). When looking at the distribution of the exponents n upstream of the measured species (Figure (Figure4F)4F) we found that after the optimization 36 out of 49 (73%) were below the ad hoc estimate (3) and only 8(16%) were above it. This was to be expected, as in order to fit three different concentration levels the model had to contain mainly slow switches (low n) and only few Boolean-like switches (n → ∞). Figure Figure4E4E shows the distribution of the threshold parameters k upstream of the measured species. As explained above, we had manually set these parameters to 0.3, with the exception of four k's which had been set to k_slow= 0.8 and two k's which had been set to k_fast= 0.1. Interestingly, this structure had been preserved during the optimization. The mean of the thresholds k was 0.298 and hence well agreed with the ad hoc estimate. Also the high and low thresholds were still at least two standard deviations away from the mean, cf. the red and blue markers '+' in Figure Figure4E.4E. Only one other parameter k also had a Z-score above 2: the threshold k_L144for the stimulation of the TCR by the ligand L144, cf. marker 'o'. Possible implications hereof are discussed below.

Due to the large number of parameters the obtained parameter set was, of course, far from being unique. But we showed that a continuous model inferred from a Boolean model is able to reproduce experimental data in a quantitative way. Moreover, this transformation could enhance the explanatory power of the model in the sense that it was enabled to differentiate between more than two states.

In our example, the threshold parameters k are rather tightly centered around their mean. In principle, however, we could also have extreme outliers, i.e. very large or very small (≈ 0) values in the distribution. Mapped back to the Boolean model, this would imply a change of the network topology, as the corresponding reactions are then quasi-constant, either 'off' (for very large k) or 'on' (for very small k). Thus, a fitting of the continuous model to experimental data may also yield information about the network structure of the Boolean model.

Prediction of binding affinities of different ligands

As already mentioned, the fitted threshold k_L144for the stimulation of the TCR by the ligand L144 was significantly below the mean of the other thresholds k, cf. marker 'o' in Figure Figure4E.4E. This gave rise to the question of which affinities the model predicts for the other two ligands. To this end, we fitted the affinities of Q144 and Y144, mapped to the inverse of the Hill function's threshold parameters k_Q144and k_Y144, respectively, keeping the rest of the parameters constant. As expected, the fits themselves were far from perfect, due to parameter indeterminacies. Surprisingly however, the values of k_Q144and k_Y144we obtained were significantly above the mean of the other thresholds k, cf. markers ' [diamond with plus]

' and '□' in Figure Figure4E.4E. This suggests that the predicted relation between the parameters k_L144 [double less-than sign]

k_Y144<k_Q144is not simply an artefact of the optimization process. And indeed, it agrees well with experimental data [21].

Discussion

We now further discuss the presented transformation method and compare it to various other approaches. Also the relation between discrete and continuous models is discussed, especially with respect to their steady-state behavior and monotony properties.

Comparison of different transformation approaches

The relation between discrete and continuous models has already been investigated and various approaches to the problem of constructing the continuous homologues of the Boolean update rules have been proposed. In the following we shortly review previous work and compare the different approaches.

Piecewise linear differential equations

The idea to compare continuous and discrete models is almost as old as Boolean modeling itself. In 1973, Glass et al. [11] studied the relation between discrete models and ODE models of the form (2). Their motivation, however, was quite the opposite of ours. While we intend to enrich the dynamic behavior of discrete models, Glass et al. wanted to investigate the qualitative properties of continuous networks by studying corresponding simpler discrete models. They propose Hill functions as a suitable continuous homologue of one-variable Boolean step functions. In the case of multi-variable Boolean functions B_ia (perfect) continuous homologue An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

is constructed as follows: Note that, when building a Boolean model, one implicitly introduces a threshold 0 <θ_i< 1 for each species X_iand defines its state as 'on' if its concentration is above this threshold. The hyperplanes An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i8.gif

= θ_i, i = 1, 2, ..., N decompose the cubes An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i31.gif

into

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i32.gif

rectangular regions called domains. Each of these domains contains exactly one vertex An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i33.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i33.gif

and is denoted by An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i34.gif

accordingly (see Additional data file 4). A simple way of defining the functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

is now to set

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i35.gif

(6)

With this definition model (2) is a so-called piecewise linear ODE model which means that within each domain An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i34.gif

equation (2) is a linear ODE of the form An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i36.gif

. This kind of equation is very well understood and can be solved analytically. The functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i37.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i37.gif

defined in (6) perfectly agree with B_ion the vertices of the unit cube. Using piecewise linear ODEs Glass et al. could prove some theoretical results on the relation between discrete and continuous models, e.g. that Boolean steady-states are also steady-states of the continuous model. Some of these results are restated and generalized in the Methods section. Piecewise linear models are typically not used for quantitative simulations, as the step-like transition between the different domains is often unrealistic. Rather they are analyzed in a qualitative and semi-qualitative way, where their trajectories between the different domains are treated analogously to the state transition graphs of Boolean models [22].

Fuzzy logic

Another well studied way of generalizing Boolean models is fuzzy logic [23]. Recall that in a Boolean model one defines the state of a species as 'on' if its concentration is above a certain threshold. In fuzzy logic this concept is relaxed and a so-called degree of membership (DOM) function An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i38.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i38.gif

is introduced for each species X_i. For concentrations 0 = An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i8.gif

= 1 this function gives the degree with which we say that X_iis 'on'. There are two standard ways of generalizing the Boolean operators AND, OR, NOT:

(i) min-max logic

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i39.gif

(7)

(ii) product-sum logic

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i40.gif

(8)

When using product-sum logic, a normalization to the unit interval is necessary, since An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i41.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i41.gif

can assume values greater than 1. Both, (i) as well as (ii) are ways to construct the functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

from the Boolean functions B_i. If the DOM functions satisfy μ(0) = 0 and μ(1) = 1, the functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i42.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i42.gif

obtained by min-max logic agree with B_ion the vertices of the unit cube and hence the steady-states of the Boolean model are also steady-states of the continuous model. The major drawback of this method is that the functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i42.gif

are in general not differentiable, and do not have a 'nice' analytic representation (Figure (Figure5B).5B). Hence most analysis methods are not applicable to the resulting ODE systems. When we use product-sum logic, we encounter the problem that the resulting functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i43.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i43.gif

do not necessarily agree with B_i(Figure (Figure5C5C).

Figure 5

Comparison of different transformation techniques. Different continuous homologues of a Boolean OR gate. (A) Piecewise linear function An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i44.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i44.gif

. (B) Function An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i45.gif

obtained by min-max fuzzy logic and linear DOM functions. (C) Function An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i46.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i46.gif

obtained by product-sum fuzzy (more ...)

Standardized qualitative dynamical systems

Mendoza et al. [10] put forward a method to transform a Boolean model into a system of ODEs similar to (2) called standardized qualitative dynamical system. Their approach, however, is applicable only to a subclass of Boolean models: For each species X_iwe have a set of activators A_iand a set of inhibitors I_iof X_i. Then x_iis set to 1 at the next time step if at this time any of its activators and none of its inhibitors are acting upon it, otherwise x_iis set to 0. This corresponds to the Boolean logic An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i50.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i50.gif

. Clearly, many logic functions, like XOR gates for example, cannot be represented that way. From this Boolean model a continuous model is built up consisting of ODEs of the form

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i51.gif

Where

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i52.gif

(9)

The right-hand side of the above ODE consists of two parts: an activation function and a term for decay as in (2). The activation is given by a sigmoid shaped function of ω_i, where ω_irepresents the total input to node X_i. The steepness of the activation function is determined by the parameter h. Decay is assumed to be proportional to An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i8.gif

. Actually, a more general form of ω_iis introduced in [10], where the influence of the activators and inhibitors can be differently weighted. For the sake of better comparability we set all these weights equally to 1, as suggested in [10]. The more activators and the less inhibitors of a node are 1, the greater ω_iis. It takes its minimum (0), iff all activators are 0 or all inhibitors are 1, and its maximum (1), iff all activators are 1 and all inhibitors are 0. This, however, does not exactly correspond to the Boolean rule from above. In consequence, steady-states of the Boolean model are not necessarily steady-states of the continuous model. This will be further discussed below.

Multivariate polynomial interpolation

The aforementioned attempts only lead to functions which are either not differentiable or do not precisely generalize the Boolean logic. A straightforward approach to eliminate these drawbacks is to use multivariate polynomial interpolation [13] for the construction of the functions An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

(see Methods). This technique can be applied to any Boolean function B_i. The resulting BooleCubes An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i21.gif

are smooth and can easily be analytically differentiated and integrated. They agree with the Boolean functions B_ion the vertices of the unit cube and hence the steady-states from the Boolean model are also steady-states in the continuous model.

We define Hill functions f_ijfor all interactions and consider the HillCubes An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i26.gif

from equation (4). The idea behind this is that each interaction is described by its own Hill function with specific parameters and the different interactions are coupled by the BooleCubes. Since they are affine multilinear, the latter preserve the shape of the individual Hill functions. Thus, we can mimic single-component non-linearities which are common in switch-like regulatory systems. Additional data file 3 shows the HillCubes derived from all 16 two-variable Boolean gates. A mathematically rigorous treatment of this kind of dynamical systems can be found in [14]. If the Hill coefficient n goes to infinity the Hill function becomes more and more like a Boolean step function (see Additional data file 2). Hence for large exponents the HillCubes are very similar to the Boolean functions B_iand the continuous system will likely show an almost Boolean behavior. This is further illustrated in Figure Figure66.

Figure 6

Effect of increasing Hill exponents. We consider a simple cascade between the four species X₁, X₂, X₃, X₄as shown in the inset in (A). Each activation is modeled using a Hill function with threshold k = 0.5 and Hill coefficient n. The life-times τ (more ...)

The HillCubes do not perfectly agree with the Boolean update functions due to the asymptotic behavior of the Hill functions. By a suitable choice of the Hill parameters the difference can be reduced but not fully eliminated. An easy way to achieve a perfect agreement is to normalize the Hill functions to the unit interval, as is done in the normalized HillCubes An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i58.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i58.gif

from equation (5).

Comparison

To conclude, we illustrate the above methods applied to a simple OR gate between two species X₁and X₂. We compute

• the piecewise linear function

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i59.gif

from equation (6),

•

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i60.gif

obtained by fuzzy logic (with linear DOM functions) following (7),

•

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i61.gif

obtained by fuzzy logic (with linear DOM functions) following (8),

• the input function from equation (9) introduced by Mendoza et al.

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i62.gif

• the BooleCube An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i63.gif

from equation (3) obtained by the interpolation technique,

• the HillCube An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i64.gif

from equation (4) for Hill functions f₁and f₂with parameters n = 3, k = 0.5,

• and, finally, the normalization

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i65.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i48.gif

from equation (5).

Figures Figures5C5C and and5D5D show the product-sum fuzzy logic function An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i46.gif

and the input function ω. One can clearly see that they do not represent a pure OR gate, where the values at (x₁, x₂) = (1, 0) and (x₁, x₂) = (0, 1) should already be maximal. This is the case in Figures Figures5A5A and and5B5B which show the piecewise linear An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i44.gif

and the min-max fuzzy logic function An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i45.gif

. Here however, the problem is that the functions are not differentiable, as can easily be seen from their plots. The BooleCube An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i47.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i47.gif

shown in Figure Figure5E5E is both, smooth and maximal as soon as any concentration is equal to 1. Finally, Figures Figures5F5F and and5G5G show the (normalized) HillCubes An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i48.gif

and

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i49.gif

, respectively, which are also smooth and can be considered good transformations of the Boolean OR gate. An overview about the discussed advantages and disadvantages of the different transformation techniques is provided in Figure Figure5H5H.

Theoretical results about steady-states

A fundamental principle of biological modeling is that steady-states of a model typically correspond to the different operating modes or states of the biological system under study. This correspondence was also the motivation for Kauffman's seminal study [1], where Boolean models were introduced for the first time in biology. A critical step in the justification of any transformation method therefore is to ensure that at least the steady-states of the Boolean model are still steady-states in the homologue continuous system. In the case of a perfect agreement of the Boolean update rules B_iwith their continuous homologues An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

this can easily be shown (see Methods and [11]). This perfect agreement, however, is not a biologically plausible assumption; biological interactions, such as enzyme kinetics for example, are known to asymptotically approach but never fully reach a saturation level. Empiric evidence that also in real-world examples, the steady-states of a Boolean model correspond to steady-states of a homologue continuous model, is given by Mendoza et al. [10]. The method of transformation used therein has already been described above and we also mentioned that it does not accurately transform the Boolean update rule into a continuous activation function. This inaccuracy is due to a systematic difference between the Boolean logic and the analytic form of the activation function. It is not the result of an asymptotic sigmoid function; in fact, the used sigmoid function assumes its maximal values 0 and 1. One can easily construct an example where due to this systematic difference Boolean steady-states are not conserved under the transformation (see Additional data file 5).

In the case of the HillCube model, it is the other way round. There is no systematic difference between the Boolean update rules and the HillCube functions, the imperfect agreement is caused by the asymptotic behavior of the Hill functions. Therefore, the difference between both can be made arbitrarily small — albeit not zero — by a suitable choice of parameters. In this situation, we can show that for certain parameters, more precisely for sufficiently large exponents, there will be a steady-state of the continuous system in the neighborhood of each Boolean steady-state (see Methods). This theoretical result further justifies the presented transformation method.

Monotony properties

A nice feature of our method for converting Boolean into continuous models is that monotony properties, typically captured in the underlying interaction graph of the system, are preserved. Interaction graphs are signed directed graphs where each directed edge reflects a causal dependency, which can either be positive or negative, between its start and end node. Boolean models represented as interaction hypergraphs have a unique underlying interaction graph which can easily be derived from the logical model (by splitting the ANDs, see Klamt et al. [6]). For example, the Boolean function A = (¬B ∧ C) [logical or]

D would be translated into two positive arcs (C → A, D → A) and one negative arc (B [dashv R: dash, vertical]

A). In the interaction graph one may then compute the recently introduced dependency matrix [6,24], which determines for each pair (X, Y) of species the global effect of X on Y. This effect can — in some cases only initially — be positive, negative, ambivalent or vanishing.

For example, the dependency matrix of the T-cell model tells us, that LAT-phosp exerts purely positive effects on ERK, JNK, IKK, and NFAT, because there are only positive paths from LAT-phosp to these species and no negative feedback is involved. If we simulate a scenario with the logical model and repeat it then with e.g. fixing LAT-phosp to 1 (i.e. to the highest possible value), the resulting Boolean values in the four species mentioned above cannot decrease. In continuous systems, the interaction graph is encoded in the sign structure of the Jacobian matrix. In fact, in a continuous system obtained from a Boolean model the interaction graph is up to negative self-loops identical and monotony properties are therefore preserved. Accordingly, a positive perturbation in LAT-phosp, e.g. by permanently decreasing the threshold of ZAP-70 in the interaction activating LAT-phosp, results in a trajectory that is always above the trajectory of the non-perturbed system (Figure (Figure3F).3F). In fact, in accordance with the interaction graph of the Boolean model, we observe purely positive effects on all species downstream of LAT-phosp with the exception of ikB. Hence, important qualitative properties of the dynamics derived from the logical model are reflected in the dynamics of the continuous system.

Conclusion

With increasing amounts of quantitative data being available, the challenge arises how we can use our typically qualitative knowledge about biological systems to explain this data. For this purpose, we presented a canonical and fully standardized way of transforming qualitative discrete into continuous models. The transformation is accurate and we can show that it preserves the steady-state behavior as well as the monotony properties of the discrete model. The feasibility of the presented approach was substantiated by applying it to a logical model of T-cell receptor signaling. The resulting model is an extensive continuous model of T-cell activation. In contrast to the Boolean model it allowed to accommodate different time scales by adjusting kinetic parameters. It was competent to reproduce time courses of key signaling molecules measured for three different ligand concentration levels. Moreover, the model was able to predict the binding affinities of different ligands.

Being fully automatized [Krumsiek et al.: Odefy — From discrete to continuous models. In preparation (2009)] the presented method recommends itself to be applied to further biological systems. Future work could also aim at generalizing the approach from Boolean (binary) to s-state systems, where one no longer differentiates between two but s > 2 discrete states, e.g. 'low', 'medium' and 'high'. Finally, the relation between a discrete model and its continuous homologue needs to be further investigated, especially with respect to more complex behaviors like oscillations, which are of importance in many biological systems, such as cell cycle.

Methods

Multivariate polynomial interpolation

We now explain the technique of multivariate polynomial interpolation of a single Boolean function B_i. Therefore i [set membership]

{1, 2, ..., N} is fixed and for the sake of simplicity the subscript i is omitted. We remark that here B can be any real-valued function on the vertices of the unit cube {0, 1}^N, i.e. does not necessarily have to be a Boolean function. The idea is, to find a polynomial An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i47.gif

^N→

that is a continuation of B_iin the sense that An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i66.gif

for all

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i67.gif

One can easily see, that there is no unique solution to this problem. Therefore, we additionally require that the degree of the polynomial be minimal, where the degree of some polynomial

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i68.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i69.gif

for almost all (m₁, m₂, ..., m_N) [set membership]

^N, is defined as

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i70.gif

For

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i71.gif

we define

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i72.gif

We now show that An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i47.gif

indeed satisfies the three requirements for interpolation functions that we set out at the beginning (see section on continuous homologues of Boolean functions).

Theorem. The function An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i47.gif

has the following properties:

(i)

for all

(ii)

is the unique minimal degree polynomial interpolating B.

(iii) Let s denote the number of symmetry hyperplanes of B, i.e. the number of variables x_isatisfying

B(x₁, ..., x_i-1, 0, x_i+1, ..., x_N) = B(x₁, ..., x_i-1, 1, x_i+1, ..., x_N), for all (x₁, ..., An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i73.gif

, ..., x_N) ε {0, 1}^N-1. Then deg ( An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i74.gif

) = N - s.

(iv)

is affine multilinear. It is multilinear iff it corresponds to an AND gate, i.e.

B(x₁, x₂, ..., x_N) = An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i75.gif

(v) Let

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i76.gif

. Then it holds

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i77.gif

In particular, An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i78.gif

if B is a Boolean function.

Proof. (i) Note that for An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i8.gif

{0, 1} we have

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i79.gif

(ii) A minimal degree interpolation polynomial is of the form

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i80.gif

since exponents greater than 1 can be replaced by 1 without changing the values of f for An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i67.gif

. We order the vertices of the unit cube such that the number of 1's in the coordinates is not decreasing and denote the reordered sequence by An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i81.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i81.gif

. Then we consider the sequence of equations f(V_k) = B(V_k), k = 1, 2, ..., 2^N. For all k [set membership]

{1, 2, ..., 2^N} there is exactly one coefficient An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i82.gif

whose monomial An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i83.gif

satisfies g(V_k) = 0, k = 1, 2, ..., k - 1 and g(V_k) = 1. This allows to uniquely determine the coefficients An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i82.gif

and that way also the polynomial f. Since An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i47.gif

is of the form (*), it is the unique minimal degree interpolation polynomial.

(iii) The degree of An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i47.gif

is clearly deg ( An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i47.gif

) ≤ N. If B is symmetric with respect to the hyperplane x_i= 0.5, w.l.o.g. i = 1, we have B(0, x₂, ..., x_N) = B(1, x₂, ..., x_N) for all (x₂, ..., x_N) [set membership]

{0, 1}^N-1and consequently

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i84.gif

Hence, deg ( An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i47.gif

) ≤ N - 1. Inductively this proves (iii).

(iv) Let i [set membership]

{1, 2, ..., N}. Then at fixed An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i85.gif

the derivative An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i86.gif

is constant so An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i47.gif

is affine linear in An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i8.gif

. Moreover, we have

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i87.gif

For the last equivalence note that if An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i47.gif

is multilinear, then An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i88.gif

if any

= 0.

(v) Assume An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i89.gif

and w.l.o.g. An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i90.gif

. Then it follows from (iv) that An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i91.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i92.gif

. Inductively, we obtain x [set membership]

{0, 1}^Nsuch that An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i93.gif

, a contradiction to (i). Analogously, one proves that An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i94.gif

. □

Note that in general deg ( An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i47.gif

) = N - s does not hold, as can be seen from the Boolean function in three variables x₁, x₂, x₃

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i95.gif

This function does not have any symmetry hyperplanes, but its interpolation An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i96.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i96.gif

has degree 2.

Theoretical results about steady-states

The following theorem investigates the steady-state behavior of discrete and continuous models.

Theorem. Assume we are given a Boolean model and perfect continuous homologues An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

of the Boolean update functions. Then for any state vector An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i97.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i97.gif

the following are equivalent:

(i)

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i98.gif

is a steady-state of the Boolean model.

(ii)

is a steady-state of the model (1).

(iii)

is a steady-state of the ODE model (2).

Proof. The steady-state conditions are

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i99.gif

for the Boolean model,

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i100.gif

for model (1) and

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i101.gif

for model (2). Considering that the An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i9.gif

are perfect homologues of the B_i, these conditions are clearly equivalent. □

Remark. Note that the above theorem generally applies to transformations using perfect continuous homologues of Boolean functions. In the special case of piecewise linear ODEs, i.e. An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i102.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i102.gif

from equation (6), Glass et al. [11] could show that the above theorem also holds when (iii) is replaced by ' An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i98.gif

is a stable steady-state of the ODE model (2)'.

We now extend the above theorem to HillCube models. The problem is that we can no longer assume a perfect agreement between An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i26.gif

and B_i. The main idea is to use the implicit function theorem to prove the existence of a steady-state of the continuous model in a neighborhood of a Boolean steady-state.

Theorem. Assume we are given a Boolean model and the corresponding HillCube model. Let x^∞ [set membership]

{0, 1}^Nbe a steady-state of the Boolean model. We fix the thresholds k_ijat values in (0, 1) and the life-times τ_iat values in (0, ∞). Then there exists a neighborhood U of x^∞such that for all sufficiently large Hill exponents n = (n_ij) the HillCube model has a stable steady-state An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i103.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i103.gif

in U . It holds

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i104.gif

Proof. For a single Hill exponent n we write n = m^-2, m ≠ 0, fix some 0 <k < 1 and consider the Hill function

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i105.gif

We extend this function at m = 0 by setting

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i106.gif

Then there are open neighborhoods U₁, U₂of 0 and 1, respectively, such that f( An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i98.gif

, m) is continuous on U_1/2×

and it can be easily shown that f is even continuously differentiable on U_1/2×

. Now consider the HillCube model. It depends only on the Hill exponents n and we define m = (m_ij) by An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i107.gif

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i107.gif

. For concentrations An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i108.gif

let Φ(m,

) denote the right hand side of the HillCube ODE system. As explained above, we continuously extend the Hill functions and hence also Φ at m = 0. Then there exists an open neighborhood U [subset or is implied by]

^Nof x^∞such that Φ is continuously differentiable on

^|m|× U (as the composition of continuously differentiable functions).

For m = 0 the Hill functions become Boolean step functions and hence the HillCubes perfectly agree with the Boolean update rules. Therefore, we have Φ (0, x^∞) = 0. Now, let us compute the Jacobian DΦ of Φ in (0, x^∞).

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i109.gif

where δ_ijis the Kronecker delta. Note that for m = 0 the Hill functions in the HillCubes are step functions, i.e. constant in a neighborhood of x^∞. Hence the derivative of the HillCubes vanishes. This shows that DΦ is diagonal negative definite and hence, in particular, invertible. Therefore, the implicit function theorem guarantees that there are open neighborhoods U^'of x^∞and V' [subset or is implied by]

^|m|of 0 such that for each m [set membership]

V', i.e. for large exponents n, there is a An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i103.gif

U' such that Φ(m, An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i103.gif

) = 0, i.e. An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i103.gif

is a steady-state of the model. Since the mapping m [mapsto]

is continuous we have An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i103.gif

→ x^∞as m → 0 with respect to the euclidean norm, i.e. all m_ij→ 0 or equivalently n_ij→ ∞. Moreover, it follows that on subsets V [subset, dbl equals]

V' and U

U' the Jacobian of Φ is still negative definite and, consequently, the An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i103.gif

are stable steady-states. □

Manual parameter determination for the continuous T-cell model

The question is how to encode the information about slow and fast interactions in the numeric values of the parameters. We illustrate our approach using the subnetwork shown in Figure Figure2C,2C, where we focus on ZAP-70, LAT phosp and cCbl. The activation of LAT phosp stands for the interactions of scenario 2, whereas the activation of cCbl represents scenario 3. The ODE system for this network is given by

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i110.gif

For the activation of ZAP-70 k = 0.3 and τ_ZAP-70= 10 are used. LAT phosp and cCbl are activated at thresholds k_fast= 0.1 and k_slow= 0.8, respectively, and both their life-times are set to τ = 1. Initial condition for all species is 0 and the inputs for ZAP-70 are fixed at 1. The result of the numeric simulation is shown in Figure Figure2D.2D. From the beginning on ZAP-70 is activated but rather slowly due to the increased life-time τ_ZAP-70. The activations of LAT-phosp and cCbl occur at around the time points when the concentration of ZAP-70 crosses the thresholds k_fastand k_slow, respectively. The high threshold k_slowalso leads to a lower total activation level of cCbl. One can clearly see the time lag between the activations of LAT phosp and cCbl, i.e. between the interactions in scenario 2 and the interactions not occurring until scenario 3.

Model transformation and simulation

For the transformation of the Boolean T-cell model we choose HillCube ODEs. HillCubes are better suited to describe signaling cascades than BooleCubes. Normalization is not necessary as we will choose thresholds k [double less-than sign]

1 and consequently the Hill functions will already satisfy f(1) ≈ 1. For the transformation and simulation we developed a MATLAB toolbox called Odefy [Krumsiek et al.: Odefy — From discrete to continuous models. In preparation (2009)], which is publicly available at http://cmb.helmholtz-muenchen.de/odefy and allows the experimentalist to easily transform Boolean models into ODE models. Since Odefy can be integrated into CellNetAnalyzer [24], we were able to export the Boolean T-cell model from there. Numeric integration of the ODE system was carried out using MATLAB ode15s, a variable-order multistep solver based on the numerical differentiation formulas.

Experimental data

Kemp et al. [12] created a data set describing the dynamics of the activation of the key signaling elements ERK, JNK, IKK and NFAT upon activation of the TCR. The data was generated by stimulation of a T-cell line (1B6 T cell hybridoma) with three peptides with different affinities for the T-cell receptor, Q144, Y144 and L144. In the case of L144 experiments were conducted for three different peptide concentrations 0.04 μg/ml (low), 0.4 μg/ml (medium) and 4 μg/ml (high). In the case of Q144 and Y144 only the high concentration of 4 μg/ml was used. Concentrations of ERK, JNK, IKK and NFAT were measured at 0, 10, 30, 60, 120, 240 and 2400 minutes. Here we neglected the last time point, as on this slow time scale many interactions play a role, that are not included in the model, such as gene expression. For the parameter fitting the data were linearly rescaled to the unit interval.

Parameter fitting

The T-cell model consists of 40 species, 55 pairwise interactions and three external inputs. Hence, we have 40 life-time parameters and 58 pairs of Hill parameters, amounting to a total of 156 parameters. Due to this large number of parameters compared to the number of experimental data points, the fitting problem is obviously ill-posed as for many different parameter sets the model reproduces the data equally well. For this reason, we performed a two-step fitting process. First, we determined a parameter set for which the model fits the experimental data reasonably well. Second, we added a regularization to account for the indeterminacies. Both steps are optimization problems. The two cost functions are given below. They take a parameter set consisting of all Hill parameters (n, k) and all life-times τ as input a yield a scalar loss value, that needs to be minimized. In both steps we used a simulated annealing algorithm [25] for minimization. As the threshold parameters k have to be precisely adjusted at small values, these were fitted on a log-scale, as is also done in [26]. Parameters downstream of the measured species were, of course, not changed but fixed at their manually determined value. We used the SBPD package of the Systems Biology Toolbox [27] to create a compiled MATLAB simulation function of our ODE model for faster performance.

Least squares fitting

In the first step we determined a parameter set for which the model reproduces the data reasonably well. To this end, we employed a least squares fitting, i.e. we minimized the sum of the squared offsets of the data points from the model prediction. In a first attempt, we used only the offsets at the six time points t = 0, 10, 30, 60, 120, 240. This led to the model showing fast oscillations, which almost perfectly fitted the data, and, clearly, were an unrealistic overfitting. To avoid this, we linearly interpolated the experimental data and minimized the cost function

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i111.gif

where

An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i112.gif

denotes the time courses predicted by the model and An external file that holds a picture, illustration, etc.
Object name is 1752-0509-3-98-i113.gif

the interpolated data points. The simulated annealing algorithm was started from the manually determined parameter set and finally converged at a cost function value of 18.98.

Regularization

In the regularization step we minimized the sum of the coefficients of variation σ_n/μ_n, σ_k/μ_kand σ_τ/μ_τwithin the three parameter groups of n's, k's and τ's, respectively, under the constraint that the model's fit to the data did not deteriorate. The corresponding cost function is

A mathematical equation, expression, or formula.
Object name is 1752-0509-3-98-i114.gif

where the last term is a penalty term ensuring a constant quality of the model's fit. The simulated annealing was started from the result of the first step. Idea of this regularization is to account for parameter indeterminacies by reducing their variation and to enhance the significance of 'outliers' like the affinity k_L144of the TCR for the ligand L144, cf. Figure Figure4E4E.

Additional data file 6 shows model simulations for 5 different results of step 1. While not perfectly agreeing, their overall dynamic behavior is the same. We can reasonably assume that neither the parameter indeterminacies nor the regularization significantly influence the dynamics of the experimentally observed species.

Authors' contributions

DMW conceived and designed the methodology, performed the computational experiments and drafted the manuscript. JK helped to perform the computational experiments and developed the computational platform. JSR, DAL and SK interpreted the numerical and experimental data and contributed to the manuscript. FJT participated in the design of the methodology and of the computational platform and helped to draft the manuscript. All authors read and approved the final manuscript.

Supplementary Material

Additional file 1

Example for the hypergraph representation of a Boolean model. Supplementary text (.pdf) giving an example for the hypergraph representation of Boolean models.

Click here for file^{(42K, PDF)}

Additional file 2

Hill functions. Figure (.pdf) showing Hill functions f(x) = xⁿ/(xⁿ+ kⁿ) with Hill coefficients n = 2, 4, 8, 16, 32 and threshold k = 0.5.

Click here for file^{(7.4K, PDF)}

Additional file 3

HillCubes. Figure (.png) showing HillCubes of all 16 two-variable Boolean gates. Hill parameters are n = 3 and k = 0.5 for both inputs.

Click here for file^{(108K, PNG)}

Additional file 4

Regulatory domains. Figure (.pdf) showing the regulatory domains in a two-variable example.

Click here for file^{(8.6K, PDF)}

Additional file 5

Steady-states in discrete and continuous models. Supplementary text (.pdf) discussing a toy example where the steady-states of a discrete model are not preserved in a continuous version of the model.

Click here for file^{(317K, PDF)}

Additional file 6

Comparison of different best fit parameter sets with respect to model dynamics. Figure (.pdf) showing simulations of the continuous T-cell model for 5 different best fit parameter sets (without regularization, cf. section on parameter fitting). While not perfectly agreeing, the overall dynamic behavior is the same in all simulations.

Click here for file^{(43K, PDF)}

Acknowledgements

The authors would like to thank Melissa Kemp for providing the experimental data, Melody Morris for critical reading of the manuscript and useful discussions as well as the anonymous reviewers for their helpful comments during the preparation of this manuscript.

This research was partially supported by the Initiative and Networking Fund of the Helmholtz Association within the Helmholtz Alliance on Systems Biology (project CoReNe), by the German Federal Ministry of Education and Research (MaCS, Magdeburg Centre for Systems Biology) and by the Ministry of Education of Saxony-Anhalt (Research Center "Dynamic Systems"). JSR and DAL thank funding of Pfizer Inc., NIH grant GM68762, and DOD Institute for Collaborative Biotechnologies.

References

Kauffman SA. Metabolic stability and epigenesis in randomly constructed genetic nets. Journal of Theoretical Biology. 1969;22:437–467. doi: 10.1016/0022-5193(69)90015-0. [PubMed] [Cross Ref]
Fauré A, Naldi A, Chaouiya C, Thieffry D. Dynamical analysis of a generic Boolean model for the control of the mammalian cell cycle. Bioinformatics. 2006;22:124–131. doi: 10.1093/bioinformatics/btl210. [PubMed] [Cross Ref]
Davidich MI, Bornholdt S. Boolean network model predicts cell cycle sequence of fission yeast. PloS ONE. 2008;3 doi: 10.1371/journal.pone.0001672. [PMC free article] [PubMed] [Cross Ref]
Mendoza L, Thieffry D, Alvarez-Buylla ER. Genetic control of flower morphogenesis in Arabidopsis thaliana: a logical analysis. Bioinformatics. 1999;15:593–606. doi: 10.1093/bioinformatics/15.7.593. [PubMed] [Cross Ref]
Albert R, Othmer HG. The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster. Journal of Theoretical Biology. 2003;223:1–18. doi: 10.1016/S0022-5193(03)00035-3. [PubMed] [Cross Ref]
Klamt S, Saez-Rodriguez J, Lindquist JA, Simeoni L, Gilles ED. A methodology for the structural and functional analysis of signalling and regulatory networks. BMC Bioinformatics. 2006;7:56. doi: 10.1186/1471-2105-7-56. [PMC free article] [PubMed] [Cross Ref]
Chavez M, Albert R, Sontag E. Robustness and fragility of Boolean models for genetic regulatory networks. Journal of Theoretical Biology. 2005;235:431–449. doi: 10.1016/j.jtbi.2005.01.023. [PubMed] [Cross Ref]
Saez-Rodriguez J, Simeoni L, Lindquist JA, Hemenway R, Bommhardt U, Arndt B, Haus UU, Weismantel R, Gilles ED, Klamt S, Schraven B. A Logical Model Provides Insights into T Cell Receptor Signaling. PLoS Comput Biol. 2007;3:e163. doi: 10.1371/journal.pcbi.0030163. [PubMed] [Cross Ref]
Altan-Bonnet G, Germain R. Modeling T cell antigen discrimination based on feedback control of digital ERK responses. PLoS Biol. 2005;3:e356. doi: 10.1371/journal.pbio.0030356. [PMC free article] [PubMed] [Cross Ref]
Mendoza L, Xenarios I. A method for the generation of standardized qualitative dynamical systems of regulatory networks. Theoretical Biology and Medical Modelling. 2006;3 [PMC free article] [PubMed]
Glass L, Kauffman SA. The logical analysis of continuous, non-linear biochemical control networks. Journal of Theoretical Biology. 1973;39:103–129. doi: 10.1016/0022-5193(73)90208-7. [PubMed] [Cross Ref]
Kemp ML, Wille L, Lewis CL, Nicholson LB, Lauffenburger DA. Quantitative network signal combinations downstream of TCR activation can predict IL-2 production response. J Immunol. 2007;178:4984–4992. [PubMed]
Gasca M, Sauer T. On the history of multivariate polynomial interpolation. Journal of Computational and Applied Mathematics. 2000;122:23–35. doi: 10.1016/S0377-0427(00)00353-8. [Cross Ref]
Plahte E, Mestl T, Omholt S. A methodological basis for description and analysis of systems with complex switch-like interactions. Journal of Mathematical Biology. 1998;36:321–348. doi: 10.1007/s002850050103. [PubMed] [Cross Ref]
Hill A. The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves. J Physiol. 1910;40:4–7.
El Snoussi H, Thomas R. Logical identification of all steady states: The concept of feedback loop characteristic states. Bulletin of Mathematical Biology. 1993;55:973–991.
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. Bulletin of Mathematical Biology. 1995;57:247–276. [PubMed]
Werlen G, Hausmann B, Palmer E. A motif in the alphabeta T-cell receptor controls positive selection by modulating ERK activity. Nature. 2000;406:422–426. doi: 10.1038/35019094. [PubMed] [Cross Ref]
Werlen G. Signaling Life and Death in the Thymus: Timing Is Everything. Science. 2003;299:1859–1863. doi: 10.1126/science.1067833. [PubMed] [Cross Ref]
Daniels MA, Teixeiro E, Gill J, Hausmann B, Roubaty D, Holmberg K, Werlen G, Holländer GA, Gascoigne NRJ, Palmer E. Thymic selection threshold defined by compartmentalization of Ras/MAPK signalling. Nature. 2006;444:724–729. doi: 10.1038/nature05269. [PubMed] [Cross Ref]
Munder M, Bettelli E, Monney L, Slavik J, Nicholson L, Kuchroo V. Reduced Self-Reactivity of an Autoreactive T Cell After Activation with Cross-reactive Non-Self-Ligand. The Journal of Experimental Medicine. 2002;196:1151–1162. doi: 10.1084/jem.20020390. [PMC free article] [PubMed] [Cross Ref]
de Jong H, Gouzé J, Hernandez C, Page M, Sari T, Geiselmann J. Qualitative simulation of genetic regulatory networks using piecewise-linear models. Bulletin of Mathematical Biology. 2004;66:301–340. doi: 10.1016/j.bulm.2003.08.010. [PubMed] [Cross Ref]
Zadeh LA. Fuzzy sets. Information and Control. 1965;8:338–353. doi: 10.1016/S0019-9958(65)90241-X. [Cross Ref]
Klamt S, Saez-Rodriguez J, Gilles E. Structural and functional analysis of cellular networks with CellNetAnalyzer. BMC Systems Biology. 2007;1 [PMC free article] [PubMed]
Kirkpatrick S, Gelatt C, Vecchi M. Optimization by Simulated Annealing. Science. 1983;220:671–680. doi: 10.1126/science.220.4598.671. [PubMed] [Cross Ref]
Ma W, Lai L, Ouyang Q, Tang C. Robustness and modular design of the Drosophila segment polarity network. Molecular Systems Biology. 2006;2 [PMC free article] [PubMed]
Schmidt H, Jirstrand M. Systems Biology Toolbox for MATLAB: a computational platform for research in systems biology. Bioinformatics. 2006;22:514–515. doi: 10.1093/bioinformatics/bti799. [PubMed] [Cross Ref]
Huang C, Ferrell J. Ultrasensitivity in the mitogen-activated protein kinase cascade. Proceedings of the National Academy of Sciences. 1996;93:10078–10083. doi: 10.1073/pnas.93.19.10078. [PubMed] [Cross Ref]

Articles from BMC Systems Biology are provided here courtesy of
BioMed Central