 Methodology
 Open Access
 Open Peer Review
 Published:
Identifying timedelayed gene regulatory networks via an evolvable hierarchical recurrent neural network
BioData Mining volume 10, Article number: 29 (2017)
Abstract
Background
The modeling of genetic interactions within a cell is crucial for a basic understanding of physiology and for applied areas such as drug design. Interactions in gene regulatory networks (GRNs) include effects of transcription factors, repressors, small metabolites, and microRNA species. In addition, the effects of regulatory interactions are not always simultaneous, but can occur after a finite time delay, or as a combined outcome of simultaneous and time delayed interactions. Powerful biotechnologies have been rapidly and successfully measuring levels of genetic expression to illuminate different states of biological systems. This has led to an ensuing challenge to improve the identification of specific regulatory mechanisms through regulatory network reconstructions. Solutions to this challenge will ultimately help to spur forward efforts based on the usage of regulatory network reconstructions in systems biology applications.
Methods
We have developed a hierarchical recurrent neural network (HRNN) that identifies timedelayed gene interactions using timecourse data. A customized genetic algorithm (GA) was used to optimize hierarchical connectivity of regulatory genes and a target gene. The proposed design provides a nonfully connected network with the flexibility of using recurrent connections inside the network. These features and the nonlinearity of the HRNN facilitate the process of identifying temporal patterns of a GRN.
Results
Our HRNN method was implemented with the Python language. It was first evaluated on simulated data representing linear and nonlinear timedelayed genegene interaction models across a range of network sizes and variances of noise. We then further demonstrated the capability of our method in reconstructing GRNs of the Saccharomyces cerevisiae synthetic network for in vivo benchmarking of reverseengineering and modeling approaches (IRMA). We compared the performance of our method to TDARACNE, HCCCLINDE, TSNI and ebdbNet across different network sizes and levels of stochastic noise. We found our HRNN method to be superior in terms of accuracy for nonlinear data sets with higher amounts of noise.
Conclusions
The proposed method identifies timedelayed genegene interactions of GRNs. The topologybased advancement of our HRNN worked as expected by more effectively modeling nonlinear data sets. As a nonfully connected network, an added benefit to HRNN was how it helped to find the few genes which regulated the target gene over different time delays.
Background
New opportunities to reverse engineer the activities of different components of complex cellular systems are arising due to technologies like DNA microarrays and RNA sequencing which provide genomicscale data sets [1, 2]. Time series data may be collected either in longitudinal studies of cell or tissue samples collected over multiple time points [3], or expressional change across state space [4]. Effective models of gene regulatory networks (GRNs) have successfully identified regulatory interactions between genes and the specific functional roles of individual genes in cellular systems [5, 6].
Reverse engineering of GRNs occurs within the context of stochastic properties of the system, measurement noise, and high dimensionality [3]. There is strong nonlinearity on temporal patterns of regulatory genes [7]. Further complexity ensues given that genetic interactions among different genes can have different time delays [8, 9]. These delays are due to the transcription and translation of genes varying in composition and length, along with varying kinetics of binding and completion with respect to genes being processed by the transcriptome, the spliceosome and the ribosome. Transcribed and translated products may be further converted and are eventually degraded, with some products being more stable than others. Changing physiological conditions can impact many of the above factors of time delays. As shown in Fig. 1, there are complex combinations by which the expression level of a gene at a certain time could depend upon the expression level of another gene at a previous time point.
Diverse methods with different levels of complexity have been used to model, analyze and infer complex regulatory interactions [10–13]. Boolean networks are the simplest among them [14]. They are based only upon binary outcomes (on and off) for gene expression and therefore lack adequate dynamic resolution. Bayesian networks represent probabilistic relationships among genes and have shown some success in capturing the inherent noise and stochasticity of gene expression data [15]. Dynamic Bayesian Networks (DBN) are an extension of Bayesian networks that can unravel the feedback cycles and loops over time points [16]. However, due to their high computational cost, the application of dynamic Bayesian networks is limited to small networks. Ordinary Differential Equations (ODE) are deterministic models, where interactions among genes represent causal interactions rather than statistical dependencies [17]. They can offer continuous representations of genetic networks, but are not robust for imprecise data.
Methods such as time delay linear regression [18], correlation matrices [19], stochastic simulation algorithms [9, 20], dynamic Bayesian networks [16] and delayed differential equations [21] have been proposed to incorporate a fixed time delay in GRN models. In [22], pairwise correlations between each pair of genes have been used to address the various time delays in gene interactions. TDARACNE (Time DelayAlgorithm for the Reconstruction of Accurate Cellular Networks) has been proposed in [23]. This algorithm detects the timedelayed dependencies between the expression profiles in terms of mutual information by assuming a stationary Markov Random Field as its underlying probabilistic model. The TDARACNE algorithm does not assign any specific delay or regulatory effect on the edges of the GRN. HCCCLINDE [24] is an extension of CLINDE [25], and has been developed to infer a timedelayed GRN in the presence of hidden common causes. All directed pairs of genes in the network have possible delays up to a maximum allowed delay, which is obtained based on either a correlation test or mutual information test.
The main objective of this paper is to reconstruct a timedelayed GRN which takes into account the nonlinearity of gene interaction and the noise of temporal measurements. Recurrent Neural Networks (RNNs) are computational tools inspired by the structural and functional aspects of biological nervous systems, and are noted for their effectiveness in temporal data processing and approximating nonlinear patterns of dynamic temporal behaviors [26]. The ability of RNNs to learn from temporal data, estimate multivariate nonlinear functions, and tolerate noise in measurements makes RNNs an ideal fit for the modeling of gene regulatory interactions using gene expression profiles. Several variants of RNNs have been deployed for the modeling of GRNs including neural fuzzy recurrent networks [27], RNNs combined with particle swarm optimization [28], ensemble of RNNs and support vector machines [29], RNNs combined with differential evolution [30] and RNNs hybridized with the generalized extended Kalman filter [31]. Despite the great capabilities of RNNs for predictive modeling with high accuracy, RNNs are usually considered “black box” models whose internal structure and learned parameters are not interpretable. Due to the multiple layers, the nonlinearity of the model, and cyclic (feedback) connections in the network structure, their interpretability still remains vague [32]. This, in particular, impedes goals with GRN reconstruction to identify pairs of genes, directions of regulation, effects (i.e. up or down regulation), and time delays.
In this paper, we have proposed a hierarchical RNN (HRNN) that surmounts the interpretation difficulties of the RNNs for application of GRN modeling. The proposed design lets us use the features of hierarchical representation in addition to the capabilities of RNNs for finding temporal dependencies. In this way, timedelayed regulations can be captured through hierarchical paths between leaf nodes (regulatory genes) and a target node (regulated gene) in the HRNN. For discovering the underlying hierarchical structure among the regulatory genes and a target gene, the network topology and connection weights are encoded by a customized genetic algorithm (GA). Through the training procedure, in addition to evolving network connection weights, the GA rewires the connectivity and length of hierarchical paths between leaf nodes and the target gene of a population of candidate networks. From the trained HRNN, the direction and effect of gene regulations in the presence of time delays can be captured. Our proposed model is evaluated on a real biological system and linear/nonlinear synthetic generated data for different sizes of networks and variances of noise. The results of our HRNN method are compared with TDARACNE, HCCCLINDE, ODE implemented in the TSNI package [33] and ebdbNet package (Empirical Bayes Dynamic Bayesian Network Inference) [34].
Method
Assume that \(\{G_{1}^{l}(t), G_{2}^{l}(t),\ldots, G_{P}^{l}(t)\}\) are expression levels of P genes at time t in experiment l where t∈{1,…,T ^{l}} and l∈{1,…,L}. The aim is to capture the potential regulators for each gene G _{ i } in a decoupled hierarchical RNN. In this network, G _{ i } is the target gene and the rest of the genes are the potential regulators. At the beginning, a population of candidate hierarchical RNNs should be randomly generated. A candidate network can have 0≤c≤C context nodes in its structure. The network with c context nodes has c+1 neurons. Neurons are the processing units in the RNNs which induce nonlinearity on the inputs. Neurons have multiple inputs and one output. The maximum number of context nodes in the candidate networks is set to C. Context nodes in the hierarchical RNN are nodes without experimental measurements and assist with modeling of temporal dynamics.
Assume that x _{1},…,x _{ C } are context nodes and x _{ C+1},…,x _{ C+P } are genes. In a network with c≤C context nodes, the first c context nodes x _{1},…,x _{ c } and genes (excluded the target gene) are potential inputs of the c+1 neurons in the network. In each candidate network, the target gene is the output of the first neuron, and the context node c _{ i } is the input of n e u r o n _{ i } and output of n e u r o n _{ i+1} where i∈{1,…,c}. In addition to the genes and context node c _{ i }, other context nodes could also be the potential inputs of the neuron n e u r o n _{ i+1}, except for n e u r o n _{ c+1}, which has no context nodes as its inputs. Each input connection has a weight. If context node is the input of the neuron, the corresponding connection weight is positive. Else, it could be positive or negative. Through training, the customized GA evolves the connectivity between nodes and neurons and connection weights.
Figure 2 shows a candidate network generated from a maximum possible number of three context nodes (x _{1},x _{2},x _{3}) and five genes (x _{4},x _{5},x _{6},x _{7},x _{8}). The candidate network in this figure uses two out of the three possible context nodes; thus it has three neurons. This figure shows the regulatory interactions of target gene x _{8} with x _{4},x _{5} and x _{7}. Figure 3(a) shows the corresponding hierarchical recurrent structure obtained from Fig. 2. For Fig. 3(b), timedelayed regulations of the target gene x _{8} are captured by x _{4},x _{5},x _{7} in presence of two context nodes x _{1} and x _{2}.
For the target and context nodes x _{ i } that are outputs of the neurons, Eq. 1 shows the updated value from time t to t+1.
where x _{ j }(t) is the value of j ^{th} input node connected to x _{ i }(t+1). f is a sigmoid function in the form of Eq. 2 which is monotonically increasing in the range of [0,1] and is commonly used in the literature to induce nonlinearity.
In the case of selfregularization of the target gene, Eq. 3 is used for updating the value of the target gene:
where μ is the decay rate of the target gene’s expression over time. Estimating the decay rate for each gene helps to model the suppression effect of a gene on itself [35].
Evolutionary training algorithm
A customized GA is proposed for training the HRNN. At the beginning, the GA generates a population of random candidate networks. The structure and connection weights of the candidate networks are evolved over generations of the GA with the guidance of the fitness function. In each generation, new candidate networks (children) are formed by applying the evolutionary operators (crossover and mutation) on the old networks (parents) within the constraints of the HRNN. Parents are selected according to their fitness values, where networks with higher accuracy have more chances to reproduce. The newly generated population is used for the next generation of the GA. At each generation, an elitist evolving strategy is applied to keep the best candidate networks from the last population. The evolutionary process is repeated until the terminating conditions are satisfied. The proposed procedure is summarized in Algorithm 1.
Representation of candidate networks
Candidate networks in the GA are represented by their number of neurons (N _{ n }), number of inputs to each neuron (N _{ in }), indices of the input nodes (In), weights of the input connections (W) and the decay rate of the target gene’s expression level (μ) if it exists. Components of the first, c ^{th} and last neurons in a candidate network with c+1 neurons are represented in Fig. 4. In this candidate network, P genes and c out of C context nodes are used in the network. One of the genes is considered as a target gene. The output of each neuron (Out) does not change in the training process.
Fitness of candidate networks
The performance of the candidate networks (fitness) is evaluated by measuring the tradeoff between the goodness of fit and complexity of the model by using the Akaike information criterion (AIC) and the Akaike information criterion with correction (AICc). AIC is a model selection criterion which estimates the quality of a model relative to other models (Eq. 4).
where \(x_{i}^{l}(t)\) is the expression value of the target gene i at experiment l, and \(\hat {x}_{i}^{l}\) is the corresponding estimation by the candidate HRNN at time t. k is the number of leaf nodes in the HRNN and n is the total number of temporal samples for gene expression. If n is small or k is large, the AICc is preferred rather than AIC (Eq. 5). As n gets larger, AICc converges to AIC.
Crossover operator
Crossover is an evolutionary operator for generating a new candidate network. Before applying the proposed crossover, the tournament selection reproduces a new pool of candidate networks. The tournament selection is a method of selecting a candidate among a few candidates chosen at random from the population. The winner of each tournament (the one with the best fitness) will be replaced in the new pool.
For applying the proposed crossover, first the networks in the population are shuffled and sorted by the number of neurons (N _{ n }) in their structure. Then, for each pair of selected networks with the same number of neurons (parents 1 and 2), crossover with probability of P _{ c } swaps the random neurons i∈{1,…,N _{ n }} in two parents. Figure 5(c)(d) show the crossover operating on neuron 2 of the parents in Figure 5(a)(b). Crossover creates new candidate networks (crosschildren) with new connectivity and connection weights.
Mutation operator
The mutation operator mutates the number of inputs of the neurons, rewires the connectivity of the inputs of the neurons, and evolves the connection weights with the probability P _{ m }. For a mutation site m _{ site } in the network, the mutation works as below:

If m _{ site } is on the number of inputs of a neuron (N _{ in }), it is mutated to N _{ in }=N _{ in }±1. Therefore, a new input and its corresponding weight are added or deleted.

If m _{ site } is on an input connection of a neuron (In), the selected connection is rewired to another node in the network.

If m _{ site } is on a connection weight of a neuron and input is a context node, the Gaussian mutation evolves the weight in the range of [0,w _{ max }]; else, the weight is mutated in the range of [w _{ min },w _{ max }]

If m _{ site } is on the decay rate μ, the Gaussian mutation is applied to evolve the the decay rate in the range [0,w _{ max }].
Experimental results
In order to evaluate the performance of the proposed method, we have tested our method with both synthetic data and real data against TDARACNE, HCCCLINDE, TSNI and ebdbNet. As the underlying regulatory networks for the real biological datasets are generally unknown, synthetic data are helpful for checking the efficiency of methods. The generated synthetic models in this paper have different levels of complexity and enable us to have a broadranging performance evaluation of our proposed approach in comparison to other approaches. In a real life experiment, we applied our method for finding the GRN of Saccharomyces cerevisiae.
We assess the performance of the inference algorithm on three aspects, namely Links (which is considered correct if and only if both the gene pair and the direction are correct), Delays (which is considered correct if and only if both the link and the time delay are correct), and Effects (which is considered correct if and only if both the link and the sign of an effect are correct) [24]. For each aspect, Recall \(=\frac {TP}{TP+FN}\), Precision \(=\frac {TP}{TP+FP}\) and Fscore \(=\frac {2 \times Precision \times Recall}{Precision+Recall}\) metrics are computed. In these metrics, T P,F P and FN are numbers of true positives, false positives, and false negatives respectively. The Fscores of the results have been compared as an overall measurement of performance. The HCCCLINDE method provides Fscores of Link, Delay and Effect criteria. However, TDARACNE, TSNI and ebdbNet provide information for finding the Fscore of the Link. In all simulations, algorithms have been tested by their default parameters.
Synthetic data
The generation of synthetic data has been considered in two instances of linear and nonlinear models. To compare the accuracy of the HRNN with TDARACNE, HCCCLINDE, TSNI and ebdbNet, the effects of the noise levels and number of genes in small and medium size networks are investigated. For a chosen number of genes and level of noise, ten random GRNs with random connectivity between nodes, weight of the connections, time delay and initial value of gene expressions are generated. The purpose of these experiments is to evaluate the performance of the proposed method in terms of linearity versus nonlinearity of gene expression values, network size (P∈{5,10,20,30}) and noise variance (σ ^{2}∈{0.1,0.25,0.5,1.0,1.5}). The linear and nonlinear models of random GRNs are generated by Eqs. 6 and 7 respectively:
where G _{ i }(t) is the expression value of gene i at time t, i∈{1,…,P} and t∈{1,…,50}. τ _{ ij }∈{1,…,τ _{ max }} is the time delay of the edge j→i (if τ _{ ij }≠0). a _{ ij } is the regulatory effect of gene j on gene i, where the regulatory effect is repressive if a _{ ij } is negative, activatory if positive, and absent if zero. ε _{ i }(t) is a random Gaussian noise with zero mean and variance σ ^{2}. The regulatory effects a _{ ij } are randomly selected at the beginning of each simulation run. The network generation algorithm is set in such a way that each gene could have a maximum of 3 regulators with maximum delay τ _{ max }=4. The number of regulators and time lag for each edge in the synthetic networks are also generated randomly. In Eq. 7, f is a sigmoid and monotonically increasing function in the form of Eq. 2 which adds nonlinearity to the model.
The accuracy of GRN reconstruction using synthetic gene expression data generated from the linear model are presented in Figs. 6, 7 and 8. These figures compare the Link Fscore, Delay Fscore and Effect Fscore of the HRNN with HCCCLINDE, TDARACNE, TSNI and ebdbNet for different number of genes in the network. The variance of noise in these experiments is the same, and is equal to (σ ^{2}=1). In part a of these figures, the box plot of Fscore values of Link, Delay and Effect criterion are compared. In parts b of these figures, a linear regression model is fit on the Fscore values with respect to the number of genes in the network. Rsquared and pvalues of the linear regression models in Figs. 6(b), 7(b) and 8(b) are stated in Table 1. In case of the linear data, the linear regression models in Fig. 6(b) shows that the HRNN and the HCCCLINDE are responsibly competitive for finding the correct Link between the nodes in networks.
Table 2 includes the average of the TP, FP, FN, Precision, Recall and Fscore of Link criteria for 10 independent runs of the methods and different number of genes in the network. Also, we conduct a hypothesis test for the difference between means of Fscore. For a selected number of genes in Table 2, a ttest performed on Fscore values of the HRNN and other methods. The null hypothesis is defined as two population means are equal. The nominal and adjusted pvalues are mentioned in the table. HRNN is tested multiple times for four different number of genes in the network, to obtain the adjusted pvalues, the nominal pvalues are multiplied by four. If the corresponding adjusted pvalue is less than 0.05, the null hypothesis is rejected, meaning that mean of Fscores are significantly different. pvalues are for how the Fscores of the other methods may significantly differ with the Fscores of HRNN; therefore, pvalues are not shown for the HRNN row in the tables. Table 3 includes the average of the TP, FP, FN, Precision, Recall, Fscore and also pvalues of Delay and Effect criterion for 10 independent runs of the HRNN and HCCCLINDE and different number of genes in the network. Results shows that HRNN and HCCCLINDE are not significantly different in terms of Delay and Effect Fscores in case of linearity among genes.
In the next step for testing synthetic data, we considered a more realistic scenario where the gene expression values are generated from nonlinear models. For 10 different randomly generated networks, the effect of the number of genes in the accuracy of GRN reconstruction of HRNN was compared with other methods. Figures 9, 10 and 11 compare the accuracy of Link Fscore, Delay Fscore and Effect Fscore for different number of genes in the network respectively. The variance of noise in these experiments is the same, and is equal to (σ ^{2}=1). The linear regression models between Fscore and number of genes, shown in part b of these figure. Rsquared and pvalues of the linear regression models in Figs. 9(b), 10(b) and 11(b) are stated in Table 4.
Tables 5 and 6 include the average of the TP, FP, FN, Precision, Recall and Fscore in term of Link, Delay and Effect criterion for 10 independent runs of the methods and different number of genes in the network. For a selected number of genes in these tables, a ttest performed on Fscore values. If the corresponding adjusted pvalue is less than 0.05, the null hypothesis is rejected, meaning that mean of Fscore of HRNN is significantly different from other method. The results show that the proposed HRNN works better than HCCCLINDE, TDARACNE, TSNI and ebdbNet for cases of nonlinearity among gene interactions in bigger networks. The accuracy of HCCCLINDE for finding the correct link drops significantly in case of nonlinearity between the genes in comparison to the the linear relationships.
In another simulation, the effect of noise in nonlinear models is examined. In synthetic nonlinear networks with 20 genes, the level of noise is changed (σ ^{2}∈{0.1,0.25,0.5,1.0,1.5}). Figures 12, 13 and 14 show the results of GRN reconstruction of HRNN in comparison to TDARACNE, HCCCLINDE, TSNI and ebdbNet in the forms of box plot and linear regression model and in terms of Link, Delay and Effect Fscores respectively. Rsquared and pvalues of the linear regression models in Figs. 12(b), 13(b) and 14(b) are stated in Table 7. Tables 8 and 9 include the average of the TP, FP, FN, Precision, Recall, Fscore, nominal and adjusted pvalue for 10 independent runs of the methods and different levels of the noise. HRNN is tested multiple times for five different noise levels of gene measurements, to obtain the adjusted pvalues, the nominal pvalues are multiplied by five. Figure 12 show that the accuracy of our proposed method in terms of Link Fscore is often higher than HCCCLINDE, TDARACNE, TSNI and ebdbNet. Figures 13 and 14 show that the HRNN and HCCCLINDE are competitive in terms of Delay and Effect Fscores.
Reallife biological data of Saccharomyces cerevisiae (IRMA)
In order to validate the performance of the proposed method on reallife biological GRNs, we considered a recent significant contribution to systems biology reported in [36] where the authors built a synthetic network, called IRMA, of the yeast organism Saccharomyces cerevisiae. The researchers tested transcription of network genes by culturing cells in presence of galactose and glucose. This is one of the first attempts at building a reference data set, having a fairly true table of regulations [8, 23]. The regulatory network includes five genes. It is negligibly affected by endogenous genes. Two sets of gene profiles called Switch ON and Switch OFF were provided, each containing 16 and 21 time points. The former corresponds to the shifting of the growing cells from glucose to the galactose medium; the latter corresponds to the reverse shift. Due to the lack of stimulus, reconstruction of the GRN from the Switch OFF dataset is difficult [8, 23].
The performance comparisons among various methods for the IRMA ON dataset are shown in Fig. 15 and Table 10. In Fig. 15(a), the true IRMA network is shown. Figure 15(b) displays the GRN obtained by the proposed method. Figure 15(cf) present the GRN reconstructions by TDARACNE, HCCCLINDE, TSNI and ebdbNet obtained with default parameters. In Table 10, TP, FP, FN, precision, recall and Fscore values are also compared. The proposed HRNN finds the true regulations of ASH1 by SWI5, CBF1 by ASH1, GAL80 by SWI5, GAL4 by GAL80, GAL80 by GAL4 and CBF1 by SWI5. The regulation of SWI5 by GAL4 and regulation of GAL4 by CBF1 are not found. The method finds three regulations (regulations of SWI5 by CBF1, SWI5 by ASH1 and SWI5 by GAL4) that are not in the true network of IRMA. Among the eight connections in the true network, TDARACNE finds two correct regulations. The HCCCLINDE method found one true regulation and one false regulation. Also, HCCCLINDE finds the regulations of ASH1 and SWI5 by a hidden common cause which is not reported in the actual GRN of the IRMA. Results show higher accuracy in the proposed HRNN approach for finding the regulatory interactions between the genes in comparison to the other two approaches.
Conclusions
In this study, we developed and implemented a hierarchical recurrent neural network (HRNN) approach to identify timedelayed regulatory interactions of genes. The designed HRNN facilitates capturing the paths with different lengths from the leaf nodes in the network to the target node. Hierarchy and nonlinearity in the network and the allowance for recurrent connections in HRNN provide an effective capability for modeling the temporal patterns of gene expression. Furthermore, partial connectivity of the nodes aids in finding the limited set of genes which regulate the target gene over different time delays. The proposed method outperformed TDARACNE, HCCCLINDE, TSNI and ebdbNet in terms of reconstructing small and medium size networks having nonlinearity and high levels of noise for measurement data.
Abbreviations
 GA:

Genetic algorithm
 GRN:

Gene regulatory network
 HRNN:

Hierarchical recurrent neural network
References
 1
Sima C, Hua J, Jung S. Inference of gene regulatory networks using timeseries data: a survey. Current Genomics. 2009; 10(6):416–29. doi:10.2174/138920209789177610.
 2
Wang Z, Yang F, Ho DWC, Swift S, Tucker a, Liu X. Stochastic dynamic modeling of short gene expression timeseries data. IEEE Trans Nanobiosci. 2008; 7(1):44–55. doi:10.1109/TNB.2008.2000149.
 3
BarJoseph Z, Gitter A, Simon I. Studying and modelling dynamic biological processes using timeseries gene expression data. Nat Rev Genet. 2012; 13(8):552–64.
 4
Mar JC, Quackenbush J. Decomposition of gene expression state space trajectories. PLoS Comput Biol. 2009; 5(12):1–10. doi:10.1371/journal.pcbi.1000626.
 5
Bansal M, Belcastro V, AmbesiImpiombato A, di Bernardo D. How to infer gene networks from expression profiles. Mol Syst Biol. 2007; 3:78. doi:10.1038/msb4100120.
 6
Chai LE, Loh SK, Low ST, Mohamad MS, Deris S, Zakaria Z. A review on the computational approaches for gene regulatory network construction. Comput Biol Med. 2014; 48:55–65. doi:10.1016/j.compbiomed.2014.02.011.
 7
