State-space Modelling of Replicated Dynamic GeneticNetworks

The genomic reality is highly complex and dynamic. Recent developments of high-throughput technologies have enabled researchers to measure the RNA abundance of thousands of genes simultaneously. The challenge is to unravel from such measurements genomic interactions and key biological features of cellular systems. Two common problems are the high-dimensionality of the system and the spurious correlations induced by unmeasured intermediate substrates. Furthermore most currently available models cannot deal with biological replicates. Our goal is to devise a method for inferring large transcriptional or gene regulatory networks from high-throughput data sources such as gene expression microarrays with potentially hidden states, such as unmeasured transcription


INTRODUCTION
Since the turn of the century a new scientific field has emerged: system biology has started to view biological processes as interrelated events, which ought to be understood in its entirety to make progress within the life sciences [1]. It is a biology-based, but interdisciplinary field that focuses on the systematic study of complex interactions in biological systems. The aim of this holistic approach is to discover new emergent properties that may arise from a systemic view, which are inaccessible to reductionist approaches. The concept of gene networks is central in system biology. A network is an abstract representation of a system, where the substrates or genes are seen as the nodes and the links as some kind of relationship, such as binding or some chemical reaction between them. It is an abstract representation of the stability and interconnectedness of molecular reactions.
The challenge is to give this a precise statistical interpretation in order to allow one to be able to infer the network from quantitative observations on the nodes. Nowadays, expression levels of many genes can be measured simultaneously through many techniques including DNA hybridization arrays [2,3] or RNA-seq methods [4]. A major challenge in system biology is to uncover, from such measurements, gene-protein interactions and key biological features of cellular systems.
The inverse problem of system biology requires a flexible statistical method that in a computationally efficient manner infers the complexity, the dependence structure of the network topology and the functional relationship between the genes. A lot of the statistical system biology literature only consider static networks [5,6,7].
In this paper, we will focus on the well-known linear state space models (SSM) [8,9], which consider dynamic interactions across observed variables and nonobserved states. Several authors have used Kalman filtering of SSM on gene expression data to reverse engineer transcriptional networks. [10] used a two-step approach. In the first step, factor analysis was employed to estimate the state vector and the design matrix; this resulted in the choice of the dimension of the state vector by means of BIC. In the second step, the matrix representing protein-protein translation was estimated using least squares regression. [11,12,13] have applied SSM to T-cell activation data, in which a bootstrap procedure was used to derive a classical confidence interval for parameters representing gene-gene interaction through a re-sampling technique. [14] approached the problem of inferring the model structures of the SSM using variational approximations in the Bayesian context. They used a Variational Bayes method to make the method computationally tractable and identify the dimension of the latent state. [15] also applied SSM to infer the topology of Gene regulatory networks. They introduce an empirical Bayes estimation procedure for a feedback state space model in a hierarchical Bayesian framework that is complementary to the method developed by [14] Recently, [16] used SSMs to rank observed genes in gene expression time series experiments according to their degree of regulation in a biological process.
Their technique is based on Kalman smoothing and maximum likelihood estimation to obtain estimates of the model parameters; however, little attention was paid to the dimension of the hidden state. [17] also presented a novel approach based on the state space model to identify the transcriptional modules and module-based gene networks simultaneously using SSM.
The common problem of all current statistical implementations of SSMs has been that they ignore the existence of biological replicates. [17,18,15] however used technical replicates of gene expression profiles, which are often measured in duplicate or triplicate.
In the presence of biological replicates, time-series are typically averaged out within each time point and the SSM is applied to the average profiles. It is well-known however that replicated genomic time-series typically undergo a gradual phaseshift. These diffusion-like shifts are typically stochastic and not under any genomic control. Averaging out time-series will blur the genomic control and reduce the ability of correct inference. It is our aim to build a dynamic model of replicated dynamic RNA transcripts and unobserved quantities that represent (linear combinations of) commonly unmeasured protein regulators. We infer the model structure as a biological network by estimating model interaction parameters through the EM algorithm [19] combined with the Kalman smoothing algorithm [20,21] in the context of maximum likelihood estimation. We use a bootstrap approach [22] to infer the complex transcriptional response of the network and to reveal interactions between components. Choosing SSM to model network kinetics has a number of advantages. Most importantly, it allows the inclusion of hidden regulators, which can either be unobserved gene expression values or transcription factors (TFs). It can be used to model gene-gene and gene-protein interactions. The parameter estimates obtained through the EM algorithm and the state estimates from the Kalman filter have been shown to be consistent and asymptotically normal under some general conditions [23,24]. In this paper, we demonstrate how the EM algorithm with the Kalman smoothing algorithm are used in the maximum likelihood setup to reverse engineer transcriptional networks from gene expression profiling data. By so doing, we are able to add some useful interpretations to the model. The EM algorithm itself guarantees at least a monotonically increasing likelihood. Model selection or determining a suitable dimension of the hidden state is an additional complication. [12] approached the problem of deciding on a suitable dimension of the hidden state through cross-validation. In their approach, they continuously increased the dimension of the hidden states and monitored the predictive likelihood using the test data.
One major drawback of this approach is that it is very slow and that it cannot be applied in an exploratory analysis of the data. As a result, we focus on faster information-based criteria.
The rest of the paper is organized as follows.
In section 2, we introduce the model and give it a precise mathematical and biological interpretation.
Crucially, we will focus on replicated genomic time-series that will undergo stochastic time-shifts. Section 3 describes the inference method including a model selection procedure for the regulating, but unobserved substrates. Identifiability issues of the model are also discussed and resolved through a minimal number of assumptions. In section 4 we assess via simulation the performance of our method extensively in terms of the F1-score, true positive and false positive rates under various scenarios. Section 5 consists of the application of our model to a well-studied T-cell data set through a bootstrap procedure where we identify the network kinetics, by identifying genetic regulatory networks.
Importantly, our analysis makes explicitly use of the replicated time-series, which has been ignored thus far. We summarize the results, analyze their statistical significance and their biological plausibility. We conclude with a discussion of the method, possible extensions and a summary of related work in section 6.

