Identifying dynamical time series model parameters from equilibrium samples,with application to gene regulatory networks

Abstract

Gene regulatory network reconstruction is an essential task of genomics in order to further our understanding of how genes interact dynamically with each other. The most readily available data, however, are from steady-state observations. These data are not as informative about the relational dynamics between genes as knockout or over-expression experiments, which attempt to control the expression of individual genes. We develop a new framework for network inference using samples from the equilibrium distribution of a vector autoregressive (VAR) time-series model which can be applied to steady-state gene expression data. We explore the theoretical aspects of our method and apply the method to synthetic gene expression data generated using GeneNetWeaver.

1 Introduction

With continuing advances in gene expression measurement technologies, more and more expression data are becoming available for analysis. One of the important uses for these data has been in inferring gene regulatory networks. These networks identify relationships between genes and give insight into the complex inner workings of the cell.

Even with the increasing amounts of data, constructing gene regulatory networks is difficult in practice. The large number of genes involved, generally in the thousands, makes methods developed for smaller networks computationally infeasible. Time-series genetic data, with the same population of cells observed at multiple times, is useful in inferring the direction of the relationships (Yeung et al., 2011; Young et al., 2014). Knockdown and over-expression data can also be used to infer directionality of relationships (Young et al., 2016). However, most data take the form of perturbation screens or steady-state data at a single timepoint. This makes it harder to identify the direction of the relationship between two genes without prior knowledge.

There have been many methods developed for the analysis of steady-state data. Mutual information and correlation-based methods are commonly used, although the resulting networks tend to lack directionality (Basso et al., 2005; Faith et al., 2007; Margolin et al., 2006; Meyer et al., 2007; Tusher et al., 2001). Regression-based methods have also been applied to steady-state data to infer network structure (Omranian et al., 2016; Singh and Vidyasagar, 2016). Bayesian network methods (Friedman et al., 2010) can result in directed graphs, but the networks generated are, of necessity, acyclic. This means that any cycles, such as those that are found in real biological systems, will not be captured, so the resulting networks will be unrealistic in this respect. Shojaie et al. (2014) developed a method combining information from both perturbation and steady-state data to increase the accuracy of the inferred network.

Here, we propose a new method based on an implicit multivariate time-series model that uses steady-state data and has the capability to test the existence of edges in directed networks. This method is not restricted to acyclic networks, but can be used to infer information about cycles in the network. We present proofs of consistency and efficiency in a constrained case as well as a likelihood ratio test for the existence of an edge. We also give simulation results and apply the method to a synthetic dataset, showing that the method performs well in practice.

2 First-order vector autoregressive model

To build a model for equilibrium data, first consider a first-order vector autoregres- sive, or VAR(1) model for time-series data for a system with p genes with no correlation in the error term between genes (Lütkepohl, 2005). We can write the model as

$$\begin{array}{ccc}\hfill {\mathbf{x}}_1& \hfill =\hfill & {\mathbf{\epsilon }}_1,\hfill \\ \hfill {\mathbf{x}}_t& \hfill =\hfill & \mathbf{A}{\mathbf{x}}_{t-1}+{\mathbf{\epsilon }}_t\mathit{for}t>1,\hfill \\ \hfill {\mathbf{\epsilon }}_t& \hfill \sim \hfill & N(0,\mathbf{D}),\hfill \end{array}$$

where $\mathbf{D}$ is diagonal, ${\mathbf{x}}_t$ and ${\mathbf{\epsilon }}_t$ are vectors of length $p$ , while $\mathbf{A}$ and $\mathbf{D}$ are $p\times p$ matrices. Here, $\mathbf{A}$ identifies the relationships between the genes. A non-zero $(i,j)$ -th element of $\mathbf{A}$ indicates that there is an edge from gene $j$ to gene $i$ .

Autoregressive models have been broadly used to infer gene networks from time-series data, with approaches ranging from penalized regression (LASSO, elastic net) to Bayesian model averaging (BMA) to dynamic Bayesian networks (DBN). Michailidis and dAlché Buc (2013) review many of these approaches to modelling time series using autoregressive models.

At first glance, the VAR(1) model does not appear to be applicable to steady-state data. However, if the eigenvalues of $\mathbf{A}$ are less than 1 in absolute value, then as $t\to \infty$ , ${\mathbf{x}}_t$ converges to a stable equilibrium distribution (Anderson, 2000). We can use this equilibrium distribution as the model for steady-state data. This can be applied to either a set of experiments with the same perturbation or to the control experiments.

2.1 Equilibrium distribution

Since this is a Gaussian VAR model, we can write the equilibrium distribution as

$${\mathbf{x}}_{\infty }\sim N(0,\mathbf{\Sigma }),$$

where $\mathbf{\Sigma }$ is also a $p\times p$ matrix. To find $\mathbf{\Sigma }$ , we use an iterative relationship between the variances at consecutive time points:

$$\begin{array}{ccc}\hfill \mathit{Var}[{\mathbf{x}}_1]& \hfill =\hfill & \mathbf{D},\hfill \\ \hfill \mathit{Var}[{\mathbf{x}}_{t+1}]& \hfill =\hfill & \mathbf{A}\mathit{Var}[{\mathbf{x}}_t]{\mathbf{A}}^T+\mathbf{D},\hfill \\ \hfill \mathbf{\Sigma }& \hfill =\hfill & \sum _{i=0}^{\infty }{\mathbf{A}}^i\mathbf{D}({\mathbf{A}}^T{)}^i.\hfill \end{array}$$

We will call this the summation identity (Sheppard, 2013).

Another expression for the asymptotic variance uses the identity $\mathit{vec}(\mathbf{A}\mathbf{X}\mathbf{B})=({\mathbf{B}}^T\otimes \mathbf{A})\mathit{vec}(\mathbf{X})$ :

