A non-linear reverse-engineering method for inferring genetic regulatory networks

Experimental data

In this work, we used the sub-series GSE49987 as the experimental data from the published microarray dataset GSE49991 (May et al., 2013). This dataset contains the expression profiles collected by experiments using the cell line FDCPmix. This dataset was generated with the probe name version of Agilent Whole Mouse Genome Microarray 4 × 44 K (May et al., 2013). It provides microarray gene expression profiles of hematopoietic stem cells (HSCs) differentiating into erythrocytes and neutrophils. This microarray dataset is available at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE49991. To convert all microarray probe IDs to gene names, we pre-processed this dataset based on the Ensembl BioMart and GO Enrichment Analysis (The Gene Ontology Consortium, 2017). From a previous study, the regulatory network of 18 core genes during the hematopoiesis has been curated (Moignard et al., 2013). Moreover, the same research team studied the regulatory interaction of 26 core genes during the hematopoiesis (Moignard et al., 2015). The total number of distinct genes in these two studies is 30. Thus, in our work we considered 30 genes whose names are listed in Table S1. There are three repeated experiments for each developmental process, each of which contains the expression levels of 30 genes from HSCs to differentiated cells at 30 time points spanning over 1 week. The observation time points are those starting from the HSCs/progenitors stage (1 point), then every 2 h over the first day (12 points), every 3 h over the second day (8 points), every 4 h over the third day (6 points), every 24 h until the fifth day (2 points), and the seventh day (1 point). In this study, we used the average data of these three repeated tests as the experimental data for each time point. The time points and expression data of four genes can be found in the following Figs. 1 and 2.

Selection of candidate genes

Based on our research experience (Wang et al., 2016), it is challenging to study a dynamic network with 30 genes. Thus, we conducted an extensive literature review for selecting a smaller number of important genes based on their relationship with the three genes Gata1, Gata2 and PU.1. These candidate genes should be essential for the cell-fate choice in hematopoiesis, or they significantly interact with these three genes. For example, gene Scl/Tal1 interacts with Gata1, Eto2/Cbfa2t3 and Ldb1 (Goardon et al., 2006), and is a regulator in the differentiation of hematopoietic stem cells (HSCs) (Shivdasani, Mayer & Orkin, 1995; Zhang et al., 2005; Porcher et al., 1996; Real et al., 2012). In addition, Eto2/Cbfa2t3 regulates the differentiation of HSCs by repressing the expression of target gene Scl/Tal1 (Goardon et al., 2006). Moreover, Ldb1 is a significant transcriptional factor (TF) for the differentiation of erythroid lineage (Soler et al., 2010). According to the ChIPSeq analysis, Ldb1 is necessary for HSCs to control their maintenance since it binds to the majority of enhancer elements in hematopoiesis (Li et al., 2011).

We also included a number of genes with potential regulatory relationship with the three genes Gata1, Gata2 and PU.1. For example, it was indicated that there might be unclear regulations between Gata2 and Gfi1 (Moignard et al., 2013). Gfi1 is an important TF in the regulation of HSCs differentiation (van der Meer, Jansen & van der Reijden, 2010; Lancrin et al., 2012). Gfi1 is required for the differentiation of common lymphoid progenitors (CLPs) and common myeloid progenitors (CMPs) from HSCs and exists in the majority of HSCs, CLPs and CMPs. Similar to gene Gfi1, gene Runx1 is also expressed in most HSCs and progenitor cells as well. Then, Gfi1 and/or Runx1 are expressed continually in most cells which differentiate into the granulocyte lineage (North et al., 2004). Lmo2 is a master regulator of hematopoiesis (Inouea et al., 2013). However, its specific role in regulation is still unclear. Experimental studies suggested that the knockdown of Lmo2 does not affect the expression of Gata1 and Scl/Tal1 (Inouea et al., 2013). However, the overexpression of Lmo2 gene also inhibited erythroid differentiation (Visvader et al., 1997). In addition, gene Ets1 is a suppressor in the erythrocyte differentiation. It is downregulated in erythrocyte differentiation by binding to and activating the Gata2 promoter (Lulli et al., 2006). The last candidate gene is Notch1 that inhibits the differentiation of granulocyte lineage by maintaining the expression of gene Gata2. It also enhances the HSCs differentiate to CLPs (Kumano et al., 2001; Stier et al., 2002). Therefore, in this study we considered the regulatory networks with the following 11 genes: Gata1, Gata2, PU.1/Sfpi1, Runx1, Eto2/Cbfa2t3, Ets1, Notch1, Scl/Tal1, Ldb1, Gfi1 and Lmo2. The detailed information of the references for these 11 genes is also given in Table S2 in Supplemental Information.

