Sophie Lèbre
email: s.lebre@imperial.ac.uk

Inferring dynamic genetic networks with low order independencies

Keywords: conditional independence, Dynamic Bayesian Network, Directed Acyclic Graph, networks inference, time series modelling

In this paper, we propose a novel inference method for dynamic genetic networks which makes it possible to deal with a number of time measurements n much smaller than the number of genes p. The approach is based on the concept of low order conditional dependence graph which we extend here to the case of Dynamic Bayesian Networks. Most of our results are based on the theory of graphical models associated with Directed Acyclic Graphs (DAGs). In this way, we define a DAG G which describes exactly the full order conditional dependencies given the past of the process. Then, to cope with the large p and small n estimation case, we propose to approximate DAG G by considering low order conditional independencies. We introduce partial q th order conditional dependence DAGs and analyze their probabilistic properties. In general, DAGs G (q) differ from G but still reflect relevant dependence facts for sparse networks such as genetic networks. By using this approximation, we set out a non-Bayesian inference method and demonstrate the effectiveness of this approach on both simulated and real data analysis. The inference procedure is implemented in the R package 'G1DBN' which is available from the CRAN archive.

Introduction

The development of microarray technology allows to simultaneously measure the expression levels of many genes at a precise time point. Thus it has 1 become possible to observe gene expression levels across a whole process such as the cell cycle or response to radiation or different treatments. The objective is now to recover gene regulation phenomena from this data. We are looking for simple relationships such as "gene i activates gene j". But we also want to capture more complex scenarios such as auto-regulations, feed-forward loops, multi-component loops... as described by Lee et al. [START_REF] Lee | Transcriptional regulatory networks in Saccharomyces cerevisiae[END_REF] in the case of the transcriptional regulatory network of the yeast Saccharomyces cerevisiae.

To such an aim, we both need to accurately take into account temporal dependencies and to deal with the dimension of the problem when the number p of observed genes is much higher than the number n of observation time points. Moreover we know that most of the genes whose expression has been monitored using microarrays are not taking part in the temporal evolution of the system. So we want to determine the few 'active' genes that are involved in the regulatory machinery, as well as the relationships between them. In short, we want to infer a network representing the dependence relationships which govern a system composed of several agents from the observation of their activity across short time series.

Static Modelling Such gene networks were first described using static modelling and mainly non oriented networks. One of the first tools used to describe interactions between genes is the relevance network [START_REF] Butte | Discovering functional relationships between RNA expression and chemotherapeutic susceptibility using relevance networks[END_REF] or correlation network [START_REF] Steuer | Observing and interpreting correlations in metabolomic networks[END_REF]. Better known as the covariance graph [START_REF] Cox | Multivariate dependencies: Models, analysis and interpretation[END_REF] in graphical models theory, this undirected graph describes the pair-wise correlation between genes. Its topology is derived from the covariance matrix between the gene expression levels; an undirected edge is drawn between two variables whenever they are correlated. However, the correlation between two variables may be caused by linkage with other variables. This creates spurious edges due to indirect dependence relationships.

