Identifying regulational alterations in gene regulatory networks by state space representation of vector autoregressive models and variational annealing

Vector autoregressive model

where A is a p × p autoregressive coefficient matrix, and ε _t is observation noise at time t and follows $\mathcal{N}\left(0,\mathsf{\text{diag}}\left[{\sigma }_1^2,...,{\sigma }_p^2\right]\right)$, a normal distribution with mean 0 and variance $\mathsf{\text{diag}}\left[{\sigma }_1^2,...,{\sigma }_p^2\right]$. The (i, j)th element of A, A _ij , indicates a temporal regulation from the j th gene to the i th gene, and if A _ij ≠ 0, regulation from the j th gene to the i th gene is considered. By examining whether A _ij is zero or not for all i and j, a gene network is constructed. Since equally space time points are assumed in the VAR model, it has difficulty on handling unequally spaced time series data.

State space representation of VAR model (VAR-SSM)

where ρ _t is the observation noise normally distributed with mean 0 and variance R = diag[r₁,..., r _p ]. Unequally spaced time series data are handled by ignoring observation model at non-observed time points.

Effectiveness of variational annealing

As alternatives of the variational annealing, we may consider the variational method and the EM algorithm where X is the set of hidden variables and E is handled as the set of parameters to be maximized. We show the effectiveness of variational annealing compared to the variational method and the EM algorithm under the following conditions: P(X, Θ, E) is factorized into P(X, E)P(Θ, E), and P(E _i |X, Θ, E\{E _i }) is given as a binomial distribution, where E _i is the i th element of E. If factorization of P(X, Θ, E) = Q(X)Q(Θ)Q(E) is assumed in the calculation of the variational method, arg max _E Q(E) is not the optimal solution of Equation (2) in general. In the EM algorithm, by allowing E to move around a p-dimensional continuous space ${\left(0,1\right)}^p,{Ê}_{\mathsf{\text{EM}}}=arg\underset{E\in {\left(0,1\right)}^m}{max}P\left(E\right)$ can be calculated. Let ${Ê}_{\mathsf{\text{EM}},i}$ be the i th element of ${Ê}_{\mathsf{\text{EM}}}$. Usually, ${Ê}_{\mathsf{\text{EM}}}\mathsf{\text{,}}i$ is mapped to 1 if ${Ê}_{\mathsf{\text{EM,}}\mathsf{\text{i}}}>0.5$ and 0 otherwise for discretizing ${Ê}_{\mathsf{\text{EM}}}$ to the p-dimensional binary space {0,1} ^p , but such a mapping is not guaranteed to provide $arg\underset{E\in {\left\{0,1\right\}}^m}{max}P\left(E\right)$. Although other mappings can be considered, to the best of our knowledge, no mapping is guaranteed to provide the optimal solution in polynomial time of p.

In the following, we prove a proposition in order to show that the variational annealing possibly give the optimal solution of Equation (2) even if the factorization of $Q\left(E\right)={\prod }_{i=1}^{p}Q\left(E_i\right)$ is additionally considered.

Then, the set of E _i ∈ {0,1} maximizing ${Ê}_{\mathsf{\text{EM}}}$is $arg\underset{E\in {\left\{0,1\right\}}^p}{max}P\left(E\right)$.

For the proof of the proposition, see Section 1 in Additional file 1. From Proposition 1, if the factorization P(X, Θ, E) = P(Θ, E)P(X, E) is satisfied and optimal Q functions are found, the variational annealing is guaranteed to provide the optimal solution of Equation (2) while the variational method and the EM algorithm are not. Although the factorization is not a generally satisfied property, the factorization is often assumed in approaches based on the variational method, and the assumption usually works as good approximations. Thus, the variational annealing is expected to provide the better performance than the variational method and the EM algorithm even if the factorization is not satisfied exactly.

Variational E-step

Parameters of Q(X) are mean of x _t , variance of x _t , and cross time variance of x_t-1and x _t . These parameters can be calculated via variational Kalman filter by using following terms expected with Q(Θ)Q(E): 〈E^(c)〉_Q(Θ)Q(E), 〈A〉_Q(Θ)Q(E), 〈H^-1 A ○ E^(c)〉_Q(Θ)Q(E), and 〈(A ○ E^(c))'H^-1A ○ E^(c)〉_Q(Θ)Q(E). For the details of variational Kalman filter, see [4]. From the parameters of Q(X), expectations of the following terms with Q(X) required in other steps are calculated: ${〈{\mathit{x}}_t^{\left(c\right)}〉}_{Q\left(X\right)},{〈{\mathit{x}}_t^{\left(c\right)}{\left({\mathit{x}}_t^{\left(c\right)}\right)}^{\prime }〉}_{Q\left(X\right)}$, and ${〈{\mathit{x}}_{t+1}^{\left(c\right)}{\left({\mathit{x}}_t^{\left(c\right)}\right)}^{\prime }〉}_{Q\left(X\right)}$.

