- Methodology article
- Open Access

# Stochastic Boolean networks: An efficient approach to modeling gene regulatory networks

- Jinghang Liang
^{1}and - Jie Han
^{1}Email author

**6**:113

https://doi.org/10.1186/1752-0509-6-113

© Liang and Han; licensee BioMed Central Ltd. 2012

**Received:**14 May 2012**Accepted:**6 August 2012**Published:**28 August 2012

## Abstract

### Background

Various computational models have been of interest due to their use in the modelling of gene regulatory networks (GRNs). As a logical model, probabilistic Boolean networks (PBNs) consider molecular and genetic noise, so the study of PBNs provides significant insights into the understanding of the dynamics of GRNs. This will ultimately lead to advances in developing therapeutic methods that intervene in the process of disease development and progression. The applications of PBNs, however, are hindered by the complexities involved in the computation of the state transition matrix and the steady-state distribution of a PBN. For a PBN with *n* genes and *N* Boolean networks, the complexity to compute the state transition matrix is *O*(*nN* 2^{2n}) or *O*(*nN* 2^{
n
}) for a sparse matrix.

### Results

This paper presents a novel implementation of PBNs based on the notions of stochastic logic and stochastic computation. This stochastic implementation of a PBN is referred to as a stochastic Boolean network (SBN). An SBN provides an accurate and efficient simulation of a PBN without and with random gene perturbation. The state transition matrix is computed in an SBN with a complexity of *O*(*nL* 2^{
n
}), where *L* is a factor related to the stochastic sequence length. Since the minimum sequence length required for obtaining an evaluation accuracy approximately increases in a polynomial order with the number of genes, *n*, and the number of Boolean networks, *N*, usually increases exponentially with *n*, *L* is typically smaller than *N*, especially in a network with a large number of genes. Hence, the computational efficiency of an SBN is primarily limited by the number of genes, but not directly by the total possible number of Boolean networks. Furthermore, a time-frame expanded SBN enables an efficient analysis of the steady-state distribution of a PBN. These findings are supported by the simulation results of a simplified p53 network, several randomly generated networks and a network inferred from a T cell immune response dataset. An SBN can also implement the function of an asynchronous PBN and is potentially useful in a hybrid approach in combination with a continuous or single-molecule level stochastic model.

### Conclusions

Stochastic Boolean networks (SBNs) are proposed as an efficient approach to modelling gene regulatory networks (GRNs). The SBN approach is able to recover biologically-proven regulatory behaviours, such as the oscillatory dynamics of the p53-Mdm2 network and the dynamic attractors in a T cell immune response network. The proposed approach can further predict the network dynamics when the genes are under perturbation, thus providing biologically meaningful insights for a better understanding of the dynamics of GRNs. The algorithms and methods described in this paper have been implemented in Matlab packages, which are attached as Additional files.

## Keywords

- Boolean Function
- Boolean Network
- State Transition Matrix
- Stochastic Simulation Algorithm
- Chemical Master Equation

## Background

Biological systems are inherently noisy, yet robust in the presence of noise. The function and malfunction of a system are regulated through the interactions among genes, proteins and other molecules in the cellular network. For instance, the tumour suppressor gene p53 controls cell growth and plays an important role in preventing the development and progression of tumour cells[1–4]. Therefore, it has been of great interest to understand the regulatory mechanisms of genes, and various computational models have been developed for a better understanding of gene regulatory networks (GRNs)[5].

These models can be classified into three broad categories: logical models, continuous models and stochastic models at the single-molecule level[6]. Boolean networks (BNs) are logical models that utilize discrete state levels and usually assume synchronous and discrete time steps in the evolution of a network[7], whereas continuous models, such as those using linear or ordinary differential equations[8], employ real-valued state variables over a continuous timescale. Although continuous models are in principle more accurate and may describe the dynamics of a system in more detail, they require extensive high-quality experimental data that may not always be available to modellers. As a single-molecule level model, Gillespie’s stochastic simulation algorithm (SSA)[9, 10] is based on the chemical master equation; it describes the interactions among single molecules and accounts for noise and stochasticity in the regulation of a genetic network. While the SSA provides the most accurate description of the regulatory behaviour, it requires a large number of parameters and a detailed understanding of the regulatory mechanism. Despite the development of approximate SSAs that trade off accuracy for efficiency[11, 12], algorithms using single-molecule level models are generally slow to run, especially in the modelling of large genetic networks.

Albeit simplistic, BNs have been shown to be efficient in the modelling of GRNs by taking advantages of low complexity and a minimum requirement on the quality (and quantity) of experimental data[13]. To account for the intrinsic noise in genetic and molecular interactions, probabilistic Boolean networks (PBNs) have been developed as a generalization of BNs[14–16]. In a PBN, the inherent stochastic nature of molecular and genetic interactions dictates that the next state of target genes is predicted by several BNs with various probabilities. The evolution of such a system is thus a Markov chain and the state transitions can be described by a transition probability matrix. A steady-state analysis further tells whether a PBN will evolve into a stable target state in the presence of random gene perturbations, thereby providing valuable information for developing intervention-based therapeutic approaches[17–21].

The computation of the steady-state distribution of a PBN, however, presents a challenge. In a PBN with *n* genes and *N* Boolean networks, the complexity to compute the state transition matrix is *O*(*nN* 2^{2n})[15] and it is more difficult to compute the steady-state distribution. This complexity is reduced to *O*(*nN* 2^{
n
}) for a sparse state transition matrix[22] and can further be reduced (to the same order, but with a smaller *N*) by ignoring the Boolean networks with probabilities below certain threshold[23]. Methodologies have also been developed by eliminating genes[24] and using optimal control policies[25] to increase computational efficiency. State reduction techniques have been used for network intervention[26] and to reduce the model complexity of context-sensitive PBNs[27]. Nevertheless, it remains a difficult problem to reduce the computational complexity of a PBN without a compromise on the accuracy of an evaluation.

Although synchronicity is usually assumed in the state transitions of PBNs, asynchronous PBNs have been considered by accounting for different updating periods of genes in the constituent BNs. Asynchronous PBNs are potentially more accurate in describing the regulatory behaviour of genetic networks and may provide a better vehicle for investigating intervention strategies that lead to optimal therapeutic methodologies[28, 29].

As an application of BNs, logic circuits have been used to simulate genetic networks[30]. Recently, circuit diagnosis techniques have been utilized to identify the most vulnerable molecules in cellular networks[31]. Synchronous simulation of Boolean networks has been proposed for the analysis of biological regulatory networks[32]. An unreliable logic circuit usually behaves probabilistically and thus becomes an instance of PBNs. Initially proposed for reliable circuit design[33, 34], stochastic computation has been demonstrated in several physical and biological applications[35, 36].

In this paper, a stochastic computational model is presented for an efficient representation and simulation of PBNs; this implementation of a PBN is referred to as a stochastic Boolean network (SBN). It is shown that in an SBN, the complexity to compute the state transition matrix is *O*(*nL* 2^{
n
}), where *L* is a factor related to the minimum sequence length required for obtaining an evaluation accuracy. In a network with a large number of genes, *L* is usually significantly smaller than *N*. By using a time-frame expanded structure of the SBN, the steady-state distribution can be efficiently computed. Asynchronous PBNs can also be modelled by SBNs for studying the state dynamics of GRNs. With the recent development of BN models[13, 37, 38], the SBN technique is potentially useful in the modelling of large genetic networks. The accuracy and efficiency of the proposed techniques are demonstrated through extensive simulation results. Albeit proposed on the formalism of PBNs, the SBN framework is potentially applicable in improving the simulation efficiency of continuous models (using linear or ordinary differential equations) and single-molecule level models such as those based on SSAs. These aspects are further discussed in the Results and Discussion section.

