Volume 12 Supplement 10
Proceedings of the Eighth Annual MCBIOS Conference. Computational Biology and Bioinformatics for a New Decade
Probabilistic reconstruction of the tumor progression process in gene regulatory networks in the presence of uncertainty
 Mohammad Shahrokh Esfahani^{1}Email author,
 ByungJun Yoon^{1} and
 Edward R Dougherty^{1, 2}
DOI: 10.1186/1471210512S10S9
© Esfahani et al; licensee BioMed Central Ltd. 2011
Published: 18 October 2011
Abstract
Background
Accumulation of gene mutations in cells is known to be responsible for tumor progression, driving it from benign states to malignant states. However, previous studies have shown that the detailed sequence of gene mutations, or the steps in tumor progression, may vary from tumor to tumor, making it difficult to infer the exact path that a given type of tumor may have taken.
Results
In this paper, we propose an effective probabilistic algorithm for reconstructing the tumor progression process based on partial knowledge of the underlying gene regulatory network and the steady state distribution of the gene expression values in a given tumor. We take the BNp (Boolean networks with pertubation) framework to model the gene regulatory networks. We assume that the true network is not exactly known but we are given an uncertainty class of networks that contains the true network. This network uncertainty class arises from our partial knowledge of the true network, typically represented as a set of local pathways that are embedded in the global network. Given the SSD of the cancerous network, we aim to simultaneously identify the true normal (healthy) network and the set of gene mutations that drove the network into the cancerous state. This is achieved by analyzing the effect of gene mutation on the SSD of a gene regulatory network. At each step, the proposed algorithm reduces the uncertainty class by keeping only those networks whose SSDs get close enough to the cancerous SSD as a result of additional gene mutation. These steps are repeated until we can find the best candidate for the true network and the most probable path of tumor progression.
Conclusions
Simulation results based on both synthetic networks and networks constructed from actual pathway knowledge show that the proposed algorithm can identify the normal network and the actual path of tumor progression with high probability. The algorithm is also robust to model mismatch and allows us to control the tradeoff between efficiency and accuracy.
Background
The construction of gene regulatory networks is an extremely difficult problem owing to their complexity and the difficulty of obtaining the relevant time series data, in terms of sampling rate, measurement accuracy, and quantity. For instance, microarray data usually come in samples much too small for accurate inference, have a very low sampling rate relative to most cell signaling, measure average transcript across many cells, and are susceptible to many confounding factors which adversely affect the signaltonoise ratio. In particular, for human cells, with data coming from patients, there are no timecourse data and the data come from cells that have already undergone a sequence of mutations, so that the regulatory mechanisms of the original cell are no longer intact. Rather than depend on expression data, one can use known pathway information to construct regulatory relations and thereby develop an uncertainty class of networks whose regulatory dynamics are consistent with the pathway knowledge. An algorithm for doing this has been developed in the context of Boolean networks [1]. If one could obtain wildtype timecourse data, then one could reduce this uncertainty class by standard Boolean network inference methods. Given that in practice we usually only have access to stationary patient data and that the progression of mutations leading to the cancerous state has already occurred, we would like to use the available data to reduce the uncertainty class. In fact, since all we require is that we have an uncertainty class to begin with and wish to use the tumor data, from an algorithmic perspective it does not matter whether the uncertainty class arises from prior biological knowledge, wildtype data, or a combination of both. The proposed algorithm operates on the basis of probabilistic sequential faultdetection, which views regulatory alterations, such as gene mutations, as faults in the network wiring [2]. It estimates the sequence of faults leading to the current (cancerous) regulatory structure, and from these estimates, a reduced uncertainty class for the original (healthy) network. By taking this approach the algorithm simultaneously accomplishes a dual purpose: network inference and fault detection.
The methodology is based on certain fundamental notions regarding cancer development, in particular, that the formation of a tumor is a complex process usually proceeding over a period of decades. Normal cells evolve into cells with increasingly neoplastic phenotypes. Tumor progression is driven by a sequence of randomly occurring mutations and epigenetic alterations of DNA that affect the genes controlling cell proliferation, survival, and other traits associated with malignant cell phenotype. To wit, tumor progression is a multistep process of changes in the regulatory pathways. A set of pathways must be deregulated during the tumor progression until the tissue reaches a cancer state. A wide variety of normal adult human cell types can be transformed experimentally by perturbing five pathways [3]. Certain normal human cells require a greater or lesser number of changes before they will become transformed. Moreover, the regulatory pathways can be altered in many different ways leading to the same cancer. For instance, studies on colon cancer show that the great majority (~ 90%) of colon carcinomas suffer inactivation of the APC gene on Chromosome 5q as an early step in this process, about 40% to 50% acquire a Kras mutation, 50% to 70% show an LOH of Chromosome 17p involving p53, and about 60% show an LOH on Chromosome 18q. Most colon cancers will therefore begin with a Chromosome 5 alteration, but then will take alternative genetic paths on the road toward fullfledged malignancy [3].
In sum, although some common alterations may happen in tumor progression, different patients confront with different types of alterations during the progression, thereby making it important to find a way to identify mutations in order to apply appropriate intervention. Very little work has been done on the identification of genetic alterations (e.g. mutations) in cancer progression using network models. One such example is the work by Gerstung et al. [4], where they predicted cancer progression by applying a conjunctive Bayesian network, in which the order of gene mutations is extracted.
In this paper, we use the Boolean networks with perturbation (BNp) framework to model signaling pathways and ultimately predict the gene mutations that occurred during the tumor progression process. Boolean networks (BNs) have been used in a variety of other contexts and with different objectives in biological applications. Kauffman [5] proposed that the cell types are the attractors. He introduced randomization into the networks, in terms of environmental noise (random perturbation of individual genes) and mutation (not to be confused with the notion of mutation in cancer progression), which refers to changes in the wiring of the network. Random BNs and their characteristics have been extensively studied by Aldana et. al [6]. In a random BN, the average function indegrees are constant and function outputs are assigned randomly. Serra et al. [7] investigated the effects of perturbation in the context of random BNs by knocking out a single gene. Additionally, intervention in BNp has been also studied by Dougherty et al. [8] and Qian and Dougherty [9].
In this work, we focus on BNs with perturbation owing to their role in modeling gene regulatory networks, a key point being that their dynamics can be modeled as Markov chains, thereby facilitating the modeling of genetic alterations in signaling pathways by shifting the network steady state distribution (SSD) from the normal (healthy) SSD toward the cancerous SSD. Having this tool in one hand to model signaling pathways, and the cancerous SSD extracted from the malignant tissue (e.g., based on gene expression data) on the other hand, one can test all the possible alterations on the BN satisfying the pathway information to see which one makes the SSD of the altered BN as close to the cancerous SSD as possible. This allows one to track the sequence of mutations. However, there are two concerns for using BNps to model signaling pathways: (1) the network perturbation probability should be determined, and (2) signaling pathways provide us with incomplete information, which means that there may be too many BNs that satisfy the pathway information. The first issue can be mitigated by finding a good estimate of the perturbation probability. For example, inferring a BNp from a sequence of gene expression data has been studied in [10]. In fact, the second issue is the main source of uncertainty in our problem. To the best of our knowledge, the paper by Layek et al. [1] is the only work that proposed a method to extract the BN underlying the normal tissue from a set of biological pathways. Although this paper introduced an elegant method for extracting the information needed for constructing Boolean networks from biological pathways, it yields a large number of networks since the available network knowledge is often incomplete and not enough to point out the true network. To address this issue, we define the notion of a family of Boolean networks, which is the set of all BNs that satisfy the constraints that are imposed by a given set of pathways. For instance, for a 6gene signaling pathway, the resulting family can contain 2^{12} networks, all of which satisfy the constraints imposed by the pathway.
As mentioned earlier, the main goal of this paper is twofold: (1) to infer the normal network underlying healthy cells from this family, and at the same time, (2) to find the set of alterations that have occurred during the cancer progression. Toward this goal, we propose a probabilistic sequential fault detection algorithm that can effectively identify the best candidates for the original normal (healthy) network and the accumulated gene mutations.
Methods
Boolean networks (BNs)
Each x_{ i } represents the state (expression) of gene i, where x_{ i } = 1 and x_{ i } = 0 represent gene i being expressed and not expressed, respectively. It is commonplace to refer to x_{ i } as the i th gene. The list F of Boolean functions represents the rules of regulatory interactions between genes. That is, any given gene transforms its inputs (regulatory factors that bind to it) into an output, which is the state or expression of the gene itself. All genes are assumed to update synchronously in accordance with the functions assigned to them and this process is then repeated. At any time t, the state of the network is defined by a state vector x(t) = (x_{1}(t), x_{2}(t), …,x_{ n }(t)), called a gene activity profile (GAP). Given an initial state, a BN will eventually reach a set of states, called an attractor cycle, through which it will cycle endlessly. Each initial state corresponds to a unique attractor cycle and the set of states leading to a specific attractor cycle is known as the basin of attraction (BOA) of the attractor cycle.
where, Q is the TPM of the corresponding deterministic BN and H is a 2^{ n } × 2^{ n } matrix depending only on n and p [11].
We assume there exists a “normal” BN, denoted N_{ normal }, corresponding to a healthy wildtype phenotype, and a family of BNs possessing identical predictor sets as N_{ normal } such that . We refer to this family as the “uncertainty class” relative to N_{ normal }.
Alterations in cancer progression are commonly gene mutations, and the accumulation of gene mutations is usually responsible for cancer. Gene mutation includes both oncogene activation and tumor suppressor gene deactivation, resulting in either continuous activation or deactivation of the corresponding genes. In the context of the BN model, such alteration in gene x_{ i } leads to permanently setting the boolean function to f_{ i } ≡ 1 or f_{ i } ≡ 0. We denote a gene mutation by a pair (i, k), which indicates that gene x_{ i } is stuck at k ∈ {0, 1}. For convenience, we define and . If a Boolean network N is altered by a mutation (i, k), this mutated BN is denoted as [N; {(i, k)}]. For example, for a 4gene Boolean network N, [N; {(1,0), (4,1)}] refers to a mutated version of N, where gene x_{1} is permanently deactivated and gene x_{4} is permanently activated. In this case, in the regulatory truthtable, we will have f_{1} = 0 and f_{4} = 1 for every set of predictors. Gene mutation, also called “1gene function perturbation”, has been studied [9]. It should be noted that, physically, the order of mutations can make a difference in cancer progression, since alterations affect the regulatory structure, thereby affecting subsequent cancer progression. There is, however, no way to take this into account given that we only have steadystate data and no data on transient behavior. From a mathematical perspective, the commutativity in (3) means that a path is a set of alterations rather than a sequence of alterations; however, we employ the latter terminology owing to its commonplace usage.
where , , and ρ(π_{i},π_{j}) is the KullbackLeibler divergence (KLdivergence) between the SSDs π_{ i } and π_{ j }. This distance measure can be extended to BNs by first building the corresponding BNp for each BN using (2) and a given probability of perturbation p and then computing the distance between the resulting BNps. Without any ambiguity, in what follows, we use the same notation for a BN and the induced BNp for notational simplicity.
Gene mutation effects
Effects of gene mutation in a BNp
The effect of rankone perturbations in the TPM of a Markov chain on the SSD has been studied in the context of structural intervention in gene regulatory networks [9], and more generally in the framework of Markov chain perturbation theory [12]. We can utilize these results to analyze the SSD of the altered BNp, whose TPM is given by (5).
for some vectors a_{ j } and b j satisfying , and a positive integer u ≤ 2^{n – 1}. The proof can be found in Additional file 1. Based on (6) and (7), the SSD of the altered BNp can be iteratively calculated in at most 2^{ n }^{–1} iterations.
Example: effect of mutations in a 3gene network
where π and are the SSD vectors for P and , respectively, which satisfy π^{ T }P = π^{ T } and . Z_{ j } are the corresponding fundamental matrices, as defined earlier.
Overview of the proposed algorithm
Initially, for each network , there can be 6 possible altered networks based on a single gene mutation. These altered networks are shown in the first row of Figure 1, in the middle plot. Among these networks, the algorithm only selects the networks whose SSDs are close to SSD_{ cancer }. Suppose we select the altered networks that are within β_{1}distance of cancerous network. The selected networks constitute a new (and reduced) uncertainty class of networks . Next, each network in can be altered into 4 different networks based on an additional single gene mutation. These networks are shown in the second row of Figure 1, in the middle plot. Among these networks, the algorithm selects only those that are within β_{2}distance from the cancerous network, resulting in a further reduced uncertainty class of networks . The family contains the best candidates for the normal network and the cancerous path. For example, in Figure 1 contains two candidates (N^{1} and N^{2}) for the normal network. For N^{1}, the cancerous path {(1,0), (2,1)}, and equivalently, {(2,1), (1,0)}, may lead it to the cancerous network with the given steady state distribution SSD_{ cancer }. Similarly, N^{2} may be another candidate for the normal network, which may get close to the cancerous network also through the path {(1,0), (2,1)}. Note that the actual number of networks in and that in will depend on the parameters β_{1} and β_{2}, respectively.
Details of the proposed algorithm
The detailed procedure of the proposed algorithm is shown in Algorithm 1. At each step, one additional single gene mutation is considered. Therefore, to detect all M alterations that may have led the normal network into the cancerous network, the algorithm needs to go through M sequential steps. In the first step, we consider all possible single mutations of the form (i, k) for every network N in the family , which results in altered networks. Among these networks, we select only those networks whose SSD can get close enough to SSD_{ cancer }, after M – 1 additional gene mutations. Based on this criterion, we select the networkmutation pairs, whose residual values (distance to the cancerous network measured based on SSD) are smaller than a threshold β_{1}. In the second iteration we start with a new family of BNs that contains the networks selected in the previous iteration. Since the gene that was mutated in the first step cannot go through another mutation, every network in can go through one of 2(n – 1) possible single gene mutations. Among these possible altered networks, we select only those networks whose residual values are smaller than a threshold β_{2}. After repeating these steps M times, the final class will contain the best networkpath pairs, where each pair consists of a candidate for the normal network and the cancerous path that may drive the given network into the cancerous state with the given SSD. In each iteration, the threshold β_{ j } can be used as a control parameter for trading between efficiency and accuracy. In general, we will have β_{1} ≥ β_{2} ≥ … ≥ β_{ M }.
Performance metrics
These two metrics can be used to evaluate the accuracy of the proposed algorithm.
which is exponential with respect to the number of mutations M.
residual value computations, where α_{0} = 1. A smaller β_{ i } will lead to a smaller α_{ i }, thereby reducing the overall complexity of the algorithm. However, this will also increase the probability of missing the true network, hence the parameters β_{ j } can be used to control the tradeoff between computational efficiency and the prediction accuracy of the algorithm. As we can see from (18), the computational complexity of the proposed algorithm is polynomial with respect to the number of genes n (for a fixed M), while it is exponential with respect to the number of mutations M (for a fixed n). However, the parameters α_{ j } (j = 0, … , M – 1) allow one to trade between computational efficiency and prediction accuracy. As a result, the proposed algorithm can accurately reconstruct the cancer progression path in a much more efficient manner compared to the exhaustive search, as will be demonstrated in our simulation results.
Cumulative distribution function of the distance between a random BNp and its mutated version
We estimate the CDF of the distance between a network and its altered version based on random BNs. We define a random Boolean network (RBN) as a BN: (1) whose gene predictors are randomly chosen such that every gene has k predictors, and (2) the truth table of every Boolean function f_{ i } follows an independent and identically distributed Bernoulli(p_{b}) distribution, where p_{ b } is typically called the bias of the Boolean function f_{ i }. By allowing random perturbations with probability p in the RBN, we can obtain a random BNp (RBNp). First, we generate large number of RBNps with certain properties. Second, for each RBNp N_{ p }, we randomly introduce m mutations to obtain an altered network , and measure their distance . Based on these observations, we can estimate the CDF, where m is the number of single gene mutations and p is the perturbation probability in the RBNp.
Results and discussion
Estimating the CDF of the distance between networks
To execute the algorithm, we first estimated the CDF of the distance d between an RBNp and its mutated copy. As with ensemble analysis in [15][16][17], we estimate these CDFs based on a large number of randomly generated networks with similar structural properties. The two most important parameters for generating random BNs are their bias probability p_{ b } and connectivity k. As described earlier, p_{ b } is the mean of the Bernoulli distribution used to randomly generate the predictor function for each gene in a BN, and k is the maximum number of input variables for each of these functions. We randomly generated 4,000 BNps with these properties. For each network, we introduced random gene mutations and computed the distance between the original BNp and the altered BNp. We used the MATLAB function KSDENSITY to find the CDF that best fits the observed distance distribution. We repeated the overall experiment for different numbers of genes n, different perturbation probabilities p, and different numbers of mutations m.
As we can see from Figure 2(B), for 6gene networks, increasing the perturbation probability from p = 0.001 to p = 0.01 decreases the distance. This is intuitive, since gene mutation (see (5)) only affects the regulatory part, which plays less important roles as the perturbation probability p increases. As a result, changing the regulatory structure of a BNp will have less significant effects when p is larger. Figure 2(D) shows the results for 8gene networks, which show similar tendencies.
Performance of the algorithm on synthetic network families
We evaluated the performance of the proposed algorithm based on randomly generated families of Boolean networks through Monte Carlo simulations. All random networks in each of these families have identical structural properties (i.e., p_{ b } and k). In each family, one network whose Boolean functions are canalizing functions was selected as the true normal network N_{ normal }. A canalizing Boolean function is a function in which an input with a specific value determines the output of the function regardless of the other inputs. For instance, f(x_{1}, x_{2}) = x_{1}ORx_{2} is a canalizing function, where x_{1} = 1 (and similarly, x_{2} = 1) will make the output f(x_{1},x_{2}) = 1, regardless of the value of the other input variable. We randomly chose a path of length M and altered the normal network according to the given path to obtain the cancerous network. The steady state distribution (SSD) of the cancerous network was computed, to be used as an input for the proposed algorithm. Next, we used the proposed algorithm to find out whether it was able to infer the true normal network from a given family of networks and correctly predict the actual cancer progression path, when provided with the number of mutations M and the cancerous SSD. In our simulations, we used p_{ cancer } = 0.001, pb = 0.3, and k = 2. We considered 6gene networks with M = 2 mutations and 8gene networks with M = 3 mutations. For the case of 6gene networks, we considered families of size and . For the case of 8gene networks, we considered families of size . The algorithm was implemented using MATLAB 7.9.0 (R2009b), and all simulations have been performed on a desktop computer with 2.67GHz Intel Core i7 CPU and 12GB RAM. Each SSD computation took around 9.2 × 10^{–4} sec and 5.7 × 10^{–3} sec for n = 6 and n = 8, respectively.
Performance of the proposed algorithm evaluated on 500 randomly generated network families.
, p_{ cancer } = p = 0.001, β_{2} = 0.1  

