Hybrid-Controlled Neurofuzzy Networks Analysis Resulting in Genetic Regulatory Networks Reconstruction

Reverse engineering of gene regulatory networks (GRNs) is the process of estimating genetic interactions of a cellular system from gene expression data. In this paper, we propose a novel hybrid systematic algorithm based on neurofuzzy network for reconstructing GRNs from observational gene expression data when only a medium-small number of measurements are available. The approach uses fuzzy logic to transform gene expression values into qualitative descriptors that can be evaluated by using a set of defined rules. The algorithm uses neurofuzzy network to model genes effects on other genes followed by four stages of decision making to extract gene interactions. One of the main features of the proposed algorithm is that an optimal number of fuzzy rules can be easily and rapidly extracted without overparameterizing. Data analysis and simulation are conducted on microarray expression profiles of S. cerevisiae cell cycle and demonstrate that the proposed algorithm not only selects the patterns of the time series gene expression data accurately, but also provides models with better reconstruction accuracy when compared with four published algorithms: DBNs, VBEM, time delay ARACNE, and PF subjected to LASSO. The accuracy of the proposed approach is evaluated in terms of recall and F-score for the network reconstruction task.


Introduction
Biological systems are inherently stochastic, uncertain, and fuzzy [1]. Therefore, research in bioinformatics and computational biology, where computer technologies are applied to manage and analyze biological data and make computational models, is faced with a great deal of uncertainty. For instance, growth and development as well as environmental stresses can all contribute to change in gene expression levels. In addition, under such conditions, some genes influence the expression of other genes and their functionalities.
With the advent of high-throughput technologies in transcriptomics, proteomics, and metabolomics, now, biologists have the ability to investigate the expression of genes and consequences on a genome-wide scale. Gene expression data in the form of high-throughput microarray experiments measure the amounts of RNA associated with each of thousands of genes in parallel. Time-series microarrays have attracted biologists' interests for deciphering the dynamic and complex nature of biological networks. Time-series microarrays record multiple expression profiles at discrete time points (i.e., hours or days) of a continuous cellular process. Thus, analytical methods are needed to handle many genes with uncertain functions based on discrete datasets of continuous biological processes. The methodological areas range from experimental design [2] to data normalization [3,4], missing value imputation [5], cluster analysis [6,7], classification [8], identification of differentially expressed genes [9], and network modelling [10,11].