Consequently, there has been great interest in the concentration graph [START_REF] Lauritzen | Graphical models[END_REF], also called the covariance selection model, which describes the conditional dependence structure between gene expression using Graphical Gaussian Models (GGMs). Let Y = (Y i ) 1≤i≤p be a multivariate Gaussian vector representing the expression levels of p genes. An undirected edge is drawn between two variables Y i and Y j whenever they are conditionally dependent given the remaining variables (See Figure 1B). The standard theory of estimation in GGMs [START_REF] Lauritzen | Graphical models[END_REF][START_REF] Whittaker | Graphical models in applied multivariate statistics[END_REF] can be exploited only when the number of measurements n is much higher than the number of variables p. This ensures that the sample covariance matrix is positive definite with probability one. However, in most microarray gene expression datasets, we have to cope with the opposite situation (n << p). Thus, the growing interest in "small n, large p" furthered the development of numerous alternatives (Schäfer and Strimmer [31,[START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF] , Waddell and Kishino [START_REF] Waddell | Correspondence analysis of genes and tissue types and finding genetics links from microarray data[END_REF][START_REF] Waddell | Cluster inference methods and graphical models evaluated on nci60 microarray gene expression data[END_REF], The concentration graph corresponding to the motif A. For all i ≥ 3, Y i is a Gaussian variable representing the expression level of gene G i . Some cycles cannot be represented on the concentration graph. (C) Dynamic network equivalent to the regulation motif A. Each vertex X i t represents the expression level of gene G i at time t. This graph is acyclic and allows to define a Bayesian network.

Toh and Horimoto [START_REF] Toh | Inference of a genetic network by a combined approach of cluster analysis and graphical gaussian modeling[END_REF][START_REF] Toh | System for automatically inferring a genetic network from expression profiles[END_REF], Wu et al. [START_REF] Wu | Interactive analysis of gene interactions using graphical gaussian model[END_REF], Wang et al. [START_REF] Wang | Mgraph: graphical models for microarray data analysis[END_REF]). Even though concentration graphs allow to point out some dependence relationships between genes, they do not offer an accurate description of the interactions. Firstly, no direction is given to the interactions. Secondly, some motifs containing cycles as in Figure 1A cannot be properly represented.

Contrary to the previous undirected graphs, Bayesian networks (BNs) [START_REF] Friedman | Using Bayesian networks to analyse expression data[END_REF] model directed relationships. Based on a probabilistic measure, a BN representation of a model is defined by a Directed Acyclic Graph (DAG) and the set of conditional probability distributions of each variable given its parents in the DAG [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF]. The theory of graphical models [START_REF] Whittaker | Graphical models in applied multivariate statistics[END_REF][START_REF] Edwards | Introduction to Graphical Modelling[END_REF][START_REF] Lauritzen | Graphical models[END_REF] then allows to derive conditional independencies from this DAG. However, the acyclicity constraint in static BNs is a serious restriction given the expected structure of genetic networks.

Dynamic Bayesian networks This limitation can be overcome by employing Dynamic Bayesian networks (DBNs) introduced for the analysis of gene expression time series by Friedman et al. [START_REF] Friedman | Learning the structure of dynamic probabilistic networks[END_REF] and Murphy and Mian [START_REF] Murphy | Modelling gene expression data using dynamic Bayesian networks[END_REF]. In DBNs, a gene is no longer represented by a single vertex but by as many vertices as time points in the experiment. A dynamic network (Figure 1C) can then be obtained by unfolding in time the initial cyclic motif in Figure 1A. The direction according to time guarantees the acyclicity of this dynamic network and consequently allows to define a Bayesian network. The nature of the relationships (positive/negative) does not appear in this DAG but is derived from estimates of the model parameters.

The very high number p of genes simultaneously observed raises a dimension problem. Moreover, a large majority of time series gene expression data contain no or very few repeated measurements of the expression level of the same gene at a given time. Hence, we assume that the process is homogeneous across time. This means that the system is considered to be governed by the same rules during the whole experiment. Consequently, the temporal dependencies are homogeneous: any edge is present or absent during the whole process. This is a strong assumption which is not necessarily satisfied. Nevertheless, this condition is necessary to carry out estimation unless we have several measurements of each gene expression at each time point.

Up to now, various DBN representations based on different probabilistic models have been proposed (discrete models [START_REF] Ong | Modelling regulatory pathways in e. coli from time series expression profiles[END_REF][START_REF] Zou | A new dynamic Bayesian network (DBN) approach for identifying gene regulatory networks from time course microarray data[END_REF], multivariate autoregressive process [START_REF] Opgen-Rhein | Learning causal networks from systems biology time course data: an effective model selection procedure for the vector autoregressive process[END_REF], State Space or Hidden Markov Models [START_REF] Perrin | Gene networks inference using dynamic Bayesian networks[END_REF][START_REF] Wu | Modeling gene expression from microarray expression data with state-space equations[END_REF][START_REF] Rangel | Modeling t-cell activation using gene expression profiling and state-space models[END_REF][START_REF] Beal | A Bayesian approach to reconstructing genetic regulatory networks with hidden factors[END_REF], nonparametric additive regression model [START_REF] Imoto | Estimation of genetic networks and functional structures between genes by using Bayesian networks and nonparametric regression[END_REF][START_REF] Imoto | Bayesian network and nonparametric heteroscedastic regression for nonlinear modeling of genetic network[END_REF][START_REF] Kim | Dynamic Bayesian network and nonparametric regression for nonlinear modeling of gene networks from time series gene expression data[END_REF][START_REF] Sugimoto | Inference of gene regulatory networks by means of dynamic differential Bayesian networks and nonparametric regression[END_REF]). See also Kim et al. [START_REF] Kim | Inferring gene networks from time series microarray data using dynamic Bayesian networks[END_REF] for a review of such models. Faced with so much diversity, we introduce in this paper sufficient conditions for a model to admit a DBN representation and we set out a concrete interpretation in terms of dependencies between variables by using the theory of graphical models for DAGs.

Our DBN representation is based on a DAG G (e.g. like the DAG of Fig. 1C) which describes exactly the full order conditional dependencies given all the remaining past variables (See Section 1). This approach extends the principle of the concentration graph showing conditional independencies to the dynamic case.

Dimension reduction Even under the assumption of homogeneity, which enables to use the pairs of successive time point gene expression as repeated measurements, we have to deal with the "curse of dimensionality" when inferring the structure of DAG G. The difficulty lies in coping with the large p and small n estimation case. Several inference methods have been proposed for the estimation of the topology of the DAG defining the various DBNs quoted above. To name a few, Murphy [START_REF] Murphy | The bayes net toolbox for matlab[END_REF] implemented several Bayesian structure learning procedures for dynamic models in the Matlab package BNT (Bayes Net Toolbox); Ong et al. [START_REF] Ong | Modelling regulatory pathways in e. coli from time series expression profiles[END_REF] reduce the dimension of the problem by considering prior knowledge; Perrin et al. [START_REF] Perrin | Gene networks inference using dynamic Bayesian networks[END_REF] use an extension of the linear regression; Wu et al. [START_REF] Wu | Modeling gene expression from microarray expression data with state-space equations[END_REF] use factor analysis and Beal et al. [START_REF] Beal | A Bayesian approach to reconstructing genetic regulatory networks with hidden factors[END_REF] develop a variational Bayesian method; Zou and Conzen [START_REF] Zou | A new dynamic Bayesian network (DBN) approach for identifying gene regulatory networks from time course microarray data[END_REF] limit potential regulators to the genes with either earlier or simultaneous expression changes and estimate the transcription time lag; Opgen-Rhein and Strimmer [START_REF] Opgen-Rhein | Learning causal networks from systems biology time course data: an effective model selection procedure for the vector autoregressive process[END_REF] proposed a model selection procedure based on an analytic shrinkage approach. However, a powerful approach based on the consideration of zero-and first-order conditional independencies to model concentration graphs has gained attention. When n << p, Wille et al. [START_REF] Wille | Sparse graphical gaussian modeling for genetic regulatory network inference[END_REF][START_REF] Wille | Low-order conditional independence graphs for inferring genetic networks[END_REF] propose to approximate the concentration graph by the graph G 0-1 describing zero-and first-order conditional independence. An edge between the variables Y i and Y j is drawn in the graph G 0-1 if and only if, zero-and first-order correlations between these two variables both differ from zero, that is, if

r(Y i , Y j ) = 0 and ∀k ∈ {1, ..., p}\{i, j}, r(Y i , Y j |Y k ) = 0, (1)
where

r(Y i , Y j |Y k ) is the partial correlation between Y i and Y j given Y k .
Hence, whenever the correlation between two variables Y i and Y j can be entirely explained by the effect of some variable Y k , no edge is drawn between them. This procedure allows a drastic dimension reduction: by using first order conditional correlations, estimation can be carried out accurately even with a small number of observations. Even if the graph of zero-and first-order conditional independence differs from the concentration graph in general, it still reflects some measure of conditional independence. Wille et al. show through simulations that the graph G 0-1 offers a good approximation of sparse concentration graphs and demonstrate that both graphs coincide exactly if the concentration graph is a forest ([47], Corollary 1). This approach has also been used by Magwene and Kim [START_REF] Magwene | Estimating genomic coexpression networks using first-order conditional independence[END_REF] and de la Fuente et al. [START_REF] De La Fuente | Discovery of meaningful associations in genomic data using partial correlation coefficients[END_REF] for estimating undirected gene networks from microarray gene expression of the yeast Saccharomyces cerevisiae. Castelo and Roverato [START_REF] Castelo | Graphical model search procedure in the large p and small n paradigm with applications to microarray data[END_REF] investigate such undirected q th order partial independence graphs for q ≥ 1 and present a thorough analysis of their properties. In this paper, we extend this approach by defining q th order order conditional dependence DAGs G (q) for DBN representations. Then, by basing our results on these low order conditional dependence DAGs, we propose a novel inference method for dynamic genetic networks which makes it possible to deal with the "small n, large p" problem.

The remainder of the paper is organized as follows. In Section 1, we provide sufficient conditions for a DBN modelling of time series describing temporal dependencies. In particular, we show the existence of a minimal DAG G which allows such a DBN representation. To reduce the dimension of the estimation of the topology of G, we propose to approximate G by q th order conditional dependence DAGs G (q) and analyze their probabilistic properties in Section 2. From conditions on the topology of G and the faithfulness assumption, we establish inclusion relationships between both DAGs G and G (q) . In Section 3, we exploit our results on DAGs G (q) Finally, validation is obtained on both simulated and real data in Section 4. We use our inference procedure for the analysis of two microarray time course data sets: the Spellman's yeast cell cycle data [START_REF] Spellman | Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization[END_REF] and the diurnal cycle data on the starch metabolism of Arabidopsis Thaliana collected a DAG whose vertices are defined by X and edges by E(G) ⊆ X × X, G the "true" DAG describing the set of full order conditional dependencies, G (q) q th order conditional dependence DAG, by Smith et al. [START_REF] Smith | Diurnal Changes in the Transcriptome Encoding Enzymes of Starch Metabolism Provide Evidence for Both Transcriptional and Posttranscriptional Regulation of Starch Metabolism in Arabidopsis Leaves[END_REF].

A minimal DBN representation

Let P = {1 ≤ i ≤ p} describe the set of observed genes and N = { 1 ≤ t ≤ n } the set of observation times. In this paper, we consider a discrete-time stochastic process X = {X i t ; i ∈ P, t ∈ N} taking real values and assume the joint probability distribution P of the process X has density f with respect to Lebesgue measure on R p×n . We denote by X t = {X i t ; i ∈ P } the set of the p random variables observed at time t and X 1:t = {X i s ; i ∈ P, s ≤ t} the set of the random variables observed before time t.

The main result of this section is set out in Proposition 3; we show that process X admits a DBN representation according to a minimal DAG G whose edges describe exactly the set of direct dependencies between successive variables X j t-1 , X i t given the past of the process. For an illustration, the minimal DAG G is given in the case of an AR(1) model in Subsection 1.2. Most of our results are derived from the theory of graphical models associated with DAGs [START_REF] Lauritzen | Graphical models[END_REF]. Note that, even though we need to consider a homogeneous DBN for the inference of gene interaction networks, the theoretical results introduced in Sections 1 and 2 are valid without assuming homogeneity across time.

Background

Theory of graphical models associated with DAGs Let G = (X, E(G)) be a DAG whose vertices are the variables X = {X i t ; i ∈ P, t ∈ N} and whose set of edges E(G) is a subset of X × X. We quickly recall here elements of the theory of graphical models associated with DAGs [START_REF] Lauritzen | Graphical models[END_REF]. A characterization of a Bayesian Network (BN) representation for a process X is given in Proposition 1.

Definition 1 (Parents, Lauritzen [START_REF] Lauritzen | Graphical models[END_REF]) The parents of a vertex X i t in G, denoted by pa(X i t , G), are the variables having an edge pointing towards the vertex X i t in G, pa(X i t , G) := {X j s such that (X j s , X i t ) ∈ E(G); j ∈ P, s ∈ N}.

Proposition 1 (BN representation, Pearl [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF]) The probability distribution P of process X admits a Bayesian Network (BN) representation according to DAG G whenever its density f factorizes as a product of the conditional density of each variable

X t i given its parents in G, f (X) = i∈P t∈N f (X i t |pa(X i t , G)).
Throughout this paper, a central notion is that of conditional independence of random variables. Two random variables U and V are conditionally independent given a third variable W (and we write U ⊥ ⊥ V | W ) if they are independent in the joint probability distribution P U,V,W of the three random variables (U, V, W ). In other words, U and V are conditionally independent given W if for any possible value w of W , variables U and V are independent given the variable W = w. This result generalizes to disjoint sets of variables. Such conditional independence relationships can be obtained from a BN representation by using graphical theory associated with DAGs, which is essentially based on the directed global Markov property recalled in Proposition 2.

Definition 2 (Moral graph, Lauritzen [START_REF] Lauritzen | Graphical models[END_REF]) The moral graph G m of DAG G is obtained from G by first 'marrying' the parents (draw an undirected edge between each pair of parents of each variable X i t ) and then deleting the directions of the original edges of G. For an illustration, Figure 2A displays the moral graph of the DAG in Figure 1C. Definition 3 (Ancestral set, Lauritzen [START_REF] Lauritzen | Graphical models[END_REF]) The subset S is ancestral if and only if, for all α ∈ S, the parents of α satisfy pa(α, G) ⊆ S. Hence, for any subset S of vertices, there is a smallest ancestral set containing S which is denoted by An(S). Then G An(S) refers to the graph of the smallest ancestral set An(S). See Figure 2B for an illustration.

Proposition 2 (Directed global Markov property, Lauritzen [START_REF] Lauritzen | Graphical models[END_REF], Corollary 3.23) Let P admit a BN representation according to G. Then,

E ⊥ ⊥ F | S,
whenever all paths from E to F intersect S in (G An(E∪F ∪S) ) m , the moral graph of the smallest ancestral set containing E ∪ F ∪ S. We say that S separates E from F . . For all t > 1, the parents of the variable X 1 t are 'married', that is connected by an undirected edge. (B) Moral graph of the smallest ancestral set containing the variables X 1 t+1 , its parents in the DAG in Figure 1C and

3 t-1 t-1 t-1 t t t t+1 t+1 t+1 X X X X X X X X X . . . . . . . . . X X X X X X 1 2 t-1 t-1 t t t 1 2 3 t+1 1 . . . . . .
X 3 t . As the set (X 1 t , X 2 t ) blocks all paths between X 3 t and X 1 t+1 , thus {X 1 t , X 2 t } separates X 1 t+1 from X 3 t and we have X 1 t+1 ⊥ ⊥ X 3 t | (X 1 t , X 2 t ).

Sufficient conditions for DBNs representation

We recall here sufficient conditions under which the probability distribution P of process X admits a BN representation according to a dynamic network (e.g. in Figure 1C). We first assume that the observed process X t is first-order Markovian (Assumption 1). That is, the expression level of a gene at a given time t only depends on the past through the gene expression levels observed at the previous time t -1. Then we assume that the variables observed simultaneously are conditionally independent given the past of the process (Assumption 2). In other words, we consider that time measurements are close enough so that gene expression level X i t measured at time t is better explained by the previous time expression levels X t-1 than by some current expression level X j t .

Assumption 1 The stochastic process X t is first-order Markovian, that is,

∀t ≥ 3, X t ⊥ ⊥ X 1:t-2 | X t-1 .
Assumption 2 For all t ≥ 1, the random variables {X i t } i∈P are conditionally independent given the past of the process

X 1:t-1 , that is, ∀t ≥ 1, ∀i = j, X i t ⊥ ⊥ X j t | X 1:t-1 .
Assumptions 1 and 2 allow the existence of a DBN representation of the distribution P according to DAG G full = (X, {(X j t-1 , X i t )} i,j∈P,t>1 ) which contains all the edges pointing out from a variable observed at some time t-1 towards a variable observed at the next time t (See Lemma 1 in Appendix A.1). The direction of the edges according to time guarantees the acyclicity of G full .

Minimal DAG G

Existence and definition Among the DAGs included in G full , we show that the probability distribution P factorizes according to a minimal DAG, which we denote by G (See Lemma 2, Appendix A.1). The set of edges of G is exactly the set of full order conditional dependencies between successive variables given the past of the process as set up in the Proposition 3 (See Proof in Appendix A.2).

Proposition 3 (Existence of minimal DAG G, the smallest subgraph of G full allowing DBN modelling) Let P j = P \{j} and X P j t = {X k t ; k ∈ P j } refer to the set P j of p -1 variables observed at time t. Whenever Assumptions 1 and 2 are satisfied, the probability distribution P admits a DBN representation according to DAG G whose edges describe exactly the full order conditional dependencies between successive variables X j t-1 and X i t given the remaining variables X

P j t-1 observed at time t -1, G = X, (X j t-1 , X i t ); X i t ⊥ ⊥ X j t-1 |X P j t-1 i,j∈P,t∈N
Moreover, DAG G is the smallest subgraph of G full according to which P admits a DBN representation.

Thus in DAG G, the set of parents pa(

X i t , G) of a variable X i t is the smallest subset of X t-1 such that the conditional densities satisfy f (X i t |pa(X i t , G)) = f (X i t |X t-1
). The set of parents of a variable can be seen as the only variables on which this variable depends directly. So G is the DAG we want to infer in order to recover potential regulation relationships from gene expression time series. From Proposition 3, any pair of successive variables (X j t-1 , X i t ) which are non adjacent in G are conditionally independent given the parents of X i t . In short, for all i, j in P , for all t > 1, we have,

(X j t-1 , X i t ) / ∈ E( G) ⇔ X i t ⊥ ⊥ X j t-1 | pa(X i t , G
). We will make use of this result in Section 2 in order to define low order conditional dependence DAGs for the inference of G.

Minimal DAG G for an AR(1) process Consider the following first order auto-regressive model (AR(1)) with a diagonal error covariance matrix Σ,

X 1 ∼ N (µ 1 , Σ 1 ) (2) ∀t > 1, X t = AX t-1 + B + ε t , ε t ∼ N (0, Σ), (3) ∀s, t ∈ N, Cov(ε t , ε s ) = δ ts Σ, ( 4
) ∀s > t, Cov(X t , ε s ) = 0. ( 5
)
where A=(a ij ) 1≤i≤p,1≤j≤p is a real matrix of size p×p, B=(b i ) 1≤i≤p is a real column vector, Σ = Diag (σ 2 ii ) 1≤i≤p is the diagonal error covariance matrix of size p × p and for all s,t in N, δ ts =1l {s=t} . Equation [START_REF] Butte | Discovering functional relationships between RNA expression and chemotherapeutic susceptibility using relevance networks[END_REF] implies that the coefficient matrices are uniquely determined from the covariance function of X t . This modelling assumes homogeneity across time (constant matrix A) and linearity of the dependency relationships. From (3) and ( 5), the model is first order Markovian (Assumption 1). From (4), Assumption 2 is satisfied whenever the error covariance matrix Σ is diagonal. Thus from Proposition 3, the probability distribution of the AR(1) process defined by equations (2-5) factorizes according to the minimal DAG GAR(1) whose edges correspond to the non-zero coefficients of matrix A. Indeed, if matrix Σ is diagonal, each element a ij is the regression coefficient of the variable X i t on X j t-1 given X

P j t-1 , that is a ij = Cov(X i t , X j t-1 | X P j t-1 )/V ar(X j t-1 | X P j t-1 ).
As process X is Gaussian, the set of null coefficients of matrix A exactly describes the conditional independencies between successive variables, thus if Σ is diagonal, we have,

a ij = 0 ⇔ ∀t > 1, X i t ⊥ ⊥ X j t-1 |X P j t-1 .
Finally, DAG GAR(1) has an edge between two successive variables X j t-1 and X i t , for all t > 1, whenever the coefficient a ij of the matrix A differs from zero, GAR(1) := X, (X j t-1 , X i t ) such that a ij = 0; t > 1, i, j ∈ P .

As an illustration, any AR(1) process whose matrix Σ is diagonal and matrix A has the following form,

A =    a 11 a 12 0 a 21 0 0 0 a 32 0    ,
admits a BN representation according to the dynamic network of Fig. 1C (p =3).

2 Introducing q th order dependence DAGs G (q) for DBNs

In this paper, we propose to use the DBN modelling according to DAG G (introduced in Proposition 3) to model genetic regulatory networks from gene expression time series. Reverse discovery of DAG G requires to determine, for each variable X i t , the set of variables X j t-1 observed at time t -1 on which variable X i t is conditionally dependent given the remaining variables X P j t-1 . However, even under the time homogeneity assumption discussed in the introduction, standard estimation methods do not allow us to infer the parameters of a regression model for p genes (i.e. p 2 possible edges) from np measurements. We still have to face the 'curse of dimensionality' since the number of genes p, is much higher than the number of measurements n.

In order to reduce the dimension, we approximate DAG G by q th order conditional dependence DAGs G (q) (q < p). To such an end, we extend to DBNs the approach based on the consideration of low order independencies introduced by Wille et al. [START_REF] Wille | Sparse graphical gaussian modeling for genetic regulatory network inference[END_REF][START_REF] Wille | Low-order conditional independence graphs for inferring genetic networks[END_REF] for GGM approximation (See more details on low order independence graphs for GGMs in Section ). After defining q th order conditional dependence DAGs G (q) for DBNs, we investigate the manner in which they allow us to approximate the DAG G describing full order conditional dependencies.

DAG G (q) definition

Let q be smaller than p. In the q th order dependence DAG G (q) , whenever there exists a subset X Q t-1 of q variables among the set of p -1 variables X P j t-1 such that X j t-1 and X i t are conditionally independent given X Q t-1 , no edge is drawn between the two successive variables X j t-1 and X i t . In short, DAGs G (q) are defined as follows, Definition 4 q th -order conditional dependence DAG G (q) ∀q < p, G (q) = X, (X j t-1 , X i t );

∀Q ⊆ P j , |Q| = q, X i t ⊥ ⊥ X j t-1 |X Q t-1 i,j∈P,t∈N
.

DAGs G (q) offer a way of producing dependence relationships between the variables, but they are no longer associated with a BN representation which would call for more global relationships. Note that the definition of q th order partial dependence DAG G (q) is based on exact q th order independencies (not on all partial independencies lower than q as in the partial order correlation network used by Wille and Bühlmann [START_REF] Wille | Low-order conditional independence graphs for inferring genetic networks[END_REF]). Indeed, we consider that including only the q th order dependencies better reflects the true DAG G. In particular, for p variables, DAG G (p-1) is DAG G. This definition is possible for DBNs because dynamic modelling essentially differs from static correlation network modelling 1 .

In general, DAGs G (q) differ from DAG G. For instance, the approximation of the DAG of Figure 1C by the 1 st order conditional dependence 1 In particular, contrary to the case of correlation network, the " V " structures (or structures with multiple parents) do not generate spurious edges in the case of DBN since the definition of the DAG G defining full order dependencies does not allow edges between variables observed at the same time. Thus, for instance, when considering the following " V " structure X j t-1 → X i t ← X k t-1 , no spurious edge can be inferred between the variables X j t-1 and X k t-1 . (1) (obtained from the DAG in Figure 1C). The spurious dashed arrow may appear in G (1) .

DAG may give rise to the spurious edge X 3 t → X 1 t+1 , for all t < n (See Figure 3). Indeed, X 1 t (resp. X 2 t ) does not separate X 1 t+1 from X 3 t in the smallest moral graph containing the variables 2B. Nevertheless, if the vertices of G have few parents, DAGs G (q) bring relevant information about the topology of G, even for small values of q. In the following, we give characterizations of low order conditional dependence DAGs G (q) and analyze the accuracy of the approximations they offer.

X 1 t+1 ∪ X 3 t ∪ X 1 t (resp. X 1 t+1 ∪ X 3 t ∪ X 2 t ) displayed in Figure

A restricted number of parents

In some known gene regulation mechanisms, it is the case that a few genes regulate many other genes (e.g. the single input modules in the transcriptional regulatory network of S. Cerevisiae [START_REF] Lee | Transcriptional regulatory networks in Saccharomyces cerevisiae[END_REF]). However, we do not expect a single gene to be regulated by many genes at the same time. So the number of parents in gene interaction networks is expected to be relatively small. In this section, we analyze the properties of G (q) when the number of parents in G is lower than q.

Let us denote by N pa (X i t , G) the number of parents of X i t in DAG G and N Max pa ( G) the maximal number of parents of any variable

X i t in G, N pa (X i t , G) = pa(X i t , G) , N Max pa ( G) = Max i∈P,t∈N N pa (X i t , G) .
The next results hold when the number of parents in G is restricted.

Proposition 4 If N pa (X i t , G) ≤ q then we have, (X j t-1 , X i t ) / ∈ E( G) ⇒ (X j t-1 , X i t ) / ∈ E(G q ) .
Corollary 1 For all q ≥ N Max pa ( G), we have G ⊇ G (q) . Proposition 5 Let X be a Gaussian process. If N Max pa ( G) ≤ 1 then G = G (1) .

Consider a variable X i t having at most q parents in G (q < p). Let X j t-1 be a variable observed at the previous time t -1 and having no edge pointing towards X i t in G. In the moral graph of the smallest ancestral set containing X i t ∪ X j t-1 ∪ pa(X i t , G), the set of parents pa(X i t , G) separates X i t from X j t-1 . From Proposition 2, we have X i t ⊥ ⊥ X j t-1 | pa(X i t , G). The number of parents pa(X i t , G) is smaller than q, so the edge X j t-1 → X i t is not in G (q) . This establishes Proposition 4. Consequently, if the maximal number of parents in G is lower than q, then G (q) is included in G (Corollary 1). In this case, G (q) does not contain spurious edges.

The converse inclusion relationship is not true in general 2 . Nevertheless, if each variable has at most one parent, the converse inclusion G ⊆ G (1) is true if the process is Gaussian and q = 1 (Proposition 5, see proof in Appendix A.2). At a higher order, we need to assume that all conditional independencies can be derived from G, that is P is faithful to G.

Faithfulness

Definition 5 (faithfulness, Spirtes [START_REF] Spirtes | Causation, prediction and search[END_REF]) A distribution P is faithful to a DAG G if all and only the independence relationships true in P are entailed by G (as set up in Proposition 2).

Theorem 1 (Measure zero for unfaithful Gaussian (Spirtes [START_REF] Spirtes | Causation, prediction and search[END_REF]) and discrete (Meek [23]) distributions) Let π N G (resp. π D G ) be the set of linearly independent parameters needed to parameterize a multivariate normal distribution (resp. discrete distribution) P which admits a factorization according to a DAG G. The set of distributions which are unfaithful to G has measure zero with respect to Lebesgue measure over π N G (resp. over π D G ).

From Definition 5, whenever the distribution P is faithful to G, any subset X Q t-1 ⊆ X t-1 , with respect to which X i t and X j t-1 are conditionally independent, separates X i t and X j t-1 in the moral graph of the smallest ancestral set containing

X i t ∪ X j t-1 ∪ X Q t-1 .
Under this assumption, we can derive interesting properties on G from the topology of low order dependence DAGs G (q) . As there is no way to assess a probability distribution to be faithful to a DAG, this assumption has often been criticized. However, Theorem 1, established by Spirtes [START_REF] Spirtes | Causation, prediction and search[END_REF] for the Gaussian distribution and extended to discrete distributions by Meek [START_REF] Meek | Strong completeness and faithfulness in Bayesian networks[END_REF], makes this assumption 2 As an illustration, let X j t-1 → X i t be an edge of G then in essence (See Prop 3) X i t and X j t-1 are conditionally dependent given the remaining variables X Pj t-1 . There may however exist a subset of q variables X Q t-1 , where Q is a subset of P \{j} of size q, such that X i t and X j t-1 are conditionally independent with respect to this subset X Q t-1 . Indeed, even though the topology of G allows us to establish some conditional independencies, DAG G does not necessarily allow to derive all of them. Two variables can be conditionally independent given a subset of variables whereas this subset does not separate these two variables in G.

reasonable at least in a measure-theoretic sense. Moreover this assumption remains very reasonable in a modelling framework where the network to be inferred describes actual interaction relationships. The next propositions are derived from the faithfulness of the distribution P to G (See proofs in Appendix A.2).

Proposition 6 Assume P is faithful to G. For all q < p, we have G ⊆ G (q) . Corollary 2 Assume P is faithful to G. For all q ≥ N Max pa ( G), we have

G = G (q) . Proposition 7 Assume P is faithful to G. If N pa (X i t , G (q) ) ≤ q then (X j t-1 , X i t ) ∈ E(G (q) ) ⇒ (X j t-1 , X i t ) ∈ E( G).
Corollary 3 Assume P is faithful to G. For all q ≥ N M ax pa (G (q) ), G = G (q) . Whenever P is faithful to G, DAG G (q) contains DAG G (Proposition 6). Even though we expect the number of parents in a gene interaction networks to be bounded aboce, the exact maximal number of parents N Max pa ( G) remains mostly unknown. However, we show that the edges of DAG G (q) pointing towards a variable having less than q parents in G (q) are edges of G too (Proposition 7). Thus, if P is faithful to G, knowledge of the topology of DAG G (q) only allows us to ascertain some edges of DAG G. From Propositions 6 and 7, we establish that both DAG G (q) and DAG G exactly coincide if any node of G (q) has less than q parents (Corollary 3).

G1DBN, a procedure for DBN inference

We introduced and characterized the q th order dependence DAGs G (q) , for all q < p, for dynamic modelling. We now exploit our results to develop a non-Bayesian inference method for DAG G defining a DBN representation for process X. Let q max be the maximal number of parents in G. From Corollary 3, inferring G amounts to inferring G (qmax) . However, the inference of G (qmax) requires to check, for each pair (i, j), if there exists a subset Q ⊆ P j of dimension q max such that X i t ⊥ ⊥ X j t-1 |X Q t-1 for all t > 1. So, for each pair (i, j), there are qmax p-1 potential sets that can lead to conditional independence. To test each conditional independence given any possible subset of q max variables is questionable both in terms of complexity and multiple testings.

To circumvent these issues, we propose to exploit the fact that the true DAG G is a subgraph of G (1) (Proposition 6) in order to develop an inference procedure for G. Indeed, the inference of G (1) is both faster (complexity) and more accurate (number of tests). Thus we introduce a 2 step-procedure for DBN inference. In the first step, we infer the 1 st order dependence DAG G (1) , then we infer DAG G from the estimated DAG Ĝ(1) . This 2 step-procedure, summarized in Figure 4, is implemented in a R package 'G1DBN' [START_REF] G1dbn | A package performing Dynamic Bayesian Network inference[END_REF] freely available from the Comprehensive R Archive Network.

3.1

Step 1: inferring G (1) We evaluate the likelihood of an edge (X j t-1 , X i t ) by measuring the conditional dependence between the variables X j t-1 and X i t given any variable X k t-1 . Assuming linear dependencies, we consider the partial regression coefficient a ij|k defined as follows,

X i t = m ijk + a ij|k X j t-1 + a ik|j X k t-1 + η i,j,k t ,
where the rank of the matrix (X j t-1 , X k t-1 ) t≥2 equals 2 and the errors {η i,j,k t } t≥2 are centered, have same variance and are not correlated.

We measure the conditional dependence between the variables X j t-1 and X i t given any variable X k t-1 by testing the null assumption H i,j,k 0 : "a ij|k = 0". To such an aim, we use one out of three M-estimators for this coefficient: either the familiar Least Square (LS) estimator, the Huber estimator, or the Tukey bisquare (or biweight) estimator. The two latter are robust estimators [START_REF] Fox | An R and S-Plus companion to applied regression[END_REF]. Then for each k = j, we compute the estimates âij|k according to one of these three estimators and derive the p-value p ij,k from the standard significance test:

under (H i,j,k 0 ) : " a ij|k = 0 ", âij|k σ(â ij|k ) ∼ t(n -4), (7)
where t(n-4) refers to a student probability distribution with n-4 degrees of freedom and σ(â ij|k ) is the variance estimates for âij|k . Thus, we assign a score S 1 (i, j) to each potential edge (X j t-1 , X i t ) equal to the maximum Max k =j (p ij|k ) of the p -1 computed p-values, that is the most favorable result to 1 st order conditional independence. This procedure does not derive p-values for the edges but allows to order the possible edges of DAG G (1) according to how likely they are. The smallest scores point out the most significant edges for G (1) . The inferred DAG Ĝ(1) contains the edges assigned a score below a chosen threshold α 1 .

3.2

Step 2: inferring G from G (1) We use the inferred DAG Ĝ(1) as a reduction of the search space. Indeed, from faithfulness, we know that G ⊆ G (1) (Proposition 6). Moreover, when DAG G is sparse, there are far fewer edges in G (1) than in the complete DAG G full defined in Section 1.1. Consequently, the number of parents of each variable in Ĝ( 1) is much smaller than n. Then model selection can be Choose either LS, Huber or Tukey estimator and set α 1 and α 2 thresholds.

Step 1: inferring G (1) . For all i ∈ P , For all j ∈ P , for all k = j, compute the p-value p ij|k from (7), S 1 (i, j) = M ax k =j (p ij|k ). E( Ĝ(1) ) = {(X j t-1 , X i t ) t>1 ; i, j ∈ P, such that S 1 (i, j) < α 1 }.

Step 2: inferring G from Ĝ(1) . If N M ax pa ( Ĝ( 1) ) ∼ n -1, choose a higher threshold α 1 and go to Step1. For all i such that N pa (X i t , Ĝ(1) ) ≥ 1, compute the p-value p

(2) ij from [START_REF] Edwards | Introduction to Graphical Modelling[END_REF]. For each pair (i, j) such that the set of edges (X j t-1 , X i t ) t>1 is in Ĝ(1) , we denote by a

S 2 (i, j) = p (2) ij f or all i, j ∈ P such that (X j t-1 , X i t ) t>1 ∈ Ĝ(1) , 1 otherwise. E( G) = {(X j t-1 , X i t ) t>1 ; i ∈ P, (i, j) ∈ P such that S 2 (i, j) < α 2 }.
(2) ij the regression coefficient,

X i t = m i + j∈pa(X i t , Ĝ(1) ) a (2) ij X j t-1 + η i t , (8)
where the rank of the matrix (X j t-1 ) t≥2,j∈pa(X i t , Ĝ(1) ) is |pa(X i t , Ĝ(1) )| and the errors {η i t } t≥2 are centered, have the same variance, and are not correlated. We assign to each edge of Ĝ(1) a score S 2 (i, j) equal to the p-value p

(2) ij derived from the significance test, under (H i,j 0 ) : " a

(2) ij = 0 ", â(2) ij σ(â (2) ij ) ∼ t(n -1 -|pa(X i t , Ĝ(1) )|). ( 9
)
The score S 2 (i, j) = 1 is assigned to the edges that are not in Ĝ(1) . The smallest scores indicate the most significant edges. The inferred DAG for G contains those edges whose score is below a chosen threshold α 2 . When G is sparse, Step 1 of G1DBN inference procedure gives already a good estimation of G (See Precision-Recall curves obtained for simulated data in Figure 5). Even better results can be obtained with the 2 stepprocedure which requires to tune two parameters α 1 and α 2 . Parameter α 1 is the selection threshold of the edges of Ĝ(1) in Step 1 (that is the dimension reduction threshold), whereas parameter α 2 is the selection threshold for the edges of G among the edges of DAG Ĝ(1) .