MODEL
Linear Gaussian state space models, also known as linear dynamical systems [25,26], are a class of dynamic Bayesian networks that relate p temporal observations yt ∈ R p to k temporal hidden state variable θt ∈ R k . We consider a sequence (y1, ..., yT ) of p-dimensional realvalued observation vectors through time, which we shall simply denote by y1:T , representing a gene expression data matrix with p rows and T columns, where p and T are the number of genes and the measuring time points, respectively. The model assumes that the evolution of the hidden variables θt is governed by the state dynamics, which follows a first-order Markov process and is further corrupted by a Gaussian intrinsic biological noise ηt. However, these hidden variable are not directly accessible to the experimenter but rather can be noticed through their effect on the observed data vector, yt, the quantity of mRNA for each of the p genes at time t. The observation yt is a linear transformation of a k-dimensional real-valued θt with observational Gaussian noise ξt. It is assumed that the entire experiment is replicated nR times, resulting in the following model formulation: The terms ηt and ξt are zero-mean independent system noise and measurement noise, respectively with Both Q and R are assumed to be diagonal in many practical applications. The initial state θ0 is independently Gaussian distributed with mean a0 = 0 and covariance Q. This model is more complex and represents an extension of the standard SSM described in [27,28,10] as it includes various forms of feedback and can also be extended to include additional covariates.
The novelty of our approach as compared to other method such as [14] and [12] is the fact that we take biological replication into account. This is a crucial difference as can be seen in Fig. 1.
In this simple 2 gene system with a single latent state, the expression of the same gene, i.e. 1, varies dramatically between replicates, although the underlying kinetics, given by are exactly the same, just as the covariance structure and initial states, given by It is clear that averaging such profiles would lose all biological meaning. Nevertheless, this is what most methods do currently. Our method takes the inter-replicate variability explicitly into account.
It must be noted that the method proposed in [15] is also capable of dealing with replicated time-series gene expression data. Their approach is based on an iterative empirical Bayesian procedure with the introduction of hyperparameters that estimate the posterior distributions of network parameters.
We proposed a complementary method. The novelty in our paper is that we used a Maximum Likelihood inference approach which is a direct inference of the parameters and do not have to worry much about convergence problems . We also used AICc which is a data driven technique in estimating the dimension of the hidden state k.
A mathematical representation of the model is depicted in Fig. 2 indicating the latent and observed dynamics across 3 consecutive time points, where we assumed k = p = 2. The model in Fig. 2 assumes RNA-protein translation at two consecutive time points through the matrix A, and instantaneous protein-RNA transcription through Z. From a biological point of view, the model describes two fundamental stages in gene regulation which are in conformity with the central dogma which states information flows from DNA via RNA to proteins through transcription and translation. The observation-to-state matrix A, is of dimension k × p, and models the influence or the effects of the gene expression values from previous time steps on the hidden states. Matrix B is the p × p matrix indicating the direct genegene interactions. The state dynamic matrix F describes the temporal development of the regulators or the evolution of the transcription factors from previous time step t−1 to the current time step t and is of dimension k × k. It describes the influences of the hidden regulators on each other. The p × k observation dynamics matrix Z relates the transcription factors to the RNAs at a given time point. We collect the model interaction parameters into a single parameter represents our genomic network of interactions including the hidden states. We must point out that [12] focused on the matrix CB + D which is just direct gene-gene interactions. The corresponding CB + D in our case is ZA + B.

