Transforming Boolean models to continuous models: methodology and application to Tcell receptor signaling
 Dominik M Wittmann^{1, 2},
 Jan Krumsiek^{1},
 Julio SaezRodriguez^{3, 4},
 Douglas A Lauffenburger^{3},
 Steffen Klamt^{5} and
 Fabian J Theis^{1, 2, 6}Email author
DOI: 10.1186/17520509398
© Wittmann et al; licensee BioMed Central Ltd. 2009
Received: 19 January 2009
Accepted: 28 September 2009
Published: 28 September 2009
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 Tcells. We discuss how parameters for this model can be determined such that quantitative experimental results are explained and predicted, including timecourses 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 wellstudied subsystems.
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 topdown 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.
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, midtime events: Inputs are still 1, feedback loops are active.
A qualitative analysis of an expanded version of this Boolean model by SaezRodriguez et al. [8] yielded new and nonobvious signaling pathways. There are also quantitative models covering aspects of Tcell signaling in mechanistic detail, e.g. [9]. In contrast thereto, the continuous model we obtain in the following describes the Tcell signaling cascade on a larger scale yet in less mechanistic detail.
We apply the transformation method to the Tcell 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 steadystates 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 Tcell 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
giving the value of x_{ i }at the next time step for every possible combination of .
This is a socalled sumofproduct 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 sumofproduct form of B_{ i }by a directed hyperedge between a set of start nodes and the end node X_{ i }. Each pair (s, X_{ i }), s ∈ S, carries a sign — '+' or '' — depending on whether there is a factor s or ¬s in the product. All incoming hyperedges at node X_{ i }are 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_{ i }by continuous variables taking values in [0, 1], i.e. we normalize concentrations to the unit interval. Consequently, we have to 'extend' the functions B_{ i }to functions : [0, 1]^{ N }→ [0, 1]. We call the functions 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 numeric simulations the discretization of time is obviously irrelevant. It complicates, however, the detection of smallscale continuous effects and is a serious drawback in the further investigation of the model by analytical methods.
and by setting τ_{ i }= 1/γ_{ i }and 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 can be obtained from a Boolean function B_{ i }in a computationally efficient manner. A suitable transformation has to satisfy three conditions:
Accuracy: It has to be accurate, which means that and B_{ i }must agree on the vertices of the unit cube, i.e. for values in . Functions satisfying this condition are called perfect continuous homologues. As will be shown, for perfect continuous homologues discrete and continuous models exhibit a similar steadystate behavior.
Good analytical properties: The functions 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 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 is defined as a polynomial in the variables that agrees with B_{ i }on 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
by linearly interpolating the functions B_{ i }using the technique of multivariate polynomial interpolation as explained in the Methods section. These functions are called BooleCubes. By substituting for B_{ i }in equation (2), we can then define a system of ODEs that describes the temporal development of the .
HillCube
The functions are affine multilinear, i.e. for each 1 ≤ j ≤ N_{ i }and fixed , k ≠ j, there exist constants a, b ∈ ℝ such that . Molecular interactions, however, are known to show a switchlike behavior, which can be modeled using sigmoid shaped Hill functions [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.
which we call HillCubes. Plots of the HillCubes of all 16 twovariable Boolean gates can be found in Additional data file 3. Now a new system of ODEs can be defined by replacing the in equation (2) by the HillCubes .
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
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 in equation (2) are perfect homologues of the Boolean update functions B_{ i }, a steadystate of the Boolean model will also be a steadystate of the continuous system. This result can be applied to the BooleCube model as well as the normalized HillCube model, but not to the nonnormalized HillCube model. Therefore, Boolean steadystates will in general not be steadystates of the nonnormalized HillCube system. The question is if Boolean steadystates still correspond to 'similar' steadystates 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) steadystates may arise. Besides the steadystate behavior, monotony properties are a further important characteristic of a dynamical system. These properties determine the effect of a down or upregulation 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 Tcell 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 Tcell model
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 steadystate was reached.
Continuous Tcell model
We applied the transformation technique to the Boolean Tcell model. Using nonnormalized HillCubes we obtained a large quantitative model of Tcell 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 Tcells, 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 lifetimes to τ = 1, with the exception of Fyn and ZAP70. Here we knew that the interactions Fyn → PAGCsk and ZAP70 → 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 → TCRphos, ZAP70 → LATphosp, ZAP70 → PLCg (act) and ZAP70 → Itk were set to k_{fast} = 0.1. The Hill functions for the three thresholds k, k_{slow} and k_{fast} are displayed in Figure 2B. Finally, we set the lifetimes of Fyn and ZAP70 to τ_{Fyn} = τ_{ZAP70} = 10 to enlarge the time gap between the two switching points.
The numeric simulation of the continuous Tcell model with the manually determined parameters is shown in Figure 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 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 ZAP70 were activated more slowly. Nonetheless TCRphos, LATphosp, 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 ZAP70 were high enough, the feedback loops Fyn → PAGCsk and ZAP70 → cCbl became active and began to switch on PAGCsk and cCbl. While PAGCsk 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 steadystate 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 Tcell 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 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 3C) showed again a similar expression pattern.
We then analyzed this part of the cascade more deeply. Figure 3C 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 3A and 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 3C 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_{ActBeg} as the time when the first species (of the lower 29) reaches 5% of its maximal value,