Choice of the thresholds

The choice of thresholds is often something non trivial, especially when using multiple testing. However, Step 1 of the procedure is conservative by construction. Indeed, the definition of score S 1 (equal to the maximum of p -1 p-values computed for testing 1st-order conditional independence) clearly supports the acceptation of the null assumption, i.e. the absence of an edge. Standard approaches for multiple testing correction do not apply to choose α 1 threshold. Thus we introduce a heuristic approach to choose α 1 threshold which is detailed in Supplementary Material [2], Section B. Overall, α 1 threshold is chosen so that, after the Step 1, the number of genes having exactly one parent in DAG G (1) predominates.

The choice of α 2 threshold is less problematic. Indeed, the second Step of the inference procedure is a standard multivariate regression. Then the usual thresholds 1%, 5% or 10 % can be chosen or even a lower threshold when a low number of edges is wanted. However, a large number of tests are computed (as many as edges in DAG G (1) ). In such multiple testing situations, a set of the predictions are expected to be false and it is useful to control this. We control the expected proportion of false positives edges, i.e. the False Discovery Rate (FDR) with the approach introduced by Benjamini and Hochberg 3 [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF].

Complexity of the algorithm

The complexity of this algorithm is O(p 3 ). However the scores (S 1 (i, j)) j∈P of the incoming edges of each target gene i can be computed separately by using parallel run. This option is available in the R package G1DBN by specifying the target gene i in the function DBNScoreStep1 dedicated to the Step 1 computation.

