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 p <dig> controls cell growth and plays an important role in preventing the development and progression of tumour cells
<cit> .
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
<cit> .
these models can be classified into three broad categories: logical models, continuous models and stochastic models at the single-molecule level
<cit> .
boolean networks  are logical models that utilize discrete state levels and usually assume synchronous and discrete time steps in the evolution of a network
<cit> , whereas continuous models, such as those using linear or ordinary differential equations
<cit> , 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
<cit>  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
<cit> , 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  of experimental data
<cit> .
to account for the intrinsic noise in genetic and molecular interactions, probabilistic boolean networks  have been developed as a generalization of bns
<cit> .
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
<cit> .
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
<cit>  and it is more difficult to compute the steady-state distribution.
this complexity is reduced to o for a sparse state transition matrix
<cit>  and can further be reduced  by ignoring the boolean networks with probabilities below certain threshold
<cit> .
methodologies have also been developed by eliminating genes
<cit>  and using optimal control policies
<cit>  to increase computational efficiency.
state reduction techniques have been used for network intervention
<cit>  and to reduce the model complexity of context-sensitive pbns
<cit> .
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
<cit> .
as an application of bns, logic circuits have been used to simulate genetic networks
<cit> .
recently, circuit diagnosis techniques have been utilized to identify the most vulnerable molecules in cellular networks
<cit> .
synchronous simulation of boolean networks has been proposed for the analysis of biological regulatory networks
<cit> .
an unreliable logic circuit usually behaves probabilistically and thus becomes an instance of pbns.
initially proposed for reliable circuit design
<cit> , stochastic computation has been demonstrated in several physical and biological applications
<cit> .
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 .
it is shown that in an sbn, the complexity to compute the state transition matrix is o, 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
<cit> , 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  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
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
<cit> , a pbn is defined by g , where v = {x <dig>  x <dig>  … xn}, a set of binary-valued nodes, f = , a list of sets of boolean functions:
fi=f1i,f2i,…flii and l is the number of possible functions for gene i,i∈ <dig> n. each node xi represents the state of gene i, denoted by xi and xi =  <dig>  indicates that gene i is  expressed.
the set fi contains the rules that determine the next state of gene i. each
fjii: <dig> n→ <dig> , for
1≤ji≤li, 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 functions in fi, i.e., by
f1i,f2i…flii with probabilities
c1i,c2i…clii.
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  model for a pbn.
a pbn is independent if the functions from fi are independent.
this means that the selection of boolean functions for gene i has no influence on the selection of boolean functions for gene
ji≠j <cit> .
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=∏i=1nli possible bns, each of which is a possible realization of the genetic network.
for the jth bn
1≤j≤n, let
fj=fj <dig> fj22…fjnn, where
1≤ji≤li and i =  <dig>   <dig> … n. the probability that the jth bn is selected is:
pj=∏i=1ncjii,
where
cjii is the probability that the boolean function j 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:
a=p0|0p1|0……p2n−1|0p0|1p1|1……p2n−1|1…………………………p0|2n−1p1|2n−1……p2n−1|2n− <dig>
where each entry is a conditional  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
a=∑j=1npjaj, where pj is the probability that the jth bn occurs and aj is the state transition matrix due to the jth bn.
this way of computing a results in a complexity of o
<cit> .
random gene perturbation, which can occur in an open genome system, is caused by random inputs from outside under external stimuli
<cit> .
by a perturbation, a gene flips its state from  <dig> to  <dig> or vice versa.
since a pbn with perturbation is an aperiodic and irreducible homogeneous markov chain
<cit> , 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
<cit> .
usually, synchronicity is assumed in the state transitions of pbns.
however, a gene-level asynchronous model considers different updating periods of genes in the constituent bns.
in a deterministic-asynchronous boolean network , a gene is assumed to have a fixed updating period
<cit> .
a pbn that uses da-bns as constituent networks is defined as a deterministic-asynchronous probabilistic boolean network .
more rigorously, a da-pbn of n genes consists of a set of
xii=1n, where xi represents the expression level of the ith gene, denoted by xi and
xi∈ <dig>  <cit> .
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 can be used to indicate whether the state of gene i is updated or not, by a value of  <dig> or  <dig> respectively.
the next state of gene i, xi, is then determined by:
xit+1={fjiix1t,…,xntwith probabilitycjii,ifθit+1=1xit,otherwise
where
fjii is a function in the da-bn for gene i, selected with probability
cjii1≤ji≤li.
stochastic boolean networks
<dig>  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
<cit> .
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
<cit> .
figure
<dig> shows an inverter , 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 2-to- <dig> multiplexer of figure
<dig>  for example, its output takes the value of its input ‘a’ or ‘b’ when the control bit ‘c’ is  <dig> or  <dig>  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  <dig> and  <dig> encoded in the random sequences of the control bits.
for a sequence length of  <dig> bits, for example, an input probability of  <dig>  indicates that approximately  <dig> 1’s are in the random sequence of the input ‘a,’ as shown in figure
<dig>  if the random input sequences are independent, the output of the multiplexer is expected to be pa + pbpc =  <dig> , which means that approximately  <dig> 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
<cit> .
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 non-bernoulli sequences.
also see additional files
<dig> and
3: matlab programs that implements the functions of two-input and four-input stochastic multiplexers.
a general structure of the stochastic boolean network  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
<dig>
generally, if a total number of l boolean functions are needed to determine the next state of gene i, an l-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=log2li.
in fact, the number of possible boolean functions for one gene is usually small—between  <dig> and  <dig> for 93% of genes
<cit> .
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 s1 ~ sm, as shown in figure
<dig>  a function in the jth bn is selected with probability
cjii for gene i. when all the genes are accounted for, therefore, an sbn accurately implements the probabilistic functions of a pbn, as dictated by .
<dig>  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 =  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
<cit> :
x′={x⊕γwith a probability of1−1−pnfkxwith a probability of1−pn
where ⊕ is the modulo  <dig> of additions and fk 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
<cit> .
the perturbation matrix is determined by the number of genes and the gene perturbation probability p. it is usually computed by a  analytical approach.
however, the effect of perturbation can be readily accounted for in an sbn.
figure
<dig> 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  <dig> of the perturbation vector and the present state.
the probability that either a boolean function works or a perturbation works ) is computed by a stochastic 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  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:
x′=x⊕γ·1−1−pn+fkx·1−pn,
which is equivalent to .
this indicates that a pbn with perturbation can be accurately implemented by an sbn with perturbation.
<dig>  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  for the so-called deterministic-asynchronous probabilistic boolean networks .
in a da-pbn, 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 da-pbn.
to model this asynchronicity, an sbn can be constructed by considering the timing information as follows:
construct the boolean functions for each gene using the proposed sbn structure.
sort the genes by the updating period and record the sequence.
for example, a sequence can be created as
gt=gt <dig> gt <dig> …,gtn, where the updating periods of
gt <dig> gt <dig> …,gtn are in an ascending order.
consider the current first gene, i.e., the gene with the smallest updating period in gt, denoted by gt.
since the state of gt will first be updated while the states of the other genes remain unchanged, the bns at this stage consist of the boolean functions of gt 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.
delete gt from gt.
repeat steps  and  until gt is empty.
an sbn for a da-pbn is shown in figure
<dig>  since the state transition of a fast-response gene may occur several times before a slow-response gene updates its state, the boolean functions for a fast gene may appear in a number of times in the network of figure
<dig>
applications of sbns
<dig>  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 2n input states and an o 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 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
<cit> .
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, is significantly smaller than the analytical result of o, 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:
construct an sbn by inserting a multiplexer for each gene in a pbn;
for each input state, generate initial random binary streams encoding the control signal probabilities for each multiplexer;
propagate the binary streams from the present state  to the next state  and obtain a random bit sequence for each output;
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;
repeat steps ,  and  for all 2n 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 .
<dig>  estimation of the steady-state distribution
given the size of the state transition matrix of a pbn, the analysis of the steady-state 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 .
an fsm is equivalent to a sequential circuit implementation.
by a time-frame 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
<dig>  this spatial extension of an sbn can be used for the steady-state analysis and the required iterations of the sbn are determined by the number of state transitions before reaching a steady state.
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-mdm <dig> network
in a p <dig> network, signaling pathways are triggered by dna damages and external factors such as chemotherapeutic drugs and ultraviolet light.
for instance, dna double strand breaks  activate pathways that involve the p <dig> and mdm <dig> genes
<cit> .
in response to dsbs, the atm kinase is first stimulated and the chk <dig> is then stimulated by atm.
these activated kinases subsequently induce an increase in the concentration level of p <dig> and a decrease in the interactions between p <dig> and mdm <dig>  the increase in the p <dig> protein level and its transcription activity promote the expression of the mdm <dig> gene, which in turn proceeds to trigger the degradation and destruction of p <dig>  this prior knowledge enables us to come up with the transition rules for the p53-mdm <dig> interactions, as shown in table
<dig>  based on these rules, an independent pbn of the two genes p <dig> and mdm <dig> can be established: v =  with the function classes
f1=f <dig> f <dig> f <dig> f <dig> and
f2=f <dig> f <dig> f <dig> f <dig>  the state transitions of this pbn are given in the truth table of table
<dig>
in table
<dig>  the leftmost column indicates the present state of the genes p <dig> and mdm <dig>  the internal entries in the table indicate whether a function will result in a logical  <dig> or  <dig> 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 ‘ <dig> ’ for example, the next state of the genes can be ‘00’ with a probability of  ×  =  <dig> , ‘01’ with a probability of  ×  =  <dig> , ‘10’ with a probability of  ×  =  <dig>  or ‘11’ with a probability of  ×  =  <dig> .
a pbn is determined by the truth table of table
<dig> and its state transition matrix can be computed as:
apbn= <dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig>
for this pbn, an sbn can be constructed using stochastic multiplexers and random binary bit streams as information carriers, as shown in figure
<dig>  as discussed previously, the control binary sequences determine the probability that each boolean network is selected.
for example, as the boolean functions for the p <dig> gene occur with probabilities  <dig> ,  <dig> ,  <dig>  and  <dig> , the binary bit sequences for the control vectors ‘s1s2’ to the multiplexer are generated with a probability of  <dig>  to be ‘ <dig> ’ a probability of  <dig>  to be ‘01’, a probability of  <dig>  to be ‘10’ and a probability of  <dig>  to be ‘ <dig> ’ then the output bit sequences are read out and decoded into  probabilities.
with a sequence length of  <dig> bits, the state transition matrix is obtained as follows:
asbn= <dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig>
<dig> for the network in figure
<dig>  it has  <dig> genes and  <dig> boolean networks in total.
the difference between  and  is evaluated using the following norms:
.
<dig> and
.∞, which specify the maximum absolute value of the summed differences of columns and rows of the two matrices respectively, and
.
<dig>  which is a measure on the average difference of all the entries in these matrices.
for  and , we obtain
asbn-apbn1= <dig> ,
asbn-apbn2= <dig>  and
asbn-apbn∞= <dig> , 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
<dig>  if the stochastic or outputs a ‘1’ , which means that at least one of the p <dig> and mdm <dig> are perturbed, the multiplexer is then switched to the perturbation network.
if the output of the or is  <dig>  the multiplexer is switched to the original sbn and the network works as the one in figure
<dig> without perturbation.
a similar procedure can be used to compute the state transition matrix of the sbn with perturbation—the result is shown in  for a perturbation probability of  <dig> :
a˜sbn= <dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig>
compared to the analytical result by a method based on :
a˜pbn= <dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig> .
<dig>
the differences between  and  are revealed in the measures of
a˜sbn-a˜pbn1= <dig> ,
a˜sbn-a˜pbn2= <dig>  and
a˜sbn-a˜pbn∞= <dig> .
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  further confirms the p53-mdm <dig> oscillatory dynamics observed in experiments.
an sbn for an asynchronous p53-mdm <dig> network can also be constructed, as in figure
<dig> and following the aforementioned procedure.
due to space limitation, however, this is not further discussed and will be pursued in future work.
RESULTS
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  and a total number of bns .
the simulation is run on a pc with an intel core i3- <dig> cpu  and 6g memory.
the results for using sequence lengths of  <dig> and  <dig> bits are first compared to those obtained using an analytical approach, as shown in table
<dig>  while a larger sequence length of  <dig> bits produces results with a higher precision, a sequence length of  <dig> 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
<dig> and
<dig> 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  <dig> 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
<cit> , the number of boolean functions is considered no larger than  <dig> 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, pseudo-random 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
<dig> and
<dig>  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
<cit>  for several networks with more than  <dig> genes.
the results are shown in table
<dig>
.
the results are obtained from five randomly generated networks, so the standard deviations of the minimum sequence length and run time are also shown.
.
the results are obtained from five randomly generated networks, so the standard deviations of the minimum sequence length and run time are also shown.
.
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 , 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 time-frame 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
<cit> .
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
<dig>  the average runtime for simulating networks of  <dig> and  <dig> genes is shown for various accuracy requirements and perturbation rates.
since the runtime for reaching the steady state is dependent on the initial probabilities , five independent experiments with randomly-selected initial probabilities are performed to obtain an average result.
however, it should be noted that the run time of the time-frame expanded sbn technique is also dependent on the threshold value and perturbation rate.
in table
<dig>  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 20-gene network with a perturbation rate of  <dig>  can be evaluated in approximately  <dig>  seconds using the time-frame expanded sbn technique for a threshold value of  <dig>  .
these results indicate that the time-frame expanded sbn technique is potentially useful in the analysis of large grns.
experiments on a t-cell time series dataset
a network inferred from a time series gene expression dataset
<cit>  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- <dig>  cells were collected at  <dig> different time points before il- <dig> stimulation  and after il- <dig> stimulation .
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
<cit> .
this algorithm is discussed in detail next.
<dig>  inference of boolean dynamics of the grn
pbns have been inferred from steady-state data using the coefficient of determination
<cit>  and from time series data to estimate the perturbation probabilities and switching probabilities between the constituent bns
<cit> .
large amounts of data are usually required by these methods due to their computational complexity.
in
<cit> , 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
<cit> , 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 t-cell time series dataset, a total of  <dig>  networks were discovered by the inference algorithm  <cit> .
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
<cit> .
during this step, the steady states or attractors are computed to validate the inferred networks.
it was found that  <dig>   of these networks did not exhibit the fluctuations expected in the steady-state dynamics of the il- <dig> stimulated t cell network  <cit> .
therefore, these networks were discarded and  <dig>  of the networks that produced biologically meaningful attractors were left for further analysis.
the  <dig> networks were based on twelve meta-genes and yielded a consensus network as shown in figure
<dig>  the steady-state dynamics in the  <dig> networks consist of three time points .
it has also been shown that the computational complexity of this inference algorithm increases exponentially with the maximum number of inputs to a node
<cit> .
however, the maximum input number is limited by the size of a network with a power law
<cit> , so this number is expected to be smaller than  <dig> for a network with less than  <dig> nodes.
the resulting network is not unique in that the occurrence of different boolean functions results in different bns.
in figure
<dig>  the activation and inhibition relationships that occur in all  <dig> 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  <dig> networks in figure
<dig>  this function is considered to occur with a probability of  <dig>  the inferred pbn is shown in the truth tables , for which the boolean functions are assumed to occur independently in a bn.
<dig>  modeling the network with sbn
to build an sbn for the inferred network of figure
<dig>  each of the  <dig> genes is assigned a number, as shown in table
<dig>  for these  <dig> genes, there are  <dig> or  <dig> states, each of which is indexed by the state of each gene as follows:
k=∑i=112gi·2i−1+ <dig>
where i is the gene index and g is the state of gene i .
since solid arrows in figure
<dig> indicate regulatory interactions found in all  <dig> 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-myb <dig> inhibits e-stat5b in all the networks, so the activation of l-myb <dig> 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-bic <dig> and i-myc, as the upregulation of e-stat5b has been observed in the experiments.
an sbn is constructed for the genetic network of figure
<dig>  as shown in figure
<dig>  the construction is based on the following principles:
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-myb <dig> inhibits e-jun-fos, so an inverter is used to simulate this relationship between gt and gt+ <dig>  for the activation of l-foxm <dig> by l-nsbp <dig>  a buffer is applied between gt and gt+ <dig>
an or gate is applied to model multiple activations while a nor  gate is applied to model multiple inhibitions on the same gene.
for example, l-myb <dig> can be activated by any one of e-jun-fos, i-rpol-hnr, e-stat5b and l-mcmd, so in figure
<dig>  gt, gt, gt and gt are used as the four inputs to an or gate.
however, due to the inhibition of l-myb <dig> 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 gt+ <dig>  the use of the and is dictated by the priority rule of the inhibition over the activation of l-myb <dig>  as explained as follows.
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-myb <dig> in the example of , as shown in figure
<dig>
a solid arrow indicates a relationship that exists in all  <dig> networks and therefore is considered to occur with a probability of  <dig>  the corresponding function then exists in every boolean function that produces an input to a mux.
for example, e-stat5a inhibits l-myb <dig> in all the networks, so inverters are present in both of the two boolean functions that lead to gt+ <dig>
<dig>  steady-state evaluation
for this sbn, the state transition matrix at is of the size  <dig> x  <dig> and computed in about 70s.
see additional file
5: the matlab program that describes the structure of the sbn in figure
<dig> and computes its state transition matrix .
given an initial input, i <dig> = , as indicated by the vector at t = 1 h in table
<dig> of
<cit>  that corresponds to the state  <dig> ), the output response after t clock cycles can be computed by:
outputt=i0att
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
<cit> .
although these steady states, or attractors, can be computed using a bn-based method ,  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  <dig> clock cycles.
as shown in table
<dig>  the obtained stationary states with the highest probabilities perfectly match the three attractors found at the time points t <dig>  t <dig> and t <dig> in
<cit> , referred to as attractors  <dig>   <dig> and  <dig> at states  <dig>   <dig> and  <dig>
alternatively, and more efficiently, the aforementioned time-frame expansion technique can be used to estimate the attractors with a greatly reduced complexity.
the results are shown in figure
<dig> for the same sbn simulation of  <dig>   <dig> and 30 cycles and the largest runtime is only  <dig> s, compared to more than 70s by using the matrix-based analysis.
it can be seen that the steady states in figure
<dig> match the attractors in table
<dig>  this shows the effectiveness and efficiency of the time-frame expansion technique.
<dig>  perturbation and prediction
when the genes in a network are perturbed with a small probability, an sbn with perturbation can be constructed  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  <dig> , figure
<dig> shows the steady state distribution of the sbn with perturbation for the network in figure
<dig>
it can be seen that the steady states in figure
<dig> with the highest probabilities  <dig> ,  <dig>  and  <dig>  match the known attractors  <dig>   <dig> and  <dig> .
what is interesting, however, is that pseudo-attractors exist in a perturbed network.
pseudo-attractors 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 pseudo-attractors with a steady state probability equal or larger than  <dig>  are listed in table
<dig>  it can be seen that most of these pseudo-attractors differ from the closest known attractor by only one gene.
in particular, the most prominent pseudo-attractor, located at state  <dig> with a probability larger than  <dig> , differs from attractor  <dig> or state  <dig> by the expression of l-myb <dig>  l-myb <dig> is a late response gene and plays an important role in the regulation of the t-cell network, so this result confirms the sensitivity of l-myb <dig> in the regulatory behaviour.
since biological experiments are not straightforward or easy to be implemented for investigating the t-cell 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.
.
application of the time-frame expansion technique yields similar predictions for the network under perturbation.
for a perturbation rate of  <dig>  and a threshold value of  <dig>  for norm infinity, it only takes  <dig>  seconds to obtain the steady state distribution using a sequence length of  <dig>  bits, in contrast to  <dig>  seconds using the matrix-based sbn method and  <dig>  seconds using the analytical method in
<cit> .
the simulation results are shown in figure
<dig> for the initial state  <dig> , which agree with those in figure
<dig>  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  perturbation.
see additional file
6: 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  <dig> .
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
<cit> .
in figure s <dig> of the additional file
<dig>  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  method.
the time-frame expansion technique is also more efficient compared to the markov chain monte carlo  method.
in table s <dig> of the additional file
<dig>  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  in the computation of state transition matrices and the evaluation of steady state distributions.
relationship to other grn models
<dig>  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:
dg1dtdg2dt⋮dgndt=tg1g2⋮gn,
where gi, , 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 pi, , 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  can be expressed as:
dp1dtdp2dt⋮dpndt=tp1p2⋮pn.
in an sbn, the next state of genes, xt+ <dig>  is determined by the current state, xt, and the state transition matrix, a, i.e.,
xt+1=axt,
where a is a 2n × 2n matrix, as given by .
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
pt+1=gpt,
where pt+ <dig> and pt indicate the gene levels at two consecutive time steps.
further assume that
Δp=pt+1−pt.
in the limit, we obtain:
dpdt=pt+1−ptdt=g−idtpt,
where i is the identity matrix.
finally,  and  lead to
g−i=t·dt,
which describes the relationship between the transition matrices in a continuous model and an sbn.
<dig>  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  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
<cit> .
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
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
<dig>  matlab packages for applications using sbns, including both for the matrix-based analysis and the time-frame expansion technique, are provided as additional files.
CONCLUSIONS
this paper proposes a novel structure of stochastic boolean networks  for an accurate and efficient implementation of probabilistic boolean networks .
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, 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 and o. 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-mdm <dig> 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
<cit> , 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.
supplementary material
additional file 1
stochastic logic using non-bernoulli sequences.
click here for file
additional file 2
mux <dig> m. ‘mux <dig> m’ is a matlab program, which implements the function of a two-input stochastic multiplexer  for an sbn.
click here for file
additional file 3
mux <dig> m. ‘mux <dig> m’ is a matlab program, which implements the function of a four-input stochastic multiplexer  for an sbn.
click here for file
additional file 4
truth table of the pbn inferred from the t cell microarray time series data.
click here for file
additional file 5
t_cell_sbn.m.
‘t_cell_sbn.m’ is a matlab program, which describes the structure of an sbn for the t-cell genetic network and computes its state transition matrix for both without and with perturbation.
the programs ‘mux <dig> m’ and ‘mux <dig> m’ are needed to run ‘t_cell_sbn.m.’
click here for file
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 t-cell genetic network.
the programs ‘mux <dig> m’ and ‘mux <dig> m’ are needed to run ‘time_frame_expansion.m.’
click here for file
acknowledgements
funding: this work was supported by the nserc discovery grant  and the startup fund of the university of alberta .
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.

 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 p controls cell growth and plays an important role in preventing the development and progression of tumour cells. 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. these models can be classified into three broad categories: logical models, continuous models and stochastic models at the single-molecule level. boolean networks are logical models that utilize discrete state levels and usually assume synchronous and discrete time steps in the evolution of a network, whereas continuous models, such as those using linear or ordinary differential equations, 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 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, 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 of experimental data. to account for the intrinsic noise in genetic and molecular interactions, probabilistic boolean networks have been developed as a generalization of bns. 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. 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 and it is more difficult to compute the steady-state distribution. this complexity is reduced to o for a sparse state transition matrix and can further be reduced by ignoring the boolean networks with probabilities below certain threshold. methodologies have also been developed by eliminating genes and using optimal control policies to increase computational efficiency. state reduction techniques have been used for network intervention and to reduce the model complexity of context-sensitive pbns. 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. as an application of bns, logic circuits have been used to simulate genetic networks. recently, circuit diagnosis techniques have been utilized to identify the most vulnerable molecules in cellular networks. synchronous simulation of boolean networks has been proposed for the analysis of biological regulatory networks. an unreliable logic circuit usually behaves probabilistically and thus becomes an instance of pbns. initially proposed for reliable circuit design, stochastic computation has been demonstrated in several physical and biological applications. 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. it is shown that in an sbn, the complexity to compute the state transition matrix is o, 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, 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 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 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, a pbn is defined by g, where v = {x x … xn}, a set of binary-valued nodes, f =, a list of sets of boolean functions: fi=f1i,f2i,…flii and l is the number of possible functions for gene i,i∈ n. each node xi represents the state of gene i, denoted by xi and xi =  indicates that gene i is expressed. the set fi contains the rules that determine the next state of gene i. each fjii: n→, for 1≤ji≤li, 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 functions in fi, i.e., by f1i,f2i…flii with probabilities c1i,c2i…clii. 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 model for a pbn. a pbn is independent if the functions from fi are independent. this means that the selection of boolean functions for gene i has no influence on the selection of boolean functions for gene ji≠j. 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=∏i=1nli possible bns, each of which is a possible realization of the genetic network. for the jth bn 1≤j≤n, let fj=fj fj22…fjnn, where 1≤ji≤li and i =  … n. the probability that the jth bn is selected is: pj=∏i=1ncjii, where cjii is the probability that the boolean function j 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: a=p0|0p1|0……p2n−1|0p0|1p1|1……p2n−1|1…………………………p0|2n−1p1|2n−1……p2n−1|2n− where each entry is a conditional 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 a=∑j=1npjaj, where pj is the probability that the jth bn occurs and aj is the state transition matrix due to the jth bn. this way of computing a results in a complexity of o. random gene perturbation, which can occur in an open genome system, is caused by random inputs from outside under external stimuli. by a perturbation, a gene flips its state from to or vice versa. since a pbn with perturbation is an aperiodic and irreducible homogeneous markov chain, 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. usually, synchronicity is assumed in the state transitions of pbns. however, a gene-level asynchronous model considers different updating periods of genes in the constituent bns. in a deterministic-asynchronous boolean network, a gene is assumed to have a fixed updating period. a pbn that uses da-bns as constituent networks is defined as a deterministic-asynchronous probabilistic boolean network. more rigorously, a da-pbn of n genes consists of a set of xii=1n, where xi represents the expression level of the ith gene, denoted by xi and xi∈. 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 can be used to indicate whether the state of gene i is updated or not, by a value of or respectively. the next state of gene i, xi, is then determined by: xit+1={fjiix1t,…,xntwith probabilitycjii,ifθit+1=1xit,otherwise where fjii is a function in the da-bn for gene i, selected with probability cjii1≤ji≤li. stochastic boolean networks 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. 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. figure  shows an inverter, 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 2-to- multiplexer of figure  for example, its output takes the value of its input ‘a’ or ‘b’ when the control bit ‘c’ is or 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 and encoded in the random sequences of the control bits. for a sequence length of bits, for example, an input probability of indicates that approximately 1’s are in the random sequence of the input ‘a,’ as shown in figure  if the random input sequences are independent, the output of the multiplexer is expected to be pa + pbpc = , which means that approximately 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. 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 non-bernoulli sequences. also see additional files and 3: matlab programs that implements the functions of two-input and four-input stochastic multiplexers. a general structure of the stochastic boolean network 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  generally, if a total number of l boolean functions are needed to determine the next state of gene i, an l-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=log2li. in fact, the number of possible boolean functions for one gene is usually small—between and for 93% of genes. 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 s1 ~ sm, as shown in figure  a function in the jth bn is selected with probability cjii for gene i. when all the genes are accounted for, therefore, an sbn accurately implements the probabilistic functions of a pbn, as dictated by. 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 =  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: x′={x⊕γwith a probability of1−1−pnfkxwith a probability of1−pn where ⊕ is the modulo of additions and fk 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. the perturbation matrix is determined by the number of genes and the gene perturbation probability p. it is usually computed by a analytical approach. however, the effect of perturbation can be readily accounted for in an sbn. figure  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 of the perturbation vector and the present state. the probability that either a boolean function works or a perturbation works ) is computed by a stochastic 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 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: x′=x⊕γ·1−1−pn+fkx·1−pn, which is equivalent to. this indicates that a pbn with perturbation can be accurately implemented by an sbn with perturbation. 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 for the so-called deterministic-asynchronous probabilistic boolean networks. in a da-pbn, 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 da-pbn. to model this asynchronicity, an sbn can be constructed by considering the timing information as follows: construct the boolean functions for each gene using the proposed sbn structure. sort the genes by the updating period and record the sequence. for example, a sequence can be created as gt=gt gt …,gtn, where the updating periods of gt gt …,gtn are in an ascending order. consider the current first gene, i.e., the gene with the smallest updating period in gt, denoted by gt. since the state of gt will first be updated while the states of the other genes remain unchanged, the bns at this stage consist of the boolean functions of gt 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. delete gt from gt. repeat steps and until gt is empty. an sbn for a da-pbn is shown in figure  since the state transition of a fast-response gene may occur several times before a slow-response gene updates its state, the boolean functions for a fast gene may appear in a number of times in the network of figure  applications of sbns 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 2n input states and an o 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 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. 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, is significantly smaller than the analytical result of o, 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: construct an sbn by inserting a multiplexer for each gene in a pbn; for each input state, generate initial random binary streams encoding the control signal probabilities for each multiplexer; propagate the binary streams from the present state to the next state and obtain a random bit sequence for each output; 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; repeat steps, and for all 2n 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. estimation of the steady-state distribution given the size of the state transition matrix of a pbn, the analysis of the steady-state 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. an fsm is equivalent to a sequential circuit implementation. by a time-frame 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  this spatial extension of an sbn can be used for the steady-state analysis and the required iterations of the sbn are determined by the number of state transitions before reaching a steady state. 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-mdm network in a p network, signaling pathways are triggered by dna damages and external factors such as chemotherapeutic drugs and ultraviolet light. for instance, dna double strand breaks activate pathways that involve the p and mdm genes. in response to dsbs, the atm kinase is first stimulated and the chk is then stimulated by atm. these activated kinases subsequently induce an increase in the concentration level of p and a decrease in the interactions between p and mdm the increase in the p protein level and its transcription activity promote the expression of the mdm gene, which in turn proceeds to trigger the degradation and destruction of p this prior knowledge enables us to come up with the transition rules for the p53-mdm interactions, as shown in table  based on these rules, an independent pbn of the two genes p and mdm can be established: v =  with the function classes f1=f f f f and f2=f f f f the state transitions of this pbn are given in the truth table of table  in table  the leftmost column indicates the present state of the genes p and mdm the internal entries in the table indicate whether a function will result in a logical or 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 ‘ ’ for example, the next state of the genes can be ‘00’ with a probability of × =, ‘01’ with a probability of ×  = , ‘10’ with a probability of ×  =  or ‘11’ with a probability of ×  = . a pbn is determined by the truth table of table  and its state transition matrix can be computed as: apbn=........ for this pbn, an sbn can be constructed using stochastic multiplexers and random binary bit streams as information carriers, as shown in figure  as discussed previously, the control binary sequences determine the probability that each boolean network is selected. for example, as the boolean functions for the p gene occur with probabilities,, and, the binary bit sequences for the control vectors ‘s1s2’ to the multiplexer are generated with a probability of to be ‘ ’ a probability of to be ‘01’, a probability of to be ‘10’ and a probability of to be ‘ ’ then the output bit sequences are read out and decoded into probabilities. with a sequence length of bits, the state transition matrix is obtained as follows: asbn=........ for the network in figure  it has genes and boolean networks in total. the difference between and 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 and, we obtain asbn-apbn1=, asbn-apbn2= and asbn-apbn∞=, 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  if the stochastic or outputs a ‘1’, which means that at least one of the p and mdm are perturbed, the multiplexer is then switched to the perturbation network. if the output of the or is the multiplexer is switched to the original sbn and the network works as the one in figure  without perturbation. a similar procedure can be used to compute the state transition matrix of the sbn with perturbation—the result is shown in for a perturbation probability of: a˜sbn=........ compared to the analytical result by a method based on: a˜pbn=........ the differences between and are revealed in the measures of a˜sbn-a˜pbn1=, a˜sbn-a˜pbn2= and a˜sbn-a˜pbn∞=. 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 further confirms the p53-mdm oscillatory dynamics observed in experiments. an sbn for an asynchronous p53-mdm network can also be constructed, as in figure  and following the aforementioned procedure. due to space limitation, however, this is not further discussed and will be pursued in future work.  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 and a total number of bns. the simulation is run on a pc with an intel core i3- cpu and 6g memory. the results for using sequence lengths of and bits are first compared to those obtained using an analytical approach, as shown in table  while a larger sequence length of bits produces results with a higher precision, a sequence length of 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  and 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 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, the number of boolean functions is considered no larger than 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, pseudo-random 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  and 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 for several networks with more than genes. the results are shown in table . the results are obtained from five randomly generated networks, so the standard deviations of the minimum sequence length and run time are also shown.. the results are obtained from five randomly generated networks, so the standard deviations of the minimum sequence length and run time are also shown.. 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, 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 time-frame 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. 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  the average runtime for simulating networks of and genes is shown for various accuracy requirements and perturbation rates. since the runtime for reaching the steady state is dependent on the initial probabilities, five independent experiments with randomly-selected initial probabilities are performed to obtain an average result. however, it should be noted that the run time of the time-frame expanded sbn technique is also dependent on the threshold value and perturbation rate. in table  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 20-gene network with a perturbation rate of can be evaluated in approximately  seconds using the time-frame expanded sbn technique for a threshold value of. these results indicate that the time-frame expanded sbn technique is potentially useful in the analysis of large grns. experiments on a t-cell time series dataset a network inferred from a time series gene expression dataset 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- cells were collected at different time points before il- stimulation and after il- stimulation. 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. this algorithm is discussed in detail next. inference of boolean dynamics of the grn pbns have been inferred from steady-state data using the coefficient of determination and from time series data to estimate the perturbation probabilities and switching probabilities between the constituent bns. large amounts of data are usually required by these methods due to their computational complexity. in, 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, 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 t-cell time series dataset, a total of networks were discovered by the inference algorithm. 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. during this step, the steady states or attractors are computed to validate the inferred networks. it was found that of these networks did not exhibit the fluctuations expected in the steady-state dynamics of the il- stimulated t cell network. therefore, these networks were discarded and of the networks that produced biologically meaningful attractors were left for further analysis. the networks were based on twelve meta-genes and yielded a consensus network as shown in figure  the steady-state dynamics in the networks consist of three time points. it has also been shown that the computational complexity of this inference algorithm increases exponentially with the maximum number of inputs to a node. however, the maximum input number is limited by the size of a network with a power law, so this number is expected to be smaller than for a network with less than nodes. the resulting network is not unique in that the occurrence of different boolean functions results in different bns. in figure  the activation and inhibition relationships that occur in all 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 networks in figure  this function is considered to occur with a probability of the inferred pbn is shown in the truth tables, for which the boolean functions are assumed to occur independently in a bn. modeling the network with sbn to build an sbn for the inferred network of figure  each of the genes is assigned a number, as shown in table  for these genes, there are or states, each of which is indexed by the state of each gene as follows: k=∑i=112gi·2i−1+ where i is the gene index and g is the state of gene i. since solid arrows in figure  indicate regulatory interactions found in all 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-myb inhibits e-stat5b in all the networks, so the activation of l-myb 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-bic and i-myc, as the upregulation of e-stat5b has been observed in the experiments. an sbn is constructed for the genetic network of figure  as shown in figure  the construction is based on the following principles: 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-myb inhibits e-jun-fos, so an inverter is used to simulate this relationship between gt and gt+ for the activation of l-foxm by l-nsbp a buffer is applied between gt and gt+ an or gate is applied to model multiple activations while a nor gate is applied to model multiple inhibitions on the same gene. for example, l-myb can be activated by any one of e-jun-fos, i-rpol-hnr, e-stat5b and l-mcmd, so in figure  gt, gt, gt and gt are used as the four inputs to an or gate. however, due to the inhibition of l-myb 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 gt+ the use of the and is dictated by the priority rule of the inhibition over the activation of l-myb as explained as follows. 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-myb in the example of, as shown in figure  a solid arrow indicates a relationship that exists in all networks and therefore is considered to occur with a probability of the corresponding function then exists in every boolean function that produces an input to a mux. for example, e-stat5a inhibits l-myb in all the networks, so inverters are present in both of the two boolean functions that lead to gt+ steady-state evaluation for this sbn, the state transition matrix at is of the size x and computed in about 70s. see additional file 5: the matlab program that describes the structure of the sbn in figure  and computes its state transition matrix. given an initial input, i =, as indicated by the vector at t = 1 h in table  of that corresponds to the state ), the output response after t clock cycles can be computed by: outputt=i0att 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. although these steady states, or attractors, can be computed using a bn-based method, 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 clock cycles. as shown in table  the obtained stationary states with the highest probabilities perfectly match the three attractors found at the time points t t and t in, referred to as attractors and at states and alternatively, and more efficiently, the aforementioned time-frame expansion technique can be used to estimate the attractors with a greatly reduced complexity. the results are shown in figure  for the same sbn simulation of and 30 cycles and the largest runtime is only s, compared to more than 70s by using the matrix-based analysis. it can be seen that the steady states in figure  match the attractors in table  this shows the effectiveness and efficiency of the time-frame expansion technique. perturbation and prediction when the genes in a network are perturbed with a small probability, an sbn with perturbation can be constructed 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, figure  shows the steady state distribution of the sbn with perturbation for the network in figure  it can be seen that the steady states in figure  with the highest probabilities, and match the known attractors and. what is interesting, however, is that pseudo-attractors exist in a perturbed network. pseudo-attractors 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 pseudo-attractors with a steady state probability equal or larger than are listed in table  it can be seen that most of these pseudo-attractors differ from the closest known attractor by only one gene. in particular, the most prominent pseudo-attractor, located at state with a probability larger than, differs from attractor or state by the expression of l-myb l-myb is a late response gene and plays an important role in the regulation of the t-cell network, so this result confirms the sensitivity of l-myb in the regulatory behaviour. since biological experiments are not straightforward or easy to be implemented for investigating the t-cell 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.. application of the time-frame expansion technique yields similar predictions for the network under perturbation. for a perturbation rate of and a threshold value of for norm infinity, it only takes  seconds to obtain the steady state distribution using a sequence length of bits, in contrast to  seconds using the matrix-based sbn method and  seconds using the analytical method in. the simulation results are shown in figure  for the initial state, which agree with those in figure  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 perturbation. see additional file 6: 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. 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. in figure s of the additional file 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 method. the time-frame expansion technique is also more efficient compared to the markov chain monte carlo method. in table s of the additional file 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 in the computation of state transition matrices and the evaluation of steady state distributions. relationship to other grn models 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: dg1dtdg2dt⋮dgndt=tg1g2⋮gn, where gi,, 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 pi,, 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 can be expressed as: dp1dtdp2dt⋮dpndt=tp1p2⋮pn. in an sbn, the next state of genes, xt+ is determined by the current state, xt, and the state transition matrix, a, i.e., xt+1=axt, where a is a 2n × 2n matrix, as given by. 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 pt+1=gpt, where pt+ and pt indicate the gene levels at two consecutive time steps. further assume that Δp=pt+1−pt. in the limit, we obtain: dpdt=pt+1−ptdt=g−idtpt, where i is the identity matrix. finally, and lead to g−i=t·dt, which describes the relationship between the transition matrices in a continuous model and an sbn. 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 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. 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 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  matlab packages for applications using sbns, including both for the matrix-based analysis and the time-frame expansion technique, are provided as additional files.  this paper proposes a novel structure of stochastic boolean networks for an accurate and efficient implementation of probabilistic boolean networks. 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, 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 and o. 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-mdm 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, 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. supplementary material additional file 1 stochastic logic using non-bernoulli sequences. click here for file additional file 2 mux m. ‘mux m’ is a matlab program, which implements the function of a two-input stochastic multiplexer for an sbn. click here for file additional file 3 mux m. ‘mux m’ is a matlab program, which implements the function of a four-input stochastic multiplexer for an sbn. click here for file additional file 4 truth table of the pbn inferred from the t cell microarray time series data. click here for file additional file 5 t_cell_sbn.m. ‘t_cell_sbn.m’ is a matlab program, which describes the structure of an sbn for the t-cell genetic network and computes its state transition matrix for both without and with perturbation. the programs ‘mux m’ and ‘mux m’ are needed to run ‘t_cell_sbn.m.’ click here for file 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 t-cell genetic network. the programs ‘mux m’ and ‘mux m’ are needed to run ‘time_frame_expansion.m.’ click here for file acknowledgements funding: this work was supported by the nserc discovery grant and the startup fund of the university of alberta. 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.
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. 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. these findings are supported by the simulation results of a simplified p network, several randomly generated networks and a network inferred from a t cell immune response dataset. this paper presents a novel implementation of pbns based on the notions of stochastic logic and stochastic computation. 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. 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. furthermore, a time-frame expanded sbn enables an efficient analysis of the steady-state distribution of a pbn. as a logical model, probabilistic boolean networks consider molecular and genetic noise, so the study of pbns provides significant insights into the understanding of the dynamics of grns. for a pbn with n genes and n boolean networks, the complexity to compute the state transition matrix is o or o for a sparse matrix. 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. stochastic boolean networks are proposed as an efficient approach to modelling gene regulatory networks. the state transition matrix is computed in an sbn with a complexity of o, where l is a factor related to the stochastic sequence length. the algorithms and methods described in this paper have been implemented in matlab packages, which are attached as additional files. various computational models have been of interest due to their use in the modelling of gene regulatory networks. the sbn approach is able to recover biologically-proven regulatory behaviours, such as the oscillatory dynamics of the p53-mdm network and the dynamic attractors in a t cell immune response network. an sbn provides an accurate and efficient simulation of a pbn without and with random gene perturbation. this stochastic implementation of a pbn is referred to as a stochastic boolean network.
 various computational models have been of interest due to their use in the modelling of gene regulatory networks. as a logical model, probabilistic boolean networks 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 or o for a sparse matrix.  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. 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, 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 p 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.  stochastic boolean networks are proposed as an efficient approach to modelling gene regulatory networks. the sbn approach is able to recover biologically-proven regulatory behaviours, such as the oscillatory dynamics of the p53-mdm 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. 