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

Background

• ODE parameters are often fitted to expected phenotypes, and stochastic chemical-model translation can make initial conditions and computation time explode.ODEs also focus on average behavior, limiting their suitability for heterogeneous cell populations.Source: Background · p. 2
• Boolean models encode genes and proteins as 0s and 1s and assign each node a logical rule, but asynchronous updates can produce an explosion of solutions.Their transitions also function as time steps rather than characterizing biological time.Source: Background · p. 2
• BKMC generalizes asynchronous Boolean dynamics with biological time, stochastic update order, and a numerical rate for each transition.It combines advantages of discrete and continuous frameworks.Source: Background · p. 2
• BKMC is suited to influence-network models, processes with known biological timing, and heterogeneous cell populations.The passage cites cell-cycle models as an influence-network example and developmental biology as a setting with timing information.Source: Background · p. 2
• The present implementation supports models of up to 64 nodes, enabling finer descriptions of transient effects in networks already modeled with discrete-time Boolean methods.The authors describe MaBoSS, a C++ implementation, and illustrate BKMC with a toy model and a published p53-MDM2 model.Source: Background · p. 3

Continuous time in Boolean modeling: past and present

• Explicit-time discrete-step methods assign durations to node states, but noisy durations can make results sensitive to noise and initial conditions and require careful tuning.The authors note that experimental data for these durations are difficult to extract.Source: Continuous time in Boolean modeling:… · p. 3
• Implicit-time methods assign probabilities to nondeterministic transitions, but the resulting dynamics remain difficult to interpret in biological time.Priority rules can be represented by assigning zero probability to some transitions.Source: Continuous time in Boolean modeling:… · p. 3
• BKMC is proposed as an alternative to those explicit- and implicit-time approaches.The authors say their goal is to represent the time evolution of Boolean network states using real time.Source: Continuous time in Boolean modeling:… · p. 3

Properties of BKMC algorithm

• BKMC keeps Boolean node and network states and logic-based updates, while representing time as real-valued and evolution as stochastic.It uses a continuous-time Markov process with transition rates in a master equation; specified transitions follow asynchronous Boolean dynamics.Source: Properties of BKMC algorithm · p. 4

Markov process for Boolean model

• A Boolean network state is a vector of node values, and the set of possible states forms the state space over which the stochastic process evolves.At each time, the process assigns probabilities to states; these probabilities alone do not specify the full process.Source: Markov process for Boolean model · p. 4
• The Markov property means future conditional probabilities depend on the present rather than the past.The process is specified by an initial condition and conditional probabilities.Source: Markov process for Boolean model · p. 4
• Discrete-time Markov processes use transition probabilities, whereas continuous-time processes use transition rates; BKMC considers time-independent rates.For time-independent rates, an edge in the transition graph corresponds to a positive transition rate.Source: Markov process for Boolean model · p. 5

Asynchronous Boolean dynamics as a discrete time Markov process

• In asynchronous Boolean dynamics, each node has a logic based on its incoming links, and a permitted transition updates one node.With γ(S) possible outgoing asynchronous transitions from state S, each permitted transition receives probability 1/γ(S), while non-permitted transitions receive zero.Source: Asynchronous Boolean dynamics as a… · p. 5
• This construction defines a time-independent discrete-time Markov process once an initial condition is specified.Other definitions of γ(S) can yield a different Markov process with the same transition graph.Source: Asynchronous Boolean dynamics as a… · p. 5

Continuous time Markov process as a generalization of asynchronous Boolean Dynamics

• The continuous-time generalization replaces transition probabilities with rates and computes conditional probabilities by solving a master equation.To retain asynchronous dynamics, transitions are nonzero only when the two states differ at one node.Source: Continuous time Markov process as a… · pp. 5–6
• For the node that changes state, the model assigns activation or inhibition rates, which together with an initial condition define the continuous-time process.The rates are specified for node transitions.Source: Continuous time Markov process as a… · p. 6

Asymptotic behavior of continuous time Markov process

• A complete description of asymptotic behavior is a linear combination of indecomposable stationary distributions, and instantaneous probabilities in a stationary process are time-independent.Instantaneous probabilities converge to a stationary distribution.Source: Asymptotic behavior of continuous time… · p. 6

