Modelling Nonstationary Gene Regulatory Processes

Abstract

An important objective in systems biology is to infer gene regulatory networks from postgenomic data, and dynamic Bayesian networks have been widely applied as a popular tool to this end. The standard approach for nondiscretised data is restricted to a linear model and a homogeneous Markov chain. Recently, various generalisations based on changepoint processes and free allocation mixture models have been proposed. The former aim to relax the homogeneity assumption, whereas the latter are more flexible and, in principle, more adequate for modelling nonlinear processes. In our paper, we compare both paradigms and discuss theoretical shortcomings of the latter approach. We show that a model based on the changepoint process yields systematically better results than the free allocation model when inferring nonstationary gene regulatory processes from simulated gene expression time series. We further cross-compare the performance of both models on three biological systems: macrophages challenged with viral infection, circadian regulation in Arabidopsis thaliana, and morphogenesis in Drosophila melanogaster.

1. Introduction

An important objective in systems biology is to infer regulatory networks from postgenomic data. Bayesian networks have been widely applied as a popular tool to this end—see, for example, [1]—and novel fast Markov chain Monte Carlo (MCMC) algorithms can be applied to systematically search the space of network structures for those that are most consistent with the data [2, 3]. One reason for the popularity of Bayesian networks has been the tractability of the marginal likelihood of the network structure. This term describes how well the model structure explains the data. Its computation is usually challenging, as it requires the solution of an integral over the entire parameter space. To obtain a closed-form expression, two probabilistic models with their respective conjugate prior distributions have been employed: the multinomial distribution with the Dirichlet prior (BDe) and the Gaussian distribution with the normal-Wishart prior (BGe). These approaches are restricted in that they either require a data discretisation (BDe: Bayesian Dirichlet equivalence score) or can only capture linear relationships (BGe: Bayesian Gaussian equivalence score). Recently, a generalisation of the BGe model based on a combination of a mixture model and the allocation sampler has been proposed [4], the Bayesian Gaussian Mixture (BGM) model. In the BGM model, data points are assigned to different compartments (subsets of the data) with the allocation sampler [5]. Model parameters (and their distributions) are allowed to differ between compartments, while information is shared among the compartments via a common network structure. Given the network structure each compartment is modelled separately and independently with the Gaussian BGe model.

The present work proposes a modification of the BGM model, which is more suitable for dynamic gene regulatory networks and gene expression time series. The proposed model, which we refer to as the BGM_𝐷 model, replaces the free allocation model by a multiple changepoint process, as, for example, employed in [6], to divide the time series into segments and thereby take the time structure into account.

In a preliminary study [7], we focused on synthetic data from small network domains. We found that the new model avoids spurious self-feedback loops if feedback loops, such as 𝑋(𝑡)→𝑋(𝑡+1) (e.g., in molecular biology: transcription factors regulating their own transcription), are not excluded from the analysis. In this paper, we present the mathematical details of the new BGM_𝐷 model in greater depth, and we demonstrate that BGM_𝐷 yields better inference results for real biological systems. To this end, we analyse two gene expression time series from macrophages and Arabidopsis thaliana and cross-compare the inference results of the BGM model and the new changepoint process BGM_𝐷 model.

Furthermore, we combine both models (the free allocation model and the changepoint process) with discrete Bayesian network methodology and we compare the performance of both models on inferring the morphogenetical stages of muscle development in Drosophila melanogaster from binary gene expression time series.

We note that our modelling paradigm is complementary to other recently proposed approaches. A nonlinear nondiscretised model based on heteroscedastic regression has been proposed in [8]. However, this approach no longer allows the marginal likelihood to be obtained in closed-form and requires a restrictive approximation (the Laplace approximation, that is, an approximation based on a 2nd order Taylor series approximation of the log likelihood) to be adopted. Another nonlinear model based on node-specific Gaussian mixture models has been proposed in [9, 10]. Again, the marginal likelihood is intractable. The authors resort to the Bayesian information criterion (BIC) of [11] for model selection, which is only a good approximation to the marginal likelihood in the limit of very large data sets. A more flexible approach based on changing graphs between changepoints has been proposed in [12, 13]. Conceptually, the assumption of changing networks is reasonable for some biological scenarios, like morphogenesis, where gene-regulatory processes have been measured over a long period of time. However, for cellular processes on a shorter time scale, it is questionable whether it is the network structure rather than just the strength of the regulatory interactions that changes with time. For example, it is not the principle ability of a transcription factor to potentially bind to the promoter of a gene and thereby initiate transcription (i.e., not the network structure), but the extent to which this happens (i.e., the regulatory interaction strength). We therefore argue that, especially for short time series, it is more appropriate to leave the network structure invariant among segments and to allow the interaction parameters to vary with time by modelling the time segments (between changepoints) separately and independently with the Gaussian BGe model. The systematic sharing of information among segments via a common network structure also reduces potential problems with model overflexibility and overfitting, which a more flexible approach that allows for different network structures is susceptible to when the time series are short.

2. Methods

2.1. The Dynamical BGe Network

Dynamical Bayesian networks (DBNs) are flexible models for representing probabilistic relationships among interacting variables 𝑋₁,…,𝑋_𝑁. The graph 𝒢 of a DBN describes the relationships between the variables (nodes), which have been measured at equidistant time points 𝑡=1,…,𝑚, in the form of conditional probability distributions. An edge pointing from 𝑋_𝑖 to 𝑋_𝑗 means that the realisation of 𝑋_𝑗 at time point 𝑡, symbolically: 𝑋_𝑗(𝑡), is influenced by the realisation of 𝑋_𝑖 at time point 𝑡−1, symbolically: 𝑋_𝑖(𝑡−1). Effectively there is a bipartite graph structure between two time steps 𝑡 and 𝑡+1 so that the acyclicity constraint—known from static Bayesian networks—is guaranteed to be satisfied. Therefore, in principle, each node can be its own parent node in DBNs. Such self-feedback loops 𝑋_𝑛(𝑡−1)→𝑋_𝑛(𝑡) model autocorrelations and it depends on the application whether such self-feedback loops should be valid edges or ruled out altogether to focus the inference on the interactions between genes. 𝜋_𝑛=𝜋_𝑛(𝒢) denotes the parent (node) set of node 𝑋_𝑛 in 𝒢; that is, the set of all nodes from which an edge points to node 𝑋_𝑛 in 𝒢. Given a data set 𝐃, where 𝐃_𝑛,𝑡 and 𝐃_{(𝜋𝑛,𝑡)} are the 𝑡th realisations 𝑋_𝑛(𝑡) and 𝜋_𝑛(𝑡) of 𝑋_𝑛 and 𝜋_𝑛, respectively, DBNs are based on the following homogeneous Markov chain expansion:=𝑃(𝐃∣𝒢,𝜽)_{𝑁𝑚𝑛=1𝑡=2}𝑃(𝑋_𝑛(𝑡)=𝐃_𝑛,𝑡∣𝜋_𝑛(𝑡−1)=𝐃_{(𝜋𝑛,𝑡−1)},𝜽_𝑛), (1) where 𝜽 is the total parameter vector, composed of subvectors 𝜽_𝑛. 𝜽_𝑛 specifies the 𝑛th local conditional distributions 𝑃(𝑋_𝑛(𝑡)∣𝜋_𝑛(𝑡−1),𝜽_𝑛) in the factorisation. The BGe model specifies the distributional form 𝑃(𝐃∣𝒢,𝜽) as multivariate Gaussian distribution, and it assumes a normal-Wishart distribution as prior distribution 𝑃(𝜽|𝒢) [14]. The local probability distributions 𝑃(𝑋_𝑛(𝑡)|𝜋_𝑛(𝑡−1),𝜽_𝑛) are then given by conditional linear Gaussian distributions. For discrete data 𝐃 the BDe model specifies the distributional form of the likelihood 𝑃(𝐷|𝒢,𝜽) as a set of multinomial distributions, and assumes a Dirichlet distribution as the prior 𝑃(𝜽∣𝒢) for the unknown parameters [15]. In principle, the BDe model yields a higher modelling flexibility than the BGe model, but BDe requires a data discretisation that usually incurs a substantial information loss. Under fairly weak conditions imposed on the parameter vector (prior independence and modularity) and the prior distribution 𝑃(𝜽) (conjugacy), so that 𝑃(𝜽|𝒢)=∏^𝑁_𝑛=1𝑃(𝜽_𝑛|𝜋_𝑛), where 𝜋_𝑛=𝜋_𝑛(𝒢), the parameters can be integrated out analytically. For both scoring metrics BGe and BDe [16], the marginal likelihood then satisfies the same expansion rule as the Bayesian network with fixed parameters:𝑃(𝐃∣𝒢)=𝑃(𝐃∣𝒢,𝜽)𝑃(𝜽∣𝒢)𝑑𝜽=^𝑁_𝑛=1Ψ(𝐃^𝜋𝑛_𝑛), (2) whereΨ(𝐃^𝜋𝑛_𝑛)=^𝑚_𝑡=2𝑃(𝑋_𝑛(𝑡)=𝐃_𝑛,𝑡∣𝜋_𝑛(𝑡−1)=𝐃_{(𝜋𝑛,𝑡−1)},𝜽_𝑛)×𝑃(𝜃_𝑛∣𝜋_𝑛)𝑑𝜽_𝑛, (3) and 𝐃^𝜋𝑛_𝑛∶={(𝐃_𝑛,𝑡,𝐃_{𝜋𝑛,𝑡−1})∶2≤𝑡≤𝑚} denotes the subset of the data pertaining to node 𝑋_𝑛 and its parent set 𝜋_𝑛. For the Gaussian BGe model the (local) factors Ψ(𝐃^𝜋𝑛_𝑛) in (3) can be computed in closed-form, according to (15) and (24) in Geiger and Heckerman [14]. For the discrete multinomial BDe model, the (local) factors can be computed in closed form according to (19) in [16].

