Gene Regulatory Identification Based on the Novel Hybrid Time-Delayed Method.

Gene regulatory network (GRN) inference with biology data is a difficult and serious issue in the field of system biology. In order to detect the direct associations of GRN more accurately, a novel two-step GRN inference technique based on the time-delayed correlation coefficient (TDCC) and time-delayed complex-valued S-system model (TDCVSS) is proposed. First, a TDCC algorithm is utilized to construct an initial network. Second, a TDCVSS model is utilized to prune the network topology in order to delete false-positive regulatory relationships for each target gene. The complex-valued restricted additive tree and complex-valued differential evolution are proposed to approximate the optimal TDCVSS model. Finally, the overall network could be inferred by integrating the regulations of all target genes. Two real gene expression datasets from E. coli and S. cerevisiae gene networks are utilized to evaluate the performances of our proposed two-step GRN inference algorithm. The results demonstrated that the proposed algorithm could infer GRN more correct than classical methods and time-delayed methods.

Introduction

With the accomplishment of the human genome project and the emergence of high-throughput gene analysis techniques, it has been recognized that great heterogeneity in gene mutation profiles of cancer tissues has been reported (Thomas et al., 2007; Tesniere et al., 2010; Rupaimoole et al., 2016; Zhang et al., 2018). Mutant genes are not only numerous and complex but could also construct a dynamic network system. At the cell level, interrelated genes/proteins constitute complex cellular networks, including signaling pathways, gene regulatory networks (GRN), and metabolic networks (Quach et al., 2007; Ma’ayan, 2009). When it comes to the GRN issue, any gene could not execute the independent function. They must coordinate with other genes to complete certain biological processes and participate in determining the behaviors and phenotypes of cells (Hernández- Prieto et al., 2014). Thus, the research on GRN has become a hotspot in the field of system biology during the past several decades (Bracken et al., 2016; Liu et al., 2018).

The gene regulatory network contains two parts, including nodes and edges. The main work is to identify the direct regulations of all pairs of nodes. In real organisms, genetic and non-genetic factors could cause an objective phenomenon that time delays occur in the gene expression process. Recently, several machine learning methods have been utilized to identify time-delayed regulations among genes (Parmar et al., 2015; Wang et al., 2020). Lo proposed a causal network model by identifying causal-directed regulations with time delays (Lo et al., 2015). Kordmahalleh et al. presented a hierarchical recurrent neural network (HRNN) and genetic algorithm (GA) to infer the time-delayed gene regulatory network (Kordmahalleh et al., 2017). Li et al. utilized the Max–Min high-order dynamic Bayesian network (MMHO-DBN) to identify the synthetic and real time–delayed gene regulatory network (TDGRN) (Hu et al., 2020). Yu and Li also utilized the dynamic Bayesian network (DBN) to infer the TDGEN (Li et al., 2014; Yu et al., 2017). Some time-delayed information theory methods have been proposed to infer the TDGEN, such as TimeDelay-ARACNE (Zoppoli et al., 2010) and the time-lagged correlation method (Sefidmazgi et al., 2016; Abduallah and Wang, 2017). Zhao et al. proposed a kind of delay differential equation model to infer the gene regulatory networks (Zhao et al., 2016).

The S-system model is a nonlinear ordinary differential equation (ODE) composed of power-law functions, which is very suitable for simulating nonlinear biological systems with a large number of components, such as GRNs and metabolic pathways (Miyawaki-Kuwakado et al., 2020). Liu et al. proposed a multiobjective optimization method to evolve the S-system in order to infer the biochemical network (Liu and Wang, 2008). Orland et al. utilized simulated annealing (SA) to search the optimal parameters of S-system in order to simulate real biochemical networks (Gonzalez et al., 2007). Wang et al. proposed a simplified S-system and a multi-dimensional optimization method for GRN inference (Wang et al., 2010). Iwata et al. utilized the S-system to simulate metabolic reaction systems (Iwata et al., 2014). Chowdhury et al. proposed the time-delayed S-system and stochastic S-system to model time-delayed and stochastic regulations in GRN, respectively, and gained a good performance (Chowdhury et al., 2013a; Chowdhury et al., 2013b; Ji et al., 2017).