## Methods

### Probabilistic Boolean networks (PBNs)

In a PBN, genes are represented by a set of binary-valued nodes and the state transitions of genes are described by a list of Boolean functions. Following[15], a PBN is defined by *G (V, F)*, where *V* = {*X*_{
1
}, *X*_{
2
}, … *X*_{
n
}}, a set of binary-valued nodes, *F =* (*F*_{
1
}, *F*_{
2
}, … *F*_{
n
}), a list of sets of Boolean functions:${F}_{i}=\left\{{f}_{1}^{\left(i\right)},{f}_{2}^{\left(i\right)},\dots {f}_{l\left(i\right)}^{\left(i\right)}\right\}$ and *l*(*i*) is the number of possible functions for gene *i,*$i\in \left[1,n\right]$. Each node *X*_{
i
} represents the state of gene *i*, denoted by *x*_{
i
} and *x*_{
i
} = 1 (or 0) indicates that gene *i* is (or not) expressed. The set *F*_{
i
} contains the rules that determine the next state of gene *i*. Each${f}_{j\left(i\right)}^{\left(i\right)}:{\left\{0,1\right\}}^{\text{n}}\to \left\{0,1\right\}$, for$1\le j\left(i\right)\le l\left(i\right)$, is a mapping or a Boolean function determining the state of gene *i*.

Due to the noise in genetic networks, the functions in a PBN occur with certain probabilities. The next state of gene *i* is determined by all the *l(i)* functions in *F*_{
i
}, i.e., by${f}_{1}^{\left(i\right)},{f}_{2}^{\left(i\right)}\dots {f}_{l\left(i\right)}^{\left(i\right)}$ with probabilities${c}_{1}^{\left(i\right)},{c}_{2}^{\left(i\right)}\dots {c}_{l\left(i\right)}^{\left(i\right)}$. Thus, the next state of genes is totally determined by the possible functions and the present state of genes. This indicates that a PBN is modelled as a Markov chain. The fact that all genes are supposed to be updated synchronously also suggests a finite state machine (FSM) model for a PBN.

A PBN is independent if the functions from *F*_{
i
} are independent. This means that the selection of Boolean functions for gene *i* has no influence on the selection of Boolean functions for gene$j\left(i\ne j\right)$[39]. As a first study, the discussions in this paper are limited to independent PBNs. For an independent PBN of *n* genes, there are a total number of$N={\displaystyle \prod _{i=1}^{\text{n}}l\left(i\right)}$ possible BNs, each of which is a possible realization of the genetic network.

*j*th BN$\left(1\le j\le N\right)$, let${f}_{j}=\left[{f}_{j\left(1\right)}^{\left(1\right)},{f}_{j\left(2\right)}^{\left(2\right)}\dots {f}_{j\left(n\right)}^{\left(n\right)}\right]$, where$1\le j\left(i\right)\le l\left(i\right)$ and

*i*= 1, 2 …

*n.*The probability that the

*j*th BN is selected is:

*j*(

*i*) is selected for gene

*i*. By a different selection of the BNs during a state transition, the genes can reach a different state from their present state. This property of a PBN can be described by a state transition matrix as:

where each entry is a conditional (transition) probability that the genes transfer from a given present state into a next state. Since each BN results in a unique next state, the matrix A can be obtained by$\mathit{A}={\displaystyle \sum _{j=1}^{N}{P}_{j}{\mathit{A}}_{j}}$, where *P*_{
j
} is the probability that the *j* th BN occurs and A_{
j
} is the state transition matrix due to the *j* th BN. This way of computing A results in a complexity of *O*(*nN* 2^{2n})[14]. Random gene perturbation, which can occur in an open genome system, is caused by random inputs from outside under external stimuli[17]. By a perturbation, a gene flips its state from 1 to 0 or vice versa. Since a PBN with perturbation is an aperiodic and irreducible homogeneous Markov chain[15], any PBN with perturbation will reach a steady state in a long run. A variant of the state transition matrix A can be used to model the effect of perturbation; however the analysis of its steady-state distribution is complex[17].

*n*genes consists of a set of${\left\{{X}_{i}\right\}}_{i=1}^{n}$, where

*X*

_{ i }represents the expression level of the

*i*th gene, denoted by

*x*

_{ i }and${x}_{i}\in \left\{0,1\right\}$[16]. In a DA-PBN, a gene updates its state by its updating period using the DA-BN involved. At time

*t*, a binary variable

*θ*

_{ i }(t) can be used to indicate whether the state of gene

*i*is updated or not, by a value of 1 or 0 respectively. The next state of gene

*i*,

*x*

_{ i }(t + 1), is then determined by:

where${f}_{j\left(i\right)}^{\left(i\right)}$ is a function in the DA-BN for gene *i*, selected with probability${c}_{j\left(i\right)}^{\left(i\right)}\left(1\le j\left(i\right)\le l\left(i\right)\right)$.

### Stochastic Boolean networks (SBNs)

#### 1. An SBN without perturbation

In stochastic computation, real numbers are represented by random binary bit streams and information is carried in the statistics of the binary streams[40]. A stochastic processing element is typically implemented by a logic gate. Stochastic logic processes information encoded in the random binary bit streams. Probability is represented by a proportional number of bits, usually the mean number of 1’s in a bit sequence. Given independent inputs, for example, an inverter computes the complement of a probability while the multiplication of probabilities is implemented by an AND gate. Hence, stochastic computation transforms Boolean logic operations into probabilistic computation in the real domain. Signal correlations can be efficiently handled in a stochastic network by the bit-wise dependencies encoded in the random binary streams, so making it an efficient approach to computing probabilities[41].

_{a}(1 − P

_{c}) + P

_{b}P

_{c}= 0.34, which means that approximately 340 1’s are expected in the output sequence. Note that this number is only approximate due to the stochastic fluctuations inherent in the representation of the random binary bit streams. This is an important feature in stochastic computation as probabilistic values are propagated rather than deterministic ones, which results in inevitable random fluctuations in the representation of probabilities. It has been shown, however, when non-Bernoulli sequences of random permutations of fixed numbers of 1’s and 0’s are used for representing initial probabilities, these fluctuations are significantly smaller than using Bernoulli sequences, which is equivalent to a random sampling based simulation[41]. It is shown later in the Result section that these fluctuations are generally negligible when reasonably long random bit sequences are used. See Additional file1: Stochastic Logic using non-Bernoulli Sequences. Also see Additional files2 and3: Matlab programs that implements the functions of two-input and four-input stochastic multiplexers.

Generally, if a total number of *l(i)* Boolean functions are needed to determine the next state of gene *i*, an *l(i)*-input multiplexer is used to simulate the selection of functions in the PBN for gene *i*. The number of control bits is given by$m=\u2308{log}_{2}\left(l\left(i\right)\right)\u2309$. In fact, the number of possible Boolean functions for one gene is usually small—between 1 and 4 for 93% of genes[23, 42]. This indicates that one or two bits are usually sufficient to control a multiplexer in an SBN. By using a stochastic multiplexer with the control bit streams S_{1} ~ S_{m}, as shown in Figure 2, a function in the *j* th BN is selected with probability${c}_{j\left(i\right)}^{\left(i\right)}$ for gene *i*. When all the genes are accounted for, therefore, an SBN accurately implements the probabilistic functions of a PBN, as dictated by (1).