Top-down approach: extended forward search algorithm

To reduce the number of unknown parameters in our proposed mathematical model, we used the probabilistic graphical models as the top-down approach to infer the topological structure of gene regulatory networks. Probabilistic graphical model is a useful tool for inferring the network structure (Noor et al., 2013). One type of probabilistic graphical models is the Gaussian graphical model (GGM), which provides a simple and effective method to characterize the regulatory relationship between genes. The GGM is based on the calculation of the conditional dependencies among genes using the gene expression data. The edge connecting two genes in the model is neglected if they are conditionally independent given all other genes (Krämer, Schäfer & Boulesteix, 2009). In this work it is assumed that a system includes genes {G₁, …, G_m} with expression levels x_ij for gene G_i at time point j. Compared with the existing methods that study networks with genes only, this work will study gene networks that include not only genes in the form of monomers {G₁, …, G_m}, which are represented by the linear terms in the model, but also protein heterodimers and/or synergistic effect {G_k–G_l} (k, l = 1,…,m), which are represented by the non-linear terms (NLTs) in the model. There are two reasons for using the NLTs {G_k–G_l}. Firstly, we can use the product of two variables to represent the synergistic effect of these two genes. Secondly, if the NLT represents the protein heterodimer, we assumed that the binding and disassociation reactions for the heterodimer {G_k–G_l} reach an equilibrium state quickly. Thus the level of the heterodimer {G_k–G_l} can be written as C_kl × G_k × G_l, where C_kl is the equilibrium constant. We can consider this constant C_kl as a coefficient in our mathematical model. In both cases, we only need to consider the product of the expression levels of these two genes, namely y_klj = x_kjx_lj, as the level of NLT {G_k–G_l} at time t_j for our algorithm computation. Since the number of possible regulations from NLTs to genes is much larger than that of possible regulations among genes (i.e. 726 vs 110), the regulations from NLTs to genes will dominate the whole genetic regulatory system with high probability. However, the regulations among genes should be the core mechanisms rather than the regulations from NLTs to genes. To avoid the dominance of NLTs regulations, we assume that the number of regulations from NLTs to genes does not exceed that between genes.

According to the GGM (Wang, Myklebost & Hovig, 2003; Wang et al., 2016), we proposed a new algorithm, named Extended Forward Search Algorithm (EFSA), to infer the topological structure of regulatory networks that includes both genes and NLTs. Let X = (x₁, x₂, …, x_N) be a vector that consists of m genes and n NLTs (N = m + n). The following three matrices are constructed, namely a m × m covariance matrix A of m genes, a m × n covariance matrix B to measure the covariance between m genes and n NLTs, and a n × n covariance matrix C of n NLTs. The N-dimensional matrix M is defined by (1) $$\mathbf{M}=\left[\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ \mathbf{B}\mathrm{\prime }& \mathbf{C}\end{array}\right],$$where B′ is the transpose of B. An initial empty graph G is built by the N-dimensional identity matrix. This initial graph G consists of four matrices G₁, G₂, G₃ and G₄ which have the same dimensions as A, B, B′ and C, respectively, namely (2) $$\mathbf{G}=\left[\begin{array}{cc}{\mathbf{G}}_1& {\mathbf{G}}_2\\ {\mathbf{G}}_3& {\mathbf{G}}_4\end{array}\right],$$where G₁ and G₄ are identity matrix with dimensions m and n, respectively, and G₂ and G₃ are m × n and n × m zero matrices, respectively.