All the computations were performed on Redhat WS 4 AMD opteron 270 (2GHz). The computation time mostly depends on the number of TF genes, i.e. the genes allowed to be parents in the DAG to be inferred. For an illustration based on DBN inference performed from a real data set by Spellman [START_REF] Spellman | Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization[END_REF] containing 786 target genes in Section 4.3, the computation of Step 1 required 7 minutes when the set of possible TF genes was restricted to 18 genes (resp. 4 minutes with the lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and 7 seconds with the shrinkage procedure [START_REF] Opgen-Rhein | Learning causal networks from systems biology time course data: an effective model selection procedure for the vector autoregressive process[END_REF], which are two alternative approaches for DBN inference introduced in Section 4.1). When all the 786 genes can be TFs, 3 Let m be the number of remaining edges after Step 1, then Step 2 requires to compute m tests. Choose a maximal FDR level q and order the set of m observed p-values: p

(1) ≤ • • • ≤ p (i) ≤ • • • ≤ p (m)
. Then reject the null assumption (H (i) 0 : "Edge i is not DAG G") for all i ≤ k where k is defined as follows: k = max i : p (i) ≤ i m q . If no such i exists, reject no hypothesis. [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF] showed that this procedure ensures the FDR is lower than q m0 m ≤ q where m 0 is the number of true null hypotheses.