2.2. Gaussian Mixture Bayesian Network Model

In the Gaussian BGe model, the local distributions 𝑃(𝑋_𝑛(𝑡)|𝜋_𝑛(𝑡−1),𝜽_𝑛) are conditional linear Gaussian distributions. We generalise the BGe model by the introduction of a latent allocation vector 𝐕, which assigns the data points to 𝐾 different mixture components, where 𝐾 is inferred from the data by applying changepoint birth and death moves, along the lines of the changepoint model (e.g., see [6, 17]). As in the BGM model [4], conditional on the latent vector 𝐕, a separate BGe score can be computed for each of the 𝐾 mixture components.

The allocation vector 𝐕 of size 𝑚−1 describes the allocation of the time points 𝑡=2,…,𝑚 to the 𝐾 components. 𝐃^(𝑉,𝑘) denotes all time points that are allocated to component 𝑘. The posterior probability of 𝒢, 𝐕, and 𝐾 is proportional to the joint distribution𝑃(𝒢,𝐕,𝐾∣𝐃)=𝑃(𝒢,𝐕,𝐾,𝐃)𝑃(𝐃)∝𝑃(𝒢,𝐕,𝐾,𝐃), (4) and the joint distribution can be factorised as follows:𝑃(𝒢,𝐕,𝐾,𝐃)=𝑃(𝐾)𝑃(𝐕∣𝐾)𝑃(𝒢)𝑃(𝐃∣𝒢,𝐕,𝐾), (5) where𝑃(𝐃∣𝒢,𝐕,𝐾)=^𝐾_𝑘=1𝑃(𝐃^(𝐕,𝑘)∣𝒢)=^𝐾𝑁_𝑘=1𝑛=1Ψ(𝐃^{(𝐕,𝑘),𝜋𝑛}_𝑛), (6) and 𝐃^{(𝑉,𝑘),𝜋𝑛}_𝑛∶={(𝐃_𝑛,𝑡,𝐃_{𝜋𝑛,𝑡−1})∣𝑡∈{2,…,𝑚}∶𝐕(𝑡)=𝑘} denotes the realisations of node 𝑋_𝑛 and its parent set 𝜋_𝑛 for those time points that have been allocated to the 𝑘th component. It can be seen from these equations that 𝐕 acts as a filter which divides the data 𝐃 into 𝐾 different compartments 𝐃^(𝐕,𝑘) (𝑘=1,…,𝐾), for which separate independent BGe scores can be computed in closed-form using (2) and (3). The BGM counterpart of (3) is given by:Ψ(𝐃^{(𝐕,𝑘),𝜋𝑛}_𝑛)=_{𝑡∶𝐕(𝑡)=𝑘}𝑃(𝑋_𝑛(𝑡)=𝐃_𝑛,𝑡∣𝜋_𝑛(𝑡−1)=𝐃_{(𝜋𝑛,𝑡−1)},𝜽_𝑛)×𝑃(𝜽_𝑛∣𝜋_𝑛)𝑑𝜽_𝑛. (7) For instance, if we have 𝑚=11 time points and one changepoint between 𝑡₆ and 𝑡₇ so that 𝐕 assigns the time points 𝑡₂,…,𝑡₆ to the first and the remaining time points 𝑡₇,…,𝑡₁₁ to the second segment, then separate local BGe scores are computed for the data subsets 𝐃^{(𝐕,1),𝜋𝑛}_𝑛∶={(𝐃_𝑛,𝑡,𝐃_{𝜋𝑛,𝑡−1})|2≤𝑡≤6} and 𝐃^{(𝐕,2),𝜋𝑛}_𝑛∶={(𝐃_𝑛,𝑡,𝐃_{𝜋𝑛,𝑡−1})∣7≤𝑡≤11}, according to (15) and (24) in Geiger and Heckerman [14].

When a data compartment 𝐃^(𝐕,𝑘) is empty, then we set the factors Ψ(𝐃^{(𝐕,𝑘),𝜋𝑛}_𝑛) equal to 1 (𝑛=1,…,𝑁). For 𝑃(𝒢), we take a uniform distribution over all graph structures subject to a fan-in restriction of |𝜋_𝑛|≤ℱ. For 𝑃(𝐾), we take a truncated Poisson distribution with 𝜆=1 restricted to 1≤𝐾≤𝐾^⋆ as prior. We note that the MCMC inference scheme, which we will discuss in Section 2.3, does not sample 𝐕 directly, but is based on local modifications of 𝐕 based on changepoint birth, death, and reallocation moves. That is, different from the free allocation in the BGM model [4], we change the assignment of data points to components via a changepoint process (e.g., see [6, 17]). This reduces the complexity of the allocation space substantially and incorporates our prior knowledge that adjacent time points are likely to be assigned to the same component. We will refer to the new changepoint BGM model as the BGM_𝐷 model.

We identify 𝐾 with 𝐾−1 changepoints: 𝑏₁,…,𝑏_𝐾−1 on the continuous interval [2,𝑚], and for notational convenience we introduce the pseudo-changepoints 𝑏₀=2 and 𝑏_𝐾=𝑚. The observation at time point 𝑡 is assigned to the 𝑘th component, symbolically 𝐕(𝑡)=𝑘, if 𝑏_𝑘−1≤𝑡<𝑏_𝑘. Following Green [17], we assume that the changepoints are distributed as the even-numbered order statistics of 𝐿∶=2(𝐾−1)+1 points 𝑢₁,…,𝑢_𝐿 uniformly and independently distributed on the interval [2,𝑚]. The motivation for this prior, instead of taking 𝐾−1 uniformly distributed points, is to encourage a priori an equal spacing between the changepoints, that is, to discourage mixture components that contain only a short compartment of the time series. We note that the even-numbered order statistics prior on the changepoint locations induces a prior distribution on the allocation vector 𝐕. Deriving a closed-form expression is involved but not required, as our MCMC inference scheme does not sample 𝐕 directly. The inference is based on local modifications of 𝐕. In particular, we employ three different move types: (i) an additional component (𝐾→𝐾+1) can be generated by setting a new changepoint on the interval [2,𝑚] (changepoint birth move), (ii) the number of components can be decremented (𝐾→𝐾−1) by removing one of the existing changepoints (changepoint death move), and finally, (iii) the allocation vector 𝐕 can be changed without affecting the number of components 𝐾 by changing the position of one of the existing changepoints (changepoint reallocation move). For the acceptance probabilities of these local moves, only 𝑃(𝐕∣𝐾) ratios, which are straightforward to compute, are required.

2.3. Multinomial Mixture Bayesian Network Model

The discrete multinomial BDe model can be generalised in analogy to the continuous Gaussian Bge model. As before, the allocation vector 𝐕 divides the data into 𝐾 compartments, and each compartment of the data is modelled separately and independently with the multinomial BDe scoring metric. We note that closed-form solutions for the Ψ(⋅)-terms in (7) can be obtained with (19) in [16]. We will refer to the free individual allocation model BGM with the BGe score being replaced by the discrete multinomial BDe score as the Bayesian Discrete Mixture (BDM) model, and accordingly we will refer to the discrete BDe counterpart of BGM_𝐷 as the BDM_𝐷 model. We note that the BDM_𝐷 model is similar to the 𝑛𝑠𝐵𝐷𝑒 model of Robinson and Hartemink [12], except that BDM_𝐷 leaves the inferred network structures invariant in time to allow for more information sharing among segments. We will employ the discrete counterparts BDM and BDM_𝐷 for an independent comparison of the suitability of the free allocation model and the changepoint process for inferring dynamic gene-regulatory processes (from discrete data).

2.4. MCMC Inference