The proposed algorithm is given below.

Algorithm 1: extended forward search algorithm

1. Let X = (x₁,x₂,…, x_N) be a vector with N elements, and N be the number of components consist of m genes and n NLTs. An initial empty graph G is built by the N-dimensional identity matrix, which is defined by Eq. (2).

2. Substitute all covariance values from the diagonal positions of sub-matrix A into the corresponding positions of sub-matrix G₁, and then based on the updated G₁, use the Iterative Maximum Likelihood Estimates Algorithm (IMLEA) to compute the new covariance matrix (Dempster, Laird & Rubin, 1977).

3. Add an undirected edge $E_{ij}^1$ ((i, j) ∈ [1, m]²) into G₁, namely add the symmetrical covariance value between the ith gene and jth gene from the positions A(i, j) and A(j, i) into the positions G₁(i, j) and G₁(j, i), respectively. Then compute a new covariance matrix by the IMLEA. Based on the deviance difference between the new covariance matrix and that before addition, test the significance of the added edge $E_{ij}^1$ by using the Chi-square distribution with one degree of freedom. The p-value of the Chi-square test is used in the next step as the edge selection criterion. Record the p-value of this tested edge and remove it from G₁.

4. Add a new undirected edge into G₁. Then, repeat the computation in Step 3. After all possible undirected edges have been tested, sort all tested edges in ascending order by their p-values. If the smallest p-value is lower than the predefined cut-off value, add the edge with the smallest p-value into the sub-graph G₁ permanently.

5. Go back to step 3, add the second edge in the updated sub-graph G₁. Repeat the computation in steps 3 and 4 until the smallest p-value of an added edge is larger than the cutoff p-value.

6. Based on the last updated undirected graph G₁, the graph orientation rules are applied to transform the undirected graph into a directed acyclic graph (DAG) (Meek, 1995). The inferred DAG with m₁ directed edges, denoted as A_s, represents the predicted regulatory network among m genes.

7. Test the possible edges between m genes and n NLTs. Based on the latest matrix G, add an undirected edge $E_{ij}^2$ between the ith gene and the jth NLT. That is, add the symmetrical covariance value between the ith gene and jth NLT from the positions B(i, j) and B′(j, i) into the positions G₂(i, j) and G₃(j, i), respectively. Then, compute a new covariance matrix by the IMLEA. Based on the deviance difference between the new covariance matrix and that before addition, test the significance of the added edge $E_{ij}^2$ by using the Chi-square distribution with one degree of freedom. The p-value of the Chi-square test is used as the edge selection criterion. Record the p-value of this tested edge $E_{ij}^2$ and remove it from G.

8. Repeat the computation in steps 7 for the regulation between genes and NLTs. The last updated sub-graph G₃ with n₁ edges, denoted as B′_s, is the predicted directed regulatory network from n NLTs to m genes. Since we only consider regulations among genes and those from NLTs to genes, the result matrix is given as follows: (3) $${\mathbf{G}}_{\mathbf{s}}=\left[\begin{array}{c}{\mathbf{A}}_{\mathbf{s}}\\ {\mathbf{B}}_{\mathbf{s}}^{\mathrm{\prime }}\end{array}\right].$$

The output network includes m₁ directed edges among m gene and n₁ directed edges from n NLTs to m genes.