the computation was parallel run and required 19 minutes by target gene with G1DBN (resp. 8 minutes by target gene with the lasso and 5 minutes for the whole set of 786 target genes with the shrinkage procedure).

Step 2 of G1DBN is very quick and requires less than 5 seconds for the 786-TF study. Despite the need for more time, inference with G1DBN for a data set containing 800 genes is fully computable, especially when parallel running.

Validation

Comparison with two reference methods

We compare the G1DBN inference procedure with two reference methods for model selection for multivariate AR(1) process: the shrinkage approach by Opgen-Rhein and Strimmer [START_REF] Opgen-Rhein | Learning causal networks from systems biology time course data: an effective model selection procedure for the vector autoregressive process[END_REF] and the lasso (Least Absolute Shrinkage and Selection Operator) introduced by Tibshirani [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. Opgen-Rhein and Strimmer recently proposed a model selection based on an analytic approach using James-Stein-Type shrinkage. The procedure consists of first computing the partial correlation coefficients, r(X i t , X j t-1 |X P j t-1 ), from the shrinkage estimates of the partial regression coefficients, and second, selecting the edges with a local false discovery rate approach [START_REF] Efron | Local false discovery rates[END_REF]. Shrinkage inference is performed using the R code for shrinkage estimation4 by Opgen-Rhein and Strimmer.

The lasso (also called L1 shrinkage) combines shrinkage and model selection. The lasso estimates are obtained by minimizing the residual sum of squares subject to the sum of the absolute values of the coefficients being less than a constant. This approach offers the advantage that it automatically sets many regression coefficients to zero. We performed the lasso with the R package LARS developped by Efron et al. [START_REF] Efron | Least angle regression[END_REF].

Simulation study

As the discovery of genetic regulatory interaction is a field in progress, validation of predictions made on real gene expression data is only partial, which may render the estimation of true and false positive detection rate not fully reliable [START_REF] Husmeier | Sensitivity and specificity of inferring genetic regulatory interactions from microarray experiments with dynamic Bayesian networks[END_REF]. Thus we first investigate the accuracy of G1DBN, the shrinkage and the lasso inference procedures on simulated data.

Data generation

We generated 100 random time series according to a multivariate AR(1) model defined by parameters (A [p×p] ,B, Σ) for p = 50 genes. Since gene regulation networks are sparse, each matrix A contains 5% of non zero coefficients. While keeping the number of parents low, this does not prevent a vertex from having more than one parent. Non zero regression coefficients a ij , mean coefficients b i and error variances σ i were drawn from uniform distributions (a ij , b i ∼ U([-0.95; -0.05]∪ [0.05; 0.95]), i ∼ U[0.03, 0.08]). Time series were generated under the corresponding multivariate AR(1) models for n = 20 to 50.

Evaluation based on PR curves

We evaluated the performance of DBN inference procedures using the Precision-Recall (PR) curve as plotted in Figure 5. PR curves show the precision, equal to the Positive Predictive Value (PPV) on the ordinate against the recall, equal to the power, on the abscissa. PR curves are drawn by first ordering the edges by decreasing significance, and then by computing the PPV and power for the first selected edge and for each newly included edge successively. We recall the next definitions, Positive Predictive Value (PPV) = True Discovery Rate (TDR) = 1-False Discovery Rate (FDR) = T P T P + F P Recall = Sensitivity = Power = T P T P + F N where TP refers to the number of true positive edges, i.e. the number of edges which are selected by the inference procedure and actually belongs to the true DAG (used for simulating the data); FP refers to the number of false positive edges, i.e. the edges which are selected by the procedure but are not in the true DAG and FN refers to the number of false negative edges, i.e. the number of edges which are not selected by the procedure but are in the true DAG. Step 2 of the G1DBN approach drastically improves the results (threshold α 1 = 0.7). (B) Impact of noisy data, simulated using a non diagonal matrix Σ with either Gaussian or uniform noise, on the G1DBN procedure (Step 2) computed with LS estimates.

Simulation results

We show on Figure 5 the results obtained with n = 20, a length one can expect from existing gene expression time series. Figure 5A displays the average Precision Recall (PR) curves obtained with the various inference approaches when the error covariance matrix Σ is diagonal and the noise distribution is Gaussian. The Step 1 of the G1DBN procedure computed either with the LS estimator or with the Tukey estimator (dashed lines) gives a very high PPV for the very first selected edges. The Step 2 of the G1DBN procedure (solid line) drastically improves the results. It allows to maintain the PPV greater than 95 % while the power goes up to 50%. PR curves computed with the Huber estimates (not shown) led to comparable results. The lasso (dotted line) is clearly outperformed by the other approaches and the shrinkage approach (dashed-dotted line) gives results comparable to the Step 1 of the G1DBN procedure only. The results of the three methods are naturally improved for greater values of n but their relative perfomances are preserved (curves not shown).

We investigated the impact of the violation of the model assumptions. First we performed DBN inference on simulated data where the error covariance matrix Σ is not diagonal (3% of the coefficients outside the diagonal differ from 0) and the noise distribution is either Gaussian or uniform (U[-2; 2]). As shown on Figure 5B, the accuracy of the G1DBN procedure (Step 2) is not strongly affected when these assumptions on the noise distribution are not satisfied. However, it is difficult to get rid of the 1 st order Markov Assumption which was chosen in order to reduce the model dimension. When simulating an AR(2) model, the 2-order time dependencies existing in the model are missed. However, the 1-order time dependencies existing in the model are still recovered. Then, when considering a 2 nd order Markov process, an approximation can still be performed by successively inferring 1-and 2-order time dependencies. Note that the procedure also performs well when the number of parents in the true DAG G is greater than one (See Supp. Material [2], Section A).

Analysis of microarray time course data sets

Spellman's Yeast cell cycle data set We performed dynamic network inference from the Saccharomyces cerevisiae cell cycle data collected by Spellman et al. [START_REF] Spellman | Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization[END_REF]. We used the α Factor-based synchronization data (18 time points) and we focus here on a set of 786 genes which demonstrated consistent periodic changes in transcription level (See Supplementary Material [2], Section D.1 for more details). We carried out two surveys on this dataset. First, we allow only a subset of 18 genes5 identified as putative TFs to be possible parent genes (i.e. to have edges pointing out towards other genes in DAG G) and look for their target genes. Then we extend the search for parent genes to the whole dataset of 786 genes in a second survey. We set α 1 threshold for the G1DBN procedure according to guidelines detailed in Supplementary Material [2], Section B (α 1 = 0.1 for the 18 TF-survey, α 1 = 0.05 for the 786 TF-survey).

It is somehow difficult to assert the validity of the results obtained from real data as the whole regulatory machinery is not known yet. However the yeast cell cycle has been studied a lot and many regulation relationships have been recovered. We study the consistency of the first inferred edges with annotations in the Yeastract database [START_REF] Teixeira | The YEASTRACT database: a tool for the analysis of transcription regulatory associations in Saccharomyces cerevisiae [http://www.yeastract.com[END_REF], a curated repository currently listing found regulatory associations between TFs and target genes in S. cerevisiae.

In the 18 TF-survey, the first few selected edges are biologically validated. In the DAG comprising the 18 first selected edges (Figure 6A), 11 edges refer to identified regulatory relationships (thick blue edges). The first detected TFs are the genes coding for proteins FKH2, NDD1, RAP1 and SWI4. In particular, the proteins FKH2 (known as a TF with a major role in the expression of G2/M phase genes) and SWI4 (TF regulating late G1-specific transcription of targets) are pointed out as being essential TFs; they have the most target genes and the high majority (73%) of these regulatory relationships is listed in Yeastract.

As introduced in Section 3.3, we chose α 2 threshold in order to keep the False Discovery Rate (FDR) smaller than 1% with the approach by Benjamini and Hochberg [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF]. This lead to α 2 = 0.0059. The corresponding inferred DAG is shown in Figure 7. The two proteins FKH2 and SWI4 are still part of the TFs having the most targets, together with NDD1, which is an essential component of the activation of the expression of a set of late-S-phase-specific genes and TEC1, a transcription factor required for full Ty1 expression and Ty1-mediated gene activation (Ty transposableelement own for causing cell-type-dependent activation of adjacent-gene expression). The set of selected TFs is listed in Supplementary Material [2], Section D.2, Table 1, where the third column indicates the number of validated edges out of the selected ones. Except for NDD1, for which no target gene is listed in yeastract, one forth of the targets genes of the top four TFs are validated. For a comparative overview, the histogram of Figure 6B displays the percentage of validated edges out of the first 5 to 1000 selected edges inferred with each inference procedure When considering the 1000 first inferred edges, the results are very similar to what could be expected by chance only. Note that, as the Step 2 of G1DBN choose 308 edges only, it is not considered when comparing the 1000 first edges.

In the second survey including all the 786 genes as putative TFs, the dimension is far higher and the results are consequently more restricted. Indeed, the proportion of validated edges doesn't exceed 12.5%, obtained with the 2nd step of G1DBN procedure among the first selected edges. However, this is still a subtantial result as compared with the proportion of validated edges (equal to 0.26%). In order to keep the FDR smaller than 0.01, we chose α 2 = 0.0067 by following the Benjamini and Hochberg approach [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF]. The inferred DAG for the 786 TF-survey contains 437 genes and 380 edges. The display of this DAG, as well as the list of its edges and the list of the genes selected as TFs, is available in Supplementary Material [2].

Diurnal cycle on the starch metabolism of A. Thaliana We applied the G1DBN inference procedure to the expression time series data generated by Smith et al. [START_REF] Smith | Diurnal Changes in the Transcriptome Encoding Enzymes of Starch Metabolism Provide Evidence for Both Transcriptional and Posttranscriptional Regulation of Starch Metabolism in Arabidopsis Leaves[END_REF] to investigate the impact of the diurnal cycle on the starch metabolism of Arabidopsis Thaliana. We restricted our study to the 800 genes selected by Opgen-Rhein and Strimmer [START_REF] Opgen-Rhein | Learning causal networks from systems biology time course data: an effective model selection procedure for the vector autoregressive process[END_REF] as having periodic expression profiles 6 .

Using the heuristic approach detailed in Supplementary Material [2], Section B, we choose threshold α 1 = 0.02 allowing the distribution of the number of parents in the DAG G (1) having the number of 0-parent genes to dominate and the number of 1-parent genes to be half as large. We set α 2 = 0.005 in order to maintain the False Discovery Rate smaller than 0.01 by using the approach by Benjamini and Hochberg [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF] (See Section 3.3 for details). We recover the DAG in Figure 8 which has a "hub" connectivity structure. This network contains 206 edges implicating 277 different genes. We may notice that this DAG differs from the one inferred by Opgen-Rhein and Strimmer [START_REF] Opgen-Rhein | Learning causal networks from systems biology time course data: an effective model selection procedure for the vector autoregressive process[END_REF]. However the edges selected by the three inference procedures discussed in this section differ somewhat (See the proportion of edges in common by using the various inference approaches in Supplementary Material [2], Section C) and may, in fact, yield complementary information or insights.

Figure 8: DAG inferred with G1DBN from the data by Smith et al. [START_REF] Smith | Diurnal Changes in the Transcriptome Encoding Enzymes of Starch Metabolism Provide Evidence for Both Transcriptional and Posttranscriptional Regulation of Starch Metabolism in Arabidopsis Leaves[END_REF] in order to investigate starch metabolism of A. thaliana (LS estimates, α 1 = 0.1, α 2 =0.005 such that FDR< 0.01). The dark colored nodes are the 3 nodes with the most targets, 2 out of them are known for being implicated in starch metabolism. The light colored nodes are parent nodes already identified as TF or DNA binding protein (See Supp. Material[2], Section E, Table 2). This network contains 277 genes and 206 edges. See the edges Among the 'parent' nodes in the inferred DAG displayed in Figure 8, two nodes (799 and 628) out of the three having the most target refers to proteins that are known to be implicated in starch metabolism. Indeed, node 799, which has 14 'target' nodes, refers to DPE2 (DISPROPORTION-ATING ENZYME 2), which is an essential component of the pathway from starch to sucrose and cellular metabolism in plant leaves at night. Node 628 (6 targets) is a transferase (At5g24300) implicated in the starch synthase. Node 702, which is an unknown protein (At5g58220), has also 6 targets. These three nodes are dark-colored in the DAG of Figure 8. Note that there is no prior knowledge regarding the role of each gene (TF or target) in this survey. As a consequence, some edges might be inferred wrong way around7 . Thus node 799, which is a gene coding for an enzyme (DPE2), is most probably not a TF for its 14 apparent target genes. However node 799 is still the gene whose expression level best explains the expression of the 14 genes. Consequently these genes might be implicated in the same pathway as DPE2. The remaining parent nodes have from 1 to 4 targets. Among them, 9 genes, which are listed in Supplementary Material [2], Section E, Table 2, have already been identified as TFs or as DNA binding proteins. These 9 nodes are light-colored in the displayed DAG. Finally a list of 37 unknown proteins have been selected as parents in the inferred DAG. Potentially implicated in the regulation machinery of starch metabolism, these proteins represent a subset of genes which is relevant for further analyses. See more details on the inferred network displayed in Figure 8 in the Supplementary Material [2].

Discussion and conclusion

As more and more gene expression time series has become available, the need for efficient tools to analyze such data has become imperative. In this paper, we first determine sufficient conditions for Dynamic Bayesian Network modelling of gene expression time series. This type of modelling offers a straightforward interpretation: the edges of the DAG G defining the DBN exactly describe the set of conditional dependencies between successive gene expression levels. Having defined and characterized low order conditional dependence DAGs for DBNs, we point out relevant characteristics for the approximation of sparse DAGs. In particular, under faithfulness assumption, DAG G is included in the 1 st order conditional dependence DAG G (1) .

From these results, we develop G1DBN, a novel procedure for DBN inference, which makes it possible to tackle the 'small n, large p' estimation case that occurs with genetic time series data. Based on the consideration of low order conditional dependencies, the G1DBN procedure proved to be powerful on both simulated and real data analysis. With respect to other methods, the shrinkage approach considerably improves the precision of the overall estimation of the partial correlation coefficients when the number of observations n is small compared to the number of genes p. However, considering 1 st order conditional independence proved to be more efficient for DBN inference in terms of power and PPV on simulated data, and gave promising results on real data analysis. As for the lasso, one might notice that a drawback lies in the fact that the edge selection is done vertex by vertex whereas the DAG G is globally sparse but not uniformally. As a consequence, the lasso tends to uniformally reduce the number of parents of each vertex instead of only keeping the total number of edges contained. The power of the G1DBN procedure comes from the accuracy improvement of the testing made possible by the dimension reduction. Indeed, as the first step selection is based on the 1 st order conditional independence consideration, significance tests are performed in a model of dimension 4 (See Section 3.1). This represents a drastic dimension reduction compared to full order independence testing and makes the testing much more accurate. Thus, even if there are more edges in the DAG G (1) than in the true DAG G (Proposition 6), Step 1 of the procedure is already very predictive.

Throughout the analyses performed for this paper, we point out two major directions for further research. On the one hand, we noticed that the edges selected by the three inference procedures differ somewhat (See Supplementary Material [2], Section C). A further relevant study would consist of analyzing in which way these DBN inference procedures could have different strenghts and may be complementary. On the other hand, the use of robust estimators like Huber or Tukey bisquare did not allow a noticeable change of the inference approach on real data. Another interesting survey lies in the investigation of which measures of dependence, like non linear or other robust estimates, are the more pertinent to analyze gene expression data.

Proof of Lemma 2. Consider a discrete-time stochastic process X = {X i t ; i ∈ P, t ∈ N} whose joint probability P distribution has the density f with respect to Lebesgue measure on R p×n .

Let G 1 and G 2 be two different subgraphs of G full according to which the joint probability distribution P factorizes. Let i ∈ P , t ∈ N, we consider the random variable X i t . We denote as follows,

• the following subsets of P , pa 1 = {j ∈ P ; X j t-1 ∈ pa(X i t , G 1 )} pa 1 = P \{pa 1 } pa 2 = {j ∈ P ; X j t-1 ∈ pa(X i t , G 2 )} pa 2 = P \{pa 2 }

• and the densities of the joint or marginal probability distributions of (X i t , X t-1 ), g : R p+1 → R the density of the joint probability distribution of (X i t , X t-1 ), g i the density of the probability distribution of X i t , g P the density of the joint probability distribution of (X t-1 ), g i,pa 1 the density of the joint probability distribution of (X i t , X pa 1 t-1 ) where, X pa 1 t-1 = pa(X i t , G 1 ), g i,pa 2 the density of the joint probability distribution of (X i t , X pa 2 t-1 ) where X pa 2 t-1 = X t-1 \{pa(X i t , G 2 )}), etc... In the following, y ∈ R, x = (x 1 , ..., x p ) ∈ R p and we denote by x pa 1 = {x j ; j ∈ pa 1 } ∈ R |pa 1 | (Thus x = (x pa 1 , x pa 1 ) = (x pa 2 , x pa 2 ) ∈ R p ). As the probability distribution P factorizes according to G 1 , we derive from the DAG theory the conditional independence,

X i t ⊥ ⊥ X pa 1 t-1 |X pa 1 t-1 , that is, ∀y ∈ R, ∀x ∈ R p ,
g(y, x) g P (x) = g i,pa 1 (y, x pa 1 ) g pa 1 (x pa 1 ) .

Equivalent results derived from the factorization according to G 2 gives, ∀y ∈ R, x ∈ R p , N g i,pa 2 (y, x pa 2 ) = g i,pa 1 (y, x pa 1 ) g pa 1 (x pa 1 ) g pa 2 (x pa 2 ).

By taking the integral with respect to x pa 2 ∩pa 1 , we write for all y ∈ R, for all x pa 1 ∪pa 2 ∈ R |pa 1 ∪pa 2 | , g i,pa 2 (y, x pa 2 )d(x pa 2 ∩pa 1 ) = g i,pa 1 (y, x pa 1 ) g pa 1 (x pa 1 ) g pa 2 (x pa 2 )d(x pa 2 ∩pa 1 ) g i,pa 1 ∩pa 2 (y, x pa 1 ∩pa 2 ) = g i,pa 1 (y, x pa 1 ) g pa 1 (x pa 1 ) g pa 1 ∩pa 2 (x pa 1 ∩pa 2 )

Finally we have, ∀y ∈ R, ∀x ∈ R p , g(y, x) g P (x) = g i,pa 1 ∩pa 2 (y, x pa 1 ∩pa 2 ) g pa 1 ∩pa 2 (x pa 1 ∩pa 2 ) , that is the conditional density of the probability distribution of X i t given X t-1 is the conditional density of the probability distribution of X i t given X pa 1 ∩pa 2 t-1

. Then P factorizes according to G 1 ∩ G 2 .

Lemma 3 (Conditional independence between non adjacent successive variables)

Let G be a subgraph of G full according to which the probability distribution P admits a BN representation. For any pair of successive variables (X j t-1 , X i t ) which are non adjacent in G, we have X i t ⊥ ⊥ X j t-1 | pa(X i t , G) and X i t ⊥ ⊥ X j t-1 | pa(X i t , G) ∪ S, for all S subset of {X k u ; k ∈ P, u < t}. As an illustration of Lemma 3, assume P admits a BN representation according to the DAG of Figure 1C. There is no edge between X 3 t and X 1 t+1 in this DAG. Now consider in Figure 2B the moral graph of the smallest ancestral graph containing X 3 t , X 1 t+1 and the parents (X 1 t , X 2 t ) of X 1 t+1 . The set (X 1 t , X 2 t ) blocks all paths between X 3 t and X 1 t+1 . From Proposition 2, we have X 1 t+1 ⊥ ⊥ X 3 t | pa(X 1 t+1 , G).

Proof of Lemma 3. Assume P admits a BN representation according to G, a subgraph of G full . Let X j t-1 and X i t be two non adjacent vertices of G (there is no edge between them in G) and consider the moral graph (G An(X i t ∪X j t-1 ∪pa(X i t ,G)) ) m of the smallest ancestral set containing the variables X i t , X j t-1 and the parents pa(X i t , G) of X i t in G. As DAG G is a subgraph of G full , the set of parents pa(X i t , G) blocks all paths between X j t-1 and X i t in the moral graph (G An(X i t ∪X j t-1 ∪pa(X i t ,G)) ) m . From Proposition 2, this establishes the conditional independence X i t ⊥ ⊥ X j t-1 | pa(X i t , G). This result holds for the conditioning according to any subset S ⊆ {X k u ; k ∈ P, u < t}.

Figure 1 :

 1 Figure 1: (A) A biological regulation motif. (B)The concentration graph corresponding to the motif A. For all i ≥ 3, Y i is a Gaussian variable representing the expression level of gene G i . Some cycles cannot be represented on the concentration graph. (C) Dynamic network equivalent to the regulation motif A. Each vertex X i t represents the expression level of gene G i at time t. This graph is acyclic and allows to define a Bayesian network.

Figure 2 :

 2 Figure 2: (A) Moral graph of the DAG in Figure1C. For all t > 1, the parents of the variable X 1 t are 'married', that is connected by an undirected edge. (B) Moral graph of the smallest ancestral set containing the variables X 1 t+1 , its parents in the DAG in Figure1Cand X 3 t . As the set (X 1 t , X 2 t ) blocks all paths between X 3 t and X 1 t+1 , thus {X 1 t , X 2 t } separates X 1 t+1 from X 3 t and we have X 1 t+1 ⊥ ⊥ X 3 t | (X 1 t , X 2 t ).

Figure 3 :

 3 Figure3: First-order conditional dependence DAG G(1) (obtained from the DAG in Figure1C). The spurious dashed arrow may appear in G(1) .

Figure 4 :

 4 Figure 4: Outline of the 2 step-procedure G1DBN for DBN inference.

Figure 5 :

 5 Figure 5: Precision-Recall (PR) curves obtained for network inference from simulated data (n = 20). (A) Comparison of the inference procedures: G1DBN (LS or Tukey), shrinkage and lasso.Step 2 of the G1DBN approach drastically improves the results (threshold α 1 = 0.7). (B) Impact of noisy data, simulated using a non diagonal matrix Σ with either Gaussian or uniform noise, on the G1DBN procedure (Step 2) computed with LS estimates.

Figure 6 :

 6 Figure 6: Some results of the 18 TF-survey of S. cerevisiae cell cycle. (A) DAG containing the 18 first selected edges with G1DBN with LS estimates (PPV=60%). Colored nodes represent the TFs and the dark blue edges are validated by the Yeastract database. (B) Percentage of validated edges out of the first 5 to 1000 edges inferred with the G1DBN procedure, after Step 2 or after Step 1 only, the shrinkage or the lasso procedure. The dashed line shows the proportion of validated edges out of the 786×18 possible edges.

Figure 7 :

 7 Figure 7: DAG inferred by G1DBN with LS estimates, using α 1 = 0.1, α 2 = 0.0059 (ensuring FDR< 0.01), in the 18 TF-survey of the S. cerevisiae cell cycle. The 17 colored nodes represent the 16 TFs selected as parent node out of the 18 TFs under study, plus node FKH1 which is selected as a target of NDD1. The dark blue edges are validated by Yeastract. This network contains 286 genes and 308 edges. See the complete edges list in Supp. Material [2].

Table 1

 1

	: Notations

available at http://strimmerlab.org/software.html.

The 18 genes code for proteins ACE2, FKH1, FKH2, GAT3, MBP1, MCM1, MIG2, NDD1, PHD1, RAP1, RME1, STB1, SUT1, SWI4, SWI5, SWI6, TEC1 and YOX1. consist of the overlap between the 786 genes under study and the 50 genes identified as putative TFs in a recent study by Tsai et al.[START_REF] Tsai | Statistical methods for identifying yeast cell cycle transcription factors[END_REF].

The data are available in the GeneNet R package at http://strimmerlab.org/ software/genenet/html/ar th800.html or in our R package G1DBN (arth800line).

In particular if some assumption of the model is not satisfied. For instance if an essential TF is missing or if the regulation is not transcriptional, i.e. does not depend on the amount of mRNA coding for the protein.

APPENDIX A Proofs

A.1 Lemmas 1 to 3 and proofs Lemma 1 Under Assumptions 1 and 2, the probability distribution P admits a DBN representation according to a DAG whose edges only join nodes representing variables observed at two successive time points, at least according to DAG G full = (X, {(X j t-1 , X i t )} i,j∈P,t>1 ) which has edges between any pair of successive variables.

Proof of Lemma 1. From assumption 1, the density f of the joint probability distribution of process X be written as the product of conditional densities,

where f (X t |X t-1 ) refers to the density of the conditional probability distribution of X t given X t-1 .

From Assumption 2, for all t > 1, the conditional density f (X t |X t-1 ) can be written as the product of the conditional density of each variable X i t given the set of variables X t-1 observed at the previous time,

From equations [START_REF] Efron | Local false discovery rates[END_REF] and [START_REF] Efron | Least angle regression[END_REF], the density f writes as the product of the conditional density of each variable X i t given its parents in G full . From Proposition 1, the probability distribution P admits a BN representation according to G full .

Lemma 2 Assume the joint probability distribution P of process X has density f with respect to Lebesgue measure on R p×n . If P factorizes according to two different subgraphs of G full , G 1 and G 2 , then P factorizes according to

From Lemma 2, it is straightforward that, among the DAGs included in G full , there exists a minimal DAG (denoted by G in the paper) according to which the probability distribution P factorizes, thus establishing a BN representation of process X.

A.2 Proof of Propositions 3, 5, 6 and 7

Proof of Proposition 3. First, we show that P admits a BN representation according to G. Let i, j ∈ P such that X i t ⊥ ⊥ X j t-1 |X

). Under Assumptions 1 and 2, from Lemma 1 (See Appendix A.1) and Prop. 1, P admits a BN representation according to the DAG (X, E(G full ) \ (X j t-1 , X i t )) which has the edges of G full except for the edge (X j t-1 , X i t ). This holds for any pair of successive variables that are conditionally independent.

