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

Abstract

Background

In the analysis of effects by cell treatment such as drug dosing, identifying changes on gene network structures between normal and treated cells is a key task. A possible way for identifying the changes is to compare structures of networks estimated from data on normal and treated cells separately. However, this approach usually fails to estimate accurate gene networks due to the limited length of time series data and measurement noise. Thus, approaches that identify changes on regulations by using time series data on both conditions in an efficient manner are demanded.

Methods

We propose a new statistical approach that is based on the state space representation of the vector autoregressive model and estimates gene networks on two different conditions in order to identify changes on regulations between the conditions. In the mathematical model of our approach, hidden binary variables are newly introduced to indicate the presence of regulations on each condition. The use of the hidden binary variables enables an efficient data usage; data on both conditions are used for commonly existing regulations, while for condition specific regulations corresponding data are only applied. Also, the similarity of networks on two conditions is automatically considered from the design of the potential function for the hidden binary variables. For the estimation of the hidden binary variables, we derive a new variational annealing method that searches the configuration of the binary variables maximizing the marginal likelihood.

Results

For the performance evaluation, we use time series data from two topologically similar synthetic networks, and confirm that our proposed approach estimates commonly existing regulations as well as changes on regulations with higher coverage and precision than other existing approaches in almost all the experimental settings. For a real data application, our proposed approach is applied to time series data from normal Human lung cells and Human lung cells treated by stimulating EGF-receptors and dosing an anticancer drug termed Gefitinib. In the treated lung cells, a cancer cell condition is simulated by the stimulation of EGF-receptors, but the effect would be counteracted due to the selective inhibition of EGF-receptors by Gefitinib. However, gene expression profiles are actually different between the conditions, and the genes related to the identified changes are considered as possible off-targets of Gefitinib.

Conclusions

From the synthetically generated time series data, our proposed approach can identify changes on regulations more accurately than existing methods. By applying the proposed approach to the time series data on normal and treated Human lung cells, candidates of off-target genes of Gefitinib are found. According to the published clinical information, one of the genes can be related to a factor of interstitial pneumonia, which is known as a side effect of Gefitinib.

Background

Gene network estimation from time series gene expression data is a key task for elucidating cellular systems. Thus far, wide variety of approaches have been proposed based on the vector autoregressive (VAR) model [1, 3], the state space model [4–6], and the dynamic Bayesian network [7, 8]. Recently, time series gene expression data on multiple conditions aiming at analyzing effects of cell treatment such as drug dosing and heat shock are available. We here assume that some gene regulations are disrupted but many of the gene regulations do not change due to some treatment of interest, and try to find a small number of changes on regulations as keys for elucidating effects by the treatment.

A possible way for finding changes on regulations is to estimate networks from two data sets separately and then compare their structures. However, due to the limited length of time series data (usually less than 10 time points) and unignorable measurement noise, networks are estimated with high error rates and the estimation errors cause the serious failure on identifying changes on regulations. Thus, approaches using two time series data in an efficient manner are strongly demanded. Also, widely used statistical methods such as the VAR model and dynamic Bayesian network assume equally spaced time points in time series data. However, observed time points on usually available time series data are not equally spaced [5, 6, 9], and approaches that can handle unequally spaced time series data in a theoretically correct way should be considered.

We propose a new statistical model that estimates gene networks on two different conditions in order to identify changes on regulations between the conditions. As the basis of the proposed model, we employ the state space representation for VAR model (VAR-SSM), in which observation noise is considered between the measured or observed gene expressions and the true gene expressions in observation model and gene regulations between true gene expressions are considered in the system model [10]. The VAR-SSM can handle unequally spaced time series data by ignoring observation model on the non-observed time points. For considering the changes on regulations, we introduce hidden variables to the VAR-SSM in order to indicate the presence of regulations in each condition. If hidden binary variables on two conditions indicating the presence of a regulation are both estimated as one, the regulation is considered as a commonly existing regulation. On the other hand, if only one of the hidden binary variables for the regulation is estimated as one, the regulation is considered as a condition specific regulation. We also introduce a potential function between the hidden binary variables that is designed to take high probability if the hidden binary variables on two conditions take the same value. From the design of the potential function, the similarity of networks on two conditions is automatically considered. Since the time series data on both conditions are used for estimating commonly existing regulation due to the use of the hidden binary variables, an efficient data assignment is achieved. In addition, from the more accurate estimation of commonly existing regulations by the efficient data assignment, accurate identification of changes on regulations is induced.