#### 2. An SBN with perturbation

While a switch of Boolean functions may indicate a structural change in the network, a random perturbation could cause a transient change of a gene’s state under external stimuli. In a PBN with perturbation, a gene may change its value with a small probability *p* during each state transition.

*x*

_{1},

*x*

_{2}, …

*x*

_{ n }) represents the current state of an

*n*-gene network at time

*t*and γ is the vector that indicates the effect of random perturbation, the next state x

^{′}is given by[17]:

where ⊕ is the modulo 2 of additions and *f*_{
k
}(·) represents the function of the *k* th Boolean network at time *t.* The effect of perturbation to the state transition matrix can then be described by a matrix called the perturbation matrix[23]. The perturbation matrix is determined by the number of genes and the gene perturbation probability *p*. It is usually computed by a (complex) analytical approach*.*

*n*-input OR gate. This probability is then encoded into the output sequence of the OR gate and used as the control sequence of a bus multiplexer. If the perturbation vectors (‘Pert 1’ … ‘Pert n’ in Figure 3) are all 0’s, which means there is no perturbation, then the output sequence of the OR gate contains all 0’s, which subsequently determines that the next state is given by the original SBN without perturbation; otherwise, the next state is determined by the perturbation probability encoded in the output sequence of the stochastic OR gate. Per the stochastic functions of XOR, OR and the multiplexer, the next state is given as the output of the SBN with perturbation, by:

#### 3. An SBN for asynchronous PBNs

In contrast to synchronous PBNs, each gene in an asynchronous PBN has a different period of updating time. Mathematically, this is described by (3) for the so-called deterministic-asynchronous probabilistic Boolean networks (DA-PBNs). In a DA-PBN, the state of each gene is independently updated according to its own updating period.

- (1)
Construct the Boolean functions for each gene using the proposed SBN structure.

- (2)
Sort the genes by the updating period and record the sequence. For example, a sequence can be created as ${G}_{t}=\left\{{{g}_{t}}^{\left(1\right)},{{g}_{t}}^{\left(2\right)},\dots ,{{g}_{t}}^{\left(n\right)}\right\}$, where the updating periods of ${{g}_{t}}^{\left(1\right)},{{g}_{t}}^{\left(2\right)},\dots ,{{g}_{t}}^{\left(n\right)}$ are in an ascending order.

- (3)
Consider the current first gene, i.e., the gene with the smallest updating period in

*G*_{ t }, denoted by*g*_{ t }^{(i)}. Since the state of*g*_{ t }^{(i)}will first be updated while the states of the other genes remain unchanged, the BNs at this stage consist of the Boolean functions of*g*_{ t }^{(i)}and buffers for the other genes. A buffer is a logic element with a delayed input as its output. In this structure, a buffer is used to preserve the state of a gene that is not being updated. - (4)
Delete

*g*_{ t }^{(i)}from*G*_{ t }. - (5)
Repeat steps (3) and (4) until

*G*_{ t }is empty.

### Applications of SBNs

#### 1. Computation of the state transition matrix

In an SBN, each input combination yields output sequences that contain information about the transition probability from this input state to an output state. Therefore, the statistics, i.e., the proportions of the number of each state encoded in the output sequences return the transition probabilities in a row in the state transition matrix. This row corresponds to the given input state and thus all the transition probabilities from this input can be generated in a single run. For a PBN with *n* genes, the SBN needs to be run for each of the 2^{
n
} input states and an *O*(*n*) number of sequences need to be generated for the control signals of the multiplexers.

The accuracy in the computed state transition probabilities is determined by the sequence length of the random binary bit streams. In general, longer sequences are required in a larger network for achieving an evaluation accuracy. To consider the overhead incurred in the use of a larger sequence length, a factor, *L*, is introduced and therefore, a complexity of *O*(*nL* 2^{
n
}) results for computing all the entries in the state transition matrix for a desired accuracy.

It has been shown that the required sequence length is related to the reliability and thus the size of a combinational network[41]. In an SBN, the network size is typically on a polynomial order of the number of genes. This is in contrast with the number of BNs, *N*, which generally increases exponentially with the number of genes. As a result, the complexity of using an SBN to compute the transition matrix, i.e., *O*(*nL* 2^{
n
}), is significantly smaller than the analytical result of *O*(*nN* 2^{
n
}), especially for a network with a large number of genes. This is demonstrated later by simulations using several measures to determine the minimum sequence length required for certain accuracy.

- (1)
Construct an SBN by inserting a multiplexer for each gene in a PBN;

- (2)
For each input state, generate initial random binary streams encoding the control signal probabilities for each multiplexer;

- (3)
Propagate the binary streams from the present state (inputs) to the next state (outputs) and obtain a random bit sequence for each output;

- (4)
Obtain the statistics, i.e., the proportions of the number of each state encoded in the output sequences as the transition probabilities for this input state;

- (5)
Repeat steps (2), (3) and (4) for all 2

^{n}input states to compute all the entries in the state transition matrix.

For an SBN with perturbation, the state transition matrix can be similarly computed using the procedure outlined above with an exception in the construction of the SBN in step (1).

#### 2. Estimation of the steady-state distribution

A steady-state analysis using a time-frame expanded SBN starts with an initial input state, generates the random bit sequences for the inputs and control bits of multiplexers, and then propagates the stochastic signals through the expanded SBN structure. This process is equivalent to an analytical procedure of multiplying the input probabilities with the powers of the state transition matrix. Finally, a small variance threshold is used to determine whether the system has reached a steady state. The steady-state distribution is then obtained from the output sequences at the end of the operation.

In the above process, the speed of convergence to a steady state is dependent on a number of factors, including the length of random bit sequences, the variance threshold value and the perturbation rate. In practice, a sequence length that is long enough to have a resolution of at least two magnitudes smaller than the threshold value, is used to guarantee that the convergence is not dominated by stochastic fluctuations. It is shown later that the analysis using an extended SBN structure provides an alternative and efficient way of estimating the steady-state distribution of a PBN without resorting to the state transition matrix.

### Example: the p53-Mdm2 network

*X*

_{1},

*X*

_{2}) with the function classes${\text{F}}_{1}=\left\{{f}_{1}^{\left(1\right)},{f}_{2}^{\left(1\right)},{f}_{3}^{\left(1\right)},{f}_{4}^{\left(1\right)}\right\}$ and${\text{F}}_{2}=\left\{{f}_{1}^{\left(2\right)},{f}_{2}^{\left(2\right)},{f}_{3}^{\left(2\right)},{f}_{4}^{\left(2\right)}\right\}$. The state transitions of this PBN are given in the truth table of Table 2.

**State transition probabilities of the p53-Mdm2 network**

Present State | Next State Probability | |||
---|---|---|---|---|

p53, Mdm2 | p53 | Mdm2 | ||