$$\begin{array}{ccc}\hfill \mathbf{\Sigma }& \hfill =\hfill & \sum _{i=0}^{\infty }{\mathbf{A}}^i\mathbf{D}({\mathbf{A}}^T{)}^i,\hfill \\ \hfill \mathit{vec}(\mathbf{\Sigma })& \hfill =\hfill & \mathit{vec}\left(\sum _{i=0}^{\infty }{\mathbf{A}}^i\mathbf{D}({\mathbf{A}}^T{)}^i\right),\hfill \\ \hfill & \hfill =\hfill & \sum _{i=0}^{\infty }\mathit{vec}({\mathbf{A}}^i\mathbf{D}({\mathbf{A}}^T{)}^i),\hfill \\ \hfill & \hfill =\hfill & \sum _{i=0}^{\infty }\left({\mathbf{A}}^i\otimes {\mathbf{A}}^i\right)\mathit{vec}(\mathbf{D}),\hfill \\ \hfill & \hfill =\hfill & \sum _{i=0}^{\infty }{\left(\mathbf{A}\otimes \mathbf{A}\right)}^i\mathit{vec}(\mathbf{D}),\hfill \\ \hfill & \hfill =\hfill & {\left(\mathbf{I}-\mathbf{A}\otimes \mathbf{A}\right)}^{-1}\mathit{vec}(\mathbf{D}).\hfill \end{array}$$

This Kronecker identity is convenient for calculating the asymptotic variance from a given $\mathbf{A}$ and $\mathbf{D}$ since it eliminates the infinite sum (Sheppard, 2013). We can also write out a recursive identity for $\mathbf{\Sigma }$ since, at the equilibrium, the variance does not change from one time step to the next:

$$\mathbf{\Sigma }=\mathbf{A}\mathbf{\Sigma }{\mathbf{A}}^T+\mathbf{D}.$$

(2.1)

3 Identifiability

When looking at steady-state data as opposed to time-series data, we must rely on the relationship between $(\mathbf{A},\mathbf{D})$ and $\mathbf{\Sigma }$ in order to infer anything about $\mathbf{A}$ and $\mathbf{D}$ . This raises an issue of identifiability. As a simple example, substituting $-\mathbf{A}$ for $\mathbf{A}$ in any of the earlier variance identities does not change anything, so the sign of $\mathbf{A}$ cannot be identified from equilibrium data only. More generally, let $\mathbf{Q}$ be any $p\times p$ orthogonal matrix. An orthogonal matrix is defined as any matrix such that $\mathbf{Q}{\mathbf{Q}}^T=\mathbf{I}$ . For any such $\mathbf{Q}$ , if $\tilde{\mathbf{A}}$ is defined as

$$\tilde{\mathbf{A}}=\mathbf{A}{\mathbf{\Sigma }}^{1/2}\mathbf{Q}{\mathbf{\Sigma }}^{-1/2}$$

if will satisfy the recursive identity (Tong et al., 1992).

We can state the identifiability problem as follows. Under what conditions is the map

$$(\mathbf{A},\mathbf{D})⟷\mathbf{\Sigma }$$

bijective? As a necessary condition we can count parameters.

Proposition A necessary condition for identifiability in the VAR equilibrium model of dimension $p$ , with connected nodes and non-zero diagonal elements of $\mathbf{A}$ , is that the number of non-zero off-diagonal elements of $\mathbf{A}$ be no more than $p(p-3)/2$ and that $p\ge 5$ .

Proof The number of parameters in $\mathbf{\Sigma }$ for a $p$ -dimensional model is $p(p+1)/2$ . For identifiability, the total number of independently estimable parameters in $\mathbf{A}$ and $\mathbf{D}$ must not exceed that number. $\mathbf{D}$ is diagonal and thus accounts for $p$ parameters. Similarly, the diagonal of $\mathbf{A}$ is assumed to be non-zero, implying that the expression level of each gene at timepoint $t+1$ is dependent on the expression level of that gene at timepoint $t$ for all genes. This accounts for an additional $p$ parameters. Thus, a maximum $p(p-3)/2$ off-diagonal elements of $\mathbf{A}$ may be non-zero. Now, the graph defined by $\mathbf{A}$ must be connected. If it is not, then it can be decomposed into two smaller, independent components, and each of those must be identifiable. We need at least $p-1$ off-diagonal elements of $\mathbf{A}$ for it to be connected. Thus, $p$ must be at least 5, because if $p=4$ , only a maximum of 2 off-diagonal elements is allowed, which is not enough to fully connect the 4 nodes.

Parameter counting yields a necessary condition, but it is not a sufficient condition. We have not found sufficient conditions in the most general case. We have, however, been able to examine identifiability in two specific examples, one acyclic and one cyclic. To assess the identifiability of a specific model, we used the recursive identity,

$$\mathbf{\Sigma }=\mathbf{A}\mathbf{\Sigma }{\mathbf{A}}^T+\mathbf{D}.$$

If a given model is identifiable, then there will be only one $(\mathbf{A},\mathbf{D})$ pair which satisfies this equation.

The first example is a five-dimensional model with a simple network structure, given in Figure 1. This model has four off-diagonal non-zero elements in $\mathbf{A}$ . The maximum allowable by parameter counting is 5, so we are close to but not at that limit. Writing out the recursive identity as a system of equations yields equations, the first two of which are:

$$\begin{array}{ccc}\hfill {\sigma }_{11}& \hfill =\hfill & a_{11}^2{\sigma }_{11}+2a_{11}{a}_{12}{\sigma }_{12}+a_{22}^2{\sigma }_{22}+d_1,\hfill \\ \hfill {\sigma }_{12}& \hfill =\hfill & a_{11}{a}_{22}{\sigma }_{12}+a_{12}{a}_{22}{\sigma }_{22}+a_{11}{a}_{23}{\sigma }_{13}+a_{12}{a}_{23}{\sigma }_{23},\hfill \\ \hfill \dots & \hfill \hfill & \hfill \end{array}$$

where $a_{\mathit{ij}}$ is the entry of $\mathbf{A}$ in the $i$ ’th row and $j$ ’th column, and $d_i$ is the $i$ -th diagonal entry in $\mathbf{D}$ . Next, we solve the system of equations to eliminate as many variables as possible. Then we perform a grid search on the free variables to find all potential solutions of the recursive identity. In this example, we were able to solve the system down to a single free parameter, $a_{55}$ .

To test the parameter solutions, we used the following model parameters:

$$\begin{array}{ccc}\hfill \mathbf{A}& \hfill =\hfill & \left[\begin{array}{ccccc}\hfill 0.6\hfill & \hfill 0.3\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0.4\hfill & \hfill -0.4\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0.8\hfill & \hfill 0.5\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.5\hfill & \hfill -0.3\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.7\hfill \end{array}\right],\hfill \\ \hfill \mathbf{D}& \hfill =\hfill & {\mathbf{I}}_5.\hfill \end{array}$$

The grid search was performed by allowing $a_{55}$ to vary between 0 and 1 in increments of 0.0001. This is the allowable range for $a_{55}$ since $\mathbf{A}$ and $-\mathbf{A}$ both result in the same $\mathbf{\Sigma }$ . For each prospective value of $a_{55}$ , the other elements of $\mathbf{A}$ and $\mathbf{D}$ were calculated. These values were then used to evaluate the objective function $\text{log}||\mathbf{\Sigma }-\mathbf{A}\mathbf{\Sigma }{\mathbf{A}}^T-\mathbf{D}|{|}_2^2$ . Figure 2 shows the results. From this it is clear that the only value of $a_{55}$ that satisfies the recursive identity is the true value, 0.7, where the objective function is minimized. The flat region of the objective curve to the left of $a_{55}\approx 0.65$ reflects the fact that there is no valid solution for $\mathbf{A}$ and $\mathbf{D}$ with that starting value for $a_{55}$ , either because some element of $\mathbf{A}$ is not a real value or because a diagonal element of $\mathbf{D}$ is negative. Also note that there is a degenerate solution at $a_{55}=1$ , such that $\mathbf{A}=\mathbf{I}$ and $\mathbf{D}=0$ . Since $\mathbf{D}$ represents the variance of the noise added to the VAR system, its diagonal elements must be positive. So in this example the model is identifiable.

Figure 1:

The network structure for the cyclic identifiability example

Figure 2:

Grid search for solutions of the acyclic identifiability example. The objective function is minimized at the true value of $a_{55}=0.7$ , showing that the model is identifiable in this example

Our second example involves a cyclic network. We added a single edge to the previous network, making a 3-node cycle, as shown in Figure 3. We used the following model parameters:

$$\begin{array}{ccc}\hfill \mathbf{A}& \hfill =\hfill & \left[\begin{array}{ccccc}\hfill 0.3\hfill & \hfill 0.6\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0.1\hfill & \hfill -0.3\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0.7\hfill & \hfill 0.2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0.3\hfill & \hfill 0\hfill & \hfill 0.4\hfill & \hfill -0.6\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.5\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \end{array}\right],\hfill \\ \hfill \mathbf{D}& \hfill =\hfill & {\mathbf{I}}_5.\hfill \end{array}$$

For this model, we were able to solve down to two free parameters, $a_{45}$ and $a_{55}$ . We again performed a grid search over the appropriate parameter space ( $a_{45}\in (-1,1)$ , $a_{55}\in (0,1)$ ) with an initial grid increment of 0.001 and evaluated the objective function at each point. Figure 4 shows the results of the grid search, both in the larger space as well as zoomed in around the truth, where a much smaller grid increment of ${10}^{-11}$ was used. There are large areas where a valid solution is not achieved, indicated by the white areas. Again, however, the original solution is the only valid one, indicating that this model is also identifiable.

Figure 3:

The network structure for the cyclic identifiability example

Figure 4:

Grid search for solutions of the cyclic identifiability example showing the value of the objective function at each point searched. The left plot shows the grid search over the entire model space, while the right plot is zoomed in on the true parameter value. The white areas indicate parameter values that result in invalid solutions. The true solution at ( $-$ 0.6, 0.5) is the only valid solution

3.1 Simple case

Although the parameter counting condition is helpful and results in identifiability in at least a couple of examples, it is important to find a more general result for identifiability.

We now find a more general condition for identifiablity within a specific subset of VAR(1) models. We do this by restricting the model considerably to make the analysis more tractable. First, we constrain $\mathbf{D}$ to be known. Second, we require that the model be acyclic, and thus the nodes can be ordered such that $\mathbf{A}$ is lower-triangular. Finally, we require that the diagonal entries of $\mathbf{A}$ be less than 1 in absolute value and that their signs be known. With these constraints, the model is identifiable and has good asymptotic properties.

Consider the model as constrained earlier. We order the nodes in the network topologically, meaning that there are no edges from node $i$ to node $j$ , if $i>j$ . This creates an $\mathbf{A}$ that is lower-triangular. We will use the recursive identity for the equilibrium variance:

$$\mathbf{A}\mathbf{\Sigma }{\mathbf{A}}^T+\mathbf{D}=\mathbf{\Sigma },$$

(3.1)

$$\mathbf{A}\mathbf{\Sigma }{\mathbf{A}}^T=\mathbf{\Sigma }-\mathbf{D}.$$

(3.2)