ISRN Bioinformatics
As a challenging concept, reverse engineering can be employed to estimate gene regulatory networks from highthroughput expression data. Genetic regulatory network reconstruction provides a concise representation of the interactions between multiple genes at the system level. In addition, it confers a broader insight for biologists about the manner in which genes interact with one another, and about the roles that they play in various biological functions.
Viability of an organism down to the cells is essentially controlled by gene expression regulation at the transcript level. This concept leads one to ponder how the changes in the expression patterns of genes during the ordinary and the stressful conditions of cells may infer the way genes are affected by environmental conditions (e.g., lack of nutrients) [12]. Such studies of gene expression patterns can improve our understanding of biological systems, and may enhance our ability to combat undesired situations (e.g., diseases such as cancer) with the hope of improving human life quality.
There are several limitations in the study of time-series gene expression data such as small sample size (due to the time-consuming nature in which samples are produced and the high costs associated with microarray experiments, especially in clinical studies), genes with low level expression and noisy data structure [13,14]. Problems related to high dimensionality accompanied by a small sample size, such as matrix singularity, model over-fitting and model overparameterization become more pronounced in the case of the most available data [15]. Also, the unavoidable presence of noise has more influence on the analysis of short term rather than long term data series. This enhances the difficulty in distinguishing actual patterns from random data, thereby raising the potential of misleading analyses [16].
Reconstruction of gene regulatory networks based on expression data should therefore: be able to handle constrained data; should be robust to, and compensate for noise and incompleteness of data; and should be capable of providing interpretable results.
Modelling based on Boolean Networks is one of the common methods employed in GRNs inference [32]. The goal of these models is basically to infer rules based on their computational simplicity and ability to handle noisy experimental data [17]. Even though these models can be easily applied, much information is lost in binary encoding and, in practice, the derived models have insufficient dynamic resolution because they depend on arbitrary discretizations of the gene expression values [33][34][35][36][37]. In general, Boolean networks are limited by their definition.
Bayesian network modelling is based on probabilistic transitions between network states and assumes that there is no feedback in a network; in spite of the fact that cycles of events are the major mechanism to ensure robustness of the biological systems [21].
HMM also has been applied for analyzing time-series gene expression data [24]. However, there are several problems with HMMs. The number of parameters that need to be set in an HMM is quite high. As a result, the amount of data required to train an HMM is very large. Also, concepts learnt by an HMM are framed in terms of emission and transition probabilities. If one is trying to understand the concept learnt by the HMM, then this concept representation is difficult to understand.
Dynamic Bayesian Network (DBN) as another method combines the features of hidden Markov models to incorporate feedback [38]. Models based on Bayesian networks, despite attractiveness due to their ability to deal with stochastic aspects of gene expression and noisy measurements, have the disadvantage of minimizing the dynamical aspects of gene regulation [20]. A graphical Gaussian model (GGM) is an undirected probabilistic graphical model [25]. This model allows the identification of conditional independence relations among genes, under the assumption of a multivariate Gaussian distribution of gene expression data. The GGM does not identify the direction of gene relationships, but rather only calculates the correlations between their gene expression data.
The Kalman filter (KF) [28] is only applicable to linear models and the Gaussian posterior density probability. Since position information is linear, standard Kalman filtering can be easily applied to the tracking problem without much difficulty. However, gene regulatory networks pose nonlinear information, requiring a modification to the KF. To overcome this problem, much research has been reported on nonlinear filtering methods such as extended Kalman filter (EKF) [29], unscented Kalman filter (UKF), and Particle filter (PF) [30]. Currently, considerable research is being devoted to introduce improvements in the working of these algorithms and enhance our understanding about gene interactions.
In this paper, we describe a novel algorithm that benefits from using rule-based neurofuzzy networks (RBNFNs) to extract information from time-series gene expression data. The suggested algorithm combines neural networks with fuzzy systems and allows for mapping the dynamics of gene expression data. Fuzzy logic [39][40][41] and artificial neural networks [42,43] are complementary technologies in the design of an intelligent system, and their combination appears to be a promising path, since neural networks are essentially lowlevel, computational algorithms that sometimes offer a good performance in pattern-recognition tasks; whilst fuzzy logic provides a structural framework that uses and exploits those low-level capabilities of neural networks. Thus, the combination seems to offer potential for capturing subtle effects of genes. Neural networks can learn from data sets while fuzzy logic solutions are easy to verify and optimize. Fuzzy logic and neural networks generally approach the design   of intelligent systems from quite different angles, and the combined system can have advantages from both sides: neural networks are implicit; although the system is not easily interpreted or modified, it trains itself by data sets. Fuzzy logic is explicit; thus the system verification and optimization is more efficient.
Here, we present a two-level algorithm based on RBNFNs for extracting gene interaction networks. This algorithm is advantageous as the network is never overparameterized, avoids redundant fuzzy sets, decreases the redundancy of the model, thus it is simpler and drastically reduces the computational cost. Also, in this network, the training process does not depend on the number of inputs and the sample size. The proposed networks are trained with a subset of experimental samples and tested with the remaining samples. Then, the constructed fuzzy functions in the gene interaction reconstruction (referred as edges) depend on the dynamics of input-output patterns of the networks. We applied the network to yeast cell-cycle regulation data and the resulted interactions from our model were compatible with previous experimentally verified interactions.

Algorithm
The method proposed in this paper is a quantitative computational approach consisting of five main stages shown in Figure 1. Our methodology includes two levels. In the first level (stage 1), RBNFNs are used for time series prediction of gene expression; in the second level (stage 2 to 5), the rules created by the RBNFNs for reconstructing gene regulatory networks are employed. The detailed descriptions of these stages are outlined as follows: (1) Training RBNFNs by gene expression data and creating IF-THEN rules (Section 2.1),  . Expression data for each gene is normalized (input layer) before being fuzzified using a three-state Gaussian matrix (fuzzier layer), then the preceding part of fuzzy rules are combined using a product operator (rule base layer), before the rules are normalized (rule normalization layer) and the fuzzy target value is defuzzified (combination & defuzzifier layer).
(3) Sorting the obtained fuzzy rules of each RBNFN to have further analyses (Section 2.3), (4) Extracting similarity rules as a connectivity matrix, from sorted fuzzy rules of all RBNFNs, and generating the interactions from each rule (Section 2.4), (5) Analyzing the extracted interactions of stage 4 to reconstruct the final model (Section 2.5).