(or, x | 0 | 1 | 0 | 1 |

00 | 0.01 | 0.99 | 0.99 | 0.01 |

01 | 0.1 | 0.9 | 0.9 | 0.1 |

10 | 0.9 | 0.1 | 0.1 | 0.9 |

11 | 0.5 | 0.5 | 0.5 | 0.5 |

**Truth table of the PBN for the p53-Mdm2 network**

x | ${\mathit{f}}_{\mathit{1}}^{\left(1\right)}$ | ${f}_{2}^{\left(1\right)}$ | ${f}_{3}^{\left(1\right)}$ | ${f}_{4}^{\left(1\right)}$ | ${f}_{1}^{\left(2\right)}$ | ${f}_{2}^{\left(2\right)}$ | ${f}_{3}^{\left(2\right)}$ | ${f}_{4}^{\left(2\right)}$ |
---|---|---|---|---|---|---|---|---|

00 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 |

01 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 |

10 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 |

11 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 0 |

${c}_{j}^{\left(i\right)}$ | 0.5 | 0.4 | 0.09 | 0.01 | 0.5 | 0.4 | 0.09 | 0.01 |

The difference between (6) and (7) is evaluated using the following norms:${\Vert .\Vert}_{1}$ and${\Vert .\Vert}_{\infty}$, which specify the maximum absolute value of the summed differences of columns and rows of the two matrices respectively, and${\Vert .\Vert}_{2}$, which is a measure on the average difference of all the entries in these matrices. For (6) and (7), we obtain$\Vert {\mathit{A}}_{\mathit{SBN}},-,{\mathit{A}}_{\mathit{PBN}}\Vert {}_{1}=0.0018$,$\Vert {\mathit{A}}_{\mathit{SBN}},-,{\mathit{A}}_{\mathit{PBN}}\Vert {}_{2}=0.0024$ and$\Vert {\mathit{A}}_{\mathit{SBN}},-,{\mathit{A}}_{\mathit{PBN}}\Vert {\infty}_{}=0.0044$, which indicate that the SBN structure accurately computes the state transition matrix of the PBN.

the differences between (8) and (9) are revealed in the measures of$\Vert {\tilde{\mathit{A}}}_{\mathit{SBN}},-,{\tilde{\mathit{A}}}_{\mathit{PBN}}\Vert {}_{1}=0.0032$,${\Vert {\tilde{\mathit{A}}}_{\mathit{SBN}},-,{\tilde{\mathit{A}}}_{\mathit{PBN}}\Vert}_{2}=0.0030$ and$\Vert {\tilde{\mathit{A}}}_{\mathit{SBN}},-,{\tilde{\mathit{A}}}_{\mathit{PBN}}\Vert {}_{\infty}=0.0042$. This shows that the proposed approach using an SBN can accurately and efficiently compute the state transition matrix. The differences in these results come from the stochastic fluctuation, which is an intrinsic property of stochastic computation. More simulation results are presented in the Results and Discussion section, which show that the fluctuations are generally small. A steady state analysis using (8) further confirms the p53-Mdm2 oscillatory dynamics observed in experiments.

An SBN for an asynchronous p53-Mdm2 network can also be constructed, as in Figure 4 and following the aforementioned procedure. Due to space limitation, however, this is not further discussed and will be pursued in future work.

## Results and discussion

### Simulations with randomly generated networks

*n*) and a total number of BNs (

*N*). The simulation is run on a PC with an Intel Core i3-2100 CPU (@3.10 GHz) and 6G memory. The results for using sequence lengths of 10000 and 1000 bits are first compared to those obtained using an analytical approach, as shown in Table 3. While a larger sequence length of 10000 bits produces results with a higher precision, a sequence length of 1000 bits also provides highly accurate results for networks of such size.

**Errors in the state transition matrices obtained using SBNs without perturbation, compared to the results by using the analytical approach in**[22]

Number of genes (n) | 2 | 3 | 4 | 5 | 6 | |
---|---|---|---|---|---|---|

Error | Length (bits) | |||||

${\text{Error}}_{{\Vert \xb7\Vert}_{1}}$ | 1000 | 0.0070 | 0.0330 | 0.0420 | 0.0477 | 0.0649 |

10000 | 0.0027 | 0.0052 | 0.0105 | 0.0179 | 0.0186 | |

${\text{Error}}_{{\Vert \xb7\Vert}_{2}}$ | 1000 | 0.0100 | 0.0314 | 0.0408 | 0.0287 | 0.0405 |

10000 | 0.0038 | 0.0047 | 0.0102 | 0.0109 | 0.0099 | |

${\text{Error}}_{{\Vert \xb7\Vert}_{\infty}}$ | 1000 | 0.0160 | 0.0640 | 0.0908 | 0.0735 | 0.1293 |

10000 | 0.0056 | 0.0096 | 0.0248 | 0.0303 | 0.0248 |

**Minimum sequence length and run time required in the computation of state transition matrices for a given accuracy, measured by Norm 2**

n | N | SBN (Norm 2 = 0.04) | SBN (Norm 2 = 0.02) | Method[22] | |||||||
---|---|---|---|---|---|---|---|---|---|---|---|

Sequence length | Std. deviation | Avg. time (s) | Std. deviation | Sequence length | Std. deviation | Avg. time (s) | Std. deviation | Avg. time (s) | Std. deviation | ||

2 | 6 | 150 | 46 | 0.006324 | 0.003315 | 480 | 84 | 0.013655 | 0.007568 | 0.005468 | 0.004100 |

3 | 8 | 460 | 89 | 0.019755 | 0.008942 | 800 | 122 | 0.017634 | 0.009536 | 0.011655 | 0.007036 |

4 | 16 | 520 | 109 | 0.024337 | 0.009108 | 1120 | 84 | 0.043844 | 0.010102 | 0.031391 | 0.009388 |

5 | 32 | 860 | 134 | 0.052112 | 0.017356 | 1540 | 182 | 0.118927 | 0.036943 | 0.157794 | 0.020922 |

6 | 64 | 1240 | 270 | 0.209416 | 0.030298 | 2460 | 241 | 0.548156 | 0.042366 | 0.532971 | 0.037483 |

7 | 128 | 1340 | 167 | 0.453192 | 0.048960 | 3680 | 239 | 1.208252 | 0.060325 | 2.441066 | 0.163347 |

8 | 256 | 2260 | 378 | 2.030217 | 0.171125 | 5480 | 335 | 4.110083 | 0.326308 | 9.368184 | 0.863544 |

9 | 512 | 2580 | 303 | 4.751360 | 0.421918 | 6820 | 471 | 12.81050 | 2.061854 | 39. 26049 | 4.208466 |

10 | 1024 | 3920 | 923 | 16.06112 | 4.252810 | 8760 | 1135 | 38.60258 | 6.377620 | 201.5433 | 10.90932 |

11 | 2048 | 4700 | 836 | 40.44380 | 5.742303 | 10400 | 1140 | 95.40610 | 7.547263 | 811.6358 | 15.88395 |

12 | 4096 | 5660 | 882 | 118.3426 | 9.031772 | 13000 | 1000 | 286.5043 | 12.37633 | 3501.744 | 86.66141 |

n | N | SBN (Norm 2 = 0.04) | SBN (Norm 2 = 0.02) | Method[22] | |||||||
---|---|---|---|---|---|---|---|---|---|---|---|

Sequence length | Std. deviation | Avg. time (s) | Std. deviation | Sequence length | Std. deviation | Avg. time (s) | Std. deviation | Avg. time (s) | Std. deviation | ||

2 | 6 | 180 | 45 | 0.008052 | 0.005219 | 340 | 55 | 0.017285 | 0.010683 | 0.050477 | 0.010140 |

3 | 8 | 460 | 114 | 0.020473 | 0.011034 | 920 | 130 | 0.027358 | 0.010944 | 0.026389 | 0.014326 |

4 | 16 | 660 | 152 | 0.032089 | 0.023041 | 1220 | 148 | 0.055602 | 0.022138 | 0.053726 | 0.021034 |