p _{1}  P _{ miss } 


 AVG of  AVG of # of paths  AVG of # of SSD calculations 
p_{1} = 0.1  0.81  0.66  0.58  0.57  24.3  3.71  3,421 
p_{1} = 0.3  0.49  0.45  0.41  0.39  45.11  4.17  4,677 
p_{1} = 0.5  0.25  0.24  0.21  0.20  64.27  4.41  5,918 
p_{1} = 0.7  0.09  0.06  0.04  0.04  94.84  4.60  9,208 
Performance of the proposed algorithm evaluated on 200 randomly generated network families.
, p_{ cancer } = p = 0.001, β_{2} = 0.1  

p _{1}  P _{ miss } 


 AVG of  AVG of # of paths  AVG of # of SSD calculations 
p_{1} = 0.1  0.81  0.65  0.57  0.56  68.3  4.45  13,289 
p_{1} = 0.3  0.49  0.46  0.40  0.39  138.17  4.70  17,319 
p_{1} = 0.5  0.25  0.23  0.17  0.17  206.7  4.97  21,846 
p_{1} = 0.7  0.09  0.05  0.03  0.03  293.4  4.98  36,348 
Performance of the proposed algorithm evaluated on 100 randomly generated network families
, p_{ cancer } = p = 0.001, β_{3} = 0.1  