From Lemma 2 (See Appendix A.1), P admits a BN representation according to the intersection of the DAG (X, E(G full ) \ (X j t-1 , X i t )) for any pair (X i t , X j t-1 ) such that X i t ⊥ ⊥ X j t-1 |X

, that is DAG G. Also, DAG G cannot be reduced. Indeed, let (X l t-1 , X k t ) be an edge of G and assume that P admits a BN representation according to G\(

Proof of Proposition 5.

First, from Corollary 1, G ⊇ G (1) . Second, let X be a Gaussian process and (X j t-1 , X i t ) ∈ E( G), then according to Proposition 3, X i t ⊥ ⊥ X j t-1 | X 1) ). We are going to prove that this contradicts the nullity of covariance Cov(X i t , X j t-1 |X

) can be written,

.

However both terms in the latter expression of Cov(ij|k, l) are null:

• as N Max pa ( G) ≤ 1, X j t-1 is the only parent of X i t in G. So the variable X j t-1 and thus also the set (X j t-1 , X k t-1 ) blocks all paths between X l t-1 and X i t in the moral graph of the smallest ancestral set containing X i t ∪ X j,k,l t-1 . Then we have,

Then Cov(ij|k, l) = 0. By induction, we obtain Cov(X i t , X j t-1 |X P j t-1 ) = 0 leading to a contradiction with (X j t-1 , X i t ) ∈ E( G). Therefore (X j t-1 , X i t ) ∈ G (1) and we have G ⊆ G (1) .

) then there exists a subset of q variables X Q t-1 with respect to which X j t-1 and X i t are conditionally independent. From faithfulness, the subset X Q t-1 separates X j t-1 and X i t in the moral graph of the smallest ancestral set containing X i t ∪ X j t-1 ∪ X Q t-1 . This contradicts the presence of the edge (X j t-1 , X i t ) in G.

Proof of Prop 7 .

From faithfulness, G ⊆ G (q) . Then for all i in P , for all t > 1, we have N pa (X i t , G) ≤ N pa (X i t , G (q) ) ≤ q. From Proposition 4, (X j t-1 , X i t ) / ∈ E( G) ⇒ (X j t-1 , X i t ) / ∈ E(G (q) ), that is (X j t-1 , X i t ) ∈ E(G (q) ) ⇒ (X j t-1 , X i t ) ∈ E( G).