5 | 32 | 880 | 130 | 0.071256 | 0.020862 | 1620 | 130 | 0.162794 | 0.047719 | 0.161462 | 0.039981 |

6 | 64 | 1320 | 228 | 0.235628 | 0.038845 | 2460 | 288 | 0.443522 | 0.056302 | 0.613840 | 0.047252 |

7 | 128 | 1480 | 130 | 0.574352 | 0.062129 | 4240 | 261 | 1.540875 | 0.071316 | 2.663523 | 0.180211 |

8 | 256 | 2420 | 319 | 2.124709 | 0.228612 | 5620 | 319 | 4.411751 | 0.413352 | 11.90834 | 1.412206 |

9 | 512 | 3220 | 650 | 7.248265 | 2.301722 | 6940 | 498 | 14.36077 | 3.253704 | 61.45203 | 6.881528 |

10 | 1024 | 4140 | 882 | 18.09032 | 4.112405 | 9400 | 1140 | 41.41356 | 5.289815 | 261.3189 | 12.29343 |

11 | 2048 | 4860 | 606 | 47.37403 | 5.822926 | 11800 | 837 | 120.3839 | 6.107372 | 975.3821 | 33.25207 |

12 | 4096 | 5820 | 782 | 132.8137 | 10.90686 | 14600 | 1342 | 373.7601 | 13.64551 | 4022.140 | 78.42531 |

**Run time and errors in the computation of state transition matrices for SBN and the approximate method in**[23]