p_{1}, p_{2}  P _{ miss } 


 AVG of  AVG of # of paths  AVG of # of SSD calculations 
p_{1} = 0.1, p_{2} = 0.1  0.95  0.74  0.68  0.68  28.6  8.74  8,158 
p_{1} = 0.3, p_{2} = 0.3  0.66  0.42  0.36  0.34  57.5  10.2  13,999 
p_{1} = 0.5, p_{2} = 0.5  0.34  0.16  0.13  0.13  111.3  13.04  35,039 
p_{1} = 0.7, p_{2} = 0.7  0.11  0.05  0.03  0.03  123.1  11.45  89,788 
Performance of the proposed algorithm in case of model mismatch. Evaluated on 500 randomly generated network families.
, p_{ cancer } = 0.001, p = 0.003, β_{2} = 0.1  

p _{1} 


 AVG of  AVG of # of paths  AVG of # of SSD calculations 
p_{1} = 0.1  0.88  0.77  0.76  8.65  2.3  3177.1 
p_{1} = 0.3  0.49  0.43  0.42  46.8  4.29  4693.9 
p_{1} = 0.5  0.25  0.20  0.19  61.31  4.35  5802.1 
p_{1} = 0.7  0.18  0.15  0.14  77.82  4.41  6870.7 
Performance of the proposed algorithm in case of model mismatch. Evaluated on 200 randomly generated network families.
, p_{ cancer } = 0.001, p = 0.003, β_{2} = 0.1  