We solve this identity recursively by looking at the leading principal submatrices, an approach inspired by Drton et al. (2011). First, define ${\mathbf{S}}_i$ to be the leading $i\times i$ principal submatrix of $\mathbf{\Sigma }$ . That is the matrix made of elements in the first $i$ rows and columns of $\mathbf{\Sigma }$ . Similarly, define ${\mathbf{B}}_i$ to be the leading principal submatrix of $\mathbf{B}=\mathbf{\Sigma }-\mathbf{D}$ and ${\mathbf{\Gamma }}_i$ to be the leading principal submatrix of $\mathbf{A}$ . Also, define ${\mathbf{s}}_{i+1}$ to be the vector containing the first $i$ elements of the $(i+1{)}^{\mathit{th}}$ column of $\mathbf{\Sigma }$ or $\mathbf{B}$ . That is, ${\mathbf{s}}_{i+1}=[{\sigma }_{1,i+1},\dots ,{\sigma }_{i,i+1}{]}^{\prime }$ . Similarly, ${\mathbf{\gamma }}_{i+1}^T$ is the vector containing the first $i$ elements of the $(i+1{)}^{\mathit{th}}$ row of $\mathbf{\Gamma }$ . So,

$$\begin{array}{ccc}\hfill {\mathbf{S}}_{i+1}& \hfill =\hfill & \left[\begin{array}{cc}\hfill {\mathbf{S}}_i\hfill & \hfill {\mathbf{s}}_{i+1}\hfill \\ \hfill {\mathbf{s}}_{i+1}^T\hfill & \hfill s_{i+1,i+1}\hfill \end{array}\right],s_{\mathit{ii}}={\sigma }_{\mathit{ii}}\hfill \\ \hfill {\mathbf{B}}_{i+1}& \hfill =\hfill & \left[\begin{array}{cc}\hfill {\mathbf{B}}_i\hfill & \hfill {\mathbf{s}}_{i+1}\hfill \\ \hfill {\mathbf{s}}_{i+1}^T\hfill & \hfill b_{i+1,i+1}\hfill \end{array}\right],b_{\mathit{ii}}={\sigma }_{\mathit{ii}}-d_{\mathit{ii}}\hfill \\ \hfill {\mathbf{\Gamma }}_{i+1}& \hfill =\hfill & \left[\begin{array}{cc}\hfill {\mathbf{\Gamma }}_i\hfill & \hfill 0\hfill \\ \hfill {\mathbf{\gamma }}_{i+1}^T\hfill & \hfill {\gamma }_{i+1,i+1}\hfill \end{array}\right],{\gamma }_{\mathit{ii}}=a_{\mathit{ii}}.\hfill \end{array}$$

Now, for $i=1$ , we get

$$\begin{array}{ccc}\hfill {\gamma }_{\mathit{ii}}^2{s}_{\mathit{ii}}& \hfill =\hfill & b_{\mathit{ii}},\hfill \\ \hfill a_{\mathit{ii}}& \hfill =\hfill & \sqrt{\frac{{\sigma }_{\mathit{ii}}-d_{\mathit{ii}}}{{\sigma }_{\mathit{ii}}}}.\hfill \end{array}$$

This gives us ${\mathbf{A}}_1$ . Now, suppose we know ${\mathbf{A}}_i$ . Then, for $i+1$ , we have

$$\begin{array}{ccc}\hfill {\mathbf{B}}_{i+1}& \hfill =\hfill & {\mathbf{\Gamma }}_{i+1}{\mathbf{S}}_{i+1}{\mathbf{\Gamma }}_{i+1}^T,\hfill \\ \hfill \left[\begin{array}{cc}\hfill {\mathbf{B}}_i\hfill & \hfill {\mathbf{s}}_{i+1}\hfill \\ \hfill {\mathbf{s}}_{i+1}^T\hfill & \hfill b_{i+1,i+1}\hfill \end{array}\right]& \hfill =\hfill & \left[\begin{array}{cc}\hfill {\mathbf{\Gamma }}_i{\mathbf{S}}_i{\mathbf{\Gamma }}_i^T\hfill & \hfill {\mathbf{\Gamma }}_i{\mathbf{S}}_i{\mathbf{\gamma }}_{i+1}+{\gamma }_{i+1,i+1}{\mathbf{\Gamma }}_i{\mathbf{s}}_{i+1}\hfill \\ \hfill -\hfill & \hfill {\mathbf{\gamma }}_{i+1}^T{\mathbf{S}}_i{\mathbf{\gamma }}_{i+1}+2{\gamma }_{i+1,i+1}{\mathbf{s}}_{i+1}^T{\mathbf{\gamma }}_{i+1}+{\gamma }_{i+1,i+1}^2{s}_{i+1,i+1}\hfill \end{array}\right].\hfill \end{array}$$

We solve for ${\gamma }_{i+1,i+1}$ and ${\mathbf{\gamma }}_{i+1}$ to get

$$\begin{array}{ccc}\hfill {\gamma }_{i+1,i+1}& \hfill =\hfill & \sqrt{\frac{b_{i+1,i+1}-{\mathbf{s}}_i^T({\mathbf{\Gamma }}_i{\mathbf{S}}_i{\mathbf{\Gamma }}_i^T{)}^{-1}{\mathbf{s}}_i}{s_{i+1,i+1}-{\mathbf{s}}_i^T{\mathbf{\Gamma }}^{-1}{\mathbf{s}}_i}},\hfill \\ \hfill {\mathbf{\gamma }}_{i+1}& \hfill =\hfill & ({\mathbf{\Gamma }}_i{\mathbf{S}}_i{)}^{-1}(\mathbf{I}-{\gamma }_{i+1,i+1}{\mathbf{\Gamma }}_i){\mathbf{s}}_i.\hfill \end{array}$$