Hasty J, McMillen D, Isaacs F, Collins JJ. Computational studies of gene regulatory networks: in numero molecular biology. Nat Rev Genet. 2001; 2(4):268–79.
 8
Morshed N, Chetty M, Xuan Vinh N. Simultaneous learning of instantaneous and timedelayed genetic interactions using novel information theoretic scoring technique. BMC Syst Biol. 2012; 6(1):62. doi:10.1186/17520509662.
 9
Bratsun D, Volfson D, Tsimring LS, Hasty J. Delayinduced stochastic oscillations in gene regulation. Proc Natl Acad Sci. 2005; 102(41):14593–8. doi:10.1073/pnas.0503858102.
 10
Schlitt T, Brazma A. Current approaches to gene regulatory network modelling. BMC Bioinforma. 2007; 8(6):9. doi:10.1186/147121058S6S9.
 11
Karlebach G, Shamir R. Modelling and analysis of gene regulatory networks. Nat Rev Mol Cell Biol. 2008; 9(10):770–80.
 12
Lee WP, Tzou WS. Computational methods for discovering gene networks from expression data. Brief Bioinform. 2009; 10(4):408–23. doi:10.1093/bib/bbp028.
 13
Ristevski B. A survey of models for inference of gene regulatory networks. Nonlinear Anal Modell Control. 2013; 18(4):444–65.
 14