We now describe an MCMC inference algorithm that can be used to obtain a sample {𝒢^𝑖,𝐕^𝑖,𝐾^𝑖}_{𝑖=1,…,𝐼} from the posterior distribution 𝑃(𝒢,𝐕,𝐾|𝐃). Our algorithm combines the structure MCMC algorithm for Bayesian networks [18] with the changepoint model (e.g., see [6, 17]), and draws on the fact that conditional on the allocation vector 𝐕, separate BGe scores 𝑃(𝐃^(𝐕,𝑘)|𝒢) can be computed for the 𝐾 data compartments. Note that this approach is equivalent to the idea underlying the allocation sampler [5]. The resulting algorithm is effectively an RJMCMC sampling scheme in the discrete space of network structures and latent allocation vectors, where the Jacobian in the acceptance criterion is always 1 and can be omitted. With probability 𝑃=.5, we perform a traditional structure MCMC move on the current graph 𝒢^𝑖 and leave the latent vector 𝐕 and the number of mixture components 𝐾 unchanged, symbolically: 𝐕^𝑖+1=𝐕^𝑖 and 𝐾^𝑖+1=𝐾^𝑖. A new candidate graph 𝒢^𝑖+1 is randomly drawn out of the set of graphs 𝒩(𝒢^𝑖) that can be reached from the current graph 𝒢^𝑖 by deletion or addition of one single edge. The proposed graph 𝒢^𝑖+1 is accepted with probability𝐴(𝒢^𝑖+1∣𝒢^𝑖)=min{1,𝑃(𝐃∣𝒢^𝑖+1,𝐕^𝑖,𝐾^𝑖)𝑃(𝐃∣𝒢^𝑖,𝐕^𝑖,𝐾^𝑖)⋅𝑃(𝒢^𝑖+1)𝑃(𝒢^𝑖)⋅|𝒩(𝒢^𝑖)||𝒩(𝒢^𝑖+1)|}, (8) where |⋅| is the cardinality, and the marginal likelihood terms have been specified in (6). The graph is left unchanged, symbolically 𝒢^𝑖+1∶=𝒢^𝑖, if the move is not accepted. We note that the network reconstruction will be based on the marginal posterior probabilities of the individual edges, which can be estimated for each edge from the MCMC sample 𝒢¹,…,𝒢^𝐼 by the fraction of graphs in the sample that contain the edge of interest𝐸_𝑖,𝑗=1𝐼^𝐼_𝑘=1𝐼_𝑖,𝑗(𝒢^𝑘), (9) where 𝐼_𝑖,𝑗(⋅) is the indicator function with 𝐼_𝑖,𝑗(𝒢^𝑘)=1 if there is an edge from 𝑋_𝑖 to 𝑋_𝑗 in 𝒢^𝑘.

With the complementary probability 1−𝑃 we leave the graph unchanged: 𝒢^𝑖+1=𝒢^𝑖, and perform a move on (𝐕^𝑖,𝐾^𝑖). We change the current number of components 𝐾^𝑖 via a changepoint birth or death move, or the allocation vector 𝐕^𝑖 by a changepoint reallocation move along the lines of the Reversible Jump Markov chain Monte Carlo algorithm (RJMCMC) algorithm [17].

The changepoint birth (death) move increases (decreases) 𝐾^𝑖 by 1 and may also have an effect on 𝐕^𝑖. The changepoint reallocation move leaves 𝐾^𝑖 unchanged and may have an effect on 𝐕^𝑖. If with probability (1−𝑃) a changepoint move on (𝐾^𝑖,𝐕^𝑖) is performed, we randomly draw the move type. Under fairly mild regularity conditions (ergodicity), the MCMC sampling scheme converges to the desired posterior distribution if the equation of detailed balance is fulfilled [17]. The condition of detailed balance implies that for each move a complementary move is defined, and that the acceptance probability depends on the proposal probability of the complementary move. The moves presented below are designed such that there is a unique complementary death move for each birth move and vice versa. Moreover, each reallocation move can be reversed by a single (complementary) reallocation move. The acceptance probabilities for these three changepoint moves (𝐾^𝑖,𝐕^𝑖)→(𝐾^𝑖+1,𝐕^𝑖+1) are of the following form [17]:𝐴=min{1,𝑃(𝐃∣𝒢^𝑖,𝐕^𝑖+1,𝐾^𝑖+1)𝑃(𝐃∣𝒢^𝑖,𝐕^𝑖,𝐾^𝑖)×𝑅×𝐵}, (10) where 𝑅=𝑃(𝐕^𝑖+1|K^𝑖+1)𝑃(𝐾^𝑖+1)/𝑃(𝐕^𝑖|𝐾^𝑖)𝑃(𝐾^𝑖) is the prior probability ratio, and 𝐵 is the inverse proposal probability ratio. The exact form of these factors depends on the move type.

(i) For a changepoint reallocation (𝑟), we randomly select one of the existing changepoints 𝑏_𝑗∈{𝑏₁,…,𝑏_𝐾−1}, and the replacement value 𝑏^†_𝑗 is drawn from a uniform distribution on [𝑏_𝑗−1,𝑏_𝑗+1] where 𝑏₀=2 and 𝑏_𝐾=𝑚. Hence, the proposal probability ratio is one, the prior probabilities 𝑃(𝐾^𝑖+1)=𝑃(𝐾^𝑖) cancel out, and the remaining prior probability ratio 𝑃(𝐕^𝑖+1|𝐾^𝑖+1)/𝑃(𝐕^𝑖|𝐾^𝑖) can be obtained from page 720 in Green's RJMCMC paper [17]𝑅_𝑟=(𝑏_𝑗+1−𝑏^†_𝑗)(𝑏^†_𝑗−𝑏_𝑗−1)(𝑏_𝑗+1−𝑏_𝑗)(𝑏_𝑗−𝑏_𝑗−1),𝐵_𝑟=1. (11) If there is no changepoint (𝐾^𝑖=1) the move is rejected and the Markov chain is left unchanged.