Rule-Based Neurofuzzy Networks (RBNFNs).
Let us first assume that there are N genes in the expression dataset we are studying. In the first stage, N independent RBNFNs are generated with the following structure: in each RBNFN, the expression data for all genes are loaded as inputs except for one gene, which is considered as an output. Training each RBNFN with expression data gives rise to a network capable of predicting the expression pattern of the output as the result of input genes. RBNFNs are employed in this stage because of their ability to overcome the drawbacks of pure neural networks. By incorporating elements of fuzzy reasoning processes, the RBNFN gives meaning and function as part of a fuzzy rule to each node via an associated weight. We have utilized the least square (LS) technique to learn the parameters of these RBNFNs. The architecture of a typical RBNFN (illustrated in Figure 2) includes the following steps:

Input Layer (Step 1).
In this layer, the input expression dataset is normalized using (1). Each node, corresponding to each input variable, normalizes each input value to the scale of [0, 1] for the next layer, to facilitate in fuzzification: where G i is the expression values of the ith gene, and O (1) i is the ith output of input layer.

Fuzzier Layer (
Step 2). In layer 2, the normalized expression values are fuzzified using a three-state Gaussian matrix, with linguistic values of "Low," "Medium," and "High." These linguistic labels represent the data in fuzzy logic terms. Figure 3 depicts the Gaussian membership functions (MFs) for the three-state model employed to find Table 1: Membership matrix of a three-state fuzzy logic model.

Labels
Genes Entries within the columns are the normalized expression value of genes in the kth time point applied to the membership functions of low, medium, and high expressed. the maximum degree that each sample of normalized input belongs to the respective label. In (2), the membership functions are represented in Gaussian form: where i, k, andjindicate the label of the input variable, the time point and the MF, respectively, m and σ are the mean and variance of the membership functions, while O (1) i,k is the value of the ith input variable in the kth time point.
We have considered constant values for the means and variances of MFs in a manner that they cover the scale of [0, 1] in equal partitions. The constant means and variances are assumed in order to maintain an ability to compare and analyze the extracted rules in the next layers. In other words, if these values are different for each input, the concept of low, medium, and high expressed for each gene will not be the same as for the others. In addition, employing constant means and variances allows for fewer parameters to be trained in the network; an advantage of which is that, when encountering small sample size temporal genetic data, the network will not be overparameterized during training.
As described earlier, the algorithm in step 2 partitions the input/output space. This space is an n-dimensional unit hypercube [0, 1] n similar to Table 1. This hypercube results when membership functions map the points in the input/ output space to a degree of membership between 0 and 1. Table 1 presents the membership value of the normalized expression data of each gene assigned to low, medium, and high labels by (2). In this table, O j i,k represents the membership value of kth data point for ith gene assigned to jth label.

Rule Base Layer (
Step 3). Upon the fuzzified output of layer 2, a rule-based profile is set up in layer 3. This profile consists of fuzzy IF-THEN rules for determining the expression value of a target gene according to the expression levels of input genes. For instance, as illustrated in Figure 4 for a sample three-state model, a fuzzy IF-THEN rule may be presented as follows: IF "the expression level of G 1 , G 2 and G 3 are respectively Low, High and Medium as input genes," THEN "the corresponding expression level for the target gene, G T , will be Low." The method used in our fuzzy system is based on a singleton output function, which assigns a single value to each of the N fuzzy states in the model.
In order to reduce the search space of possible fuzzy rule combinations from the rule-based layer, we only create new rules associated with output states observed in the training sample set. As such, each new rule corresponds to interaction clusters discovered in the input/output space. This indicates that the number of created rules in the RBNFN depends on the number of sample sets of input/output genes. Thus, the more samples of genes presented, the greater number of rules created. The expected benefit of this approach is that there would be less contamination in knowledge extracted by not considering irrelevant states for gene interaction.
When a new set of input/output data is applied to the RBNFN, the network determines whether to generate a new rule for describing the incoming pattern (x, y) or not. This occurs after checking the similarities between the new rule and the previous ones. This process leads to a reduction in the number of fuzzy sets and avoids the existence of redundant rules.
Each node in layer 3 combines the antecedent part of a fuzzy rule using a T-norm operator. In this study, the Tnorm operator is considered as a product operation. The output of each node represents the firing strength of the corresponding fuzzy rule, which is calculated by 6 ISRN Bioinformatics where r runs through all the selected nodes of step 2 corresponding to the kth rule.