Note that we have initially applied the GGM in our previous work to the whole matrix M directly (Wang et al., 2016). However, since the number of NLTs is much larger than that of genes, numerical results showed that the majority of selected edges connect NLTs, but few edges are selected to connect genes. This result is not appropriate because the regulations between genes should be the primary mechanisms of the network. Then we conducted another test, in which we did not consider the regulations between NLTs by changing matrix C into an identity matrix I_m. Matrix M now is (4) $${\mathbf{M}}_1=\left[\begin{array}{cc}{\mathbf{A}}_{\mathbf{s}}& \mathbf{B}\\ \mathbf{B}\mathrm{\prime }& {\mathbf{I}}_{\mathbf{m}}\end{array}\right].$$However, when we applied the GGM to M₁ directly, the singular problem arose during the computation of IMLEA. To satisfy our intention and make the algorithm stable, we proposed EFSA which is executed in two steps. The first step selects regulations between genes and the second step finds regulations from NLTs to genes. The EFSA can be used to predict the gene-gene interactions and the effect from NLTs to genes based on the time-course experimental data.

Bottom-up approach: mathematical model

For a regulatory network with m genes, the expression levels of the i-th gene at time t is denoted as x_i(t). We used the following ordinary differential equation (ODE) model to describe the dynamics of the network (de Jong, 2002) (5) $${\displaystyle \frac{d\mathbf{x}}{dt}=F(t,\mathbf{x})},$$where x = (x₁,…,x_m) is a vector representing the expression levels of m genes. A number of mathematical formalisms have been proposed to describe the dynamical interactions between different genes in the network, such as the models with linear functions (de Jong, 2002) (6) $$F_i(t,\mathbf{x})=\sum _{j=1,j\ne i}^n{a}_{ij}{x}_j-k_i{x}_i$$or the models with non-linear functions (Olariu & Peterson, 2018) (7) $$F_i(t,\mathbf{x})={\displaystyle \frac{{\sum }_{j=1}^n{a}_{ij}{x}_j}{1+{\sum }_{j=1}^n{b}_{ij}{x}_j}-k_i{x}_i}$$

The advantage of the model (Eq. 5) with the linear functions (Eq. 6) is that it has a much smaller number of unknown parameters than the non-linear functions (Eq. 7). However, the non-linear model is able to describe the non-linear dynamics more precisely. Therefore, we proposed a method that combines the feature of additive terms in the linear model and the advantages of non-linear model. We applied the second truncated Taylor series approach to approximate the non-linear function (Eq. 7). Here the Taylor series is a mathematical formula to approximate a function by using a polynomial function (Stewart, 2018). Thus, we proposed an ODE model (Eq. 5) with the following functions (8) $$F_i(t,\mathbf{x})=\sum _{j=1,j\ne i}^m{\mathrm{\alpha }}_{ij}{x}_j+\sum _{1\le j<k\le n}{\mathrm{\beta }}_{ijk}{x}_j{x}_k-k_i{x}_i$$where k_i is the degradation rate of x_i. This proposed model (Eq. 5) with the non-linear function (Eq. 8) is based on the following assumptions:

1. The regulations from different genes to a particular gene are additive. Similarly, the regulations from non-linear terms (NLTs) to a particular gene are also additive.

2. The regulations from gene j to gene i is represented by α_ijx_j, where α_ij is the coefficient of regulation strength.

3. The regulation of NLT x_jx_k to gene i is represented by β_ijkx_jx_k, where β_ijk consists of the regulation strength and equilibrium constant C_ij, as we discussed in the sub-section Top-down Approach.

4. The auto-regulation is not considered, namely α_ii = 0, to avoid confusion between auto-regulation term α_iix_i and degradation term k_ix_i. Note that the issue of auto-regulation may be addressed using a model with non-linear function (Eq. 7). In addition, we just consider the effect of NLTs x_jx_k for j ≠ k since the expression levels of x_j may be highly correlated to that of $x_j^2$. Therefore, we assume that β_ijj = 0.

5. If the value of α_ij is positive (negative or zero), it means that gene x_j activates (represses or has no regulation to) the expression of gene x_i. Similar assumption is applied to the value of β_ijk.

