Volume 5 Supplement 2
Selected Proceedings of the 6th International Symposium on Bioinformatics Research and Applications (ISBRA'10)
Constraintbased analysis of gene interactions using restricted boolean networks and timeseries data
 Carlos HA Higa^{1, 2},
 Vitor HP Louzada^{1},
 Tales P Andrade^{1} and
 Ronaldo F Hashimoto^{1}Email author
DOI: 10.1186/175365615S2S5
© Higa et al; licensee BioMed Central Ltd. 2011
Published: 28 April 2011
Abstract
Background
A popular model for gene regulatory networks is the Boolean network model. In this paper, we propose an algorithm to perform an analysis of gene regulatory interactions using the Boolean network model and timeseries data. Actually, the Boolean network is restricted in the sense that only a subset of all possible Boolean functions are considered. We explore some mathematical properties of the restricted Boolean networks in order to avoid the full search approach. The problem is modeled as a Constraint Satisfaction Problem (CSP) and CSP techniques are used to solve it.
Results
We applied the proposed algorithm in two data sets. First, we used an artificial dataset obtained from a model for the budding yeast cell cycle. The second data set is derived from experiments performed using HeLa cells. The results show that some interactions can be fully or, at least, partially determined under the Boolean model considered.
Conclusions
The algorithm proposed can be used as a first step for detection of gene/protein interactions. It is able to infer gene relationships from timeseries data of gene expression, and this inference process can be aided by a priori knowledge available.
Background
One of the goals of Systems Biology is to study the various cellular mechanisms and components. In many cases these mechanisms are complex, where some of the interactions between the proteins are still unknown. To represent these interactions it is common to use gene regulatory networks (GRN). There are several models of GRN, both discrete and continuous. The simplest discrete model was introduced by Kauffman [1] and its known as Boolean network. Later, this model was modified to express uncertainty giving rise to the probabilistic Boolean network[2, 3]. Friedman introduced Bayesian networks[4] as a probabilistic tool for the identification of regulatory data and showed that they can reproduce certain known regulatory relationships. Among the continuous models we can cite the ordinary differential equations which was suggested several decades ago [5]. For a more detailed review about models of gene regulatory networks see [6].
Models of gene regulatory networks help us study biological phenomena (e.g. cell cycle) and diseases (e.g. cancer). Therefore, revealing such networks, or at least some of its connections, is an important problem to address. The ability to uncover the mechanisms of GRN has been possible due to developments in highthroughput technologies, allowing scientists to perform analysis on the DNA and RNA levels. The most common type of data provided by these technologies are gene expression data (microarray). The biological systems are notoriously complex. Determining how the pieces of this puzzle come together to create living systems is a hard challenge known as reverse engineering, which is the process of elucidating the structure of a system by reasoning backwards from observations of its behavior [7]. However, in many cases, GRN cannot be precisely unraveled due to measurement noise and the limited number of data sets compared to the number of genes that are involved.
The most common approach to reverse engineering GRN is to use gene expression data. For example, Marshall et al. [8] considered the approach of using timeseries data and probabilistic Boolean networks (PBNs) as a model of GRN. In fact, there are many other works about inferring PBNs [9–13].Computational algebra approaches were also proposed in [14, 15]. Jarrah et al. [16, 17] used polynomial dynamical systems for reverse engineering of GRN. One good survey for inferring GRN from timeseries data can be found in [18]. Some algorithms use additional information from heterogeneous data sources, e.g. genome sequence and proteinDNA interaction data, to assist the inference process. Hecker et al. [19] presents a good review of GRN inference and data integration.
Usually, an inference algorithm aims to construct one single network which is believed to be the true network. The issue is that the inverse problem is illposed, meaning that several networks could explain (or generate) the data set given as the input for the algorithm. In fact, a study for validation of GRN inference procedures can be found in [20]. The problem becomes more complicated if we take into account the noise that may be present in the data and the small amount of samples. For this reason, our approach aims to analyze several networks that could explain the data. By analyzing the similarities among these networks, we propose a confidence measure of the regulatory relationship between the genes.
In this paper, we present an algorithm for analysis of gene interactions. Although this analysis is directly connected to the process of inference of gene regulatory networks, the main goal of this work is not the inference. The idea is that the algorithm could be used as a first step of an inference process, that is, a preprocessing of the data, in order to support an inference process. To perform the analysis, the algorithm generates a limited number of consistent networks (to be explained in the next section). Unlike any inference algorithm, our algorithm does not take these networks as the final result (the true network). It uses these networks to perform the analysis of gene interactions.
The algorithm is based on Boolean networks and timeseries gene expression. Actually, the Boolean networks are called restricted in the sense that not all Boolean functions are allowed in the model. Restricting the network reduces the search space, which can be significant, since the inverse problem is very complex. This restricted model allows us to find constraints that turn our problem into what can be seen as a Constraint Satisfaction Problem (CSP) and CSP techniques can be used to find feasible solutions, that is, networks. The timeseries data allows us to observe part of the dynamics of the system. These observations are used to generate the constraints of the CSP.
A challenge always presented in any gene regulatory model is its usefulness. It would be interesting if a model could help biological experiments in understanding gene interactions. The model here presented is capable of inferring some of these connections from timeseries data of gene expressions, and this inference process is aided by all a priori knowledge available. What we envisage with our method is a model that points out which connections should be inspected in the wet lab that would constrain as many other connections as possible and consequently could facilitate some biological experiments.
Methods
Restricted Boolean network model
A Boolean network (BN) is defined by a set X = {x_{1}, x_{2},…,x_{ n }} of n Boolean variables and a set F = {f_{1}, f_{2},…,f_{ n }} of n Boolean functions. In the case of GRN the variables are called genes. Each gene x_{ i }, i = 1,…,n, can assume only two possible values: 0 (OFF) or 1 (ON). The value of the gene x_{ i } at time t + 1 is determined by genes at time t through a Boolean function . Given that, there are k_{ i } genes assigned to gene x_{ i }, and the mapping j_{ k } : {1,…,n} → {1,…,n}, k = 1,…,k_{ i } determines the “wiring” of x_{ i }[21]. This way,
We assume that all genes are updated synchronously by the functions in F assigned to them and this process is repeated. The artificial synchrony simplifies computation while preserving the qualitative, generic properties of global network dynamics [22, 23]. A state of the network at time t is a binary vector s(t) = (x_{1}(t),…,x_{ n }(t)). Therefore, the number of states is 2^{ n }, labeled by s_{0}, s_{1},…,s_{2}^{ n }–1. The dynamics of the network is represented by the transition between states. This model is deterministic given that there is a single Boolean function to regulate each gene. Because of the finite number of states and the deterministic behavior, some of the states may be visited cyclically. These states form what is known by the attractor of the BN. The states outside the attractor are called transient states. The transient states together with the corresponding attractor states form the basin of attraction of that attractor.
We call the summation ∑_{ j }a_{ ij }x_{ j }(t) the input of x_{ i } at time t. Besides the regulatory relationships of the matrix A, each gene can have a selfdegradative behavior implying that its value is set to 0 whenever its input is null. Observe that not all Boolean functions can be represented using (2) and that is why the Boolean network is called “restricted”. It is also worth to notice here that each gene x_{ i } depends only on the ith row of A.
Constraint satisfaction problem
A constraint satisfaction problem (CSP) is defined by a set of variables X = {x_{1}, x_{2},…,x_{ n }}; a collection of finite sets D = {D_{1}, D_{2},…,D_{ n }}, where D_{ i } is the domain set for x_{ i }; and a collection of set of constraints C = {C_{1}, C_{2},…,C_{ m }} restricting the values that all variables can assume simultaneously, where each set C_{ i } involves constraints of a subset of variables and specifies the allowable combinations of values for that subset. A solution of a CSP is an assignment, to every variable x_{ i }, a value from its domain D_{ i } such that all constraints in C are satisfied [24, 25].
CSPs defined on finite domains are usually solved by search algorithms meaning systematically assigning possible values to variables and verifying whether all constraints may be satisfied or not. The most used techniques are variants of backtracking, constraint propagation, and local search[24–27]. In these search algorithms, the assignment process requires an order in which the variables are considered. In addition, after selecting a variable, it is necessary to decide the order in which a value is picked up from its domain (to assign to it). In fact, there are many heuristics for variable and value ordering [24, 27]. If one is interesting in generating a sample of uniform solutions, a good heuristic for variable and value ordering may be, firstly, select randomly a variable (with uniform probability), and, then choose randomly (again, with uniform probability) a value from its domain. There are many CSP solvers in the literature; the one we are using is from Gecode project (http://www.gecode.org) [27].
Analysis of gene interactions using CSP
One important remark about the Boolean model presented by Equation (2) is that each row in the regulation matrix A is independent from each other. Using this property, instead of applying the CSP to find all the possible solutions for the A matrix, we can apply the CSP to find all the feasible rows r_{ i } for each gene x_{ i }. In this way, the time complexity of the algorithm is reduced by decreasing the number of possible combinations from 3^{ n }^{2} (number of possible matrices) to n · 3^{ n } (total number of possible rows). Thus, in this context, the gene interaction problem considered in this paper can be modeled as a set of n CSPs. For each gene x_{ i } (i = 1,2, …,n), the problem P_{ i } is defined by the set of variables R_{ i } = {a_{i,1}, a_{i,2},…,a_{ i,n }} (corresponding to the n entries of the ith row of the regulation matrix); a collection of domain sets D_{ i } = {D_{i,1}, D_{i,2},…,D_{ i,n }} (each D_{ i,j } = {–1,0,1}); and a set of constraints obtained by considering all successions of states in the timeseries data (see next subsection).
In order to analyze the gene interactions, since there may be a combinatorial explosion in generating the rows, we consider some of them (a random sampling process) to perform the analysis. In fact, as the number of genes grows, the problem becomes intractable and there could be too many consistent networks for a given data set. We would like to highlight that we are not concerned in generate all the consistent networks, but a limited number of them in order to perform an analysis of the gene interactions.
In situations where a large number of genes is considered and a small number of timeseries data is available, one can reduce the number of genes to perform the gene interaction analysis by employing clustering analysis [28–31] or feature selection algorithms [32], and/or building small subnetworks by using the paradigm of growing seed genes [11, 33]. Of course, one can still use prior knowledge by selecting a small number of genes involved in a specific biological process (e.g., cell cycle division, metabolic pathway).
Constraints generation for the CSP
The algorithm was designed under the assumption that the gene expression data were generated by a biological system which can be modeled as a restricted Boolean network. Let S = {S(1), S(2),…,S(m)} be a set of m timeseries gene expression profiles, where S(t) ∈ {0,1}^{ n } for t = 1,…,m. The algorithm aims to analyze networks that produce the sequence
S(1) → S(2) → ⋯ → S(m) . (3)
When the network produces the timeseries data we say that the network is consistent with the data.
Naturally, there may exist several consistent networks for a single sequence. That is, the inverse problem is a “onetomany” or illposed problem, and this is very difficult to handle.
One naïve way to solve this illposed problem is to find all possible networks by a full search algorithm. In fact, Lau et al. [34] proposed a “smart” full search algorithm to enumerate all possible networks. Here, in this paper, we explore some mathematical properties of the restricted Boolean networks in order to avoid this full search approach.
The algorithm aims to analyze the interactions between the genes through the information provided by the timeseries sequence, which can be seen as a state transition sequence of the corresponding BN. These timeseries data and the restricted Boolean network model are used to generate the constraints of the CSP, as we show next.
First constraints set
Timeseries data
x_{1}(t)  x_{2}(t)  x_{3}(t)  x_{4}(t)  

S(1)  1  1  0  0 
S(2)  1  0  0  0 
S(3)  1  0  1  0 
S(4)  1  0  1  1 
S(5)  0  0  1  1 
Proposition 1. Let S(t – 1), S(t) and S(t + 1) be three consecutive states according to the restricted Boolean network model. If S(t – 1) and S(t) differ by a single gene x_{ k }, then for each gene x_{ i } such that x_{ i }(t) ≠ x_{ i }(t + 1) we have that x_{ k } regulates x_{ i } directly, that is, a_{ ik } ≠ 0.
Proof. Suppose that S(t – 1) and S(t) differ by a single gene x_{ k }, and that there is at least one gene x_{ i } such that x_{ i }(t) ≠ x_{ i }(t + 1). As x_{ i }(t) ≠ x_{ i }(t + 1), the summations ∑_{ j }a_{ ij }x_{ j }(t – 1) and ∑_{ j }a_{ ij }x_{ j }(t) have different signs. Given that x_{ k } is the only gene possessing different values in S(t – 1) and S(t), this difference signal must have been caused by x_{ k }. Therefore, a_{ ik } ≠ 0.
Second constraints set
To generate the second set of constraints, the algorithm takes into account two consecutive states, S(t) and S(t + 1). There is one important observation here: only the active genes at time t can possibly regulate genes at time t + 1. This fact becomes clear when we look at Equation (2). The active genes can give us an insight of which genes are regulating other gene, although the type of the regulatory relationship can not be determined. However, the input given by the summation in Equation (2) can help us to determine the regulatory relationships. For example, if we observe that a gene x_{ i } changes its value from 0 (at time t) to 1 (at time t + 1), we can deduce that the input for gene x_{ i } is positive at time t and only the active genes at time t are responsible for this positive input. Following this logic, the algorithm generates all possible combinations of regulatory relationships using the active genes such that the input of gene x_{ i } at time t is coherent to the values of x_{ i } at time t + 1. More formally, we can state the following proposition.
Possible transitions
x_{ i }(t)  x_{ i }(t + 1)  constraint for a_{ ij } 

1  0 

0  1 

0  0 

1  1 

Proof. Let us first prove the first constraint, that is, if x_{ i }(t) = 1 and x_{ i }(t + 1) =0, then . According to Equation (2), the only way to change the state of gene x_{ i } from 1 (at time t) to 0 (at time t + 1) is if its input . If we just consider the genes x_{ j } that are active at time t, that is, those which x_{ j }(t) = 1, we can rewrite this constraint as . Analogously, one can prove the other constraints given in Table 2.
To exemplify, consider the data in Table 1 where t = 3, that is, S(3). At this time, there are two active genes, x_{1} and x_{3}. These genes are the only ones that can contribute to the sign of the input for each gene for the next time. If we look at gene x_{4}, we observe that its value turned from 0 to 1. According to Equation (2), the input must be positive in this case, that is, . Considering only the active genes at time t = 3, we must have a_{41} + a_{43} > 0. Therefore, neither a_{41} or a_{43} can take the value –1, only 1 or 0 (not both). The same procedure can be applied to all genes and then, the constraints for the CSP are generated.
Third constraints set
This set is generated by analyzing any two pairs of consecutive states in the timeseries data. Let t_{1} and t_{2} be two time points in the timeseries data:
⋯S(t_{1}) → S(t_{1} + 1) → ⋯ → S(t_{2}) → S(t_{2} + 1) ⋯.
Therefore, the difference between S(t_{1} + 1) and S(t_{2} + 1) in this case must be caused by the change on x_{4}. In this step, the algorithm checks how each gene changed in the two pairs of consecutive states.
This result implies that the entry a_{14} of the matrix must have value 1.
If S(t_{1}) and S(t_{2}) differ in more than one gene, we can still generate hypotheses of regulation. In fact, this step tries to construct a system of inequalities with the inputs of each gene for every combination of two consecutive pairs. More formally, we can state the following proposition.
Possible transitions
x_{ i }(t_{ a }) → t_{ a } + 1  x_{ i }(t_{ b }) → x_{ i }(t_{ b } + 1)  Constraints for a_{ ij } 

1→0  0→1 

0→1  0→0 

1→1  0→0 

0→0  0→1 

0→1  1→0 

1→1  1→0 

0→0  1→1 

1→0  1→1 

The other constraints can be obtained in a similar way.
Bar chart of connection frequencies
For a fixed row r_{ i }, the algorithm generates a collection of consistent rows R_{ i } = {r_{i1}, r_{i2},…,r_{ im }} from the constraints generated by the timeseries data and the CSP solver. Each consistent row has n entries , each one corresponding to a one type of connection (inhibition, no connection or activation) on gene x_{ i } from gene x_{ j }. Thus, we can estimate the frequency of all possible connections for the entry a_{ ij } of the regulation matrix by computing the frequency of entries –1, 0 and 1 of all , for k = 1,2,…,m.
We can exhibit the frequency of different types of connections on gene x_{ i } from gene x_{ j }, by showing the estimated frequencies of –1, 0 and 1 for a_{ ij } using a bar chart, as we will see in the next section. Evidently, for a fixed row r_{ i }, determined connections on gene x_{ i } will appear with frequency 100% in all rows R_{ i }; while partially determined connections on gene x_{ i } will have, at least, one type of connection (inhibition, no connection or activation) with frequency 0%; and, for undetermined connections, all relationships will have nonzero frequencies.
Interactions rank
where inh(x_{ i } ← x_{ j }) and act(x_{ i } ← x_{ j }) denote the estimated frequency of –1 (inhibition) and 1 (activation), respectively, for the entry a_{ ij } obtained from the set of the rows in R_{ i }.
For a set of n genes, we rank n^{2}/2 gene interactions. Typically, the interactions with the highest rank are used in order to search for interactions already known in the literature.
Inducing a connection
An interesting application of the bar chart of connection frequencies would be the answer of the following question: which nondetermined (partially determined or undetermined) connection, once determined, would constrain as many other connections as possible? In an experimental context, a method that can point which connection would aggregate more “knowledge” to the network if determined leads to an empirical construction of optimal GRNs.
To investigate this point, a simple exhaustive search was implemented in the space of consistent rows R_{ i } for a gene x_{ i }, according to the following steps:
1. Set a nondetermined (partially determined or undetermined) connection a_{ ij } of the set R_{ i } as a determined connection with a value v of its domain;
2. Set ;
3. Construct the bar chart connection frequencies of ;
4. Compute score(a_{ ij },v) = score(x_{ i }) = ∑_{ j } determination(x_{ i },x_{ j }) on the constrained set , where determination(x_{ i }, x_{ j }) is 0 if the connection x_{ i } ← x_{ j } is undetermined, 0.5 if it is partially determined, and 1 if it is completely determined;
5. Repeat Steps 14 for all nondetermined connections and all their domain values. In the limit, 3^{ g } will be tested, where g is the number of nondetermined connections;
6. Chose a_{ ij } = v that has the highest score.
By the end of these steps, we can say that connection a_{ ij } = v determines as many other connections as possible for gene x_{ i }.
Results and discussion
Budding yeast cell cycle model
Attractors
Basin size  Cln3  MBF  SBF  Cln1  Cdh1  Swi5  Cdc20  Clb5  Sic1  Clb1  Mcm1 

1,764  0  0  0  0  1  0  0  0  1  0  0 
151  0  0  1  1  0  0  0  0  0  0  0 
109  0  1  0  0  1  0  0  0  1  0  0 
9  0  0  0  0  0  0  0  0  1  0  0 
7  0  1  0  0  0  0  0  0  1  0  0 
7  0  0  0  0  0  0  0  0  0  0  0 
1  0  0  0  0  1  0  0  0  0  0  0 
Temporal evolution
Time  Cln3  MBF  SBF  Cln1  Cdh1  Swi5  Cdc20  Clb5  Sic1  Clb1  Mcm1  Phase 

1  1  0  0  0  1  0  0  0  1  0  0  Start 
2  0  1  1  0  1  0  0  0  1  0  0  G_{1} 
3  0  1  1  1  1  0  0  0  1  0  0  G_{1} 
4  0  1  1  1  0  0  0  0  0  0  0  G_{1} 
5  0  1  1  1  0  0  0  1  0  0  0  S 
6  0  1  1  1  0  0  0  1  0  1  1  G_{2} 
7  0  0  0  1  0  0  1  1  0  1  1  M 
8  0  0  0  0  0  1  1  0  0  1  1  M 
9  0  0  0  0  0  1  1  0  1  1  1  M 
10  0  0  0  0  0  1  1  0  1  0  1  M 
11  0  0  0  0  1  1  1  0  1  0  0  M 
12  0  0  0  0  1  1  0  0  1  0  0  G_{1} 
13  0  0  0  0  1  0  0  0  1  0  0  Stationary G_{1} 
Results for the budding yeast artificial data
Rank table
Genes interacting  Rank 

Clb5 and Mcm1  1.72 
Cln3 and SBF  1.70 
Clb5 and Clb1  1.68 
Cln3 and MBF  1.66 
CLN1 and Sic1  1.66 
Sic1 and Clb1  1.66 
Swi5 and Cdc20  1.65 
Clb1 and Mcm1  1.61 
Cdh1 and Clb1  1.55 
CLN1 and Cdh1  1.54 
HeLa cells
HeLa genes
Phase  Genes 

G1/S  CCNE1, E2F1, CDC6 and PCNA 
S  RFC4, DHFR, RRM2 and RAD51 
G2  CDC2, TOP2A, CCNF and CCNA2 
G2/M  STK15, BUB1, CCNB1 and PLK1 
M/G1  PTTG1, RAD21, VEGFC and CDKN3 
Another preprocessing step was to split the data into three data sets, one for each cycle. Considering the sample S = {S(1), S(2),…,S(48)}, we identified the binary state vectors that represent the attractors of the system. For example, the sequence of states S(1),…,S(5) are very similar and we consider them as equal states and represent them as a one singleton attractor (an attractor composed by a single state). The same approach was taken regarding the sequences of states S(16), S(17), S(18); S(30), S(31), S(32) and S(43),…,S(48). Although the binary states in the sequence may be not exactly the same, we are assuming that this difference is caused by the noise in the data. These singleton attractors are similar to the G_{1} stationary state of the budding yeast cell cycle model [35].
We identified three cell cycles C_{1} = S(6),…,S(17); C_{2} = S(19),…,S(31) and C_{3} = S(33),…,S(47). Therefore, to apply the algorithm on this data, we considered the three cell cycles separately. To analyze our results, for a fixed gene x_{ i }, we considered the union of all sets of rows R_{ i }, obtained from the application of the algorithm to each cell cycle, and then compute the connection frequencies.
Results for the HeLa cells
Comparing the three cycles present in the timeseries data, we can see some effects of noise in the gene expression measurements. Supposedly, the three cycles should be equal, but there are minor differences among them. According to [36], the cells utilized in the microarray experiment, by the time of C_{3}, could not be in the same cycle phase, compromising the experiment. Therefore, we did not utilize the data from C_{3} in our analysis.
Gene set interactions
Genes in G1/S* are not regulated by genes in G2/M 
Genes in G1/S are not regulated by genes in M 
Genes in G2/M are not regulated by genes in G1/S 
Genes in G2/M are not regulated by genes in S 
Genes in M* are not regulated by genes in G1/S 
Rank table
Inducing a connection in HeLa cells data
The second highest rank produced is the connection CCNF← RAD21 determined as no relationship. As these two genes are classified as members of the same cell cycle phase, and the experimental determination of “no relationship” between two genes is difficult, we do not consider this result relevant.
Discussion
By looking at Figures 2, 3, 6 and 7, it is interesting to note that, in some cases, the frequency analysis of connections was capable of almost excluding one relationship possibility, transforming some undetermined connections into partially determined connections. These results show that the cell cycle pathway constrains some connections, therefore restricting the whole network [34].
We can attribute this phenomenon to the high dependency that the determination of a gene connection has on other connections. The proposed algorithm performs a search over the space of possibilities of the influence of a set of genes over a single gene. If one of these influences is a priori determined (or known), this result can bias other connections. For example, let us suppose that genes A and B have to produce a positive output over a gene C, according to some restriction imposed by the timeseries data. If we already know that gene A has no relationship to gene C, gene B must have a positive relationship on gene C. Therefore, this high dependency on the determination of a gene connection over the network makes the use of Figures 2, 3, 6 and 7 very restricted. If we simply use a relationship with a high weight to be our “best guess” on the connection between two genes, this choice can constrain other relationships, leading the system to a more or less determined state, or even creating a connection in a network that is not consistent with the data.
Another fact to be pointed out is the importance of the inferred partially determined connections. Although these connections can not be directly used to construct a network like the determined connections, it can guide some biological experiments, since a partially determined connection states that at least one type of relationship between two genes is not possible.
We could use the connection frequencies generated to attribute a strength of connection to the relationship of a partially determined connection, e.g., in the yeast cell cycle, the interference of Clb1 on SBF can be stated as 80% (or a probability of 0.8) of being an inhibition. In fact, we use the frequencies in Figures 2, 3, 6 and 7 to compute the rank of undirected relationships (Equation 11).
Regarding the validation, in the yeast cell cycle data, which is artificially generated, the true positive rate (Figure 5) of predicted connections is very high (75%) when 25 connections are considered, and 100% for 5 and 10 connections. These results show that the proposed algorithm can be successfully applied in an artificial data, i.e., a data set without noise and good balance between time points and number of genes. Evidently, our rank procedure is constructed in a way that determined connections and partially determined connections would be benefited. Hence, as our algorithm correctly determines 13 directed connections in the yeast cell cycle model, its true positive rate for a small number of undirected predictions is high as well.
Considering the HeLa cells data, the results are not quite optimistic. The true positive rate (Figure 8) stands between 17% and 30% for the inference procedure without biological knowledge and 25% to 35% using biological knowledge. However, if we consider the small amount of time points (12 per cell cycle) and the number of genes (20 in this simplified version), the difficulty of obtaining a higher true positive rate is clear.
Evidently, the method here proposed can only be used to aid a wet lab experiment on finding gene interactions if considerations about the network size and amount of timeseries data were made. In situations where a large set of g genes is investigated and only a small amount of timeseries data is available, as in the HeLa cells data, we would recommend that a rank of the r first interactions, with g ≪ r <g^{2}/2, to reduce the set of possible gene interactions to be tested.
To the HeLa cells, we can also explain the low true positive rate by considering the 20 genes version of the network too simplified. Maybe our algorithm predicts interactions that are not directly observed in nature, but only through a series of interactions of genes not present in our network. Therefore, our validation procedure is compromised.
It is worth to notice that adding some biological knowledge the results are improved for the HeLa cells. This fact reinforces the need for an integrated work with biologists in a network inference process, as we show that even using little pieces of biological information we can improve the whole procedure.
Regarding the example of finding a connection which most determines others, we expect to exemplify here that this use of the algorithm could substantially aid biological experiments. It is also worth noticing that the connection found, CCNF← RRM2 as an inhibition, could make biological sense, as RRM2 is classified as a gene of the S phase and CCNF is a G_{2} phase gene. As the connection CCNF← RRM2, defined as an activation, is also well ranked, we can say that this relation is worth for an empirical test.
A closer look at the frequency analysis raises another interesting question: would not the network chosen by nature be easily detectable? Or even better: would not the utilized data be enough to constrain the connection frequencies into nature’s choice? We could answer this question by pointing out a truth that unequivocally distinguishes our model from nature’s choice: the chemical interactions between proteins. Evidently, some of the connections considered on many steps of the algorithm here presented cannot exist due to chemical incompatibilities. In some sense, nature has more information to constrain its network than we do.
Conclusions
This paper proposes an algorithm to perform analyses for discovering gene regulatory interactions from timeseries data under the Boolean network model and in the context of Constraint Satisfaction Problem (CSP). In fact, the inference of gene regulatory network is a onetomany inverse problem in the sense that there may exist several networks consistent with the dataset. In order to analyze the gene interactions, we have generated several gene connections in consistent networks by using CSP solver techniques which in turn utilized constraints sets built from three algorithms provided by this work. We have applied our methodology to an artificial dataset that had been generated by a Boolean network that models the budding yeast cell cycle [35], and to an experimental dataset of HeLa cells [36]. By these applications, we have shown that our analyses could be a first step for detection of gene relationships with a high flexibility to include biological knowledge.
A challenge always presented in any gene regulatory model is its usefulness. It would be very interesting if a model could help biological experiments in understanding gene interactions. The model presented here together with the algorithm proposed is a first step to aid an inference process from timeseries data of gene expression, and it can be improved by all a priori knowledge available. As it was made clear in the HeLa cells data, the use of biological knowledge can improve the efficiency of the proposed algorithm. For future steps, an interesting feature to be improved on our method is the ability to indicate which connection should be verified in the wet lab to help determine others. As exemplified in this work, this feature could lead to important contributions on wet lab experiments. To use this method in an empirical gene connection survey, we would recommend a search over all possible connections between all genes, and then proceed with the ranking process. Evidently, biological considerations over the highest ranks produced is heavily necessary.
However, there are other characteristics to be sought that could constrain the network towards nature’s choice. On one hand, so far, only constraints built from successive states are considered. Thus, constraints constructed from considering the whole trajectory (e.g., some kind of powers of the regulation matrix in order consider succession of more than two or three states) could help to obtain more precise solutions. In fact, although we have carefully generated uniform samples to build a set of solutions (to produce representations good enough for connection frequencies) by using an appropriate heuristic for variable and value ordering, it is important to keep in mind that in order to make a more precise frequency analysis, one needs the consistent solutions in the CSP context, meaning that, in our case, the solutions obtained from considering the whole trajectory. On the other hand, one feature not explored in this paper is the dynamics of the network. There are indications, as stated by Kauffman [23], that nature would prefer networks with a small amount of attractors  the gene pattern expression that leads the system to itself and large basins of attraction  the set of gene pattern expressions that leads the system to an attractor. The network assembled by Li et al. [35] has these characteristics. Therefore, an analysis of connections computed only from networks with a few number of attractors  or other dynamical characteristic  could create a well established result. One naïve way to proceed is to build regulation matrices from the solutions of the CSPs subproblems (possible rows) and select the ones such that present the dynamical features described before (large basins of attraction and small number of attractors)
Concluding, the analysis presented here is a remarkable first step for the construction of a system to infer gene interactions. The true positive rate on the artificial data is excellent and, considering noises and lack of time points, the true positive rate for the experimental data is beyond expectation.
We understand that any inference procedure can not have success if it does not contain biological and computational expertise, therefore the future steps of this research tend to be centered on the difficulties of a wet lab, or its limitations.
Declarations
Acknowledgements
This work was supported by FAPESP (Fundação de Amparo à Pesquisa do Estado de São Paulo), CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico), and CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior).
This article has been published as part of BMC Proceedings Volume 5 Supplement 2, 2011: Proceedings of the 6th International Symposium on Bioinformatics Research and Applications (ISBRA'10). The full contents of the supplement are available online at http://www.biomedcentral.com/17536561/5?issue=S2.
Authors’ Affiliations
References
 Kauffman SA: Metabolic stability and epigenesis in randomly constructed genetic nets. Journal of Theoretical Biology. 1969, 22: 437467. 10.1016/00225193(69)900150.View ArticlePubMedGoogle Scholar
 Shmulevich I, Dougherty ER, Zhang W: Probabilistic Boolean networks: a rulebased uncertainty model for gene regulatory networks. Bioinformatics. 2002, 18 (2): 261274. 10.1093/bioinformatics/18.2.261.View ArticlePubMedGoogle Scholar
 Zhang SQ, Ching WK, Ng MK, Akutsu T: Simulation study in Probabilistic Boolean Network models for genetic regulatory networks. International Journal on Data Mining and Bioinformatics. 2007, 1 (3): 217240. 10.1504/IJDMB.2007.011610.View ArticleGoogle Scholar
 Friedman N, Linial M, Nachman I, Pe’er D: Using Bayesian networks to analyze expression data. Journal of Computational Biology. Journal of Computational Biology. 2000, 7: 601620. 10.1089/106652700750050961.View ArticlePubMedGoogle Scholar
 Goodwin BC: Temporal Organization in Cells; A Dynamic Theory of Cellular Control Process. 1963, New York: Academic PressGoogle Scholar
 Karlebach G, Shamir R: Modelling and analysis of gene regulatory networks. Nature. 2008, 9: 770780.Google Scholar
 Hartemink AJ: Reverse engineering gene regulatory networks. Nature. 2005, 5: 554555.Google Scholar
 Marshall S, Yu L, Xiao Y, Dougherty ER: Inference of a probabilistic Boolean network from a single observed temporal sequence. EURASIP Journal on Bioinformatics and Systems Biology. 2007, 32454Google Scholar
 Dougherty ER, Xiao Y: Design of Probabilistic Boolean Networks under the Requirement of Contextual Data Consistency. IEEE Transactions on Signal Processing. 2006, 54 (9): 36033613. 10.1109/TSP.2006.877641.View ArticleGoogle Scholar
 Xiao Y, Dougherty ER: Optimizing ConsistencyBased Design of ContextSensitive Gene Regulatory Networks. IEEE Transactions on Circuits and Systems I. 2006, 53 (11): 24312437. 10.1109/TCSI.2006.883883.View ArticleGoogle Scholar
 Hashimoto RF, Kim S, Shmulevich I, Zhang W, Bittner ML, Dougherty ER: Growing Genetic Regulatory Networks from Seed Genes. Bioinformatics. 2004, 20 (8): 12411247. 10.1093/bioinformatics/bth074.View ArticlePubMedGoogle Scholar
 Vahedi G, Ivanov I, Dougherty E: Inference of Boolean networks under constraint on bidirectional gene relationships. Systems Biology. 2009, IET, 3 (3): 191202.
 Liu W, Lahdesmaki H, Dougherty ER, Shmulevich I: Inference of Boolean Networks Using Sensitivity Regularization. EURASIP Journal on Bioinformatics and Systems Biology. 2008, 2008 (780541): 112. 10.1155/2008/780541.View ArticleGoogle Scholar
 Laubenbacher R, Stigler B: A computational algebra approach to the reverse engineering of gene regulatory networks. Journal of Theoretical Biology. 2004, 229 (4): 523537. 10.1016/j.jtbi.2004.04.037.View ArticlePubMedGoogle Scholar
 Laubenbacher R, Jarrah ASS: Algebraic models of biochemical networks. Methods in Enzymology. 2009, 467: 163196. full_text.View ArticlePubMedGoogle Scholar
 Jarrah A, Laubenbacher R, Stigler B, Stillman M: Reverseengineering of polynomial dynamical systems. Advances in Applied Mathematics. 2007, 39: 477489. 10.1016/j.aam.2006.08.004.View ArticleGoogle Scholar
 Stigler B, Jarrah A, Stillman M, Laubenbacher R: Reverse engineering of dynamic networks. Annals of the New York Academy of Sciences. 2007, 1115: 168177. 10.1196/annals.1407.012.View ArticlePubMedGoogle Scholar
 Sima C, Hua J, Jung S: Inference of Gene Regulatory Networks Using TimeSeries Data: A Survey. Current Genomics. 2009, 10 (14): 416429. 10.2174/138920209789177610.PubMed CentralView ArticlePubMedGoogle Scholar
 Hecker M, Lambeck S, Toepfer S, van Someren E, Guthke R: Gene regulatory network inference: Data integration in dynamic models  A review. BioSystems. 2009, 96: 86103. 10.1016/j.biosystems.2008.12.004.View ArticlePubMedGoogle Scholar
 Dougherty ER: Validation of Inference Procedures for Gene Regulatory Networks. Current Genomics. 2007, 8 (6): 351359. 10.2174/138920207783406505.PubMed CentralView ArticlePubMedGoogle Scholar
 Shmulevich I, Dougherty ER, Zhang W: From Boolean to Probabilistic Boolean Networks as Models of Genetic Regulatory Networks. Proceedings of the IEEE. 2002, 90: 17781792. 10.1109/JPROC.2002.804686.View ArticleGoogle Scholar
 Huang S: Gene expression profiling, genetic networks, and cellular states: An integrating concept for tumorigenesis and drug discovery. Journal of Molecular Medicine. 1999, 77: 469480. 10.1007/s001099900023.View ArticlePubMedGoogle Scholar
 Kauffman SA: The Origins of Order: SelfOrganization and Selection in Evolution. 1993, USA: Oxford Univ. PressGoogle Scholar
 Tsang E: Foundations of Constraint Satisfaction. 1993, Academic PressGoogle Scholar
 Russell S, Norvig P: Artificial Intelligence: A Modern Approach. 2002, Prentice Hall, secondGoogle Scholar
 Tack G: Constraint Propagation  Models, Techniques, Implementation. PhD thesis. 2009, Saarland UniversityGoogle Scholar
 Schulte C, Tack G, Lagerkvist MZ: Modeling. 2010Google Scholar
 Jiang D, Tang C, Zhang A: Cluster Analysis for Gene Expression Data: A Survey. IEEE TKDE. 2004, 16 (11): 13701386.Google Scholar
 D’haeseleer P, Liang S, Somogyi R: Genetic network inference: from coexpression clustering to reverse engineering. Bioinformatics. 2000, 16 (8): 707726. 10.1093/bioinformatics/16.8.707.View ArticlePubMedGoogle Scholar
 Dougherty ER, Barrera J, Brun M, Kim S, Junior RMC, Chen Y, Bittner ML, Trent JM: Inference from Clustering with Application to GeneExpression Microarrays. Journal of Computational Biology. 2002, 9: 105126. 10.1089/10665270252833217.View ArticlePubMedGoogle Scholar
 Zhou X, Wang X, Dougherty ER, Russ D, Suh E: Gene Clustering Based on Clusterwide Mutual Information. Journal of Computational Biology. 2004, 11: 147161. 10.1089/106652704773416939.View ArticlePubMedGoogle Scholar
 Hashimoto RF, Dougherty ER, Brun M, Zhou Z, Bittner ML, Trent JM: Efficient Selection of Feature Sets Possessing High Coefficients of Determination Based on Incremental Determinations. Signal Processing. 2003, 83 (4): 695712. 10.1016/S01651684(02)004681.View ArticleGoogle Scholar
 Peña J, Björkegren J, Tegnér J: Growing Bayesian network models of gene networks from seed genes. Bioinformatics. 2005, 21 (Suppl. 2): ii224ii229.PubMedGoogle Scholar
 Lau KY, Ganguli S, Tang C: Function Constrains Network Architecture and Dynamics: A Case Study on the Yeast Cell Cycle Boolean Network. Physics Review E. 2007, 75 (5): 19. 10.1103/PhysRevE.75.051907.View ArticleGoogle Scholar
 Li F, Long T, Lu T, Ouyang Q, Tang C: The Yeast CellCycle Network is Robustly Designed. PNAS of the USA. 2004, 101 (14): 47814786. 10.1073/pnas.0305937101.View ArticleGoogle Scholar
 Whitfield ML, Sherlock G, Saldanha AJ, Murray JI, Ball CA, Alexander KE, Matese JC, Perou CM, Hurt MM, Brown PO, Botstein D: Identification of genes periodically expressed in the human cell cycle and their expression in tumors. Molecular biology of the cell. 2002, 13 (6): 19772000. 10.1091/mbc.02020030..PubMed CentralView ArticlePubMedGoogle Scholar
 Koff A, Giordano A, Desai D, Yamashita K, Harper JW, Elledge S, Nishimoto T, Morgan DO, Franza BR, Roberts JM: Formation and activation of a cyclin Ecdk2 complex during the G1 phase of the human cell cycle. Science. 1992, New York, N.Y., 257 (5077): 16891694. 10.1126/science.1388288.
 Shanahan F, Seghezzi W, Parry D, Mahony D, Lees E: Cyclin E associates with BAF155 and BRG1, components of the mammalian SWISNF complex, and alters the ability of BRG1 to induce growth arrest. Molecular and Cellular Biology. 1999, 19 (2): 14601469.PubMed CentralView ArticlePubMedGoogle Scholar
 Chan AK, Persad S, Litchfield DW, Wright JA: Ribonucleotide reductase R2 protein is phosphorylated at serine20 by P34cdc2 kinase. Biochimica et Biophysica Acta (BBA)  Molecular Cell Research. 1999, 1448 (3): 363371. 10.1016/S01674889(98)001153.View ArticleGoogle Scholar
 DeGregori J, Kowalik T, JR N: Cellular targets for activation by the E2F1 transcription factor include DNA synthesis and G1/Sregulatory genes. Molecular Cell Biology. 1995, 15 (8): 42154224.View ArticleGoogle Scholar
 Chiang SY, Azizkhan JC, Beerman TA: A comparison of DNAbinding drugs as inhibitors of E2F1and Sp1DNA complexes and associated gene expression. Biochemistry. 1998, 37 (9): 31093115. 10.1021/bi9721142.View ArticlePubMedGoogle Scholar
 Kong M, Barnes EA, Ollendorff V, Donoghue DJ: Cyclin F regulates the nuclear localization of cyclin B1 through a cyclincyclin interaction. The EMBO Journal. 2000, 19 (6): 13781388. 10.1093/emboj/19.6.1378.PubMed CentralView ArticlePubMedGoogle Scholar
 Holt LJ, Krutchinsky AN, Morgan DO: Positive feedback sharpens the anaphase switch. Nature. 2008, 454 (7202): 353357. 10.1038/nature07050.PubMed CentralView ArticlePubMedGoogle Scholar
 Cai J, Gibbs E, Uhlmann F, Phillips B, Yao N, O’Donnell M, J H: A complex consisting of human replication factor C p40, p37, and p36 subunits is a DNAdependent ATPase and an intermediate in the assembly of the holoenzyme. Journal of Biological Chemistry. 1997, 272 (30): 1897418981. 10.1074/jbc.272.30.18974.View ArticlePubMedGoogle Scholar
 Ohta S, Shiomi Y, Sugimoto K, Obuse C, Tsurimoto T: A proteomics approach to identify proliferating cell nuclear antigen (PCNA)binding proteins in human cell lysates. Identification of the human CHL12/RFCs25 complex as a novel PCNAbinding protein. Journal of Biological Chemistry. 2002, 277 (43): 4036240367. 10.1074/jbc.M206194200.View 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.