p _{1} 


 AVG of  AVG of # of paths  AVG of # of SSD calculations 
p_{1} = 0.1  0.94  0.82  0.80  14.33  2.52  12,454 
p_{1} = 0.3  0.43  0.37  0.36  140.5  4.9  17,850 
p_{1} = 0.5  0.28  0.22  0.21  172.48  4.54  21,175 
p_{1} = 0.7  0.19  0.13  0.13  234.8  4.60  26,060 
Performance of the proposed algorithm in case of model mismatch. Evaluated on 100 randomly generated network families.
, p_{ cancer } = 0.001, p = 0.003, β_{3} = 0.1  

p_{1},p_{2} 


 AVG of  AVG of # of paths  AVG of # of SSD calculations 
p_{1} = 0.1, p_{2} = 0.1  0.95  0.76  0.75  25.75  7.88  5,724 
p_{1} = 0.3, p_{2} = 0.3  0.48  0.44  0.43  46  10.79  12,720 
p_{1} = 0.5, p_{2} = 0.5  0.21  0.17  0.16  97.8  14.34  30,146 
p_{1} = 0.7, p_{2} = 0.7  0.19  0.14  0.14  100.8  9.69  47,755 
Performance on cancerous networks involving the p53 gene
Case1: deactiviation of p53
Performance of the proposed algorithm in the case when p 53 is deactivated.
p  # networks  # networkpath pairs  # paths  result 