With several decades of efforts, complex-valued methods have been proposed to solve the real prediction and classification issues. Compared with real-valued methods, complex-valued methods have stronger modeling and noise tolerance abilities (Yang and Bao, 2019; Yuan et al., 2021). Fink et al. proposed a complex-valued multilayer feedforward neural network to forecast the degradation of railway track turnouts (Fink et al., 2014). Chen et al. utilized a complex-valued radial basis function network to solve a nonlinear signal processing problem (Chen et al., 1994). Goh et al. utilized complex-valued recurrent neural networks (RNNs) to predict Santa Fe and chaotic Mackey–Glass time series data (Goh et al., 2006). Savitha et al. proposed a complex-valued version of the extreme learning machine to solve real-valued classification problems (Savitha et al., 2012). Rashid utilized the complex-valued neural network to solve classification problems in the bioinformatics field (Rashid et al., 2016). Bakbak et al. presented the complex wavelet neural network to classify the sonar signal (Bakbak and Peker, 2020).

In order to enhance the accuracy of GRN inference, this study presents a time-delayed complex-valued S-system model (TDCVSS) to identify time-delayed and nonlinear relationships among genes. Compared with the S-system, TDCVSS contains time-delayed and complex-valued parameters, and the variables are complex-valued. The time-delayed correlation coefficient (TDCC) algorithm is first utilized to construct a TDCC matrix and the optimal time delay vector between genes. According to the TDCC matrix, the initial network is constructed. A complex-valued hybrid swarm intelligent algorithm based on the restricted additive tree and differential evolution is utilized to search for the optimal TDCVSS model in order to prune the network topology further.

Methods

Time-Delayed Correlation Coefficient

The time-delayed correlation coefficient (TDCC) is utilized to evaluate the linear relationship between two genes under the conditions of a constant system time delay, which is described as follows.

where and are gene expression profiles of gene and gene at the time points and , respectively. is the number of sample points, and and are the means of gene expression levels of gene and gene , respectively.

Time-Delayed Complex-Valued Time-Delayed S-System

The time-delayed complex-valued S-system (TDCVSS) is the time-delayed and complex-valued version of the S-system. Compared with a real-valued S-system, TDCVSS contains two improvements. Input variables and rate constants ( and ) are complex-valued. Also, in a TDCVSS model, a time-delayed factor is included. The time-delayed and complex-valued differential equation is given in Eq. 2.

where and are real-valued kinetic orders. is the time delay between variable and variable .

Complex-Valued Restricted Additive Tree

The TDCVSS model contains complex-valued variables and coefficients. In the GRN, each target gene corresponds to a small number of regulatory factors. Thus, for each dependent variable in TDCVSS, the proper independent variables need to be selected. A complex-valued restricted additive tree algorithm (CVRAT) is utilized to evolve the structure of the model. An example of the chromosome of TDCVSS can be demonstrated in Figure 1. The node in the first layer is fixed to subtraction (-). Two operator sets ( and ) were utilized to create the nodes randomly in other layers. is the product of complex-valued input variables. In order to represent the parameters of the TDCVSS model, a real-valued parameter ( or ) is given to each variable node and a complex-valued parameter ( or ) is given to each branch of the root node. The time delay vector between variable and other variables is given as . The corresponding TDCVSS model is .

FIGURE 1

Example of the chromosome of CVRAT.

In the CVRAT algorithm, three genetic operators (selection, crossover, and mutation) are used to evolve the chromosome populations, which are the same as some structure-based evolutionary algorithms, such as genetic programming (GP). The detail crossover and mutation operators are shown in Figures 2, 3, respectively (Yang et al., 2020).

FIGURE 2

Crossover operator.

FIGURE 3