The hidden binary variables are estimated by searching the configuration of binary variables that maximizes the marginal likelihood of the model. However, searching the optimal configuration is computationally intractable. Thus, as an alternative approach, we derive a new variational annealing method based on [11] in order to estimate the hidden binary variables. We also give a proof for the effectiveness of the variational annealing compared to other candidate alternatives, the variational annealing and the EM algorithm, in order to show the validity of using the variational annealing.

For the performance evaluation, we generate two regulatory networks in such a way that most of the regulations commonly exist and some exist only on one of the networks. We then apply our proposed approach and existing var model based and dynamic Bayesian network based approaches to two equally spaced time series data drawn separately from the generated networks. From the comparisons of true positive rates and false positive rates of these approaches, we confirm the effectiveness of our approach. We also generate unequally spaced time series data from these networks, and show that our approach works correctly on unequally spaced time series data while the performance of the existing approaches assuming equally spaced time points is drastically worsened.

Our proposed approach is used to analyze changes on regulations in gene networks between normal Human lung cells and Human lung cells treated by stimulating EGF-receptors and dosing an anticancer drug termed Gefitinib. A lung cancer condition is simulated by the stimulation of EGF-receptors in the treated cells. Since Gefitinib is known as a selective inhibitor of EGF-receptors, the stimulation of EGF-receptors would be counteracted by Gefitinib, and hence the treated cells are expected to be the same condition as normal cells. However, gene expression profiles from normal and treated cells are actually different, and off-targets of Gefitinib causing unexpected positive or negative effects are implied. We focus on genes with changes on regulations between the networks estimated by our approach and find possible off-target genes of Gefitinib. According to the published clinical information, one of the possible off-target genes is suggested as one of factors of interstitial pneumonia, which is known as a side effect of Gefitinib.

Methods

Vector autoregressive model and its state space representation

Vector autoregressive model

Given gene expression profile vectors of p genes during T time points {y₁, ..., y _T }, the first order vector autoregressive (VAR(1)) model at time point t is given by

{\mathit{y}}_t=A{\mathit{y}}_{t-1}+{\epsilon }_1,

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)

Let \mathcal{T} be the set of equally spaced entire T time points and {\mathcal{T}}_{obs} the set of time points where gene expressions are observed. Note that {\mathcal{T}}_{obs}\subseteq \mathcal{T} holds. VAR-SSM is comprised of two models: system model and observation model. Let x _t be hidden variable vector representing true gene expression at time t. The system model is given as the VAR model of x _t :

{\mathit{x}}_t=A{\mathit{x}}_{t-1}+{\mathit{\eta }}_t,t\in \mathcal{T},

where η _t is the system noise normally distributed with mean 0 and variance H = diag[h₁,...,h _p ]. The observation model represents measurement error of observed gene expression y _t and true gene expression x _t at observed time point t\in {\mathcal{T}}_{obs}:

{\mathit{y}}_t={\mathit{x}}_t+{\mathit{\rho }}_t,t\in {\mathcal{T}}_{obs},

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.

Joint model of VAR-SSM for two time series data

Let {\left\{{\mathit{y}}_t^{\left(c\right)}\right\}}_{t\in {\mathcal{T}}_{obs}^{\left(c\right)}} be time series gene expression data on cell condition c, where {\mathcal{T}}_{obs}^{\left(c\right)} is the set of observed time points on cell condition c. We also let {\mathcal{T}}^{\left(c\right)} be the set of time points from 1 to T^(c), where T^{\left(c\right)}=max\left\{t\in \underset{obs}{\overset{\left(c\right)}{\mathcal{T}}}\right\}. Given time series data on two types of cell conditions c = 1 and 2, we propose a new VAR-SSM model to estimate gene networks in the two conditions as well as identify changes on regulations between them. The model is comprised of the following two equations:

\begin{array}{c}{\mathit{x}}_t^{\left(c\right)}=A\circ E^{\left(c\right)}{\mathit{x}}_{t-1}^{\left(c\right)}+{\eta }_t^{\left(c\right)},t\in {\mathcal{T}}^{\left(c\right)},\\ {\mathit{y}}_t^{\left(c\right)}={\mathit{x}}_t^{\left(c\right)}+{\rho }_t^{\left(c\right)},t\in {\mathcal{T}}_{obs}^{\left(c\right)},\end{array}

where ∘ denotes the Hadamard product, E^(c)is a p × p binary matrix, and {\mathit{\eta }}_t^{\left(c\right)} and {\mathit{\rho }}_t^{\left(c\right)} are respectively system and observation noises from \mathcal{N}\left(0,H\right) and \mathcal{N}\left(0,R\right). In this model, the (i, j)th element E_{ij}^{\left(c\right)} takes one if regulation from gene j to gene i exists on condition c and zero otherwise, i.e., the presence of regulations is controlled by E^(c)and the AR coefficient matrix A is commonly used in conditions 1 and 2. Changes on regulations are identified when regulations exist only in a condition.