Rule Normalization Layer (
Step 4). The number of nodes in this step is the same as step 3. The nodes in layer 4 calculate the ratio of the ith rule's firing power to the sum of all rules' firing strengths, which can be formulated by (4). Indeed, the rules are normalized in this layer in the scale of where m is the number of all rule sets. The output of this layer is also called normalized firing power.

Combination and Defuzzification Layer (Step 5).
After applying the decision matrix to the fuzzified expression levels in step 4, the membership degree of a target to a statement is determined. At step 5, the fuzzy target value is transformed back into a value between 0 and 1 via the process of defuzzification; indeed, the output is the predicted value of the RBNFN.
Equation (5) formulates the output of layer 5, as the overall output of RBNFN. In fact, this equation combines and defuzzifies all the rules from previous steps where w i is the weight multiplied to the ith rule. Figure 5 illustrates this process as an example. In this figure, each row indicates a rule calculated from normalized expression values of input genes at a sample time point.
As described above, in our methodology, we have an RBNFN for each gene. In a particular RBNFN, the expression value of the output gene is assumed to be the outcome of all other genes. In other words, all remaining genes are considered to be either activators or repressors for the target gene, and the predicted expression pattern for this gene is deduced from the expression levels of all other ones.
Also we compare the prediction of the RBNFN with the real expression pattern of a target gene in order to calculate the weights of the network. A greater weight specifies that the corresponding input gene has a greater affect on the output, so it can be considered as a greater activator or repressor effect on the target gene.

Rule
Gene

RBNFNs Learning
Procedure. The parameters of proposed RBNFN must be trained to obtain a prediction system for our data. Therefore, in the learning procedure, the related parameters of rules (weights) are trained by LS the least squares (LS) learning algorithm. The simplicity of LS algorithm makes it a common and suitable method for training and tuning parameters of neural or fuzzy networks. An error criterion is defined for each RBNFN in order to test the accuracy of a given data set. The error criterion is an overall error measurement based on the differences between the predictions and target values. We used mean square error (MSE) as our error criterion which is a very common criterion described by where N is the number of samples in the data set. Intuitively, an increase in the number of rules results in a decrease in the error measure of (6).
In the next step, we attempt to extract the most of available information from provided genetic data in order to make the most possible understanding of behaviour and interaction among genes.

Extracting Fuzzy IF-THEN Rules.
As stated above, in the proposed HRBNF algorithm, we construct network models for as many genes as possible from the data available. In each model, one of the genes is considered as the output and the others as the input nodes. The model attempts to describe the behaviour of the output gene based on the possible interactions of input genes.
Finally, it stores the derived knowledge from gene-gene interactions in the fuzzy rules and their corresponding weights. Table 2 presents the extracted rules of ith network. In this table, i, k run through all the input variables and extracted rules, respectively, for n genes. L i k, j indicates the fuzzy label of jth gene in kth rule of ith RBNFN. Table 2, the rules are sorted in a way that the ith gene takes the ith position in the matrix; so we shift it between i − 1th and i + 1th column of the matrix. This leads us to have similar rule matrices in all RBNFNs and a decision agent to find gene interactions. The interactions are delineated according to the weight values and upon comparing them in similar rules of different RBNFNs. As stated in previous sections, to extract more precise relationships 8 ISRN Bioinformatics Table 3: The networks effects on jth rule.

Rule
Gene · · · L n j, n w n j between genes, a full capacity of all available data is used in HRBNF algorithm. The proposed analysis of obtained rules and the way interactions among genes occurs are described in the next stage.

