BACKGROUND
a key area in the field of translational genomics is to derive therapeutic intervention that can beneficially alter cell dynamics in such a way as to avoid cancerous phenotypes.
the first step to derive therapeutic interventions is to understand the regulatory relationships among genes.
the interactions among genes are studied in the context of gene regulatory networks .
grn models often possess high uncertainty.
this inherent uncertainty might be due to many factors such as the complex nature of biological phenomena, lack of enough training data, etc.
uncertainty in grns can affect the accuracy and performance of the therapeutic interventions.
therefore, biologists aim at reducing uncertainty of grn model via conducting additional biological experiments.
however, there are many limitations for carrying out experiments.
biological experiments are usually expensive and time consuming as they need to be done on living organisms.
on the other hand, the resources are limited which prevents researchers from conducting all experiments they need for identification of grn.
thus, it is prudent to prioritize potential experiments based on the information they provide and then conduct the most informative.
this process is called experimental design.
boolean networks  have been extensively used for the study of grns.
they have been shown to be effective in capturing the multivariate relationships among entities within a cell.
bns also facilitate the use of the well-studied markov decision theory for deriving beneficial interventions.
to date, intervention methods for bns have been categorized in two different groups: structural interventions, and dynamical interventions.
while structural interventions  <cit>  permanently change the behavior of a network via one-time change in the underlying regulatory structure, the goal in dynamical interventions  <cit>  is to interfere with the signaling of grn through flipping or not flipping the expression state of genes over time.
the assumption behind most of these intervention methods is that the network being studied is perfectly known.
thus, presence of uncertainty can degrade the performance of the intervention methods.
from a translational perspective it is crucial to reduce uncertainty pertinent to the objective, such as intervention.
hence, uncertainty should be quantified in such a way that the objective for modeling grn is taken into account.
mean objective cost of uncertainty   <cit>  quantifies the uncertainty of model in terms of the expected increased cost due to the presence of uncertainty.
an experimental design method based on mocu has been proposed in  <cit> .
the experimental design method in  <cit>  evaluates the effect of each potential experiment in reducing the model uncertainty, which is measured in terms of mocu, and suggests that the experiment which results in the minimum expected remaining mocu should be conducted first.
the long-run performance of the experimental design method in  <cit>  is guaranteed to be optimal in terms of reducing the error of interventions obtained after conducting the chosen experiments.
although, the method in  <cit>  is optimal, it is computationally expensive.
because our final objective is to improve the performance of the therapeutic interventions, method in  <cit>  involves finding optimal interventions for all networks which are compatible with the prior knowledge.
finding optimal interventions is computationally expensive whose complexity grows exponentially with the number of genes in network.
therefore, the computational complexity of finding optimal experiment can be prohibitively high for large networks.
thus, it is inevitable to construct a smaller network via deleting some genes from the original large size network and then estimate the optimal interventions using the resulting reduced network.
generally the goal in network reduction methods is to produce networks of smaller size while the dynamical behavior of the original network is preserved.
there have been some efforts for network reduction to reduce the complexity of designing interventions  <cit> .
in this paper, we propose a novel cost function for the gene deletion process which takes into account the disruption in the order of potential experiments when they are ranked according to the experimental design method in  <cit> .
since experiments are ranked based upon the expected remaining mocu or the mocu that is expected to remain after performing the experiment, we desire that the network reduction step has a low effect on the expected remaining mocu corresponding to the potential experiments.
when the gene  suggested by the cost function are deleted from network, the optimal  interventions are estimated using the reduced networks and then they are used for calculating expected remaining mocu for prioritizing potential experiments.
we show the effectiveness of our proposed cost effective experimental design method through simulations on synthetic and real networks.
the simulation results verify that our method can perform comparable to the optimal experimental design method in  <cit>  with much lower computations.
mocu-based optimal experimental design is very general and does not even require a markovian network  <cit> .
as we will see, finding the best gene to delete is also very general; however, once the genes are deleted, the regulatory structure of the original network must be mapped onto a corresponding regulatory structure on the reduced network, an optimal intervention must be found on the reduced network, and that intervention must be induced to the full network.
reduction and inducement are nontrivial and depend on the nature of the regulatory structure.
the problem has been addressed for boolean networks in  <cit> , to which we refer, and a theoretical analysis is given in  <cit> , where it is noted that the methodology applies to probabilistic boolean networks   <cit>  by applying the reduction to each of constituent bn of the pbn.
moreover, whereas we will restrict intervention to rank-one perturbations  <cit> , which provide a one-time alteration of the regulatory logic, the reduction-inducement paradigm applies to other forms of intervention  <cit> .
methods
boolean networks
gene regulatory network models are increasingly used as a tool to study interactions among genes  <cit> .
boolean networks   <cit>  and probabilistic boolean networks   <cit>  are widely used models for grns that have been shown to be effective in capturing these interactions  <cit> .
a boolean network on n genes is defined by a pair .
v = {x <dig>  x <dig> ..., xn} is a set of binary variables that represent the expression state of genes with xi =  <dig> and xi =  <dig> corresponding to gene i being off and on, respectively.
the on state means that a gene is translated into a protein whereas the off state means that the gene is not translated.
it has been shown that significant biological information can be extracted from binarized gene expression data  <cit> .
the gene values at each time step are updated according to the list of boolean functions f = {f <dig>  f <dig> ..., fn} with xi = fi where xi <dig>  xi <dig> ..., xiki are ki predictor genes for xi.
the state of the bn at time t is the vector x = , x <dig> ..., xn), which is called gene activity profile .
a bn with n genes possesses 2n different states.
the state of a bn for the next time instant t +  <dig> is updated as x = f).
in a boolean network with perturbation , there is a perturbation probability of each gene flipping its value independently from other genes.
thus, in the presence of perturbation x = f) with probability n and with probability  <dig> − n at least one gene flips its value.
the sequence of states of a bnp over time can be viewed as the states of a markov chain with transition probability matrix  p = piji,j=12n where pij is the probability that state i transitions into state j. positive perturbation probability p > <dig> makes the markov chain ergodic and irreducible, thereby possessing a unique steady-state distribution  πt = πt p, where πi is the steady-state probability of state i and t is the transpose operator.
more details on how to compute tpm can be found in  <cit> .
the long-run behavior of a bnp is characterized by its ssd, which is indicative of phenotypes.
from a translational perspective, the states of a bnp can be partitioned into the set of desirable states d containing those states associated with healthy phenotypes and undesirable states u containing states associated with cancerous phenotypes.
the goal of therapeutic interventions is to drive the network away from undesirable states and consequently reduce the steady-state probability mass of undesirable states, πu = ∑uπi.
without loss of generality regarding the experimental design analysis, in this paper we will assume that the undesirable states are determined by a single gene, called the "target gene" with the aim of intervention being to flip an undesirable value of the target gene to a desirable value.
two intervention approaches are commonly considered:  structural interventions change the long-run behavior of the network by altering its underlying rule-based structure  <cit> , and  dynamical interventions affect the dynamical evolution of the network by manipulating the expression states of some "control genes"  <cit> .
this paper considers optimal structural intervention via rank-one function perturbations  <cit> .
in a rankone function perturbation intervention, the relation between the tpm before intervention, p, and after intervention, p˜, is p˜=p+abt where a and b are two arbitrary vectors satisfying bt e =  <dig>  where e is an all unity column vector.
single-gene perturbation, a special case of rank-one function perturbation, changes the regulatory function f for only one state.
if originally f = w, then the single-gene perturbation  means that for the regulatory function f˜after intervention, f˜=v and f˜=f for i ≠ u. for this intervention, all entries of p and p˜ would be identical except for two entries.
p˜uw=puw-n and p˜uv=puv+n.
the steady-state probability of state i after intervention  is found using the following equation:
π˜i=πi+nπu1-n,
where π˜i is the steady-state probability mass of state i after single-gene perturbation  and zvi, zwi, zvu, and zwu belong to the fundamental matrix z = − <dig> where i is the n × n identity matrix.
the undesirable steady state probability mass after single-gene function perturbation  is π˜u=∑i∈uπ˜i.
to find the optimal intervention, one needs to search among all possible 2n × 2n state pairs .
optimal experimental design
this section briefly reviews the optimal experimental design method in  <cit> .
suppose θ <dig>  θ <dig> ..., θk are k uncertain parameters in the network model.
they compose an uncertainty vector θ = , the set of all such vectors being denoted Θ. the collection of all networks corresponding to θ ∈ Θ constitute an uncertainty class, denoted by Θ. let Ψ be a class of potential interventions .
let ξθ  denote the error of applying intervention ψ ∈ Ψ to the network with uncertainty vector θ ∈ Θ. in case of intervention, ξθ  can  be the steady-state probability mass of undesirable states after intervention, denoted π˜u,θ.
for a network with uncertainty vector θ, the optimal intervention ψ is defined by
ψ=argminψ∈Ψξθ.
in the presence of uncertainty, we desire an intervention that performs optimally on average across the uncertainty class.
the intrinsic bayesian robust  intervention ψibr is defined by
ψibr=argminψ∈Ψeθξθ,
where eθ  denotes an expectation relative to the probability density function f  governing the uncertainty class Θ. note that although ψibr has the best expected performance across the uncertainty class, it is not guaranteed to have the best performance for each network inside the uncertainty class.
the mean objective cost of uncertainty  relative to Θ and Ψ is defined by  <cit>
mΨ=eθ,
the expression inside the expectation being called the objective cost of uncertainty.
mocu is the expected difference between the performance of the robust intervention and the optimal intervention for each network inside Θ. mocu is an uncertainty quantification scheme which measures uncertainty by taking into account the objective of operator design, herein, network intervention.
mocu can effectively capture the amount of uncertainty and tries to quantify the uncertainty based on a specific objective  to gauge how much it is going to affect the objective  <cit> .
mocu has been used to design experiments that can maximally enhance the performance of the interventions  <cit> .
suppose that there are k potential experiments e <dig>  e <dig> ..., ek corresponding to k uncertain parameters θ <dig>  θ <dig> ..., θk .
it is assumed that experiment ei, which might be a complex experiment consisting of several sub-experiments, fully identifies θi.
the goal of experimental design is to find out which experiment ei,  <dig> ≤ i ≤ k should be conducted first, or how to rank potential experiments effectively.
although we do not know the actual outcome of an experiment before conducting it, we do know the possible outcomes.
suppose the outcome of experiment ei is ϕ. we define the remaining mocu, given θi = ϕ, as
mψ=eθϕξθϕ)-ξθϕ)),
where θϕ, called the conditional uncertainty vector, has ith parameter equal to  ϕ and other parameters are still random, and Θi,ϕ is the remaining uncertainty class when θi = ϕ, i.e., Θi,ϕ = {θ|θ ∈ Θ, θi = ϕ}.
the expectation is taken over the conditional density function f)=f that governs the remaining uncertainty class Θi,ϕ.
ψibr is the robust intervention for Θi,ϕ and is found in a similar way as .
corresponding to each experiment ei, we define the expected remaining mocu after conducting that experiment as:
mΨ=eθimΨ=eθieθθiξθθi)-ξθθi)),
where eθi is the expectation relative to the marginal density function f of uncertain parameter θi.
then we find the minimum expected remaining mocu:
i*=argmini= <dig> ,...,kmΨ.
ei* is the optimal experiment to be conducted first  <cit> .
approximate experimental design method
according to  and , calculating the expected remaining mocu requires finding the optimal intervention ψ for each θ ∈ Θ and the robust intervention ψibr for each possible remaining uncertainty class.
the complexity of finding optimal interventions grows exponentially with network size n. for finding an optimal single-gene structural intervention, we need to search among all possible 2n × 2n state pairs and calculate the new steady-state probability π˜i for each state i in the set of undesirable states u. thus, the complexity is  o. this heavy computational cost motivates us to reduce the size of network in order to reduce the complexity of finding optimal interventions, thereby reducing the complexity of the experimental design.
assuming that gene g is deleted from a network with regulatory function f, we define a regulatory function fred for the reduced network.
doing this for each network with uncertainty vector θ in Θ produces the uncertainty class, Θg, of reduced networks via the mapping θ → θg.
to approximate the optimal intervention for a network in Θ, we use the corresponding network in Θg, find the optimal intervention for the reduced network ψ, and then induce the intervention to the original network in Θ. this approximate optimal intervention denoted by ψ is called the induced optimal intervention.
also, to find the induced robust intervention, ψibrind, for Θ, first we find the robust intervention, ψibr, for Θg using  and then find the induced robust intervention ψibrind from ψibr.
as illustrated in figure  <dig>  in the proposed approximate experimental design method, we find the best gene g∗ for deletion via a novel cost function c and then obtain the induced optimal and robust interventions needed for the mocu calculations in the experimental design step by inducing interventions from uncertainty class of reduced networks Θg∗ to the original uncertainty class Θ.
we now aim to find a gene whose deletion results in minimum degradation in the experimental design process.
keeping in mind that the experimental design is based on the expected remaining mocu for potential experiments, let mΨg be the remaining mocu when uncertain parameter i has value θi and we delete gene g,
mΨg=eθθiξθθi)-ξθθi)).
we define the cost of deleting gene g by
c=∑i=1k|mΨg-mΨ|,
where
mΨg=eθimΨg.
the gene g∗ minimizing the cost function in  is selected for deletion:
g*=argming∈ <dig> ,...,nc.
the intuition behind this cost function is that our choice of optimal experiment is based upon the expected remaining mocu corresponding to each experiment.
therefore, we desire that the network reduction step has minimum effect on these quantities.
deleting genes increases the inherent uncertainty of the network because the induced robust intervention cannot perform better than the original robust intervention on average.
we want to reduce this increase in the uncertainty of model caused by network reduction.
since
eθθiξθθi)≥eθθiξθθi),
mΨg≥mΨ.
hence, we can omit the absolute value operator in  to obtain
g*=argming∑i=1k-mΨ)=argming∑i=1kmΨg,
where the second equality follows from the fact that mΨ does not depend on the gene being deleted.
expanding mΨg yields
g*=argming∑i=1keθieθθiξθθi)-ξθθi))=argming∑i=1keθieθθiξθθi).
the minimization problem in  is equivalent to the one in .
based on the cost function in , for each gene g, we find the expected performance of the induced robust intervention ψibrind across the remaining uncertainty class Θi,θi, then take the expectation of this average performance relative to the marginal distribution of the uncertain parameter θi, and finally sum all values found for the k uncertain parameters.
algorithm  <dig> approximate experimental design
1: input: Θ, Ψ, f , θ =
2: output: ei*, i∗ ∈ { <dig> ..., k}: the estimated optimal experiment to be conducted first
3: for g =  <dig> : n do
4:   cost ← 0
5:   for i =  <dig> : k do
6:      for all θi do
7:         build remaining uncertainty class of reduced networks Θi,θig
8:         compute conditional density function f)
9:         find induced robust intervention Ψibrind
10:         hg←eθθiξθθi)
11:      cost ← cost + eθi
12: g*←argming= <dig> ,...,ncost
13: for i =  <dig> : k do
14:   for all θi do
15:      build remaining uncertainty class Θi,θi
16:      compute conditional density function f)
17:      compute mΨg* via equation  using Ψibrind and Ψind;g*)
18:      mΨg*←eθimΨg*
19: i*←argmini= <dig> ,...,kmΨg*
20: return i∗
after removing gene g∗, we find the expected remaining mocu corresponding to each experiment using equation  by replacing ψ) with ψind;g*) and ψibr with ψibrind.
an abstract form of the proposed experimental design method has been given in algorithm  <dig>  a step by step toy example illustrating algorithm  <dig> is also provided in the additional file  <dig>
this procedure for estimating optimal experiment via deleting one gene can be easily extended to the deletion of two or more genes.
for example, to delete two genes, we need to evaluate the cost function in  for all possible two-gene combinations and delete the pair whose cost is minimum.
reduction mappings and induced interventions
if we want to delete gene g from network, we need to find the regulatory function fred for the reduced network.
following  <cit> , every two states of the original network that differ only in the value of gene g can be collapsed to find the transition rule of the reduced network.
let s1g and s0g be two states with value  <dig> and  <dig> for gene g, respectively, and identical values for other genes.
state sg¯ can be obtained from either s1g or s0g by removing the value of gene g. if for the original network, the transition rules for these two states are f = p and f = q, then for the reduced network, fred=pg¯ if πs1g>πs0g and otherwise fred=qg¯, where pg¯ and qg¯ are found from states p and q via removing the value of gene g, respectively.
following this procedure, we find the regulatory function fred for all states in the reduced network.
algorithm  <dig> finding induced optimal interventions
1: input: ψ=,π
2: output: ψind=
3: u˜1g← place  <dig> in the gth coordinate of u˜
4: u˜0g← place  <dig> in the gth coordinate of u˜
5: if πu˜1g≥πu˜0g then
6:   ugind←u˜1g
7: else
8:   ugind←u˜0g
9: υ˜1g← place  <dig> in the gth coordinate of υ^
10: υ˜0g← place  <dig> in the gth coordinate of υ^
11: if πv^1g≥πv^0g then
12:   υgind←υ˜1g
13: else
14:   υgind←υ˜0g
as illustrated in algorithm  <dig>  we find the induced optimal intervention from the optimal intervention for the reduced network.
suppose that the optimal intervention for the reduced network θg is ψ=.
the two corresponding states to u^ in the original network are u˜1g and u˜0g, which are found by placing  <dig> and  <dig> in the gth coordinate of u^, respectively.
similarly, there are two states v˜1g and v˜0g in the original network corresponding to state υ^ the induced optimal intervention for the original network is ψind=, where ugind is the one among u˜1g and u˜0g having larger steady-state probability in the original network and vgind is the one among υ˜1g and υ˜0g with larger steady-state probability in the original network.
analogous to the induced optimal intervention, the induced robust intervention ψibrind is found from the robust intervention ψibr according to algorithm 3; however, here we choose the two states possessing larger expected steady-state probability across Θ using the expected ssd, π = eθ , where π is the ssd of the network with uncertainty vector θ in uncertainty class Θ. we can use this procedure to find the induced robust intervention for each remaining uncertainty class Θi,ϕ.
algorithm  <dig> finding induced robust interventions
1: input: ψibr=,π∀θ∈Θ
2: output: ψibrind=
3: π←eθ
4: u˜1g← placing  <dig> in the gth coordinate of u^
5: u˜0g← placing  <dig> in the gth coordinate of u^
6: if πu˜1g≥πu˜0g then
7:   ugind←u˜1g
8: else
9:   ugind←u˜0g
10: v˜1g← placing  <dig> in the gth coordinate of υ^
11: v˜0g← placing  <dig> in the gth coordinate of υ^
12: if πυ˜1g≥πυ˜0g then
13:   υgind←υ˜1g
14: else
15:   υgind←υ˜0g
preliminary gene elimination via the coefficient of determination
to further reduce the computational cost of the experimental design, we utilize the coefficient of determination   <cit>  to eliminate some genes from the optimization problem without evaluating the cost function and then search among the remaining genes for choosing genes to be removed using the cost function .
the cod measures the strength of relationship between a target gene y and a vector x of predictor genes as the difference between the error of the best estimation of gene y in the absence of other genes and in the presence of genes in x. the cod is between  <dig> and  <dig> and a larger cod means a stronger connection between the target and predictor genes, in our case the target gene being the aim of intervention.
we use the intuition that genes possessing large cod in relation to the target gene are not likely among the genes that should be deleted because they have strong connection to the target gene.
the cod of the target gene y, relative to a vector x =  of predictor genes is defined by
codx=εy-εx,yεy
where εy is the error of the best estimation of y without any predictors,
εy=min,
and εx,y is the error of the optimal estimation of y upon observing x. by assuming that the value of the binary vector x of predictor genes changes from  <dig> to 2m, εx,y can be calculated by
εx,y=∑j=12mpmin.
if codx denotes the cod of y relative to x in a network with uncertainty vector θ, then given the uncertainty class Θ the expected cod of y relative to x is given by
codx=eθ.
genes possessing strong connection with the target gene in terms of codx are not considered for deletion.
when excluding genes using the cod it is important to recognize the possibility of intrinsic multivariate prediction  <cit> , where a set of genes may have low individual cods with respect to the target gene but may have significant cod when used together for multivariate prediction.
first we calculate codx for all 3-gene combinations and pick the one with largest cod.
we compute cod for 3-gene predictors because it has been shown in  <cit>  that the average connectivity of the model cannot be too high providing that the model is not chaotic and it is commonplace to assume 3-gene predictivity in bns.
if we want to exclude less than  <dig> genes from the search space, then among the 3-gene combination with the largest expected cod, we choose those genes that have larger expected individual cod.
if we want to exclude more than  <dig> genes, then in addition to the three genes in the combination with the largest cod, we choose those genes in the 3-gene combination with the second largest cod that have larger expected individual cod and do not belong to the first 3-gene combination.
we repeat this process until we reach the desired number of genes to exclude.
if there are initially n genes and we want to delete  <dig> genes, then we need to evaluate cost function  for all c 3-gene combinations, where c denotes the number of combinations of n objects taken k at a time; however, if we exclude s genes from search space then the number of evaluations of  decreases to c.
having performed the cod-based exclusion process and excluded s genes, x1′,x2′,...,xs′, we search for the genes to be deleted using the cost function in  among the remaining genes, {x <dig>  x <dig> ..., xn} − {x1′,x2′,...,xs′}.
computational complexity analysis
the first step for the optimal experimental design in  <cit>  is estimating optimal interventions ψ for each network in Θ. we also need to compute the robust intervention ψibr for each possible remaining uncertainty class Θi,θi.
most of the computations are devoted to this step.
finding robust interventions does not require additional calculations because we can store the error of each intervention ψ∈Ψ for the network θ when finding optimal interventions and later use these errors to find robust interventions.
therefore, complexity analysis requires computing the complexity of estimating the optimal interventions.
with n genes that take on binary expression levels, the network has 2n states.
finding an optimal single-gene function intervention requires searching among all possible 22n state pairs  according to .
assuming without loss of generality that states 2n− <dig> to 2n are undesirable,  must be evaluated 2n− <dig> times for each state pair.
thus, the complexity of finding the optimal intervention ψ is  o. if there are k uncertain parameters and each can take on l different values, then the uncertainty class Θ contains lk different networks for which an optimal intervention must be found.
hence, the complexity of the optimal experimental design method in  <cit>  is  o.
to analyze the complexity of the proposed approximate method, suppose p genes are to be deleted.
then the cost function in  must be evaluated for all c p-gene combinations, n −  <dig> instead of n because the target gene cannot be deleted.
the complexity of finding an induced optimal intervention for each network after deleting p genes is  o).
therefore, the complexity of the approximate method is  o × lk × 23n−3p).
for large n, it is possible that for small p the complexity of the approximate method can exceed that of the original method; however, by deleting more genes the complexity of the approximate method drops sharply because by deleting each additional gene the complexity of estimating the optimal intervention decreases by eight-fold.
by incorporating the cod-based gene exclusion step in the approximation method and excluding s genes we are able to decrease the number of p-gene combinations from c to c, which reduces the complexity of the approximate method to  o × lk × 23n−3p).
define the computational gain λ by
λ=lk×23nc×lk×23n-3p=23pc,
which is the ratio of the complexity of the optimal method in  <cit>  to the complexity of the approximate method when deleting p genes using the cost function in  and excluding s genes from the search space using the cod-based gene exclusion step.
RESULTS
this section evaluates the performance of the approximate method for both synthetic and real grns where the majority vote rule is used as the transition rule.
majority vote rule  <cit>  is popular in systems biology, especially when we are interested in the overall dynamics of the network.
for example, majority vote is used in  <cit>  to model the dynamics of yeast cell-cycle network.
for the majority vote rule, a regulatory matrix r is defined component-wise by
rij=1genejactivatesgenei-1genejsuppressesgenei0norelationfromgenejtogenei.
according to this rule, gene i takes value  <dig> if the number of genes that are on and activate it is more than the number of genes that are on and suppress it:
xi=fi)=1if∑jrijxj>00if∑jrijxj<0xiif∑jrijxj= <dig>
uncertainty is introduced by assuming that the exact values of some of the nonzero components of r are unknown; that is, for some regulations it is not known whether they are activating or suppressive.
each uncertain parameter θi can be − <dig> or  <dig>  conducting experiment ei determines the value of parameter θi.
let μ =  denote the true value for the uncertainty vector θ = .
conducting experiment ei results in a remaining uncertainty class Θi,μi consisting of networks with θi = μi and other uncertain parameters being − <dig> or  <dig>  for Θi,μi we can determine a robust intervention ψibr.
we evaluate the effectiveness of experiment ei in terms of the error of the resulting robust intervention obtained after experiment on the underlying true network, ξμ).
we define the gain of conducting the chosen experiment ei* over a random experiment ernd  by
ρ=ξμ)-ξμ).
if ρ >  <dig>  then the chosen experiment outperforms the random experiment; if ρ <  <dig>  then the random experiment outperforms the chosen experiment; and if ρ =  <dig>  then they perform the same.
simulation results based on the synthetic bnps
for the performance evaluation based on synthetic bnps, we generated  <dig> networks randomly and chose  <dig> different sets of k regulations in each to be unknown -  <dig> simulations in total.
we assigned  <dig> random predictor genes to each gene where each one can be randomly activating or suppressive.
the gene perturbation probability was set to  <dig> .
without loss of generality, we assume that states with up-regulated x <dig> are undesirable.
we removed the regulatory type of those regulations that have been assumed to be uncertain and retained other regulatory information of the network.
we assume that all uncertain parameters are independent from each other and have uniform marginal distribution.
the analysis can be easily extended to other distributions.
because x <dig> is the target gene, it was excluded from the reduction process.
hence, we look for the best p-gene set to be deleted among {x <dig> ..., xn}.
to consider larger networks, we generated  <dig> random 13-gene networks and in each chose one set of  <dig> regulations to be unknown.
we used the approximate experimental design method and deleted  <dig> genes.
for this size of network it is not possible to perform the optimal experimental design method or compute original optimal and robust interventions to calculate the gain of the chosen experiment over a randomly selected experiment.
hence, we use the induced robust intervention obtained by deleting the set of  <dig> genes having the minimum cost of deletion as in .
therefore, here gain ρ is defined as
ρ=ξμ;g*)-ξμ;g*)
where g* is the set of  <dig> genes with minimum cost of deletion.
for this set of simulations, the average gain ρ¯ is  <dig> .
note that here the average gain might not be very accurate owing to the small number of simulations.
the approximate run time for each simulation was around  <dig> seconds.
we now evaluate the proportion of times that we obtain the optimal experiment found by  <cit>  when using the approximate method.
figure  <dig> shows the percentage of finding the optimal experiment when using the approximate method and deleting different number of genes from 7-gene networks.
in this figure, there are  <dig>   <dig>  and  <dig> values for deleting  <dig>   <dig>  or  <dig> genes, respectively.
we observe that deleting the set of genes corresponding to the minimum of the cost function yields the highest likelihood of obtaining an optimal experiment, which is what we would hope for from an efficient approximate method.
according to figure  <dig>  when we delete the gene which attains the minimum cost,  <dig> % of the simulations yield an optimal experiment, whereas this percentage is  <dig> % when deleting the gene with the largest value of the cost function.
similar behavior is observed when deleting  <dig> or  <dig> genes.
a salient reason that the largest average gain of the approximate method over random experiments is when we delete genes corresponding to the minimum cost function is that it is more likely to get an optimal experiment.
an issue that arises when evaluating experimental design on synthetic networks, as opposed to real biological networks, which typically manifest substantial controllability on account of their need to maintain functionality within changing contexts, a large portion of randomly generated networks might not be controllable and therefore not be responsive to intervention.
hence, intervention has negligible effect on their ssds and including them in the analysis masks the effect of optimality.
to address this issue, we define controllability Δ as the percentage decrease of undesirable probability mass after applying optimal intervention:
Δ=πu-π˜uπu×100%,
where πu and π˜u are the undesirable probability masses before and after applying optimal intervention to the network, respectively.
a larger Δ means that a network is more controllable.
figure  <dig> considers the effect of controllability on the performance of experimental design when networks have n =  <dig> genes and k =  <dig> uncertain regulations.
the figure shows the average gain ρ¯ for the optimal method and the proposed approximate method for networks possessing controllability greater than a given threshold.
we observe that ρ¯ increases when networks are more controllable, regardless of the number of genes deleted from network.
note that as controllability increases, the difference between the performance of different methods increases.
for example, for all networks the average gain for the optimal method and the proposed method when deleting one, two, three, and four genes is  <dig> ,  <dig> ,  <dig> ,  <dig> , and  <dig> , respectively; but for networks with Δ ≥ 40% the average gains are  <dig> ,  <dig> ,  <dig> ,  <dig> , and  <dig> , respectively.
to evaluate the effectiveness of the cod-based gene exclusion algorithm, we compare the average gain of the approximate method when excluding genes from the search space using the cod-based exclusion algorithm against the average gain when excluding randomly selected genes from the search space.
figure  <dig> shows the average gain ρ¯ for networks with n =  <dig> genes and k  =  <dig>   <dig>   <dig>   <dig> uncertain regulations.
for deleting p genes, we exclude up to 6− p − <dig> genes from the search space so that for the largest number of genes excluded, the search space contains at least p +  <dig> genes.
for example, when deleting p =  <dig> gene, we exclude  <dig>   <dig>   <dig>  and  <dig> genes; when deleting p =  <dig> genes we exclude  <dig>   <dig>  and  <dig> genes; and so on.
for each number of uncertain regulations, we observe that the average gain when excluding genes using the cod-based algorithm is always larger than random gene exclusion, regardless of the number of deleted genes.
for example, when k =  <dig>  for deleting one gene and excluding  <dig>   <dig>   <dig>  and  <dig> genes randomly, ρ¯= <dig> ,  <dig> ,  <dig> , and  <dig>  respectively but using the cod-based scheme and excluding the same number of genes, ρ¯= <dig> ,  <dig> ,  <dig> , and  <dig>  respectively.
if we delete three genes ρ¯= <dig>  without gene exclusion, and if we exclude  <dig> and  <dig> genes, then ρ¯= <dig>  and ρ¯= <dig> , respectively, when we exclude genes randomly and ρ¯= <dig>  and ρ¯= <dig> , respectively, when we exclude genes based on cod.
note that when deleting more genes, the difference between random exclusion and cod-based exclusion increases because as more genes are deleted, exclusion has a larger impact on the number of candidate sets for evaluating the cost function.
for example, when deleting  <dig> gene, if we exclude one gene, then the number of candidate sets decreases from  <dig> to  <dig>  but when deleting  <dig> genes, if we exclude one gene, then the number of candidate sets decreases from c =  <dig> to c =  <dig>
in figure  <dig>  we also show the box plot for the gain of conducting the chosen experiment if we delete  <dig> genes and exclude genes from the search space either randomly or via the proposed cod-based method for  <dig> gene networks possessing  <dig> uncertain regulations.
we observe that the first quartile, median, and third quartiles are higher when excluding genes using cod.
for example, when randomly excluding  <dig> genes from the search space, the first quartile, median, and third quartile are − <dig>  × 10− <dig>   <dig>  × 10− <dig>  and  <dig> , respectively; however, when excluding genes using cod they are − <dig>  × 10− <dig>   <dig>  × 10− <dig>  and  <dig> , respectively.
we have also evaluated the performance of the approximate experimental design method when a sequence of experiments is conducted.
suppose there are k =  <dig> uncertain regulations and we conduct five experiments to identify all unknown regulations.
for each set of unknown regulations, at each step we utilize the experimental design to choose one of the possible experiments, conduct the chosen experiment, and measure the performance  of the robust intervention obtained after the experiment on the underlying true network.
continuing, for the remaining uncertain regulations we use the experimental design method and repeat the previous procedure until there is no more unknown regulation remaining in the network.
we also do sequential experiments randomly where at each step we choose an experiment randomly, measure the performance of its corresponding robust intervention, and again choose one experiment randomly among the remaining ones.
the difference between the undesirable probability mass after applying the robust interventions corresponding to the randomly chosen experiment and the chosen experiment through experimental design at each step is the gain of conducting the chosen experiment at that step.
figure  <dig> shows the average gain over random selection for the optimal method and the approximate design method deleting up to three genes for 7-gene networks and up to four genes for 8-gene networks.
the figure indicates that the approximate design method has reliable performance compared to the optimal method.
moreover, similar to the optimal method, the average gain increases sharply in the beginning for the approximate method.
this is very important in real applications owing to the cost and time required for conducting experiments.
note that when we conduct all five experiments the average gain is zero because after five experiments the network is fully known and we can exactly calculate the optimal intervention regardless of the approach taken to choose experiments.
performance evaluation based on the colon cancer pathways
in this section, we analyze the performance of the proposed experimental design method on the colon cancer pathways used in  <cit> .
we focus on the pathways formed by  <dig> genes extracted from the complete pathway set, as used in  <cit> .
these are shown in figure 12: stat <dig>  ras, il <dig>  hgf, pik3ca, egf, tsc1/tsc <dig>  mtor, spry <dig>  pkc, and mek 1/ <dig>  normal and blunt arrows represent activating and suppressive regulations, respectively.
we modeled the pathways as a bnp with perturbation probability  <dig> .
genes are named as they have been introduced.
for example, stat <dig> is x <dig> and mek 1/ <dig> is x <dig>
egf, hgf, and il <dig> are three stimulation factors that carry the external signals generated by neighboring cells to downstream genes and activate downstream pathways.
signal transducers and activators of transcription  constitute a family of transcription factors that can be activated via extracellular signaling proteins such as cytokins and growth factors.
these play a major role in regulating downstream processes such as cell growth, survival, and apoptosis  <cit> .
stat <dig> is an oncogene observed to be highly activated in many cancers, in particular, colon cancer  <cit> .
hence, stat <dig> has been recognized as a legitimate target for cancer therapy  <cit> .
we considered states with up-regulated stat <dig>  as undesirable states, so that the set of undesirables states is u = { <dig> ..., 2047}.
before intervention the probability mass of undesirable states πu is  <dig> .
the optimal intervention for this network is transitioning state  <dig> to state 01011001101; that is, f˜= <dig> for the regulatory function after intervention.
the undesirable probability mass after intervention π˜u is  <dig> .
to evaluate the proposed approximate method, we randomly selected  <dig> different sets of k =  <dig> regulations and assumed that they are uncertain, meaning that their regulatory information is unknown.
if experiments are chosen according to the optimal experimental design method in  <cit> , then ρ¯= <dig> .
table  <dig> compares the average gain ρ¯ of the experiments chosen by the approximate method when deleting p =  <dig> genes and excluding s genes from the search space using the cod-based algorithm.
the table shows that we can obtain meaningful gain when the approximate experimental design method is used to select the experiment to be conducted first.
CONCLUSIONS
we have proposed a computationally effective experimental design method for reducing uncertainty in gene regulatory networks.
this method can effectively approximate the optimal experimental design method in  <cit> , which is based on the mean objective cost of uncertainty .
to reduce computational complexity, we use network reduction to estimate the optimal and robust interventions needed for finding an optimal experiment.
we introduced a novel cost function for gene deletion that takes into account the effect of gene deletion on the ranking of potential experiments.
because potential experiments are ranked based on the mocu in  <cit> , the proposed cost function is also based on the mocu.
simulation results on both synthetic and real networks show that while our proposed method can greatly reduce computations, its performance is comparable to the optimal method in  <cit>  and much better than random gene deletion.
greater computational reduction is achieved by excluding genes from the search space based on their cod with the target gene whose expression the intervention is aimed at altering.
we have assumed a uniform distribution over the uncertainty class.
if one has relevant prior knowledge, perhaps it can be used to construct a distribution reflecting it.
care must be taken because concentrating the mass of the distribution in the wrong place can lead to poorer results.
in bayesian terminology, the distribution on the uncertainty class is called a prior distribution.
putting a non-uniform prior on Θ does not change the reduction procedure introduced in the current paper; however, some calculations are altered by including the weights.
prior construction is a difficult problem and has been considered in the context of gene regulation, but not in the context of network construction.
rather, pathway knowledge has been used to construct prior distributions governing uncertainty classes of feature-label distributions for optimal bayesian classification  <cit> .
prior construction is particular to each application, examples being gene/protein signaling pathways in discrete phenotype classification  <cit>  and model-based rna-seq classification  <cit> .
prior construction for uncertainty classes of the kind considered in this paper constitutes an important issue for future study - and not just in relation to the specific problem considered herein.
competing interests
the authors declare that they have no competing interests.
authors' contributions
conceived the method: rd, bjy, erd.
developed the algorithm and performed the simulations: rd.
analyzed the results and wrote the paper: rd, bjy, erd.
supplementary material
additional file 1
supplementary materials.
this is a file in pdf format that contains an illustrative step by step example providing details about the proposed experimental design method.
this file also includes box plots for the simulations on 8-gene networks.
click here for file
acknowledgements
this work was supported in part by the national science foundation through the nsf award ccf- <dig>  it was also funded in part by a grant from the los alamos national laboratory.
publication costs for this article were funded by the corresponding author's institution.
this article has been published as part of bmc bioinformatics volume  <dig> supplement  <dig>  2015: proceedings of the 12th annual mcbios conference.
the full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/16/s <dig>

 a key area in the field of translational genomics is to derive therapeutic intervention that can beneficially alter cell dynamics in such a way as to avoid cancerous phenotypes. the first step to derive therapeutic interventions is to understand the regulatory relationships among genes. the interactions among genes are studied in the context of gene regulatory networks. grn models often possess high uncertainty. this inherent uncertainty might be due to many factors such as the complex nature of biological phenomena, lack of enough training data, etc. uncertainty in grns can affect the accuracy and performance of the therapeutic interventions. therefore, biologists aim at reducing uncertainty of grn model via conducting additional biological experiments. however, there are many limitations for carrying out experiments. biological experiments are usually expensive and time consuming as they need to be done on living organisms. on the other hand, the resources are limited which prevents researchers from conducting all experiments they need for identification of grn. thus, it is prudent to prioritize potential experiments based on the information they provide and then conduct the most informative. this process is called experimental design. boolean networks have been extensively used for the study of grns. they have been shown to be effective in capturing the multivariate relationships among entities within a cell. bns also facilitate the use of the well-studied markov decision theory for deriving beneficial interventions. to date, intervention methods for bns have been categorized in two different groups: structural interventions, and dynamical interventions. while structural interventions permanently change the behavior of a network via one-time change in the underlying regulatory structure, the goal in dynamical interventions is to interfere with the signaling of grn through flipping or not flipping the expression state of genes over time. the assumption behind most of these intervention methods is that the network being studied is perfectly known. thus, presence of uncertainty can degrade the performance of the intervention methods. from a translational perspective it is crucial to reduce uncertainty pertinent to the objective, such as intervention. hence, uncertainty should be quantified in such a way that the objective for modeling grn is taken into account. mean objective cost of uncertainty quantifies the uncertainty of model in terms of the expected increased cost due to the presence of uncertainty. an experimental design method based on mocu has been proposed in. the experimental design method in evaluates the effect of each potential experiment in reducing the model uncertainty, which is measured in terms of mocu, and suggests that the experiment which results in the minimum expected remaining mocu should be conducted first. the long-run performance of the experimental design method in is guaranteed to be optimal in terms of reducing the error of interventions obtained after conducting the chosen experiments. although, the method in is optimal, it is computationally expensive. because our final objective is to improve the performance of the therapeutic interventions, method in involves finding optimal interventions for all networks which are compatible with the prior knowledge. finding optimal interventions is computationally expensive whose complexity grows exponentially with the number of genes in network. therefore, the computational complexity of finding optimal experiment can be prohibitively high for large networks. thus, it is inevitable to construct a smaller network via deleting some genes from the original large size network and then estimate the optimal interventions using the resulting reduced network. generally the goal in network reduction methods is to produce networks of smaller size while the dynamical behavior of the original network is preserved. there have been some efforts for network reduction to reduce the complexity of designing interventions. in this paper, we propose a novel cost function for the gene deletion process which takes into account the disruption in the order of potential experiments when they are ranked according to the experimental design method in. since experiments are ranked based upon the expected remaining mocu or the mocu that is expected to remain after performing the experiment, we desire that the network reduction step has a low effect on the expected remaining mocu corresponding to the potential experiments. when the gene suggested by the cost function are deleted from network, the optimal interventions are estimated using the reduced networks and then they are used for calculating expected remaining mocu for prioritizing potential experiments. we show the effectiveness of our proposed cost effective experimental design method through simulations on synthetic and real networks. the simulation results verify that our method can perform comparable to the optimal experimental design method in with much lower computations. mocu-based optimal experimental design is very general and does not even require a markovian network. as we will see, finding the best gene to delete is also very general; however, once the genes are deleted, the regulatory structure of the original network must be mapped onto a corresponding regulatory structure on the reduced network, an optimal intervention must be found on the reduced network, and that intervention must be induced to the full network. reduction and inducement are nontrivial and depend on the nature of the regulatory structure. the problem has been addressed for boolean networks in, to which we refer, and a theoretical analysis is given in, where it is noted that the methodology applies to probabilistic boolean networks by applying the reduction to each of constituent bn of the pbn. moreover, whereas we will restrict intervention to rank-one perturbations, which provide a one-time alteration of the regulatory logic, the reduction-inducement paradigm applies to other forms of intervention. methods boolean networks gene regulatory network models are increasingly used as a tool to study interactions among genes. boolean networks and probabilistic boolean networks are widely used models for grns that have been shown to be effective in capturing these interactions. a boolean network on n genes is defined by a pair. v = {x x..., xn} is a set of binary variables that represent the expression state of genes with xi = and xi = corresponding to gene i being off and on, respectively. the on state means that a gene is translated into a protein whereas the off state means that the gene is not translated. it has been shown that significant biological information can be extracted from binarized gene expression data. the gene values at each time step are updated according to the list of boolean functions f = {f f..., fn} with xi = fi where xi xi..., xiki are ki predictor genes for xi. the state of the bn at time t is the vector x =, x..., xn), which is called gene activity profile. a bn with n genes possesses 2n different states. the state of a bn for the next time instant t + is updated as x = f). in a boolean network with perturbation, there is a perturbation probability of each gene flipping its value independently from other genes. thus, in the presence of perturbation x = f) with probability n and with probability − n at least one gene flips its value. the sequence of states of a bnp over time can be viewed as the states of a markov chain with transition probability matrix p = piji,j=12n where pij is the probability that state i transitions into state j. positive perturbation probability p > makes the markov chain ergodic and irreducible, thereby possessing a unique steady-state distribution πt = πt p, where πi is the steady-state probability of state i and t is the transpose operator. more details on how to compute tpm can be found in. the long-run behavior of a bnp is characterized by its ssd, which is indicative of phenotypes. from a translational perspective, the states of a bnp can be partitioned into the set of desirable states d containing those states associated with healthy phenotypes and undesirable states u containing states associated with cancerous phenotypes. the goal of therapeutic interventions is to drive the network away from undesirable states and consequently reduce the steady-state probability mass of undesirable states, πu = ∑uπi. without loss of generality regarding the experimental design analysis, in this paper we will assume that the undesirable states are determined by a single gene, called the "target gene" with the aim of intervention being to flip an undesirable value of the target gene to a desirable value. two intervention approaches are commonly considered: structural interventions change the long-run behavior of the network by altering its underlying rule-based structure, and dynamical interventions affect the dynamical evolution of the network by manipulating the expression states of some "control genes". this paper considers optimal structural intervention via rank-one function perturbations. in a rankone function perturbation intervention, the relation between the tpm before intervention, p, and after intervention, p˜, is p˜=p+abt where a and b are two arbitrary vectors satisfying bt e = where e is an all unity column vector. single-gene perturbation, a special case of rank-one function perturbation, changes the regulatory function f for only one state. if originally f = w, then the single-gene perturbation means that for the regulatory function f˜after intervention, f˜=v and f˜=f for i ≠ u. for this intervention, all entries of p and p˜ would be identical except for two entries. p˜uw=puw-n and p˜uv=puv+n. the steady-state probability of state i after intervention is found using the following equation: π˜i=πi+nπu1-n, where π˜i is the steady-state probability mass of state i after single-gene perturbation and zvi, zwi, zvu, and zwu belong to the fundamental matrix z = − where i is the n × n identity matrix. the undesirable steady state probability mass after single-gene function perturbation is π˜u=∑i∈uπ˜i. to find the optimal intervention, one needs to search among all possible 2n × 2n state pairs. optimal experimental design this section briefly reviews the optimal experimental design method in. suppose θ θ..., θk are k uncertain parameters in the network model. they compose an uncertainty vector θ =, the set of all such vectors being denoted Θ. the collection of all networks corresponding to θ ∈ Θ constitute an uncertainty class, denoted by Θ. let Ψ be a class of potential interventions. let ξθ denote the error of applying intervention ψ ∈ Ψ to the network with uncertainty vector θ ∈ Θ. in case of intervention, ξθ can be the steady-state probability mass of undesirable states after intervention, denoted π˜u,θ. for a network with uncertainty vector θ, the optimal intervention ψ is defined by ψ=argminψ∈Ψξθ. in the presence of uncertainty, we desire an intervention that performs optimally on average across the uncertainty class. the intrinsic bayesian robust intervention ψibr is defined by ψibr=argminψ∈Ψeθξθ, where eθ denotes an expectation relative to the probability density function f governing the uncertainty class Θ. note that although ψibr has the best expected performance across the uncertainty class, it is not guaranteed to have the best performance for each network inside the uncertainty class. the mean objective cost of uncertainty relative to Θ and Ψ is defined by mΨ=eθ, the expression inside the expectation being called the objective cost of uncertainty. mocu is the expected difference between the performance of the robust intervention and the optimal intervention for each network inside Θ. mocu is an uncertainty quantification scheme which measures uncertainty by taking into account the objective of operator design, herein, network intervention. mocu can effectively capture the amount of uncertainty and tries to quantify the uncertainty based on a specific objective to gauge how much it is going to affect the objective. mocu has been used to design experiments that can maximally enhance the performance of the interventions. suppose that there are k potential experiments e e..., ek corresponding to k uncertain parameters θ θ..., θk. it is assumed that experiment ei, which might be a complex experiment consisting of several sub-experiments, fully identifies θi. the goal of experimental design is to find out which experiment ei, ≤ i ≤ k should be conducted first, or how to rank potential experiments effectively. although we do not know the actual outcome of an experiment before conducting it, we do know the possible outcomes. suppose the outcome of experiment ei is ϕ. we define the remaining mocu, given θi = ϕ, as mψ=eθϕξθϕ)-ξθϕ)), where θϕ, called the conditional uncertainty vector, has ith parameter equal to  ϕ and other parameters are still random, and Θi,ϕ is the remaining uncertainty class when θi = ϕ, i.e., Θi,ϕ = {θ|θ ∈ Θ, θi = ϕ}. the expectation is taken over the conditional density function f)=f that governs the remaining uncertainty class Θi,ϕ. ψibr is the robust intervention for Θi,ϕ and is found in a similar way as. corresponding to each experiment ei, we define the expected remaining mocu after conducting that experiment as: mΨ=eθimΨ=eθieθθiξθθi)-ξθθi)), where eθi is the expectation relative to the marginal density function f of uncertain parameter θi. then we find the minimum expected remaining mocu: i*=argmini=,...,kmΨ. ei* is the optimal experiment to be conducted first. approximate experimental design method according to and, calculating the expected remaining mocu requires finding the optimal intervention ψ for each θ ∈ Θ and the robust intervention ψibr for each possible remaining uncertainty class. the complexity of finding optimal interventions grows exponentially with network size n. for finding an optimal single-gene structural intervention, we need to search among all possible 2n × 2n state pairs and calculate the new steady-state probability π˜i for each state i in the set of undesirable states u. thus, the complexity is  o. this heavy computational cost motivates us to reduce the size of network in order to reduce the complexity of finding optimal interventions, thereby reducing the complexity of the experimental design. assuming that gene g is deleted from a network with regulatory function f, we define a regulatory function fred for the reduced network. doing this for each network with uncertainty vector θ in Θ produces the uncertainty class, Θg, of reduced networks via the mapping θ → θg. to approximate the optimal intervention for a network in Θ, we use the corresponding network in Θg, find the optimal intervention for the reduced network ψ, and then induce the intervention to the original network in Θ. this approximate optimal intervention denoted by ψ is called the induced optimal intervention. also, to find the induced robust intervention, ψibrind, for Θ, first we find the robust intervention, ψibr, for Θg using and then find the induced robust intervention ψibrind from ψibr. as illustrated in figure in the proposed approximate experimental design method, we find the best gene g∗ for deletion via a novel cost function c and then obtain the induced optimal and robust interventions needed for the mocu calculations in the experimental design step by inducing interventions from uncertainty class of reduced networks Θg∗ to the original uncertainty class Θ. we now aim to find a gene whose deletion results in minimum degradation in the experimental design process. keeping in mind that the experimental design is based on the expected remaining mocu for potential experiments, let mΨg be the remaining mocu when uncertain parameter i has value θi and we delete gene g, mΨg=eθθiξθθi)-ξθθi)). we define the cost of deleting gene g by c=∑i=1k|mΨg-mΨ|, where mΨg=eθimΨg. the gene g∗ minimizing the cost function in is selected for deletion: g*=argming∈,...,nc. the intuition behind this cost function is that our choice of optimal experiment is based upon the expected remaining mocu corresponding to each experiment. therefore, we desire that the network reduction step has minimum effect on these quantities. deleting genes increases the inherent uncertainty of the network because the induced robust intervention cannot perform better than the original robust intervention on average. we want to reduce this increase in the uncertainty of model caused by network reduction. since eθθiξθθi)≥eθθiξθθi), mΨg≥mΨ. hence, we can omit the absolute value operator in to obtain g*=argming∑i=1k-mΨ)=argming∑i=1kmΨg, where the second equality follows from the fact that mΨ does not depend on the gene being deleted. expanding mΨg yields g*=argming∑i=1keθieθθiξθθi)-ξθθi))=argming∑i=1keθieθθiξθθi). the minimization problem in is equivalent to the one in. based on the cost function in, for each gene g, we find the expected performance of the induced robust intervention ψibrind across the remaining uncertainty class Θi,θi, then take the expectation of this average performance relative to the marginal distribution of the uncertain parameter θi, and finally sum all values found for the k uncertain parameters. algorithm approximate experimental design 1: input: Θ, Ψ, f, θ = 2: output: ei*, i∗ ∈ {..., k}: the estimated optimal experiment to be conducted first 3: for g =: n do 4:   cost ← 0 5:   for i =: k do 6:      for all θi do 7:         build remaining uncertainty class of reduced networks Θi,θig 8:         compute conditional density function f) 9:         find induced robust intervention Ψibrind 10:         hg←eθθiξθθi) 11:      cost ← cost + eθi 12: g*←argming=,...,ncost 13: for i =: k do 14:   for all θi do 15:      build remaining uncertainty class Θi,θi 16:      compute conditional density function f) 17:      compute mΨg* via equation using Ψibrind and Ψind;g*) 18:      mΨg*←eθimΨg* 19: i*←argmini=,...,kmΨg* 20: return i∗ after removing gene g∗, we find the expected remaining mocu corresponding to each experiment using equation by replacing ψ) with ψind;g*) and ψibr with ψibrind. an abstract form of the proposed experimental design method has been given in algorithm a step by step toy example illustrating algorithm is also provided in the additional file this procedure for estimating optimal experiment via deleting one gene can be easily extended to the deletion of two or more genes. for example, to delete two genes, we need to evaluate the cost function in for all possible two-gene combinations and delete the pair whose cost is minimum. reduction mappings and induced interventions if we want to delete gene g from network, we need to find the regulatory function fred for the reduced network. following, every two states of the original network that differ only in the value of gene g can be collapsed to find the transition rule of the reduced network. let s1g and s0g be two states with value and for gene g, respectively, and identical values for other genes. state sg¯ can be obtained from either s1g or s0g by removing the value of gene g. if for the original network, the transition rules for these two states are f = p and f = q, then for the reduced network, fred=pg¯ if πs1g>πs0g and otherwise fred=qg¯, where pg¯ and qg¯ are found from states p and q via removing the value of gene g, respectively. following this procedure, we find the regulatory function fred for all states in the reduced network. algorithm finding induced optimal interventions 1: input: ψ=,π 2: output: ψind= 3: u˜1g← place in the gth coordinate of u˜ 4: u˜0g← place in the gth coordinate of u˜ 5: if πu˜1g≥πu˜0g then 6:   ugind←u˜1g 7: else 8:   ugind←u˜0g 9: υ˜1g← place in the gth coordinate of υ^ 10: υ˜0g← place in the gth coordinate of υ^ 11: if πv^1g≥πv^0g then 12:   υgind←υ˜1g 13: else 14:   υgind←υ˜0g as illustrated in algorithm we find the induced optimal intervention from the optimal intervention for the reduced network. suppose that the optimal intervention for the reduced network θg is ψ=. the two corresponding states to u^ in the original network are u˜1g and u˜0g, which are found by placing and in the gth coordinate of u^, respectively. similarly, there are two states v˜1g and v˜0g in the original network corresponding to state υ^ the induced optimal intervention for the original network is ψind=, where ugind is the one among u˜1g and u˜0g having larger steady-state probability in the original network and vgind is the one among υ˜1g and υ˜0g with larger steady-state probability in the original network. analogous to the induced optimal intervention, the induced robust intervention ψibrind is found from the robust intervention ψibr according to algorithm 3; however, here we choose the two states possessing larger expected steady-state probability across Θ using the expected ssd, π = eθ, where π is the ssd of the network with uncertainty vector θ in uncertainty class Θ. we can use this procedure to find the induced robust intervention for each remaining uncertainty class Θi,ϕ. algorithm finding induced robust interventions 1: input: ψibr=,π∀θ∈Θ 2: output: ψibrind= 3: π←eθ 4: u˜1g← placing in the gth coordinate of u^ 5: u˜0g← placing in the gth coordinate of u^ 6: if πu˜1g≥πu˜0g then 7:   ugind←u˜1g 8: else 9:   ugind←u˜0g 10: v˜1g← placing in the gth coordinate of υ^ 11: v˜0g← placing in the gth coordinate of υ^ 12: if πυ˜1g≥πυ˜0g then 13:   υgind←υ˜1g 14: else 15:   υgind←υ˜0g preliminary gene elimination via the coefficient of determination to further reduce the computational cost of the experimental design, we utilize the coefficient of determination to eliminate some genes from the optimization problem without evaluating the cost function and then search among the remaining genes for choosing genes to be removed using the cost function. the cod measures the strength of relationship between a target gene y and a vector x of predictor genes as the difference between the error of the best estimation of gene y in the absence of other genes and in the presence of genes in x. the cod is between and and a larger cod means a stronger connection between the target and predictor genes, in our case the target gene being the aim of intervention. we use the intuition that genes possessing large cod in relation to the target gene are not likely among the genes that should be deleted because they have strong connection to the target gene. the cod of the target gene y, relative to a vector x = of predictor genes is defined by codx=εy-εx,yεy where εy is the error of the best estimation of y without any predictors, εy=min, and εx,y is the error of the optimal estimation of y upon observing x. by assuming that the value of the binary vector x of predictor genes changes from to 2m, εx,y can be calculated by εx,y=∑j=12mpmin. if codx denotes the cod of y relative to x in a network with uncertainty vector θ, then given the uncertainty class Θ the expected cod of y relative to x is given by codx=eθ. genes possessing strong connection with the target gene in terms of codx are not considered for deletion. when excluding genes using the cod it is important to recognize the possibility of intrinsic multivariate prediction, where a set of genes may have low individual cods with respect to the target gene but may have significant cod when used together for multivariate prediction. first we calculate codx for all 3-gene combinations and pick the one with largest cod. we compute cod for 3-gene predictors because it has been shown in that the average connectivity of the model cannot be too high providing that the model is not chaotic and it is commonplace to assume 3-gene predictivity in bns. if we want to exclude less than genes from the search space, then among the 3-gene combination with the largest expected cod, we choose those genes that have larger expected individual cod. if we want to exclude more than genes, then in addition to the three genes in the combination with the largest cod, we choose those genes in the 3-gene combination with the second largest cod that have larger expected individual cod and do not belong to the first 3-gene combination. we repeat this process until we reach the desired number of genes to exclude. if there are initially n genes and we want to delete genes, then we need to evaluate cost function for all c 3-gene combinations, where c denotes the number of combinations of n objects taken k at a time; however, if we exclude s genes from search space then the number of evaluations of decreases to c. having performed the cod-based exclusion process and excluded s genes, x1′,x2′,...,xs′, we search for the genes to be deleted using the cost function in among the remaining genes, {x x..., xn} − {x1′,x2′,...,xs′}. computational complexity analysis the first step for the optimal experimental design in is estimating optimal interventions ψ for each network in Θ. we also need to compute the robust intervention ψibr for each possible remaining uncertainty class Θi,θi. most of the computations are devoted to this step. finding robust interventions does not require additional calculations because we can store the error of each intervention ψ∈Ψ for the network θ when finding optimal interventions and later use these errors to find robust interventions. therefore, complexity analysis requires computing the complexity of estimating the optimal interventions. with n genes that take on binary expression levels, the network has 2n states. finding an optimal single-gene function intervention requires searching among all possible 22n state pairs according to. assuming without loss of generality that states 2n− to 2n are undesirable, must be evaluated 2n− times for each state pair. thus, the complexity of finding the optimal intervention ψ is  o. if there are k uncertain parameters and each can take on l different values, then the uncertainty class Θ contains lk different networks for which an optimal intervention must be found. hence, the complexity of the optimal experimental design method in is  o. to analyze the complexity of the proposed approximate method, suppose p genes are to be deleted. then the cost function in must be evaluated for all c p-gene combinations, n − instead of n because the target gene cannot be deleted. the complexity of finding an induced optimal intervention for each network after deleting p genes is  o). therefore, the complexity of the approximate method is  o × lk × 23n−3p). for large n, it is possible that for small p the complexity of the approximate method can exceed that of the original method; however, by deleting more genes the complexity of the approximate method drops sharply because by deleting each additional gene the complexity of estimating the optimal intervention decreases by eight-fold. by incorporating the cod-based gene exclusion step in the approximation method and excluding s genes we are able to decrease the number of p-gene combinations from c to c, which reduces the complexity of the approximate method to  o × lk × 23n−3p). define the computational gain λ by λ=lk×23nc×lk×23n-3p=23pc, which is the ratio of the complexity of the optimal method in to the complexity of the approximate method when deleting p genes using the cost function in and excluding s genes from the search space using the cod-based gene exclusion step.  this section evaluates the performance of the approximate method for both synthetic and real grns where the majority vote rule is used as the transition rule. majority vote rule is popular in systems biology, especially when we are interested in the overall dynamics of the network. for example, majority vote is used in to model the dynamics of yeast cell-cycle network. for the majority vote rule, a regulatory matrix r is defined component-wise by rij=1genejactivatesgenei-1genejsuppressesgenei0norelationfromgenejtogenei. according to this rule, gene i takes value if the number of genes that are on and activate it is more than the number of genes that are on and suppress it: xi=fi)=1if∑jrijxj>00if∑jrijxj<0xiif∑jrijxj= uncertainty is introduced by assuming that the exact values of some of the nonzero components of r are unknown; that is, for some regulations it is not known whether they are activating or suppressive. each uncertain parameter θi can be − or conducting experiment ei determines the value of parameter θi. let μ = denote the true value for the uncertainty vector θ =. conducting experiment ei results in a remaining uncertainty class Θi,μi consisting of networks with θi = μi and other uncertain parameters being − or for Θi,μi we can determine a robust intervention ψibr. we evaluate the effectiveness of experiment ei in terms of the error of the resulting robust intervention obtained after experiment on the underlying true network, ξμ). we define the gain of conducting the chosen experiment ei* over a random experiment ernd by ρ=ξμ)-ξμ). if ρ > then the chosen experiment outperforms the random experiment; if ρ < then the random experiment outperforms the chosen experiment; and if ρ = then they perform the same. simulation results based on the synthetic bnps for the performance evaluation based on synthetic bnps, we generated networks randomly and chose different sets of k regulations in each to be unknown - simulations in total. we assigned random predictor genes to each gene where each one can be randomly activating or suppressive. the gene perturbation probability was set to. without loss of generality, we assume that states with up-regulated x are undesirable. we removed the regulatory type of those regulations that have been assumed to be uncertain and retained other regulatory information of the network. we assume that all uncertain parameters are independent from each other and have uniform marginal distribution. the analysis can be easily extended to other distributions. because x is the target gene, it was excluded from the reduction process. hence, we look for the best p-gene set to be deleted among {x..., xn}. to consider larger networks, we generated random 13-gene networks and in each chose one set of regulations to be unknown. we used the approximate experimental design method and deleted genes. for this size of network it is not possible to perform the optimal experimental design method or compute original optimal and robust interventions to calculate the gain of the chosen experiment over a randomly selected experiment. hence, we use the induced robust intervention obtained by deleting the set of genes having the minimum cost of deletion as in. therefore, here gain ρ is defined as ρ=ξμ;g*)-ξμ;g*) where g* is the set of genes with minimum cost of deletion. for this set of simulations, the average gain ρ¯ is. note that here the average gain might not be very accurate owing to the small number of simulations. the approximate run time for each simulation was around seconds. we now evaluate the proportion of times that we obtain the optimal experiment found by when using the approximate method. figure shows the percentage of finding the optimal experiment when using the approximate method and deleting different number of genes from 7-gene networks. in this figure, there are and values for deleting or genes, respectively. we observe that deleting the set of genes corresponding to the minimum of the cost function yields the highest likelihood of obtaining an optimal experiment, which is what we would hope for from an efficient approximate method. according to figure when we delete the gene which attains the minimum cost, % of the simulations yield an optimal experiment, whereas this percentage is % when deleting the gene with the largest value of the cost function. similar behavior is observed when deleting or genes. a salient reason that the largest average gain of the approximate method over random experiments is when we delete genes corresponding to the minimum cost function is that it is more likely to get an optimal experiment. an issue that arises when evaluating experimental design on synthetic networks, as opposed to real biological networks, which typically manifest substantial controllability on account of their need to maintain functionality within changing contexts, a large portion of randomly generated networks might not be controllable and therefore not be responsive to intervention. hence, intervention has negligible effect on their ssds and including them in the analysis masks the effect of optimality. to address this issue, we define controllability Δ as the percentage decrease of undesirable probability mass after applying optimal intervention: Δ=πu-π˜uπu×100%, where πu and π˜u are the undesirable probability masses before and after applying optimal intervention to the network, respectively. a larger Δ means that a network is more controllable. figure considers the effect of controllability on the performance of experimental design when networks have n = genes and k = uncertain regulations. the figure shows the average gain ρ¯ for the optimal method and the proposed approximate method for networks possessing controllability greater than a given threshold. we observe that ρ¯ increases when networks are more controllable, regardless of the number of genes deleted from network. note that as controllability increases, the difference between the performance of different methods increases. for example, for all networks the average gain for the optimal method and the proposed method when deleting one, two, three, and four genes is,,,, and, respectively; but for networks with Δ ≥ 40% the average gains are,,,, and, respectively. to evaluate the effectiveness of the cod-based gene exclusion algorithm, we compare the average gain of the approximate method when excluding genes from the search space using the cod-based exclusion algorithm against the average gain when excluding randomly selected genes from the search space. figure shows the average gain ρ¯ for networks with n = genes and k = uncertain regulations. for deleting p genes, we exclude up to 6− p − genes from the search space so that for the largest number of genes excluded, the search space contains at least p + genes. for example, when deleting p = gene, we exclude and genes; when deleting p = genes we exclude and genes; and so on. for each number of uncertain regulations, we observe that the average gain when excluding genes using the cod-based algorithm is always larger than random gene exclusion, regardless of the number of deleted genes. for example, when k = for deleting one gene and excluding and genes randomly, ρ¯=, , , and respectively but using the cod-based scheme and excluding the same number of genes, ρ¯=, , , and respectively. if we delete three genes ρ¯= without gene exclusion, and if we exclude and genes, then ρ¯= and ρ¯=, respectively, when we exclude genes randomly and ρ¯= and ρ¯=, respectively, when we exclude genes based on cod. note that when deleting more genes, the difference between random exclusion and cod-based exclusion increases because as more genes are deleted, exclusion has a larger impact on the number of candidate sets for evaluating the cost function. for example, when deleting gene, if we exclude one gene, then the number of candidate sets decreases from to but when deleting genes, if we exclude one gene, then the number of candidate sets decreases from c = to c = in figure we also show the box plot for the gain of conducting the chosen experiment if we delete genes and exclude genes from the search space either randomly or via the proposed cod-based method for gene networks possessing uncertain regulations. we observe that the first quartile, median, and third quartiles are higher when excluding genes using cod. for example, when randomly excluding genes from the search space, the first quartile, median, and third quartile are − × 10− × 10− and, respectively; however, when excluding genes using cod they are − × 10− × 10− and, respectively. we have also evaluated the performance of the approximate experimental design method when a sequence of experiments is conducted. suppose there are k = uncertain regulations and we conduct five experiments to identify all unknown regulations. for each set of unknown regulations, at each step we utilize the experimental design to choose one of the possible experiments, conduct the chosen experiment, and measure the performance of the robust intervention obtained after the experiment on the underlying true network. continuing, for the remaining uncertain regulations we use the experimental design method and repeat the previous procedure until there is no more unknown regulation remaining in the network. we also do sequential experiments randomly where at each step we choose an experiment randomly, measure the performance of its corresponding robust intervention, and again choose one experiment randomly among the remaining ones. the difference between the undesirable probability mass after applying the robust interventions corresponding to the randomly chosen experiment and the chosen experiment through experimental design at each step is the gain of conducting the chosen experiment at that step. figure shows the average gain over random selection for the optimal method and the approximate design method deleting up to three genes for 7-gene networks and up to four genes for 8-gene networks. the figure indicates that the approximate design method has reliable performance compared to the optimal method. moreover, similar to the optimal method, the average gain increases sharply in the beginning for the approximate method. this is very important in real applications owing to the cost and time required for conducting experiments. note that when we conduct all five experiments the average gain is zero because after five experiments the network is fully known and we can exactly calculate the optimal intervention regardless of the approach taken to choose experiments. performance evaluation based on the colon cancer pathways in this section, we analyze the performance of the proposed experimental design method on the colon cancer pathways used in. we focus on the pathways formed by genes extracted from the complete pathway set, as used in. these are shown in figure 12: stat ras, il hgf, pik3ca, egf, tsc1/tsc mtor, spry pkc, and mek 1/ normal and blunt arrows represent activating and suppressive regulations, respectively. we modeled the pathways as a bnp with perturbation probability. genes are named as they have been introduced. for example, stat is x and mek 1/ is x egf, hgf, and il are three stimulation factors that carry the external signals generated by neighboring cells to downstream genes and activate downstream pathways. signal transducers and activators of transcription constitute a family of transcription factors that can be activated via extracellular signaling proteins such as cytokins and growth factors. these play a major role in regulating downstream processes such as cell growth, survival, and apoptosis. stat is an oncogene observed to be highly activated in many cancers, in particular, colon cancer. hence, stat has been recognized as a legitimate target for cancer therapy. we considered states with up-regulated stat as undesirable states, so that the set of undesirables states is u = {..., 2047}. before intervention the probability mass of undesirable states πu is. the optimal intervention for this network is transitioning state to state 01011001101; that is, f˜= for the regulatory function after intervention. the undesirable probability mass after intervention π˜u is. to evaluate the proposed approximate method, we randomly selected different sets of k = regulations and assumed that they are uncertain, meaning that their regulatory information is unknown. if experiments are chosen according to the optimal experimental design method in, then ρ¯=. table compares the average gain ρ¯ of the experiments chosen by the approximate method when deleting p = genes and excluding s genes from the search space using the cod-based algorithm. the table shows that we can obtain meaningful gain when the approximate experimental design method is used to select the experiment to be conducted first.  we have proposed a computationally effective experimental design method for reducing uncertainty in gene regulatory networks. this method can effectively approximate the optimal experimental design method in, which is based on the mean objective cost of uncertainty. to reduce computational complexity, we use network reduction to estimate the optimal and robust interventions needed for finding an optimal experiment. we introduced a novel cost function for gene deletion that takes into account the effect of gene deletion on the ranking of potential experiments. because potential experiments are ranked based on the mocu in, the proposed cost function is also based on the mocu. simulation results on both synthetic and real networks show that while our proposed method can greatly reduce computations, its performance is comparable to the optimal method in and much better than random gene deletion. greater computational reduction is achieved by excluding genes from the search space based on their cod with the target gene whose expression the intervention is aimed at altering. we have assumed a uniform distribution over the uncertainty class. if one has relevant prior knowledge, perhaps it can be used to construct a distribution reflecting it. care must be taken because concentrating the mass of the distribution in the wrong place can lead to poorer results. in bayesian terminology, the distribution on the uncertainty class is called a prior distribution. putting a non-uniform prior on Θ does not change the reduction procedure introduced in the current paper; however, some calculations are altered by including the weights. prior construction is a difficult problem and has been considered in the context of gene regulation, but not in the context of network construction. rather, pathway knowledge has been used to construct prior distributions governing uncertainty classes of feature-label distributions for optimal bayesian classification. prior construction is particular to each application, examples being gene/protein signaling pathways in discrete phenotype classification and model-based rna-seq classification. prior construction for uncertainty classes of the kind considered in this paper constitutes an important issue for future study - and not just in relation to the specific problem considered herein. competing interests the authors declare that they have no competing interests. authors' contributions conceived the method: rd, bjy, erd. developed the algorithm and performed the simulations: rd. analyzed the results and wrote the paper: rd, bjy, erd. supplementary material additional file 1 supplementary materials. this is a file in pdf format that contains an illustrative step by step example providing details about the proposed experimental design method. this file also includes box plots for the simulations on 8-gene networks. click here for file acknowledgements this work was supported in part by the national science foundation through the nsf award ccf- it was also funded in part by a grant from the los alamos national laboratory. publication costs for this article were funded by the corresponding author's institution. this article has been published as part of bmc bioinformatics volume supplement 2015: proceedings of the 12th annual mcbios conference. the full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/16/s 