Bornholdt S. Boolean network models of cellular regulation: prospects and limitations. J R Soc Interface. 2008; 5 Suppl 1:85–94. doi:10.1098/rsif.2008.0132.focus.
 15
Friedman N, Linial M, Nachman I, Pe’er D. Using Bayesian networks to analyze expression data. J Comput Biol; 7(34):601–20.
 16
Zou M, Conzen SD. A new dynamic Bayesian network (DBN) approach for identifying gene regulatory networks from time course microarray data. Bioinformatics. 2005; 21(1):71–9. doi:10.1093/bioinformatics/bth463.
 17
Sakamoto E, Iba H. Inferring a system of differential equations for a gene regulatory network by using genetic programming. In: Proceedings of the 2001 Congress on Evolutionary Computation, vol. 1.Seoul: IEEE: 2001. p. 720–261. doi:10.1109/CEC.2001.934462.
 18
Huang T, Liu L, Qian Z, Tu K, Li Y, Xie L. Using GeneReg to construct time delay gene regulatory networks. BMC Res Notes. 2010; 3(1):142. doi:10.1186/175605003142.
 19
van Someren EP, Vaes BLT, Steegenga WT, Sijbers AM, Dechering KJ, Reinders MJT. Least absolute regression network analysis of the murine osteoblast differentiation network. Bioinformatics. 2006; 22(4):477–84. doi:10.1093/bioinformatics/bti816.
 20