Similar Rules
Extraction. Assuming that we have N genes, the outputs of previous stages are N matrices with K similar rules in each one. In this stage, we create a new set of matrices from previous tables in order to extract information. To achieve this target, N matrices are defined with K rules in them. A new table in this stage consists of K similar rules gained from all RBNFNs. Regarding the weights of a specified rule in RBNFNs, we can find the effect of that particular rule in every RBNFN; therefore, its effect on every output gene can be evaluated. This guides us to extract the relationships between genes by exploiting these matrices. In Table 3, an example to the jth rule is illustrated. We call this table as the table of network effects on the similar rule (NESR) matrix.
In order to estimate gene interactions from NESR matrix, the weights have to be analyzed. Therefore, we define some thresholds to identify the effect of a rule on all genes.
It is not easy to determine an optimal threshold value since a large value may result in the elimination of true connections with small effects during the procedure. On the other hand, a small threshold is not able to retrieve true connections from false connections in an inherently sparse genetic dataset. Our proposed criterion to obtain an appropriate threshold is called "effectivity threshold (ET)." This threshold is different for each RBNFN and is defined as the average of all weights calculated in that RBNFN, as shown in (7). We defined different thresholds for RBNFNs because the range of weights in each RBNFN is different from others; thus, by having different thresholds, a more realistic view of the effect of similar rules in different networks is provided In (7) i, k denote ith RBNFN and the number of rules in ith RBNFN, respectively. Also, w i p is the calculated weight value for pth rule in ith RBNFN.
After defining the thresholds, ETs are compared with thresholds. We also define a parameter named effectivity symbol (ES), which describes the state of each ET. According to (8), an ES equals to "+1" when the related weight is positive with absolute value greater than the defined threshold, or equals to "−1" when the related weight is negative with absolute value bigger than the defined threshold. An ES is considered to be "0" if the absolute value of related weight is less than the threshold In the above equations, i, j denotes ith RBNFN for jth rule; and ES i j shows whether the jth rule has an effect more than the threshold on the ith gene or not.
By defining ESs, new matrices are obtained from NESR matrices with ES values instead of weights. We named these matrices as "ES matrices." Table 4 is an instance of an ES matrix related to jth rule. This matrix contains similar jth rule extracted from different RBNFNs. This stage is finalized by a decision-making process, in which we extract the connections from ES matrices. The connection of two genes is displayed by a directed graph in which an edge from one gene is directed to the other one, showing that the former gene has effect on the latter one. Our proposed procedure of extracting the connections among genes from the ES matrix is outlined in the following.
As stated before, a rule with ES equals to "0" implies that the weight of this rule in the related RBNFN is less than defined threshold; this means that the input genes of the RBNFN have effects less than the threshold on the output gene. Thus, we can consider that the output gene is not affected by the input genes. In other words, a rule with ES = "0" in an RBNFN means there is no inferred connection from the input gene to output gene. On the other hand, a rule with ES = "1" in an RBNFN means that the input genes have effects on the output gene as activators. Finally, a rule with ES = "−1" in an RBNFN indicates that the input genes are repressors of the output gene.
The decision-making process is formalized by (9):

Final Genetic Interaction Network Extraction.
In the former section, we presented a procedure for extracting the connections among genes from ES tables. This procedure leads us to create two connection networks (CNs) for each ES matrix: one for activator connections and the other for repressor connections. As previously discussed, ES matrices are obtained from the same rules extracted from all RBNFNs. Thus, the number of CNs attained from a genetic dataset is two times that of the extracted rules. In this section, we attempt to achieve a final genetic interaction network from these CNs. This goal is acquired by analyzing and interpreting the activator and repressor connections of all rules. We define two criteria for all connections of different CNs: activator criterion (AC) and repressor criterion (RC). Activator criterion, as shown in (10), has been employed to survey the activator interactions; as a result of which, Table 6 is obtained from all connections in CNs For Rule j ⇒ if G(i) connect to G p p = 1, . . . , n, p / = i and connection is activator In (10), AC = 1 denotes the existence of an activator connection and AC = 0 does not indicate whether an activator connection exists or not. The interaction is estimated by summation of ACs as shown in Table 5. In this estimation, first we calculate the summation of ACs that confirm an interaction from ith gene to pth gene ( k j=1 AC j i → p ), and the summation of ACs that confirm an interaction from pth gene to ith gene ( k j=1 AC j p → i ); then, the interaction between ith gene and pth gene is estimated as the subtraction of these two values. By estimating the interactions between all pairs of genes, the final activator interaction network can be ascertained. Similar to the procedure described above, we define a repressor criterion in order to obtain the repressor connections and to achieve the final repression interactions network. The repressor criterion is presented in (11) and the repressor interactions are illustrated in Table 6 For Rule j ⇒ if G(i) connect to G p p = 1, . . . , n, p / = i and connection is repressor where RC denotes the repressor connections. In Table 6, similar to Table 5, the repressor interaction between ith gene and pth gene is estimated as the subtraction of k j=1 RC j i → p and k j=1 RC j p → i . The final repressor interaction network will be provided by estimating the repressor interactions between all pairs of genes.