n | N | SBN (s) (Length = 10000 bits) | Method[23] (s) (lower bound = 10 | ${\mathbf{\text{Error}}}_{{\Vert \xb7\Vert}_{\mathbf{2}}}$ (SBN) | ${\mathbf{\text{Error}}}_{{\Vert \xb7\Vert}_{\mathbf{2}}}$[23] | ||||
---|---|---|---|---|---|---|---|---|---|

${\mathbf{\text{Error}}}_{{\Vert \xb7\Vert}_{\mathbf{1}}}$ | ${\mathbf{\text{Error}}}_{{\Vert \xb7\Vert}_{\mathbf{2}}}$ | ${\mathbf{\text{Error}}}_{{\Vert \xb7\Vert}_{\infty}}$ | ${\mathbf{\text{Error}}}_{{\Vert \xb7\Vert}_{\mathbf{1}}}$ | ${\mathbf{\text{Error}}}_{{\Vert \xb7\Vert}_{\mathbf{2}}}$ | ${\mathbf{\text{Error}}}_{{\Vert \xb7\Vert}_{\infty}}$ | ||||

11 | 2048 | 92.367577 | 183.617225 | 0.2031 | 0.0268 | 0.1209 | 0.2416 | 0.0463 | 0.0221 |

12 | 4096 | 221.849183 | 1125.969347 | 0.3448 | 0.0301 | 0.1540 | 0.6387 | 0.0929 | 0.0386 |

13 | 8192 | 489.265478 | 4395.954714 | 0.4581 | 0.0552 | 0.2249 | 1.6583 | 0.1414 | 0.0874 |

14 | 16384 | 1063.892415 | 9415.812415 | 1.0152 | 0.0825 | 0.4287 | 2.1642 | 0.2283 | 0.1895 |

As revealed in the tables, while an analytical approach is fast in computing the state transition matrices of small networks, it becomes cumbersome to use for larger networks. This is because an analytical approach is limited by the number of BNs (*N*), which generally increases exponentially with the number of genes in a PBN. In an SBN, however, all the state transition probabilities for each input state are encoded in the output sequences, so the computation of the state transition matrix is very efficient. Although a longer stochastic sequence length is required to meet a higher accuracy, the proposed SBN approach still outperforms an analytical approach for networks with a large number of genes and BNs, because its efficiency is not directly limited by the number of BNs.

**Run time of the time frame expansion technique for randomly-generated networks**

Number of genes | Sequence length (bits) | Threshold value (Norm infinity) | Perturbation rate | SBN (results from five experiments with different initial values) | |||
---|---|---|---|---|---|---|---|

Average number of periods for convergence | Standard deviation of the number of periods | Run time (s) | Standard deviation of the run time (s) | ||||

20 | 100,000 | 0.001 | 0.0001 | 1648 | 278.3 | 1477.8 | 125.1 |

100,000 | 0.001 | 0.001 | 202 | 21.3 | 182.65 | 18.62 | |

10,000 | 0.01 | 0.01 | 20 | 4.9 | 2.5878 | 0.5663 | |

30 | 1,000,000 | 0.01 | 0.0001 | 1128 | 173.2 | 15536 | 3121 |

1,000,000 | 0.1 | 0.001 | 66 | 26.2 | 904.04 | 395.6 |

### Experiments on a T-cell time series dataset

A network inferred from a time series gene expression dataset[43] is further modelled using SBNs. The dataset was taken from an IL-2-stimulated immune response experiment using a murine T cell line called CTLL-2. Cells were collected at 12 different time points before IL-2 stimulation (0 h) and after IL-2 stimulation (15, 30 mins, 1, 2, 4, 6, 8, 10, 12, 16 and 24 h). The dataset was then normalized to the same expression level and clustered based on the similarities in the regulatory behaviour of the genes. This produced simplified networks of gene groups, referred to as meta-genes, instead of actual genes. This result has significantly reduced the complexity of the analysis and interpretation of the inferred networks. Finally, the dataset was discretized for the implementation of a Boolean network inference algorithm[43]. This algorithm is discussed in detail next.

#### 1. Inference of Boolean dynamics of the GRN

PBNs have been inferred from steady-state data using the coefficient of determination[15] and from time series data to estimate the perturbation probabilities and switching probabilities between the constituent BNs[44]. Large amounts of data are usually required by these methods due to their computational complexity. In[43], the Boolean inference is based on the activation and inhibition functions of a target gene and its control genes. This is similar to the qualitative inference method used in[45], but it considers all possible networks rather than a single most likely one. While the number of possible inputs to a Boolean function is limited in this method, the restriction on the amount of data required to perform an inference is released. The number of possible networks is then counted and all networks are enumerated.

The resulting network is not unique in that the occurrence of different Boolean functions results in different BNs. In Figure 9, the activation and inhibition relationships that occur in all 901 networks are indicated by solid arrows, while the relationships that occur in a fraction of the networks are indicated by dashed arrows. The value associated with a dashed arrow indicates the fraction of networks having that relationship. To infer a PBN, this fractional occurrence of a function is considered probabilistic and its associated value is taken as the occurrence probability of a Boolean function in the network. These probabilities are then utilized to obtain the switching probabilities between the constituent BNs in the PBN. Since a solid arrow indicates a relationship that exists in all 901 networks in Figure 9, this function is considered to occur with a probability of 1. The inferred PBN is shown in the truth tables (see Additional file4: Truth table of the PBN inferred from the T cell microarray time series data), for which the Boolean functions are assumed to occur independently in a BN.

#### 2. Modeling the network with SBN

^{12}or 4096 states, each of which is indexed by the state of each gene as follows:

*i*is the gene index and

*g(i)*is the state of gene

*i*(i.e., 1 or 0).

**Code of the 12 genes in the T cell immune response network**

Gene | E-Jun-Fos | L-Nsbp1 | L-Foxm1 | I-BIc3 | I-Myc | L-Myb12 | E-Cdkn2c | E-Stat1-6 | I-Rpol-hnr | E-stat5a | E-stat5b | L-Mcmd |
---|---|---|---|---|---|---|---|---|---|---|---|---|

Symbol | g(1) | g(2) | g(3) | g(4) | g(5) | g(6) | g(7) | g(8) | g(9) | g(10) | g(11) | g(12) |

Since solid arrows in Figure 9 indicate regulatory interactions found in all 901 networks, they are considered to have a priority over other interactions, i.e., any other relationships are overruled by a solid-line interaction if they occur simultaneously. For the dashed arrows, the priority is determined according to the observations in the experiments. Take ‘E-stat5b’ for example; the solid arrow indicates that L-Myb12 inhibits E-stat5b in all the networks, so the activation of L-Myb12 overrules any other function applied on E-stat5b. When the state of E-stat5b is only affected by the dashed arrows, the activation by E-Cdkn2c is considered to take precedence over the inhibitions by I-BIc3 and I-Myc, as the upregulation of E-stat5b has been observed in the experiments.

- (1)
An inhibited signal is considered logical “low” while an activated signal is considered logical “high.” Therefore, an inverter or a buffer is applied to represent an inhibition or an activation relationship between genes. For example, L-Myb12 inhibits E-Jun-Fos, so an inverter is used to simulate this relationship between

*g*_{ t }(6) and*g*_{ t+1 }(1). For the activation of L-Foxm1 by L-Nsbp1, a buffer is applied between*g*_{ t }(2) and*g*_{ t+1 }(3). - (2)
An OR gate is applied to model multiple activations while a NOR (inverted OR) gate is applied to model multiple inhibitions on the same gene. For example, L-Myb12 can be activated by any one of E-Jun-Fos, I-Rpol-Hnr, E-stat5b and L-Mcmd, so in Figure 10,

*g*_{ t }(1),*g*_{ t }(9),*g*_{ t }(11) and*g*_{ t }(12) are used as the four inputs to an OR gate. However, due to the inhibition of L-Myb12 by E-stat5a, an inverter is applied and its output is ANDed with the output of the 4-input OR gate to produce the output of*g*_{ t+1 }(6). The use of the AND is dictated by the priority rule of the inhibition over the activation of L-Myb12, as explained as follows. - (3)
When an inhibition and activation occur on the same gene, the logic gate is determined by the priority of the two functions: an AND gate is applied if the inhibition has a higher priority, whereas an OR gate is used if the activation has a higher priority. For instance, an AND gate is used to model the relationship between the activation and inhibition of L-Myb12 in the example of (2), as shown in Figure 10.

- (4)
A solid arrow indicates a relationship that exists in all 901 networks and therefore is considered to occur with a probability of 1. The corresponding function then exists in every Boolean function that produces an input to a MUX. For example, E-stat5a inhibits L-Myb12 in all the networks, so inverters are present in both of the two Boolean functions that lead to

*g*_{ t+1 }(6).

#### 3. Steady-state evaluation

For this SBN, the state transition matrix A_{
T
} is of the size 4096 x 4096 and computed in about 70s. See Additional file5: The Matlab program that describes the structure of the SBN in Figure 10 and computes its state transition matrix (for both without and with perturbation).

_{ 0 }= [0, 0… 0, 1, 0… 0], as indicated by the vector at

*T = 1 h*in Table 3 of[43] that corresponds to the state 1730 (by (10)), the output response after

*t*clock cycles can be computed by:

A clock cycle here corresponds to the time interval between two discrete time points as a period of biological response. It has been shown that the network exhibits a steady-state dynamics consisting of three time points[43]. Although these steady states, or attractors, can be computed using a BN-based method (e.g.[47]), (11) is used here to estimate the attractors as a means to validate the constructed T-cell SBN. In this evaluation, a periodic behaviour of state transitions has been observed after 20 clock cycles.

*t1*,

*t2*and

*t3*in[43], referred to as Attractors 1, 2 and 3 at states 1224, 1768 and 711.

**Attractors found by the SBN approach, compared to the experimental results in**[43]

Number of cycles | States with highest probabilities | Attractors found in[43] |
---|---|---|

28 | 1224 | Attractor 1 |

711 | Attractor 3 | |

1768 | Attractor 2 | |

29 | 1768 | Attractor 2 |

1224 | Attractor 1 | |

711 | Attractor 3 | |

30 | 711 | Attractor 3 |

1768 | Attractor 2 | |

1224 | Attractor 1 |

#### 4. Perturbation and prediction

**Pseudo-attractors with a steady state probability no smaller than 0.01, as found in the SBN with perturbation**

State Number | Probability | Closest attractor | Difference |
---|---|---|---|

1736 | 0.1099 | Attractor 2 | g(6) (L-Myb12) |

967 | 0.0203 | Attractor 3 | g(9) (I-Rpol-hnr) |

199 | 0.0164 | Attractor 3 | g(10) (E-stat5a) |

3816 | 0.0147 | Attractor 2 | g(12) (L-Mcmd) |

3866 | 0.0135 | Different from all the attractors by more than 3 genes | |

743 | 0.0120 | Attractor 3 | g(6) (L-Myb12) |

1352 | 0.0101 | Attractor 1 | g(8) (E-Stat1-6) |

g(9) (I-Rpol-hnr) | |||

1256 | 0.0100 | Attractor 1 | g(6) (L-Myb12) |

Application of the time-frame expansion technique yields similar predictions for the network under perturbation. For a perturbation rate of 0.01 and a threshold value of 0.01 for Norm infinity, it only takes 3.7 seconds to obtain the steady state distribution using a sequence length of 10,000 bits, in contrast to 212.1 seconds using the matrix-based SBN method and 2532.9 seconds using the analytical method in[22]. The simulation results are shown in Figure 12(b) for the initial state 1730 (as considered in[43]), which agree with those in Figure 12(a). As the speed of convergence of the time frame expansion technique is dependent on the initial state of the network, several different initial states have been randomly selected and all of them have resulted in a runtime less than 100 seconds. Therefore, the time-frame expansion technique provides a highly efficient tool for analysing the dynamics of a network with (and without) perturbation. See Additional file6: The Matlab program that evaluates the steady state distribution using the time frame expansion technique for the T-cell genetic network with a perturbation rate of 0.01.

The proposed SBN technique is more efficient than a random sampling approach, due to the use of non-Bernoulli sequences of random permutations of fixed numbers of 1’s and 0’s in the representation of initial probabilities[41]. In Figure S3 of the Additional file1, it is shown that smaller variations generally result in the state transition matrices computed using the SBN technique compared to those obtained using the Monte Carlo (MC) method. The time-frame expansion technique is also more efficient compared to the Markov Chain Monte Carlo (MCMC) method. In Table S1 of the Additional file1, it is shown that the time-frame expanded SBN technique converges faster to a steady state than the MCMC method, because it requires a fewer number of clock cycles or time frames to converge and generates less pseudo-random numbers at each time frame. These indicate that the proposed SBN approach is more accurate and more efficient than a simple random sampling approach (such as the MC simulation) in the computation of state transition matrices and the evaluation of steady state distributions.

### Relationship to other GRN models

#### 1. Continuous models

*n*genes is modelled by:

*g*

_{ i }, (

*i =*1, 2, …

*n*), indicates the level of a gene and T is a matrix of

*n*rows and

*n*columns. The entries in T are determined by factors such as the reaction rate constants. If the gene level can be expressed as the occurrence rate of a gene, denoted by

*p*

_{ i }, (

*i =*1, 2, …

*n*), which, for example, can be obtained by the ratio between the number of a particular type of genes and the total number of genes, then (12) can be expressed as:

_{ t +1 }, is determined by the current state, X

_{ t }, and the state transition matrix, A, i.e.,

^{n}× 2

^{n}matrix, as given by (2). Then a new transition matrix of

*n*rows and

*n*columns, denoted by

*G*, can be obtained by summarizing the entries in the rows and columns of A, such that

_{ t +1 }and P

_{ t }indicate the gene levels at two consecutive time steps. Further assume that

which describes the relationship between the transition matrices in a continuous model and an SBN.

#### 2. Single-molecule level models

In a single-molecule level model, significant stochastic effects of biochemical reactions are accounted for each molecular species. The stochastic simulation algorithm (SSA) tracks the number of molecular species in a biochemical system, so it accurately simulates the discrete, random biochemical reactions specified by the chemical master equation (CME)[9, 10]. Essentially, the SSA follows a discrete Markov process, in which two values are generated from two independent random variables at each time step. The first value predicts when the next reaction will occur and the second decides which reaction will occur. In order to characterize the evolution of the system, repeated trials are required to perform, which leads to a significant run time for simulating a large network.

Due to the same underlying Markov models in the SSA and PBNs, the SSA can, in principle, be implemented using SBNs. However, this implementation is not straightforward as the SSA simulates the function of the CME while the SBN implements the state transitions of Boolean functions. A challenge is therefore to formulate the underlying principles of the CME in the form of state transition matrices. Nevertheless, it is possible for the SSA and SBN to be used in a hybrid method. In this method, a logical model is first used to simulate a large network and to identify the sensitive nodes in the network. Then, a single-molecule level model such as the SSA can be used to find out more details of the identified sensitive genes. In this way, this hybrid method leverages the efficiency of a logical model and the accuracy of a single-molecule level model, so it may provide an effective means to model large gene regulatory networks.

### Application on GRN analysis

## Conclusions

This paper proposes a novel structure of stochastic Boolean networks (SBNs) for an accurate and efficient implementation of probabilistic Boolean networks (PBNs). The application of an SBN is demonstrated through the computation of the state transition matrix and the steady-state analysis of a PBN. The state transition matrix can be accurately and efficiently computed in an SBN with a complexity of *O*(*nL* 2^{
n
}), where *n* is the number of genes in a PBN and *L* is a factor determined by the stochastic sequence length. Since the required minimum sequence length for a given evaluation accuracy usually increases slower with *n* than the number of Boolean networks, i.e., *N*, *L* is typically smaller than *N*, especially in a network with a large number of genes. This result is an improvement compared to the previous results of *O*(*nN* 2^{2n}) and *O*(*nN* 2^{
n
}). The steady state distribution can be estimated using the obtained state transition matrix or a time-frame expansion technique. The latter approach has shown a significant speedup in the computation of the steady state distribution.

SBNs have been constructed for the p53-Mdm2 network and an inferred T cell immune response network. Simulations of the SBNs have recovered state dynamics that have been experimentally demonstrated for these two networks. The proposed approach is able to discover network dynamics when the genes are under perturbation, which is a difficult task to implement in experiments or by other modeling approaches due to its complexity. So in this case, the SBN technique can be used to provide biologically meaningful insights for a first understanding of the dynamics of a GRN. The relationship between an SBN and continuous/stochastic models has also been discussed and a hybrid approach may be useful in a more efficient modelling of a large GRN. Finally, the SBN approach is able to account for signalling pathway information[48], so it may provide an effective solution to the modeling of complex genetic networks.

## Declarations

### Acknowledgements

*Funding*: This work was supported by the NSERC Discovery Grant (JL) and the Startup fund of the University of Alberta (JH). The authors would like to thank Lukasz Kurgan of the University of Alberta for his helpful advices during the preparation of the manuscript and the reviewers for their constructive and insightful comments.

## Authors’ Affiliations

## References

- Weinberg RA: The Biology of Cancer. 2006, New York,Garland Science, 1Google Scholar
- Vogelstein B, Lane D, Levine AJ: Surfing the p53 network. Nature. 2000, 408: 307-310. 10.1038/35042675.View ArticleGoogle Scholar
- Lahav G, Rosenfeld N, Sigal A, Geva-Zatorsky N, Levine AJ, Elowitz MB, Alon U: Dynamics of the p53-Mdm2 feedback loop in individual cells. Nat Genet. 2004, 36: 147-150. 10.1038/ng1293.View ArticleGoogle Scholar
- Batchelor E, Loewer A, Lahav G: The ups and downs of p53: understanding protein dynamics in single cells. Nat Rev Cancer. 2009, 9: 371-377. 10.1038/nrc2604.View ArticleGoogle Scholar
- de Jong H: Modeling and simulation of genetic regulatory systems: a literature review. J Comput Biol. January 2002, 9 (1): 67-103. 10.1089/10665270252833208.View ArticleGoogle Scholar
- Karlebach G, Shamir R: Modelling and analysis of gene regulatory networks. Nat Rev Mol Cell Biol. 2008, 9: 770-780. 10.1038/nrm2503.View ArticleGoogle Scholar
- Kauffman SA: Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol. 1969, 22: 437-467. 10.1016/0022-5193(69)90015-0.View ArticleGoogle Scholar
- Klipp E: Systems Biology In Practice: Concepts, Implementation And Application. Wiley-VCH, Weinheim. 2005Google Scholar
- Gillespie DT: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comput Phys. 1976, 22: 403-10.1016/0021-9991(76)90041-3.View ArticleGoogle Scholar
- Gillespie DT: Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977, 81: 2340-2361. 10.1021/j100540a008.View ArticleGoogle Scholar
- Gibson M, Bruck J: Efficient exact stochastic simulation of chemical systems with many species and many channels. J Phys Chem. 1999, 104: 1876-1889.View ArticleGoogle Scholar
- Gillespie DT: Approximate accelerated stochastic simulation of chemically reacting systems. J Chem Phys. 2001, 115: 1716-1733. 10.1063/1.1378322.View ArticleGoogle Scholar
- Pandey S, Wang R, Wilson L, Li S, Zhao Z, Gookin T, Assmann S, Albert R: Boolean modeling of transcriptome data reveals novel modes of heterotrimeric G-protein action. Mol Syst Biol. 2010, 372: 10.1038/msb.2010.28.Google Scholar
- Shmulevich I, Dougherty ER, Zhang W: From Boolean to probabilistic Boolean networks as models of genetic regulatory networks. Proc IEEE. 2002, 90: 1778-1792. 10.1109/JPROC.2002.804686.View ArticleGoogle Scholar
- Shmulevich I, Dougherty ER, Kim S, Zhang W: Probabilistic Boolean networks: a rule-based uncertainty model for gene regulatory networks. Bioinformatics. 2002, 18: 261-274. 10.1093/bioinformatics/18.2.261.View ArticleGoogle Scholar
- Shmulevich I, Dougherty ER: Probabilistic Boolean Networks: The Modeling and Control of Gene Regulatory Networks. 2010, U.S., Society for Industrial & Applied MathematicsView ArticleGoogle Scholar
- Shmulevich I, Dougherty ER, Zhang W: Gene perturbation and intervention in probabilistic Boolean networks. Bioinformatics. 2002, 18 (10): 1319-1331. 10.1093/bioinformatics/18.10.1319.View ArticleGoogle Scholar
- Shmulevich I, Gluhovsky I, Hashimoto RF, Dougherty ER, Zhang W: Steady-state analysis of genetic regulatory networks modelled by probabilistic Boolean networks. Comp Funct Genom. 2003, 4: 601-608. 10.1002/cfg.342.View ArticleGoogle Scholar
- Dougherty ER, Pal R, Qian X, Bittner ML, Datta A: Stationary and Structural Control in Gene Regulatory Networks: Basic Concepts. Int J Syst Sci. 2010, 41 (1): 5-16. 10.1080/00207720903144560.View ArticleGoogle Scholar
- Faryabi B, Vahedi G, Datta A, Chamberland JF, Dougherty ER: Recent advances in intervention in Markovian regulatory networks. Curr Genomics. 2009, 10 (7): 463-477. 10.2174/138920209789208246.View ArticleGoogle Scholar
- Karlebach G, Shamir R: Minimally perturbing a gene regulatory network to avoid a disease phenotype: the glioma network as a test case. BMC Syst Biol. 2010, 4: 15-10.1186/1752-0509-4-15.View ArticleGoogle Scholar
- Zhang S, Ching W, M Ng, Akutsu T: Simulation study in probabilistic Boolean network models for genetic regulatory networks. Int J Data Min. 2007, 1: 217-240. 10.1504/IJDMB.2007.011610.View ArticleGoogle Scholar
- Ching W, Zhang S, Ng M, Akutsu T: An approximation method for solving the steady-state probability distribution of probabilistic Boolean networks. Bioinformatics. 2007, 23: 1511-1518. 10.1093/bioinformatics/btm142.View ArticleGoogle Scholar
- Ivanov I, Pal R, Dougherty ER: Dynamics Preserving Size Reduction Mappings for Probabilistic Boolean Networks. IEEE Trans Signal Process. 2007, 55 (5): 2310-2322.View ArticleGoogle Scholar
- Qian X, Ivanov I, Ghaffari N, Dougherty ER: Intervention in gene regulatory networks via greedy control policies based on long-run behavior. BMC System Biology. 2009, 3: 61-10.1186/1752-0509-3-61.View ArticleGoogle Scholar
- Qian X, Ghaffari N, Ivanov I, Dougherty ER: State reduction for network intervention in probabilistic Boolean networks. Bioinformatics. 2010, 26 (24): 3098-3104. 10.1093/bioinformatics/btq575.View ArticleGoogle Scholar
- Pal R: Context-sensitive probabilistic Boolean networks: Steady-state properties, reduction, and steady-state approximation. IEEE T Signal Proces. 2010, 58 (2): 879-890.View ArticleGoogle Scholar
- Faryabi B, Chamberland J-F, Vahedi G, Datta A, Dougherty ER: IEEE J Sel Top Signa. 2008, 2 (3): 412-423.View ArticleGoogle Scholar
- Faryabi B, Vahedi G, Chamberland JF, Datta A, Dougherty ER: Intervention in context-sensitive probabilistic Boolean networks revisited. EURASIP Journal on Bioinformatics Systems Biology. 2009, 1-13. Special section:Google Scholar
- McAdams HH, Shapiro L: Circuit simulation of genetic networks. Science. 1995, 269 (5224): 650-10.1126/science.7624793.View ArticleGoogle Scholar
- Abdi A, Tahoori MB, Emamian ES: Fault diagnosis engineering of digital circuits can identify vulnerable molecules in complex cellular pathways. Sci Signal. 2008, 1 (42): ra10-10.1126/scisignal.2000008.View ArticleGoogle Scholar
- Kervizic G, Corcos L: Dynamical modeling of the cholesterol regulatory pathway with Boolean networks. BMC Syst Biol. 2008, 2: 99-10.1186/1752-0509-2-99.View ArticleGoogle Scholar
- von Neumann J: Probabilistic logics and the synthesis of reliable organisms from unreliable components. Automata Studies. Edited by: Shannon CE, McCarthy J. 1956, Princeton, Princeton University Press, 43-98.Google Scholar
- Gaines BR: Stochastic Computing Systems. Adv Inf Syst Sci. 1969, 2: 37-172.View ArticleGoogle Scholar
- Adar R, Benenson Y, Linshiz G, Rosner A, Tishby N, Shapiro E: Stochastic computing with biomolecular automata. PNAS. 2004, 101 (27): 9960-9965. 10.1073/pnas.0400731101.View ArticleGoogle Scholar
- Benenson Y, Gil B, Ben-Dor U, Adar R, Shapiro E: An autonomous molecular computer for logical control of gene expression. Nature. 2004, 429: 423-429. 10.1038/nature02551.View ArticleGoogle Scholar
- Wittmann D, Krumsiek J, Saez-Rodriguez J, Lauffenburger D, Klamt S, Theis F: Transforming Boolean models to continuous models: methodology and application to T-cell receptor signaling. BMC System Biology. 2009, 3: 98-10.1186/1752-0509-3-98.View ArticleGoogle Scholar
- Duarte N, Becker S, Jamshidi N, Thiele I, Mo M, Vo T, Srivas R, Palsson B: Global reconstruction of the human metabolic network based on genomic and bibliomic data. PNAS. 2007, 104 (6): 1777-1782. 10.1073/pnas.0610772104.View ArticleGoogle Scholar
- Lähdesmäki H, et al.: Relationships between probabilistic Boolean networks and dynamic Bayesian networks as models of gene regulatory networks. Signal Process. 2006, 86: 814-834. 10.1016/j.sigpro.2005.06.008.View ArticleGoogle Scholar
- Brown B, Card H: Stochastic neural computation I: Computational elements. IEEE Tran. Computers. 2001, 50: 891-905. 10.1109/12.954505.View ArticleGoogle Scholar
- Han J, Chen H, Liang J, Zhu P, Yang Z, Lombardi F: A Stochastic Computational Approach for Accurate and Efficient Reliability Evaluation. IEEE Trans Comput. 2013Google Scholar
- Guelzim N, Bottani S, Bourgine P, Kepes F: Topological and causal structure of the yeast transcriptional regulatory network. Nat Genet. 2002, 31: 60-63. 10.1038/ng873.View ArticleGoogle Scholar
- Martin S, Zhang Z, Martino A, Faulon J-L: Boolean dynamics of genetic regulatory networks inferred from microarray time series data. Bioinformatics. 2007, 23 (7): 866-874. 10.1093/bioinformatics/btm021.View ArticleGoogle Scholar
- Marshall S, Yu L, Xiao Y, Dougherty ER: Inference of a probabilistic boolean network from a single observed temporal sequence. EURASIP J Bioinforma Syst Biol. 2007, 2007: 32454-Google Scholar
- Akutsu T, et al.: Inferring qualitative relations in genetic networks and metabolic pathways. Bioinformatics. 2000, 16: 727-734. 10.1093/bioinformatics/16.8.727.View ArticleGoogle Scholar
- Basso K, et al.: Reverse engineering of regulatory networks in human B cells. Nat Genet. 2005, 37: 382-390. 10.1038/ng1532.View ArticleGoogle Scholar
- Garg A, Di Cara A, Xenarios I, Mendoza L, De Micheli G: Synchronous versus asynchronous modeling of gene regulatory networks. Bioinformatics. 2008, 24: 1917-1925. 10.1093/bioinformatics/btn336.View ArticleGoogle Scholar
- Layek RK, Datta A, Dougherty ER: From biological pathways to regulatory networks. Mol Biosyst. 2011, 7: 843-851. 10.1039/c0mb00263a.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.

## Comments

View archived comments (1)