We emphasize that the proposed method in this work is substantially different from our previous work (Wang et al., 2016). The first difference is that the proposed non-linear model (Eq. 8) is different from the non-linear model in Wang et al. (2016). This new model not only can study the regulations from genes to genes, as we considered in our previously proposed model (Wang et al., 2016), but also can investigate the effects of heterodimers and/or synergistic effect in genetic regulation. This new model also leads to the second difference compared with our previous top-down approach, namely the proposed Extended Forward Search Algorithm (EFSA) not only includes the probabilistic graphical model in our previous work (Wang et al., 2016) but also can predict the possible regulations from NLTs to genes. In addition, in this work, we will infer a medium-sized network first by using EFSA and then reduce the network size by removing regulations from the network in the Results section, rather than inferring a core network first and then adding regulations to the core network in our previous approach (Wang et al., 2016).

Parameter inference

When considering the full connected graph among m genes and n non-linear terms (NLTs), we have an ordinary differential equation (ODE) system with m differential equations. The total number of all unknown coefficients is m(m + n). After applying the Extended Forward Search Algorithm (EFSA), we have an inferred regulatory network which contains only m₁ edges among genes and n₁ edges from NLTs to genes. Thus, the numbers of coefficients α_ij and β_ijk are reduced from m(m − 1) to m₁ and from mn to n₁, respectively. It is easier to estimate the parameters for the inferred network than for the fully connected network.

In this work, we used a MATLAB toolbox of Genetic Algorithm to estimate the parameters in the proposed mathematical model (Chipperfield, Fleming & Fonseca, 1994). The algorithm begins by generating a population of initial parameter values, for example, 100 values. Each initial value is called an individual and the whole population is called one generation. Then it calculates the fitness value for each individual of current generation. Based on the fitness values, the algorithm next creates new values for each individual and thus forms a population of the next generation. This process is repeated until a pre-defined number of generations have been calculated. In this work, we used the following functions, namely function crtbp to generate initially binary populations, function reins to effect fitness-based reinsertion, function select to give a convenient interface to the selection routines, function recombine to conduct crossover operators, and function mut to conduct binary and integer mutations. The detailed information of these functions and their alternatives can be found in the relevant reference (Chipperfield, Fleming & Fonseca, 1994).

To ensure the accuracy of estimates, we set the number of generations as 1,000 and the number of individuals for each generation as 300. For the parameter vector (α_ij, β_ijk, k_i), we used the uniform distribution over the interval (W_min, W_max) to generate the initial estimates. Here W_min and W_max are the minimal value and maximal value, respectively, for choosing the samples of the parameters. The values of W_min and W_max are adjusted by computation. For example, if the majority of estimated parameters all are close to W_min, then we will further decrease the value of W_min. However, if the majority of estimated values are well above W_min, then we need to increase the value of W_min accordingly. The similar consideration is applied to W_max. In this study, for the erythroid lineage pathway, numerical results suggest that the values of W_min and W_max for (α_ij, β_ijk, k_i) are ( −3, −3, 0) and (3, 3, 1), respectively. In addition, for the neutrophil lineage pathway, numerical results suggest that the values of W_min and W_max for (α_ij, β_ijk, k_i) are (−2.5, −2.5, 0) and (2.5, 2.5, 1), respectively. We run the algorithm using an initial random number to generate an initial set of model parameters, which leads to a set of estimated parameters. For each model, we used 200 different initial random numbers, which lead to 200 different sets of estimated model parameters. Denote x_i(t_j) and x*_i(t_j) as the observation data and numerical simulations at time point t_j for j = {1,2,…,M}, respectively. The simulation error is calculated by (9) $$E=\sqrt{\sum _{i=1}^m\sum _{j=1}^M(x_i(t_j)-x_i^{\ast }(t_j){)}^2}.$$

We selected the top ten sets with the minimal estimated errors out of 200 estimates for further analysis and comparison.

Inference of regulatory network