Results
In this section, we present the results of HRBNF algorithm for an experimental data. In Section 3.1, we will introduce yeast (Saccharomyces cerevisiae) cell cycle microarray time series data sets presented in [47,48]. This data has been extensively exploited for both practical and academic applications. Researchers frequently use these data sets to demonstrate and validate statistical and clustering analysis (e.g., [49][50][51][52]), mathematical modelling [37,38], and reverse engineering methods [36,53].
In Section 3.2, we present the results of RBNFNs performance in predicting the time series data sets. Finally, in Table 6: Resultant of repressor interactions.

Gene
Gene

Data: Yeast Cell Cycle Dataset.
In order to evaluate HRBNF algorithm, we generated fuzzy gene networks based on yeast cell cycle microarray time series data sets. We focused on twelve yeast genes playing key roles in the control of cell cycle as listed in Table 7 with descriptions taken from the Yeast Proteome Database [46]. The protein-protein and the regulatory interaction of the coded proteins for these genes which are involved in yeast cell-cycle are well-studied. The high-throughput techniques such as microarray provided us with the time-series expression of these genes reflecting the dynamic behaviour of these genes in cell cycle.
Although the new techniques such as RNA-seq are already out there which give us better resolution of expression profile, the main pattern of genes expression profile can be extracted from microarray expression data. Consequently, HRBNF algorithm is used as a "reverse engineering" method to find all possible genetic interaction network models that fit the data for the set of twelve genes and to demonstrate its ability to handle other similarly noisy data sets. In addition to the vast studies on the yeast data sets specially in systems biology that allow us to have a better assessment on the acquired results, these data sets are advantageous because of having adequate number of genes and time points for testing results.
As described by Cho et al. [48], gene expression profiles for yeast cell cycle have been studied through four microarray time series data sets: Alpha, Cdc15, Cdc28, and Elu with 18, 24, 17, and 14 time points, respectively. We have used Alpha, Cdc15, and Cdc28 data sets, in which still some samples were missed. The missing values are computed using an estimation method based on the K-Nearest Neighbourhood (KNN) algorithm [54].
S. cerevisiae cell cycle regulatory protein-DNA interactions were also the subject of a recent extensive experimental study [55] for which a great deal of information has been compiled in KEGG pathway database [44,45]. Figure 6, shows the interactions of the cell cycle regulatory protein subset shown experimentally so far [56]. We applied HRBNF algorithm on gene expression values corresponding to these twelve genes, and compared to the estimated directed network with the pathway depicted in KEGG.

Prediction of Gene Expression Level.
In order to constrain the number of parameters, we have to restrict the structure of our RBNFNs according to the size of employed data set. As a result, we consider the following hypothesis for the RBNFNs: (i) focusing on 12-key yeast cell cycle genes as the input/output nodes of the RBNFNs; (ii) setting the number  Figure 6: KEGG yeast cell cycle interactions-Schematic of known yeast cell cycle interactions between protein products of twelve studied genes regulatory (Table 7). An arrow indicates a positive interaction and a closed circle indicates negative interaction [44,45].
of MFs to three, as the minimum meaningful resolution. The centers and variances of search space with three MFs are considered as follows (identified empirically): The centers are selected according to the considered three MFs and linguistic values of "Low," "Medium," and "High"; the variances are obtained from results of Table 8. This table presents the accuracy of extracted interactions compared with real interactions, and shows that variance of 0.25 can extract the best result. Also this variance represents a distribution over the entire range of our search space. The RBNFNs are trained using Cdc15 dataset, and tested by Cdc28 and Alpha datasets. The test results are shown in Table 9.

Genetic Interaction Network Reconstruction for Yeast Cell
Cycle Data. To test the capability of our proposed method, we used the RBNFNs to predict time-series gene expression values and then, extract gene regulatory network (GRN) structures from inferred rules. The results of GRN reconstruction were evaluated by a part of yeast cell cycle  [44,45]. Figure 7 illustrates the extracted genetic network of activator and repressor interactions.  Figure 7: Inferred genetic interaction network for twelve yeast cell cycle regulatory genes. Each node represents a gene and the presence of an edge between the two nodes represents the existence of interaction between the two genes. Symbols " → " and " ", shown by blue and red edges, illustrate activator and repressor interactions, respectively. Dashed edges represent interactions that have been verified. In contrast, dotted edges are incorrect extracted interactions.