Mutation operators.

Complex-Valued Differential Evolution

Differential evolution (DE) is an efficient and global evolutionary algorithm, which is based on the continuous variable optimization (Das and Suganthan, 2011). Its idea comes from the genetic algorithm (GA), which also contains the crossover, mutation, and reproduction. However, the mutation vector of differential evolution is generated by the difference vector of the parent generation, and the new individual could be generated by the crossover of the parent individuals. Considering its simple structure, easy implementation, and strong robustness, DE is widely utilized in many fields, such as bioinformatics, image processing, document extraction, artificial neural network, and electromagnetics. Complex-valued differential evolution (CVDE) is the complex-valued version of DE. In CVDE, a complex-valued individual includes the real part and the imaginary part, which need to be evolved simultaneously. CVDE could improve the diversity of population and premature convergence of DE. The optimization process of parameters of TDCVSS with CVDE is introduced in Algorithm 1.

Parameter optimization of TDCVSS with complex-valued differential evolution.

Gene Regulatory Network Inference Algorithm

The proposed network inference algorithm contains two steps, whose flowchart is demonstrated in Figure 4.

FIGURE 4

Flowchart of gene regulatory network inference.

Construction of the Initial Gene Regulatory Network

1) According to the given gene expression data, the TDCC method is utilized to obtain the TDCC matrix and the optimal time delay vector between regulatory factors and target genes. The maximum time delay is set, and the CC values are calculated among genes with the different time delays. The maximum CC value and time delay are selected as TDCC and the optimal time-delayed value between the regulatory factor and target gene separately.

2) According to the TDCC matrix, a directed weighted graph is obtained. The initial network is constructed with the selected threshold. If the weight of an edge is less than the threshold, the edge is deleted. If the weight of an edge is higher than the threshold, the edge is retained.

Pruning of the Gene Regulatory Network

In order to improve the false-positive rate, the TDCVSS model is utilized to select the regulatory factors of each target gene. According to the initial network, the optimal time factor vector, and gene expression data, the optimal TDCVSS model of each target gene is found by CVRAT and CVDE. If a regulatory factor is not included in the TDCVSS, this regulatory factor could not regulate the target gene. With such an approach, some false-positive regulatory relationships could be deleted in order to prune the network.

Experiments

In this part, two real gene regulatory networks from E. coli and Saccharomyces cerevisiae are utilized. The true-positive rate (TPR), false-positive rate (FPR), positive predictive value (PPV), accuracy (ACC), and F-score are utilized to evaluate the performance of our method. In order to test our method well, the TDCVSS model is utilized to infer two real GRNs without the TDCC algorithm. Some classical GRN inference methods such as DBN (Perrin et al., 2003) [MMHO-DBN (Hu et al., 2020), DBN-ZC (Zou and Conzen, 2005) and DBmcmc (Husmeier, 2003)], RNN (Xu et al., 2007; Kordmahalleh et al., 2017), ODE (Chen et al., 2011), time-delayed methods [TDARACNE (Zoppoli et al., 2010), and TDLASSO (Mundra et al., 2013)] are also utilized.

SOS Repair Network

The first real biological gene expression data were derived from the SOS (Save Our Souls) DNA repair system. SOS DNA repair is a kind of DNA repair method induced by the serious damage of DNA and the cell in a crisis state under the action of a variety of enzymes, in order to maintain the integrity of the genome. The SOS reaction in DNA of E. coli is controlled by , and the network structure is depicted in Figure 5, which contains eight genes: uvrD, lexA, umuD, recA, uvrA, uvrY, ruvA, and polB (Ronen et al., 2002).

FIGURE 5

SOS DNA repair network of E. coli.