To reduce the complexity of regulatory networks, we first used the Extended Forward Search Algorithm (EFSA) to predict the topological structure of genetic networks. The algorithm controls the number of edges by adjusting a pre-defined cut-off value. This value is equivalent to the significant value in statistics. If the threshold is too low, we may miss some significant regulations. However, if the threshold is relatively high, it is quite possible to select insignificant regulations. This work considers the networks including 11 genes and 55 non-linear terms (NLTs). For the sub-network of 11 genes only (i.e. matrix A_s in Eq. 3), to ensure the statistically significant, we set a specific threshold as 0.1 for both the erythroid regulatory network and neutrophil regulatory network. The selection of this threshold value (i.e. 0.1) is based on the balance between neither selecting much insignificant regulations nor choosing a small number of candidate regulations. Then we had 46 and 40 directed edges for the erythroid regulatory network and neutrophil network, respectively.

For the regulations from NLTs to genes (i.e. matrix B′_s in Eq. 3), the size of matrix B′_s is much larger than that of A_s. To avoid the dominance of the regulations from NLTs to genes, we also set the cut-off value as 0.1 for the two networks, or take the first 46 and 40 directed edges from NLTs to genes for the erythrocyte and neutrophil differentiation, respectively, if more edges are selected when using the cut-off value 0.1. The reason we still applied threshold 0.1 here is that the number of selected edges that satisfy this value is much larger than the required number (i.e. 46 for the erythroid regulatory network and 40 for the neutrophil regulatory network). Since the edges are selected and ranked by their significance, we can simply select the top 46 edges and 40 edges for the erythroid and neutrophil pathway, respectively, without conducting any further numerical tests.

Figures S1 and S2 in Supplemental Information present the inferred regulatory networks for the erythroid and neutrophil networks, respectively. Note that there are 11 and 17 isolated NLTs in the erythroid and neutrophil networks, respectively, since no significant edges have been selected from these NLTs by our algorithm. All these isolated NLTs are listed in the “Isolated NLTs Table”. Moreover, all arrows in these figures only represent the direction of regulations, rather than the types of regulations (i.e. positive or negative regulation). We will study the detailed regulatory mechanisms in the next subsection. We found that the targeted gene of the protein heterodimer is a component of that heterodimer in all situations. The possible explanation of this observation is that the expression levels of a heterodimer are the product of the expression levels of the two corresponding genes (namely x_ix_j for genes i and j with expression levels x_i and x_j, respectively). Thus, the expression data of the NLTs {x_ix_j} may be highly correlated to those of the component genes, namely {x_i} or {x_j}.

Inference of mathematical model

After the success of constructing regulatory networks in the previous sub-section, we next study the detailed dynamics of genetic networks in fate determination of hematopoietic stem cells (HSCs) by using our proposed mathematical model. The major step is to infer the values of unknown parameters in the model (Eq. 8). If we consider the fully connected model, there should be 11 × (11 + 55) = 726 parameters. However, after the application of EFSA, the number of unknown parameters is reduced to 103 (including 46 directed edges between genes, 46 directed edges from non-linear terms (NLTs) to genes and 11 self-degradation rate constants) for the differentiation of erythrocytes and 91 (including 40 directed edges between genes, 40 directed edges from NLTs to genes and 11 self-degradation rate constants) for the differentiation of neutrophils. We next applied the Genetic-Algorithm to estimate these unknown parameters for two networks. We used 200 different random numbers to obtain different initial values of rate constants (α_ij, β_ijk, k_i) over the defined range (W_min, W_max), which was discussed in the Methods section. This leads to 200 different sets of estimated parameters. Then, we chose the top ten sets of estimated results for each differentiated lineage with the smallest estimation errors for further robustness analysis. According to the definition of estimation error (Eq. 9), the optimal inferred network for the erythrocyte differentiation in our tests has estimation error 0.9902. In addition, the robust average (Eq. 12) and robust standard deviation (Eq. 13) are 0.3977 and 0.1066, respectively. For the neutrophil differentiation, the optimal inferred network has estimation error 0.8726, robust average 0.3983 and robust standard deviation 0.1275.