Notice here that the expression for ${\gamma }_{i+1,i+1}$ has two solutions, the positive and negative square root. These elements are the diagonal entries in $\mathbf{A}$ , thus requiring that the signs of the diagonal elements be specified. If they are not specified, there are $2^p$ possible solutions. Once the signs of these elements are specified, we have a unique solution for $\mathbf{A}$ given $\mathbf{\Sigma }$ and $\mathbf{D}$ . However, there still exists positive definite $\mathbf{\Sigma }$ that do not map back to a valid $\mathbf{A}$ given $\mathbf{D}$ . As a simple example, if ${\sigma }_{\mathit{ii}}<d_{\mathit{ii}}$ , then the model will not work since that fact means that the equilibrium variance is less than the error variance added at each step. Now we show that the maximum likelihood estimator (MLE) for $\mathbf{A}$ in this model is both consistent and efficient:

Theorem Suppose $p$ -dimensional ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_n$ are independent and are observed from the distribution

$$X|\mathbf{A}\sim N\left(0,\sum _{i=0}^{\infty }{\mathbf{A}}^i\mathbf{D}({\mathbf{A}}^T{)}^i\right),$$

where $\mathbf{D}$ is known and $\mathbf{A}$ is lower-triangular with the signs of the diagonal elements known, and the diagonal elements less than 1 in absolute value. Then $\widehat{\mathbf{A}}$ , the MLE of $\mathbf{A}$ , is consistent and asymptotically efficient with

$$\sqrt{n}(\widehat{\mathbf{a}}-\mathbf{a})\to N(0,\mathcal{F}(\mathbf{a}{)}^{-1}),$$

where $\mathbf{a}$ is the vectorization of the non-zero elements of $\mathbf{A}$ and $\mathcal{F}(\mathbf{a})$ is the Fisher information matrix, whose elements are

$$\mathcal{F}(\mathbf{a}{)}_{\mathit{ij}}=E\left[\left(\frac{\partial }{\partial a_i}\text{log}f(\text{x}|\text{A})\right)\left(\frac{\partial }{\partial a_j}\text{log}f(\text{x}|\text{A})\right)\right].$$

Proof We will show that the conditions for Cramér's theorem, as found in Ferguson (1996, p. 121), are satisfied.

The appendix contains derivations of the first- and second-partial derivatives and show that they exist and are continuous. For example, the first partial derivatives take the form

$$\begin{array}{ccc}\hfill \frac{\partial {\sigma }_{\mathit{ij}}}{\partial a_{\mathit{mn}}}& \hfill =\hfill & \frac{1}{1-a_{\mathit{ii}}{a}_{\mathit{jj}}}\left[\begin{array}{c}\hfill (\text{A}\mathbf{\Sigma }{)}_{\mathit{in}}{1}_{[j=m,j\ge n]}+(\text{A}\mathbf{\Sigma }{)}_{\mathit{jn}}{1}_{[i=m,i\ge n]}\hfill \\ \hfill +\sum _{k=1}^i\sum _{l=1}^j{a}_{\mathit{ik}}{a}_{\mathit{jl}}\frac{\partial {\sigma }_{\mathit{kl}}}{\partial a_{\mathit{mn}}}{1}_{[(k,l)\ne (i,j)]}\hfill \end{array}\right].\hfill \end{array}$$

Both the first- and second-partial derivatives can be calculated iteratively.

Now let

$$\begin{array}{ccc}\hfill \mathbf{\Psi }(x,\mathbf{a})=\frac{\partial }{\partial \mathbf{a}}\text{log}f(x|\mathbf{a}),& \hfill \hfill & \mathit{a}k\mathit{vector},\hfill \\ \hfill \dot{\mathbf{\Psi }}(x,\mathbf{a})=\frac{{\partial }^2}{\partial {\mathbf{a}}^2}\text{log}f(x|\mathbf{a}),& \hfill \hfill & \mathit{a}k\times k\mathit{matrix}.\hfill \end{array}$$

Note that the components of $\dot{\mathbf{\Psi }}(x,\mathbf{a})$ , as shown in the appendix, are all of the form

$$\sum _{i,j}{c}_{\mathit{ij}}(\mathbf{x}{\mathbf{x}}^T{)}_{\mathit{ij}},$$

where the $c_{\mathit{ij}}$ ’s are constant. Further, the expectations of the absolute values of the components are finite since $E[|\mathbf{x}{\mathbf{x}}^T{|}_{\mathit{ij}}]<\infty$ (Li and Wei, 2012). Now lets

$$K(\mathbf{x})=\sum _{i,j}{m}_{\mathit{ij}}|\mathbf{x}{\mathbf{x}}^T{|}_{\mathit{ij}},$$

where $m_{\mathit{ij}}$ is the maximum of all the $|c_{\mathit{ij}}|$ ’s across all components. Then each component of $\dot{\mathbf{\Psi }}(x,\mathbf{a})$ is bounded by $K(\mathbf{x})$ in absolute value and $K(\mathbf{x})$ has finite expectation.

These results, along with the earlier proof of identifiability and the form of the constraints on $\mathbf{A}$ , satisfy Cramér's theorem. Thus, the MLE of $\mathbf{A}$ is asymptotically consistent and efficient.

In addition to asymptotic consistency and efficiency, we also can form a likelihood-ratio test for the existence of an edge.

Corollary Let ${\mathcal{A}}_1$ be the set of all valid $\mathbf{A}$ in the VAR equilibrium model, possibly with some entries constrained to be 0. Let ${\mathcal{A}}_0\subset {\mathcal{A}}_1$ with a difference in dimension between the two spaces of $c$ . This corresponds to $c$ additional entries of $\mathbf{A}$ constrained to be 0. Then the log-likelihood ratio test statistic