(ii) If a changepoint birth move (𝑏) on 𝐾^𝑖 is proposed, the location of the new changepoint 𝑏^† is randomly drawn from a uniform distribution on the interval [2,𝑚]; the proposal probability for this move is 𝑏_𝐾^𝑖/(𝑚−2), where 𝑏_𝐾^𝑖 is the (𝐾^𝑖-dependent) probability of selecting a birth move. The reverse death move, which is selected with probability 𝑑_(𝐾^𝑖₊₁₎, consists in discarding randomly one of the (𝐾^𝑖−1)+1=𝐾^𝑖 changepoints. The inverse proposal probability ratio is thus given by 𝐵=𝑑_(𝐾^𝑖₊₁₎(𝑚−2)/(𝑏_𝐾^𝑖𝐾^𝑖). The prior probability ratio is given by the expression at the bottom of page 720 in Green's RJMCMC paper [17] slightly modified to allow for the fact that 𝐾 components correspond to 𝐾−1 changepoints, and we obtain𝑅_𝑏=𝑃(𝐾^𝑖+1)𝑃(𝐾^𝑖)2𝐾^𝑖(2𝐾^𝑖+1)(𝑚−2)²(𝑏_𝑗+1−𝑏^†)(𝑏^†−𝑏_𝑗)(𝑏_𝑗+1−𝑏_𝑗),𝐵_𝑏=𝑑_(𝐾^𝑖₊₁₎(𝑚−2)𝑏_𝐾^𝑖𝐾^𝑖. (12) For 𝐾^𝑖=𝐾^⋆ the birth of a new changepoint is invalid and the Markov chain is left unchanged.

We note that the proposal probabilities 𝑏_𝐾 and 𝑑_(𝐾+1) for birth and death moves can be chosen as follows:𝑏_𝐾=𝑐⋅min{1,𝑃(𝐾+1)𝑑𝑃(𝐾)},_(𝐾+1)=𝑐⋅min{1,𝑃(𝐾)𝑃(𝐾+1)}, (13) where 𝑐 is a constant that can be chosen as large as possible subject to the constraint 𝑏_𝐾+𝑑_𝐾≤0.9 for 𝐾=1,…,𝐾^⋆. This choice yields both a certain acceptance rate of the MCMC sampling scheme [17] and a simple prior probability (Hastings) ratio. From 𝑏_𝐾^𝑖⋅𝑃(𝐾^𝑖)=𝑑_(𝐾^𝑖₊₁₎⋅𝑃(𝐾^𝑖+1) it follows that the ratio 𝑑_(𝐾^𝑖₊₁₎/𝑏_𝐾^𝑖 in the expression 𝑅_𝑏 cancels out against the prior ratio 𝑃(𝐾^𝑖+1)/𝑃(𝐾^𝑖) in the expression 𝐵_𝑏, and the prior probability ratio simplifies to𝑅_𝑏𝐵_𝑏=2(2𝐾^𝑖+1)(𝑚−2)(𝑏_𝑗+1−𝑏^†)(𝑏^†−𝑏_𝑗)(𝑏_𝑗+1−𝑏_𝑗). (14)

(iii) A changepoint death move (𝑑) is the reverse of the birth move, and we obtain𝑅_𝑑𝐵_𝑑=(𝑚−2)2(2𝐾^𝑖−1)(𝑏_𝑗+1−𝑏_𝑗)(𝑏_𝑗+1−𝑏^†)(𝑏^†−𝑏_𝑗). (15)

3. Data

We have evaluated the proposed BGM_𝐷 model on various synthetic data sets. For illustration purposes, we present results obtained for two studies with small networks where the nonlinearity was implemented by a sinusoidal transformation. In a second study, we focus on the two-gene expression time series from macrophages and Arabidopsis thaliana, which have been used for evaluating the BGM model. More details on the experimental settings can be found in the paper on BGM [4]. In a third study, we switch from the Gaussian BGe score to the discrete multinomial BDe score to infer discrete nonstationary time series from Drosophila melanogaster.

3.1. Small Synthetic Networks

The first synthetic network consists of two nodes 𝑋 and 𝑌 and possesses two edges. Node 𝑋 has a recurrent feedback loop, symbolically𝑋(𝑡)→𝑋(𝑡+1), and 𝑋 acts as a regulator of node 𝑌, symbolically 𝑋(𝑡)→𝑌(𝑡+1). We consider the scenario of a nonlinear regulatory influence that 𝑋 exerts on 𝑌, whereby we implement the nonlinearity by a sinusoid transformation of 𝑋(𝑡)'s signal on 𝑌(𝑡+1). The state-space equations are given by𝑋(𝑡)=𝑋(𝑡−1)+𝑐+𝑐_𝑋⋅𝜀_𝑋(𝑡),𝑌(𝑡)=sin(𝑋(𝑡−1))+𝑐_𝑌⋅𝜀_𝑌(𝑡), (16) where 𝑐,𝑐_𝑋, and 𝑐_𝑌 are constants, 𝜀_𝑋(.),𝜀_𝑌(.) are iid normally distributed random variables.

The second synthetic network is a generalisation of the two node domain where three nodes 𝑌₁, 𝑌₂, and 𝑌₃ are regulated by 𝑋. The relationships are again realised by sinusoids, whereby we shift the periods: 𝑌_𝑖(𝑡)=sin(𝑋(𝑡−1)+𝜏_𝑖⋅𝜋)+𝑐_𝑌⋅𝜀_𝑌,𝑖(𝑡) with 𝜏₁=0, 𝜏₂=2/3, and 𝜏₃=4/3. For both networks we set the drift term 𝑐=2𝜋/40 to ensure that (on average) the complete period [0,2𝜋] of the sinusoid is involved, and we generate 𝑚=41 observations for four different combinations of the coefficients 𝑐_𝑋 and 𝑐_𝑌.

3.2. Bone Marrow Derived Macrophages

Interferons (IFNs) play a pivotal role in the innate and adaptive mammalian immune response against infection, and central research efforts therefore aim to elucidate their regulatory interactions [19]. For the present study, we have analysed gene expression time series from bone marrow-derived macrophages, which were sampled at 𝑚=25×30 minute time intervals. The macrophages were subjected to three external conditions: (1) infection with Cytomegalovirus (CMV), (2) treatment with Interferon Gamma (IFN𝛾), and (3) infection with Cytomegalovirus after pretreatment with IFN𝛾 (CMV+IFN𝛾). Samples derived from the macrophages were hybridised to Agilent mouse genome arrays. We focus on 𝑁=3 Interferon-regulatory factors Irf1, Irf2, and Irf3. These factors are the key regulators in the response of the macrophage cell to pathogens. The macrophages data sets used in the study are publicly available from http://www.bioss.ac.uk/associates/marco/supplement/.

3.3. Circadian Regulation in Arabidopsis thaliana

We have also analysed two-gene expression time series from Arabidopsis thaliana cells, which were sampled at 𝑚 = 13 × 2 hour time intervals with Affymetrix microarray chips. The expressions were measured twice independently under experimentally generated constant light condition, but differed with respect to the prehistories. In the first experimental scenario, 𝑇₂₀, the plants were entrained in a 10 h : 10 h light/dark-cycle, while the plants in the second experimental setting, 𝑇₂₈, were entrained in 14 h : 14 h light/dark-cycle. The analysis focuses on 𝑁=9 genes, namely LHY, CCA1, TOC1, ELF4, ELF3, GI, PRR9, PRR5, and PRR3, which are known to be involved in circadian regulation [20, 21]. The Arabidopsis data sets used in the study are publicly available from http://www.bioss.ac.uk/associates/marco/supplement/.

3.4. Muscle Development in Drosophila melanogaster

The gene expressions in Drosophila melanogaster cells were sampled at 𝑚=67 time-steps during four different morphological stages of life: embryonic, larval, pupal, and adult stages. Since these phases cover time periods of different lengths, gene expression profiles were collected at nonequidistant time-points. The true morphological transitions occur at time points 𝑡₃₂ (embryonic to larval), 𝑡₄₂ (larval to pupal), and 𝑡₆₀ (pupal to adult) [22]. Like other researchers [12], we focus our analysis on 𝑁=13 genes involved in muscle growth and muscle development: EVE, GFL, TWI, MLC1, SLS, MHC, PRM, ACTN, UP, MYO61F, and MSP300. The quantile-discretised binary data set of these 13 Drosophila genes is available from Robinson and Hartemink [12].

4. Simulations and Evaluation

We impose a fan-in restriction of size ℱ=3 on the maximal number of parent nodes per node as done in related articles [2, 4]. For the Gaussian BGe model the hyperparameters of the normal-Wishart prior were chosen maximally uninformative subject to certain regularity constraints [14]. In the second study, the hyperparameters of the discrete multinomial BDe model were chosen as explained in Heckerman and Geiger [16], whereby the total prior precision parameter was set to 1.

We set 𝐾^⋆=10, and the burn-in and the sampling-phase of MCMC runs were set to 500,000 iterations each, and we sampled every 1,000 iterations during the sampling-phase. For each data set, we started 5 independent MCMC simulations from different initialisations, and we computed the potential scale reduction factor PSRF [23] based on the marginal edge posterior probabilities. As we observed, a sufficient degree of convergence for all data sets (PSRF<1.2), we report only the results of the empty-seeded (graphs without any edges) MCMC runs. We note that each single MCMC simulation (even for the bigger domains Arabidopsis with 𝑁=9 genes and Drosophila with 𝑁=13 genes) was accomplished within a few hours using Matlab^© code on a SunFire X4100M2 machine with MAD Opteron 2224 SE dual-core processor.

More generally, we note that the computational complexity in network structure space is 𝑁∑^ℱ_𝑘=0(^𝑁_𝑘)∼𝑂(𝑁^ℱ+1), where 𝑁 is the number of nodes and ℱ is the maximum number of parent nodes. The complexity in changepoint configuration space is of order 𝑚^𝐾, where 𝐾 is the number of changepoints, and 𝑚 is the number of time points. Hence, the problem is of polynomial complexity in 𝑁 and 𝑚, and not exponential complexity (i.e., it is not NP-hard). Polynomial complexity—as opposed to exponential complexity—does not impose any principled restrictions on the network size. However, the practical feasibility will depend on varying factors, like the efficiency of the software implementation and the capacity of parallel clusters. It also depends on the information content of the data, as increasing the number of nodes with limited number of experimental replications will increase the intrinsic uncertainty of inference. The practical decision on how many nodes and how large a network to consider will therefore usually be based on some preliminary data analysis.

We note that we allowed for self-feedback loops for the synthetic data only. For the three real applications to biological systems, we ruled out self-feedback loops altogether to enable direct comparability with the results reported for the BGM model [4].

If two independent data sets 𝐃 and 𝐃_‡ are available for a network domain, predictive probabilities of the models BGM and BGM_𝐷 can be estimated straightforwardly [4]. For example, having sampled {𝒢^𝑖,𝐕^𝑖,𝐾^𝑖}_{𝑖=1,…,𝐼} from the posterior distribution 𝑃(𝒢,𝐕,𝐾|𝐃) via the MCMC inference scheme described above, the predictive probability of BGM_𝐷 can be estimated by𝑃(𝐃_‡∣𝐃,BGM_𝐷1)=𝐼^𝐼𝐾_𝑖=1^𝑖_𝑘=1Ψ(𝒢^𝑖,𝐃^(𝐕𝑖_‡^,𝑘)∣𝐃^{(𝐕𝑖,𝑘)}), (17) where Ψ(𝒢^𝑖,𝐃^(𝐕𝑖_‡^,𝑘)∣𝐃^{(𝐕𝑖,𝑘)})=^𝑁_{𝑛=1𝑡∶𝐕𝑖(𝑡)=𝑘}𝑃(𝑋_𝑛(𝑡)=𝐃_‡𝑛,𝑡∣𝜋_𝑛(𝑡−1)=𝐃_{‡(𝜋𝑛,𝑡−1)},𝜽_𝑛)×𝑃(𝜽_𝑛∣𝜋_𝑛,𝐃^{(𝐕𝑖,𝑘)})𝑑𝜽_𝑛. (18) That is, when computing the BGe score for the compartment 𝐷^{(𝐕𝑖,𝑘)} with (6) and (7), the local prior distributions 𝑃(𝜽_𝑛|𝜋_𝑛) in (7) are replaced by posterior distributions 𝑃(𝜽_𝑛|𝜋_𝑛,𝐃^{(𝐕𝑖,𝑘)}). This results in a straightforward modification of the BGe score as follows. In (13) in Geiger and Heckerman [14], those training data that have been allocated to component 𝑘, symbolically 𝐃^(𝐕,𝑘), are included in the conditioning part of the distribution, and the sufficient statistics are adjusted accordingly.

5. Results

5.1. Inference on Synthetic Data

For the synthetic sinusoid data, the true underlying network topologies are known so that the network reconstruction accuracy can be assessed via the area under the ROC (receiver operator characteristic) curve: AUC; this is a standard criterion that has been applied in numerous related articles (e.g., [24]). Figure 1 shows histograms of the average AUC scores obtained from the marginal edge posterior probabilities for the sinusoid networks with 𝑁=2 and 𝑁=4 nodes. Figure 1 is laid out as matrix, in which the rows and columns correspond to different 𝑐_𝑋 and 𝑐_𝑌 coefficients (noise levels). We note that an increase of 𝑐_𝑋 reduces the autocorrelation of 𝑋, while increasing 𝑐_𝑌 blurs the functional dependence of 𝑌(𝑡+1) on 𝑋(𝑡). The autocorrelation of 𝑌 is jointly influenced by both parameters. From the histograms it can be seen that the novel BGM_𝐷 model leads to a better network reconstruction accuracy in terms of average AUC values than the standard BGe model and the BGM model in the majority of scenarios. Further investigations showed that BGe and BGM yield lower AUC scores, since they tend to infer spurious self-feedback loops on node 𝑌 (𝑁=2) or on the nodes 𝑌_𝑖 (𝑁=4), respectively.

(a)

(b)

Figure 1 

AUC histograms—Cross-method comparison on synthetic sine data. AUC scores for the synthetic sine network data with 𝑁=2 nodes (a) and 𝑁=4 nodes (b). The figure is laid out as a matrix, where rows and columns correspond to different noise levels 𝑐𝑋 (rows) and 𝑐𝑌 (columns). In each histogram, the white bar shows the average AUC score for the BGe model. The grey bar shows the average AUC score of the BGM model, and the black bar shows the AUC score for the proposed BGM𝐷 model. Each histogram shows averages and standard deviations obtained from 50 independent data instantiations.

This trend can be visualised by histograms of the average edge posterior probabilities. As an example, Figure 2 shows the average marginal edge posterior probabilities of the four possible edges for the 𝑁=2 nodes sinusoid network with 𝑐_𝑋=0.5 and 𝑐_𝑌=0.5. Consistently, all three models under comparison assign the highest posterior probability to the true self loop 𝑋→𝑋 and the lowest posterior probability to the false edge 𝑌→𝑋. But BGe and BGM favour the spurious feed-back loop 𝑌→𝑌 over the true edge 𝑋→𝑌 while the proposed BGM_𝐷 suppresses the false self-feedback loop and assigns a higher edge posterior probability to the true edge 𝑋→𝑌. This shows that BGM_𝐷 yields a higher network reconstruction accuracy (see Figure 1), as it is less susceptible to inferring spurious self-feedback loops (see Figure 2).

(a)

(b)

(c)

Figure 2 

Edge Posterior Probabilities—Cross-method comparison on synthetic sine data. The figure shows three histograms of the inferred marginal edge posterior probabilities in the sinusoid network with 𝑁=2 nodes and 𝑐𝑋=0.5 and 𝑐𝑌=0.5 as obtained with BGe (a), BGM (b), and BGM𝐷 (c). In each histogram, the four bars represent the four possible edges: Left: self-loop 𝑋→𝑋 (true); centre left: 𝑋→𝑌 (true); centre right: self-loop 𝑌→𝑌 (false); right: 𝑌→𝑋 (false). Each bar shows the average marginal posterior probability, averaged over 50 independent data instantiations. It is seen that BGe and BGM have a high propensity for learning the spurious feedback loop 𝑌→𝑌 (centre right white bars). BGM𝐷 (right histogram) assigns a higher posterior probability to the correct edge 𝑋→𝑌 (centre left black bar) and suppresses the spurious feedback loop 𝑌→𝑌 (centre right white bar)

5.2. Inference on Macrophages Data

For the macrophages data the BGM model inferred a biologically plausible state change in the host macrophage brought about by infection (CMV) or immune activation (IFN𝛾), and a less pronounced state change in the combined condition CMV+IFN𝛾 [4]. We compare these findings with results obtained with the novel BGM_𝐷 model. The fraction of sampled states for which two time points 𝑡_𝑖 and 𝑡_𝑗 are allocated to the same component 𝑘 (1≤𝑘≤𝐾) can be used as a connectivity measure 𝐶(𝑡_𝑖,𝑡_𝑗), and the resulting temporal connectivity structures are displayed graphically as heat maps in Figure 3. All six heatmaps in Figure 3 reflect the two-stage nature of the gene-regulatory processes in the host macrophages: the first part (time points (𝑡₂,…,𝑡₆)) and the last part of the three time series (time points (𝑡₇,…,𝑡₂₅)) are allocated to different components. For all three conditions, a stronger separation between the two regulatory states is inferred by the BGM_𝐷 model (see Figures 3(d), 3(e), and 3(f)). It appears that the BGM_𝐷 inference results are more consistent, as even for the combined condition (CMV+IFN𝛾) a clear trend towards a dichotomous regulatory process can be found (see Figure 3(d)). This finding (stronger separation) is consistent with our conjecture that the novel BGM_𝐷 assigns neighbouring time-points to the same compartment more likely a priori. Interestingly, the BGM inference outlier at time point 𝑡₉ in Figure 3(b) yields a certain trend for a subdivision of the second compartment (𝑡₇,…,𝑡₂₅) by the BGM_𝐷 model. Instead of one outlying time point two substages (𝑡₇,…,𝑡₁₀) and (𝑡₁₁,…,𝑡₂₅) are inferred (see Figure 3(e)). To provide statistical evidence that the new BGM_𝐷 model does not overfit the data, we compute predictive probabilities for the BGM_𝐷 model and compare them with those reported for the BGM model [4].

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

Heat maps—macrophages data. Graphical heat map presentation of the temporal connectivity structure for the macrophage gene expression time series. (a), (b), and (c): Heat matrices for experiments CMV (a), IFNG𝛾 (b), and CMV+IFN𝛾 (c) inferred with the BGM model. (d), (e), and (f): Heat matrices for experiments CMV (d), IFNG𝛾 (e), and CMV+IFN𝛾 (f) inferred with the novel BGM𝐷 model. Each heat map indicates the estimated posterior probability of two time points being assigned to the same compartment (mixture component). The probabilities are represented by a grey shading, where white corresponds to a probability of 1, and black corresponds to a probability of 0. The numbers on the axes represent the time points of the time course experiment.

To this end, we treat the three experiments as independent replications, and check via predictive probabilities whether the superiority of the proposed BGM_𝐷 model can be confirmed statistically. Table 1 gives the predictive probabilities reported in the BGM paper [4] and those obtained with the proposed BGM_𝐷 model. As the predictive probabilities for BGM_𝐷 are systematically better than those of BGM, we conclude that the BGM_𝐷 model yields more stable inference results, that is, a better generalisation performance.

Table 1 

Predictive Probabilities: Macrophages data. Logarithmic predictive probabilities for the macrophage data: log_𝑒(𝑃(𝐃_‡|𝐃)) for BGe, BGM (as reported earlier [4]) and the new BGM_𝐷 model. The standard deviations of the logarithmic probabilities are given in brackets. We note that the BGM values could be confirmed by our reanalysis of the data: the deviations were smaller than one standard deviation.

5.3. Inference on Arabidopsis thaliana Data

For the Arabidopsis thaliana data, the BGM model also inferred a biologically plausible two-stage process [4]. In this application, the two stages are likely to be related to the diurnal nature of the dark-light cycle influencing the circadian genes. The plants were subjected to different prehistories, related to different lengths of the artificial, experimentally controlled light-dark cycle. The plants in experimental scenario 𝑇₂₈ were entrained in an increased day length of 14 hours light followed by 14 hours darkness and in experiment 𝑇₂₀ the plants were entrained in a decreased day length of 10 hours light followed by 10 hours darkness. As an effect of these two entrainments, a phase shift in the gene-regulatory processes between these two experiments was expected [4]. The BGM model inferred a certain trend for a phase shift of the changepoint (subjective day to subjective night) of about 4–6 hours as a consequence of the increased day length. The heat maps in Figures 4(a) and 4(b) show that the connected blocks (compartments) of the time series are shifted along the diagonal by 2-3 time-points (4–6 hours). The BGM_𝐷 model infers the same trend but with a stronger separation score between these compartments (see Figures 4(c) and 4(d)). We note that the BGM_𝐷 model is based on changepoints so that compartments once left cannot be revisited. That is, while the BGM model tends to allocate the first time points (𝑡₂,𝑡₃) and the last time points (𝑡₉,…,𝑡₁₃) in experiment T₂₈ to one single component (light grey shading in the top right and bottom left area of the heat map in the top centre panel of Figure 4), the BGM_𝐷 model has to allocate the last time points (𝑡₉,…,𝑡₁₃) to an additional third component, as the first compartment (𝑡₂,𝑡₃) cannot be reused after the transition to the second compartment (𝑡₄,…,𝑡₈).

(a)

(b)

(c)

(d)

Figure 4 

Heat maps. Arabidopsis data. Graphical heat map representations of the temporal connectivity structures for the Arabidopsis thaliana data. (a) and (b): heat matrices for experiments 𝑇20 (a) and 𝑇28 (b) inferred with the BGM model. (c) and (d): heat matrices for experiments 𝑇20 (c) and 𝑇28 (d) inferred with the novel BGM𝐷 model. Each heat map indicates the posterior probability of two time points being assigned to the same compartment (mixture component). The probabilities are represented by a grey shading, where white corresponds to a probability of 1, and black corresponds to a probability of 0. The numbers on the axes represent the time points of the time course experiment.

As for the macrophages data, predictive probabilities can be computed by treating the two experiments 𝑇₂₀ and 𝑇₂₈ as independent replications. Table 2 gives the predictive probabilities reported for the BGM model [4] and those obtained with the new BGM_𝐷 model. In consistency with the results for the macrophages data, the resulting two predictive probabilities for BGM_𝐷 are better than those of the BGM model. A scatter plot 𝐸_𝑖,𝑗(𝑇₂₀) versus 𝐸_𝑖,𝑗(𝑇₂₈) of the inferred (marginal) posterior probabilities of the individual edges for the BGM_𝐷 model inference is given in Figure 5. As the individual edge posterior probabilities 𝐸_𝑖,𝑗 do not differ drastically, we extract a network structure from the averaged probabilities: 𝐸_𝑖,𝑗=(𝐸_𝑖,𝑗(𝑇₂₀)+𝐸_𝑖,𝑗(𝑇₂₈))/2.

Table 2 

Predictive probabilities: Arabidopsis thaliana data. Logarithmic predictive probabilities for the Arabidopsis thaliana data: log_𝑒(𝑃(𝐃_‡|𝐃)) for BGe, BGM (as earlier reported [4]) and the new BGM_𝐷 model. The standard deviations of the logarithmic probabilities are given in brackets. We note that the BGM values could be confirmed by our reanalysis of the data: the deviations were smaller than one standard deviation.

Figure 5 

BGM𝐷 scatter plot. Arabidopsis data. Scatter plot of the marginal edge posterior probabilities inferred with the proposed BGM𝐷 model. In the plot the marginal edge posterior probabilities for time series 𝑇20: 𝐸𝑖,𝑗(𝑇20) (horizontal axis) are plotted versus the marginal edge posterior probabilities for time series 𝑇28: 𝐸𝑖,𝑗(𝑇28) (vertical axis).

Figure 6 shows the gene interaction network that is predicted when keeping all edges with marginal posterior probability above 0.5. There are two groups of genes. Empty circles in the figure represent morning genes (i.e., genes whose expression peaks in the morning), shaded circles represent evening genes (i.e., genes whose expression peaks in the evening). There are several directed edges pointing from the group of morning genes to the evening genes such that each evening gene is regulated by at least one morning gene. Moreover, the two genes LHY and CCA1 seem to play a central role. This result is consistent with biological findings, where the morning genes were found to activate the evening genes, with LHY and CCA1 being central regulators [25]. Our reconstructed network also contains edges pointing into the opposite direction, from the evening genes back to the morning genes. This finding is consistent with biological observations [25], where the evening genes were found to inhibit the morning genes via negative feedback loops. In the reconstructed network, there are 9 edges (drawn in black) originating from the four morning genes while only 7 edges (drawn in grey) originate from the group of five evening genes. Biologically, this means that the activity of the morning genes is stronger than the activity of the group of evening genes and that the regulatory mechanisms are dominated by the morning genes in the network topology. This finding is consistent with the fact that following the light-dark cycle entrainment, the experiments were carried out in constant-light condition, resulting in a higher activity of the morning genes overall. Within the group of evening genes, the reconstructed network contains an edge between GI and TOC1. This interaction has been confirmed independently [26]. Hence, while a proper evaluation of the reconstruction accuracy is currently unfeasible—like many related studies, we lack a gold-standard owing to the unknown nature of the true interaction network—our study suggests that the essential features of the reconstructed network are biologically plausible and consistent with the literature.

Figure 6 

Predicted regulatory network of nine circadian genes in Arabidopsis thaliana. From the averaged marginal edge posterior probabilities (average of 𝐸𝑖,𝑗(𝑇20) and 𝐸𝑖,𝑗(𝑇28)) of the proposed BGM𝐷 model inference results a regulatory network can be extracted. Empty circles represent morning genes. Shaded circles represent evening genes. Edges indicate predicted interactions with an inferred marginal posterior probability greater than 0.5. Edges are black (grey) if they refer to a morning (evening) gene as regulator.

5.4. Inference on Drosophila melanogaster Data

For an independent comparison of the free allocation model and the changepoint model, we carried out an analysis similar to Robinson and Hartemink [12] on the binary Drosophila muscle development gene expression time series. This time series can be analysed with the discrete counterparts BDM and BDM_𝐷 of our Gaussian Mixture models; see Section 2.3 for details. The graphical heat map representations in Figure 7 show that the BDM model does not infer the morphological stages of Drosophila melanogaster. Almost all time-points are strongly connected (white shading) and no separated blocks of connected time points have been inferred. Only a few time-points are allocated (as outliers) to other mixture components (black vertical and horizontal lines).

(a)

(b)

Figure 7 

Heat maps for the discrete Drosophila data. Graphical heat map representations of the temporal connectivity structures for the Drosophila melanogaster data. (a) Heat matrix inferred with the BDM model. (b) Heat matrix inferred with the novel changepoint variant BDM𝐷. Each heat map indicates the posterior probability of two time points being assigned to the same compartment (mixture component). The probabilities are represented by a grey shading, where white corresponds to a probability of 1, and black corresponds to a probability of 0. The numbers on the axes represent the time points of the time course experiment. Both axes have been ticked at the three real morphological stage transitions: embryonic to larval (𝑡31→𝑡32), larval to pupal (𝑡41→𝑡42), and pupal to adult (𝑡59→𝑡60).

The novel BDM_𝐷 infers a clear block structure with different separated compartments. The connectivity structure corresponds well to the first two morphological transitions (i) embryonic to larval (𝑡₃₁→𝑡₃₂) and (ii) larval to pupal (𝑡₄₁→𝑡₄₂), whereby the separation between the embryonic and the larval stage is less pronounced (grey shading) than the separation between the larval and the pupal stage (black shading). The exact third morphological transition pupal to adult (𝑡₅₉→𝑡₆₀) has not been inferred but it can be seen that two changepoints occur during the pupal stage before the third transition to the adult stage. Figure 8 shows a graphical presentation of the changepoint location posterior probabilities inferred with BDM_𝐷 and the same trends become obvious: the first two stage transitions have been inferred correctly while two further changepoints occur before the morphological stage transition from pupal to adult. We note that Robinson and Hartemink did find the same trends with their 𝑛𝑠𝐵𝐷𝑒 model and they conclude that the third premature transition in the gene regulatory process can be explained biologically, since the gene expression program governing the transition from pupal to adult morphology should be active well before the time of the real morphological transition [12].

Figure 8 

BDM𝐷 transition posterior probabilities for Drosophila data. Graphical representation of the posterior probabilities of the changepoint locations inferred with the novel BDM𝐷. The transition posterior probabilities (vertical axis) are plotted against the time axis (horizontal axis). The time axis has been ticked at the three real morphological stage transitions: embryonic to larval (𝑡31→𝑡32), larval to pupal (𝑡41→𝑡42), and pupal to adult (𝑡59→𝑡60). We note that the BDM model is based on a free individual allocation of time points so that a transition posterior probability plot cannot be interpreted properly.

6. Discussion

Our empirical results have shown that the proposed BGM_𝐷 model performs consistently better than the BGM model. A possible explanation could be related to the latent space complexity of the models. For two components, we have 𝑚−2 changepoint locations, but 2^𝑚−1 free allocations. This difference in latent space complexity gets aggravated for more components. The difference in performance between BGM_𝐷 and BGM could therefore be a consequence of the different degrees of convergence of the MCMC simulations. However, our convergence diagnostics based on potential scale reduction factors [23] did not indicate any significant difference in the convergence. It therefore appears that the higher latent space complexity of the BGM model does not adversely affect the MCMC convergence for gene expression time series of the length investigated in our study. This suggests that another explanation for the better performance of BGM_𝐷 over BGM has to be found.

We will discuss that the performance difference is most likely a consequence of the different prior probabilities intrinsic to the models, which determine the factor 𝑅 in (10). Since both models BGM and BGM_𝐷 employ the same Poisson distribution with parameter 𝜆=1 truncated to the interval [1,𝐾^⋆] as prior for the number of mixture components, the difference between the two models BGM and BGM_𝐷 is imposed by the prior distribution of the allocation vector 𝑃(𝐕|𝐾). While the BGM model is based on a free allocation, the BGM_𝐷 model takes the time structure into account and employs a changepoint process. In this section we describe the differences between the two models in detail by three theoretical examinations, and these examinations reveal trends that appear to be immediately related to some of our empirical findings.

In the first study, we compare the (temporal) connectivity structures that are introduced by the prior 𝑃(𝐕|𝐾). To this end, we infer the prior distribution 𝑃(𝐾)𝑃(𝐕|𝐾) of both models BGM and BGM_𝐷 via MCMC simulations that are purely prior-driven. That is, we employ an empty data set (without any data points) so that all acceptance probabilities depend on the prior probability ratios only, as 𝑃(𝐃=∅|(𝒢,𝐕,𝐾))=1 for all combinations of 𝒢, 𝐕, and 𝐾. Note that sampling network structures via edge operation moves on the graph 𝒢 becomes obsolete, because (i) the graphs do not have any effect on the likelihoods for an empty data set, and (ii) the graph priors cancel out in the prior probability ratio if a uniform graph prior 𝑃(𝒢)=const is used. We set 𝑚=26; note that in this theoretical consideration based on empty data, 𝑚 only determines the length of the allocation vector.

After running 10 independent MCMC simulations with 𝑚=26 to infer the prior distribution 𝑃(𝐾)⋅𝑃(𝐕|𝐾), we can compute the average prior connectivity strengths from the sampled allocation vectors. As before, the fraction of sampled allocation vectors for which two time points 𝑡_𝑖 and 𝑡_𝑗 are allocated to the same component𝑘 (1≤𝑘≤𝐾) can be used as a connectivity measure 𝐶(𝑡_𝑖,𝑡_𝑗). Figure 9 shows heat maps of the inferred prior connectivity structure for the free allocation BGM model and the proposed changepoint process BGM_𝐷 model. The heat maps confirm our earlier conjecture that the proposed BGM_𝐷 model, which takes the time structure of the data into account, allocates neighbouring time points to the same compartment more likely a priori than the BGM model. More precisely, it can be seen from Figure 9(a) that the prior connectivity strength 𝐶(𝑡_𝑖,𝑡_𝑗) is the same for all 𝑡_𝑖 and 𝑡_𝑗 with 𝑡_𝑖≠𝑡_𝑗 in the BGM model. On the contrary, for the BGM_𝐷 model (Figure 9(b)), the connectivity strength 𝐶(𝑡_𝑖,𝑡_𝑗) decreases with the temporal distance between 𝑡_𝑖 and 𝑡_𝑗: for three time points 𝑡_𝑖, 𝑡_𝑗 and 𝑡_𝑘 with 𝑡_𝑖<𝑡_𝑗<𝑡_𝑘 we have 𝐶(𝑡_𝑖,𝑡_𝑗)<𝐶(𝑡_𝑖,𝑡_𝑘). This finding explains why the proposed BGM_𝐷 model yields a stronger separation between the light : darkness induced stages in Arabidopsis thaliana (see Figure 4): the two-stage structure of gene-regulation in Arabidopsis is of a temporal form that is supported by the allocation vector prior 𝑃(𝐕|𝐾) of the BGM_𝐷 model.

(a)

(b)

Figure 9 

Prior connectivity structure for time series of length 𝑚=26. Graphical heat map representations of the temporal connectivity structures imposed by the prior distribution 𝑃(𝐾)⋅𝑃(𝐕∣𝐾) for 𝑚=26 time points. (a) Free allocation model (BGM). (b) Changepoint process (BGM𝐷). Each heat map indicates the prior probability of two time points being assigned to the same compartment. The probabilities are represented by a grey shading, where white corresponds to a probability of 1, and black corresponds to a probability of 0.5. The connectivity strengths were estimated from 10 independent MCMC simulations. In these simulations, an empty data set (without any data points) was used so that the inference was driven exclusively by the prior probability distribution 𝑃(𝐾)⋅𝑃(𝐕∣𝐾).

In the second and third theoretical study, we cross-compare the BGM model and the proposed BGM_𝐷 model in terms of prior probability ratios between the heterogeneous (nonstationary) and the homogeneous (stationary) state. That is, both models can be in the homogeneous state (i.e., the complete time series is modelled as one single compartment (𝐾=1)) or in a heterogeneous state (i.e., the time series is divided into 𝐾 different compartments (𝐾>1)). For the BGM as well as the BGM_𝐷 model the prior probability ratio 𝑅 between (i) the heterogeneous state consisting of two segments (𝐾=2) and (ii) the homogeneous state (𝐾=1) is given by𝑅=𝑃(𝐾=2)⋅𝑃(𝐾=1)𝑃(𝐕∣𝐾=2)1, (19) since 𝑃(𝐕|𝐾=1)=1 for both models. We consider the scenario where the time series is of length 𝑚 and the allocation vector 𝐕 divides the time series into two (non-empty) connected segments 𝑡₂,…,𝑡_𝑗 and 𝑡_𝑗+1,…,t_𝑚 (2≤𝑗≤𝑚). The prior probability ratio of the BGM model is then given by𝑅_BGM=𝑃(𝐾=2)⋅𝑃(𝐾=1)Γ(𝛼₀)Γ(𝛼₀⋅+(𝑚−1))²_𝑘=1Γ(𝛼_𝑘+𝑚_𝑘)Γ(𝛼_𝑘), (20) where 𝛼₁=𝛼₂=1, 𝛼₀=𝛼₁+𝛼₂=2, and 𝑚_𝑘 is the number of time points that have been allocated to the 𝑘th segment (𝑘=1,2); that is, 𝑚₁=𝑗−1 and 𝑚₂=𝑚−𝑗 in our scenario. See [4] for details. The proposed BGM_𝐷 model requires a changepoint in the interval [𝑗,𝑗+1] and it can be derived straightforwardly from (12):𝑅_BGM𝐷=𝑃(𝐾=2)𝑃(𝐾=1)⋅_𝑗^𝑗+16(𝑚−𝑏₁)(𝑏₁−2)(𝑚−2)³𝑑𝑏₁. (21) In the second theoretical study we vary the length of the time series 𝑚=3,5,7,…,25, and we consider a heterogeneous time series consisting of two equally-spaced segments 𝑡₂,…,𝑡_{⌊𝑚/2⌋+1} and 𝑡_{⌊𝑚/2⌋+2},…,𝑡_𝑚. This corresponds to 𝑚₁=𝑚₂=0.5⋅(𝑚−1) in the BGM model. For the BGM_𝐷 model, we obtain with 𝑗=⌊𝑚/2⌋+1 that the changepoint has to be located in the interval 𝑏₁∈[𝑡_{⌊𝑚/2⌋+1},𝑡_{⌊𝑚/2⌋+2}]. Figures 10(a) and 10(b) show the resulting (logarithmic) prior probability ratios in dependence on 𝑚. It can be seen that the prior ratio 𝑅 for the BGM model is considerably lower than for the BGM_𝐷 model. Moreover, the logarithmic plot in Figure 10(b) shows that the prior ratio of the BGM model shows a much stronger decrease with the sample size 𝑚 than the BGM_𝐷 model. This suggests that the BGM model imposes a more severe penalty for complexity (non-stationarity), which increases with increasing sample size 𝑚. This tendency may explain the finding in [4] for the macrophage gene expression time series, which we have reproduced in the present study (Figure 3(c)): the BGM model does not infer a clear two-phase nature of the time series under simultaneous immune activation (with IFN𝛾) and viral infection (with CMV). A possible biological explanation was offered in [4]. However, the novel BGM model does not support the hypothesis of a decreased probability for the two-phase nature (Figure 3(f)). Moreover, the previous analysis has revealed that a strong penalty against the two-phase process is inherent in the BGM model. This suggests that the results reported in [4], which we have reproduced in our study, might be an artefact of the BGM model rather than of genuine biological nature.

(a)

(b)

(c)

(d)

Figure 10 

Prior probability ratios between the heterogeneous and the homogeneous state for (i) varying time series length 𝑚 ((a) and (b)) and (ii) varying segment length proportions ((c) and (d)). (a) and (b): prior probability ratio 𝑅 between (i) the heterogeneous state that consists of two equally-spaced segments 𝑡2,…,𝑡⌊𝑚/2⌋+1 and 𝑡⌊𝑚/2⌋+2,…,𝑡𝑚 and (ii) the homogeneous state consisting of one single segment 𝑡2,…,𝑡𝑚. The prior probability ratios (vertical axis) are plotted in dependence on the time series length 𝑚=3,5,7,…,25 (horizontal axis). The prior probability ratio was defined in (19). For clarity the logarithmic prior probability ratios are plotted in (b). (c) and (d): prior probability ratio 𝑅 for a time series of length 𝑚=26 between (i) the heterogeneous state with two segments 𝑡2,…,𝑡𝑗 and 𝑡𝑗+1,…,𝑡𝑚 and (ii) the homogeneous state consisting of one single segment 𝑡2,…,𝑡𝑚 only. The prior probability ratios (vertical axis) are plotted in dependence on the changepoint location (horizontal axis). For the sake of clarity the logarithmic prior probability ratios are plotted in (b). See text for further details.

In the third theoretical study, we fix the length of the time series 𝑚=26 and we vary the last time point 𝑗=2,3,…,25 of the first segment in the heterogeneous state to illustrate the effect of unequal segment lengths 𝑗−1 (𝑡₂,…,𝑡_𝑗) and 𝑚−𝑗 (𝑡_𝑗+1,…,𝑡_𝑚). For the BGM model, this scenario yields 𝑚₁=𝑗−1 and 𝑚₂=𝑚−𝑗, and for the BGM_𝐷 model we have the changepoint location intervals 𝑏₁∈[𝑗,𝑗+1], 𝑗=2,3,…,25. Figures 10(c) and 10(d) show the resulting prior probability ratios in dependence on the changepoint location 𝑗. The figure reveals contrary trends for the BGM and the BGM_𝐷 model. Figure 10(c) shows that the prior ratio of the BGM_𝐷 model peaks in the middle, that is for equally long time series segments, whereas it decreases monotonically with increasing asymmetry of the segment lengths. The BGM model (Figure 10(d)) exhibits the converse behaviour: the more disproportionate the segment lengths, the higher the prior ratio 𝑅. This behaviour becomes even more obvious when varying the time series length 𝑚 and the breakpoint location 𝑗 simultaneously. Figure 11 shows the logarithmic prior probabilities in dependence on both: the time series length 𝑚 and the segment length proportions. For the BGM model (Figure 11(a)) the penalty term 𝑃(𝐕|𝐾) increases drastically with the sample size 𝑚 if both segments are of a certain length. The higher the sample size the stronger the penalty for symmetric segment lengths. Only if 𝑗 is either very low or very high, that is if the segment lengths are strongly disproportionate, does the size of the prior probability ratio 𝑅 not change drastically with the sample size 𝑚. Figure 11(b) shows that the prior probability ratio of the proposed BGM_𝐷 model is less sensitive to both: the sample size 𝑚 and the segment length proportions. We can conclude that it is the equability of the prior ratio of the proposed BGM_𝐷 that renders it superior when modelling nonstationary behaviour in long time series. That is, the penalty for dividing a time series into segments does not change drastically with the length of the time series 𝑚, and symmetric segment lengths are supported by the prior distribution 𝑃(𝐕|𝐾). On the contrary, the penalty term of the BGM model for dividing a time series into segments increases substantially with the length of the time series: the longer the time series, the lower the prior probability 𝑃(𝐕|𝐾). Only if the segment lengths are strongly asymmetric such that one segment is very long and the other very short, is the prior probability 𝑃(𝐕|𝐾) of a comparable size for different time series lengths 𝑚. This tendency provides a possible explanation for the failure of the BGM model on the Drosophila melanogaster gene expression time series. The heat map in Figure 7 shows that the proposed BGM_𝐷 model divides the time series into segments that are consistent with the morphogenesis of Drosophila melanogaster. The BGM model also tends to detect the correct segment boundaries. However, it then erroneously infers short segments consisting of a only few time points around these boundaries. This pattern is not consistent with the morphogenetical findings. The segmentation inferred with the BGM model thus suffers from artefacts that are an immediate consequence of what has been discussed above. Namely, that the BGM model penalises symmetric partitions much more strongly than asymmetric partitions and thereby encourages the formation of segments that consist of a single or only a few time points. Note that the BGM model becomes more susceptible to these artefacts as the length of the time series increases, as demonstrated in Figure 11(a). This explains why the difference between the BGM and the BGM_𝐷 model is less pronounced for the Arabidopsis and macrophage expression time series (Figures 3 and 4), which are considerably shorter (Arabidopsis: 13 time points, macrophages: 25 time points, Drosophila: 67 time points). The susceptibility of the BGM model to short-segment artefacts is not completely avoided here either, though, as can be seen from the heat map in Figure 3(b): While BGM_𝐷 infers three segments for the macrophages data under condition IFN_𝛾 (see Figure 3(e)), BGM tends to allocate a single data point to a separate segment.

(a)

(b)

Figure 11 

Prior probability ratios in dependence on both: the time series length m and the segment lengths proportions. Graphical representation of the prior probability ratio 𝑅 between (i) the heterogeneous state consisting of two segments 𝑡2,…,𝑡𝑗 and 𝑡𝑗+1,…,𝑡𝑚 and (ii) the homogeneous state consisting of one single segment 𝑡2,…,𝑡𝑚 only. The prior probability ratio was defined in (19). For three different lengths of the time series 𝑚=13,25,67 the logarithmic prior probability ratios (vertical axis) between these two states are plotted in dependence on the changepoint location (horizontal axis). (a) BGM model. (b) BGM𝐷 model.

On the synthetic sinusoid data, the proposed BGM_𝐷 model yields a higher network reconstruction accuracy than the BGM model, as the latter is more susceptible to inferring spurious self-feedback loops. This tendency can be explained from the previous mathematical analysis. Both the BGM and the BGM_𝐷 model effectively approximate the nonlinear function by a piecewise linear function. A good approximation of the sine wave requires three segments of approximately the same length, corresponding to the ascending, stationary and descending phase. As opposed to BGM_𝐷, the prior inherent to the BGM model heavily penalises against this equal-length segmentation (Figures 10 and 11); see the previous discussion. Now, from (16) it becomes clear that the data are strongly autocorrelated. More precisely, the 𝑌(𝑡)’s tend to exhibit a strong autocorrelation by virtue of the autocorrelation of the 𝑋(𝑡)’s and the influence of 𝑋(𝑡) on 𝑌(𝑡+1). Given that the prior implicit in the BGM model impedes the proper piecewise linear model with equal-length segments, the BGM model tends to infer the second-best explanation of the data: explaining the realisation of the 𝑌(𝑡)’s via a direct modelling of the autocorrelation between the 𝑌_𝑡’s themselves. This corresponds to (spurious) self-feedback loops.

We note that the novel BGM_𝐷 model has been particularly designed for dynamic data with a temporal structure. The BGM model is not restricted to such data, and can equally be applied to both static (steady-state) and dynamic (time series) data. However, the greater flexibility of the BGM model and the intrinsic implications for the effective prior probability on segment lengths and numbers renders its application to time series data suboptimal. This suggests that the proposed BGM_𝐷 offers a useful new tool for the analysis of dynamic processes.

7. Conclusions

Two paradigms for relaxing the nonhomogeneity/nonlinearity restriction of dynamic Bayesian networks have been proposed in the literature: the changepoint process and the free allocation mixture model. The latter provides the proper approximation of a nonlinear regulation process by a piecewise linear process. The former provides a similar approximation, but under the assumption that the temporal processeses are sufficiently smooth, as the assignment of observations to mixture components of the model is done in the temporal domain rather than the domain of regulatory variables. It is obvious that inference in the free allocation model has a considerably higher computational complexity than the changepoint process. However, we have additionally discussed several principled shortcomings that are intrinsic to the methodology per se. We have proposed a new variation of the BGM model [4], which has turned out to be more suitable for the reconstruction of regulatory networks from nonstationary gene expression time series. Like the BGM model the new BGM_𝐷 model is based on a mixture model dividing the data points into compartments. The network structures are kept fixed among time series segments, and each segment is modelled separately and independently with the Gaussian BGe model for Bayesian networks. The methodological difference is that the BGM_𝐷 model employs a changepoint process to divide time series into segments instead of a free unrestricted allocation of data points. The practical inference follows the Bayesian paradigm and samples the network, the number of changepoints and the changepoint locations from the posterior distribution with Markov chain Monte Carlo (MCMC).

In a first step, the inference problem was based on synthetic data from small network domains possessing self-feedback loops. Our empirical results show that the proposed BGM_𝐷 model suppresses spurious self-feedback loops and yields a higher network reconstruction accuracy than the standard BGe model or the BGM model. We also cross-compared the performance of the models on three real biological systems. On gene expression, time series related to (i) viral challenge of macrophages and (ii) circadian rhythms in Arabidopsis heat maps of the connectivity scores revealed that the new BGM_𝐷 model infers the biologically expected two-stage structures ((i) dichotomy between the healthy and diseased state of the cell and (ii) the diurnal contrast between light and darkness) more clearly than the BGM model. For assessing the statistical significance of the improvement we focused on predictive probabilities, and the proposed BGM_𝐷 model yields consistently higher scores than the BGM model. We extracted a gene regulatory network for the circadian clock-regulated genes in the Arabidopsis thaliana domain from the BGM_𝐷 inference results, and the reconstructed network shows features that are consistent with the biological literature.

Furthermore, for an independent comparison we combined the free allocation model and the changepoint process model with the discrete multinomial BDe scoring metric for Bayesian networks. Empirical results on binary gene expression time series related to muscle development in Drosophila melanogaster were consistent with the results from the first study on continuous data. The (discrete) changepoint process model (BDM_𝐷) infers a time series segmentation that is more consistent with the morphogenetical stages in Drosophila melanogaster than the free allocation model (BDM).

We note that the ideal approach—from a biological point of view—for these three applications would be to use a supervised approach, for example, as described in Werhli and Husmeier [27], and to exploit the biological knowledge we have about the experimental conditions (e.g., the morphogenetical stages of Drosophila) for data segmentation before inference. However, we elected to use these data as a test case for evaluating the efficiency of unsupervised learning for the proposed changepoint-process DBN model.

Finally, we cross-compared both models from a theoretical point of view. The theoretical study was based on a comparison of (i) the a priori imposed temporal connectivity structures and (ii) the prior probability ratios between the heterogeneous and the homogeneous states of the models. These results are consistent with our empirical findings, and lead to a deeper insight into the intrinsic difference between the two models.

Acknowledgments

M. Grzegorcyzk is supported by the Graduate School “Statistische Modellbildung” of the Department of Statistics, TU Dortmund University. D. Husmeier is supported by the Scottish Government Rural and Environment Research and Analysis Directorate (RERAD).