BACKGROUND
diverse biological functions are regulated through the interactions among genes, proteins and other molecules in a cell.
gene expressions are however affected by the intrinsic and extrinsic noise in a gene network  <cit> .
a major source of the noise is the stochastic fluctuations in gene regulatory interactions  <cit> .
the genetic interactions are also context dependent, that is, certain regulatory functions are active in some cellular states, but inactive in others  <cit> .
this indicates the necessity to consider noise in a context-sensitive manner in the study of gene regulatory networks .
various methods have been proposed to model grns; these include logical models  <cit> , continuous models using differential equations  <cit>  and stochastic models at the single-molecule level  <cit> .
as a classic logical model, boolean networks  have been widely used to qualitatively model the interactions among genes  <cit> .
probabilistic boolean networks  have been proposed to consider noise in a bn model  <cit> .
in a pbn, the next state of a gene is determined by its current state and a boolean function.
if the boolean function is randomly selected, a pbn is referred to as an instantaneous pbn  <cit> .
as a general model, a context-sensitive pbn  considers the feature of context dependence in a biological network  <cit> .
in a cspbn, a context is a combination of boolean functions and each function determines the next state of a gene.
a context remains unchanged until a switching occurs.
this switching of contexts, possibly caused by external stimuli, is considered to occur randomly in a network.
the study of pbns has focused on the analysis of steady state distributions  under gene perturbation and intervention.
a markov chain monte carlo method is used in  <cit>  to analyze the long run behavior of a pbn; however, this method generally requires a large number of simulations to reach a steady state, due to the slow convergence typically encountered in a monte carlo method  <cit> .
an analysis is performed in  <cit>  for finding the ssd of a pbn through the computation of the state transition matrix .
however, the application of an analytical approach is generally limited to small networks due to the exponential increase in the size of an stm with gene numbers.
the analysis of cspbns presents even a greater challenge due to its significantly increased computational complexity.
analytical expressions have been derived for analyzing cspbns  <cit> ; however, this method is only applicable to the steady state analysis of a network with small perturbation and switching probabilities.
the method in  <cit>  ignores some bns with very small probabilities for reducing the size of the stm and thus provides a more efficient but approximate solution for computing the ssd of a cspbn.
in  <cit> , gene intervention is investigated for avoiding undesirable states associated with certain diseases .
due to external stimuli, the stm is changed by external control variables, so desirable states can be obtained with larger probabilities in the ssd.
in a context-sensitive network with n genes and k contexts, however, a  ×   <cit>   matrix is required for an accurate  analysis of the ssd; this results in a computational complexity of o ) for an accurate  computation of the stm.
hence, the application of current gene network analysis is limited to those of less than a dozen of genes.
recently, stochastic boolean networks  have been proposed for an efficient computation of the stm and ssd of an instantaneous pbn  <cit> .
the sbn approach can recover biologically-proven regulatory behaviors, such as the oscillatory dynamics of a simplified p53-mdm <dig> network  <cit>  and the dynamic attractors in a t cell immune response network  <cit> .
to further exploit the simplicity of logical models, stochastic multiple-valued networks  have been utilized to investigate the dynamics of gene networks with multiple-valued gene states  <cit> .
asynchronous sbns have also been developed to investigate the asynchronous state-update behavior of genes  <cit> .
in this paper, the notion of sbns is extended for the general model of pbns, i.e., the cspbns.
this presents a novel structure of context-sensitive sbns  for modeling the stochasticity in a context-dependent grn.
in particular, gene perturbation and intervention are considered in cssbns.
through an efficient simulation of a cssbn, the computational complexity in the evaluation of a cspbn is reduced from o to o for computing the stm, where l is a factor related to the required sequence length in cssbn for achieving a desired accuracy.
the use of non-bernoulli sequences of random permutations of fixed number of 1’s and 0’s further increases computational efficiency and allows for a tunable tradeoff between accuracy and efficiency.
a time-frame expanded cssbn can further simulate the stationary behavior of a cspbn and produce results that are more accurate than an approximate analysis, thus making the cssbn useful in modeling complex context-sensitive grns.
the cssbn models are applied to the study of a simplified p53-mdm <dig> network and the results using an ssd control policy are reported on the effect of external gene intervention in a glioma network  <cit> .
methods
context-sensitive probabilistic boolean networks
for a network of n genes, a probabilistic boolean network  is defined by g, where v = {x <dig>  x <dig>  …, x
n
}, 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>   <dig>  …, n <cit> .
a node x
i
represents the state of gene i; x
i
=  <dig>  indicates that gene i is  expressed.
the set f
i
contains the rules that determine the next state of gene i. each fjii: <dig> n→ <dig>  for 1 ≤ j ≤ l is a mapping or predictor function determining the state of gene i. in a context-sensitive pbn  with k contexts, a network function is defined as fj=fj <dig> fj <dig> …,fjn, where fji: <dig> n→ <dig>  is the predictor function for gene i, i =  <dig>   <dig>  …, n, in context j   <cit> .
due to the stochasticity in genetic networks, the next state of gene i is determined by all the l functions in f
i
, i.e., by f1i,f2i,…,flii, with probabilities c1i,c2i,…,clii.
for an instantaneous pbn of n genes, there are a total number of ∏i=1nli possible boolean networks , each of which is a possible realization of the genetic network.
each bn can be considered as the network function for a context  <cit> .
in the jth context, assume the network function is given by fj=fj <dig> fj <dig> ⋯,fjnn, where each fjii: <dig> n→ <dig>  for 1 ≤ j ≤ l is a mapping or predictor function determining the state of gene i. the probability that the jth context is selected, is obtained as cj=∏i=1ncjii for j =  <dig>   <dig>  …, k, where k is the number of contexts in a context-sensitive network  <cit> .
the state of a gene is updated in the selected context; thus the next state depends on both the present state and the selected context.
in a context-sensitive network, a context may remain for certain time until a random event occurs.
a context switching usually occurs with probability q. for q =  <dig>  the cspbn becomes an instantaneous pbn.
for q <  <dig>  if a new context is to be selected, it is randomly chosen from the set of network functions: {f <dig>  f <dig>  …, f
k
}, with a set of context selection probabilities: {c <dig>  c <dig>  …, c
k
}  <cit> .
if noise is considered in a cspbn, it is often referred to as a perturbation, by which a gene flips its state with a probability p .
it has been shown that a pbn with perturbation is an ergodic markov chain in that all the states are connected in the pbn  <cit> .
the transition probability for any two states is determined by the values of p and q pairs.
following  <cit> , one of four mutually exclusive events occurs at time t in a cspbn with perturbation:
∅1: the predictor functions in the currently selected context are applied to update the gene expressions and this context remains for the next transition.
∅2: the predictor functions in the currently selected context are applied to update the states of the genes and then a new context is selected for the next transition.
∅3: a random perturbation occurs and the currently selected context remains for the next transition.
∅4: a random perturbation occurs and a new context is selected for the next transition.
the effect of a switching order, i.e., whether a network switches its context before or after its state transition, is considered in  <cit> , whereas in this paper, we focus on the transition rules, i.e., the network function is applied first and then the context switches.
a gene activity profile  is defined as a vector for describing the state of a network at time t,
x
= , x <dig>  …, x
n
), where x
i
∈ { <dig>  1} for i =  <dig>   <dig>  …, n. the state of a cspbn can be represented as a combination of a context and a gap, i.e.,
s
= .
by this definition, the number of states in a cspbn increases from 2
n
to 2
n
k for a network with n genes and k contexts.
a gap can also be represented by its decimal index, i.e., d=∑j=1n2n-jxjt+ <dig>  a state of a cspbn is then given by
s
=  for i ∈ { <dig>   <dig>  …, k} and d ∈ { <dig>   <dig>  …, 2
n
}, where d is the decimal index of a gap.
for convenience, a state  is referred to as  in the following analysis.
in a cspbn with a perturbation probability p and a switching probability q, the transition probability for any two states
a
and
b
is given by  <cit> :
pst+1=b|st=a=1-pnfr <dig> x <dig> x2+1-pn-hphsh⋅1-q+qcr1ga,b+qcr21-ga,b
where
fr <dig> x <dig> x2= <dig> ,ifx1directlytransitionstox2incontextr1otherwise
ga,b= <dig> ,ifr1=r2=rotherwise
sh= <dig> ,ifh=0otherwise
and h is the hamming distance between two gaps with decimal indices x <dig> and x <dig>  x <dig>  x <dig> ∈ { <dig>   <dig>  …, 2
n
}.
this distance indicates the number of genes with different expressions in the two gaps x <dig> and x <dig>
the state transition matrix  of a cspbn without perturbation is of the size 2
n
∙ k × 2
n
∙ k, given by  <cit> :
a=1-q+qc1p1⋯qc1pkqc2p1⋯qc2pk⋯⋯⋯qckp1⋯1-q+qckpk,
and the stm of a cspbn with perturbation is given by:
a=1-q+qc1p˜1⋯qc1p˜kqc2p˜1⋯qc2p˜k⋯⋯⋯qckp˜1⋯1-q+qckp˜k,
where
p
i
and p˜i denote the stms for the boolean network i, i ∈ { <dig>   <dig>  …, k}, without and with perturbation respectively.
context-sensitive stochastic boolean networks
a cssbn without perturbation
because of the similarities between logic circuits and biological networks, digital circuits have been used to simulate genetic networks  <cit> , as well as to determine the node vulnerability in cellular networks  <cit> .
stochastic logic has been demonstrated in several biological applications  <cit> .
stochastic computation implements probabilistic analysis using boolean logic by encoding a real number or probability into a random binary bit sequence  <cit> .
a probability is usually represented by a proportional number of 1’s in a stochastic sequence.
the complement of a probability can be computed by an inverter and the multiplication of probabilities can be implemented by an and gate for independent inputs.
a multiplexer computes a weighted sum of its input probabilities, with the weights given by the selection inputs.
figure  <dig> shows an inverter , an and gate, a buffer, an or gate, an xor gate and a multiplexer  <cit> .
due to stochastic fluctuations, the computational results by stochastic logic are not deterministic but probabilistic.
however, this stochastic fluctuation can be reduced through the use of non-bernoulli sequences of random permutations of fixed numbers of 1’s and 0’s as initial inputs, thus producing more accurate results than using bernoulli sequences  <cit> .
it has been shown that a probabilistic distribution of genes’ states can be accurately determined in a stochastic boolean network  with a reasonably long stochastic sequence length  <cit> .
in a cspbn, the selection of a context is determined by the context switching probability.
this probability indicates the likelihood to maintain the current context i or to select a new context from the k contexts .
based on the present state and the selected context, a gene’s state is updated.
if no perturbation occurs in an n-gene cspbn with a switching probability of q, the transition probability from state  to  is given by  <cit> :
pr,x|s,y=fs,y,x1-q+qcsfs,y,xqcrifr=sifr≠s
where
fs,y,x= <dig> ,ifydirectlytransitionstoxincontextsotherwise
s and r denote the sth and rth contexts, y and x represent two gene activity profiles , and c
s
and c
r
indicate the probability of selecting the sth and rth contexts respectively.
an sbn structure is proposed in  <cit>  to implement an instantaneous pbn.
for a cspbn, a context-sensitive sbn  is constructed to consider the switching of contexts, as shown in figure  <dig>  in this cssbn model, the probabilistic switching is implemented using a multiplexer in stochastic computation and the switching probability q is encoded as a random binary sequence q that serves as the control sequence of a 2-to- <dig> multiplexer .
if the jth bit in the sequence q is  <dig>  a new context will be selected for the next transition.
otherwise, the current context will remain.
the selection probability of the new context is determined by the original and current context selection probabilities.
as shown in the lower section of figure  <dig>  this process is implemented by 2-to- <dig> multiplexers with the original and current context selection sequences as inputs and q as the control sequence.
the selection probability of the new context is then obtained as encoded in the output stochastic sequences of the multiplexers.
if q =  <dig>  the cssbn functions as a fixed boolean network .
if q =  <dig>  the cssbn is simplified to an instantaneous sbn.
if a switch does occur, the context selection process is implemented by another multiplexer for choosing one from the k contexts, according to predefined selection probabilities.
as a network function is a combination of each gene’s predictor function, a combination of the m control sequences of s1 ~ s
m
is used to encode the predefined selection probabilities; this in turn determines the selection probability of each new context.
for a cspbn with n genes, the cssbn needs to be run for each of the 2
n
input states and sequences need to be generated for the m control signals of the multiplexer.
for a switching probability of q in the proposed cssbn, the currently selected context i remains at time t with a probability of q and switches to one of the k contexts with a probability of 1 - q. if the network transitions from a gap y to x in the sth context , then context s will remain with probability 1 - q + qc
s
at the next time step.
otherwise, the network moves into a new context r with probability qc
r
.
from this analysis, it can be seen that the cssbn in figure  <dig> computes the transition probability from state  to  as given by .
this indicates that the proposed cssbn model accurately implements the function of a cspbn.
a cssbn with perturbation
in a pbn with random perturbation, a gene may change its state with a probability p during a state transition.
following  <cit> , the effect of perturbation is considered to flip a gene’s state.
assume in an n-gene cspbn, the current gap at time t is given by
x
=  and
γ
is the perturbation vector, the gap at t +  <dig>
x′,
is given by:
x'=x⊕γfjxwithaprobabilityof1-1-pnwithaprobabilityof1-pn
where ⊕ is the addition modulo  <dig> and
f
j
is the network function for the jth context at time t. to account for the effect of perturbation, a cssbn with perturbation is constructed, as shown in figure  <dig>  in this cssbn, xor gates are used to implement the addition modulo  <dig> of the perturbation vector and the present state, while an n-input or gate is used to compute the probability that a perturbation occurs.
the output of the or gate is then used as the control sequence of a bus  multiplexer to decide the selection of sequences with or without perturbation.
if the output sequence of the or gate contains all 0’s, which means that there is no perturbation, then the next state is given by the predictor functions in the currently selected context in the original cssbn without perturbation; otherwise, the next state is determined by the perturbation effect.
a stochastic analysis of the function of the cssbn with perturbation shows that the next state of the network is given by:
x'=x⊕γ⋅1-1-pn+fjx⋅1-pn
which is equivalent to .
this indicates that a cspbn with perturbation can be accurately implemented by a cssbn with perturbation.
intervention in a cssbn
in contrast to perturbation, gene intervention refers to the process of deliberately changing the states of some genes to guide a network into a desired state  <cit> .
external stimuli are applied to a network to avoid undesirable states that might be associated with certain diseases.
for an effective intervention, control policies are developed for different intervention strategies; as several genes may affect the state of the target gene, a single gene that is most influential on the network state is usually identified as the control gene, due to its simplistic biological implications  <cit> .
the stm is then changed by an external intervention  <cit> .
in a cspbn, the context information is usually hidden and thus difficult to obtain in practice.
however, gene expressions can readily be observed, so a gap within all possible contexts is considered to be undesirable if the gap is an undesirable state of the network  <cit> .
hence, gene intervention refers to changing the state of a network as represented by a gap.
for an n-gene network, a vector
s
=  with s
i
∈ { <dig>  1} for any i ∈ { <dig>   <dig>  …, n}, is defined as the control gene vector, that is, if s
i
=  <dig>  then gene i is selected as a control gene  <cit> .
an intervention vector is defined as u=u <dig> u <dig> u <dig> …,u2n, u
j
∈ { <dig>  1} for any j ∈ { <dig>   <dig>  …, 2
n
}, where u
j
=  <dig>  indicates a flipping  of the state of a control gene at gap j. if gene i is selected to be a control gene, for example, then s
i
=  <dig>  the expression level of the control gene i is then determined by its current state and the status of
u
.
if u
j
=  <dig>  i.e., s
i
u
j
=  <dig>  the state of control gene i is flipped by the external intervention when gap j emerges; otherwise, the state of control gene i remains unchanged and thus the current state is preserved in gap j <cit> .
an intervention vector can be obtained by various methods; in this paper, a control policy using the steady state distribution   <cit>  is used to obtain
u
.
the state of a network under intervention, i.e., a gap, is determined by the control signal, the state prior to the intervention and the intervention vector.
let x^t and
x
t
be the gaps before and after intervention at time t; if the gth gene is the control gene , the state of the network under intervention is given by  <cit> :
xt=x^t⊕sg⋅1ux^t=1+x^t⋅1ux^t= <dig>
where
s
g
is a vector of n bits with  <dig> at every bit except the gth bit.
<dig> is an indicator function: 1ux^t=1= <dig> if ux^t= <dig> and 1ux^t=1= <dig> otherwise.
when ux^t= <dig>  the gap x^t is maintained; otherwise, the gap x^t is an undesirable state, for which the control bit is to be toggled, thereby improving the occurring probability of a desirable state.
the cssbn model under intervention is given in figure  <dig>  if gene i is selected as the control gene, i.e., s
i
=  <dig>  the output of the and gate is determined by the state of
u
.
if state j is an undesirable state, i.e., u
j
=  <dig>  the state of the control gene i is toggled.
this is implemented by an xor gate: if s
i
u
j
=  <dig>  the state of gene i, given by the output of an xor gate, is flipped at the state, or gap, j; otherwise, the state is not changed by the intervention.
from this analysis, the state of a network after intervention is given by:
xt=x^t⊕si,x^t,ifux^t=1otherwise
which is equivalent to .
this indicates that the cssbn model in figure  <dig> accurately implements the process of external gene intervention.
state transition matrix  and steady state analysis
as discussed previously, the state of a cspbn can be represented as
s
=  by considering the selected context and gap information.
using a cssbn, the stm can be obtained through the statistics encoded in the output sequences.
given an input state, each simulation of a cssbn produces the transition probabilities from this input to all output states, i.e., the row in the stm for this input.
for a cspbn with n genes and k contexts, the cssbn needs to be run for each of the 2
n
k input states and an o number of sequences need to be generated for the control signals of the multiplexers.
as in  <cit> , a factor, l, is used to account for the computational overhead required by using a longer stochastic sequence.
therefore, a complexity of o results for computing the stm of a context-sensitive network.
in a cssbn, the required minimum sequence length is typically on a polynomial order of the numbers of genes and contexts, as shown in the simulation results in table  <dig>  since the number of possible bns, k, generally increases exponentially with the number of predictor functions for each gene, the complexity of using a cssbn to compute the stm, i.e., o, is smaller than o of an accurate analysis for a cspbn  <cit> .
as indicated by the shorter average run time in table  <dig>  this difference becomes significant for a network with a large number of gene states.
in table  <dig>  the simulation results are also provided for cssbns with a single context, i.e., k =  <dig>  in this case, the cssbn degenerates into a deterministic bn, for which an analytical approach is more efficient, due to the use of random binary bit sequences in the stochastic approach.
however, the cssbn approach is more efficient in computing the stm for a gene network with a large number of states.
see additional file 1: the pseudocode for computing the stm of a cssbn.
n
k
num
l
l
n: the number of genes; k: the number of contexts.
perturbation and context switch probabilities: p =  <dig>  and q =  <dig>  .
l: the required minimum sequence length for a given accuracy measured by norm  <dig>  num: the number of states for the cssbn.
the results that at least one case of the cssbn study is more efficient than the cspbn, are highlighted in bold.
in general, a  ×  stm is required for an accurate cspbn analysis  <cit>  while a 2
n
× 2
n
stm is needed by an approximate method  <cit> .
as the number of genes increases in a network, a matrix-based analysis becomes infeasible due to a significant increase in the required computational resources.
the analysis of the ssd is even more challenging for a cspbn.
for a cssbn, however, the stm can be accurately and efficiently computed.
furthermore, the ssd can be evaluated through an iterative simulation in the temporal domain .
by this technique, an iterative structure of the cssbn is used to simulate the temporal evolution of a grn.
the required number of iterations is determined by the number of state transitions before reaching a steady state.
given an initial state distribution
x
, let
x
and
x
be the state distributions after m and m +  <dig> transitions respectively.
assume that the stm of the network is given by
a
, then the state transitions from
x
to
x
are described by:
xm+1=xm⋅a
if ||
x
-
x
||∞ is used to compute the maximum absolute value of the summed difference of each row in
x
and
x
, the condition for reaching a steady state is given by  <cit> :
xm+1-xm∞<ϵ
where ϵ indicates a threshold for determining whether the steady state has been reached or not.
if  is met, a network is considered to have reached a steady state, i.e.,
x
=
x
; thus
x
is considered as the stationary distribution.
alternatively, a reasonable number of network transitions  can be run to obtain the ssd in practice.
as shown in the results and discussion section, a time-frame extended simulation of cssbns is more efficient than an analytical approach while producing more accurate results than an approximate analysis.
see additional file 2: the pseudocode for computing the ssd of a cssbn.
in a time-frame expanded cssbn, random binary bit streams are generated for predictor function selection, gene perturbation and context switching probabilities.
as in  <cit> , non-bernoulli sequences of random permutations of fixed numbers of 1’s and 0’s are used for encoding initial probabilities.
these sequences then propagate through the cssbns and the statistics encoded in the output sequences are used to obtain the ssd.
compared to the use of bernoulli sequences, as shown in  <cit> , the use of the non-bernoulli sequences produces more accurate results for a given sequence length or requires a shorter sequence length for a desired evaluation accuracy.
as an efficient alternative to an stm-based analysis, the stochastic simulation of cssbns provides flexibility in achieving a tunable accuracy-efficiency tradeoff by using stochastic sequences of different lengths.
hence, the proposed cssbn approach is potentially useful in the analysis of large grns.
RESULTS
simulation with a p53-mdm <dig> network
as a tumour suppressor gene, p <dig> plays an important role in preventing the development and progression of tumour cells  <cit> .
in a p <dig> network, signaling pathways are triggered by external stimuli.
for example, dna damages activate pathways that involve the genes p <dig> and mdm <dig>  it has been shown that the expression level of the p <dig> protein is reversely related with that of the mdm <dig> gene, which leads to an oscillatory behavior of the p53-mdm <dig> network  <cit> .
various boolean models have been developed to simulate the dynamics of a p <dig> network  <cit> .
in this section, the two-gene probabilistic boolean network  model developed in  <cit>  for the p53-mdm <dig> network is used to illustrate the applicability of context-sensitive stochastic boolean networks .
this  pbn of the two genes p <dig> and mdm <dig> consists of v = {x <dig>  x2} and the predictor function sets f1=f <dig> f <dig> f <dig> f <dig> and f2=f <dig> f <dig> f <dig> f <dig>  the truth table for the state transitions of this pbn is given as table  <dig> <cit> .
x
1
,
x
2
are the present states of p <dig> and mdm2
x
1
x
2
an internal entry indicates whether a logical  <dig> or  <dig> would result at the next state of a gene from a selected boolean function.
cji indicates the transition probability by a function fji <cit> .
an internal entry indicates whether a logical  <dig> or  <dig> would result at the next state of a gene from a selected boolean function.
cji indicates the transition probability by a function fji <cit> .
in  <cit> , a stochastic boolean network  of the p53-mdm <dig> network is proposed.
based on the general cssbn model in figure  <dig>  a cssbn for the p53-mdm <dig> network is shown in figure  <dig>  in this model, a network function ) is selected at time t -  <dig> for context i and at time t, the network remains at context i with a probability of q or switches to one of the k contexts  with a probability of 1 - q. this context switching behavior is modeled by a 2-to- <dig> multiplexer  for each gene and the switching probability q is encoded in the control bit sequence q. at the same time, the selection probabilities of the new context are determined by the 2-to- <dig> multiplexers with the original and current context selection sequences as inputs and q as the control sequence, as shown in the upper section of figure  <dig>
the context-sensitive p53-mdm <dig> network has  <dig> genes and  <dig> contexts.
if both context and gene activity profile  are considered, there are 22 × 16 =  <dig> states for the p53-mdm <dig> cssbn.
given the transition probability for each predictor function , the probability for selecting each context can readily be computed for independent functions, as shown in table  <dig> with random perturbation to the genes, a cssbn with perturbation is constructed as shown in figure  <dig>  for this two-gene network, a two-input multiplexer is used for each gene to probabilistically select a perturbed state or the original cssbn state without perturbation.
the control bits to select a predictor function in figure  <dig> are also listed.
s0s1: the control bits for p53; s2s3: the control bits for mdm <dig>
the cssbn is then used to obtain the state transition matrix  for the p53-mdm <dig> network.
see additional file 3: the matlab program that describes the structure of the cssbn in figure  <dig> and computes its stm for the p53-mdm <dig> network with perturbation.
the norms || ⋅ || <dig>  || ⋅ || <dig>  and || ⋅ ||∞ are used to measure the differences of the stms obtained for the cssbn and cspbn.
|| ⋅ || <dig> and || ⋅ ||∞ indicate the maximum absolute values of the summed differences of the columns and rows respectively, while || ⋅ || <dig> measures the average difference of all entries.
assume
a
cssbn
and
a
cspbn
are the obtained stms for the cssbn and cspbn respectively.
let Δ
a
=
a
cssbn
-
a
cspbn
; the norms of the differences of the computed matrices  are then shown in table  <dig> for different values of the switching probability q and perturbation rate p. the average run time is also shown for using the cssbn.
l
q
= 1
p
q
=  <dig>
p
q
=  <dig>
p
q
=  <dig>
p
||Δ
a
||Δ
a
||Δ
a
the average run time of  <dig> simulations using the cssbn is also shown for different perturbation rates, p, and context switching probabilities, q. l: sequence length.
as revealed in table  <dig>  the difference in the stms computed using the cssbn and an analytical cspbn approach is significantly reduced by increasing the sequence length l. however, the inaccuracies, due to the inherent stochastic fluctuations in stochastic computation, are generally small and thus negligible.
hence, the proposed cssbn model can be used to accurately and efficiently compute the stm of a cspbn.
a steady state analysis is further performed on the proposed cssbn.
for the p53-mdm <dig> network, there are four states or gaps.
the probability of each gap is given by the sum of the probabilities for all contexts.
the simulation results for the four gaps with respect to different p and q values are shown in figure  <dig> .as can be seen in figure  <dig>  while the ssd is determined by both the perturbation and switching probabilities, the perturbation rate has a greater effect on the final distribution of the steady state compared to the switching probability.
the ssd changes drastically with the increase of the perturbation rate, whereas the effect of context switching is rather limited for a given perturbation rate.
as p, q ∈  <cit> , several p and q pairs are chosen for further analysis.
the difference in the ssds obtained by using the time-frame expanded cssbn technique, the approximate  <cit>  and accurate analysis  <cit>  are shown in table  <dig>
p
p
p
q
q
q
||Δ
ssd
ap - ac
||Δ
ssd
cssbn - ac
||Δ
ssd
ap - ac
||Δ
ssd
cssbn - ac
||Δ
ssd
ap - ac
||Δ
ssd
cssbn - ac
Δ
ssd
ap - ac
: the difference between the results obtained by using the approximate and accurate analysis.
Δ
ssd
cssbn - ac
: the difference between the results obtained by using the cssbn approach and the accurate analysis.
p: perturbation rate, q: context switching probability, l: sequence length.
as revealed in table  <dig>  the cssbn approach can compute the ssd more accurately than the approximate analysis.
in fact, the difference in the results between the cssbn and the accurate analysis is negligible when reasonably long stochastic sequences are used.
with the stm obtained for a cssbn, an ssd is evaluated and the results are very accurate compared to those obtained using the accurate approach, as shown in figure  <dig>  a further analysis shows that the relative error is limited to less than approximately  <dig> % for the cssbn approach.the ssds obtained using different methods, as well as the corresponding values of norms, are shown in figure  <dig>  it can be seen that the ssd can be accurately evaluated by the cssbn model compared to the accurate analytical approach.
it has been shown that a pbn with perturbation evolves as an ergodic markov chain  <cit> .
for a larger perturbation probability p, there is an increased randomness in the network activities, thus the steady states of the network are more evenly distributed  <cit> .
this can be seen in the simulation results in figure  <dig>  figures  <dig> and  <dig> further show that context switching has little impact on the ssd of the network with perturbation.
this is due to the fact that the context switching activity only affects the selection probability of a context, but not the predictor functions.
also shown in figure  <dig> is that an oscillation of the states of p <dig> and mdm <dig> exists, indicated by the higher probabilities of the network states   <dig> and  <dig> compared with the other states.
however, this two-state oscillation is less evident when the perturbation rate is significantly higher.
on the other hand, a modest perturbation to the network only has a minor effect on the oscillatory behavior, confirming the stability and robustness of the p53-mdm <dig> network in a dynamic environment  <cit> .
experiments on a glioma network
a network in  <cit>  obtained from a human glioma gene expression data set  <cit>  is further used to illustrate the efficiency of the cssbn model and the time-frame expansion technique.
based on this data set, a pbn was inferred by a method using the coefficients of determination  and an ssd analysis was performed in  <cit> .
an approximate method was developed in  <cit>  as an extension to estimate the ssd of a cspbn with perturbation.
the gene top2a was not considered in either study as it is an input gene with an in-degree of zero.
in our study, the network setting is considered the same as in  <cit>  with the gene top2a removed; this leads to a total of  <dig> gaps.
figure  <dig> shows a detailed structure of the considered glioma network with double -headed arrows indicating the bi -directional relationships of gene pairs.
the selection probabilities for the predictor functions are shown in table  <dig> <cit>  and the boolean functions for each gene are obtained through an analysis of the data in  <cit> .
f
i
,
i
∈ { <dig>   <dig>  …, 14}, for each gene in the glioma network in figure <dig> <cit>
f
1
f
2
f
3
f
4
f
5
f
6
f
7
f
8
f
9
f
10
f
11
f
12
f
13
f
14
in table  <dig>  each column for f
i
, i ∈ { <dig>   <dig>  …, 14}, contains the selection probabilities for the boolean functions of gene i, with the value in the jth row as the probability for fji, thus the probabilities sum to  <dig> in each column.
table  <dig> also shows that six genes have only a single predictor function; only the state of gene  <dig> is determined by three functions, while the states of the other genes are determined by two predictor functions.
steady state analysis
for the 14-gene glioma network, there are a total of  <dig> contexts, as can be determined from table  <dig>  it requires an stm with 214 × 384 =  <dig> columns and rows for an accurate analysis.
this makes it infeasible to estimate the steady states of a cspbn using a matrix-based analysis.
the approximate analysis in  <cit>  would require the computation of an stm of the size 214 ×  <dig>  thus, it is difficult in general to use an analytical approach to evaluate a large network, due to the demanding computational resources required.
however, a cssbn model can be constructed for the glioma network; this cssbn is based on the constituent sbn, as shown in figure  <dig>  with the cssbn, the ssd can be estimated using the aforementioned time-frame expansion technique with a greatly reduced complexity.
the obtained ssd is then compared with that obtained from the approximate analysis in  <cit> .
in this paper, a network is considered to have reached a steady state if the discrepancy between two adjacent iterations is smaller than a required threshold ϵ ) or the number of simulation iterations has reached a maximum constant value.
see additional file 4: the matlab program that evaluates the steady state distribution using the time-frame expanded cssbn technique for the glioma network with context switching and perturbation.
the state or gap of the glioma network can be represented by a binary vector as , x
i
∈ { <dig>  1} for i ∈ { <dig>   <dig>  …, 14}, or its decimal index.
for example, the state  can be represented as state  <dig>  the ssds of the context-sensitive glioma network with all  <dig> states, obtained using the cssbn approach and the approximate analysis  <cit> , are shown in figure  <dig>
the norms of the differences of ssds, obtained using the cssbn with different sequence lengths and the approximate method in  <cit> , are shown in table  <dig>
p
=  <dig>
q
p
=  <dig>
q
l
||Δ
ssd
ap - cssbn
||Δ
ssd
ap - cssbn
||Δ
ssd
ap - cssbn
a maximum of  <dig> iterations are used for obtaining the ssd.
p: perturbation rate; q: context switching probability; l: sequence length.
Δ
ssd
ap - cssbn
: the difference between the ssds obtained by the approximate analysis  <cit>  and the context-sensitive stochastic boolean network  approach.
as shown in table  <dig>  the cssbn time-frame expansion technique efficiently evaluates the ssd of the glioma network and produces results comparable to those obtained by the approximate analysis  <cit> .
evaluation accuracy is further improved by using longer stochastic sequences with yet a shorter runtime compared to the approximate method.
since it is difficult, if not impossible, to compute the stm or ssd of a large grn by using an accurate or approximate analysis, the cssbn time-frame expansion approach provides an alternative means to accurately and efficiently estimate the ssd of a large network with a tunable accuracy by using stochastic sequences of different lengths.
intervention analysis
in an n-gene network, if gene x
i
is the target gene, the states of all genes can be divided into a set of desirable states, d, and a set of undesirable states, u, by the expression level of the target gene  <cit> .
as can be seen in figure  <dig>  gene x <dig>  the ; bfl <dig> protein; grs protein is one of the most influential genes that interacts closely with others, thus it is chosen as the target gene.
the desirable and undesirable states are then separated by the expression level of x <dig>  assume that the inactivation of x <dig> is preferred; the cumulative probabilities of the desirable and undesirable states are given by ∑x14=0πx and ∑x14=1πy, respectively, where π
x
and π
y
are elements in the desirable and undesirable ssds respectively.
as discussed previously, a grn can be intervened by applying external stimuli to minimize the likelihood of being in an undesirable state.
various methods have been proposed for deriving control vectors for an effective intervention  <cit> ; in this study the ssd algorithm proposed in  <cit>  is used to determine an intervention vector.
for simplicity,  <dig> most significant contexts that account for a total selection probability of  <dig> % are chosen from the total  <dig> contexts of the glioma network for an intervention analysis.
for these  <dig> contexts, the expression levels of gene x <dig> and x <dig> have no effect on the states of other genes, so these two genes are removed; this yields a simplified 12-gene glioma network.
the cumulative distributions of the desirable steady states of the context-sensitive 12-gene glioma network are shown in table  <dig> for using a different gene as the control gene.
p
=  <dig>  and context switching probability
q
=  <dig>
x
1
x
2
x
3
x
4
x
6
x
8
x
9
x
10
x
11
x
12
x
13
the cumulative probability of the desirable states with no intervention is  <dig> %.
sequence length l = 800k bits.
when a different gene is selected as the control gene, the effect of intervention varies with respect to improving the probabilities of desirable states in the ssd.
for the glioma network, as revealed in table  <dig>  gene x <dig>  tie- <dig>  is the most effective for maximizing the percentage of the desirable states.
the cumulative distribution of the desirable states is increased from  <dig> % to  <dig> % by an intervention using tie- <dig> as the control gene.
also revealed in the table is that an intervention via gene x <dig>  ; , has nearly no impact on the distribution of the desirable states, as the percentage of the desirable states has not been changed much by the intervention .
a modest improvement in the desirable state distribution is obtained by an intervention with another gene as the control gene.
for any target gene, this process can be applied to find the most significant gene that maximizes the probabilities of desirable states.
the effects of intervention can also be seen in figure  <dig>  where three different scenarios are considered:  no intervention,  an indirect intervention via gene x <dig>  i.e., tie- <dig>  and  a direct intervention via the target gene x <dig>  the ; bfl <dig> protein; grs protein.
the cumulative probabilities of the desirable states are  <dig> % for no intervention,  <dig> % and  <dig> % for the two intervention strategies.
as revealed in these results, a direct intervention of the target gene is perhaps optimal for avoiding the undesirable states and almost certainly taking the network into a desirable state.
however, when an intervention on the target gene is not possible, an intervention through tie- <dig> is the most effective means to maximize the probability of the down-regulation of the target gene x <dig>
CONCLUSIONS
context-sensitive stochastic boolean networks  are proposed as an efficient approach to modeling the effects of gene perturbation and intervention in gene regulatory networks .
in a cssbn, the state transition matrix can be accurately and efficiently computed with a complexity of o, where n is the number of genes in a context-sensitive probabilistic boolean network , k is the number of contexts and l is a factor determined by the stochastic sequence length.
this result is an improvement compared to the previous result of o for an accurate analysis.
the use of non-bernoulli stochastic sequences further increases computational efficiency and allows for a tunable tradeoff between accuracy and efficiency.
a steady state analysis using a time-frame expansion technique has shown a significant speedup and produced more accurate results than an approximate analysis in the computation of the steady state distribution.
cssbns are constructed for the analysis of gene perturbation in a p53-mdm <dig> network and gene intervention in a glioma network.
biologically meaningful insights are gained into the oscillatory dynamics of the p53-mdm <dig> network in a context-switching environment with random gene perturbation.
it has been shown that the steady state distribution  changes drastically with the increase of the perturbation rate, whereas the effect of context switching is rather limited for a given perturbation rate.
therefore, random gene perturbation may have a greater effect on the final distribution of the steady state compared to context switching activities.
by predicting the ssd, the cssbn approach can further evaluate the effectiveness of external gene intervention.
a case study of the glioma network shows that the cssbns are useful in modeling the effects of gene perturbation and intervention in a complex context-sensitive grn.
this will eventually help drug discovery for an effective drug intervention therapy.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
pz carried out the grn studies, performed the statistical analysis and drafted the manuscript.
jl developed the initial models and participated in discussion.
jh participated in the grn studies, revised the manuscript and supervised the research.
all authors read and approved the final manuscript.
supplementary material
additional file 1
the pseudocode for computing the state transition matrix  of a context-sensitive stochastic boolean network .
click here for file
additional file 2
the pseudocode for computing the steady-state distribution  of a context-sensitive stochastic boolean network .
click here for file
additional file 3
the matlab program that describes the structure of the context-sensitive stochastic boolean network  in figure 6and computes its state transition matrix for the p53-mdm <dig> network with perturbation.
click here for file
additional file 4
the matlab program that evaluates the steady state distribution using the time-frame expanded context-sensitive stochastic boolean network  technique for the glioma network with context switching and perturbation.
click here for file
acknowledgements
this work was supported in part by the natural sciences and engineering research council  of canada in a discovery grant and a startup fund of the university of alberta.

 diverse biological functions are regulated through the interactions among genes, proteins and other molecules in a cell. gene expressions are however affected by the intrinsic and extrinsic noise in a gene network. a major source of the noise is the stochastic fluctuations in gene regulatory interactions. the genetic interactions are also context dependent, that is, certain regulatory functions are active in some cellular states, but inactive in others. this indicates the necessity to consider noise in a context-sensitive manner in the study of gene regulatory networks. various methods have been proposed to model grns; these include logical models, continuous models using differential equations and stochastic models at the single-molecule level. as a classic logical model, boolean networks have been widely used to qualitatively model the interactions among genes. probabilistic boolean networks have been proposed to consider noise in a bn model. in a pbn, the next state of a gene is determined by its current state and a boolean function. if the boolean function is randomly selected, a pbn is referred to as an instantaneous pbn. as a general model, a context-sensitive pbn considers the feature of context dependence in a biological network. in a cspbn, a context is a combination of boolean functions and each function determines the next state of a gene. a context remains unchanged until a switching occurs. this switching of contexts, possibly caused by external stimuli, is considered to occur randomly in a network. the study of pbns has focused on the analysis of steady state distributions under gene perturbation and intervention. a markov chain monte carlo method is used in to analyze the long run behavior of a pbn; however, this method generally requires a large number of simulations to reach a steady state, due to the slow convergence typically encountered in a monte carlo method. an analysis is performed in for finding the ssd of a pbn through the computation of the state transition matrix. however, the application of an analytical approach is generally limited to small networks due to the exponential increase in the size of an stm with gene numbers. the analysis of cspbns presents even a greater challenge due to its significantly increased computational complexity. analytical expressions have been derived for analyzing cspbns; however, this method is only applicable to the steady state analysis of a network with small perturbation and switching probabilities. the method in ignores some bns with very small probabilities for reducing the size of the stm and thus provides a more efficient but approximate solution for computing the ssd of a cspbn. in, gene intervention is investigated for avoiding undesirable states associated with certain diseases. due to external stimuli, the stm is changed by external control variables, so desirable states can be obtained with larger probabilities in the ssd. in a context-sensitive network with n genes and k contexts, however, a  ×  matrix is required for an accurate analysis of the ssd; this results in a computational complexity of o ) for an accurate computation of the stm. hence, the application of current gene network analysis is limited to those of less than a dozen of genes. recently, stochastic boolean networks have been proposed for an efficient computation of the stm and ssd of an instantaneous pbn. the sbn approach can recover biologically-proven regulatory behaviors, such as the oscillatory dynamics of a simplified p53-mdm network and the dynamic attractors in a t cell immune response network. to further exploit the simplicity of logical models, stochastic multiple-valued networks have been utilized to investigate the dynamics of gene networks with multiple-valued gene states. asynchronous sbns have also been developed to investigate the asynchronous state-update behavior of genes. in this paper, the notion of sbns is extended for the general model of pbns, i.e., the cspbns. this presents a novel structure of context-sensitive sbns for modeling the stochasticity in a context-dependent grn. in particular, gene perturbation and intervention are considered in cssbns. through an efficient simulation of a cssbn, the computational complexity in the evaluation of a cspbn is reduced from o to o for computing the stm, where l is a factor related to the required sequence length in cssbn for achieving a desired accuracy. the use of non-bernoulli sequences of random permutations of fixed number of 1’s and 0’s further increases computational efficiency and allows for a tunable tradeoff between accuracy and efficiency. a time-frame expanded cssbn can further simulate the stationary behavior of a cspbn and produce results that are more accurate than an approximate analysis, thus making the cssbn useful in modeling complex context-sensitive grns. the cssbn models are applied to the study of a simplified p53-mdm network and the results using an ssd control policy are reported on the effect of external gene intervention in a glioma network. methods context-sensitive probabilistic boolean networks for a network of n genes, a probabilistic boolean network is defined by g, where v = {x x …, x n }, 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. a node x i represents the state of gene i; x i  =  indicates that gene i is expressed. the set f i contains the rules that determine the next state of gene i. each fjii: n→ for 1 ≤ j ≤ l is a mapping or predictor function determining the state of gene i. in a context-sensitive pbn with k contexts, a network function is defined as fj=fj fj …,fjn, where fji: n→ is the predictor function for gene i, i =  …, n, in context j. due to the stochasticity in genetic networks, the next state of gene i is determined by all the l functions in f i, i.e., by f1i,f2i,…,flii, with probabilities c1i,c2i,…,clii. for an instantaneous pbn of n genes, there are a total number of ∏i=1nli possible boolean networks, each of which is a possible realization of the genetic network. each bn can be considered as the network function for a context. in the jth context, assume the network function is given by fj=fj fj ⋯,fjnn, where each fjii: n→ for 1 ≤ j ≤ l is a mapping or predictor function determining the state of gene i. the probability that the jth context is selected, is obtained as cj=∏i=1ncjii for j =  …, k, where k is the number of contexts in a context-sensitive network. the state of a gene is updated in the selected context; thus the next state depends on both the present state and the selected context. in a context-sensitive network, a context may remain for certain time until a random event occurs. a context switching usually occurs with probability q. for q =  the cspbn becomes an instantaneous pbn. for q <  if a new context is to be selected, it is randomly chosen from the set of network functions: {f f …, f k }, with a set of context selection probabilities: {c c …, c k }. if noise is considered in a cspbn, it is often referred to as a perturbation, by which a gene flips its state with a probability p. it has been shown that a pbn with perturbation is an ergodic markov chain in that all the states are connected in the pbn. the transition probability for any two states is determined by the values of p and q pairs. following, one of four mutually exclusive events occurs at time t in a cspbn with perturbation: ∅1: the predictor functions in the currently selected context are applied to update the gene expressions and this context remains for the next transition. ∅2: the predictor functions in the currently selected context are applied to update the states of the genes and then a new context is selected for the next transition. ∅3: a random perturbation occurs and the currently selected context remains for the next transition. ∅4: a random perturbation occurs and a new context is selected for the next transition. the effect of a switching order, i.e., whether a network switches its context before or after its state transition, is considered in, whereas in this paper, we focus on the transition rules, i.e., the network function is applied first and then the context switches. a gene activity profile is defined as a vector for describing the state of a network at time t, x  =, x …, x n ), where x i ∈ { 1} for i =  …, n. the state of a cspbn can be represented as a combination of a context and a gap, i.e., s  =. by this definition, the number of states in a cspbn increases from 2 n to 2 n k for a network with n genes and k contexts. a gap can also be represented by its decimal index, i.e., d=∑j=1n2n-jxjt+ a state of a cspbn is then given by s  =  for i ∈ { …, k} and d ∈ { …, 2 n }, where d is the decimal index of a gap. for convenience, a state is referred to as in the following analysis. in a cspbn with a perturbation probability p and a switching probability q, the transition probability for any two states a and b is given by: pst+1=b|st=a=1-pnfr x x2+1-pn-hphsh⋅1-q+qcr1ga,b+qcr21-ga,b where fr x x2=,ifx1directlytransitionstox2incontextr1otherwise ga,b=,ifr1=r2=rotherwise sh=,ifh=0otherwise and h is the hamming distance between two gaps with decimal indices x and x x x ∈ { …, 2 n }. this distance indicates the number of genes with different expressions in the two gaps x and x the state transition matrix of a cspbn without perturbation is of the size 2 n ∙ k × 2 n ∙ k, given by: a=1-q+qc1p1⋯qc1pkqc2p1⋯qc2pk⋯⋯⋯qckp1⋯1-q+qckpk, and the stm of a cspbn with perturbation is given by: a=1-q+qc1p˜1⋯qc1p˜kqc2p˜1⋯qc2p˜k⋯⋯⋯qckp˜1⋯1-q+qckp˜k, where p i and p˜i denote the stms for the boolean network i, i ∈ { …, k}, without and with perturbation respectively. context-sensitive stochastic boolean networks a cssbn without perturbation because of the similarities between logic circuits and biological networks, digital circuits have been used to simulate genetic networks, as well as to determine the node vulnerability in cellular networks. stochastic logic has been demonstrated in several biological applications. stochastic computation implements probabilistic analysis using boolean logic by encoding a real number or probability into a random binary bit sequence. a probability is usually represented by a proportional number of 1’s in a stochastic sequence. the complement of a probability can be computed by an inverter and the multiplication of probabilities can be implemented by an and gate for independent inputs. a multiplexer computes a weighted sum of its input probabilities, with the weights given by the selection inputs. figure  shows an inverter, an and gate, a buffer, an or gate, an xor gate and a multiplexer. due to stochastic fluctuations, the computational results by stochastic logic are not deterministic but probabilistic. however, this stochastic fluctuation can be reduced through the use of non-bernoulli sequences of random permutations of fixed numbers of 1’s and 0’s as initial inputs, thus producing more accurate results than using bernoulli sequences. it has been shown that a probabilistic distribution of genes’ states can be accurately determined in a stochastic boolean network with a reasonably long stochastic sequence length. in a cspbn, the selection of a context is determined by the context switching probability. this probability indicates the likelihood to maintain the current context i or to select a new context from the k contexts. based on the present state and the selected context, a gene’s state is updated. if no perturbation occurs in an n-gene cspbn with a switching probability of q, the transition probability from state to is given by: pr,x|s,y=fs,y,x1-q+qcsfs,y,xqcrifr=sifr≠s where fs,y,x=,ifydirectlytransitionstoxincontextsotherwise s and r denote the sth and rth contexts, y and x represent two gene activity profiles, and c s and c r indicate the probability of selecting the sth and rth contexts respectively. an sbn structure is proposed in to implement an instantaneous pbn. for a cspbn, a context-sensitive sbn is constructed to consider the switching of contexts, as shown in figure  in this cssbn model, the probabilistic switching is implemented using a multiplexer in stochastic computation and the switching probability q is encoded as a random binary sequence q that serves as the control sequence of a 2-to- multiplexer. if the jth bit in the sequence q is a new context will be selected for the next transition. otherwise, the current context will remain. the selection probability of the new context is determined by the original and current context selection probabilities. as shown in the lower section of figure  this process is implemented by 2-to- multiplexers with the original and current context selection sequences as inputs and q as the control sequence. the selection probability of the new context is then obtained as encoded in the output stochastic sequences of the multiplexers. if q =  the cssbn functions as a fixed boolean network. if q =  the cssbn is simplified to an instantaneous sbn. if a switch does occur, the context selection process is implemented by another multiplexer for choosing one from the k contexts, according to predefined selection probabilities. as a network function is a combination of each gene’s predictor function, a combination of the m control sequences of s1 ~ s m is used to encode the predefined selection probabilities; this in turn determines the selection probability of each new context. for a cspbn with n genes, the cssbn needs to be run for each of the 2 n input states and sequences need to be generated for the m control signals of the multiplexer. for a switching probability of q in the proposed cssbn, the currently selected context i remains at time t with a probability of q and switches to one of the k contexts with a probability of 1 - q. if the network transitions from a gap y to x in the sth context, then context s will remain with probability 1 - q + qc s at the next time step. otherwise, the network moves into a new context r with probability qc r. from this analysis, it can be seen that the cssbn in figure  computes the transition probability from state to as given by. this indicates that the proposed cssbn model accurately implements the function of a cspbn. a cssbn with perturbation in a pbn with random perturbation, a gene may change its state with a probability p during a state transition. following, the effect of perturbation is considered to flip a gene’s state. assume in an n-gene cspbn, the current gap at time t is given by x  =  and γ is the perturbation vector, the gap at t +  x′, is given by: x'=x⊕γfjxwithaprobabilityof1-1-pnwithaprobabilityof1-pn where ⊕ is the addition modulo and f j is the network function for the jth context at time t. to account for the effect of perturbation, a cssbn with perturbation is constructed, as shown in figure  in this cssbn, xor gates are used to implement the addition modulo of the perturbation vector and the present state, while an n-input or gate is used to compute the probability that a perturbation occurs. the output of the or gate is then used as the control sequence of a bus multiplexer to decide the selection of sequences with or without perturbation. if the output sequence of the or gate contains all 0’s, which means that there is no perturbation, then the next state is given by the predictor functions in the currently selected context in the original cssbn without perturbation; otherwise, the next state is determined by the perturbation effect. a stochastic analysis of the function of the cssbn with perturbation shows that the next state of the network is given by: x'=x⊕γ⋅1-1-pn+fjx⋅1-pn which is equivalent to. this indicates that a cspbn with perturbation can be accurately implemented by a cssbn with perturbation. intervention in a cssbn in contrast to perturbation, gene intervention refers to the process of deliberately changing the states of some genes to guide a network into a desired state. external stimuli are applied to a network to avoid undesirable states that might be associated with certain diseases. for an effective intervention, control policies are developed for different intervention strategies; as several genes may affect the state of the target gene, a single gene that is most influential on the network state is usually identified as the control gene, due to its simplistic biological implications. the stm is then changed by an external intervention. in a cspbn, the context information is usually hidden and thus difficult to obtain in practice. however, gene expressions can readily be observed, so a gap within all possible contexts is considered to be undesirable if the gap is an undesirable state of the network. hence, gene intervention refers to changing the state of a network as represented by a gap. for an n-gene network, a vector s  =  with s i ∈ { 1} for any i ∈ { …, n}, is defined as the control gene vector, that is, if s i  =  then gene i is selected as a control gene. an intervention vector is defined as u=u u u …,u2n, u j ∈ { 1} for any j ∈ { …, 2 n }, where u j  =  indicates a flipping of the state of a control gene at gap j. if gene i is selected to be a control gene, for example, then s i  =  the expression level of the control gene i is then determined by its current state and the status of u. if u j  =  i.e., s i u j  =  the state of control gene i is flipped by the external intervention when gap j emerges; otherwise, the state of control gene i remains unchanged and thus the current state is preserved in gap j. an intervention vector can be obtained by various methods; in this paper, a control policy using the steady state distribution is used to obtain u. the state of a network under intervention, i.e., a gap, is determined by the control signal, the state prior to the intervention and the intervention vector. let x^t and x t be the gaps before and after intervention at time t; if the gth gene is the control gene, the state of the network under intervention is given by: xt=x^t⊕sg⋅1ux^t=1+x^t⋅1ux^t= where s g is a vector of n bits with at every bit except the gth bit. is an indicator function: 1ux^t=1= if ux^t= and 1ux^t=1= otherwise. when ux^t= the gap x^t is maintained; otherwise, the gap x^t is an undesirable state, for which the control bit is to be toggled, thereby improving the occurring probability of a desirable state. the cssbn model under intervention is given in figure  if gene i is selected as the control gene, i.e., s i  =  the output of the and gate is determined by the state of u. if state j is an undesirable state, i.e., u j  =  the state of the control gene i is toggled. this is implemented by an xor gate: if s i u j  =  the state of gene i, given by the output of an xor gate, is flipped at the state, or gap, j; otherwise, the state is not changed by the intervention. from this analysis, the state of a network after intervention is given by: xt=x^t⊕si,x^t,ifux^t=1otherwise which is equivalent to. this indicates that the cssbn model in figure  accurately implements the process of external gene intervention. state transition matrix and steady state analysis as discussed previously, the state of a cspbn can be represented as s  =  by considering the selected context and gap information. using a cssbn, the stm can be obtained through the statistics encoded in the output sequences. given an input state, each simulation of a cssbn produces the transition probabilities from this input to all output states, i.e., the row in the stm for this input. for a cspbn with n genes and k contexts, the cssbn needs to be run for each of the 2 n k input states and an o number of sequences need to be generated for the control signals of the multiplexers. as in, a factor, l, is used to account for the computational overhead required by using a longer stochastic sequence. therefore, a complexity of o results for computing the stm of a context-sensitive network. in a cssbn, the required minimum sequence length is typically on a polynomial order of the numbers of genes and contexts, as shown in the simulation results in table  since the number of possible bns, k, generally increases exponentially with the number of predictor functions for each gene, the complexity of using a cssbn to compute the stm, i.e., o, is smaller than o of an accurate analysis for a cspbn. as indicated by the shorter average run time in table  this difference becomes significant for a network with a large number of gene states. in table  the simulation results are also provided for cssbns with a single context, i.e., k =  in this case, the cssbn degenerates into a deterministic bn, for which an analytical approach is more efficient, due to the use of random binary bit sequences in the stochastic approach. however, the cssbn approach is more efficient in computing the stm for a gene network with a large number of states. see additional file 1: the pseudocode for computing the stm of a cssbn. n k num l l n: the number of genes; k: the number of contexts. perturbation and context switch probabilities: p =  and q = . l: the required minimum sequence length for a given accuracy measured by norm num: the number of states for the cssbn. the results that at least one case of the cssbn study is more efficient than the cspbn, are highlighted in bold. in general, a  ×  stm is required for an accurate cspbn analysis while a 2 n  × 2 n stm is needed by an approximate method. as the number of genes increases in a network, a matrix-based analysis becomes infeasible due to a significant increase in the required computational resources. the analysis of the ssd is even more challenging for a cspbn. for a cssbn, however, the stm can be accurately and efficiently computed. furthermore, the ssd can be evaluated through an iterative simulation in the temporal domain. by this technique, an iterative structure of the cssbn is used to simulate the temporal evolution of a grn. the required number of iterations is determined by the number of state transitions before reaching a steady state. given an initial state distribution x, let x and x be the state distributions after m and m +  transitions respectively. assume that the stm of the network is given by a, then the state transitions from x to x are described by: xm+1=xm⋅a if || x  -  x ||∞ is used to compute the maximum absolute value of the summed difference of each row in x and x, the condition for reaching a steady state is given by: xm+1-xm∞<ϵ where ϵ indicates a threshold for determining whether the steady state has been reached or not. if is met, a network is considered to have reached a steady state, i.e., x  =  x; thus x is considered as the stationary distribution. alternatively, a reasonable number of network transitions can be run to obtain the ssd in practice. as shown in the results and discussion section, a time-frame extended simulation of cssbns is more efficient than an analytical approach while producing more accurate results than an approximate analysis. see additional file 2: the pseudocode for computing the ssd of a cssbn. in a time-frame expanded cssbn, random binary bit streams are generated for predictor function selection, gene perturbation and context switching probabilities. as in, non-bernoulli sequences of random permutations of fixed numbers of 1’s and 0’s are used for encoding initial probabilities. these sequences then propagate through the cssbns and the statistics encoded in the output sequences are used to obtain the ssd. compared to the use of bernoulli sequences, as shown in, the use of the non-bernoulli sequences produces more accurate results for a given sequence length or requires a shorter sequence length for a desired evaluation accuracy. as an efficient alternative to an stm-based analysis, the stochastic simulation of cssbns provides flexibility in achieving a tunable accuracy-efficiency tradeoff by using stochastic sequences of different lengths. hence, the proposed cssbn approach is potentially useful in the analysis of large grns.  simulation with a p53-mdm network as a tumour suppressor gene, p plays an important role in preventing the development and progression of tumour cells. in a p network, signaling pathways are triggered by external stimuli. for example, dna damages activate pathways that involve the genes p and mdm it has been shown that the expression level of the p protein is reversely related with that of the mdm gene, which leads to an oscillatory behavior of the p53-mdm network. various boolean models have been developed to simulate the dynamics of a p network. in this section, the two-gene probabilistic boolean network model developed in for the p53-mdm network is used to illustrate the applicability of context-sensitive stochastic boolean networks. this pbn of the two genes p and mdm consists of v = {x x2} and the predictor function sets f1=f f f f and f2=f f f f the truth table for the state transitions of this pbn is given as table . x 1, x 2 are the present states of p and mdm2 x 1 x 2 an internal entry indicates whether a logical or would result at the next state of a gene from a selected boolean function. cji indicates the transition probability by a function fji. an internal entry indicates whether a logical or would result at the next state of a gene from a selected boolean function. cji indicates the transition probability by a function fji. in, a stochastic boolean network of the p53-mdm network is proposed. based on the general cssbn model in figure  a cssbn for the p53-mdm network is shown in figure  in this model, a network function ) is selected at time t -  for context i and at time t, the network remains at context i with a probability of q or switches to one of the k contexts with a probability of 1 - q. this context switching behavior is modeled by a 2-to- multiplexer for each gene and the switching probability q is encoded in the control bit sequence q. at the same time, the selection probabilities of the new context are determined by the 2-to- multiplexers with the original and current context selection sequences as inputs and q as the control sequence, as shown in the upper section of figure  the context-sensitive p53-mdm network has genes and contexts. if both context and gene activity profile are considered, there are 22 × 16 =  states for the p53-mdm cssbn. given the transition probability for each predictor function, the probability for selecting each context can readily be computed for independent functions, as shown in table  with random perturbation to the genes, a cssbn with perturbation is constructed as shown in figure  for this two-gene network, a two-input multiplexer is used for each gene to probabilistically select a perturbed state or the original cssbn state without perturbation. the control bits to select a predictor function in figure  are also listed. s0s1: the control bits for p53; s2s3: the control bits for mdm the cssbn is then used to obtain the state transition matrix for the p53-mdm network. see additional file 3: the matlab program that describes the structure of the cssbn in figure  and computes its stm for the p53-mdm network with perturbation. the norms || ⋅ || || ⋅ || and || ⋅ ||∞ are used to measure the differences of the stms obtained for the cssbn and cspbn. || ⋅ || and || ⋅ ||∞ indicate the maximum absolute values of the summed differences of the columns and rows respectively, while || ⋅ || measures the average difference of all entries. assume a cssbn and a cspbn are the obtained stms for the cssbn and cspbn respectively. let Δ a  =  a cssbn  -  a cspbn; the norms of the differences of the computed matrices are then shown in table  for different values of the switching probability q and perturbation rate p. the average run time is also shown for using the cssbn. l q  = 1 p q  =  p q  =  p q  =  p ||Δ a ||Δ a ||Δ a the average run time of simulations using the cssbn is also shown for different perturbation rates, p, and context switching probabilities, q. l: sequence length. as revealed in table  the difference in the stms computed using the cssbn and an analytical cspbn approach is significantly reduced by increasing the sequence length l. however, the inaccuracies, due to the inherent stochastic fluctuations in stochastic computation, are generally small and thus negligible. hence, the proposed cssbn model can be used to accurately and efficiently compute the stm of a cspbn. a steady state analysis is further performed on the proposed cssbn. for the p53-mdm network, there are four states or gaps. the probability of each gap is given by the sum of the probabilities for all contexts. the simulation results for the four gaps with respect to different p and q values are shown in figure .as can be seen in figure  while the ssd is determined by both the perturbation and switching probabilities, the perturbation rate has a greater effect on the final distribution of the steady state compared to the switching probability. the ssd changes drastically with the increase of the perturbation rate, whereas the effect of context switching is rather limited for a given perturbation rate. as p, q ∈, several p and q pairs are chosen for further analysis. the difference in the ssds obtained by using the time-frame expanded cssbn technique, the approximate and accurate analysis are shown in table  p p p q q q ||Δ ssd ap - ac ||Δ ssd cssbn - ac ||Δ ssd ap - ac ||Δ ssd cssbn - ac ||Δ ssd ap - ac ||Δ ssd cssbn - ac Δ ssd ap - ac: the difference between the results obtained by using the approximate and accurate analysis. Δ ssd cssbn - ac: the difference between the results obtained by using the cssbn approach and the accurate analysis. p: perturbation rate, q: context switching probability, l: sequence length. as revealed in table  the cssbn approach can compute the ssd more accurately than the approximate analysis. in fact, the difference in the results between the cssbn and the accurate analysis is negligible when reasonably long stochastic sequences are used. with the stm obtained for a cssbn, an ssd is evaluated and the results are very accurate compared to those obtained using the accurate approach, as shown in figure  a further analysis shows that the relative error is limited to less than approximately % for the cssbn approach.the ssds obtained using different methods, as well as the corresponding values of norms, are shown in figure  it can be seen that the ssd can be accurately evaluated by the cssbn model compared to the accurate analytical approach. it has been shown that a pbn with perturbation evolves as an ergodic markov chain. for a larger perturbation probability p, there is an increased randomness in the network activities, thus the steady states of the network are more evenly distributed. this can be seen in the simulation results in figure  figures  and further show that context switching has little impact on the ssd of the network with perturbation. this is due to the fact that the context switching activity only affects the selection probability of a context, but not the predictor functions. also shown in figure  is that an oscillation of the states of p and mdm exists, indicated by the higher probabilities of the network states and compared with the other states. however, this two-state oscillation is less evident when the perturbation rate is significantly higher. on the other hand, a modest perturbation to the network only has a minor effect on the oscillatory behavior, confirming the stability and robustness of the p53-mdm network in a dynamic environment. experiments on a glioma network a network in obtained from a human glioma gene expression data set is further used to illustrate the efficiency of the cssbn model and the time-frame expansion technique. based on this data set, a pbn was inferred by a method using the coefficients of determination and an ssd analysis was performed in. an approximate method was developed in as an extension to estimate the ssd of a cspbn with perturbation. the gene top2a was not considered in either study as it is an input gene with an in-degree of zero. in our study, the network setting is considered the same as in with the gene top2a removed; this leads to a total of gaps. figure  shows a detailed structure of the considered glioma network with double -headed arrows indicating the bi -directional relationships of gene pairs. the selection probabilities for the predictor functions are shown in table  and the boolean functions for each gene are obtained through an analysis of the data in. f i, i ∈ { …, 14}, for each gene in the glioma network in figure f 1 f 2 f 3 f 4 f 5 f 6 f 7 f 8 f 9 f 10 f 11 f 12 f 13 f 14 in table  each column for f i, i ∈ { …, 14}, contains the selection probabilities for the boolean functions of gene i, with the value in the jth row as the probability for fji, thus the probabilities sum to in each column. table  also shows that six genes have only a single predictor function; only the state of gene is determined by three functions, while the states of the other genes are determined by two predictor functions. steady state analysis for the 14-gene glioma network, there are a total of contexts, as can be determined from table  it requires an stm with 214 × 384 =  columns and rows for an accurate analysis. this makes it infeasible to estimate the steady states of a cspbn using a matrix-based analysis. the approximate analysis in would require the computation of an stm of the size 214 ×  thus, it is difficult in general to use an analytical approach to evaluate a large network, due to the demanding computational resources required. however, a cssbn model can be constructed for the glioma network; this cssbn is based on the constituent sbn, as shown in figure  with the cssbn, the ssd can be estimated using the aforementioned time-frame expansion technique with a greatly reduced complexity. the obtained ssd is then compared with that obtained from the approximate analysis in. in this paper, a network is considered to have reached a steady state if the discrepancy between two adjacent iterations is smaller than a required threshold ϵ ) or the number of simulation iterations has reached a maximum constant value. see additional file 4: the matlab program that evaluates the steady state distribution using the time-frame expanded cssbn technique for the glioma network with context switching and perturbation. the state or gap of the glioma network can be represented by a binary vector as, x i ∈ { 1} for i ∈ { …, 14}, or its decimal index. for example, the state can be represented as state the ssds of the context-sensitive glioma network with all states, obtained using the cssbn approach and the approximate analysis, are shown in figure  the norms of the differences of ssds, obtained using the cssbn with different sequence lengths and the approximate method in, are shown in table  p  =  q p  =  q l ||Δ ssd ap - cssbn ||Δ ssd ap - cssbn ||Δ ssd ap - cssbn a maximum of iterations are used for obtaining the ssd. p: perturbation rate; q: context switching probability; l: sequence length. Δ ssd ap - cssbn: the difference between the ssds obtained by the approximate analysis and the context-sensitive stochastic boolean network approach. as shown in table  the cssbn time-frame expansion technique efficiently evaluates the ssd of the glioma network and produces results comparable to those obtained by the approximate analysis. evaluation accuracy is further improved by using longer stochastic sequences with yet a shorter runtime compared to the approximate method. since it is difficult, if not impossible, to compute the stm or ssd of a large grn by using an accurate or approximate analysis, the cssbn time-frame expansion approach provides an alternative means to accurately and efficiently estimate the ssd of a large network with a tunable accuracy by using stochastic sequences of different lengths. intervention analysis in an n-gene network, if gene x i is the target gene, the states of all genes can be divided into a set of desirable states, d, and a set of undesirable states, u, by the expression level of the target gene. as can be seen in figure  gene x the; bfl protein; grs protein is one of the most influential genes that interacts closely with others, thus it is chosen as the target gene. the desirable and undesirable states are then separated by the expression level of x assume that the inactivation of x is preferred; the cumulative probabilities of the desirable and undesirable states are given by ∑x14=0πx and ∑x14=1πy, respectively, where π x and π y are elements in the desirable and undesirable ssds respectively. as discussed previously, a grn can be intervened by applying external stimuli to minimize the likelihood of being in an undesirable state. various methods have been proposed for deriving control vectors for an effective intervention; in this study the ssd algorithm proposed in is used to determine an intervention vector. for simplicity, most significant contexts that account for a total selection probability of % are chosen from the total contexts of the glioma network for an intervention analysis. for these contexts, the expression levels of gene x and x have no effect on the states of other genes, so these two genes are removed; this yields a simplified 12-gene glioma network. the cumulative distributions of the desirable steady states of the context-sensitive 12-gene glioma network are shown in table  for using a different gene as the control gene. p  =  and context switching probability q  =  x 1 x 2 x 3 x 4 x 6 x 8 x 9 x 10 x 11 x 12 x 13 the cumulative probability of the desirable states with no intervention is %. sequence length l = 800k bits. when a different gene is selected as the control gene, the effect of intervention varies with respect to improving the probabilities of desirable states in the ssd. for the glioma network, as revealed in table  gene x tie- is the most effective for maximizing the percentage of the desirable states. the cumulative distribution of the desirable states is increased from % to % by an intervention using tie- as the control gene. also revealed in the table is that an intervention via gene x;, has nearly no impact on the distribution of the desirable states, as the percentage of the desirable states has not been changed much by the intervention. a modest improvement in the desirable state distribution is obtained by an intervention with another gene as the control gene. for any target gene, this process can be applied to find the most significant gene that maximizes the probabilities of desirable states. the effects of intervention can also be seen in figure  where three different scenarios are considered: no intervention, an indirect intervention via gene x i.e., tie- and a direct intervention via the target gene x the; bfl protein; grs protein. the cumulative probabilities of the desirable states are % for no intervention, % and % for the two intervention strategies. as revealed in these results, a direct intervention of the target gene is perhaps optimal for avoiding the undesirable states and almost certainly taking the network into a desirable state. however, when an intervention on the target gene is not possible, an intervention through tie- is the most effective means to maximize the probability of the down-regulation of the target gene x  context-sensitive stochastic boolean networks are proposed as an efficient approach to modeling the effects of gene perturbation and intervention in gene regulatory networks. in a cssbn, the state transition matrix can be accurately and efficiently computed with a complexity of o, where n is the number of genes in a context-sensitive probabilistic boolean network, k is the number of contexts and l is a factor determined by the stochastic sequence length. this result is an improvement compared to the previous result of o for an accurate analysis. the use of non-bernoulli stochastic sequences further increases computational efficiency and allows for a tunable tradeoff between accuracy and efficiency. a steady state analysis using a time-frame expansion technique has shown a significant speedup and produced more accurate results than an approximate analysis in the computation of the steady state distribution. cssbns are constructed for the analysis of gene perturbation in a p53-mdm network and gene intervention in a glioma network. biologically meaningful insights are gained into the oscillatory dynamics of the p53-mdm network in a context-switching environment with random gene perturbation. it has been shown that the steady state distribution changes drastically with the increase of the perturbation rate, whereas the effect of context switching is rather limited for a given perturbation rate. therefore, random gene perturbation may have a greater effect on the final distribution of the steady state compared to context switching activities. by predicting the ssd, the cssbn approach can further evaluate the effectiveness of external gene intervention. a case study of the glioma network shows that the cssbns are useful in modeling the effects of gene perturbation and intervention in a complex context-sensitive grn. this will eventually help drug discovery for an effective drug intervention therapy. competing interests the authors declare that they have no competing interests. authors’ contributions pz carried out the grn studies, performed the statistical analysis and drafted the manuscript. jl developed the initial models and participated in discussion. jh participated in the grn studies, revised the manuscript and supervised the research. all authors read and approved the final manuscript. supplementary material additional file 1 the pseudocode for computing the state transition matrix of a context-sensitive stochastic boolean network. click here for file additional file 2 the pseudocode for computing the steady-state distribution of a context-sensitive stochastic boolean network. click here for file additional file 3 the matlab program that describes the structure of the context-sensitive stochastic boolean network in figure 6and computes its state transition matrix for the p53-mdm network with perturbation. click here for file additional file 4 the matlab program that evaluates the steady state distribution using the time-frame expanded context-sensitive stochastic boolean network technique for the glioma network with context switching and perturbation. click here for file acknowledgements this work was supported in part by the natural sciences and engineering research council of canada in a discovery grant and a startup fund of the university of alberta.