$$\lambda =-2\text{log}\left(\frac{{\text{sup}}_{\text{A}\in {\mathcal{A}}_0}L(\text{A}|\text{x})}{{\text{sup}}_{\text{A}\in {\mathcal{A}}_1}L(\text{A}|\text{x})}\right)$$

has a ${\chi }_c^2$ distribution, asymptotically.

Proof The conditions which satisfy Cramér's theorem also satisfy Wilk's theorem concerning the likelihood-ratio test statistic (Ferguson, 1996, p. 145).

3.2 Simulation results: Asymptotic variance

We have an asymptotic distribution for the MLE, and now we assess the quality of its approximation in finite samples via simulation. Because of the complexity of the calculations involved in finding the Fisher information matrix, we will restrict ourselves to a two-dimensional model. This model has three parameters for $\mathbf{A}$ ,

$$\mathbf{A}=\left[\begin{array}{cc}\hfill a_{11}\hfill & \hfill 0\hfill \\ \hfill a_{21}\hfill & \hfill a_{22}\hfill \end{array}\right],$$

so $\mathbf{a}=[a_{11},a_{21},a_{22}{]}^T$ and the Fisher information matrix will be 3 by 3, with elements

$$\mathcal{F}(\mathbf{a}{)}_{\mathit{ij}}=E\left[\left(\frac{\partial }{\partial a_i}\text{log}f(\text{x}|\text{A})\right)\left(\frac{\partial }{\partial a_j}\text{log}f(\text{x}|\text{A})\right)\right].$$

For our simulations, we used

$$\begin{array}{ccc}\hfill \mathbf{A}& \hfill =\hfill & \left[\begin{array}{cc}\hfill 0.9\hfill & \hfill 0\hfill \\ \hfill 0.7\hfill & \hfill 0.8\hfill \end{array}\right],\hfill \\ \hfill \text{D}& \hfill =\hfill & {\text{I}}_2.\hfill \end{array}$$

Using these values, the asymptotic variance of the MLE for $\mathbf{A}$ is calculated to be

$$\mathcal{F}(\mathbf{a}{)}^{-1}=\left[\begin{array}{ccc}\hfill 0.022\hfill & \hfill -0.126\hfill & \hfill 0.029\hfill \\ \hfill -0.126\hfill & \hfill 3.094\hfill & \hfill -0.850\hfill \\ \hfill 0.029\hfill & \hfill -0.850\hfill & \hfill 0.251\hfill \end{array}\right].$$

To test this, we performed a simulation where we generated $n$ samples from the equilibrium distribution. We then found the value of $\widehat{\text{A}}$ that maximized the likelihood. To do this, we initialized $\mathbf{A}$ using the recursive algorithm described earlier and used the function optim in R to optimize the log-likelihood of the model, utilizing the BFGS method developed by Broyden (1970); Fletcher (1970); Goldfarb (1970). We repeated the initialization and optimization procedure 1 000 times and used the resulting $\widehat{\text{A}}$ values to calculate an empirical variance. We used three different sample sizes $n$ , resulting in the following empirical variances:

$$\begin{array}{cc}n=1000,& \\ & Var\left(\widehat{\mathbf{A}}\right)=\left[\begin{array}{ccc}0.021& -0.077& 0.009\\ -0.077& 2.633& -0.717\\ 0.009& -0.717& 0.219\end{array}\right]\end{array},$$

$$\begin{array}{cc}n=10000,& \\ & Var\left(\widehat{\mathbf{A}}\right)=\left[\begin{array}{ccc}0.023& -0.074& 0.007\\ -0.074& 2.705& -0.746\\ 0.007& -0.746& 0.230\end{array}\right]\end{array},$$

$$\begin{array}{cc}n=100000,& \\ & Var\left(\widehat{\mathbf{A}}\right)=\left[\begin{array}{ccc}0.023& -0.091& 0.012\\ -0.091& 2.757& -0.737\\ 0.012& -0.737& 0.220\end{array}\right]\end{array},$$

These results show that the asymptotic variance approximates the finite sample variances quite well, and so provide support for the asymptotic theory.

3.3 Simulation results: Coverage

Another way to assess the accuracy of the asymptotic approximation for finite samples is to look at the coverage of the resulting confidence intervals. To do this, we used the same two-dimensional model as earlier. In each single simulation, $a_{11}$ and $a_{22}$ were chosen at random from a Uniform(0.1, 0.9) distribution and $a_{21}$ was chosen at random from a Uniform( $-$ 0.9, 0.9) distribution. We then generated $n$ samples from the equilibrium distribution and used those samples to obtain an MLE for $\text{A}$ by maximizing the log-likelihood as before. We then used the empirical asymptotic variance of the MLE to obtain confidence intervals and see if the true values are covered. We repeat this 1 000 times with the same $\text{A}$ to get coverage results.

We carried out this procedure 10 times each for $n$ equal to 100, 1 000 and 10 000, generating a different $\text{A}$ each time. The coverage results are shown in Figures 5–8. The x-axis is the true value of $a_{11}$ . This value matters more than the other $\text{A}$ values because it propagates to all the values of $\mathbf{\Sigma }$ . As $n$ increases, the coverage results clearly improve. The coverage is essentially correct for large $n$ , both marginally and jointly.

We note that with small $n$ and $a_{11}$ , the observed variance of $x_1$ can be smaller than that of the known $\text{D}$ . This can cause problems with optimizing $\text{A}$ .

Figure 5:

Joint coverage for $\mathbf{A}$ using the asymptotic distribution of the MLE]Joint coverage for all parameters of $\text{A}$ using the 95% confidence interval from the multivariate normal asymptotic distribution of the MLE