0.001  2048  2048  1  S 
0.003  2048  2048  1  S 
0.005  512  512  1  S 
0.007  0  0  0  F 
Performance of the proposed algorithm in the case when p 53 is deactivated.
p  # networks  # networkpath pairs  # paths  result 

0.001  2048  2048  1  S 
0.003  2048  2048  1  S 
0.005  1904  1904  1  S 
0.007  832  832  1  S 
Case2: amplification of Mdm2
Performance of the proposed algorithm in the case when Mdm 2 is amplified.
p  # networks  # networkpath pairs  # paths  result 

0.001  1088  1174  3  S 
0.003  520  540  2  S 
0.005  0  0  0  F 
0.007  0  0  0  F 
Performance of the proposed algorithm in the case when Mdm 2 is amplified.
p  # networks  # networkpath pairs  # paths  result 

0.001  1088  1184  3  S 
0.003  832  894  3  S 
0.005  520  544  2  S 
0.007  0  0  0  F 
Case3: simultaneous deactivation of p53 and amplification of Mdm2
Performance of the proposed algorithm when p 53 is deactivated and Mdm 2 is amplified.
p _{1}  # networks  # networkpath pairs  # paths  result  # SSD calculations 

p_{1} = 0.1  730  1333  4  F  38,684 
p_{1} = 0.3  2458  3810  4  F  69,050 
p_{1} = 0.5  3937  7208  4  S  93,650 
p_{1} = 0.7  4096  9009  4  S  113,612 
Performance of the proposed algorithm when p 53 is deactivated and Mdm 2 is amplified.
p _{1}  # networks  # networkpath pairs  # paths  result  # SSD calculations 