Identifiability Issues
A parameter of a dynamic system is said to be identifiable if given some data only one value of this parameter maximizes the observed likelihood. The identifiability property is important because it guarantees that the model parameter can be determined uniquely and with a unique interpretation from the available data. We further assume that the errors {ηt, t = 1, ..., T } and {ξt, t = 1, ..., T } are jointly normal and uncorrelated. Also the number of time points or biological replicates in microarray data are typically much smaller than the number of genes. This fundamental problem of highdimensional statistical modelling of micro array data demands additional care in the estimation of the model parameters in the state space model. This problem is avoided by requiring that the number of observations exceed the total number of parameters to be estimated, i.e., Recall that p is the number of genes, T is the measuring time points, nR is the number of replicates hence we have pT nR as total number of observations. Next, according to our model, we have B, Z, F, and A as parameters to estimate but B is a matrix of dimension for which the system is still identifiable especially for large number of replicates.

The Likelihood Function
With the identifiability constraints from the previous section, we can now write the model parameters as φ ={G, R}. As can be seen from Fig. 2, the observations at time t, ytr are conditioned on the past observations, y (t−1)r and on the regulators θtr. The latent state θtr depends on θ (t−1)r and y (t−1)r . To that effect, we can assume and N d (µ, Σ) is the d-dimensional normal distribution with mean µ and covariance matrix Σ. We write the marginal likelihood function L r y (φ) of the data is given by Maximizing (2.5) across the parameters is extremely challenging, as it involves an integral across the latent state. Therefore, we first consider the complete log-likelihood function of the augmented data (ytr, θtr), given as where the complete log-likelihood of the r th replicate l r yr θr (F, A, Z, B) is given by ignoring constant term.

Joint Parameter Estimation via EM Algorithm
The EM algorithm used in [12] is developed for a single replicate. Extension of this to multiple replicates is non-trivial and non-automatic. There are various ways to extend the SSM to multiple replicates. In our case we assume that each replicate has its own associated latent space governed by the same global parameters (i.e. F and A). The reason is that biological replicates have their own internal clock, governed by universal biological constraints. Our aim is to estimate the model parameters G, which represents the underlying directed genomic network, by maximizing the marginal likelihood function l m y (φ) given in Equation 2.5. Due to the intractability of the integral, we resort to using the EM algorithm [29,30] to learn the parameters of the model.

The EM-algorithm
The main method of inference used in this paper is the Expectation-Maximization (EM) algorithm. Frequentist estimation by and large relies on maximum likelihood (ML) estimators. It consists of maximizing the likelihood across the parameter space. Special cases of the EM algorithm were developed before it was formally introduced by [19]. The EM algorithm has become a popular method of inference in statistical estimation problems involving incomplete data, i.e, data with some missing or latent or hidden observations or problems that can be posed in a similar form, such as mixture models.
The EM algorithm is an iterative tool to compute the maximum likelihood estimate in data characterized by the presence of missing, or hidden or latent observations. This optimization can be difficult especially if the data consist of missing or latent parts. The intuition behind ML is to estimate the parameter(s) for which the observed sample is most likely. It possesses some optimality properties as discussed in [31]. Each iteration of the EM algorithm consists of an expectation step (E-step) followed by a maximization step (M-step). In the E-step, the hidden variables are "estimated" as conditional expectations given the observed data and current estimates of the model parameters. In our SSM, the Kalman filtering algorithm is precisely the Estep. The later is achieved by computing the conditional expectation of the (log) likelihood of the "complete" data. The M-step maximizes the complete likelihood function across the parameter space given the estimate of the missing data from the E-step.
To this effect the algorithm requires the computation of the conditional expectation of the log-likelihood given the complete data. The algorithm is a two-stage procedure, which alternates by calculating the Kalman smoother in the E-step and updating the model parameters in the M-step. The algorithm alternates recursively between an expectation and maximization steps until convergence is obtained.
The procedure to obtain the maximum likelihood estimator of the parameter vector φ via the EMalgorithm is summarized below: 1. Select initial values of (φ0) that is, start with initial guess for the parametersφ0 2. At the k th step, calculate the conditional expectation of the log likelihood (E-step) 3. Determine the next iterative estimated parameters (φ k+1 ) that maximizes conditional expectation of the log likelihood. (M-step) and compute the corresponding log likelihood.
4. Iterate step 2 and 3 until the log likelihood is converged.