Barrio M, Burrage K, Leier A, Tian T. Oscillatory regulation of Hes1: discrete stochastic delay modelling and simulation. PLoS Comput Biol. 2006; 2(9):117. doi:10.1371/journal.pcbi.0020117.
 21
Kim S, Kim J, Cho KH. Inferring gene regulatory networks from temporal expression profiles under timedelay and noise. Comput Biol Chem. 2007; 31(4):239–45. doi:10.1016/j.compbiolchem.2007.03.013.
 22
ElBakry O, Ahmad MO, Swamy MNS. Inference of gene regulatory networks with variable time delay from timeseries microarray data. IEEE/ACM Trans Comput Biol Bioinforma. 2013; 10(3):671–87. doi:10.1109/TCBB.2013.73.
 23
Zoppoli P, Morganella S, Ceccarelli M. TimeDelayARACNE: reverse engineering of gene networks from timecourse data by an information theoretic approach. BMC Bioinforma. 2010; 11(1):154. doi:10.1186/1471210511154.
 24
Lo LY, Wong ML, Lee KH, Leung KS. Time delayed causal gene regulatory network inference with hidden common causes. PLoS ONE. 2015; 10(9):1–47. doi:10.1371/journal.pone.0138596.
 25
Lo LY, Leung KS, Lee KH. Inferring timedelayed causal gene network using timeseries expression data. IEEE/ACM Trans Comput Biol Bioinformatics. 2015; 12(5):1169–82. doi:10.1109/TCBB.2015.2394442.
 26