this method incorporates a network reduction scheme by introducing a novel cost function that takes into account the disruption in the ranking of potential experiments. experimental designgene regulatory networksmean objective cost of uncertaintyobjective-based network reductionboolean networksstructural intervention13- march 12th annual mcbios conference little rock, ar, usa in the process, one must find the optimal intervention for every gene regulatory network compatible with the prior knowledge, which can be prohibitively expensive when the size of the network is large. the experimental design method utilizes the concept of mean objective cost of uncertainty. simulation results based on synthetic and real gene regulatory networks show that the proposed approximate method has close performance to that of the optimal method but at lower computational cost. the proposed approximate method also outperforms the random selection policy significantly. the authors have already proposed an optimal experimental design method based upon the objective for modeling gene regulatory networks, such as deriving therapeutic interventions. in this paper, we propose a computationally efficient experimental design method. a matlab software implementing the proposed experimental design method is available at http://gsp.tamu.edu/publications/supplementary/roozbeh15a/. an accurate understanding of interactions among genes plays a major role in developing therapeutic intervention methods. we then estimate the approximate expected remaining mocu at a lower computational cost using the reduced networks. the optimal experiment to be conducted first is the one which leads to the minimum expected remaining mocu subsequent to the experiment. the process of prioritizing biological experiments to reduce the uncertainty of gene regulatory networks is called experimental design. gene regulatory networks often contain a significant amount of uncertainty.
 an accurate understanding of interactions among genes plays a major role in developing therapeutic intervention methods. gene regulatory networks often contain a significant amount of uncertainty. the process of prioritizing biological experiments to reduce the uncertainty of gene regulatory networks is called experimental design. under such a strategy, the experiments with high priority are suggested to be conducted first.  the authors have already proposed an optimal experimental design method based upon the objective for modeling gene regulatory networks, such as deriving therapeutic interventions. the experimental design method utilizes the concept of mean objective cost of uncertainty. mocu quantifies the expected increase of cost resulting from uncertainty. the optimal experiment to be conducted first is the one which leads to the minimum expected remaining mocu subsequent to the experiment. in the process, one must find the optimal intervention for every gene regulatory network compatible with the prior knowledge, which can be prohibitively expensive when the size of the network is large. in this paper, we propose a computationally efficient experimental design method. this method incorporates a network reduction scheme by introducing a novel cost function that takes into account the disruption in the ranking of potential experiments. we then estimate the approximate expected remaining mocu at a lower computational cost using the reduced networks.  simulation results based on synthetic and real gene regulatory networks show that the proposed approximate method has close performance to that of the optimal method but at lower computational cost. the proposed approximate method also outperforms the random selection policy significantly. a matlab software implementing the proposed experimental design method is available at http://gsp.tamu.edu/publications/supplementary/roozbeh15a/. experimental designgene regulatory networksmean objective cost of uncertaintyobjective-based network reductionboolean networksstructural intervention13- march 12th annual mcbios conference little rock, ar, usa 