The expected log-likelihood function: The E-step
The E-step step of the EM algorithm involves the calculation of the first two moments of the hidden states θt. Let Q denote the expected loglikelihood. Then from Equation 2.6, Q becomes where φ * = (Z * , B * , F * , A * ) is the estimate obtained from the previous M-step The calculation of Q(φ|φ * ) in Equation 2.8 involves finding E(θ) and E(θ ′ θ) for each replicate r. These forms are readily found: for each replicate we run the Kalman smoothing algorithm to obtain the expected hidden states and their variance-covariance components and these are joined together to get Q(φ|φ * ).

The Kalman-Filtering Algorithm
The Kalman filter has been considered as one of the optimal solutions to many data prediction, filtering and smoothing problems. In this context, it is used to estimate the hidden or latent states in the E-step of the EM algorithm. We describe here the basic concepts that one needs to know to design and implement a Kalman filter algorithm. Given our original model from equation (2.1) The predictive step equations are given by: wherePt represents the corresponding predicted or prior state estimate error covariance.
Then the observation prediction equation step becomes:ỹ where Σt represent the observation prediction covariance.
with vt = yt −ỹt (2.13) denoting the measurement innovation or the residual and reflects the discrepancy between the predicted measurementỹt and the actual observation yt The filtered equations are also given bŷ is the Kalman gain matrix and is chosen to be the gain or blending factor that minimizes the posterior error covariance in Equation (2.15).
Finally the smoothing step is given bŷ where H is the Kalman smoothing matrix. More information on K and H can be obtained from [20,21] Equation (2.16) represents the expected hidden states needed at the E-step.

The update equations: The Mstep
A new parameter set φ * is computed by estimating the parameters that maximize the two quadratic forms in (2.8). We solve ∂ ∂φ Q1 = 0 and ∂ ∂φ Q2 = 0 to obtain estimates for Z, B and F, A, respectively. This can be solved in closed form. For a full derivation, see the appendix. The entire EM algorithm can be regarded as alternating between Kalman smoothing and least squares minimization given by the update equations.

Choice of Hidden State Dimension: AIC c
Model selection or the determination of the optimum dimension of the hidden state k is important in the application of SSM to network reconstruction. Popular model selection criteria include Akaike's Information Criterion (AIC) [32] and the Bayesian Information Criterion (BIC) [33]. We apply a corrected Akaike's Information Criterion (AICc) method in our scenario. The AICc has good model estimation properties, especially for small sample time-series data [34,35]. Furthermore, the CV approach used in [12] tends to be slow and unstable for small number of replicates.
The AICc is aimed at finding the best approximating model to the unknown data generating process via minimizing the estimated expected Kulback-Leibler divergence. Given the log-likelihood function l, the AIC for a model with k-dimensional state vector is given by: with P the number of estimated parameters, and l(yt|φ k ) the log-likelihood of the observed data. [36] recommends the use of the corrected AIC, which corrrects for small sample size bias.
In the framework of normal linear regression models, the penalty term of AICc provides an exact expression for the bias adjustment. The AICc is given by where N = pT nR represents total number of observations and P = p 2 + 2kp + k 2 is the total number of estimated parameters. We select the hidden state dimension that has the minimum AICc, i.e we find k such that (2.20) We successively increase the number of hidden states and monitor the behavior of AICc as a function of k.