Haykin SS. Neural Networks and Learning Machines. New Jersey: Pearson; 2009.
 27
Maraziotis IA, Dragomir A, Bezerianos A. Gene networks reconstruction and timeseries prediction from microarray data using recurrent neural fuzzy networks. IET Syst Biol. 2007; 1(1):41–50. doi:10.1049/ietsyb:20050107.
 28
Xu R, Venayagamoorthy GK, Wunsch DC. Modeling of gene regulatory networks with hybrid differential evolution and particle swarm optimization. Neural Netw. 2007; 20(8):917–27. doi:10.1016/j.neunet.2007.07.002.
 29
Ao SI, Palade V. Ensemble of Elman neural networks and support vector machines for reverse engineering of gene regulatory networks. Appl Soft Comput. 2011; 11(2):1718–26. doi:10.1016/j.asoc.2010.05.014.
 30
Noman N, Palafox L, Iba H. Inferring genetic networks with a recurrent neural network model using differential evolution. Berlin: Springer; 2014, pp. 355–73.
 31
Raza K, Alam M. Recurrent neural network based hybrid model for reconstructing gene regulatory network. Comput Biol Chem. 2016; 64:322–34. doi:10.1016/j.compbiolchem.2016.08.002.
 32
Castelvecchi D. Can we open the black box of AI?Nature. 2016; 538(7623):20.
 33