the end of activation t_{ActEnd} as the time when the last species reaches 95% of its maximal value,

the beginning of deactivation t_{DeactBeg} as the time after t_{ActEnd} when the first species drops under 95% of its maximal value,

the end of deactivation t_{DeactEnd} as the time when the last species drops under 5% of its maximal value,
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 lifetimes τ_{Fyn} and τ_{ZAP70} on the one hand and the concentration thresholds k_{fast} and k_{slow} on the other hand. The remaining parameters were set to n = 3, k = 0.3 and τ = 1.
First, we computed ρ for fixed k_{fast} and k_{slow} and different τ_{Fyn} = τ_{ZAP70}. The result is shown in Figure 3D. For τ_{Fyn} = τ_{ZAP70} = 1 the cascade was not activated properly. For larger values we observed a decrease in ρ implying that an increase of τ_{Fyn} and τ_{ZAP70} prolonged the deactivation phase. This was to be expected — longer lifetimes resulted in a lessened increase of the decisive elements Fyn and ZAP70 in the regulatory loop.
Second, we analyzed the effect of k_{fast} and k_{slow} for fixed τ_{Fyn} = τ_{ZAP70} = 10. The result is shown in Figure 3E. Only for parameters k_{fast} ≪ 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 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 3D and 3E were = 0.24 and = 0.27, respectively, which were both significantly smaller then ρ_{b} = 1. This suggests that ρ ≪ 1 may be an invariant of the dynamical system of biological importance.
Explanation of experimental data using the continuous Tcell model
Best fit parameter set for ligand L144.
Species  Input  n  k  Input  n  k  Input  n  k  τ 