Oscillations and cycles

• A cycle with no outgoing edge corresponds to an indecomposable stationary distribution, regardless of transition-rate values.A cycle is defined topologically as a loop in the transition graph.Source: Oscillations and cycles · p. 6
• Damped oscillations require a transition matrix with at least one non-real eigenvalue, which in turn requires a cycle, but a cycle is not sufficient.When a non-real eigenvalue exists, some initial condition produces damped oscillations; instantaneous probabilities cannot be perfectly periodic.Source: Oscillations and cycles · p. 6
• In the single-cycle toy model, the same transition graph may or may not produce non-real eigenvalues depending on the transition rates.Source: Oscillations and cycles · p. 6

• The transition matrix has size 2^n × 2^n, making exact master-equation computation impractical when n is large.The master equation describes the evolution of conditional and joint probabilities.Source: BKMC: Kinetic Monte-Carlo (Gillespie… · p. 7
• Kinetic Monte-Carlo generates stochastic trajectories that represent the Markov process, allowing probabilities to be estimated from a finite sample.Each step produces a transition time and new state from the current state using two uniform random numbers, and trajectories continue to the specified maximum time.Source: BKMC: Kinetic Monte-Carlo (Gillespie… · p. 7
• BKMC applies this Gillespie procedure to continuous-time Markov processes in Boolean network state space.At each step, possible successor states are ordered with their transition rates, and a rate-weighted selection determines the next state.Source: BKMC: Kinetic Monte-Carlo (Gillespie… · p. 7

Practical use of BKMC, through MaBoSS tool

• MaBoSS applies BKMC to influence networks with node-specific logical rules and on/off transition rates, using network and rate inputs plus simulation parameters.The rates illustrate relative reaction speeds without requiring the exact parameter values expected in nonlinear ODE models.Source: Practical use of BKMC, through MaBoSS… · p. 7

Transition rates with MaBoSS language: modeling biological processes

• MaBoSS rate expressions can distinguish input-specific activation speeds and encode a separate synergistic rate when both inputs are active.The example assigns C activation rates $ka, $kb, and $kab for A alone, B alone, and both together, respectively.Source: Transition rates with MaBoSS language:… · p. 8
• Multi-level activity can be represented with separate Boolean nodes, while differing activation and inhibition speeds can model persistence after an input shuts down.In the two-level example, B_h activates when A and B are active and is inactivated when A turns off.Source: Transition rates with MaBoSS language:… · p. 8

Simulation input parameters in MaBoSS

• MaBoSS parameters require tuning to balance estimate convergence and runtime; larger time windows and sample counts improve probability-estimate convergence.Longer maximum times improve stationary-distribution estimates but increase computing time; the maximum time can be chosen from known process duration or set slightly above the inverse of the smallest transition rate.Source: Simulation input parameters in MaBoSS · p. 9
• Identifying indecomposable stationary distributions involves repeated simulations that expand initial states or settings until the distributions stabilize.The suggested first run uses a similarity threshold around 0.8 and about 1000 trajectories; thresholds set too high can make clustering inefficient.Source: Simulation input parameters in MaBoSS · pp. 9–10

Comparison with biological data

• Biological data must be discretized, transition rates can be estimated as inverse mean transition times, and probability outputs can represent cell populations and subpopulations.States with non-zero probability in an indecomposable stationary distribution form a subpopulation whose cells evolve among those states.Source: Comparison with biological data · p. 10

Examples

• The examples include a toy single-cycle model and a published p53–Mdm2 response model; BKMC is suited to networks too large for standard state-space computation.For both examples, the authors constructed transition graphs to aid reasoning, although BKMC does not construct them explicitly.Source: Examples · p. 10

Toy model of a single cycle

• At rate up = 10, the toy model’s probability of [ABC]=[000] converges to 1, while entropy and transition entropy converge to zero, marking a fixed point.The entropy peak between times 0 and 0.6 corresponds to transiently activated states before the fixed point.Source: Toy model of a single cycle · pp. 10–11, Figure 3 · p. 12
• The slow-escape toy model is a damped oscillatory process, whereas the large-escape model is not; a transition-graph cycle alone does not ensure oscillation.A small transition entropy alongside entropy that does not converge to zero can indicate a transient cycle.Source: Toy model of a single cycle · pp. 11–12

p53-Mdm2 signaling

• The p53 model has two cycles and a fixed point, and its p53 activation appears transient under the stated initial condition.The initial state is [p53 Mdm2C Mdm2N Dam]=[0*11].Source: p53-Mdm2 signaling · p. 12
• BKMC models the p53–Mdm2 response at the cell-population level, where the results show no damped oscillations unlike the earlier article.The authors attribute this difference to faster loss of synchronization under Markovian noise; damped oscillations were observed with fast p53 activation and slow degradation, in results not shown.Source: p53-Mdm2 signaling · pp. 12–13

Toy model of a single cycle

• At rate down = 10^-5, the apparent cycle is transient: transition entropy is 10^-4, and the longer-time evolution resembles the fast-escape case.The cycle is not stable; the authors conclude stable cyclic behavior only when transition entropy is exactly zero.Source: Toy model of a single cycle · pp. 11–12, Figure 4 · p. 13

p53-Mdm2 signaling

• MaBoSS represents the system as a population and can model different contexts within one model by changing parameters; transient cycles can be detected through entropy trajectories.The authors also identify a flat segment in p53 standard deviation as an indirect effect of transient cycles.Source: p53-Mdm2 signaling · p. 13

Conclusions

• BKMC generalizes asynchronous Boolean dynamics to continuous-time signaling simulations and was applied to a toy cycle and a published p53 DNA-damage model.Source: Conclusions · p. 13
• MaBoSS supports BKMC networks up to 64 nodes and provides global and semi-global dynamics outputs with tunable convergence parameters.Model construction uses a grammar for logical rules and node activation or inhibition rates.Source: Conclusions · pp. 13–14
• In the slow-escape toy model, long-time entropy and transition-entropy behavior is similar to the fast-escape case.Source: Figure 5 · p. 14
• Future work will apply BKMC to larger signaling networks and develop sensitivity analysis, refined approximation methods, and a generalized Markov property.Translating existing Boolean models into MaBoSS requires adding transition rates.Source: Conclusions · pp. 14–15
• The p53 model transition graph contains two cycles and a fixed point, with node colors distinguishing those states.Source: Figure 7 · p. 15

Network state probabilities on time window

• State probabilities for each time window are estimated from each trajectory’s fraction of time spent in each state, then averaged across trajectories.Source: Network state probabilities on time… · p. 16

Entropies

• State entropy summarizes distributional disorder, while transition entropy measures trajectory-level transition determinism; together they provide two-number characterizations regardless of network size.The entropy uses base 2, and zero transition entropy means only one destination state is possible, though transition timing remains stochastic.Source: Entropies · pp. 16–17
• Asymptotically, zero entropy and transition entropy indicate a fixed point, whereas nonzero values indicate a cyclic stationary distribution.Source: Entropies · p. 17

Hamming distance distribution

• The Hamming distance distribution measures how far state probabilities lie from a reference state, giving a semi-global summary sized by node count rather than 2^#nodes.Source: Hamming distance distribution · p. 17

Initial conditions and outputs

• Input nodes are fixed at initialization while other nodes start randomly; internal nodes generate trajectories but are marginalized when probabilities are computed over output nodes.Source: Initial conditions and outputs · pp. 17–18
• When internal nodes are present, transition entropy ignores internal-only flips and computes transition probabilities only for events that flip non-internal nodes.Source: Initial conditions and outputs · p. 18

Stationary distribution characterization

• Finite BKMC trajectories estimate indecomposable stationary distributions by time-averaging state occupancy, then grouping estimates with the same support using a similarity threshold.Source: Stationary distribution characterization · pp. 18–19
• The stationary-distribution clustering procedure is meaningful only for Markov processes, so it uses no internal nodes.Source: Stationary distribution characterization · p. 19

• G. Stoll organized the project, set up algorithms and examples, contributed to the software, and wrote the article.Source: Author’s contributions · p. 19
• E. Viara wrote software and helped set up algorithms; E. Barillot contributed to writing the article.Source: Author’s contributions · p. 19
• L. Calzone organized the project, set up examples, and wrote the article.Source: Author’s contributions · p. 19