Variational M-step

Parameters for the above functions are calculated by using ${〈{\mathit{x}}_t^{\left(c\right)}〉}_{Q\left(X\right)},{〈{\mathit{x}}_t^{\left(c\right)}{\left(x_t^{\left(c\right)}\right)}^{\prime }〉}_{Q\left(X\right)},{〈{\mathit{x}}_{t+1}^{\left(c\right)}{\left({\mathit{x}}_t^{\left(c\right)}\right)}^{\prime }〉}_{Q\left(X\right)},$, and ${〈E_{ij}^{\left(c\right)}〉}_{Q\left(E\right)}$.

Variational A-step

For the calculation of Q(E), we assume the factorization of Q(E) to ${\prod }_c{\prod }_{ij}Q\left(E_{ij}^{\left(c\right)}\right)$ in order to make the computation tractable. $Q\left(E_{ij}^{\left(c\right)}\right)$ follows a binomial distribution that takes one with probability $e_{ij}^{\left(c\right)}$ and zero with probability $1-e_{ij}^{\left(c\right)}$, and thus the expectation of $E_{ij}^{\left(c\right)}$ with Q(E) is given by $e_{ij}^{\left(c\right)}$. For the preparation, we calculate ${〈A〉}_{Q\left(\mathrm{\Theta }\right)},{〈H^{-1}A〉}_{Q\left(\mathrm{\Theta }\right)},{〈A^{\prime }{H}^{-1}A〉}_{Q\left(\mathrm{\Theta }\right)},{〈{\mathit{x}}_t^{\left(c\right)}〉}_{Q\left(X\right)},{〈{\mathit{x}}_t^{\left(c\right)}{\left({\mathit{x}}_t^{\left(c\right)}\right)}^{\prime }〉}_{Q\left(X\right)}$, and ${〈{\mathit{x}}_{t+1}^{\left(c\right)}{\left({\mathit{x}}_t^{\left(c\right)}\right)}^{\prime }〉}_{Q\left(X\right)}$. $Q\left(E_{ij}^{\left(c\right)}\right)$ is then iteratively calculated by using these expected terms as well as ${〈E_{ik}^{\left(c\right)}〉}_{Q\left(E\right)}$ for k ≠ j. A few iterations are enough for the convergence.

Update and selection of hyperparameters

The proposed model contains u₀, k₀, v₀, l₀, ζ_{i 0}, ζ_{i 1}, α₀, and α₁ as hyperparameters. α₀ and α₁ should be α₀ <α₁ as in the model setting, but this condition can be violated in the update step of the variational method. Thus, we select α₀ and α₁ by cross validation, and update other hyperparameters as in the variational method.

${\widehat{\zeta }}_{i0}$ and ${\widehat{\zeta }}_{i1}$ can also be obtained by the Newton-Raphson method.

For the selection of α₀ and α₁, we set α₁ to some large value and select α₀ by a leave one out cross validation procedure. For condition c ∈ {1, 2} and time point $t\in {\mathcal{T}}_{obs}^{\left(c\right)}$, we remove ${\mathit{y}}_t^{\left(c\right)}$ from data set y and use the data set to train the model. We calculate square sum of residues $r_{t,c}^2$ between ${\mathit{y}}_t^{\left(c\right)}$ and the prediction of ${\mathit{x}}_t^{\left(c\right)}$ estimated from the variational Kalman filter on the trained model given by $r_{t,c}^2={\left({\mathit{y}}_t^{\left(c\right)}-{\mathit{x}}_t^{\left(c\right)}\right)}^{\prime }\left({\mathit{y}}_t^{\left(c\right)}-{\mathit{x}}_t^{\left(c\right)}\right)$. By grid search on parameter space of α₀, we select α₀ that minimizes ${\sum }_{c=1}^2{\sum }_{t\in {\mathcal{T}}_{obs}^{\left(c\right)}}{r}_{t,c}^2$.

Summary of procedures

In our setting, α₁ is set to 1,000. For the initialization of τ and other hyperparameters, we use the following settings: $\tau =2.5,E_{ij}^{\left(c\right)}=0.5,u_i=1,k_i=1,v_i=1,l_i=1,{\zeta }_{i0}=10,\text{ and }{\zeta }_{i1}=10$. If ${\mathit{y}}_t^{\left(c\right)}$ is observed, we initialize ${\mathit{x}}_t^{\left(c\right)}$ with ${\mathit{y}}_t^{\left(c\right)}$. Otherwise, we use the linearly interpolated one.

Comparison between variational annealing and EM algorithm