Figure 6:

Marginal coverage for $a_{11}$ using the asymptotic distribution of the MLE]Marginal coverage for $a_{11}$ using the 95% confidence interval from the asymptotic distribution of the MLE

Figure 7:

Marginal coverage for $a_{21}$ using the asymptotic distribution of the MLE]Marginal coverage for $a_{21}$ using the 95% confidence interval from the asymptotic distribution of the MLE

Figure 8:

Marginal coverage for $a_{22}$ using the asymptotic distribution of the MLE]Marginal coverage for $a_{22}$ using the 95% confidence interval from the asymptotic distribution of the MLE

4 Likelihood ratio test

We now assess the validity of the likelihood-ratio test for the simple VAR equilibrium model. To test this, we used a three-dimensional system with

$$\text{A}=\left[\begin{array}{ccc}\hfill 0.6\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -0.3\hfill & \hfill 0.4\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -0.5\hfill & \hfill 0.8\hfill \end{array}\right].$$

For this system, ${\mathcal{A}}_0$ has the correct constraint that the lower-left element is zero, while ${\mathcal{A}}_1$ does not have that constraint. For a dataset generated from the model, we find the MLE under both hypotheses and construct the likelihood ratio test statistic.

To obtain the empirical distribution of the test statistic, we generated 10 000 samples from the equilibrium distribution and calculated the likelihood ratio for those data. We repeated this 1 000 times to get a distribution values of the statistic. Figure 9 shows the distribution compared with a ${\chi }^2$ distribution with one degree of freedom. The distributions are close, indicating that the asymptotic distribution of the likelihood-ratio chi-squared test statistic provides a good approximation in this simulated example.

Figure 9:

Comparison of the empirical distribution of the likelihood ratio statistic to test for an additional edge in the VAR equilibrium problem with a ${\chi }^2$ distribution with one degree of freedom. The left-hand plot shows the density of the test statistic compared with the ${\chi }^2$ distribution, while the right-hand plot shows the empirical CDF comparison

5 Application to synthetic data

We have shown by simulation that the asymptotic distributions of estimators and likelihood-ratio test statistics provide good approximations to finite-sample distributions in several example scenarios. What happens when we try to apply the method in situations where the model is not known to be correct? The matrix $\mathbf{D}$ is not likely to be known a priori in a real data setting, and in gene networks the graph is not necessarily restricted to be acyclic. That the general VAR equilibrium model does not require an acyclic graph is a potential advantage of our method over Bayesian network methods, which do require a directed acyclic graph.

5.1 GeneNetWeaver data

To look at the method in a more realistic setting, we need a reasonable size network where the truth is known. For this we used networks from the DREAM4 in silico network challenge competition (http://dreamchallenges.org/project/dream4-in-[silico-network-challenge/]). The 10-gene networks from the competition are a good size for analysis, include cycles, and are known. The original data used in the competition include time-series data simulating a perturbation to the network. The first and last time points do correspond to a stable equilibrium for the network, but these cannot be used for our purposes because the original data do not include enough separate time series from which to get samples.

However, the software used to generate the data, GeneNetWeaver, is available online at http://gnw.sourceforge.net/(Marbach et al., 2009; Schaffter et al., 2011). The software can be configured to generate any number of independent time series from a given network model. The data are generated according to a dynamical model using ordinary differential equations with added stochastic noise. For the time-series data, the network begins at equilibrium, and then the network is artificially perturbed for the first half of the time course. At the midpoint of the time course, the perturbation is removed and the network is allowed to return to its equilibrium state by the end of the time course. The GeneNetWeaver software is designed to provide a rich simulation of the dynamics of biological regulatory networks. The ODEs used include terms simulating transcription and translation rates as well as protein degradation rates, and simulates mRNA levels and protein levels simultaneously. Noise is added to the system to simulate both random fluctuations in the actual levels of mRNA and protein as well as measurement error.

We used GeneNetWeaver to generate time-series data using the network model used in the first 10-gene network from the DREAM4 competition. The GeneNetWeaver software has this network available as a pre-configured model from which data can be generated. The network has 10 nodes and 15 non-self edges, including a cycle involving 3 genes. All settings for data generation were the same as those used in the DREAM4 competition, including the amount of added noise to the system and the model of noise in microarrays. We used GeneNetWeaver to generate 100 random time series of 21 points each from the model and took the first and last time points as samples from the equilibrium distribution.

5.2 Method application

To see how the VAR equilibrium method works on these data, we used the known network structure as our initial model for $\mathbf{A}$ , which includes cycles and is not lower-triangular. We did not assume a known $\mathbf{D}$ , but included it as additional parameters for the model. This resulted in a total of 35 parameters for the model: 10 from $\mathbf{D}$ , 10 diagonal elements of $\mathbf{A}$ and 15 off-diagonal elements of $\mathbf{A}$ , reflecting the 15 non-self edges in the network. The matrix $\mathbf{\Sigma }$ has 55 parameters, so this is well within the parameter counting requirement. This does not guarantee identifiability, but does suggest that it is likely to hold.

To find the MLE for the model, we randomly initialized $\mathbf{A}$ and $\mathbf{D}$ . The diagonal elements of $\mathbf{A}$ were initialized randomly from a Uniform(0, 0.9) distribution to reflect the belief that the autoregressive parameter for the effect of any gene's current expression level on its level at the next time point is positive. The off-diagonal elements were initialized from a Uniform( $-$ 0.9, 0.9) distribution. A requirement on $\mathbf{D}$ is that each element be less than the variance of the corresponding gene in the equilibrium data, so each element was initialized to be the variance of the corresponding gene times, a draw from a Uniform(0.1, 0.9) distribution.

After $\mathbf{A}$ and $\mathbf{D}$ were initialized, the log-likelihood of the data was optimized using the optim function in R with the BFGS method. The optimization scheme may get stuck in a local mode since the likelihood is not convex in $\mathbf{A}$ and $\mathbf{D}$ , so the random initialization and optimization was repeated 1 000 times to provide a better search of the model space. The best $\mathbf{A}$ and $\mathbf{D}$ were kept as the approximate MLE.

This gave an estimate of the parameters for the true model, but in order to test individual edges, we needed to find the MLE for models adjacent to the true model. These models are defined by taking the original network structure and adding or removing a single edge. When removing edges, we needed to ensure that the connectivity of the graph was not broken. Three of the edges in the original network, if removed, would isolate a single node. This leaves 12 of the original 15 edges available to be tested. Adding an edge does not create problems, and thus we were able to test all 75 extra edges. With the results from the true model and the 87 testable models, we could see which edges the likelihood ratio test identified as significant.

Running the above optimization scheme for each of the 88 models resulted in an approximate MLE for each model. Because these optimizations were run independently for each model and the model space is not necessarily fully searched, this resulted in some mis-ordering of the models. That is, the MLE found for one model resulted in a likelihood that was worse than that of a model in which some of the parameters of the first model are constrained to be zero. This is never allowed since the parameters of the smaller model can be used for the larger model, resulting in the same likelihood.

To fix the ordering issue, we performed a second optimization step for each target model. We initialized optim with the parameter values from each of the models. In some cases, this involved throwing out parameter values for extra edges or using zeros from smaller models. The model was then optimized from that starting point. Further, if there are zeroes from a smaller model, we re-randomized just the zeros a small number of times to give better search coverage. We iterated over this step until the models were appropriately ordered, which happened within a few iterations.

5.3 Results

We applied the optimization scheme described earlier to both the initial and final time point data from the generated data from GeneNetWeaver. For each model explored, we obtained the best log-likelihood value for comparison. We then compared the true model with each other model, either with one less or one more edge, and computed the likelihood ratio statistic for testing the edge. From the likelihood ratio statistic we computed a p-value by comparing the value to a ${\chi }^2$ distribution with one degree of freedom.

Ideally, the p-values corresponding to the 12 edges which are in the true network would be low and the p-values for the 75 edges added to the true model would be high. Figure 10 shows the results for both the first and the last time point. Looking at the circles, we find that 7 of the 150 edges added to the true network were identified as true edges at $p=0.05$ . This is a Type I error rate of 5%, showing that the test is producing the correct level of false positive results. Five of the 24 true edges which were tested were identified as true edges, yielding a power of 21%, which is quite promising for this area of research.

Figure 10:

p-values for testing edges in GeneNetWeaver data]p-values for each tested edge in the GeneNetWeaver data. Circles correspond to edges added to the true network (false edges) and triangles correspond to edges removed from the true network (true edges). The left figure is using data from the first time point and the right figure is using data from the last time point. The dashed line corresponds with a $p$ -value of 0.05