Discussion
We compared genetic interaction networks constructed based on five algorithms from the same dataset. The results show that the inferred network from our proposed algorithm has the best performance.
As shown in Table 10, the extracted network from HRBNF algorithm demonstrates a more complete match with the KEGG pathway than proposed methods from literature. Indeed, we observe that our HRBNF algorithm is capable to extract 13 true connections out of 33 available experimentally illustrated connections, while only 4, 5, 10, and 7 true connections are captured by DBN [57], VBEM [58], Time delay-ARACNE [59], and PF subjected to LASSO [30], respectively. It is clear that we reduced the misdirected edges.
Some criteria are useful to evaluate the goodness of fit of the inferred network [60]: the proportion of recovered true edges in the target network that is called Recall and precision corresponds to the expected success rate in the experimental validation of the predicted interactions.
We can compute sensitivity and precision from following equations: where true positive (TP) is the inferred number of edges identified correctly, false negative (FN) is the number of edges that were not identified, and false positives (FP) is the number of edges identified incorrectly.
Also, we calculated the F-score by using (14) to further quantify the performance of the algorithms [60]: where α is a weighting factor and here we consider α = 0.5 that is called the harmonic mean of precision and sensitivity because the importance of precision and sensitivity is even. Consequently, the goodness of fit of the results based on HRBNF and the other structure learning approaches in predicting the connectivity network of the KEGG pathway were compared using estimates of above criteria.
The results, presented in Table 11, show that evaluation criteria, sensitivity, and F-score, are distinctively higher for our approach compared to previous methods indicating the efficiency of our approach. In a good system, precision decreases as sensitivity increases.
Considering the results obtained from proposed method, it can be concluded that there are considerable agreements between the findings of HRBNF algorithm and experimental results reported in the literature. The main advantage of proposed method is that HRBNF algorithm searches for the relationships that fit to a logical understanding of how a set of genes should interact. By using the same criteria that biologists would use to describe the gene regulatory function (i.e., framed as an "IF-THEN-ELSE" relationship in terms of expression level), HRBNF algorithm based on fuzzy logic approximates the thought process that an expert would use to analyze these kinds of data; however, in contrast to an expert, HRBNF algorithm is automated and unbiased. If gene expression data is not analyzed properly, it can be difficult to interpret and can easily be misconstrued. In comparison to other methods, our computational algorithm analyses the data more efficient, unbiased, and fast. Our proposed model and Time delay-ARACNE [59] used less than one second against PF subjected to LASSO [30] which in turn used 23 minutes and 37 seconds on a Core i7 PC with 8 GB main memory, respectively. We speculate that this is due to the reduction in search space that results from by applying relevant rule selection (Section 2.1.3) prior to neurofuzzy network learning.
Algorithm based on neurofuzzy networks found a disproportionately large number of interactions for the roles of activators and repressors due to available datasets with small 14 ISRN Bioinformatics size samples. Of course, if expression profiling technologies become more sensitive and faster, our HRBNF algorithm would detect a more precise effect of gene expression on the activator and repressor roles, providing more information to decrease the number of interactions, which leads to a better understanding.

Conclusions
In this paper, we described a novel algorithm based on RBNFNs for gene regulatory network reconstruction. We demonstrated our approach by developing RBNFN models that more accurately predict gene expression data from typically noisy microarray experiments. Looking at the gene regulatory networks provides us a systemic view at the "gene interaction" level.
Our HRBNF algorithm was successfully validated using yeast cell cycle data. The rules of inferred networks reflect most interactions previously identified by genome-scale analysis and match the existing literature. At the network level, the inferred rules provide more detailed information about genes and the interactions among them. Potential new interesting interactions were identified, which provide novel hypotheses for new lines of further research. As a consequence, common rules among all the RBNFNs and the plausible model were identified, giving rise to better understanding of the system.
RBNFNs can simultaneously extract both quantitative and qualitative information. However, the insufficiency in both data and biological understanding of gene interactions limit the results obtained from RBNFN models.
We arrive to the conclusion that the presented HRBNF algorithm provides a comprehensive method with the capability of capturing meaningful results. It should be noted that the goal of HRBNF algorithm is not only to provide quantitative predictions, but also to extract knowledge for interactions among genes.
Finally, among these five algorithms, our proposed algorithm has the increased performance in terms of execution time, sensitivity, precision, and F-score. This study found that proposed method is a good strategy to more readily infer gene networks due to its better performance and short execution time.
Despite great progress made with different algorithms, many problems remain unsolved, and much improvement is still required.