Network Reconstruction by Bootstrapping
We use a bootstrap approach to find confidence intervals for the parameters defined in our model. By so doing we compute the bootstrap distribution of the estimator of φ. Letφ denote the MLE of the parameters defined in our model that come from using the EM algorithm described in previous section. In the following we will use the notation yr ∈ R P ×T with r ∈ {1, ..., nR} to represent each of the biological time series. The bootstrap procedure adopted is outlined below: Whereas we bootstrap the residuals in our method, [12] bootstrapped the original observations. As our approach is nonparametric, the two methods result in the same bootstrapping procedure. Given each new data we estimate, among other things, the bootstrap set of parameters through the EM algorithm. Stated differently, for each bootstrap sample the parameters that maximize the likelihood of the bootstrap data are found. We then obtain the sampling distributions of the estimators of the elements of φ. The results of the bootstrapping are the distribution of the parameters and we proceed to make statistical inferences about those underlying parameters by computing confidence interval for each of them [37,20].

SIMULATION STUDIES
In order to illustrate the performance of our method for analyzing gene expression data, we simulate artificial data and applied our proposed method to these data according to the model described in Equation 2.1 with T = 10 time points, p = 3 genes and k = 2 latent states. The true newtork is depicted on the left in Fig. 3.
In the initialization step of the EM-algorithm, Z and F are assumed to be identity matrices, whereas A to set to zero.
For B we perform a simple linear regression where we regress all observed genes on the previous time point. The diagonal variance matrix R is given the usual variance estimate coming from the regression. We apply the bootstrap procedure to the data and identify the significant and non-significant parameters through the use of bootstrap confidence intervals on the element Gij of the network. For this decision problem where we formulate two hypotheses, namely,

H0
: where rejecting H0 indicates the presence of a connection among the gene i to j. With k equals 2, we obtained network shown on the right in Fig.  3. Table 1 depicts the performance of our method as the size of the network increases. Comparison is also made to the method of [12], in which no replicates are considered. In order to deal with the replicates, an average profile across 50 biological replicates is calculated. When p increases, our method is able to detect the network more accurately. The reason is that the number of latent states k = 2 for larger network in terms of p becauses relatively speaking smaller, which makes network detection easier. Ignoring the inter-replicate variability, however, means that, unlike it is the case for our method, there is no gradual improvement in network detection for the other SSM method.