According to the gene expression levels, the time-delayed correlation coefficient is utilized to construct the network topology, which contains 29 regulatory relationships. Through our proposed evolutionary algorithm, the corresponding optimal CVTDSS models of eight target genes have been obtained in order to reduce the false-positive regulations and prune the network structures, which are described in Eq. 3. The final SOS repair network is obtained in Figure 6B, according to our proposed method. Figure 6A is the network obtained by CVTDSS without TDCC. The solid lines represent true-positive regulations, while dotted lines show the false-positive relationships. Comparing Figures 6A,B, it could be seen clearly that our proposed method infer less false-positive regulations.

FIGURE 6

Network by CVTDSS (A) and the network by TDCC + CVTDSS (B).

DBN, RNN, S-system, and ODE are also utilized to infer the SOS repair network, and the results are listed in Table 1. From Table, we can figure out that in terms of TPR, our method and CVSS have the same performance, which is 59.99% higher than the S-system and RNN, 100% higher than DBN, and 33.27% higher than ODEs, which reveal that our method could infer more real regulations. In terms of FPR, RNN has the best performance. Our method only performs better than S-system, ODE, and CVSS, which show that the network obtained by our method does not have the least false-positive edges, but it has fewer wrong regulations. In terms of PPV, RNN has the highest performance, while our method has the second best performance. In terms of ACC, our method is 31.2% higher than the S-system, 16% higher than DBN, 8.4% higher than ODE, and 33.3% higher than CVSS. The results of F-score show that our method performs best as a whole.

TABLE 1

Performance comparison of six methods for the SOS network.

	TPR	FPR	PPV	ACC	F-Score
S-system	0.5556	0.20833	0.3333	0.6667	0.41667
DBN	0.4444	0.10417	0.44444	0.75439	0.44444
RNN	0.55556	0.041667	0.71429	0.80702	0.625
ODEs	0.6667	0.3125	0.28571	0.57895	0.4
TDCVSS	0.88889	0.3818	0.27586	0.65625	0.42105
Our method	0.88889	0.12723	0.5333	0.875	0.6667

IRMA Network

The second real gene expression data are from the IRMA network, which is extracted from the switch process of galactose and glucose in Saccharomyces cerevisiae. According to the on and off of galactose creation, two kinds of gene expression datasets (on dataset and off dataset) are collected (Cantone et al., 2009). The real IRMA network is depicted in Figure 7.

FIGURE 7

IRMA network.

With the on dataset, the TDCC could select 13 potential regulations for the initial IRMA network. According to the initial network and on dataset, five optimal CVTDSS models are found for five target genes (Eq. 4), which are utilized to determine the regulations of target genes and construct the final network (Figure 8B). The network inferred only by CVTDSS with the on dataset is depicted in Figure 8A. The solid lines represent true-positive regulations, while dotted lines show the false-positive relationships. Comparing two networks, it could be seen that our method could infer six real regulations, while CVSS can gain seven real relationships. But our method could infer less false-positive regulators and obtain a more accurate network.

FIGURE 8

IRMA network by CVTDSS (A) and the IRMA network by TDCC + CVTDSS (B) with the on dataset.

HRNN, MMHO-DBN, TDARACNE, TDLASSO, DBmcmc, and DBN-ZC are also utilized to infer the IRMA network with the on dataset. Also, the results are described in Table 2. In terms of TPR, the CVSS model obtains the best performance, which proves that CVSS could infer more true-positive edges. Our method infers one edge less than CVSS, which may be because the TDCC method deletes this edge in the previous step. In terms of FPR, MMHO-DBN could gain zero, which reveals that the network inferred has no false-positive relationships. Our method could infer one false-positive edge and gain the second best FPR performance, which is 0.05882. In terms of ACC, our method could gain the highest accuracy, which shows that our method could infer more true-false edges and less false-positive edges. As a whole, our method has the highest F-score values. From the table, we could see that our method performs better than classical methods (HRNN, MMHO-DBN, DBmcmc, and DBN-ZC) and time-delayed methods (TDCVSS, TDARACNE, and TDLASSO).

TABLE 2

Performance comparison of eight methods for IRMA network inference with the on dataset.

