Continuous time Boolean modeling for biological signaling: application of Gillespie algorithm

Abstract

Mathematical modeling is used as a Systems Biology tool to answer biological questions, and more precisely, to validate a network that describes biological observations and predict the effect of perturbations. This article presents an algorithm for modeling biological networks in a discrete framework with continuous time.

Background: There exist two major types of mathematical modeling approaches: (1) quantitative modeling, representing various chemical species concentrations by real numbers, mainly based on differential equations and chemical kinetics formalism; (2) and qualitative modeling, representing chemical species concentrations or activities by a finite set of discrete values. Both approaches answer particular (and often different) biological questions. Qualitative modeling approach permits a simple and less detailed description of the biological systems, efficiently describes stable state identification but remains inconvenient in describing the transient kinetics leading to these states. In this context, time is represented by discrete steps. Quantitative modeling, on the other hand, can describe more accurately the dynamical behavior of biological processes as it follows the evolution of concentration or activities of chemical species as a function of time, but requires an important amount of information on the parameters difficult to find in the literature.

Results: Here, we propose a modeling framework based on a qualitative approach that is intrinsically continuous in time. The algorithm presented in this article fills the gap between qualitative and quantitative modeling. It is based on continuous time Markov process applied on a Boolean state space. In order to describe the temporal evolution of the biological process we wish to model, we explicitly specify the transition rates for each node. For that purpose, we built a language that can be seen as a generalization of Boolean equations. Mathematically, this approach can be translated in a set of ordinary differential equations on probability distributions. We developed a C++ software, MaBoSS, that is able to simulate such a system by applying Kinetic Monte-Carlo (or Gillespie algorithm) on the Boolean state space. This software, parallelized and optimized, computes the temporal evolution of probability distributions and estimates stationary distributions.

Conclusions: Applications of the Boolean Kinetic Monte-Carlo are demonstrated for three qualitative models: a toy model, a published model of p53/Mdm2 interaction and a published model of the mammalian cell cycle. Our approach allows to describe kinetic phenomena which were difficult to handle in the original models. In particular, transient effects are represented by time dependent probability distributions, interpretable in terms of cell populations.

Figure 1

Comparison of tools for discrete modeling, biological implication. Comparison table of the following tools: MaBoSS, GINsim, CellNetAnalyzer, BoolNet, GNA, SQUAD. Technical aspects are provided, along with the inputs/outputs relations between a model and data. The last row illustrates graphically the typical outputs that can be obtained from each tool.

Figure 2

Toy model. Toy model of a single cycle. (a) Influence network. (b) Logical rules and transition rates of the model. (c) Simulation parameters.

Figure 3

Transition graph of the toy model. Transition graph for the toy model (generated by GINsim). The node states should be read as [ABC] = [^∗∗∗]. [ABC]=[100] corresponds to a state in which only A is active. The nodes in green belong to a cycle, the node in red is the fixed point and the other nodes are in blue.

Figure 4

MaBoSS outputs of the toy model with fast escape rate. BKMC algorithm is applied to the toy model, with a fast escape rate. Trajectory of the network state probabilities [ABC]=[000] and [ABC]=[1^∗∗] (where ^∗can be either 0 or 1), the entropy (H) and the transition entropy (TH) are plotted. Because the probability of [ABC]=[000] converges to 1, [ABC]=[000] is a fixed point. The asymptotic behavior of both the entropy and the transition entropy is also the signature of a fixed point.

Figure 5

MaBoSS outputs of the toy model with slow escape rate. BKMC algorithm is applied to the toy model, with a slow escape rate. Trajectory of the network state probabilities [ABC]=[000] and [ABC]=[1**], the entropy (H) and the transition entropy (TH) are plotted. The asymptotic behavior of both the entropy and the transition entropy seems to be the signature of a cycle.

Figure 6

MaBoSS outputs of toy model with slow escape rate, large time scale. BKMC algorithm is applied to the toy model, with a slow escape rate, plotted on a larger time scale. Trajectory of probabilities ([ABC]=[000] and [ABC]=[1**]), the entropy (H) and the transition entropy (TH) are plotted. On a large time scale, the asymptotic behavior of both the entropy and the transition entropy is similar to the case of the fast escape rate (Figure 3).

Figure 7

Model of p53 response to DNA damage. Model of p53 response to DNA damage. (a) Influence network. (b) Logical rules and transition rates of the model. (c) Simulation parameters.

Figure 8

Transition graph of the model of p53 response to DNA damage. Transition graph of the p53 model (generated by GINsim). The node states should be read as [p53 Mdm2C Mdm2N Dam] = [^∗∗∗∗] (where ^∗can be either 0 or 1). For instance, [p53 Mdm2C Mdm2N Dam]=[1000] corresponds to a state in which only p53 (at its level 1) is active. The nodes in green and the nodes in light blue belong to two cycles, the node in red is the fixed point and the other nodes are in dark blue.

Figure 9

MaBoSS outputs of the model of p53 response to DNA damage. Trajectories of the network state probabilities of [p53 Mdm2C Mdm2N Dam] = [1^∗∗∗] and of [p53 Mdm2C Mdm2N Dam] = [2^∗∗∗], the entropy (H) and the transition entropy (TH) are plotted.

Figure 10

MaBoSS outputs of the model of the mammalian cell cycle: trajectories of probabilities. BKMC algorithm is applied to the mammalian cell cycle model, with an initial condition corresponding to a G1 state in the presence of growth factors (CyclinD is on). Trajectories of the cyclins probabilities, the entropy (H), transition entropy (TH) are plotted. The asymptotic behavior corresponds to the first indecomposable stationary distribution identified in Figure 10.

Figure 11

MaBoSS outputs of the model of the mammalian cell cycle: stationary distributions. BKMC algorithm is applied to the mammalian cell cycle model, with random initial conditions. Results of the clustering algorithm that associates a cluster to each indecomposable stationary distribution. (a) Probability of reaching each identified cluster; these probabilities are estimated by the proportion of trajectories that belong to each cluster. (b) First estimated cluster that can be interpreted as a desynchronized population of cells that are dividing. (c) Second estimated cluster, corresponding to a fixed point, that can be interpreted as a G1 cell cycle arrest with no growth factors.