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 steadystate distribution of a PBN. For a PBN with n genes and N Boolean networks, the complexity to compute the state transition matrix is O(nN2^{2n}) or O(nN2^{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(nL2^{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 timeframe expanded SBN enables an efficient analysis of the steadystate 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 singlemolecule 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 biologicallyproven regulatory behaviours, such as the oscillatory dynamics of the p53Mdm2 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.
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 [14]. 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 singlemolecule 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 realvalued 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 highquality experimental data that may not always be available to modellers. As a singlemolecule 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 singlemolecule 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 [1416]. 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 steadystate 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 interventionbased therapeutic approaches [1721].
The computation of the steadystate 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(nN2^{2n}) [15] and it is more difficult to compute the steadystate distribution. This complexity is reduced to O(nN2^{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 contextsensitive 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(nL2^{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 timeframe expanded structure of the SBN, the steadystate 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 singlemolecule 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 binaryvalued 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 binaryvalued nodes, F = (F_{1}, F_{2}, … F_{n}), a list of sets of Boolean functions: and l(i) is the number of possible functions for gene i,. 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 , for , 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 with probabilities . 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 [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 possible BNs, each of which is a possible realization of the genetic network.
For the jth BN , let , where and i = 1, 2 … n. The probability that the jth BN is selected is:
where is the probability that the Boolean function 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 , where P_{j} is the probability that the jth BN occurs and A_{j} is the state transition matrix due to the jth BN. This way of computing A results in a complexity of O(nN2^{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 steadystate distribution is complex [17].
Usually, synchronicity is assumed in the state transitions of PBNs. However, a genelevel asynchronous model considers different updating periods of genes in the constituent BNs. In a deterministicasynchronous Boolean network (DABN), a gene is assumed to have a fixed updating period [16]. A PBN that uses DABNs as constituent networks is defined as a deterministicasynchronous probabilistic Boolean network (DAPBN). More rigorously, a DAPBN of n genes consists of a set of , where X_{i} represents the expression level of the ith gene, denoted by x_{i} and [16]. In a DAPBN, a gene updates its state by its updating period using the DABN 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 is a function in the DABN for gene i, selected with probability .
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 bitwise dependencies encoded in the random binary streams, so making it an efficient approach to computing probabilities [41].
Figure 1 shows an inverter (NOT), an AND, a buffer, an OR, an XOR gate and a multiplexer. While an XOR gate performs a controlled inversion, a multiplexer takes one of its inputs as output according to the values of the control bits. For the 2to1 multiplexer of Figure 1(f), for example, its output takes the value of its input ‘a’ or ‘b’ when the control bit ‘c’ is 0 or 1. Similarly, a stochastic multiplexer chooses one of its inputs as output according to the distributions of 0’s and 1’s and thus the probability of 0 and 1 encoded in the random sequences of the control bits. For a sequence length of 1000 bits, for example, an input probability of 0.4 indicates that approximately 400 1’s are in the random sequence of the input ‘a,’ as shown in Figure 1(f). If the random input sequences are independent, the output of the multiplexer is expected to be P_{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 nonBernoulli 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 file 1: Stochastic Logic using nonBernoulli Sequences. Also see Additional files 2 and 3: Matlab programs that implements the functions of twoinput and fourinput stochastic multiplexers.
Figure 1. Stochastic logic. (a) a NOT gate, (b) an AND gate, (c) a buffer, (d) an OR gate, (e) an XOR gate and (f) a multiplexer. Stochastic logic performs arithmetic operations on the input probabilities encoded in the random binary bit streams. A probability is represented by a proportional number of bits, i.e., the mean number of 1’s in a binary sequence. For illustration, a sequence length of 10 bits is used from (a) to (d); however longer sequences are typically needed in a practical application, as shown in (e) and (f).
Additional file 1. Stochastic Logic using NonBernoulli Sequences.
Format: PDF Size: 141KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 2. mux2.m. ‘mux2.m’ is a Matlab program, which implements the function of a twoinput stochastic multiplexer (MUX, with one control input) for an SBN.
Format: M Size: 1KB Download file
Additional file 3. mux4.m. ‘mux4.m’ is a Matlab program, which implements the function of a fourinput stochastic multiplexer (MUX, with two control inputs) for an SBN.
Format: M Size: 1KB Download file
A general structure of the stochastic Boolean network (SBN) is defined as follows. As shown previously, the next state of genes in a PBN is updated by a set of Boolean functions according to their occurring probabilities. In an SBN, these probabilities are represented by random binary bit sequences and the selection of the Boolean functions is implemented by a stochastic multiplexer with properly generated control sequences. A general structure of an SBN for a single gene is shown in Figure 2.
Figure 2. A stochastic Boolean network (SBN) without perturbation (for a single gene). Only the Boolean functions for a single gene i are shown. The control sequences ‘S_{1} ~ S_{m}’ to the multiplexer (MUX) probabilistically determine the selection of Boolean functions for gene i. When an SBN is constructed for all the genes in a GRN, it accurately implements the probabilistic function of (1).
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 . 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 jth BN is selected with probability 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.
Assume x = (x_{1}, x_{2}, … x_{n}) represents the current state of an ngene 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 kth 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.
However, the effect of perturbation can be readily accounted for in an SBN. Figure 3 illustrates a general model of SBNs with perturbation. As perturbation introduces a probabilistic inversion to the state of a gene, XOR gates are used to implement the addition modulo 2 of the perturbation vector and the present state. The probability that either a Boolean function works or a perturbation works (given in (4)) is computed by a stochastic ninput 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:
which is equivalent to (4). This indicates that a PBN with perturbation can be accurately implemented by an SBN with perturbation.
Figure 3. An SBN with perturbation. A perturbation network is implemented by the stochastic XOR logic of the perturbation vector and the present state. The probability that either a Boolean function works or a perturbation works is given by the output sequence of a stochastic ninput OR gate, which in turn determines the selection of a BN (without perturbation) or a perturbed network by a bus multiplexer.
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 socalled deterministicasynchronous probabilistic Boolean networks (DAPBNs). In a DAPBN, the state of each gene is independently updated according to its own updating period.
While the deterministic asynchronicity changes the temporal sequence of state transitions, it has no impact on the logic relationships among genes, so the Boolean functions are preserved for each gene in a DAPBN. To model this asynchronicity, an SBN can be constructed by considering the timing information as follows:
(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 , where the updating periods of 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.
An SBN for a DAPBN is shown in Figure 4. Since the state transition of a fastresponse gene may occur several times before a slowresponse gene updates its state, the Boolean functions for a fast gene may appear in a number of times in the network of Figure 4.
Figure 4. An SBN for a deterministic asynchronous PBN. Buffers are used to preserve the states of the genes that are not being updated.
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(nL2^{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(nL2^{n}), is significantly smaller than the analytical result of O(nN2^{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.
The procedure of computing the state transition matrix using an SBN is summarized as follows:
(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 steadystate distribution
Given the size of the state transition matrix of a PBN, the analysis of the steadystate distribution is challenging for using both analytical and simulative approaches. The Markovian nature of a PBN makes its analysis similar to that of a finite state machine (FSM). An FSM is equivalent to a sequential circuit implementation. By a timeframe expansion, a sequential circuit can be unrolled into a series of identical combinational modules connected in the spatial domain. Using a similar technique, the temporal operation of an SBN can be transformed into a spatial operation of identical SBNs connected in series. This is shown in Figure 5. This spatial extension of an SBN can be used for the steadystate analysis and the required iterations of the SBN are determined by the number of state transitions before reaching a steady state.
A steadystate analysis using a timeframe 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 steadystate 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 steadystate distribution of a PBN without resorting to the state transition matrix.
Example: the p53Mdm2 network
In a p53 network, signaling pathways are triggered by DNA damages and external factors such as chemotherapeutic drugs and ultraviolet light. For instance, DNA double strand breaks (DSBs) activate pathways that involve the p53 and Mdm2 genes (Figure 6) [3,4]. In response to DSBs, the ATM kinase is first stimulated and the Chk2 is then stimulated by ATM. These activated kinases subsequently induce an increase in the concentration level of p53 and a decrease in the interactions between p53 and Mdm2. The increase in the p53 protein level and its transcription activity promote the expression of the Mdm2 gene, which in turn proceeds to trigger the degradation and destruction of p53. This prior knowledge enables us to come up with the transition rules for the p53Mdm2 interactions, as shown in Table 1. Based on these rules, an independent PBN of the two genes p53 and Mdm2 can be established: V = (X_{1}, X_{2}) with the function classes and . The state transitions of this PBN are given in the truth table of Table 2.
Figure 6. The p53Mdm2 network (adapted from [3]). Under DNA damage, p53 promotes the expression of the Mdm2 gene, which in turn causes the degradation and destruction of p53.
Table 1. State transition probabilities of the p53Mdm2 network
Table 2. Truth table of the PBN for the p53Mdm2 network
In Table 2, the leftmost column indicates the present state of the genes p53 and Mdm2. The internal entries in the table indicate whether a function will result in a logical 1 or 0 at the next state of each gene. The row on the bottom shows the probability of each transition by a function. Given an initial state of ‘01,’ for example, the next state of the genes can be ‘00’ with a probability of (0.09 + 0.01) × (0.5 + 0.4) = 0.09, ‘01’ with a probability of (0.09 + 0.01) × (0.09 + 0.01) = 0.01, ‘10’ with a probability of (0.5 + 0.4) × (0.5 + 0.4) = 0.81 or ‘11’ with a probability of (0.5 + 0.4) × (0.09 + 0.01) = 0.09. A PBN is determined by the truth table of Table 2 and its state transition matrix can be computed as:
For this PBN, an SBN can be constructed using stochastic multiplexers and random binary bit streams as information carriers, as shown in Figure 7. As discussed previously, the control binary sequences determine the probability that each Boolean network is selected. For example, as the Boolean functions for the p53 gene occur with probabilities 0.5, 0.4, 0.09 and 0.01, the binary bit sequences for the control vectors ‘S1S2’ to the multiplexer are generated with a probability of 0.5 to be ‘00,’ a probability of 0.4 to be ‘01’, a probability of 0.09 to be ‘10’ and a probability of 0.01 to be ‘11.’ Then the output bit sequences are read out and decoded into (transition) probabilities. With a sequence length of 10000 bits, the state transition matrix is obtained as follows:
The difference between (6) and (7) is evaluated using the following norms: and , which specify the maximum absolute value of the summed differences of columns and rows of the two matrices respectively, and , which is a measure on the average difference of all the entries in these matrices. For (6) and (7), we obtain , and , which indicate that the SBN structure accurately computes the state transition matrix of the PBN.
With random gene perturbation, an SBN with perturbation can be constructed, as shown in Figure 8. If the stochastic OR outputs a ‘1’ (indicated by S5 in Figure 8), which means that at least one of the p53 and Mdm2 are perturbed, the multiplexer is then switched to the perturbation network. If the output of the OR is 0, the multiplexer is switched to the original SBN and the network works as the one in Figure 7 without perturbation.
Figure 8. An SBN for the p53Mdm2 network (with perturbation).
A similar procedure can be used to compute the state transition matrix of the SBN with perturbation—the result is shown in (8) for a perturbation probability of 0.01:
Compared to the analytical result by a method based on (4):
the differences between (8) and (9) are revealed in the measures of , and . 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 p53Mdm2 oscillatory dynamics observed in experiments.
An SBN for an asynchronous p53Mdm2 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
The state transition matrices of several randomly generated PBNs have been computed using the proposed SBN structure. The Boolean functions of each network are generated for a given number of genes (n) and a total number of BNs (N). The simulation is run on a PC with an Intel Core i32100 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.
In general, a smaller sequence length leads to a shorter run time in the computation of state transition matrices. However, the error incurred due to stochastic fluctuations increases with the size of the network under evaluation. Subsequently, therefore, a minimum accuracy requirement is given and the length of the stochastic sequence is increased for a larger network in order to meet this requirement. Tables 4 and 5 show the minimum sequence lengths and run time required for two different accuracy values, given by the aforementioned “Norm 2” that measures the average difference of all the entries in two matrices. In this experiment, networks of various sizes with up to 12 genes are considered. For each size, five random networks are generated as follows. Given the number of genes in a network, the number of Boolean functions for each gene is initially randomly determined; the specific functions and their associated probabilities are then randomly generated; finally, the input genes are randomly selected for each function. Since a gene’s state is usually determined by no more than four Boolean functions [42], the number of Boolean functions is considered no larger than 4 for each gene. For simplicity, each Boolean function is selected from a set of basic functions: the buffer, NOT, AND, NAND, OR, NOR, XOR and XNOR. In this process, pseudorandom numbers are generated and used in the random selections. For these networks, the standard deviations of the minimum sequence lengths and run time are also shown in Tables 4 and 5. It can be seen that the SBN approach requires a significantly shorter runtime than the analytical approach, especially in the evaluation of large networks. Next, the efficiency of the SBN technique is compared to that of an approximate analytical approach [23] for several networks with more than 10 genes. The results are shown in Table 6.
Table 4. Minimum sequence length and run time required in the computation of state transition matrices for a given accuracy, measured by Norm 2
Table 5. Minimum sequence length and run time required in the computation of state transition matrices for a given accuracy, measured by Norm 2
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.
The state transition matrix computed using an SBN can be used to obtain the steady state distribution of a network. However, the size of the network that can be evaluated is restricted due to the exponential increase of the size of the matrix. As an alternative and efficient approach, the timeframe expansion technique can be used to evaluate much larger networks under perturbation. Recently, several BN models have been developed for GRNs with tens of genes [13,37,38]. Although the parameters for use in a PBN have not been obtained, the time frame expansion technique is well suited for simulating a network of such size, once the necessary parameters become available. In Table 7, the average runtime for simulating networks of 20 and 30 genes is shown for various accuracy requirements and perturbation rates. Since the runtime for reaching the steady state is dependent on the initial probabilities (as applicable in the general Markov chain theory), five independent experiments with randomlyselected initial probabilities are performed to obtain an average result. However, it should be noted that the run time of the timeframe expanded SBN technique is also dependent on the threshold value and perturbation rate. In Table 7, therefore, the average number of periods and run time for convergence, as well as their standard deviations, are shown for several threshold values and perturbation rates. It can be seen, for example, that a 20gene network with a perturbation rate of 0.01 can be evaluated in approximately 2.6 seconds using the timeframe expanded SBN technique for a threshold value of 0.01 (Norm infinity). These results indicate that the timeframe expanded SBN technique is potentially useful in the analysis of large GRNs.
Table 7. Run time of the time frame expansion technique for randomlygenerated networks
Experiments on a Tcell 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 IL2stimulated immune response experiment using a murine T cell line called CTLL2. Cells were collected at 12 different time points before IL2 stimulation (0 h) and after IL2 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 metagenes, 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 steadystate 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.
For the Tcell time series dataset, a total of 161,558 networks were discovered by the inference algorithm [43]. The inference algorithm further explores the dynamics of the inferred networks. This is based on the fact that finite BNs are expected to exhibit a cyclic pattern of expression [7]. During this step, the steady states or attractors are computed to validate the inferred networks. It was found that 160,657 (99.4%) of these networks did not exhibit the fluctuations expected in the steadystate dynamics of the IL2 stimulated T cell network [43]. Therefore, these networks were discarded and 901 (0.6%) of the networks that produced biologically meaningful attractors were left for further analysis. The 901 networks were based on twelve metagenes and yielded a consensus network as shown in Figure 9. The steadystate dynamics in the 901 networks consist of three time points (shown in Table 3 of [43]). It has also been shown that the computational complexity of this inference algorithm increases exponentially with the maximum number of inputs to a node [43]. However, the maximum input number is limited by the size of a network with a power law [46], so this number is expected to be smaller than 5 for a network with less than 100 nodes.
Figure 9. A T cell immune response network inferred from a time series gene expression dataset (adapted from [43]). Solid arrows indicate relationships occurring in all of the 901 networks, while the numbers associated with the dashed arrows indicate the fraction of networks having that relationship. The green lines represent activation relationships and the red lines represent inhibition relationships.
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 file 4: 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.
Additional file 4. Truth Table of the PBN Inferred from the T Cell Microarray Time Series Data.
Format: PDF Size: 75KB Download file
This file can be viewed with: Adobe Acrobat Reader
2. Modeling the network with SBN
To build an SBN for the inferred network of Figure 9, each of the 12 genes is assigned a number, as shown in Table 8. For these 12 genes, there are 2^{12} or 4096 states, each of which is indexed by the state of each gene as follows:
where i is the gene index and g(i) is the state of gene i (i.e., 1 or 0).
Table 8. Code of the 12 genes in the T cell immune response network
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 solidline interaction if they occur simultaneously. For the dashed arrows, the priority is determined according to the observations in the experiments. Take ‘Estat5b’ for example; the solid arrow indicates that LMyb12 inhibits Estat5b in all the networks, so the activation of LMyb12 overrules any other function applied on Estat5b. When the state of Estat5b is only affected by the dashed arrows, the activation by ECdkn2c is considered to take precedence over the inhibitions by IBIc3 and IMyc, as the upregulation of Estat5b has been observed in the experiments.
An SBN is constructed for the genetic network of Figure 9, as shown in Figure 10. The construction is based on the following principles:
(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, LMyb12 inhibits EJunFos, so an inverter is used to simulate this relationship between g_{t}(6) and g_{t+1}(1). For the activation of LFoxm1 by LNsbp1, 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, LMyb12 can be activated by any one of EJunFos, IRpolHnr, Estat5b and LMcmd, 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 LMyb12 by Estat5a, an inverter is applied and its output is ANDed with the output of the 4input 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 LMyb12, 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 LMyb12 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, Estat5a inhibits LMyb12 in all the networks, so inverters are present in both of the two Boolean functions that lead to g_{t+1}(6).
Figure 10. An SBN for the GRN in Figure9.
3. Steadystate evaluation
For this SBN, the state transition matrix A_{T} is of the size 4096 x 4096 and computed in about 70s. See Additional file 5: 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).
Additional file 5. T_cell_SBN.m. ‘T_cell_SBN.m’ is a Matlab program, which describes the structure of an SBN for the Tcell genetic network and computes its state transition matrix for both without and with perturbation. The programs ‘mux2.m’ and ‘mux4.m’ are needed to run ‘T_cell_SBN.m.’
Format: M Size: 5KB Download file
Given an initial input, I_{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 steadystate dynamics consisting of three time points [43]. Although these steady states, or attractors, can be computed using a BNbased method (e.g. [47]), (11) is used here to estimate the attractors as a means to validate the constructed Tcell SBN. In this evaluation, a periodic behaviour of state transitions has been observed after 20 clock cycles.
As shown in Table 9, the obtained stationary states with the highest probabilities perfectly match the three attractors found at the time points t1, t2 and t3 in [43], referred to as Attractors 1, 2 and 3 at states 1224, 1768 and 711.
Alternatively, and more efficiently, the aforementioned timeframe expansion technique can be used to estimate the attractors with a greatly reduced complexity. The results are shown in Figure 11 for the same SBN simulation of 28, 29 and 30 cycles and the largest runtime is only 0.22s, compared to more than 70s by using the matrixbased analysis. It can be seen that the steady states in Figure 11 match the attractors in Table 9. This shows the effectiveness and efficiency of the timeframe expansion technique.
Figure 11. State distributions of the SBN in Figure10after 28, 29 and 30 clock cycles obtained using the timeframe expansion technique. Our simulation shows that the output distribution starts to oscillate after 20 clock cycles.
4. Perturbation and prediction
When the genes in a network are perturbed with a small probability, an SBN with perturbation can be constructed (as in Figure 3) for analyzing the stability of the network under perturbation. Since biological networks are usually robust and stable, the same attractors are often expected to be among the steady states with the highest probabilities for the same network by a small perturbation. Assume that each gene is independently perturbed by a probability 0.01, Figure 12(a) shows the steady state distribution of the SBN with perturbation for the network in Figure 9.
Figure 12. Steady state distribution of the T cell network with perturbation rate of 0.01: (a) computed using state transition matrices and (b) obtained using the time frame expansion technique.
It can be seen that the steady states in Figure 12(a) with the highest probabilities 0.1901, 0.1804 and 0.1750 match the known Attractors 1, 2 and 3 (or, states 1224, 1768 and 711). What is interesting, however, is that pseudoattractors exist in a perturbed network. Pseudoattractors are the steady states with relatively large probabilities due to random gene perturbation, but they are not the attractors in a network without perturbation. The pseudoattractors with a steady state probability equal or larger than 0.01 are listed in Table 10. It can be seen that most of these pseudoattractors differ from the closest known attractor by only one gene. In particular, the most prominent pseudoattractor, located at state 1736 with a probability larger than 0.1, differs from Attractor 2 or state 1768 by the expression of LMyb12. LMyb12 is a late response gene and plays an important role in the regulation of the Tcell network, so this result confirms the sensitivity of LMyb12 in the regulatory behaviour. Since biological experiments are not straightforward or easy to be implemented for investigating the Tcell network under perturbation, such study may provide insights into the understanding of potential physiological implications in a perturbed network. In a long run, this may be helpful in the development of genetic therapeutic methodologies.
Table 10. Pseudoattractors with a steady state probability no smaller than 0.01, as found in the SBN with perturbation
Application of the timeframe 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 matrixbased 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 timeframe expansion technique provides a highly efficient tool for analysing the dynamics of a network with (and without) perturbation. See Additional file 6: The Matlab program that evaluates the steady state distribution using the time frame expansion technique for the Tcell genetic network with a perturbation rate of 0.01.
Additional file 6. time_frame_expansion.m. ‘time_frame_expansion.m’ is a Matlab program, which evaluates the steady state distribution using the time frame expansion technique for the Tcell genetic network. The programs ‘mux2.m’ and ‘mux4.m’ are needed to run ‘time_frame_expansion.m.’
Format: M Size: 5KB Download file
The proposed SBN technique is more efficient than a random sampling approach, due to the use of nonBernoulli 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 file 1, 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 timeframe expansion technique is also more efficient compared to the Markov Chain Monte Carlo (MCMC) method. In Table S1 of the Additional file 1, it is shown that the timeframe 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 pseudorandom 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
Continuous models based on linear or ordinary differential equations can potentially be implemented using SBNs, provided that the underlying principles of the differential equations can be formulated in state transition matrices. In this case, a network of n genes is modelled by:
where 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:
In an SBN, the next state of genes, X_{t+1}, is determined by the current state, X_{t}, and the state transition matrix, A, i.e.,
where A is a 2^{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
where P_{t+1} and P_{t} indicate the gene levels at two consecutive time steps. Further assume that
In the limit, we obtain:
where I is the identity matrix. Finally, (13) and (17) lead to
which describes the relationship between the transition matrices in a continuous model and an SBN.
2. Singlemolecule level models
In a singlemolecule 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 singlemolecule 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 singlemolecule level model, so it may provide an effective means to model large gene regulatory networks.
Application on GRN analysis
In summary, for a GRN inferred from microarray time series data, an SBN can be constructed to analyze the dynamics of the network with or without gene perturbation. This provides the biologists an efficient tool to evaluate the steady state distribution of a genetic network. A general procedure for applying the proposed SBN approach in a GRN analysis is given in the flowchart of Figure 13. Matlab packages for applications using SBNs, including both for the matrixbased analysis and the timeframe expansion technique, are provided as Additional files.
Figure 13. A flowchart for the application of the SBN approach in 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 steadystate analysis of a PBN. The state transition matrix can be accurately and efficiently computed in an SBN with a complexity of O(nL2^{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(nN2^{2n}) and O(nN2^{n}). The steady state distribution can be estimated using the obtained state transition matrix or a timeframe expansion technique. The latter approach has shown a significant speedup in the computation of the steady state distribution.
SBNs have been constructed for the p53Mdm2 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.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
JL and JH conceived the study and participated in its design. JL carried out the GRN studies, performed the statistical analysis and drafted the manuscript. JH participated in the GRN studies and revised the manuscript. All authors read and approved the final manuscript.
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.
References

Weinberg RA: The Biology of Cancer. 1st edition. New York,Garland Science; 2006.

Vogelstein B, Lane D, Levine AJ: Surfing the p53 network.
Nature 2000, 408:307310. PubMed Abstract  Publisher Full Text

Lahav G, Rosenfeld N, Sigal A, GevaZatorsky N, Levine AJ, Elowitz MB, Alon U: Dynamics of the p53Mdm2 feedback loop in individual cells.
Nat Genet 2004, 36:147150. PubMed Abstract  Publisher Full Text

Batchelor E, Loewer A, Lahav G: The ups and downs of p53: understanding protein dynamics in single cells.
Nat Rev Cancer 2009, 9:371377. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

de Jong H: Modeling and simulation of genetic regulatory systems: a literature review.
J Comput Biol January 2002, 9(1):67103. Publisher Full Text

Karlebach G, Shamir R: Modelling and analysis of gene regulatory networks.
Nat Rev Mol Cell Biol 2008, 9:770780. PubMed Abstract  Publisher Full Text

Kauffman SA: Metabolic stability and epigenesis in randomly constructed genetic nets.
J Theor Biol 1969, 22:437467. PubMed Abstract  Publisher Full Text

Klipp E: Systems Biology In Practice: Concepts, Implementation And Application.

Gillespie DT: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.
J Comput Phys 1976, 22:403. Publisher Full Text

Gillespie DT: Exact stochastic simulation of coupled chemical reactions.
J Phys Chem 1977, 81:23402361. Publisher Full Text

Gibson M, Bruck J: Efficient exact stochastic simulation of chemical systems with many species and many channels.

Gillespie DT: Approximate accelerated stochastic simulation of chemically reacting systems.
J Chem Phys 2001, 115:17161733. Publisher Full Text

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 Gprotein action.
Mol Syst Biol 2010., 372 Publisher Full Text

Shmulevich I, Dougherty ER, Zhang W: From Boolean to probabilistic Boolean networks as models of genetic regulatory networks.
Proc IEEE 2002, 90:17781792. Publisher Full Text

Shmulevich I, Dougherty ER, Kim S, Zhang W: Probabilistic Boolean networks: a rulebased uncertainty model for gene regulatory networks.
Bioinformatics 2002, 18:261274. PubMed Abstract  Publisher Full Text

Shmulevich I, Dougherty ER: Probabilistic Boolean Networks: The Modeling and Control of Gene Regulatory Networks. U.S., Society for Industrial & Applied Mathematics; 2010.

Shmulevich I, Dougherty ER, Zhang W: Gene perturbation and intervention in probabilistic Boolean networks.
Bioinformatics 2002, 18(10):13191331. PubMed Abstract  Publisher Full Text

Shmulevich I, Gluhovsky I, Hashimoto RF, Dougherty ER, Zhang W: Steadystate analysis of genetic regulatory networks modelled by probabilistic Boolean networks.
Comp Funct Genom 2003, 4:601608. Publisher Full Text

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):516. Publisher Full Text

Faryabi B, Vahedi G, Datta A, Chamberland JF, Dougherty ER: Recent advances in intervention in Markovian regulatory networks.
Curr Genomics 2009, 10(7):463477. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

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:217240. Publisher Full Text

Ching W, Zhang S, Ng M, Akutsu T: An approximation method for solving the steadystate probability distribution of probabilistic Boolean networks.
Bioinformatics 2007, 23:15111518. PubMed Abstract  Publisher Full Text

Ivanov I, Pal R, Dougherty ER: Dynamics Preserving Size Reduction Mappings for Probabilistic Boolean Networks.

Qian X, Ivanov I, Ghaffari N, Dougherty ER: Intervention in gene regulatory networks via greedy control policies based on longrun behavior.
BMC System Biology 2009, 3:61. BioMed Central Full Text

Qian X, Ghaffari N, Ivanov I, Dougherty ER: State reduction for network intervention in probabilistic Boolean networks.
Bioinformatics 2010, 26(24):30983104. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Pal R: Contextsensitive probabilistic Boolean networks: Steadystate properties, reduction, and steadystate approximation.

Faryabi B, Chamberland JF, Vahedi G, Datta A, Dougherty ER:

Faryabi B, Vahedi G, Chamberland JF, Datta A, Dougherty ER: Intervention in contextsensitive probabilistic Boolean networks revisited.
EURASIP Journal on Bioinformatics Systems Biology 2009, 113.
Special section:

McAdams HH, Shapiro L: Circuit simulation of genetic networks.
Science 1995, 269(5224):650. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text

Kervizic G, Corcos L: Dynamical modeling of the cholesterol regulatory pathway with Boolean networks.
BMC Syst Biol 2008, 2:99. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

von Neumann J: Probabilistic logics and the synthesis of reliable organisms from unreliable components. In Automata Studies. Edited by Shannon CE, McCarthy J. Princeton, Princeton University Press; 1956:4398.

Adar R, Benenson Y, Linshiz G, Rosner A, Tishby N, Shapiro E: Stochastic computing with biomolecular automata.
PNAS 2004, 101(27):99609965. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Benenson Y, Gil B, BenDor U, Adar R, Shapiro E: An autonomous molecular computer for logical control of gene expression.
Nature 2004, 429:423429. PubMed Abstract  Publisher Full Text

Wittmann D, Krumsiek J, SaezRodriguez J, Lauffenburger D, Klamt S, Theis F: Transforming Boolean models to continuous models: methodology and application to Tcell receptor signaling.
BMC System Biology 2009, 3:98. BioMed Central Full Text

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):17771782. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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:814834. Publisher Full Text

Brown B, Card H: Stochastic neural computation I: Computational elements.
IEEE Tran. Computers 2001, 50:891905. Publisher Full Text

Han J, Chen H, Liang J, Zhu P, Yang Z, Lombardi F: A Stochastic Computational Approach for Accurate and Efficient Reliability Evaluation.

Guelzim N, Bottani S, Bourgine P, Kepes F: Topological and causal structure of the yeast transcriptional regulatory network.
Nat Genet 2002, 31:6063. PubMed Abstract  Publisher Full Text

Martin S, Zhang Z, Martino A, Faulon JL: Boolean dynamics of genetic regulatory networks inferred from microarray time series data.
Bioinformatics 2007, 23(7):866874. PubMed Abstract  Publisher Full Text

Marshall S, Yu L, Xiao Y, Dougherty ER: Inference of a probabilistic boolean network from a single observed temporal sequence.

Akutsu T, et al.: Inferring qualitative relations in genetic networks and metabolic pathways.
Bioinformatics 2000, 16:727734. PubMed Abstract  Publisher Full Text

Basso K, et al.: Reverse engineering of regulatory networks in human B cells.
Nat Genet 2005, 37:382390. PubMed Abstract  Publisher Full Text

Garg A, Di Cara A, Xenarios I, Mendoza L, De Micheli G: Synchronous versus asynchronous modeling of gene regulatory networks.
Bioinformatics 2008, 24:19171925. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Layek RK, Datta A, Dougherty ER: From biological pathways to regulatory networks.
Mol Biosyst 2011, 7:843851. PubMed Abstract  Publisher Full Text