a cssbn enables an efficient simulation of a cspbn with a complexity of o for computing the state transition matrix, where l is a factor related to the required sequence length in cssbn for achieving a desired accuracy. a time-frame expanded cssbn can further efficiently simulate the stationary behavior of a cspbn and allow for a tunable tradeoff between accuracy and efficiency. a cssbn analysis provides biologically meaningful insights into the oscillatory dynamics of the p53-mdm network in a context-switching environment. context-sensitive stochastic boolean networks are proposed as an efficient approach to modeling the effects of gene perturbation and intervention in gene regulatory networks. this yields a novel structure of context-sensitive sbns for modeling the stochasticity in a grn. gene regulatory networksboolean networksstochastic boolean networkscontext dependentgene perturbationinterventioncontext switchsteady state distributionp networkglioma network in a gene regulatory network, gene expressions are affected by noise, and stochastic fluctuations exist in the interactions among genes. in a cspbn with n genes and k contexts, however, a computational complexity of o ) is required for an accurate computation of the state transition matrix of the size ×. it is shown that random gene perturbation has a greater effect on the final distribution of the steady state of a network compared to context switching activities. the cssbn approach can further predict the steady state distribution of a glioma network under gene intervention. the evaluation of a steady state distribution is more challenging. these stochastic interactions are context dependent, thus it becomes important to consider noise in a context-sensitive manner in a network model. the notion of stochastic boolean networks is extended for the general model of pbns, i.e., cspbns. ultimately, this will help drug discovery and develop effective drug intervention strategies. the cssbn approach is more efficient than an analytical method and more accurate than an approximate analysis. recently, stochastic boolean networks have been proposed as an efficient implementation of an instantaneous pbn.
 in a gene regulatory network, gene expressions are affected by noise, and stochastic fluctuations exist in the interactions among genes. these stochastic interactions are context dependent, thus it becomes important to consider noise in a context-sensitive manner in a network model. as a logical model, context-sensitive probabilistic boolean networks account for molecular and genetic noise in the temporal context of gene functions. in a cspbn with n genes and k contexts, however, a computational complexity of o ) is required for an accurate computation of the state transition matrix of the size ×. the evaluation of a steady state distribution is more challenging. recently, stochastic boolean networks have been proposed as an efficient implementation of an instantaneous pbn.  the notion of stochastic boolean networks is extended for the general model of pbns, i.e., cspbns. this yields a novel structure of context-sensitive sbns for modeling the stochasticity in a grn. a cssbn enables an efficient simulation of a cspbn with a complexity of o for computing the state transition matrix, where l is a factor related to the required sequence length in cssbn for achieving a desired accuracy. a time-frame expanded cssbn can further efficiently simulate the stationary behavior of a cspbn and allow for a tunable tradeoff between accuracy and efficiency. the cssbn approach is more efficient than an analytical method and more accurate than an approximate analysis.  context-sensitive stochastic boolean networks are proposed as an efficient approach to modeling the effects of gene perturbation and intervention in gene regulatory networks. a cssbn analysis provides biologically meaningful insights into the oscillatory dynamics of the p53-mdm network in a context-switching environment. it is shown that random gene perturbation has a greater effect on the final distribution of the steady state of a network compared to context switching activities. the cssbn approach can further predict the steady state distribution of a glioma network under gene intervention. ultimately, this will help drug discovery and develop effective drug intervention strategies. gene regulatory networksboolean networksstochastic boolean networkscontext dependentgene perturbationinterventioncontext switchsteady state distributionp networkglioma network 