Figures 1 and 2 present the simulation results based on the optimal estimated parameters for the expression levels of four genes, namely genes Gata1, PU.1, Ets1 and Tal1, for the differentiation of erythrocyte and neutrophil, respectively. The expression levels of Gata1 increase continuously in both simulated and experimental data during the erythrocyte differentiation. However, during the neutrophil differentiation, experimental data of Gata1 keep fluctuations and then turn to slightly decreasing at the end of differentiation, which is matched by our simulation. For gene PU. 1, both microarray and simulated data decline in the differentiation of erythrocyte but climb during the differentiation of neutrophil. Similarly, the expression levels of Ets1 in microarray data increase during erythrocyte differentiation but decrease during neutrophil differentiation. Simulation results also fit the trends for both differentiation pathways. The experimental data of Tal1 increase with fluctuations during the first 60 h of erythrocyte differentiation, but then rises rapidly after the first 60 h. Our simulated results are consistent with the expression levels of Tal1 with the same trend in expression levels. Thus, our simulation results fit the trend of expression levels of these genes very well during two developmental processes. Figures S3 and S4 in Supplemental Information give the simulation results of the other six genes for the differentiation of erythrocyte and neutrophil, respectively.

Reduction of network model—edge deletion

We have obtained two regulatory networks with 92 directed edges and 80 directed edges for erythroid and neutrophil differentiation, respectively. Next we tested the possibility to delete the potential insignificant edges from our predicted regulatory networks. In the first step, we tested the deletion of regulations from non-linear terms (NLTs) to genes. We removed one edge in each test to form a temporary system model, and then examined the simulation error and robustness property of the new model. Afterwards, we removed one specific edge permanently if the corresponding new system has the minimal change in simulation error and robustness property, and then formed an updated model. This test is repeated until both the simulation error and robustness property of the updated model are much worse than the original network without any removal. In the second step, we evaluated the regulatory interactions between 11 genes using the same method in the first step.

For the erythrocyte differentiation, Table 1 suggests that after removing 3 regulations from NLTs to genes, the estimation error (Eq. 9) is improved (shown in DEL1). Then, we tested the regulation reduction from gene to gene. The final result suggests that, after we deleted (Ldb1 → Lmo2), (Notch1 → Lmo2), (Cbfa2t3 → Lmo2) and (Runx1 → Lmo2) edges, the estimation error (Eq. 9) is slightly increased. However, the robustness property is better than that of the DEL1 model since the robust average (Eq. 12) is decreased. Thus, for the erythroid differentiation, numerical tests recommended to remove total seven edges from our predicted regulatory network. We stopped the deletion test after obtaining the DEL5 model. If we proceed further deletion, both simulation error and robustness property of the temporary network are much worse than the original network without removal.

Table 2 shows that, for the neutrophil differentiation, there are no insignificant regulations from NLTs to genes, because the removal of any edge from NLTs to genes will increase the simulation error (Eq. 9) substantially and/or decrease the robustness property by increasing the values of robust average (Eq. 12) and robust standard deviation (Eq. 13). For the regulations between genes, we have removed the following four regulations, namely (Gata2 → Ldb1), (Runx1 → Cbfa2t3), (Ldb1 → Lmo2) and (Tal1 → Lmo2), and formed an updated system. Table 2 shows that the simulation error and robustness property of the updated system are close to those of the original system without any removal of edges. Thus, for the neutrophil differentiation, numerical tests recommended to remove only four edges from our predicted regulatory network. Coincidentally, we stopped the deletion test after obtaining the DEL5 model because of the same reason for the erythrocyte differentiation.

Figures 3 and 4 present the inferred regulatory networks after edge deletion test for erythroid and neutrophil differentiation, respectively. Initially, we have 92 directed edges for the erythrocyte pathway and 80 directed edges for the neutrophil pathway. After the edges deletion, seven and four directed edges have been taken away from the erythrocyte network and neutrophil network, respectively, since the removal of these edges has not much negative influence on simulation error (Eq. 9), robust average (Eq. 12) and robust standard deviation (Eq. 13). Thus, there are 85 and 76 directed edges left for the erythrocyte and neutrophil pathways, respectively.