p_{1} = 0.1  1063  1629  4  F  40,016 
p_{1} = 0.3  1661  2582  4  F  48,038 
p_{1} = 0.5  2704  4089  4  F  69,146 
p_{1} = 0.7  4096  8839  4  S  101,426 
Conclusions
We proposed an effective probabilistic algorithm for reconstructing the tumor progression process. Given an uncertainty class of networks, which arises from our partial knowledge of the true gene regulatory network represented by biological pathways, and the steady state distribution of a cancerous network, the proposed algorithm tries to simultaneously infer the true gene regulatory network that underlies healthy cells and to predict the sequence of gene mutations that occurred during the tumor progression process. As demonstrated by our experiments, based on both randomly generated networks and realistic networks constructed from known biological pathways that involve the tumor suppressor gene p 53, our algorithm can effectively cope with the uncertainty present in gene regulatory networks and accurately infer the normal (healthy) network and the actual path of tumor progression with high probability. Furthermore, the proposed algorithm is robust to model mismatch and provides us with effective means for trading prediction accuracy for computational efficiency.
The computational complexity of the algorithm depends on the number of genes in the network, the number of mutations, and the number of networks in the initial uncertainty class, and increasing any of these numbers will increase the computational overhead. Based on the mathematical representation of Boolean networks, increasing the number of genes will exponentially increase the number of possible networks. However, this rapid increase does not necessarily mean that the size of the uncertainty class of networks that we need to deal with will increase at the same rate. For example, many of the mathematically possible networks may not be considered biologically viable, hence may be omitted in practice. Moreover, although the total number of states in a Boolean network with n genes is 2^{ n }, many states may be eliminated via state reduction, and the reduced network may consist of considerably fewer states [20]. In fact, the whole idea of network reduction is relevant to the present problem, just as it is to determining control polices for gene regulatory networks, where computational intractability prohibits the design of control policies without constraining the network size [21]. For example, even when the gene expression levels are restricted to be binary, a network with 15 genes, absent some form of state reduction, cannot be considered, because the size of the resulting transition probability matrix would be 2^{15} × 2^{15}, making any kind of dynamic or control analysis intractable. Another possible way to reduce the complexity of the algorithm is to restrict the possible gene mutations via the use of prior knowledge. For example, we may restrict the possible mutant genes only to a smaller subset of genes that are known to be susceptible to mutation. Furthermore, prior knowledge concerning the expected type of mutation for a susceptible gene (e.g., “amplification” for oncogenes and “deactivation” for tumor suppressor genes) can be taken into account. Although we did not constrain the possible gene mutations nor applied any network reduction technique in this study, such modifications are fairly straightforward and may be used to enhance the overall computational efficiency of the proposed algorithm.
Declarations
Acknowledgements
This work was supported by the W. M. Keck Foundation.
This article has been published as part of BMC Bioinformatics Volume 12 Supplement 10, 2011: Proceedings of the Eighth Annual MCBIOS Conference. Computational Biology and Bioinformatics for a New Decade. The full contents of the supplement are available online at http://www.biomedcentral.com/14712105/12?issue=S10.
Authors’ Affiliations
References
 Layek R, Datta A, Dougherty E: From biological pathways to regulatory networks. Mol. BioSyst 2011, 7: 843–851. 10.1039/c0mb00263aView ArticlePubMedGoogle Scholar
 Layek R, Datta A, Bittner M, Dougherty E: Cancer therapy design based on pathway logic. Bioinformatics 2011, 27(4):548. 10.1093/bioinformatics/btq703View ArticlePubMedGoogle Scholar
 Weinberg R: The biology of cancer. Garland Science New York; 2007.Google Scholar
 Gerstung M, Baudis M, Moch H, Beerenwinkel N: Quantifying cancer progression with conjunctive Bayesian networks. Bioinformatics 2009, 25(21):2809. 10.1093/bioinformatics/btp505PubMed CentralView ArticlePubMedGoogle Scholar
 Kauffman S: The origins of order: Self organization and selection in evolution. Oxford University Press, USA; 1993.Google Scholar
 Aldana M, Coppersmith S, Kadanoff L: Boolean dynamics with random couplings. Perspectives and Problems in Nonlinear Science 2003, 23–89.View ArticleGoogle Scholar
 Serra R, Villani M, Semeria A: Genetic network models and statistical properties of gene expression data in knockout experiments. Journal of theoretical biology 2004, 227: 149–157. 10.1016/j.jtbi.2003.10.018View ArticlePubMedGoogle Scholar
 Dougherty E, Pal R, Qian X, Bittner M, Datta A: Stationary and structural control in gene regulatory networks: basic concepts. International Journal of Systems Science 2010, 41: 5–16. 10.1080/00207720903144560View ArticleGoogle Scholar
 Qian X, Dougherty E: On the longrun sensitivity of probabilistic Boolean networks. Journal of theoretical biology 2009, 257(4):560–577. 10.1016/j.jtbi.2008.12.023PubMed CentralView ArticlePubMedGoogle Scholar
 Marshall S, Yu L, Xiao Y, Dougherty E: Inference of a probabilistic Boolean network from a single observed temporal sequence. EURASIP Journal on Bioinformatics and Systems Biology 2007, 2007: 5.View ArticleGoogle Scholar
 Ivanov I, Simeonov P, Ghaffari N, Qian X, Dougherty E: Selection policyinduced reduction mappings for Boolean networks. Signal Processing, IEEE Transactions 2010, 58(9):4871–4882.View ArticleGoogle Scholar
 Hunter J: Stationary distributions and mean first passage times of perturbed Markov chains. Linear Algebra and its Applications 2005, 410: 217–243.View ArticleGoogle Scholar
 Schweitzer P: Perturbation theory and finite Markov chains. Journal of Applied Probability 1968, 5(2):401–413. 10.2307/3212261View ArticleGoogle Scholar
 Qian X, Dougherty E: Effect of function perturbation on the steadystate distribution of genetic regulatory networks: optimal structural intervention. IEEE Trans. Signal Process 2008, 56(10–1):4966–4976.View ArticleGoogle Scholar
 Kauffman S, Peterson C, Samuelsson B, Troein C: Genetic networks with canalyzing Boolean rules are always stable. Proceedings of the National Academy of Sciences of the United States of America 2004, 101(49):17102. 10.1073/pnas.0407783101PubMed CentralView ArticlePubMedGoogle Scholar
 Shmulevich I, Dougherty E: Genomic Signal Processing (Princeton Series in Applied Mathematics). Princeton University Press; 2007.Google Scholar
 Shmulevich I, Lahdesmaki H, Dougherty E, Astola J, Zhang W: The role of certain Post classes in Boolean network models of genetic networks. Proceedings of the National Academy of Sciences of the United States of America 2003, 100(19):10734. 10.1073/pnas.1534782100PubMed CentralView ArticlePubMedGoogle Scholar
 Batchelor E, Loewer A, Lahav G: The ups and downs of p53: understanding protein dynamics in single cells. Nature Reviews Cancer 2009, 9(5):371–377. 10.1038/nrc2604PubMed CentralView ArticlePubMedGoogle Scholar
 Karnaugh M: The map method for synthesis of combinational logic circuits. Trans. AIEE. pt. I 1953, 72(9):593–599.Google Scholar
 Qian X, Ghaffari N, Ivanov I, Dougherty E: State reduction for network intervention in probabilistic Boolean networks. Bioinformatics 2010, 26(24):3098. 10.1093/bioinformatics/btq575PubMed CentralView ArticlePubMedGoogle Scholar
 Ghaffari N, Ivanov I, Qian X, Dougherty E: A CoDbased reduction algorithm for designing stationary control policies on Boolean networks. Bioinformatics 2010, 26(12):1556. 10.1093/bioinformatics/btq225View ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.