Method	TPR	FPR	PPV	ACC	F-Score
Our method	0.75	0.05882	0.857,143	0.88	0.8
TDCVSS	0.875	0.411,765	0.5	0.68	0.636,364
HRNN	0.75	0.176,471	0.667	0.8	0.706,069
MMHO-DBN	0.5	0	1	0.84	0.666,667
TDARACNE	0.625	0.117,647	0.7142	0.8	0.666,629
TDLASSO	0.25	0.176,471	0.4	0.64	0.307,692
DBmcmc	0.25	0.117,647	0.5	0.68	0.333,333
DBN-ZC	0.375	0.117,647	0.6	0.72	0.461,538

With the off dataset, the TDCC could select 15 potential regulations for the initial IRMA network. According to the initial network and off dataset, five optimal CVTDSS models are found for five target genes (Eq. 5), which are utilized to determine the regulations of target genes and construct the final network (Figure 9B). The network inferred only by TDCVSS with the off dataset is depicted in Figure 9A. The solid lines represent true-positive regulations, while dotted lines show the false-positive relationships. Compared with two networks, our method and TDCVSS could infer the same number of true-positive edges, which is six. But our method could infer less false-positive regulators.

FIGURE 9

IRMA network by CVTDSS (A) and the IRMA network by TDCC + CVTDSS (B) with the off dataset.

With the off dataset MMHO-DBN, TDARACNE, TDLASSO, and DBmcmc are also utilized to infer the IRMA network. The results are listed in Table 3. In terms of TPR, our method and the TDCVSS model can obtain the best performance, which is 0.75. Compared with other methods, these two methods could infer more true-positive edges. In terms of FPR, MMHO-DBN and TDARACNE have better performance, which reveals that the networks inferred by the two methods have fewer false-positive relationships than other methods. Our method and TDLASSO could gain the second best FPR performance, which is 0.176,471. TDCVSS has the worst FPR value. In terms of PPV, our method and MMHO-DBN obtain the best performance, which is 77.79% higher than TDCVSS, 33.34% higher than TDARACNE, 166.7% higher than TDLASSO, and 292.2% higher than DBmcmc. In terms of ACC and F-score, our method could gain the best performance, which shows that our method could infer more true-false edges and fewer false-positive edges than other methods as a whole.

TABLE 3

Performance comparison of six methods for IRMA network inference with the off dataset.

Method	TPR	FPR	PPV	ACC	F-Score
Our method	0.75	0.176,471	0.6667	0.8	0.705,901
TDCVSS	0.75	0.588,235	0.375	0.52	0.5
MMHO-DBN	0.25	0.058824	0.6667	0.72	0.363,641
TDARACNE	0.125	0.058824	0.5	0.68	0.2
TDLASSO	0.125	0.176,471	0.25	0.6	0.166,667
DBmcmc	0.12	0.294,118	0.17	0.52	0.14069

Conclusion

In order to improve the accuracy of time-delayed GRN inference, a novel GRN inference method is proposed. In our method, the time-delayed correlation coefficient algorithm is first utilized to construct the TDCC matrix and the optimal time delay vector between genes. According to the TDCC matrix, the initial gene regulatory network topology is inferred. In order to improve the false-positive rate of GRN obtained, the time-delayed complex-valued S-system model is proposed to identify the regulations of each target gene, which could delete many false-positive relationships. When tested on two real gene expression datasets from E. coli and S. cerevisiae gene networks, in terms of F1, our method could make the 13.3–406.5% improvements, which show that our method as a whole has better performances than DBN (MMHO-DBN, DBN-ZC, and DBmcmc), RNN, ODE, and time-delayed methods (TDARACNE and TDLASSO).

From the results, it could be seen that our method could infer more true-positive regulations and fewer false-positive relationships than other classical GRN inference methods. However, each gene needs to identify the optimal CVTDSS model through an evolutionary algorithm, so the runtime of the proposed algorithm may be high. In future research, the parallel computing framework is planned to improve the time efficiency of the algorithm.