zap70  tcrphosp  2  0.29  lck  1  0.33  ccbl  1  0.30  2.77 
tcrphosp  tcr  1  0.34  lck  1  0.34  fyn  1  0.13  0.03 
tcr  TCRlig  1  0.20  ccbl  1  0.41  0.59  
slp76  gads  1  0.35  0.24  
sek  pkcth  1  0.26  0.01  
rsk  erk  3  0.30  1.00  
rlk  lck  1  0.36  1.02  
rasgrp  pkcth  1  0.40  dag  1  0.39  2.76  
ras  rasgrp  3  0.33  grb2sos  17  0.39  7.50  
raf  ras  3  0.31  11.80  
plcgbind  lat  5  0.36  0.01  
plcgact  zap70  1  0.10  slp76  1  0.29  rlk  1  0.24  0.01 
plcgbind  1  0.31  itk  2  0.34  
pkcth  dag  2  0.31  8.66  
pagcsk  tcr  1  0.32  fyn  1  0.69  1.33  
nfkb  ikb  3  0.30  1.00  
nfat  calcen  1  0.37  5.03  
mek  raf  1  0.28  0.01  
TCRlig  ext. sig.  1  0.12  1.28  
lck  pagcsk  2  0.30  cd4  2  0.37  cd45  1  0.33  0.15 
lat  zap70  3  0.14  0.01  
jun  jnk  3  0.30  1.00  
jnk  sek  4  0.57  0.08  
itk  zap70  2  0.10  slp76  1  0.31  0.01  
ip3  plcgact  2  0.38  0.01  
ikkbeta  pkcth  18  0.33  10.55  
ikb  ikkbeta  3  0.30  1.00  
grb2sos  lat  8  0.46  5.26  
gads  lat  1  0.36  0.12  
fyn  tcr  6  0.27  lck  2  0.30  cd45  4  0.37  0.46 
fos  erk  3  0.30  1.00  
erk  mek  2  0.28  0.10  
dag  plcgact  3  0.28  0.01  
creb  rsk  3  0.30  1.00  
cre  creb  3  0.30  1.00  
cd4  ext. sig.  1  0.32  2.25  
cd45  ext. sig.  1  0.35  1.64  
ccbl  zap70  1  0.70  0.29  
calcen  ca  11  0.30  4.53  
ca  ip3  3  0.33  0.20  
ap1  jun  3  0.30  fos  3  0.30  1.00 
Best fit parameter set
Subsequently, we analyzed the distribution of the Hill parameters within the best fit parameter set for ligand L144 (Table 1). When looking at the distribution of the exponents n upstream of the measured species (Figure 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 Booleanlike switches (n → ∞). Figure 4E 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 4E. Only one other parameter k also had a Zscore above 2: the threshold k_{L144} for 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 quasiconstant, 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_{L144} for the stimulation of the TCR by the ligand L144 was significantly below the mean of the other thresholds k, cf. marker 'o' in Figure 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_{Q144} and 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_{Q144} and k_{Y144} we obtained were significantly above the mean of the other thresholds k, cf. markers '◇' and '□' in Figure 4E. This suggests that the predicted relation between the parameters k_{L144} ≪ k_{Y144} <k_{Q144} is 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 steadystate 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 onevariable Boolean step functions. In the case of multivariable Boolean functions B_{ i }a (perfect) continuous homologue is constructed as follows: Note that, when building a Boolean model, one implicitly introduces a threshold 0 <θ_{ i }< 1 for each species X_{ i }and defines its state as 'on' if its concentration is above this threshold. The hyperplanes = θ_{ i }, i = 1, 2, ..., N decompose the cubes into rectangular regions called domains. Each of these domains contains exactly one vertex and is denoted by accordingly (see Additional data file 4). A simple way of defining the functions is now to set
With this definition model (2) is a socalled piecewise linear ODE model which means that within each domain equation (2) is a linear ODE of the form . This kind of equation is very well understood and can be solved analytically. The functions defined in (6) perfectly agree with B_{ i }on 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 steadystates are also steadystates 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 steplike transition between the different domains is often unrealistic. Rather they are analyzed in a qualitative and semiqualitative 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 socalled degree of membership (DOM) function is introduced for each species X_{ i }. For concentrations 0 = = 1 this function gives the degree with which we say that X_{ i }is 'on'. There are two standard ways of generalizing the Boolean operators AND, OR, NOT:
Standardized qualitative dynamical systems
The righthand 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 ω_{ i }represents 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 . Actually, a more general form of ω_{ i }is 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 ω_{ i }is. 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, steadystates of the Boolean model are not necessarily steadystates 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 (see Methods). This technique can be applied to any Boolean function B_{ i }. The resulting BooleCubes are smooth and can easily be analytically differentiated and integrated. They agree with the Boolean functions B_{ i }on the vertices of the unit cube and hence the steadystates from the Boolean model are also steadystates in the continuous model.
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 from equation (5).
Comparison
To conclude, we illustrate the above methods applied to a simple OR gate between two species X_{1} and X_{2}. We compute
from equation (6),

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

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

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

the BooleCube from equation (3) obtained by the interpolation technique,

the HillCube from equation (4) for Hill functions f_{1} and f_{2} with parameters n = 3, k = 0.5,
of from equation (5).
Figures 5C and 5D show the productsum fuzzy logic function and the input function ω. One can clearly see that they do not represent a pure OR gate, where the values at (x_{1}, x_{2}) = (1, 0) and (x_{1}, x_{2}) = (0, 1) should already be maximal. This is the case in Figures 5A and 5B which show the piecewise linear and the minmax fuzzy logic function . Here however, the problem is that the functions are not differentiable, as can easily be seen from their plots. The BooleCube shown in Figure 5E is both, smooth and maximal as soon as any concentration is equal to 1. Finally, Figures 5F and 5G show the (normalized) HillCubes and , 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 5H.
Theoretical results about steadystates
A fundamental principle of biological modeling is that steadystates 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 steadystates of the Boolean model are still steadystates in the homologue continuous system. In the case of a perfect agreement of the Boolean update rules B_{ i }with their continuous homologues 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 realworld examples, the steadystates of a Boolean model correspond to steadystates 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 steadystates 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 steadystate of the continuous system in the neighborhood of each Boolean steadystate (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) ∨ D would be translated into two positive arcs (C → A, D → A) and one negative arc (B ⊣ 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 Tcell model tells us, that LATphosp exerts purely positive effects on ERK, JNK, IKK, and NFAT, because there are only positive paths from LATphosp 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 LATphosp 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 selfloops identical and monotony properties are therefore preserved. Accordingly, a positive perturbation in LATphosp, e.g. by permanently decreasing the threshold of ZAP70 in the interaction activating LATphosp, results in a trajectory that is always above the trajectory of the nonperturbed system (Figure 3F). In fact, in accordance with the interaction graph of the Boolean model, we observe purely positive effects on all species downstream of LATphosp 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 steadystate 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 Tcell receptor signaling. The resulting model is an extensive continuous model of Tcell 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 sstate 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 ∈ {1, 2, ..., N} is fixed and for the sake of simplicity the subscript i is omitted. We remark that here B can be any realvalued 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 : ℝ^{ N }→ ℝ that is a continuation of B_{ i }in the sense that for all .
We now show that 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 has the following properties:
 (i)
 (ii)
 (iii)
Let s denote the number of symmetry hyperplanes of B, i.e. the number of variables x _{ i } satisfying
B(x_{1}, ..., x_{i1}, 0, x_{i+1}, ..., x_{ N }) = B(x_{1}, ..., x_{i1}, 1, x_{i+1}, ..., x_{ N }), for all (x_{1}, ..., , ..., x_{ N }) ε {0, 1}^{N1}. Then deg ( ) = N  s.
 (iv)
is affine multilinear. It is multilinear iff it corresponds to an AND gate, i.e.B(x_{1}, x_{2}, ..., x_{ N }) = .
 (v)
In particular, if B is a Boolean function.
(ii) A minimal degree interpolation polynomial is of the form
since exponents greater than 1 can be replaced by 1 without changing the values of f for . 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 . Then we consider the sequence of equations f(V_{ k }) = B(V_{ k }), k = 1, 2, ..., 2^{ N }. For all k ∈ {1, 2, ..., 2^{ N }} there is exactly one coefficient whose monomial satisfies g(V_{ k }) = 0, k = 1, 2, ..., k  1 and g(V_{ k }) = 1. This allows to uniquely determine the coefficients and that way also the polynomial f. Since is of the form (*), it is the unique minimal degree interpolation polynomial.
(iii) The degree of is clearly deg ( ) ≤ 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_{2}, ..., x_{ N }) = B(1, x_{2}, ..., x_{ N }) for all (x_{2}, ..., x_{ N }) ∈ {0, 1}^{N1}and consequently
Hence, deg ( ) ≤ N  1. Inductively this proves (iii).
(iv)Let i ∈ {1, 2, ..., N}. Then at fixed the derivative is constant so is affine linear in . Moreover, we have
For the last equivalence note that if is multilinear, then if any = 0.
(v) Assume and w.l.o.g. . Then it follows from (iv) that or . Inductively, we obtain x ∈ {0, 1}^{ N }such that , a contradiction to (i). Analogously, one proves that . □
This function does not have any symmetry hyperplanes, but its interpolation has degree 2.
Theoretical results about steadystates
The following theorem investigates the steadystate behavior of discrete and continuous models.
Theorem. Assume we are given a Boolean model and perfect continuous homologues of the Boolean update functions. Then for any state vector the following are equivalent:
(i) is a steadystate of the Boolean model.
(ii) is a steadystate of the model (1).
(iii) is a steadystate of the ODE model (2).
for model (2). Considering that the 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. from equation (6), Glass et al. [11] could show that the above theorem also holds when (iii) is replaced by ' is a stable steadystate 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 and B_{ i }. The main idea is to use the implicit function theorem to prove the existence of a steadystate of the continuous model in a neighborhood of a Boolean steadystate.
Then there are open neighborhoods U_{1}, U_{2} of 0 and 1, respectively, such that f( , 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 . For concentrations 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 ⊂ ℝ^{ N }of x^{∞} such that Φ is continuously differentiable on ℝ^{m}× U (as the composition of continuously differentiable functions).
where δ_{ ij }is 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' ⊂ ℝ^{m}of 0 such that for each m ∈ V', i.e. for large exponents n, there is a ∈ U' such that Φ(m, ) = 0, i.e. is a steadystate of the model. Since the mapping m ↦ is continuous we have → 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 ⊆ V' and U ⊆ U' the Jacobian of Φ is still negative definite and, consequently, the are stable steadystates. □
Manual parameter determination for the continuous Tcell model
For the activation of ZAP70 k = 0.3 and τ_{ZAP70}= 10 are used. LAT phosp and cCbl are activated at thresholds k_{fast} = 0.1 and k_{slow} = 0.8, respectively, and both their lifetimes are set to τ = 1. Initial condition for all species is 0 and the inputs for ZAP70 are fixed at 1. The result of the numeric simulation is shown in Figure 2D. From the beginning on ZAP70 is activated but rather slowly due to the increased lifetime τ_{ZAP70}. The activations of LATphosp and cCbl occur at around the time points when the concentration of ZAP70 crosses the thresholds k_{fast} and k_{slow}, respectively. The high threshold k_{slow} also 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 Tcell 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 ≪ 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.helmholtzmuenchen.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 Tcell model from there. Numeric integration of the ODE system was carried out using MATLAB ode15s, a variableorder 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 Tcell line (1B6 T cell hybridoma) with three peptides with different affinities for the Tcell 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 Tcell model consists of 40 species, 55 pairwise interactions and three external inputs. Hence, we have 40 lifetime 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 illposed as for many different parameter sets the model reproduces the data equally well. For this reason, we performed a twostep 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 lifetimes τ 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 logscale, 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
where denotes the time courses predicted by the model and 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
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_{L 144}of the TCR for the ligand L144, cf. Figure 4E.
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.
Declarations
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 SaxonyAnhalt (Research Center "Dynamic Systems"). JSR and DAL thank funding of Pfizer Inc., NIH grant GM68762, and DOD Institute for Collaborative Biotechnologies.
Authors’ Affiliations
References
 Kauffman SA: Metabolic stability and epigenesis in randomly constructed genetic nets. Journal of Theoretical Biology. 1969, 22 (3): 437467. 10.1016/00225193(69)900150View ArticlePubMedGoogle Scholar
 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 (14): 124131. 10.1093/bioinformatics/btl210View ArticleGoogle Scholar
 Davidich MI, Bornholdt S: Boolean network model predicts cell cycle sequence of fission yeast. PloS ONE. 2008, 3 (2):
 Mendoza L, Thieffry D, AlvarezBuylla ER: Genetic control of flower morphogenesis in Arabidopsis thaliana: a logical analysis. Bioinformatics. 1999, 15 (7/8): 593606. 10.1093/bioinformatics/15.7.593View ArticlePubMedGoogle Scholar
 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: 118. 10.1016/S00225193(03)000353View ArticlePubMedGoogle Scholar
 Klamt S, SaezRodriguez 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 10.1186/14712105756PubMed CentralView ArticlePubMedGoogle Scholar
 Chavez M, Albert R, Sontag E: Robustness and fragility of Boolean models for genetic regulatory networks. Journal of Theoretical Biology. 2005, 235: 431449. 10.1016/j.jtbi.2005.01.023View ArticleGoogle Scholar
 SaezRodriguez 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 (8): e163 10.1371/journal.pcbi.0030163PubMed CentralView ArticlePubMedGoogle Scholar
 AltanBonnet G, Germain R: Modeling T cell antigen discrimination based on feedback control of digital ERK responses. PLoS Biol. 2005, 3 (11): e356 10.1371/journal.pbio.0030356PubMed CentralView ArticlePubMedGoogle Scholar
 Mendoza L, Xenarios I: A method for the generation of standardized qualitative dynamical systems of regulatory networks. Theoretical Biology and Medical Modelling. 2006, 3 (13):
 Glass L, Kauffman SA: The logical analysis of continuous, nonlinear biochemical control networks. Journal of Theoretical Biology. 1973, 39: 103129. 10.1016/00225193(73)902087View ArticlePubMedGoogle Scholar
 Kemp ML, Wille L, Lewis CL, Nicholson LB, Lauffenburger DA: Quantitative network signal combinations downstream of TCR activation can predict IL2 production response. J Immunol. 2007, 178 (8): 49844992.View ArticlePubMedGoogle Scholar
 Gasca M, Sauer T: On the history of multivariate polynomial interpolation. Journal of Computational and Applied Mathematics. 2000, 122 (12): 2335. 10.1016/S03770427(00)003538.View ArticleGoogle Scholar
 Plahte E, Mestl T, Omholt S: A methodological basis for description and analysis of systems with complex switchlike interactions. Journal of Mathematical Biology. 1998, 36 (4): 321348. 10.1007/s002850050103View ArticlePubMedGoogle Scholar
 Hill A: The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves. J Physiol. 1910, 40: 47.Google Scholar
 El Snoussi H, Thomas R: Logical identification of all steady states: The concept of feedback loop characteristic states. Bulletin of Mathematical Biology. 1993, 55 (5): 973991.View ArticleGoogle Scholar
 Thomas R, Thieffry D, Kaufman M: Dynamical behaviour of biological regulatory networksI. Biological role of feedback loops and practical use of the concept of the loopcharacteristic state. Bulletin of Mathematical Biology. 1995, 57 (2): 247276.View ArticlePubMedGoogle Scholar
 Werlen G, Hausmann B, Palmer E: A motif in the alphabeta Tcell receptor controls positive selection by modulating ERK activity. Nature. 2000, 406 (6794): 422426. 10.1038/35019094View ArticlePubMedGoogle Scholar
 Werlen G: Signaling Life and Death in the Thymus: Timing Is Everything. Science. 2003, 299 (5614): 18591863. 10.1126/science.1067833View ArticlePubMedGoogle Scholar
 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 (7120): 724729. 10.1038/nature05269View ArticlePubMedGoogle Scholar
 Munder M, Bettelli E, Monney L, Slavik J, Nicholson L, Kuchroo V: Reduced SelfReactivity of an Autoreactive T Cell After Activation with Crossreactive NonSelfLigand. The Journal of Experimental Medicine. 2002, 196 (9): 11511162. 10.1084/jem.20020390PubMed CentralView ArticlePubMedGoogle Scholar
 de Jong H, Gouzé J, Hernandez C, Page M, Sari T, Geiselmann J: Qualitative simulation of genetic regulatory networks using piecewiselinear models. Bulletin of Mathematical Biology. 2004, 66 (2): 301340. 10.1016/j.bulm.2003.08.010View ArticlePubMedGoogle Scholar
 Zadeh LA: Fuzzy sets. Information and Control. 1965, 8 (3): 338353. 10.1016/S00199958(65)90241X.View ArticleGoogle Scholar
 Klamt S, SaezRodriguez J, Gilles E: Structural and functional analysis of cellular networks with CellNetAnalyzer. BMC Systems Biology. 2007, 1 (2):
 Kirkpatrick S, Gelatt C, Vecchi M: Optimization by Simulated Annealing. Science. 1983, 220 (4598): 671680. 10.1126/science.220.4598.671View ArticlePubMedGoogle Scholar
 Ma W, Lai L, Ouyang Q, Tang C: Robustness and modular design of the Drosophila segment polarity network. Molecular Systems Biology. 2006, 2 (70):
 Schmidt H, Jirstrand M: Systems Biology Toolbox for MATLAB: a computational platform for research in systems biology. Bioinformatics. 2006, 22 (4): 514515. 10.1093/bioinformatics/bti799View ArticlePubMedGoogle Scholar
 Huang C, Ferrell J: Ultrasensitivity in the mitogenactivated protein kinase cascade. Proceedings of the National Academy of Sciences. 1996, 93 (19): 1007810083. 10.1073/pnas.93.19.10078.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.