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

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.

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 results and to design and optimize further experiments; the results of these experiments, in turn, help to refine the model.

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 scenario 3 (activation and feedback). (D) Numeric simulation of the subnet from (C). Time is plotted on a log-scale. The two dotted vertical lines indicate the time points when the concentration of ZAP-70 exceeds k_fast and k_slow, respectively.

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 factors, the activation of the feedback loops and the total deactivation of the transcription factors. (A) Boolean simulation using synchronous updates. (B) Continuous simulation using HillCubes. (C) Upper figure: difference between the time course from (A) and the time course from (B). Lower figure: correlation between the discrete and the continuous expression pattern of the lower 29 species at each time point 0 <t < 100. The curve was smoothened using a moving-average filter with 5 time-units window length. (D) Ratio of activation and deactivation time, ρ, for τ_Fyn = τ_ZAP-70 = 1, 2, ..., 20. Other parameters are set to n = 3, k = 0.3, τ = 1, k_fast = 0.1 and k_slow = 0.8. For τ_Fyn = τ_ZAP-70 = 1 the cascade was not activated properly. (E) Ratio of activation and deactivation time, ρ, for 20 different k_fast and k_slow. log10 (k_fast) and log₁₀ (k_slow) are uniform from [-2, 0]. Other parameters are set to n = 3, k = 0.3, τ = 1, τ_Fyn = τ_ZAP-70 = 10. In the white areas the cascade was not activated properly. (F) Comparison between time course y(t) with the manually determined parameters and time course y³(t) with the decreased threshold in the activation of LAT-phosp by ZAP-70 as described in the section about monotony properties. Plotted is the difference y'(t) - y(t).

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 '○'. (A) ERK. (B) JNK. (C) NFAT. (D) IKK. (E, F) Distribution of the Hill parameters in the best fit parameter set for ligand L144 (Table 1). (E) Distribution of the thresholds k. The markers '+' indicate the position of the parameters which we set to k_fast (blue) and k_slow (red) at the beginning of the optimization. The marker '○' indicates the position of k_L144; for comparability the thresholds k_Q144 and k_Y144 are also indicated, cf. the markers '◇' and '□', respectively. (F) Distribution of the Hill exponents n.

Figure 5

Comparison of different transformation techniques. Different continuous homologues of a Boolean OR gate. (A) Piecewise linear function . (B) Function obtained by min-max fuzzy logic and linear DOM functions. (C) Function obtained by product-sum fuzzy logic and linear DOM functions. (D) Input function ω introduced by Mendoza et al. [10]. (E) BooleCube obtained by multivariate polynomial interpolation. (F) HillCube . (G) normalized HillCube . In the last two figures parameters n = 3 and k = 0.5 were chosen for both inputs. (H) Overview of the different transformation techniques with respect to their analytical properties and transformation accuracy.

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 τ_i are set to 1. As initial conditions we take = c > 0, = 0, = 0, = 0, for some constant input concentration c. The input node X₁ remains constant and the other concentrations change accordingly to the ODE , i = 2, 3, 4. We simulate the model for different Hill coefficients n = 1, 4, 16 and input level c = 1; the results are shown in (A), (B) and (C). All three time courses show qualitatively the same cascade-like pattern. With growing n, however, the onset of activation of X₃ and X₄ comes closer and closer to the time point at which their activators X₂ and X₃, respectively, cross the threshold k. (D) shows the input-output curve. Plotted is the (constant) input concentration c of node X₁ against the steady-state concentration of node X₄. For n > 1, we observe the typical sigmoid stimulus-response behavior of signaling cascades, see e.g. [28]. With increasing n the steepness of the input-output curve increases, leading to an almost discrete (Boolean) output in the case n = 16.