APPLICATION
For this study, to demonstrate the application of our network inference method, we used publicly available data. Two separate experiments investigated the expression response of human T-cells to PMA and ionomicin treatment. The entire data set is a combination of the datea from these two experiments. The first data set (tcell.34) contains the temporal expression levels of 58 genes for 10 unequally spaced time points.
At each time point there are 34 separate measurements. The second data set (tcell.10) comes from a related experiment considering the same genes and identical time points, and contains 10 further measurements per time point. Therefore, at each time point there are 44 separate measurements or replicates. Corresponding to each gene expression ytr, we also assumed the existence of generative replicates for the hidden variables θtr. With p = 58 genes, R = 44 as replicates and T = 10, the constraint represented by Equation (2.3) is satisfied, indicating that we have enough data to estimate our parameters. The dimension of the hidden variables was determined using AICc as explained in section 2.4. Table 2 shows the behavior of AICc with corresponding k's. It turns out that k = 4 is the optimum number of the hidden states. This is fewer than [12] and [14] who obtained 9, 14 respectively under different criteria.
In essence, we treated the data as a time series measurement data yt r , t = 1, 2, ..., 10 and r = 1, 2, ..., 44. For each replicate, yt and θt consist of 58 genes and 4 transcriptions factors respectively, each, measured at 10 different time points, i.e for each replicate r, yt and θt are of dimension (58 × 10), (4 × 10) respectively. Some of these genes include RB1, CCNG1, TRAF5, CLU.... The parameters Q were fixed.  Based on 95% confidence intervals to detect significant interactions, we plot the connectivity matrix of the directed genomic networkĜ. The output is a directed graph showing connections from one gene expression variable at a given time point t to another gene expression variable whose expression it influences at the next time point, t + 1. The arrows indicates the direction of the regulation. The entire directed grapĥ G gives 350 genomic interactions. Fig. 4 represents a portion of the interaction network φ where we indicate genes that have at least 3 outwards connections. These genes include the FYN-binding protein gene FYB, the JUND protooncogene, the CD69 antigen p60, early T-cell activation antigen to mention but a few. Fig. 5 is the sub-network produced at 95% confidence level and it represents the interaction between, two Jun proteins family namely JUNB and JUND and various genes involved in programmed cell death. The results of our method in Fig. 5 support the hypothesis of the anti-proliferation and antiapoptotic role of JUND.
According to our method, the following genes were mostly seen as regulatory genes. These genes include the JUND proto-oncogene, the CLU gene, the cell division cycle 2 CDC2, the FYN-binding protein gene FYB, TRAF5, the CD69 and the GATA-binding protein 3. The latent variables were also seen to regulate the expression level of most genes as can be seen in Fig. 4. Our approach has revealed interesting features in the family of Jun genes. The network in Fig. 5 provides support for several hypotheses that were also confirmed in [12] and [14]. However, we also found new connections. Our results support the interaction between the protooncogene JUNB, the apoptosis-related cysteine protease genes CASP4 and CASP8. The implication is that JUNB is clearly modelled as a pro-apoptotic gene by activating CASP4 and CASP8. This interaction was also was recovered by [14].
We, however, found no evidence for interaction between JUNB and MAP3K8. Also Fig. 5 reveals that the proto-oncogene JUND activates the GATA-binding protein 3, but represses the expression level of the cell division cycle 2 (CDC2). This further supports the anti-proliferative JUND. Furthermore, in our model, the survival of motor neuron 1 gene SMN1 and the cell division cycle CDC2 influence the expression level of JUNB and MAPK8 respectively. JUNB activates the expression level of CDC2. A critical comparison of our Fig. 5 to that of similar sub-networks found in the work of [15] and [14] shows that in all the 3 subnetworks, JUND regulates the expression level of CDC2. JUNB activates CASP8 in the subnetwork found by [14] and indirectly regulates CASP8 through CASP4 in the sub-network found by [15]. However we found interaction between JUNB and both CASP8 and CASP4.
The gene FYN-binding protein FYB has been found to occupy one of the most crucial positions in the network recovered by [12] also it has a high degree of connectivity in our work. Fig. 6 reveals some crucial genes that are found to be directly connected to FYB. Most importantly, in our model FYB influences the expression level of genes such as the early T-cell activation marker CD69, the JUNB proto-oncogen. FYB is also seen to be connected to genes such as APC, API2, and CIR. Clearly, these results support the fact that FYB mRNA levels are predictive of the expression level of a number of genes. The hidden state dimensionality was found to be 4, a result similar to the work of [15] in which they developed an iterative empirical Bayesian procedure with a Kalman filter to estimate the posterior distributions of network parameters. [12] found the dimension of the hidden state to be 9 through cross validation, while [14] obtained the value of 14 through a variational Bayesian approach. At a 95% confidence level, we found no significant interactions among the hidden variables or transcription factors. However their role in the transcription process can not be ignored as the inferred matrix Z representing instantaneous protein-RNA transcription was not sparse signifying that the transcription factors regulate the expression level of most mRNAs.

CONCLUSION
In this paper, we have developed a state space model with biological replication and applied it to the T-cell data. We used the EM algorithm combined with the bootstrap to infer the structure of the underlying genomic network. The proposed method offers significant advantages over other methods that have recently appeared in the literature. For example, [14] used a variational Bayesian methodology which is an approximation of the posterior distribution of the parameters, whereas we obtained much faster results through direct inference of the parameters. [12] used cross validation as a model selection technique which is quite slow as compared to AIC. [16] used an ad hoc method for selecting the hidden state dimensionality k, while our method uses a data-driven approach. Also our model allows for dynamic correlation over time, as each observation and hidden state depend explicitly on some function of previous observations as opposed to the model described by [27,28,10]. Their model does not allow for RNA-protein translation and RNA-RNA interactions through the matrix A and B respectively in our model.
One fundamental assumption in our proposed model is the first-order linear dynamics in the state and observation equations of the SSM. This assumption can only be an approximation to the true nature of a complex biological system since more realistic models of gene regulatory interactions surely include complex interactions or nonlinear relationships. Our linear dynamics assumption is a stepping stone upon which a future model with non-linear dynamics will be explored. With application to the t-cell data, we have discovered new interactions that have not yet been reported in the current literature; as as part of our ongoing work we are investigating these interactions further.
Furthermore the AIC approach is prone to overfitting, especially in high-dimensional data. A natural way to avoid this over-fitting is through regularization.
Given that Most of existing time-series gene expression data have much fewer replicates, e.g., 3, 2, 1 or no replicates our bound condition (2.4) may be broken down making our system non-identifiable. We can overcome this situation by imposing a penalty on the parameters. In future work, we plan to employ a penalized maximum likelihood strategy in the context of the EM algorithm in the state space model.

APPENDIX
We outline here the derivations of the update equations. We write Q(Z, B) as