From the comparison, the results of the proposed approach contain more true positives than those of the EM algorithm based approach except for identifying changes on regulations for time points 50. For identifying changes on regulations, the EM algorithm based approach estimates bit more true positives than the proposed approach, but the difference is so small that it can be ignored. On the other hand, the results of the EM algorithm based approach contain more false positives than those of the proposed approach, and hence the precision of the results by the EM algorithm is worse than that of the proposed approach. Therefore, the effectiveness of the variational annealing is confirmed in the computational experiment as well.

Comparison between proposed approach and existing approaches

We employ the elastic net based VAR model approach [2] and the dynamic Bayesian network based approach termed G1DBN [8] as existing approaches for estimating networks from time series data. For the experiments, two versions of these approaches are considered: ENet1 and ENet2 from the elastic net based VAR model approach, and G1DBN1 and G1DBN2 from G1DBN. These approaches are different in the following point: ENet1 and G1DBN1 estimate G₁ and G₂ independently, i.e., for the estimation of G₁, only time series data from G₁ is used, while ENet2 and G1DBN2 assume that G₁ and G₂ have the same network structure and estimate a network by using two time series data. Thus, ENet2 and G1DBN2 are considered to more use data sample than ENet1 and G1DBN1 for network estimation, but changes on regulations between G₁ and G₂ are not considered. For selection of hyperparameters in ENet1 and ENet2, AICc is used [2], and for hyperparameters α₁ and α₂ of G1DBN1 and G1DBN2, a setting of α₁ = 0.1 and α₂ = 0.0059 considered in [8] is used.

are also provided. The results are averaged on ten data sets. The number of regulations in the true network models of G₁ and G₂ are in total 305. For the latter point, we consider the estimated regulations existing only in one of two estimated networks as changed regulations, and check if they correctly exist only in the corresponding true network. The numbers of true positives and false negatives of the estimated changes on regulations and the precisions are summarized in Table 2(b). The results are also averaged on ten data sets. The number of true changes on regulations between in G₁ and G₂ are 47, i.e., the number of true positives on this case is at most 47.

For the estimation of the regulations in Table 2(a), the proposed approach outperforms other approaches in terms of true positives. The proposed approach contains more false positives than ENet1, G1DBN1, and G1DBN2 on the data of 25 time points. We further consider these cases in terms of the precision. The precisions of the proposed approach, ENet1, G1DBN1, and G1DBN2 are given by 0.77, 0.62, 0.68, and 0.74, respectively. From this analysis, the proposed approach shows better performance than ENet2, G1DBN1, and G1DBN2 on the whole. More true positives are estimated by ENet2 than ENet1, and the precisions of ENet2 tend to be better than those of ENet1. This type of relationship is also observed between G1DBN1 and G1DBN2. Also, the elastic net based VAR model approaches estimate more true positives than the approaches from G1DBN. However, in the precision, the approaches from G1DBN are better than the elastic net based VAR model approaches. From the results in Table 2(b), we see that the proposed approach estimates more true changes on regulations than ENet1 and G1DBN1 on both data of 25 and 50 time points. In addition, the results of the proposed approach contain less false positives than those of ENet1 and G1DBN1. No changes on regulations are detected in ENet2 and G1DBN2 as topological differences are ignored in these approaches. Since the proposed approach considers differences on network structures as well as uses two time series data efficiently, it can provide better results than other approaches on both estimating regulations and identifying changes on regulations.

One may think it is strange that false positives in ENet1 and G1DBN1 in Table 2(b) is more than those in Table 2(a), but this case can occur from the following reason. If a regulation exists in both of the true network models of G₁ and G₂, but is estimated only for G₁, then the case is not counted as a false positive in Table 2(a) while it is counted as a false positive in Table 2(b). Thus, the number of false positives in Table 2(b) can be greater than those in Table 2(a).

In order to show the performance in unequally spaced time series data, we generate unequally spaced time series data of 25 and 50 observed time points. For time series data of 25 observed time points, we first generate equally spaced time series data of 40 time points and divide it into three blocks: 15 time points, 10 time points, and 15 time points. We then remove time points in the following manner: no time point is removed in the first block; one of every two time points are removed in the second block; and two of every three time points are removed in the third block. Figure 3 shows the time point schedule of time series data obtained in this process. For time series data of 50 observed time points, we first generate equally space time series data of 80 time points, divide it into three blocks: 30 time points, 20 time points, and 30 time points. Then, some time points are removed in a similar manner. We apply the proposed approach, ENet1, ENet2, G1DBN1, and G1DBN2 to the unequally spaced time series data. Results for the dataset are also summarized in Tables 2(a) and 2(b). From the comparison of results on equally and spaced time series data, the results of the proposed approach from unequally spaced time series data are worse than those from equally spaced one even with the same number of observed time points. However, results of ENet1, ENet2, G1DBN1, and G1DBN2 are worsened more than those of the proposed approach. This is probably because unequally spaced time points break their assumption, and their estimation process is misled.