Bansal M, Gatta GD, di Bernardo D. Inference of gene regulatory networks and compound mode of action from time course gene expression profiles. Bioinformatics. 2006; 22(7):815. doi:10.1093/bioinformatics/btl003.
 34
Rau A, Jaffrézic F, Foulley JL, Doerge RW, et al.An empirical bayesian method for estimating biological networks from temporal microarray data. Stat Appl Genet Mol Biol. 2010; 9(1):1544–6115.
 35
Zhang Y, Xuan J, de los Reyes BG, Clarke R, Ressom HW. Reverse engineering module networks by PSORNN hybrid modeling. BMC Genomics. 2009; 10(1):15. doi:10.1186/1471216410S1S15.
 36
Cantone I, Marucci L, Iorio F, Ricci MA, Belcastro V, Bansal M, Santini S, di Bernardo M, di Bernardo D, Cosma MP. A yeast synthetic network for in vivo assessment of reverseengineering and modeling approaches. Cell. 2009; 137(1):172–81. doi:10.1016/j.cell.2009.01.055.
Acknowledgements
This work is partially supported by the Expeditions in Computing by the National Science Foundation under Award Number: CCF1029731, and the Air Force Research Laboratory and OSD under agreement number FA87501520116. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes, notwithstanding any copyright notation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of NSF, Air Force Research Laboratory, OSD or the U.S. Government. The first, second and the forth authors would like to acknowledge the support from the NSF, and the fourth author would also like to acknowledge the support from the Air Force and the OSD.
Funding
This study was funded by the Expeditions in Computing by the National Science Foundation under award number CCF1029731.
Availability of data and materials
The Python source code providing the HRNN algorithm (files ‘mainGRN.ipynb’ and ’UtilsGRN.py’) and the datasets supporting the conclusions of this article are available in the Github repository [https://github.com/mmoradik/GRN].
Author information
Affiliations
Contributions
MMK contributed to the design of the study and experiments, invention of the HRNN algorithm, data analysis, and writing of the paper; MMK and MGS prepared samples and performed the code implementation; SHH and AH helped design the study and edited the paper. All authors reviewed the paper. All authors read and approved the final manuscript.
Corresponding author
Correspondence to Abdollah Homaifar.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Gene regulatory network
 Hierarchical recurrent neural network
 Genetic algorithm
 Time delay