The complete likelihood of our model, P(Y, X, Θ, E), where Y, X, Θ, and E are respectively the sets of {\mathit{y}}_t^{\left(c\right)},{\mathit{x}}_t^{\left(c\right)}, parameters, and E^(c), is given by the following equation:

\begin{array}{lll}\hfill P\left(Y,X,\mathrm{\Theta },E\right)& =\prod _{c=1}^2\prod _{t\in {\mathcal{T}}^{\left(c\right)}}\frac{|H{|}^{-1∕2}}{{\sqrt{2\pi }}^p}exp\left\{-\frac{1}{2}{\left({\mathit{x}}_t^{\left(c\right)}-A\circ E^{\left(c\right)}{\mathit{x}}_{t-1}^{\left(c\right)}\right)}^{\prime }{H}^{-1}\left({\mathit{x}}_t^{\left(c\right)}-A\circ E^{\left(c\right)}{\mathit{x}}_{t-1}^{\left(c\right)}\right)\right\}& \hfill \text{(1)}\\ \times \prod _{t\in {\mathcal{T}}_{obs}^{\left(c\right)}}\frac{|R{|}^{-1∕2}}{{\sqrt{2\pi }}^p}exp\left\{-\frac{1}{2}{\left({\mathit{y}}_t^{\left(c\right)}-{\mathit{x}}_t^{\left(c\right)}\right)}^{\prime }{R}^{-1}\left({\mathit{y}}_t^{\left(c\right)}-{\mathit{x}}_t^{\left(c\right)}\right)\right\}P\left(\mathrm{\Theta },E\right),& \hfill \text{(2)}\\ \hfill \text{(3) }\end{array}

where the prior distribution P(Θ, E) is assumed to be factorized as

P\left(\mathrm{\Theta },E\right)=\prod _{i}P\left(h_i\right)P\left(r_i\right)\prod _{j}P\left(A_{ij}\right)P\left(E_{ij}^{\left(1\right)},E_{ij}^{\left(2\right)},z_{ij}\right).

Here, z _ij is a parameter for a potential function of E_{ij}^{\left(1\right)} and E_{ij}^{\left(2\right)} defined later. The prior distributions of A _ij is given by

P\left(A_{ij}|h_i,E_{ij}^{\left(1\right)},E_{ij}^{\left(2\right)}\right)=\mathcal{N}{\left(A_{ij};0,h_i\cdot {\alpha }_1\right)}^{F_{ij}}\mathcal{N}{\left(A_{ij};0,h_i\cdot {\alpha }_0\right)}^{1-F_{ij}},

where α₀ and α₁ are parameters controlling the shrinkage of coefficients A and F _ij is a binary variable that takes 1 if E_{ij}^{\left(1\right)} or E_{ij}^{\left(2\right)} takes 1 and 0 otherwise, i.e., F _ij is given by 1-{\prod }_{c=1}^2\left(1-E_{ij}^c\right). α₁ is set to a large value, while α₀ is set to smaller than α₁. From the design of the prior for A _ij , if E_{ij}^{\left(1\right)} or E_{ij}^{\left(2\right)} takes one, i.e., there exists regulation from gene j to i in condition 1 or 2, weaker prior \mathcal{N}\left(A_{ij};0,h_i\cdot {\alpha }_1\right) is selected and the shrinkage of the coefficients are avoided. Otherwise stronger prior \mathcal{N}\left(A_{ij};0,h_i\cdot {\alpha }_0\right) is selected and the sparsity of the network structures is promoted. The prior distributions of h _i , and r _i are given by

\begin{array}{c}P\left(h_i\right)=\mathcal{I}\mathcal{G}\left(h_i;u_0,k_0\right),\\ P\left(r_i\right)=\mathcal{I}\mathcal{G}\left(r_i;v_0,l_0\right),\end{array}

where \mathcal{I}\mathcal{G} represents the density function of inverse gamma distribution, u₀ and v₀ are the shape parameters, and k₀ and l₀ are the inverse scaling parameters. Under the assumption that a small number of regulations change between two conditions, we design the prior distribution for E_{ij}^{\left(1\right)} and E_{ij}^{\left(2\right)},P\left(E_{ij}^{\left(1\right)},E_{ij}^{\left(2\right)},z_{ij}\right) by using the following potential function between E_{ij}^{\left(c\right)} for c = 1 or 2:

\frac{1}{2}\varphi \left(E_{ij}^{\left(1\right)},E_{ij}^{\left(2\right)};z_{ij}\right)=\left\{\begin{array}{ccc}\hfill z_{ij}\hfill & \hfill \mathsf{\text{if}}\hfill & \hfill E_{ij}^{\left(1\right)}\ne E_{ij}^{\left(2\right)}\hfill \\ \hfill 1-z_{ij}\hfill & \hfill \mathsf{\text{otherwise}}\hfill \end{array}\right..

In this setting, if z _ij is small, E_{ij}^{\left(1\right)} and E_{ij}^{\left(2\right)} tend to take the same value and thus most of the regulations exist in both two conditions. We also introduce a prior distribution for z _ij by beta distribution with parameters ζ_{i 0}and ζ_{i 1}:

\mathcal{B}\left(z_{ij};{\zeta }_{i0},{\zeta }_{i1}\right)=\frac{1}{B\left({\zeta }_{i0},{\zeta }_{i1}\right)}{z}_{ij}^{{\zeta }_{i0}-1}{\left(1-z_{ij}\right)}^{{\zeta }_{i1}-1},

where B(·) is the beta function. Thus, the prior distribution of E_{ij}^{\left(c\right)} is given by

P\left(E_{ij}^{\left(1\right)},E_{ij}^{\left(2\right)},z_{ij}\right)=\varphi \left(E_{ij}^{\left(1\right)},E_{ij}^{\left(2\right)}\right)B\left(z_{ij};{\zeta }_{i0},{\zeta }_{i1}\right).

Figure 1 shows a graphical representation of the proposed model, where dependency of the parameters and variables are indicated. The hyperparameters are omitted and the observed data y _t is represented by gray nodes. In the observed data y_t^{\left(1\right)} and y_t^{\left(2\right)} are propagated mainly via hidden variables E_{ij}^{\left(1\right)} and E_{ij}^{\left(2\right)}. Due to the data propagation, more accurate estimation is expected in the proposed model than the approaches considering data on two conditions independently.

Figure 1

A graphical representation of the proposed model. Hyperparameters are omitted from this representation. The nodes in gray denote observed data.

For the parameter estimation, we search the configuration of E maximizing the following marginal likelihood:

Ê=arg\underset{E}{max}\int dX\int d\mathrm{\Theta }P\left(Y,X,\mathrm{\Theta },E\right).

(1)

Finding the optimal configuration of E is computationally intractable, and heuristics approaches such as the EM algorithm and the variational method are used in practice. Here, we use the variational annealing, an extension of the deterministic annealing for discrete variables [11]. In the next section, we give a small explanation of the variational annealing and show its effectiveness compared to the EM algorithm and the variational method.

Parameter estimation by variational annealing

In the deterministic annealing, optimization problem is solved while gradually changing temperature in a some schedule, and maximum likelihood estimator is obtained like the EM algorithm [12–14]. Yoshida and West proposed to use the deterministic annealing to find the configuration of the binary variables that maximizes the likelihood of factor models with sparseness priors [11]. We derive a new variational annealing method by extending Yoshida and West's approach to find the configuration of the binary variables on marginal likelihood function, which can be applied for searching E that maximizes Equation (1).

Let E, X, and Θ be p dimensional binary variables, unobserved variables, and parameters, respectively, and consider to search E maximizing the following marginal likelihood:

\underset{E\in {\left\{0,1\right\}}^p}{max}log\int dX\int d\mathrm{\Theta }P\left(X,\mathrm{\Theta },E\right).

(2)

The maximum of the marginal likelihood on E is bounded by the following formula:

\underset{E\in {\left\{0,1\right\}}^p}{max}log\int dX\int {d\mathrm{\Theta }P\left(X,\mathrm{\Theta },E\right)\ge \tau log{\int }_{{\left(0,1\right)}^p}dE\left[\int dX\int d\mathrm{\Theta }P\left(X,\mathrm{\Theta },E\right)\right]}^{1∕\tau },

(3)

where τ is called temperature and the equality holds for τ → +0. Hereafter, the integral range of E is omitted if no confusion occurs. Let Q(E) be a normalized non-negative function, i.e., Q(E) ≥ 0 and ∫ Q(E)dE = 1.

From the Gibbs inequality, the right side of Equation (3) is also bounded:

\tau log\int dE{\left[\int dX\int d\mathrm{\Theta }P\left(X,\mathrm{\Theta },E\right)\right]}^{1∕\tau }\ge \int dEQ\left(E\right)log\frac{\int dX\int d\mathrm{\Theta }P\left(x,\mathrm{\Theta },E\right)}{Q{\left(E\right)}^{\tau }}.

Here, Q(E) is considered as an approximation function of \frac{P{\left(E\right)}^{1∕\tau }}{\int P{\left(E^{\prime }\right)}^{1∕\tau }d{E}^{\prime }}. Under the assumption that P(X,Θ,E) ∝ Q(X)Q(Θ), where Q(X) and Q(Θ) are normalized non-negative functions, we have the following inequality

\int dEQ\left(E\right)log\frac{\int dX\int d\mathrm{\Theta }P\left(X,\mathrm{\Theta },E\right)}{Q{\left(E\right)}^{\tau }}\ge \int dE\int dX\int d\mathrm{\Theta }Q\left(X\right)Q\left(\mathrm{\Theta }\right)Q\left(E\right)log\frac{P\left(X,\mathrm{\Theta },E\right)}{Q\left(X\right)Q\left(\mathrm{\Theta }\right)Q{\left(E\right)}^{\tau }}.

(4)

Thus, as an approximation of E maximizing Equation (2), we try to find Ê=arg\underset{E}{max}Q\left(E\right), where Q(E) is the function maximizing the lower bound in Equation (4). Since higher values on P(E) are weighed more in the approximation by Q(E) for τ < 1, the better approximation is expected for the higher values. This property on the limited case is shown in Proposition 1. As is in the variational method and the EM algorithm, the maximum of the lower bound in Equation (4) is searched by a hill climbing from high temperature τ > 1.

In the hill climbing, Q(X), Q(Θ), and Q(E) are alternately updated from the following equations:

\begin{array}{c}Q\left(E\right)\propto exp\left(\frac{1}{\tau }\int dX\int d\mathrm{\Theta }Q\left(X\right)Q\left(\mathrm{\Theta }\right)logP\left(E|X,\mathrm{\Theta }\right)\right),\\ Q\left(X\right)\propto exp\left(\int dE\int d\mathrm{\Theta }Q\left(E\right)Q\left(\mathrm{\Theta }\right)logP\left(X|\mathrm{\Theta },E\right)\right),\\ Q\left(\mathrm{\Theta }\right)\propto exp\left(\int dE\int dXQ\left(E\right)Q\left(X\right)logP\left(\mathrm{\Theta }|X,E\right)\right).\end{array}

Gradually converging τ to 0, local optimum of the lower bound and corresponding Q(E), Q(X), and Q(Θ) are obtained.

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.

Proposition 1. P(X,Θ,E) is factorized into P(Θ,E)P(X,E), and P(E _i |X,Θ,E\{E _i }) is given as a binomial distribution. Let \widehat{Q}\left(E_i\right)be Q(E _i ) maximizing the lower bound of the variational annealing for τ → +0 given by

\int d{E}_1...\int d{E}_p\int dX\int d\mathrm{\Theta }Q\left(X\right)Q\left(\mathrm{\Theta }\right)\prod _{i}Q\left(E_i\right)log\frac{P\left(X,\mathrm{\Theta },E\right)}{Q\left(X\right)Q\left(\mathrm{\Theta }\right){\prod }_{i}Q{\left(E_i\right)}^{\tau }}.

(5)

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.

Procedures of variational annealing on proposed model

In the variational annealing on the proposed model, we calculate Q functions for hidden variables X, parameters Θ, and binary variables E iteratively while cooling temperature τ to zero gradually at each iteration cycle. In the following, we show the calculation procedures of Q(X), Q(Θ), and Q(E) on the proposed model as variational E-step, variational M-step, and variational A-step, respectively. More details of the procedures are given in Additional file 1. For the notational brevity, we denote the expectation of a value x with a probability distribution Q(y) as 〈x〉_Q(y).

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

Q(Θ) is factorized into ∏ _i Q(A _i |h _i ) Q(h _i )Q(r _i ) ∏ _j Q(z _ij ), where A _i is a vector given by (A _{i 1}, ..., A _ip )'. From the design of the proposed model, Q(A _i |h _i ), Q(h _i ), Q(r _i ), and Q(z _ij ) are given in the following form:

\begin{array}{ll}\hfill Q\left(A_i|h_i\right)& =\mathcal{N}\left(A_i;{\mu }_{A_i},h_i{T}_{A_i}^{-1}\right),\\ \hfill Q\left(h_i\right)& =\mathcal{I}\mathcal{G}\left(h_i;u_i,k_i\right),\\ \hfill Q\left(r_i\right)& =\mathcal{I}\mathcal{G}\left(r_i;v_i,l_i\right),\\ \hfill Q\left(z_{ij}\right)& =\mathcal{B}\left({{\zeta }_{ij}}_{,0};{\zeta }_{ij,1}\right).\end{array}

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.

We first consider update of hyperparameters u₀, k₀, v₀, l₀, ζ_{i 0}, and ζ_{i 1}to increase the lower bound of marginal probability. u₀ and k₀ are updated by maximizing the following equation:

\begin{array}{ll}\hfill \left({û}_0,{\widehat{k}}_0\right)& =arg\underset{\left(u_0,k_0\right)}{max}\sum _i\int Q\left(h_i\right)log\mathcal{I}\mathcal{G}\left(h_i;u_0,k_0\right)d{h}_i\\ =arg\underset{\left(u_0,k_0\right)}{max}\left(u_0-1\right)\frac{{\sum }_i{〈log{h}_i〉}_{Q\left(h_i\right)}}{p}+u_{0}log{k}_0-k_0\frac{{\sum }_i{〈1∕h_i〉}_{Q\left(h_i\right)}}{p}-log\Gamma \left(u_0\right).\end{array}

{û}_0 and {\widehat{k}}_0 are obtained by numerical optimization methods such as the Newton-Raphson method. v₀ and l₀ are also updated in a similar manner to u₀ and k₀. ζ_{i 0}and ζ_{i 1}are updated by solving the following equation:

\begin{array}{ll}\hfill \left({\widehat{\zeta }}_{i0},{\widehat{\zeta }}_{i1}\right)& =arg\underset{{\zeta }_{i0},{\zeta }_{i1}}{max}\sum _j\int Q\left(z_{ij}\right)log\mathcal{B}\left(z_{ij};{\zeta }_{i0},{\zeta }_{i1}\right)d{z}_{ij}\\ =arg\underset{{\zeta }_{i0},{\zeta }_{i1}}{max}\left({\zeta }_{i1}-1\right)\frac{{\sum }_j{〈log\left(1-z_{ij}\right)〉}_{Q\left(z_{ij}\right)}}{p}+\left({\zeta }_{i0}-1\right)\frac{{\sum }_j{〈log{z}_{ij}〉}_{Q\left(z_{ij}\right)}}{p}-logB\left({\zeta }_{i0},{\zeta }_{i1}\right).\end{array}

{\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

The procedures for estimating parameters in the proposed model are summarized as follows:

• 1. Set τ to some large value. Also set α₀ to a small value and α₁ to a large value satisfying that α₀ <α₁.

• 2. Initialize other hyperparameters and hidden variables.

• 3. Perform the following procedures:

  • (a) Calculate variational M-step.

  • (b) Update hyperparameters.

  • (c) Calculate variational E-step.

  • (d) Calculate variational A-step.

  • (e) Go back to step (a) until some convergence criterion is satisfied.

• 4. Divide τ by some value > 1 such as 1.05.

• 5. Go back to step 3 if τ is larger than some very small value > 0.

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.

Results and discussion

Performance evaluation by Monte Carlo experiments

For the evaluation of the proposed approach, we generate two linear regulatory network models with similar topological structures G₁ and G₂ based on a linear regulatory network model G₀. G₀ is prepared in the following manner: (i) a scale free network of 100 nodes and 150 edges is generated; (ii) edge directions are assigned randomly; (iii) autoloop edges are added to root nodes of the directed network; and (iv) AR coefficients for the directed edges are chosen randomly from {-0.9, -0.8, -0.7, -0.6, -0.5, 0.5, 0.6, 0.7, 0.8, 0.9}. We then generate G₁ and G₂ from G₀ as follows: (i) autoloop edges and 70% of non-autoloop edges in G₀ are used for commonly existing edges in G₁ and G₂; and (ii) the other 30% of non-autoloop edges are randomly assigned as either G₁ or G₂ specific edges. Note that AR coefficients on edges of G₁ and G₂ are preserved, i.e., if a regulation from gene j to gene i exists in G₁ or G₂, then its coefficient is the same as that of the regulation from gene j to gene i in G₀.

Figure 2 gives graph structure of G₁ and G₂, where commonly regulations, G₁ specific regulations, and G₂ specific regulations are represented with black, red, and green arrows, respectively.

Figure 2

The graph structures of G ₁ and G ₂ . Commonly regulations, G₁ specific regulations, and G₂ specific regulations are represented with black, red, and green arrows, respectively.

From each of G₁ and G₂, we obtain two equally spaced time series data of 25 time points and 50 time points. For system noise and observation noise, normally distributed values with mean 0 and standard deviation 1 and mean 0 and standard deviation 0.1 and 1 are used, respectively. The signal-noise ratio for system noise with standard deviation 1 and observation noise with standard deviation 0.1 is 0.03dB, and the signal is bit stronger than the noise. On the other hand, for observation noise with standard deviation 1, the signal-noise ratio is -0.26dB. In the condition, the noise is stronger than the signal, and the noise level is quite high.

Comparison between variational annealing and EM algorithm

We first compare the performances of the proposed approach and the approach that is based on the proposed approach but uses the EM algorithm instead of the variation annealing using the equally spaced time series data of 50 and 25 time points on the system noise with standard deviation 1 and observation noise with standard deviation 0.1. From the comparison, we verify the effectiveness of the variational annealing, compared to the EM algorithm. Table 1 summarizes the results of the proposed approach and the EM algorithm based approach. For the EM algorithm, the regulation from j to i on condition c is considered to exist if the estimated E_{ij}^{\left(c\right)} is more than 0.5.

Table 1 Comparison of the variation annealing (Proposed) and EM algorithm (EM) based on the proposed model

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.

For the comparison of these approaches, we focus on the following two points: the number of correctly estimated regulations and the number of correctly estimated changes on regulations. The former is usually considered for evaluating the performance of gene network estimation methods. The numbers of true positives and false negatives of the estimated regulations are summarized in Table 2(a). The precisions of the results given by

Table 2 A summary of results for system noise with standard deviation 1 and observation noise with standard deviation 0.1

\frac{\mathsf{\text{The}}\mathsf{\text{number}}\mathsf{\text{of}}\mathsf{\text{true}}\mathsf{\text{positives}}}{\mathsf{\text{The}}\mathsf{\text{number}}\mathsf{\text{of}}\mathsf{\text{true}}\mathsf{\text{positives + The}}\mathsf{\text{number}}\mathsf{\text{of}}\mathsf{\text{false}}\mathsf{\text{positives}}}

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.

Figure 3

A time point schedule on unequally spaced time series data in the Monte Carlo experiment. Observed points in the time schedule are indicated by arrows. 15 time points are equally spaced in first block, every second point is observed in second block comprised of 5 observed time points, and every third point is observed in third block comprised of 5 observed time points.

We also consider the time series data with the high level noise: system noise with standard deviation 1 and observation noise with standard deviation 1. The results for the case are summarized in Tables 3(a) and 3(b). The proposed approach shows the better performance on the number of true positives and precisions than other approaches except for the identification of the changes on regulations from the equally spaced time series data of 50 time points. For equally spaced time series data of 50 time points, the number of true positives on the changes on regulations estimated by the proposed approach is more than that of G1DBN1, but less than that of ENet1. However, the precisions on the estimated changes by both ENet1 and G1DBN are much worse than the proposed approach. Thus, overall, the proposed approach is more effective than other methods. Although the proposed approach provides the better performance than other methods, the results of all the approaches are worsened due to the high level noise, and the differences on the performance among the approaches get smaller, compared to the case of observation noise with standard deviation 0.1.

Table 3 A summary of results for system noise with standard deviation 1 and observation noise with standard deviation 1

Analysis of time series microarray data from Human small airway epithelial cells

We apply the proposed approach to two time series microarray gene expression data from normal Human small airway epithelial cells (SAECs) and SAECs treated by stimulating EGF-receptors and dosing an anticancer drug termed Gefitinib. EGF-receptors are often overexpressed in lung cancer cells such as tumoral SAECs, and a lung cancer condition is simulated in the treated SAECs by stimulating EGF-receptors. Since Gefitinib is known as a selective inhibiter of EGF-receptors, the stimulation of EGF-receptors would be counteracted by Gefitinib, and the condition of treated SACEs should be the same as that of normal SAECs in theory. However, since some gene expression patterns are different between the two conditions in practice, some unknown effects by Gefitinib may be involved in the phenomenon. Thus, we focus on changed regulations between the gene networks estimated from gene expression data in these two conditions in order to find some insights on the unknown effects of Gefitinib.

For gene set selection, we first screen 500 genes from the ranking of the gene list sorted by coefficient variation [5]. We then select 100 genes with highly varied expression profiles between normal SAECs and treated SAECs. The time series gene expression data from both types of cells are comprised of spaced 14 time points in 48 hours. The time schedule of 14 time points are {0 h, 6 h, 9 h, 12 h, 15 h, 18 h, 21 h, 24 h, 27 h, 30 h, 33 h, 36 h, 39 h, 48 h}. For the analysis on the proposed approach, we set interval on system model to three hours.

The estimated networks from time series gene expression data in normal SAECs and treated SAECs are summarized and given in Figure 4. Black arrows, red arrows, and green arrows indicate regulations in both conditions, only in normal SAECs, and only in treated SAECs, respectively. Table 4 gives a list of the estimated regulations only in normal SAECs or in treated SAEC.

Figure 4

An estimated gene network by the proposed approach from time series gene expression data on normal SAECs and treated SAECs. In the estimated network, regulations in both conditions are in black, and regulations only in normal SAECs and treated SAECs are in red and green, respectively.

Table 4 Changes on regulations between normal and treated SAECs

From Table 4, we see that Prss22 is involved in most of the regulations only in treated SAECs. Prss22 is a tryptase, one of serine proteases, and its relationship with the airways is suggested by a report about its expression in the airways in a developmentally regulated manner. Tryptase is a potent mitogen of fibroblast [15], and it is reported that the increase and activation of fibroblast are promoted in lung cells under the condition of interstitial pneumonia. Interstitial pneumonia is known as a side effect of Gefitinib, and these findings suggest that Prss22 is an off-target of Gefitinib and is possibly related to the side effect of Gefitinib. Also, FILIP1L, a gene estimated as one of targets of Prss22 in treated SAECs, is reported as an inhibitor of cell proliferation of fibroblast [16]. This observation supports the relation between Prss22 and fibroblast as well.

We also focus on other several genes related to changes on regulations in normal and treated SAECs. LIF, leukemia inhibitory factor, is known to affect cell growth and development. Gefitinib is also known to be effective for acute myelogenous leukemia via Sky, which is an off-target gene of Gefitinib [17]. In addition, Wang et al. reported that LIF prolongs the cell cycle of stem cells on acute myelogenous leukemia lines [18]. Although, to the best of our knowledge, no direct influence from Gefitinib to LIF is reported, the above facts suggest some relation between Gefitinib and LIF, and support changes on regulations of LIF between normal and treated SAECs.

Heikema et al. reported that Human Siglecs, Siglecs-14, -15, and -16 interact with transmembrain adaptor proteins containing the immunoreceptor tyrosine-based activation motif such as DAP12 [19], and therefore they potentially mediate the activation of intracellular signaling. Gefitinib is a selective inhibitor of tyrosine kinase, and the inhibition of tyrosine kinase is considered to affect the regulations around Siglec-15.

Although the stimulation of EGF-receptors in the treated SAECs is considered to be counteracted by Gefitinib, the expressions of some genes may be affected by the stimulation in practical conditions. HAS3 is related to synthesis of the unbranched glycosaminoglycan hyaluronic acid and is reported to be up-regulated by EGF [20]. The stimulation of EGF-receptors affects the amount of EGF taken into cells, and hence the stimulation is expected to cause the changes on the regulations around HAS3. Foxn2 is a member of family of Fox proteins. Fox proteins are known to play important roles on controlling the expressions of genes related to cell growth, proliferation, and differentiation. Some members of Fox family are related to EFG-receptors, e.g., the expression of Foxn1 is suppressed by EGF-receptor signaling [21] although no direct relation between EGF-receptors and Foxn2 is found. FOXA2 is also a member of family of Fox proteins. EGF-receptor signaling is known to decrease the expression of FOXA, which prevents the mucus production [22]. [23] reported the relation between EGF-receptors and FOXA2 in the airways of asthmatic patients. Thus, it appears that the change on the regulation related to FOXA2 is caused by the stimulation of EGF-receptors.

Conclusions

We proposed the new computational model that is based on VAR-SSM and estimates gene networks from time series data on normal and treated conditions as well as identifies changes regulations by the treatment. Unlike many of existing gene network estimation approaches assuming equally spaced time points, our approach can handle unequally spaced time series data. The efficient use of time series data is achieved by representing the presence of regulations on each condition with hidden binary variables. Since finding the optimal configuration of the hidden binary variables on the proposed model is computationally in tractable, we derive the extended variational annealing method in order to address the problem as the alternative method.

In the Monte Carlo experiments, we use equally and spaced time series data from synthetically generated two regulatory networks whose structures are different in several regulations, and verified the effectiveness of the proposed model in both estimation of regulations and changes on regulations between the two conditions, compared to existing methods.

As the real data application, we use the proposed approach to analyze two time series data from normal SAECs and SAECs treated by stimulating EGF-receptors and dosing Gefitinib. From genes related to changes on regulations by the treatment, we find possible off-target genes of Gefitinib, and one of these genes is suggested to be related to a factor of interstitial pneumonia, which is known as a side effect of Gefitinib. In this study, we consider changes on regulations in two conditions, but the proposed approach can be extended to identifying changes among more than two conditions.