6 Discussion

Steady-state gene expression data present a challenge in inferring directional edges. We have presented a model which can be used to test for the existence of edges in such data. We have derived asymptotic properties of estimators and tests for this model in a constrained, but still reasonably large subset of the possible models, and we found that the resulting asymptotic approximations performed well in some finite-sample settings. The derivation of necessary and sufficient conditions for identifiability in the fully general case is a topic for future research.

The identifiability problem is also found in structural equation modeling (SEM) (Bentler and Weeks, 1980). SEM is similar to a VAR model in that it specifies relationships among a set of variables. SEM, however, relates the variables without considering the time element. If $\mathbf{X}$ is the collection of random variables observed, then a linear structural equation model with Gaussian noise can be written as

$$\begin{array}{ccc}\hfill \mathbf{X}& \hfill =\hfill & \mathbf{A}\mathbf{X}+\mathbf{\epsilon },\hfill \\ \hfill \mathbf{\epsilon }& \hfill \sim \hfill & N(0,\mathbf{D}).\hfill \end{array}$$

Seeking identifiability conditions has been the subject of a number of papers (Brito and Pearl, 2012; Drton andWeihs, 2016; Drton et al., 2011), and general conditions for this class of models are not known.

We have shown that even in cases where the theory has not been validated, the VAR equilibrium method can still identify true relationships among genes. This is evident in the analysis of the GeneNetWeaver data, where true edges were identified more consistently than extraneous edges.

If the network is partially known, then the likelihood ratio statistic can be used to test for extra edges and thus build up a more comprehensive view of the network. As an example, we may want to learn about a network perturbed with a certain drug. If the unperturbed network is known, we can test edges in and around that network using the perturbation data to learn about the changes induced by the drug. Since we may not expect the entire network to change, the VAR equilibrium method takes advantage of prior knowledge about the network structure.

To apply the VAR equilibrium-based likelihood ratio test method to steady-state data without a known network to start from will take further development. One way to go about this would be to take a subset of genes and look at all possible network structures for those genes. For each model, the steady-state data would be used to maximize the likelihood. One could then combine the models using a BMA approach. To do this, the integrated likelihood needs to be approximated for each model. As a starting point, the Bayesian Information Criteria (BIC) could be used. The BIC is equal to $-2\text{log}\widehat{L}+k\text{log}n$ , where $\widehat{L}$ is the maximized likelihood for the model, $k$ is equal to $2p$ plus the number of edges in the network, and $n$ is the number of observations. Averaging over all models would result in a posterior probability for each edge in the network.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research,authorship and/or publication of this article.

Funding

This research was supported by NIH grants U54-HL127624, R01-HD054511 and R01-HD070936. Raftery's research was also partly supported by the Center for Advanced Study in the Behavioral Sciences at Stanford University.