Océane Cassan

Charles-Henri Henri Lecellier

Antoine Martin

Laurent Brehelin
email: brehelin@lirmm.fr

Sophie Lèbre
email: sophie.lebre@umontpellier.fr

Laurent Bréhélin

Optimizing data integration improves Gene Regulatory Network inference in Arabidopsis thaliana

Keywords: Integrative regression-based Gene Regulatory Network inference, data integration optimization, simulated null hypothesis, weighted Random Forests, weighted LASSO, Arabidopsis thaliana

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Gene Regulatory Network (GRN) inference has the objective of deciphering the relationships between genes in the context of transcription, which can provide invaluable insight into environmental adaptation or developmental processes in living organisms. Statistical inference methods usually leverage highthroughput genomics to reconstruct those networks, in which nodes represent genes, and edges represent a relation of regulation between those genes. Because transcriptomic data are increasingly common and less costly, they are the input of choice for most statistical approaches to GRN inference.

Regression-based techniques mainly differ in their choice of regression function to link the expression of a target gene to the expression of its regulators. For example, TIGRESS [START_REF] Haury | TIGRESS: Trustful inference of gene REgulation using stability selection[END_REF], MERLIN [START_REF] Roy | Integrated Module and Gene-Specific Regulatory Inference Implicates Upstream Signaling Networks[END_REF] or The INFERELATOR [START_REF] Bonneau | The Inferelator: an algorithm for learning parsimonious regulatory networks from systems-biology data sets de novo[END_REF][START_REF] Madar | The Inferelator 2.0: a scalable framework for reconstruction of dynamic regulatory network models[END_REF][START_REF] Gibbs | High-performance single-cell gene regulatory network inference at scale: the inferelator 3.0[END_REF] techniques implement linear parametric models for this task, while GENIE3 [START_REF] Vân | Inferring regulatory networks from expression data using tree-based methods[END_REF] and inspired works [START_REF] Geurts | dyngenie3: dynamical genie3 for the inference of gene networks from time series expression data[END_REF][START_REF] Petralia | Integrative random forest for gene regulatory network inference[END_REF][START_REF] Cirrone | OutPredict: multiple datasets can improve prediction of expression and inference of causality[END_REF][START_REF] Cassan | Inferring and analyzing gene regulatory networks from multi-factorial expression data: a complete and interactive suite[END_REF] model non-linear relations via Random Forests (RFs) or, more broadly, ensembles of trees. Once regression models are fit, they allow the extraction of the influence of each regulator over each gene, and the strongest pairs are assembled to form a final sparse GRN.

Given the under-determined nature of GRN inference from expression alone, using additional sources of data can guide the choice between several regulators explaining expression data equally well. Complementary omics have already been used in addition to gene expression to enhance GRN inference, such as TF binding experiments (mostly ChIP-Seq) or Transcription Factor Binding Motifs (TFBM) [START_REF] Kundaje | Learning regulatory programs that accurately predict differential expression with MEDUSA[END_REF][START_REF] Marbach | Predictive regulatory models in drosophila melanogaster by integrative inference of transcriptional networks[END_REF][START_REF] Aibar | SCENIC: single-cell regulatory network inference and clustering[END_REF][START_REF] De Clercq | Integrative inference of transcriptional networks in arabidopsis yields novel ROS signalling regulators[END_REF][START_REF] Qin | Inferring gene regulatory networks by integrating ChIP-seq/chip and transcriptome data via LASSOtype regularization methods[END_REF][START_REF] Cirrone | OutPredict: multiple datasets can improve prediction of expression and inference of causality[END_REF][START_REF] Gibbs | High-performance single-cell gene regulatory network inference at scale: the inferelator 3.0[END_REF], knock-outs and protein-protein interactions [START_REF] Petralia | Integrative random forest for gene regulatory network inference[END_REF] or chromatin accessibility [START_REF] De Clercq | Integrative inference of transcriptional networks in arabidopsis yields novel ROS signalling regulators[END_REF][START_REF] Emily R Miraldi | Leveraging chromatin accessibility for transcriptional regulatory network inference in T helper 17 cells[END_REF].

In a linear context, prior information can be integrated to GRN inference by modulating the penalty strength for each TF during the estimation of regularized models with a weighted version [START_REF] Cecilie Bergersen | Weighted lasso with data integration[END_REF] of the LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and many variations (e.g. ElasticNet [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]) [START_REF] Christley | Incorporating existing network information into gene network inference[END_REF][START_REF] Greenfield | Robust data-driven incorporation of prior knowledge into the inference of dynamic regulatory networks[END_REF][START_REF] Studham | Functional association networks as priors for gene regulatory network inference[END_REF][START_REF] Emily R Miraldi | Leveraging chromatin accessibility for transcriptional regulatory network inference in T helper 17 cells[END_REF][START_REF] Gibbs | High-performance single-cell gene regulatory network inference at scale: the inferelator 3.0[END_REF] or by making use of a Bayesian prior [START_REF] Greenfield | Robust data-driven incorporation of prior knowledge into the inference of dynamic regulatory networks[END_REF][START_REF] Alireza | A prior-based integrative framework for functional transcriptional regulatory network inference[END_REF][START_REF] Gibbs | High-performance single-cell gene regulatory network inference at scale: the inferelator 3.0[END_REF]. An in-depth approach explored a resolutive range of data integration strengths with a weighted ElasticNet, and then choose a value maximizing effective data integration [START_REF] Greenfield | Robust data-driven incorporation of prior knowledge into the inference of dynamic regulatory networks[END_REF]. More recently, the StARS approach proposing to increase the stability of non-oriented GRN inference in Gaussian Graphical Models [START_REF] Liu | Stability approach to regularization selection (stars) for high dimensional graphical models[END_REF] was adapted for the LASSO in order to select a small subset of robust regulators for each gene in oriented GRN inference [START_REF] Emily R Miraldi | Leveraging chromatin accessibility for transcriptional regulatory network inference in T helper 17 cells[END_REF]. In that work, the integrated priors are TFBMs in accessible chromatin, and three values of prior reinforcement modulating penalty strength are investigated. The moderate one is chosen to maximise the area under the precision and recall curve against a CHIP-Seq gold standard.

Regarding non-linear regression, iRafNet [START_REF] Petralia | Integrative random forest for gene regulatory network inference[END_REF] proposed a Random Forest (RF) based procedure. It consists in weighting the random sub-sampling of regulators during trees elongation so that the chance of regulators supported by prior knowledge to get chosen at decision nodes is increased. This has the effect of inflating the importance metric of interactions supported by the chosen prior [START_REF] Petralia | Integrative random forest for gene regulatory network inference[END_REF]. The weights controlling the contribution of prior data to expression is provided by a predefined function, specific to each type of prior, but without specific tuning. This strategy was further adapted to time series expression data in the OutPredict method [START_REF] Cirrone | OutPredict: multiple datasets can improve prediction of expression and inference of causality[END_REF] as an extension of both a dynamic version of GENIE3 [START_REF] Geurts | dyngenie3: dynamical genie3 for the inference of gene networks from time series expression data[END_REF] and iRafNet.

Existing integrative regression models have a great potential to predict GRN from several types of omics. However, fine tuning the contribution of prior data to expression is rarely explored. When it is, the choice of prior integration strength is set using of a gold standard that is either identical to the integrated prior [START_REF] Greenfield | Robust data-driven incorporation of prior knowledge into the inference of dynamic regulatory networks[END_REF], or of a related nature [START_REF] Emily R Miraldi | Leveraging chromatin accessibility for transcriptional regulatory network inference in T helper 17 cells[END_REF] (in that latter paper, the prior information of binding motifs is necessarily correlated to ChIP-seq validation data).

In this work, we propose DIOgene (Data Integration Optimization for gene networks) for a more robust and independent calibration of prior data integration strength, based on effective data integration, gene expression prediction accuracy, and a simulated null hypothesis. Moreover, to reflect the specificities of gene regulation, we propose to tune data integration strength specifically for each gene, in contrast to previous works that traditionally enforce the same integration intensity for all genes [START_REF] Alireza | A prior-based integrative framework for functional transcriptional regulatory network inference[END_REF][START_REF] Greenfield | Robust data-driven incorporation of prior knowledge into the inference of dynamic regulatory networks[END_REF][START_REF] Petralia | Integrative random forest for gene regulatory network inference[END_REF][START_REF] Emily R Miraldi | Leveraging chromatin accessibility for transcriptional regulatory network inference in T helper 17 cells[END_REF][START_REF] Cirrone | OutPredict: multiple datasets can improve prediction of expression and inference of causality[END_REF][START_REF] Gibbs | High-performance single-cell gene regulatory network inference at scale: the inferelator 3.0[END_REF].

In order to represent the most common methods in the field of integrative GRN inference both in the linear and non-linear cases, we illustrate our results using a weightedLASSO and weightedRF model. The proposed approaches are applied to the modelling of the transcriptomic response to nitrate induction in the roots of Arabidopsis thaliana [START_REF] Varala | Temporal transcriptional logic of dynamic regulatory networks underlying nitrogen signaling and use in plants[END_REF] using TFBM in gene promoters. as prior information. With this study, we hope to open a reflection about data integration and evaluation practices in the field.

Material and Methods

We present in the following sections weightedLASSO and weightedRF, two integrative GRN inference procedures based on TFBM prior information (Figure 1). The TFBM prior matrix Π gives, for each regulator-target pair (r, t) a prior value Πr,t ∈ [0, 1] defined as: Πr,t = 0 : if the motif of r is not in the promoter of t 1 : if the motif of r is in the promoter of t 1 2 : if the motif of r is unknown [START_REF] Aibar | SCENIC: single-cell regulatory network inference and clustering[END_REF] Throughout this study, the parameter α will be used to tune data integration strength in GRN inference : its value ∈ [0, 1] controls the contribution of TFBM information to expression data. When α = 0, expression alone is used, while at α = 1, gene expression is used to chose only between regulators possessing a TFBM in the target gene. The integration of α in the two inference algorithms is detailed below, as well as our original permutation based-procedure to determine the optimal α, in a gene-specific manner. The Implementation of this integration tuning in weightedLASSO and weightedRF is detailed below.

Generalized Linear Model with Weighted LASSO (weightedLASSO)

As RNA-Seq experiments generate count data, we model the expression of a target gene t in the condition i as a Poissondistributed variable Yt,i ∼ P(µt,i). The parameter of the Poisson distribution µt,i is estimated on the log scale as a linear combination of the expression values of the regulator genes, with xr,i the expression level of regulator r in condition i:

ln(µt,i) = βt,0 + Σ R r=1 βt,rxr,i (2)
We employ a LASSO penalty in order to overcome the highdimensional setting and to select the most predictive regulators. In addition, we propose to use differential shrinkage in order to favor the selection of TFBM-supported variables. Differential shrinkage for the LASSO (weighted LASSO) allows to modulate the penalty strength of each variable individually, in a way that regulators with a binding motif in the target's promoter are less penalized during model adjustment. We model differential shrinkage using specific penalty coefficients wt,r ∈ [0, 1] defined as a linear function of the TFBM prior Πt,r and α (Figure 2a):

w LASSO t,r = 1 -Πt,rα (3)
For each target gene t, the function to minimize for model estimation is thus:

argmin β t -logL(βt; X, Yt) + λΣ R r=1 w LASSO t,r |βt,r| , (4)
where λ controls the overall strength of the penalty, R is the total number of regulator genes, N is the total number of expression measurements for t logL(βt; X, Yt) is the log-likelihood function. The value of λ is learned from 5-fold cross validation. We relied on the glmnet [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF] implementation of the LASSO, with the penalty.factor argument specifying differential shrinkage weights. In practice, for genes with a number P of TFBMsupported regulators exceeding the total number of experiments considered N , the model cannot be estimated at α = 1. We thus set α to 1 -ϵ instead of 1, with ϵ = 10 -4 .

To further reduce over-fitting problems and improve robustness, we also included a bootstrap procedure to the LASSO [START_REF] Francis R Bach | Bolasso: model consistent lasso estimation through the bootstrap[END_REF][START_REF] Meinshausen | Stability selection[END_REF][START_REF] Liu | Stability Approach to Regularization Selection (StARS) for High Dimensional Graphical Models[END_REF], as already used in some GRN inference approaches using a linear Input data for a given target gene is the vector of expression of this gene Y , an expression matrix for the regulators X, and a TFBM scoring matrix Π. Π contains information about the presence or absence of the regulators PWM in target gene's promoter. α is a parameter controlling the force of TFBM integration to expression. For each target gene, a regression model is fit to Y with weightedRF or weightedLASSO using X as predictive variables, and favoring the attribution of high influences to regulators with high values in Π. For weightedRF, this prioritization is executed by a weighted subsampling of the regulators space when elongating the regression trees, while it is achieved by differential shrinkage combined to stability selection is weightedLASSO. Once all regulator-gene pairs have been ranked based on their influence in the regression models, final GRN are built by selecting the number of strongest interactions providing a desired network density.

model [START_REF] Haury | TIGRESS: Trustful inference of gene REgulation using stability selection[END_REF][START_REF] Emily R Miraldi | Leveraging chromatin accessibility for transcriptional regulatory network inference in T helper 17 cells[END_REF]. Hence, instead of fitting a single generalized linear model, S models are adjusted on a bootstrapped version of the data as follows:

1. N observations are sampled with replacement from the N available experimental conditions (bootstrapping). 2. The N bootstrapped observations are randomly partitioned into 5-cross validation folds. We ensure that duplicated observations during bootstrapping are grouped within the same fold. 3. A model is fitted by minimizing Equation (4) during crossvalidation, allowing to learn the value of the sparsity, λ1se. λ1se is the largest value of λ in the λ grid less than one standard deviation away from the value of λ minimising prediction error on the cross-validation test folds.

In this study, results are presented for S = 50.

Weighted Random Forests (weightedRF)

Non-linear regressors such as ensembles of regression trees model combinatorics of regulators and complex relations between a target gene and the expression of its regulators. As inspired by iRafNet [START_REF] Petralia | Integrative random forest for gene regulatory network inference[END_REF], we model data integration in RF by increasing the use of regulators supported by a binding motif in the decision nodes of the regression trees. A weightedRF is inferred for each target gene t. At each decision node, the most discriminating regulator is chosen among a subset of √ R regulators. This subset, traditionally equiprobably sampled from all the regulators, is submitted here to a weighted sampling where the weights encode prior knowledge about the regulators, growing with data integration strength α and prior value Πr,t. Regulators with a high prior value are more likely to get chosen among the √ R regulators tested to create a decision node. More formally, we define that the chance of a regulator r to get picked in the decision node for the target gene t is proportional to the weight w RF r,t (Figure 2b):

w RF r,t = -1 -(α -1) 2 + 1 if Πr,t = 0 1 -α if Πr,t = 1 2 1 -(α -1) 2 + 1 if Πr,t = 1 (5)
In practice, these functions allow that in both cases, when α = 1, variable selection is performed among regulators with Πrt = 1 only (Figure S1). One weighted tree is learned for each of S bootstrapped samples. In this study, results are presented for S = 2000. The S trees are then aggregated into one RF per target gene. We implemented the new weight function upon the iRafNet R and C++ code, along with the possibility to restrict variables to regulators in regressions, which was not permitted initially, and Fig. 2: Functions linking prior TFBM scores (Πrt) to regulator-specific weights during model estimation, depending on integration strength α. a. In weightedLASSO, the penalty strength of the LASSO decreases with α when Πrt > 0 (Equation 3). b. In weightedRF, the sampling weight at regression tree nodes increases with α when Πrt = 1, and decreases otherwise (Equation 5). Weights are normalized between 0 and 1. In practice, these functions allow that in both cases, when α = 1, variable selection is performed among regulators with Πrt = 1 only (Figure S1).

with the addition of a new feature importance metric (Equation 8).

Error in predicting gene expression (MSE)

We assess the ability of weightedRF and weightedLASSO to accurately predict gene expression via the Mean Squared Error (MSE) metric. It is measured on Out Of Bag conditions (OOB) that were left out of the bootstrapped samples (and consequently not used in weightedRF and weightedLASSO models training). For weightedLASSO, the MSE for a target gene t, and a given α is

M SEtα = 1 S Σ S m=1 1 N OOB m Σ i∈OOB m (ytiα -ŷmtiα) 2 (6)
where m corresponds to one of the S bootstrapped LASSO models, OOBm refers to all conditions i that are OOB for m, N OOB m refers to the cardinal of OOBm, ŷmtiα is the prediction of m on the condition i ∈ N OOB m .

For weightedRF, the MSE for a target gene t and α is

M SEtα = 1 N Σ N i=1 1 S OOB i Σ m∈OOB i (ytiα -ŷmtiα) 2 (7)
where i corresponds to one of the N experimental conditions, OOBi refers to all trees m for which i is OOB, S OOB i refers to the cardinal of OOBi, ŷmtiα is the prediction of m ∈ S OOB i on the condition i.

Importance of a regulator

Previous RF-based models [START_REF] Vân | Inferring regulatory networks from expression data using tree-based methods[END_REF], including iRafNet [START_REF] Petralia | Integrative random forest for gene regulatory network inference[END_REF] which inspired weightedRF, used the traditional node purity metric. However, it is tailored for tree-based approaches only, and was shown to be not suitable to interpret variable importance in RFs in the presence of dependencies and interactions [START_REF] Scornet | Trees, forests, and impurity-based variable importance[END_REF][START_REF] Nicodemus | Predictor correlation impacts machine learning algorithms: implications for genomic studies[END_REF]. We therefore define a common importance metric for both for weightedLASSO and weightedRF. This importance measure between a regulator r and its target t is given by the relative increase of MSE measured on the OOB (Equation 6 and 7), induced by shuffling the expression values of r when making the prediction ŷmtiα :

Importance rtα = M SE tα,shuffle(r) -M SEtα M SE tα,shuffle(r) . (8)
We normalize the MSE difference by M SE tα,shuffle(r) to ensure that this statistic is comparable between different target genes and is included in the [0, 1] interval as similarly proposed in previous studies [START_REF] Greenfield | Robust data-driven incorporation of prior knowledge into the inference of dynamic regulatory networks[END_REF][START_REF] Gibbs | High-performance single-cell gene regulatory network inference at scale: the inferelator 3.0[END_REF].

Effective Data Integration (EDI)

In order to measure the direct consequence of modulating data integration through α, we introduce the notion of Effective Data Integration (EDI), that reflects the importance of TFBMsupported regulators in the predictions of a regression model. For a target gene t, regulators are ranked by increasing values of importance, and the EDI is the average position in this ranking of TFBM-supported regulators, i.e the regulators for which Πr,t = 1.

EDItα = Σ Π r,t =1 Rank(Importance rtα ) # Π r,t =1 (9)
EDI is close to 1 (resp. R, the total number of regulators) when all regulators with a motif have low (resp. high) importance. We expect that increasing α will increase the importance values of TFBM-supported regulators, and thus increase EDI.

Gene-specific optimisation of α (DIOgene)

Recall that both in linear and non linear models, α is the parameter controlling the extent to which TFBM contribute to model estimation relatively to gene expression, with possible values between 0 and 1. Choosing the value of α and consequently acting upon EDI is instrumental: it reflects strong modelling assumptions and has tangible impacts on inferred GRN. Given that enforcing data integration interferes with model estimation based solely on the minimization of the error, we can expect that this may also deteriorate the predictions accuracy. The foundation of DIOgene is that we should integrate the prior information of TFBM only if the prediction performance is not deteriorated too much.

In order to define what is an acceptable loss of prediction performance, we create a synthetic null hypothesis that will provide a reference for comparison. In this simulated null dataset, we break the link between gene expression and TFBM scores by randomly unmatching the expression profiles between regulators. A regulator then keeps its correct TFBM, but is attributed the wrong expression profile. In such a synthetic baseline, there is theoretically no joint information to be learned from the combination of expression and TFBM, and increasing data integration strength can only provide uninformative TF-gene interactions.

We can then assess prediction error (Equations 6, 7) on true data relatively to this synthetic dataset at comparable EDI (Equation 9) values for various value of integration strength α. In order to identify the appropriate amount of TFBM knowledge to inform GRN inference, we propose that the optimal value of α for target t (hereafter denoted as αt,opt), is chosen where true prediction error is most reduced as compared to the error committed under the simulated null hypothesis. This corresponds to a level of data integration where TFBM incorporation in the model provides a sufficient improvement of prediction over the shuffled baseline.

In order to build these curves, we run weightedRF and weightedLASSO for values of α ranging from 0 to 1 with a step of 0.1, and collect the variations of EDI and MSE depending on α. In order to capture the whole range of possible variations in the simulated null data but also the stochasticity inherent to the models, we run a large number of repetitions of weightedRF and weightedLASSO (resp. 100 and 50 repetitions) on both the true and simulated datasets. The curves M SEt(α) and EDIt(α) can then be assembled into one M SEt = ft(EDIt) curve to optimize α (see several gene examples in Figure 3 for weightedRF, and Figure S2 for weightedLASSO).

Formally, the normalized difference in MSE between true and shuffled datasets for a value of α is measured by ∆tα

∆tα = µ MSE,shuffle (EDItα) -µ MSE,true (EDItα) σ MSE,true (EDItα) (10)
with µ MSE,true (EDItα) the mean MSE at EDItα, σ MSE,true (EDItα) the standard deviation of MSE at EDItα, µ MSE,shuffle (EDItα) the mean MSE on the null dataset interpolated at EDItα, αt,opt is then the value of α that maximizes ∆tα on a certain condition:

αt,opt =      0 if max α∈[0,1] (∆tα) ≤ 1 argmax α∈[0,1] (∆ tα) otherwise. ( 11
)
When max

α∈[0,1]
(∆tα) ≤ 1, no level of data integration seems appropriate with regard to the shuffled null data, hence αt,opt is set to 0. As an illustration, the vertical blue line in Figure 3 (c-f) indicates αt,opt in 4 gene examples.

Edges selection

For each target gene t, a regression model is learned with αt,opt (Equation 11). Then, a final sparse GRN is built by selecting the edges associated with the strongest regulator-target interactions as given by the importance metric (Equation 8). Density is a common topological descriptor of biological networks: the sparsity of GRN confers them low density values, typically between 0.001 and 0.1 [START_REF] Leclerc | Survival of the sparsest: Robust gene networks are parsimonious[END_REF][START_REF] Koutrouli | A Guide to Conquer the Biological Network Era Using Graph Theory[END_REF][START_REF] Campos | Evolutionary constraints on the complexity of genetic regulatory networks allow predictions of the total number of genetic interactions[END_REF][START_REF] Hayes | Graphletbased measures are suitable for biological network comparison[END_REF]. A classical strategy in GRN inference is to select edges satisfying a biologically relevant user-specified network density defined as

D = E E total
, with E being the number of edges in the inferred network, and E total = R(T -1) being the total number of edges in a complete oriented GRN containing R regulators and T genes [START_REF] Cassan | Inferring and analyzing gene regulatory networks from multi-factorial expression data: a complete and interactive suite[END_REF]. The number of top-ranked edges to select in order to satisfy a density D is thus

E = DR(T -1). ( 12
)

Evaluation against experimental interactions

Let's define G as the set of experimental regulatory interactions (gold standard) restricted to interactions involving genes given as input for GRN inference. E is the set of inferred oriented interactions restricted to TFs studied in the gold standard. The other inferred interactions can neither be confirmed nor falsified, and are thus not taken into account here. Two standard metrics for assessing the quality of a GRN are:

1. Precision, the fraction of edges in E present in G:

Precision = |E ∩ G| |E| .
2. Recall, the fraction of edges in G retrieved by GRN inference:

Recall = |E ∩ G| |G| .

Experimental data analysis in Arabidopsis

Gene expression dataset

As a case study for GRN inference, we chose the transcriptomic root response to nitrate induction in the model plant Arabidopsis thaliana [START_REF] Varala | Temporal transcriptional logic of dynamic regulatory networks underlying nitrogen signaling and use in plants[END_REF]. This dynamic response has the advantage of being already well characterized, and used in other previous developments to chart regulatory networks [START_REF] Varala | Temporal transcriptional logic of dynamic regulatory networks underlying nitrogen signaling and use in plants[END_REF][START_REF] Brooks | Network walking charts transcriptional dynamics of nitrogen signaling by integrating validated and predicted genome-wide interactions[END_REF][START_REF] Cirrone | OutPredict: multiple datasets can improve prediction of expression and inference of causality[END_REF]. Continuing efforts to uncover these regulatory mechanisms is of great agricultural interest, as nitrate is the main source of nitrogen used by most plants. Gene expression was measured in seedling roots at times 0, 5, 20, 30, 45, 60, 90, and 120 minutes after nitrate or control treatments1 . Each combination of time point and treatment was measured in three replicates, resulting in a total of N = 45 samples. The RNA-Seq raw counts were normalized via the TMM method [START_REF] Mark | A scaling normalization method for differential expression analysis of RNA-seq data[END_REF], and lowly expressed genes were removed prior to differential expression analysis. We selected differentially expressed genes responding to nitrate induction in time by testing the interaction terms between nitrate treatment and time modelled as natural splines, as proposed in the original article that produced this dataset [START_REF] Ritchie | limma powers differential expression analyses for RNA-sequencing and microarray studies[END_REF][START_REF] Varala | Temporal transcriptional logic of dynamic regulatory networks underlying nitrogen signaling and use in plants[END_REF]. A total of T = 1, 426 genes had FDR adjusted p-values under 0.05. Among those 1,426 genes, R = 201 are annotated as transcriptional regulators. These nitrate-responsive genes and regulators are taken as input for GRN inference (Tables S1 andS2).

TFBM dataset

We define a promoter sequence in Arabidopsis as the sequence spanning -1,000 bp and +200 bp around gene start as defined in TAIR10, because this interval has been estimated to contain 86% of binding sites in plants [START_REF] Yu | Positional distribution of transcription factor binding sites in arabidopsis thaliana[END_REF]. TFBM information, encoded by Position Weight Matrices (PWM), was retrieved from the JASPAR database [START_REF] Castro-Mondragon | JASPAR 2022: the 9th release of the openaccess database of transcription factor binding profiles[END_REF] and The Plant Cistrome Database [START_REF] Ronan | Cistrome and epicistrome features shape the regulatory DNA landscape[END_REF]. Among the 201 nitrate-responsive regulators, 70 regulators were associated with a known PWM in the union of these databases. FIMO [START_REF] Grant | FIMO: scanning for occurrences of a given motif[END_REF] was used to find occurences of these 70 PWMs with a significance threshold of 1e -4 in promoter sequences. When several occurences were identified in one promoter, the maximum score was kept.

The TFBM occurrences within nitrate-responsive gene promoters form a prior TFBM network between TFs and targets Fig. 3: Gene-specific data integration with DIOgene is tuned by monitoring model performance variation relatively to a synthetic null hypothesis. For the target gene AT5G48970 : a. the EDI depending on α, b the MSE depending on α, c. the MSE depending on EDI (from panels a and b). The proposed gene-specific αt,opt is the value for which the MSE is most reduced as compared to the shuffled baseline, represented as a vertical blue line (Equations 10 and 11). d,e,f : the MSE depending on EDI for three other gene examples, representing different scenarios of data integration and thus different values of αt,opt. The trends are shown for weightedRF on true data (green) and shuffled datasets where TF expression profiles were randomly unmatched from their motif (grey). For each value of α, 100 models were run and the standard deviation around the mean is represented. The MSE is normalized by the variance of the target gene expression. Similar scenarios emerge in the linear model weightedLASSO (Figure S2). in which a target promoter harbors in average 23 TFBMs (Figure S4a) and a given TFBM is found in approximately 500 target promoters in average (Figure S4b).

Results

Optimal TFBM integration strength differs strongly between target genes

We ran weightedRF and weightedLASSO for values of α ranging from 0 to 1. We first observe that both weightedRF and weightedLASSO effectively incorporate TFBM information during their estimation, attributing higher importance measures to TFBM-supported variables as α increases. This is supported by EDI curves smoothly increasing with α (Figure S5). When applying a density threshold to build sparse GRN, we also observe that increasing α leads to the selection of edges with more and more TFBM support. At α = 1, TFBM support equals 1 meaning that, at the maximal level of data integration, GRN are restricted only to interactions supported by a TFBM (Figure S1).

An overview of the MSE profiles depending on α for all nitrate-responsive genes reveals a lot of diversity in how model performance can be driven by data integration strength, foreshadowing the usefulness of a gene-level procedure (Figure 4a). We thus applied DIOgene to optimize TFBM integration at the gene level (Equations 10 and 11). This confirmed that depending on the target genes, enforcing data integration has different effects on the predictive capabilities of the regression models, both in absolute error and relatively to the simulated null hypothesis. Very interestingly, for several genes like AT5G48970, increasing the EDI leads to a reduced MSE on test samples (Figure 3c). This illustrates that data integration can effectively guide the choice of variables toward more robust and meaningful regulators, allowing the model to better predict gene expression in unseen conditions. In this case, data integration can often be pushed to its maximal intensity, given that the maximal divergence from the simulated null data occur at αt,opt = 1. For several other genes, for example AT5G60670 (Figure 3d), the strongest improvement over the shuffled baseline is achieved for an intermediate value of α (0.5 in Figure 3d). For genes like AT3G20320, there is no reduction of MSE due to data integration, however DIOgene sets αt,opt to 1 because the MSE reduction in comparison to the shuffled baseline is sufficient (Figure 3e). Finally, the MSE of target genes can be increased by TFBM incorporation while showing no improvement over the simulated null data, like for instance AT1G30270, where αt,opt is set to 0 by our procedure (Figure 3f).

The application of DIOgene to all nitrate-responsive genes led to one αt,opt value per target gene. Among the 1,426 input genes, the number of genes for which TFBM information is integrated to expression, i.e αt,opt > 0, was 962 for weightedLASSO, and 760 for weightedRF. Overall, 572 genes have αt,opt > 0 in both weightedLASSO and weightedRF. This significant intersection (Figure S6) indicates that the two models mostly agree on a group of genes for which data integration is beneficial, even though specificities remain. The distribution of αt,opt for the 1,426 nitrate-responsive genes reveals that, similarly for the two models, a large pool of genes do not benefit from data integration according to our criterion (464 and 666 for weightedLASSO and weightedRF, respectively, Figure 4b). This suggests that data integration can often lead to a significant deterioration of model predictive capabilities as compared to our permuted control : in this case, our approach will leverage gene expression alone. These genes are further studied in the Discussion.

Properties of GRN inferred with gene-specific optimisation of data integration

Moving from the individual behavior of genes in response to data integration, we now measure global properties of GRN inferred using DIOgene. In order to evaluate the added value of tuning TFBM contributions in a gene-specific manner, we compare DIOgene GRN to GRN inferred with a parameter α identical for all genes, as done in previous approaches [START_REF] Greenfield | Robust data-driven incorporation of prior knowledge into the inference of dynamic regulatory networks[END_REF][START_REF] Petralia | Integrative random forest for gene regulatory network inference[END_REF][START_REF] Emily R Miraldi | Leveraging chromatin accessibility for transcriptional regulatory network inference in T helper 17 cells[END_REF][START_REF] Gibbs | High-performance single-cell gene regulatory network inference at scale: the inferelator 3.0[END_REF]. All GRN were built with a target density of 0.005, resulting in a total of 1,432 edges (Equation 12).

DIOgene provides a trade-off between accurate gene expression prediction and agreement with a DAP-Seq experiments

We rely on three metrics to assess inferred GRN qualities: the median MSE across the 1,426 nitrate-responsive genes together with GRN precision and recall against an experimental gold standard of in-vitro binding interactions (DAP-Seq) [START_REF] Ronan | Cistrome and epicistrome features shape the regulatory DNA landscape[END_REF].

First, the median MSE of GRN optimized with a global α displays a marked increase as the contribution of TFBM is reinforced (Figure 5a). This is in agreement with the previous observation that, for a majority of genes, TFBM deteriorate model predictions (Figure 4). Second, reinforcing TFBM contributions in GRN models equally for all genes increases both precision and recall against DAP-Seq interactions (Figures 5b, S7 and S8). Noteworthily, both models display a strong increase in precision with α, especially between α = 0 and α = 0.1, and weightedRF demonstrates a clear advantage over weightedLASSO, with a precision as high as 0.45. Both models outperform the precision of the prior PWM network for α > 0, indicating that using expression data to choose relevant links from all TFBM-supported interactions helps predicting actual TF binding, even in an invitro context. Recall values between the two models exhibit no clear differences. Thus, in globally optimised GRN, increasing data integration strength improves precision and recall, but necessarily comes with a deterioration of model predictions of the target gene expression (Figure 5a).

In contrast, gene-specific optimization of α with DIOgene offers a trade-off between these indicators. In fact, for both models, it provides a median MSE lower than any median MSE obtained with a global α (Figure 5a, blue). At the same time, GRN obtained with DIOgene achieve near-optimal precision and recall, as compared to global α curves. In this context, we actually argue that it is desirable to tolerate sub-obtimal precision and recall results while prioritizing low MSE, for three reasons: First, the MSE is specific to the conditions and cell lines used for transcriptome collection, which is not the case of DAP-Seq, that is agnostic of environmental and cellular context. Second, it is important to note that precision and recall also increase with α in shuffled datasets, where the wrong expression profile are randomly attributed to TFs (Figures S7 andS8). This illustrates that these statistics can be increased simply by enforcing data integration, even when gene expression data is uninformative. Finally, DAP-Seq and other external gold standards are often scarce and costly. Unlike the MSE, that can be measured for all genes, gold standards usually contain only a fraction of all TFs present in an inferred GRN, making it impossible to evaluate a large part of predicted interactions. We thus think that, in the context of this study and similar ones, precision and recall are unfit to properly tune, alone, the amount of a complementary omic source to incorporate into GRN inference.

DIOgene outperforms the minimal MSE approach

Because our TFBM integration approach is based on the MSE metric, we evaluated the benefit of optimising the MSE divergence from a shuffled baseline, over a simpler approach that would minimize the MSE directly.

First, the two approaches agree on setting α > 0 for a large group of genes (790 and 611 for weightedLASSO and weightedRF, respectively). They also both set α = 0 for 224 and 538 genes. Thus, our scheme and the minimal MSE approach perform data integration on globally similar sets of genes.

In contrast, some genes are set to α = 0 by our approach but not by the minimal MSE (240 and 128): these genes reach a minimal MSE for α > 0, but do not diverge sufficiently from the synthetic null hypothesis, and are thus removed from the data integration set by our approach (see two examples in Figure S9).

Other genes are set to α > 0 by our approach and to α = 0 by the minimal MSE (172 and 149). These genes typically display an increasing MSE that remains sufficiently lower than the shuffled control (see the example in Figure 3e). Our approach thus considers them for TFBM integration, tolerating that their MSE increases as long as it remains reasonably lower than the simulated control.

We then focused on the sets of genes for which we specifically integrate TFBM in one approach but not the other, and computed precision and recall curves of the corresponding sub-networks. The results (Figures S10c andS10d) showed better performance of our scheme, which suggests that the sets of genes considered for data integration on the basis of a comparison to a null hypothesis are more relevant than simply minimizing the MSE. Concordant results were obtained when computing the global qualities of inferred GRN. We observed a slight median MSE increase in our approach as compared to the minimal MSE approach (Figure S10a), while precision and recall against DAP-Seq are both globally improved in the two models by our integration scheme, especially recall (Figure S10b).

DIOgene improves the modelling of nitrate signalling

Finally, we assessed the ability of the inferred GRN to model nitrate induction pathways in Arabidopsis roots by comparing them to state of the art knowledge about this well documented response [START_REF] Bellegarde | Signals and players in the transcriptional regulation of root responses by local and systemic n signaling in arabidopsis thaliana[END_REF][START_REF] Vidal | Nitrate in 2020: Thirty years from transport to signaling networks[END_REF]. In order to identify the regulators predicted as important players in nitrate response by our models, we ranked regulators by out-degree in the inferred GRN. This was done for both weightedLASSO and weightedRF, either in GRN inferred with a global value of α = 0, α = 1, or with the proposed gene-specific optimisation of α (Figure S11).

A first observation is that, regardless of the chosen model or data integration strategy, the 25 TF with highest out-degree contain previously known master regulators of nitrate response. This includes DIV1 [START_REF] Cheng | Evolutionarily informed machine learning enhances the power of predictive gene-to-phenotype relationships[END_REF], TGA1 and TGA4 [START_REF] José | Systems approach identifies TGA1 and TGA4 transcription factors as important regulatory components of the nitrate response ofiarabidopsis thaliana/iroots[END_REF], as well as the homologs HHO2 and HHO3, belonging to the NIGT family and identified as repressing the expression of crucial nitrate transport genes [START_REF] Kiba | Repression of nitrogen starvation responses by members of the arabidopsis garp-type transcription factor nigt1/hrs1 subfamily[END_REF][START_REF] Safi | GARP transcription factors repress arabidopsis nitrogen starvation response via ROS-dependent and -independent pathways[END_REF]. Interestingly, we also uncover VRN1 and CRF4 as connectivity hubs in all inferred GRN. These regulators were respectively proposed as candidate and validated actors in nitrate signalling pathways in the studies that generated the transcriptomic data used here [START_REF] Varala | Temporal transcriptional logic of dynamic regulatory networks underlying nitrogen signaling and use in plants[END_REF][START_REF] Brooks | Network walking charts transcriptional dynamics of nitrogen signaling by integrating validated and predicted genome-wide interactions[END_REF]. Overall, whole-GRN measures of gene connectivity showed that genes involved in the regulation of nitrate pathways, nitrate uptake, transport and metabolism (Table S3) have a significantly higher total degree than other genes, in both globally optimized (at α = 0 and α = 1) and gene-specifically optimized GRN (Figure S12).

On another hand, we noticed that gene-specific calibration of data integration uniquely retrieves important regulators of nitrate nutrition, that were not present in the 25 most connected TF of the inferred GRN with a global α (α = 0 or α = 1). In the case of weightedLASSO, only the proposed gene-specific data integration strategy retrieves NLP7, which has been intensively 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0. 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0. .9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 documented as one of the main orchestrator of the early nitrate response [START_REF] Marchive | Nuclear retention of the transcription factor nlp7 orchestrates the early response to nitrate in plants[END_REF][START_REF] José | Transient genome-wide interactions of the master transcription factor NLP7 initiate a rapid nitrogen-response cascade[END_REF]. This is also the case of PHL1, a TF involved in the links between nitrate and phosphate signalling via NIGTmediated regulations [START_REF] Ueda | Nitrate-inducible NIGT1 proteins modulate phosphate uptake and starvation signalling via transcriptional regulation of iSPX/i genes[END_REF]. In the case of weightedRF, the proposed gene-specific optimization of data integration puts forward new TF as interesting candidates for nitrate response regulation. This includes HHO6, a member of the NIGT family not yet characterized for its role in the response to nitrate [START_REF] Kiba | Repression of nitrogen starvation responses by members of the arabidopsis garp-type transcription factor nigt1/hrs1 subfamily[END_REF][START_REF] Safi | GARP transcription factors repress arabidopsis nitrogen starvation response via ROS-dependent and -independent pathways[END_REF], but also BZIP53, a TF involved in the regulation of several facets of metabolism [START_REF] Garg | Targeted manipulation of bZIP53 DNA-binding properties influences arabidopsis metabolism and growth[END_REF] and JKD, a TF involved in shoot and root morphogenesis [START_REF] Kobayashi | INDETERMINATE DOMAIN PROTEIN binding sequences in the 5′-untranslated region and promoter of the SCARECROW gene play crucial and distinct roles in regulating SCARECROW expression in roots and leaves[END_REF]. Thus, this analysis reveals that this method of inference, via the optimization of data integration in a genespecific manner, not only recovers the information previously returned in the literature, but also brings to light new factors likely to be involved in this response.

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Discussion

The helpfulness of data integration is very often taken as granted in systems biology. Our work shows that it can in fact have very diverse effects on the modelling of gene expression, and that TFBM incorporation can be at the expense of model predictive capabilities for a significant number of genes. We thus propose to replace bulk data integration by a finely tuned hypothesisdriven data integration, calibrated individually for each gene. Our optimisation scheme, DIOgene, leverages TFBM in a way that their joint use with gene expression improves the target gene expression prediction over a simulated null hypothesis. On a plant biology case study, GRN inferred with this approach retrieve more DAP-Seq interactions through the use of TFBM and relevant known nitrate players, while preserving a near-optimal predictive performance on gene expression. Moreover, such conclusions hold for both the linear and non-linear regression cases, showing some general applicability to the most common models in the field. We also outlined some specificities in the tuning of TFBM integration between weightedLASSO and weightedRF (Figures S6,S11). Exploring these differences and the structure of the corresponding GRN would be a great way to test the impact of linearity and parametric assumptions in the modelling of multi-omics GRN.

The reason why some genes do not benefit from TFBM integration could stem from various factors, either technical or biological. Mining the consensual lists of genes for which α = 0 or α > 0 in both models revealed that genes for which TFBM are not integrated have in average a lower expression variance and more TFBM in their promoter (Figure S13). Further work would be needed to formulate hypotheses about the underlying general regulatory mechanisms, and also to assess the role of other forms of regulations like post-transcriptional and posttranslational modifications in these results.

Several limitations of this study should be reminded to the reader. First of all, as in all works inferring GRN from expression data alone, the expression of the regulators is taken as a proxy for their activity. This assumption is not always valid, which motivated the estimation of TF activities in other studies, typically leveraging motifs or binding experiments combined to gene expression [START_REF] Li | Regression analysis of combined gene expression regulation in acute myeloid leukemia[END_REF][START_REF] Mario | An experimentally supported model of the Bacillus subtilis global transcriptional regulatory network[END_REF][START_REF] Gibbs | High-performance single-cell gene regulatory network inference at scale: the inferelator 3.0[END_REF]. Our form of data integration, where TFBM-supported regulators have a stronger contribution in the estimated model, is another way to move away from this limitation. Even though this is a step toward more causality, challenges remain. One of them is the lack of a significant number of PWM, a problem amplified in non-model organisms. This limitation should be further reduced as PWM databases are completed and maintained by the community in the years to come, or as new computational methods are developed to predict binding affinities directly from DNA and protein sequences [START_REF] Barissi | DNAffinity: a machine-learning approach to predict DNA binding affinities of transcription factors[END_REF]. In the meantime, almost two thirds of regulators are attributed a neutral prior value, 1 2 , in our study. The choice for this neutral prior could however be different and reflect other modelling strategies that we did not explore. Another limitation of the TFBM prior is the high chance of false positives when scanning PWM: the presence of a binding motif can occur by chance, or may not necessarily cause binding nor regulation in a cellular context. In addition, TFBM with low complexity in PWM databases can result in hits in almost all the promoters of an organism (Figure S4b). In this analysis, regulators with such a widespread PWM were included, but their questionable biological relevance could lead to their exclusion or to additional weighting based on PWM complexity during GRN inference. Moreover, the PWM themselves may not be accurate priors for TF binding, due to experimental and technical biases in their identification. Finally, non canonical binding events driven by features like DNA shape, structure or repeat sequences without the need for a motif have been reported [START_REF] Md | Noncanonical binding of transcription factors: time to revisit ispecificity/i?[END_REF], but cannot be directly modelled in our approach. Another limitation is the computation time of the proposed optimization procedure. Because it requires the estimation of a large number of models over a complete range of α values, computation time could be reduced by several ways such as a re-implementation in C++, or through analyses estimating the minimal number of repetitions to properly assess the MSE and EDI curves (Equations 6, 7, 9). Finally, strong levels of correlation in the input data are still hindering accurate GRN inference. A lot of pairs of regulators have correlated expression profiles. When their TFBM are different, our integration scheme can identify the relevant regulator; however when their TFBM are very similar or both unknown, identifying the meaningful one is not guaranteed. Correlation between variables also impacts the design of simulated null datasets: in our study we randomly unmatched the expression profiles of the regulators. This has the advantage of preserving the correlation structure in the expression data, thus creating a realistic null dataset. However, when many regulators have similar expression profiles, which is, to some extent, true in our case study, the simulated null data may sometimes partly resemble the original data only by chance. Bringing more diverse expression profiles into the simulated datasets, or constraining it so that it does not contain similar versions of the true data could be envisioned.

In addition to the aforementioned perspectives, the application of the proposed data integration strategy to other complex organisms is a promising lead. In this work, TFBM influencing gene expression were assumed to be located in the promoter regions of their target genes because very few distal regulations have been reported in Arabidopsis, and are still poorly understood [START_REF] Lin | Enhancer-mediated reporter gene expression in iarabidopsis thaliana/i : a forward genetic screen[END_REF]. In organisms where regulation by distant enhancers is well documented and responsible for tissue-specificity [START_REF] Andersson | An atlas of active enhancers across human cell types and tissues[END_REF][START_REF] Kamal | scpGRaNIE/scp and scpGRaNPA/scp: inference and evaluation of enhancermediated gene regulatory networks[END_REF], delineating enhancer regions may be achieved through the use of additional molecular layers such as chromatin accessibility, chromatin contacts, or eQTLs. Enhancers and promoters could then be scanned for TFBM, linked to their target genes and further guide GRN inference.

In our case study, we favored the use of model prediction performance as a quality metric because it is a condition specific metric available for all genes and orthogonal to the integrated TFBM priors, which is often not the case of current experimental gold standards. Our results indicate that instead of directly minimizing prediction error as a function of TFBM contribution, the comparison to a shuffled baseline improved inferred GRN (Figure S10). In essence, any inference method where data integration is tuned by a parameter could be optimized based on such a simulated null dataset. As a general guideline, we believe that both the monitored quality metric and the simulated baseline should be carefully designed in order to test a clear and relevant hypothesis for the problem at hand. More generally, the concept of synthetic null datasets for in silico negative controls is gaining interest in genomic analyses. For example, scDEED [START_REF] Xia | scDEED: a statistical method for detecting dubious 2d single-cell embeddings[END_REF] optimizes 2D single-cell embeddings by simulating a null data where cells similarities are broken via permutations. Similarly, clusterDE [START_REF] Song | ClusterDE: a post-clustering differential expression (DE) method robust to false-positive inflation caused by double dipping[END_REF] reduces false discoveries in marker gene identification by generating a realistic null single-cell data via scDesign3 [START_REF] Song | scDesign3 generates realistic in silico data for multimodal single-cell and spatial omics[END_REF], and then contrasting pipeline outputs between true and null datasets. Such synthetic controls, bearing similarities with the methodology proposed in this article, are likely to enhance rigor and causal discoveries in the field. Fig. S10: Comparison of the proposed gene-specific data integration tuning approach, DIOgene (blue), to the simpler minimal MSE approach (red). a. Median MSE, computed over all nitrate-responsive genes. The MSE is normalized by the variance of the target gene expression. b. Precision and recall computed against DAP-Seq interactions, on two GRN densities : 0.005 (1432 edges), 0.01 (2864 edges). Globally, the proposed tuning approach provides GRN closer to the DAP-Seq network. Each dot is a replicate of GRN inference, randomness stemming from different bootstrap sampling. c-d. Precision and recall curves for densities ranging from 0.001 to 1, computed on sub-GRN made of the edges concerning only target genes for which we integrate TFBM exclusively in one of the compared methods. In red, the genes for which the minimal MSE sets αt,opt > 0 but not the proposed approach. In blue, the genes for which our approach sets αt,opt > 0 but not the minimal MSE approach. This highlights the precision and recall gain from changing the sets of genes for which TFBM are integrated between the two methods. In inferred GRN of a 0.005 density, regulators were ranked by their out-degree (number of regulated target genes or other regulators) for weightedLASSO and weightedRF, and for a global value of α (0 or 1), or using the proposed gene-specific approach (DIOgene). The out degree is averaged across the different runs of the same method to account for inherent stochasticity, and the top 25 TFs with the highest out-degree are shown. Regulators already identified in previous studies as important (orange) or candidate (yellow) actors of nitrate signalling and regulation are highlighted [START_REF] Bellegarde | Signals and players in the transcriptional regulation of root responses by local and systemic n signaling in arabidopsis thaliana[END_REF][START_REF] Vidal | Nitrate in 2020: Thirty years from transport to signaling networks[END_REF][START_REF] Cheng | Evolutionarily informed machine learning enhances the power of predictive gene-to-phenotype relationships[END_REF][START_REF] José | Systems approach identifies TGA1 and TGA4 transcription factors as important regulatory components of the nitrate response ofiarabidopsis thaliana/iroots[END_REF][START_REF] Kiba | Repression of nitrogen starvation responses by members of the arabidopsis garp-type transcription factor nigt1/hrs1 subfamily[END_REF][START_REF] Safi | GARP transcription factors repress arabidopsis nitrogen starvation response via ROS-dependent and -independent pathways[END_REF][START_REF] Varala | Temporal transcriptional logic of dynamic regulatory networks underlying nitrogen signaling and use in plants[END_REF][START_REF] Brooks | Network walking charts transcriptional dynamics of nitrogen signaling by integrating validated and predicted genome-wide interactions[END_REF][START_REF] Marchive | Nuclear retention of the transcription factor nlp7 orchestrates the early response to nitrate in plants[END_REF][START_REF] José | Transient genome-wide interactions of the master transcription factor NLP7 initiate a rapid nitrogen-response cascade[END_REF][START_REF] Ueda | Nitrate-inducible NIGT1 proteins modulate phosphate uptake and starvation signalling via transcriptional regulation of iSPX/i genes[END_REF][START_REF] Safi | GARP transcription factors repress arabidopsis nitrogen starvation response via ROS-dependent and -independent pathways[END_REF][START_REF] Garg | Targeted manipulation of bZIP53 DNA-binding properties influences arabidopsis metabolism and growth[END_REF][START_REF] Kobayashi | INDETERMINATE DOMAIN PROTEIN binding sequences in the 5′-untranslated region and promoter of the SCARECROW gene play crucial and distinct roles in regulating SCARECROW expression in roots and leaves[END_REF]. TFs uniquely retrieved by the proposed gene-specific approach as compared to a global optimization for a given model are reported in blue. The overall degree (in-degree and out-degree) of nitrate-related genes is higher in inferred GRN as compared to other genes. In inferred GRN (D = 0.005), the total degree of genes is shown for a global value of α (0 or 1), or using DIOgene (specific). Total degree is reported on the log scale. Only genes with at least one connection in the inferred GRN were considered. The list of nitrate-related genes (Table S3) was compiled from the literature [START_REF] Bellegarde | Signals and players in the transcriptional regulation of root responses by local and systemic n signaling in arabidopsis thaliana[END_REF][START_REF] Vidal | Nitrate in 2020: Thirty years from transport to signaling networks[END_REF][START_REF] Cheng | Evolutionarily informed machine learning enhances the power of predictive gene-to-phenotype relationships[END_REF][START_REF] José | Systems approach identifies TGA1 and TGA4 transcription factors as important regulatory components of the nitrate response ofiarabidopsis thaliana/iroots[END_REF][START_REF] Kiba | Repression of nitrogen starvation responses by members of the arabidopsis garp-type transcription factor nigt1/hrs1 subfamily[END_REF][START_REF] Safi | GARP transcription factors repress arabidopsis nitrogen starvation response via ROS-dependent and -independent pathways[END_REF][START_REF] Varala | Temporal transcriptional logic of dynamic regulatory networks underlying nitrogen signaling and use in plants[END_REF][START_REF] Brooks | Network walking charts transcriptional dynamics of nitrogen signaling by integrating validated and predicted genome-wide interactions[END_REF][START_REF] Marchive | Nuclear retention of the transcription factor nlp7 orchestrates the early response to nitrate in plants[END_REF][START_REF] José | Transient genome-wide interactions of the master transcription factor NLP7 initiate a rapid nitrogen-response cascade[END_REF][START_REF] Ueda | Nitrate-inducible NIGT1 proteins modulate phosphate uptake and starvation signalling via transcriptional regulation of iSPX/i genes[END_REF][START_REF] Safi | GARP transcription factors repress arabidopsis nitrogen starvation response via ROS-dependent and -independent pathways[END_REF][START_REF] Garg | Targeted manipulation of bZIP53 DNA-binding properties influences arabidopsis metabolism and growth[END_REF][START_REF] Kobayashi | INDETERMINATE DOMAIN PROTEIN binding sequences in the 5′-untranslated region and promoter of the SCARECROW gene play crucial and distinct roles in regulating SCARECROW expression in roots and leaves[END_REF].

Fig. 1 :

 1 Fig.1: Overview of weightedRF and weightedLASSO integrative GRN inference procedures. Input data for a given target gene is the vector of expression of this gene Y , an expression matrix for the regulators X, and a TFBM scoring matrix Π. Π contains information about the presence or absence of the regulators PWM in target gene's promoter. α is a parameter controlling the force of TFBM integration to expression. For each target gene, a regression model is fit to Y with weightedRF or weightedLASSO using X as predictive variables, and favoring the attribution of high influences to regulators with high values in Π. For weightedRF, this prioritization is executed by a weighted subsampling of the regulators space when elongating the regression trees, while it is achieved by differential shrinkage combined to stability selection is weightedLASSO. Once all regulator-gene pairs have been ranked based on their influence in the regression models, final GRN are built by selecting the number of strongest interactions providing a desired network density.

Fig. 4 :

 4 Fig.4: Gene-specific tuning of TFBM integration in the 1,426 nitrate-responsive genes with DIOgene leads to diverse MSE behaviors and integration intensities. a. Scaled MSE (z-score) in weightedLASSO and weightedRF on true data depending on α for two types of genes: genes with αopt = 0 (grey) and genes with αopt > 0 (green). b. Distribution of αopt values for the 1,426 nitrate-responsive genes in weightedLASSO and weightedRF.

0. 1 0

 1 .1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.

  4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4

Fig. 5 :

 5 Fig.5: Gene-specific integration of TFBM to gene expression with DIOgene optimises model performance and secures a satisfying intersection with binding experiments. For increasing values of α (green) and for the proposed gene-specific optimisation of α (blue), we show the median MSE (normalized by the variance of the target gene) of the nitrate-responsive genes (a) and precision as a function of recall in the inferred GRN (1,432 edges, density = 0.005) against DAP-Seq interactions (b). The precision of the prior TFBM network of nitrate-responsive genes (31956 edges, density = 0.32) is overlaid in orange. Its recall is not shown because fairly comparing the recall of GRN requires that they have the same number of edges.

Fig. S2 :

 S2 Fig.S1: TFBM support of inferred GRN. TFBM support is the average value of the inferred TF-target edges in the prior TFBM matrix Π. It is shown for three network densities : 0.005 (1432 edges), 0.01 (2864 edges) and 0.05 (14322 edges). A TFBM of 1 means that a GRN is composed only of edges supported by a TFBM. At maximal TFBM integration intensity (α = 1), GRN on both true and shuffled data achieve a TFBM support of 1.

Fig. S3 :

 S3 Fig. S3: PCA plot of the normalized counts from the RNA-Seq experiment of the dynamic root response to a nitrate treatment [59]. The three first subplots show the correlation of each sample to the first 4 principal components. Each color represents an experimental condition made of 3 replicates. C: control. N: nitrate treatment. Numbers represent time in minutes. The first component highlights the samples at 10 and 15 minutes as different from the rest of the experiment, possibly due to a batch effect. The screeplot represents the percentage of variance explained by the 4 first principal components.

  Fig.S9: Example of two genes in weightedRF for which DIOgene sets αt,opt = 0, but the minimal MSE approach would set αt,opt > 0. MSE depending on EDI for all possible values of α on true data (green) and shuffled null datasets (grey).The MSE is normalized by the variance of the target gene expression.

Fig. S11 :

 S11 Fig.S11: Most connected regulators depending on the model and TFBM integration strategy. In inferred GRN of a 0.005 density, regulators were ranked by their out-degree (number of regulated target genes or other regulators) for weightedLASSO and weightedRF, and for a global value of α (0 or 1), or using the proposed gene-specific approach (DIOgene). The out degree is averaged across the different runs of the same method to account for inherent stochasticity, and the top 25 TFs with the highest out-degree are shown. Regulators already identified in previous studies as important (orange) or candidate (yellow) actors of nitrate signalling and regulation are highlighted[START_REF] Bellegarde | Signals and players in the transcriptional regulation of root responses by local and systemic n signaling in arabidopsis thaliana[END_REF][START_REF] Vidal | Nitrate in 2020: Thirty years from transport to signaling networks[END_REF][START_REF] Cheng | Evolutionarily informed machine learning enhances the power of predictive gene-to-phenotype relationships[END_REF][START_REF] José | Systems approach identifies TGA1 and TGA4 transcription factors as important regulatory components of the nitrate response ofiarabidopsis thaliana/iroots[END_REF][START_REF] Kiba | Repression of nitrogen starvation responses by members of the arabidopsis garp-type transcription factor nigt1/hrs1 subfamily[END_REF][START_REF] Safi | GARP transcription factors repress arabidopsis nitrogen starvation response via ROS-dependent and -independent pathways[END_REF][START_REF] Varala | Temporal transcriptional logic of dynamic regulatory networks underlying nitrogen signaling and use in plants[END_REF][START_REF] Brooks | Network walking charts transcriptional dynamics of nitrogen signaling by integrating validated and predicted genome-wide interactions[END_REF][START_REF] Marchive | Nuclear retention of the transcription factor nlp7 orchestrates the early response to nitrate in plants[END_REF][START_REF] José | Transient genome-wide interactions of the master transcription factor NLP7 initiate a rapid nitrogen-response cascade[END_REF][START_REF] Ueda | Nitrate-inducible NIGT1 proteins modulate phosphate uptake and starvation signalling via transcriptional regulation of iSPX/i genes[END_REF][START_REF] Safi | GARP transcription factors repress arabidopsis nitrogen starvation response via ROS-dependent and -independent pathways[END_REF][START_REF] Garg | Targeted manipulation of bZIP53 DNA-binding properties influences arabidopsis metabolism and growth[END_REF][START_REF] Kobayashi | INDETERMINATE DOMAIN PROTEIN binding sequences in the 5′-untranslated region and promoter of the SCARECROW gene play crucial and distinct roles in regulating SCARECROW expression in roots and leaves[END_REF]. TFs uniquely retrieved by the proposed gene-specific approach as compared to a global optimization for a given model are reported in blue.

  Fig.S12: The overall degree (in-degree and out-degree) of nitrate-related genes is higher in inferred GRN as compared to other genes. In inferred GRN (D = 0.005), the total degree of genes is shown for a global value of α (0 or 1), or using DIOgene (specific). Total degree is reported on the log scale. Only genes with at least one connection in the inferred GRN were considered. The list of nitrate-related genes (TableS3) was compiled from the literature[START_REF] Bellegarde | Signals and players in the transcriptional regulation of root responses by local and systemic n signaling in arabidopsis thaliana[END_REF][START_REF] Vidal | Nitrate in 2020: Thirty years from transport to signaling networks[END_REF][START_REF] Cheng | Evolutionarily informed machine learning enhances the power of predictive gene-to-phenotype relationships[END_REF][START_REF] José | Systems approach identifies TGA1 and TGA4 transcription factors as important regulatory components of the nitrate response ofiarabidopsis thaliana/iroots[END_REF][START_REF] Kiba | Repression of nitrogen starvation responses by members of the arabidopsis garp-type transcription factor nigt1/hrs1 subfamily[END_REF][START_REF] Safi | GARP transcription factors repress arabidopsis nitrogen starvation response via ROS-dependent and -independent pathways[END_REF][START_REF] Varala | Temporal transcriptional logic of dynamic regulatory networks underlying nitrogen signaling and use in plants[END_REF][START_REF] Brooks | Network walking charts transcriptional dynamics of nitrogen signaling by integrating validated and predicted genome-wide interactions[END_REF][START_REF] Marchive | Nuclear retention of the transcription factor nlp7 orchestrates the early response to nitrate in plants[END_REF][START_REF] José | Transient genome-wide interactions of the master transcription factor NLP7 initiate a rapid nitrogen-response cascade[END_REF][START_REF] Ueda | Nitrate-inducible NIGT1 proteins modulate phosphate uptake and starvation signalling via transcriptional regulation of iSPX/i genes[END_REF][START_REF] Safi | GARP transcription factors repress arabidopsis nitrogen starvation response via ROS-dependent and -independent pathways[END_REF][START_REF] Garg | Targeted manipulation of bZIP53 DNA-binding properties influences arabidopsis metabolism and growth[END_REF][START_REF] Kobayashi | INDETERMINATE DOMAIN PROTEIN binding sequences in the 5′-untranslated region and promoter of the SCARECROW gene play crucial and distinct roles in regulating SCARECROW expression in roots and leaves[END_REF].

  Intersection between the genes with αt,opt > 0 in weightedLASSO and weightedRF with DIOgene. Enrichment pvalue : 5.9e -12 (One-sided hypergoemetric test). Precision and recall as a function of α in weightedRF, against a DAP-Seq gold standard. Precision and recall are show on true data (green) and shuffled null datasets (grey) for three network densities : 0.005 (1432 edges), 0.01 (2864 edges) and 0.05 (14322 edges). The precision of the prior TFBM network of nitrate responsive genes (31956 edges, density = 0.32) is overlaid in orange.

										weightedLASSO weightedRF
									Precision against DAP-Seq Precision against DAP-Seq
								weightedRF dataset dataset	shuffled shuffled	weightedLASSO trueData trueData
	Precision	0.00 0.20 0.25 0.30 0.35 0.40 0.45 Fig. S6: D = 0.005 0.25 0.50 0.20 0.25 0.30 0.35 0.40 D = 0.005 0.75 0.00 0.25 0.50 Precision	188 (16%) 1.00 0.00 0.75 1.00 0.00 0.25 Recall against DAPSeq 390 (34%) 572 (50%) D = 0.01 0.50 0.75 1.00 0.00 α D = 0.01 0.25 0.50 0.75 1.00 0.00 0.25 α Recall against DAPSeq	D = 0.05 0.50 D = 0.05 0.75 0.25 0.50	1.00 0.75	1.00
									dataset dataset	shuffled shuffled	trueData trueData
				D = 0.005 D = 0.005					D = 0.01 D = 0.01			D = 0.05 D = 0.05
									0.05	0.05
															0.20	0.20
		0.02	0.02						0.04	0.04
	Recall	Recall							0.03	0.03					0.15	0.15
		0.01	0.01						0.02	0.02					0.10	0.10
									0.01	0.01					0.05	0.05
		0.00	0.25 0.00	0.50 0.25	0.75 0.50	1.00 0.75	0.00 1.00	0.25 0.00	0.50 α 0.25	0.75 0.50 α	1.00 0.75	0.00 1.00	0.25 0.00	0.50 0.25	0.75 0.50	1.00 0.75	1.00
	Fig. S7: Precision and recall as a function of α in weightedLASSO, against a DAP-Seq gold standard. Precision and recall Fig. S8:
	are show on true data (green) and shuffled null datasets (grey) for three network densities : 0.005 (1432 edges), 0.01 (2864 edges) and
	0.05 (14322 edges). The precision of the prior TFBM network of nitrate responsive genes (31956 edges, density = 0.32) is overlaid in
	orange.

Although samples at 10 and 15 min were also available, we discarded them, as they were clear outliers in the first dimension of a Principal Component Analysis (PCA) (FigureS3).

Acknowledgments

This research is supported by a CNRS 80 Prime grant (BREAK), the French National Research Agency, and by the LabMUSE EpiGenMed (R-loops project). We thank Anaïs Baudot for her comment on an oral presentation of this work, which inspired a part of the discussion. We are grateful to all members of the Computational Regulatory Genomics team, Montpellier, for discussions around this project.

Data and code availability

Below are the links to the data used in the course of this study.

• The RNA-Seq data for the response to nitrate induction was downloaded from the GEO accession GSE97500 • The PWM used to build the TFBM dataset were retrieved from JASPAR and the Plant Cistrome Database. • To identify Arabidopsis TSSs and promoter regions, we relied on the TAIR10 GFF3 file. • The regulators of Arabidopsis used for GRN inference are the union between PlnTFDB and AtTFDB All results can be reproduced with the code available in the github repository: https://github.com/OceaneCsn/integrative_GRN_N_induction

Competing interests

The authors declare no competing interests.

Author contributions statement SL, LB, CHL and OC conceived the statistical inference methods and optimization scheme. OC implemented the inference methods and the optimization scheme. AM, SL and OC compared the inferred GRN to state of the art knowledge about nitrate signalling. All authors read and took part in writing the manuscript.

Supplementary figures

Table S1 : Normalized gene expression of the 1426 nitrateresponsive genes in the different treatments and time points of the experiment [START_REF] Varala | Temporal transcriptional logic of dynamic regulatory networks underlying nitrogen signaling and use in plants[END_REF]. C: control. N: nitrate induction treatment. Numbers represent time after treatment in minutes.

Table S2 : Gene identifier (AGI) of the 201 nitrate-responsive regulators.

Table S3 : 56 genes involved in the uptake, transport, metabolism and signalling (positive or negative) of nitrate in the roots of Arabidopsis thaliana. Their AGI, gene name and short description are shown. This list of nitrate-related genes was compiled from the literature [START_REF] Bellegarde | Signals and players in the transcriptional regulation of root responses by local and systemic n signaling in arabidopsis thaliana[END_REF][START_REF] Vidal | Nitrate in 2020: Thirty years from transport to signaling networks[END_REF][START_REF] Cheng | Evolutionarily informed machine learning enhances the power of predictive gene-to-phenotype relationships[END_REF][START_REF] José | Systems approach identifies TGA1 and TGA4 transcription factors as important regulatory components of the nitrate response ofiarabidopsis thaliana/iroots[END_REF][START_REF] Kiba | Repression of nitrogen starvation responses by members of the arabidopsis garp-type transcription factor nigt1/hrs1 subfamily[END_REF][START_REF] Safi | GARP transcription factors repress arabidopsis nitrogen starvation response via ROS-dependent and -independent pathways[END_REF][START_REF] Varala | Temporal transcriptional logic of dynamic regulatory networks underlying nitrogen signaling and use in plants[END_REF][START_REF] Brooks | Network walking charts transcriptional dynamics of nitrogen signaling by integrating validated and predicted genome-wide interactions[END_REF][START_REF] Marchive | Nuclear retention of the transcription factor nlp7 orchestrates the early response to nitrate in plants[END_REF][START_REF] José | Transient genome-wide interactions of the master transcription factor NLP7 initiate a rapid nitrogen-response cascade[END_REF][START_REF] Ueda | Nitrate-inducible NIGT1 proteins modulate phosphate uptake and starvation signalling via transcriptional regulation of iSPX/i genes[END_REF][START_REF] Safi | GARP transcription factors repress arabidopsis nitrogen starvation response via ROS-dependent and -independent pathways[END_REF][START_REF] Garg | Targeted manipulation of bZIP53 DNA-binding properties influences arabidopsis metabolism and growth[END_REF][START_REF] Kobayashi | INDETERMINATE DOMAIN PROTEIN binding sequences in the 5′-untranslated region and promoter of the SCARECROW gene play crucial and distinct roles in regulating SCARECROW expression in roots and leaves[END_REF]. or not integrated (αt,opt = 0, right, N = 276) in weightedRF and weightedLASSO. Expression and expression variance are reported on the log scale. The other characteristics are the number of motifs in the target gene promoter region, either among all knwon PWM in Arabidopsis (nb motifs), or only in nitrate responsive regulators (nb motifs n tfs), the number of introns (nb introns), the number of transcripts (nb transcripts). The two latter features were retrieved from the TAIR10 GFF annotation. On the last row, proportions of TFs in each group of genes are shown along an enrichment p-value result (one-sided hypergeometric test).


M Hewitt
email: mhewitt3@luc.edu

F Lehuédé

The Service Network Scheduling Problem

Keywords: Service Network Design Problem, Scheduling, Column Generation, Integer Programming

We consider the optimization problem of determining schedules for shipments on known paths within a terminal network in order to minimize vehicle transportation costs. We refer to this problem as the Service Network Scheduling Problem and present two mixed integer programming formulations of that problem. The first is based on the classical idea of a time-expanded network. The second formulation is new and is based on sets of shipment consolidations. We show both analytically and computationally that the consolidation-based formulation can be the superior of the two, but that its enumerative nature renders it ineffective for instances with large numbers of shipments. Thus, we present a column generation-based algorithm for solving the consolidation-based formulation that relies on solving relaxations that are integer programs.

We demonstrate the superior performance of this algorithm with a computational study.

INTRODUCTION

Consolidation carriers are transportation companies that transport shipments that are small relative to vehicle capacity. Consolidation carrier is an umbrella term that covers companies participating in one (or both) of the less-than-truckload freight and small parcel industries. As transportation exhibits economies of scale, profitability for such companies is driven by consolidation. Specifically, dispatching vehicles that transport multiple shipments, each potentially associated with a different customer. Such consolidation is typically enabled by routing shipments on paths through a network of terminals as opposed to directly from customer origin to destination, as is often done in full truckload transportation.

More precisely, a path for a shipment refers to the sequence of terminals it visits, with the sequence beginning at the origin terminal for the shipment and ending at its destination terminal. We also refer to these paths as consisting of one or more transportation moves, where a transportation move refers to physical transportation, either by a shipment or a capacitated vehicle, between terminals in the network. Transportation moves executed by a vehicle incur a cost that is independent of the shipments it transports. Lastly, the schedule for a path prescribes dispatch times for each of its transportation moves. These schedules in turn imply vehicle dispatch times on those same moves. Achieving high levels of consolidation, and low transportation costs, requires determining paths for shipments through such networks, and schedules for those paths, that enable multiple shipments to dispatch on the same transportation move at the same time. Jointly determining these paths and schedules is typically seen as a tactical planning problem and is modeled as a variant of the Scheduled Service Network Design Problem (SSNDP) [START_REF] Crainic | Network Design with Applications to Transportation and Logistics[END_REF].

A carrier may seek to keep shipment paths constant over some planning horizon (e.g. a month), in order to maintain consistency with the operations that support the execution of those paths. However, for some types of carriers (e.g. road-based Less-than-truckload freight transportation carriers) there is often flexibility, operationally-speaking, with respect to the scheduling of paths. In the near-term when more accurate forecasts of shipment volumes are known, it may be willing to adjust the schedules for those paths if doing so leads to greater consolidation. Thus, we consider the problem of optimally determining shipment path schedules in order to minimize total vehicle transportation costs. We refer to this problem as the Service Network Scheduling Problem (SNSP).

We propose two mixed integer programming formulations of this problem. The first is based on the classical idea of a time-expanded network [START_REF] Ford | Constructing maximal dynamic flows from static flows[END_REF][START_REF] Ford | Flows in networks[END_REF]. The second is based on enumerations of consolidations of shipments. We prove the equivalence of the formulations. We also show both analytically and computationally that the consolidation-based formulation can be the stronger of the two, but that its enumerative nature renders it computationally ineffective for instances based on large numbers of shipments.

Thus, we present an algorithm for solving the consolidation-based formulation that does not enumerate sets of consolidations a priori, but instead dynamically in the course of its execution. Clearly, one framework for such an algorithm is Branch-and-Price [START_REF] Barnhart | Branch-and-price: Column generation for solving huge integer programs[END_REF][START_REF] Desaulniers | Column generation[END_REF]. As such an algorithm involves solving linear relaxations of a formulation, many of the technical advancements of today's solvers at solving integer programs is lost. We present an algorithm that operates in a fashion that is similar to Branch-and-Price, in that at each iteration it solves an optimization problem and then uses information from the solution to that optimization problem to determine variables to add. It differs from Branch-and-Price in that the optimization problem solved is an integer program that is formulated so as to be a relaxation of the original problem. We refer to this algorithm as IP -ColGen, or, Integer Programmingbased Column Generation. We prove the correctness of this algorithm and demonstrate with a computational study its superiority to solving the static formulation. We believe this algorithm and formulation technique can be employed in formulations of and solution methods for the more general SSNDP, something we will discuss during the presentation.

Methodology

In this section, we present the consolidation-based methodology for this problem. We omit the time-expanded network formulation for brevity and because it is classical. We model the terminal network and moves within that network with the directed network D = (N , A), where N is a set of nodes that model terminals in the network and A is a set of arcs that model transportation moves within that network. Formally, the set A consists of pairs (i, j), i, j ∈ N that model physical travel from terminal i to terminal j. Associated with each arc (i, j) ∈ A is a travel time denoted by τ ij , a per-vehicle capacity denoted by u ij , and a per-vehicle cost denoted by f ij .

We model shipments that must be routed through the network with the set K of commodities. Associated with each commodity k ∈ K is an origin terminal, o k , and destination terminal, d k .

We let e k denote the earliest time at which commodity k is available at its origin terminal, o k , and l k denote the latest time at which k can be delivered to its destination terminal, d k . Lastly, associated with commodity k is a size q k quoted in the same unit as the capacity factors u ij , (i, j) ∈ A associated with arcs.

Recalling that our problem presumes that paths for shipments are known, we let

p k = v k 1 , ..., v k r k represent the sequence of nodes from N in the path of commodity k ∈ K, with v k 1 = o k and v k r k = d k .
We let P = ∪ k∈K p k denote the set of all such shipment paths. Relatedly, for commodity k we let the node set N k ⊆ N contain the nodes v k i in that path and the arc set

A k ⊆ A contain the arcs (v k i , v k i+1
) in that path. We also let K ij = {k ∈ K : (i, j) ∈ A k } denote the set of commodities with a path that contains arc (i, j) ∈ A. Lastly, we note that given the path p k and earliest available and latest due times e k , l k for commodity k, one can derive a time window [α k v , β k v ] during which k can be at each node v in its path. The proposed formulation of the SNSP that is based on sets of possible consolidations on each arc, (i, j) ∈ A. To define the formulation we let kq i ] = ∅. In words, a consolidation is feasible on an arc if and only if the time windows for all commodities in that consolidation overlap in at least one time point. That said, C ij also contains all singleton sets. Namely, if k ∈ K ij , then there exists g ∈ [1, n ij ] such that C g = {k}. In addition, we let C = ∪ (i,j)∈A C ij . Throughout this paper we will refer to "consolidation C g " as shorthand for the set of commodities {k g 1 , . . . , k g mg } contained in that consolidation. Regarding data elements associated with these consolidation sets, we let the attribute

C ij = {C 1 , . . . , C n ij }, C g ⊆ K ij ∀g = 1, . . . ,
φ k C ∈ {0, 1}, k ∈ K, C ∈ C represent whether k ∈ C. We also let s C = k∈C q k u ij
represent the number of vehicles needed to transport consolidation C.

To model the SNSP, we let the binary variable w C indicate whether the consolidation consisting of commodities in set C ∈ C ij traveling together on arc (i, j) ∈ A is chosen. Note this choice implies that the commodities must dispatch at the same time. We let the integer variable y ij represent the number of vehicles that dispatch on arc (i, j) ∈ A. We let the decision variables

γ k v k i v k i+1
, i = 1, . . . , n k -1 prescribe the time at which commodity k dispatches on the arc

(v k i , v k i+1
), i = 1, . . . , n k -1 on its path. With these sets and decision variables, we define the optimization problem Cons-SNSP(C), as

z Cons (C) = minimize (i,j)∈A f ij y ij (1) subject to C∈C ij φ k C w C = 1 ∀k ∈ K, (i, j) ∈ A k , (2)
C∈C ij s C w C ≤ y ij ∀(i, j) ∈ A, (3)
γ k ij -γ k ij ≤ M kk i (1 - C∈C ij φ k C φ k C w C ) ∀(i, j) ∈ A, k, k ∈ K ij , (4)
γ k v k i v k i+1 + τ v k i v k i+1 ≤ γ v k i+1 v k i+2 ∀k ∈ K, i = 1, . . . , n k -2, (5)
α v k i ≤ γ k v k i v k i+1 ≤ β v k i ∀k ∈ K, i = 1, . . . , n k -1, (6)
γ k v k i v k i+1 ∈ N ∀(v k i , v k i+1 ) ∈ p k , k ∈ K, (7)
w C ∈ {0, 1} ∀C ∈ C, (8)
y ij ∈ N ∀(i, j) ∈ A. (9)
The objective seeks to minimize the total costs associated with vehicle moves that transport consolidations of shipments. Constraints (2) ensure that a consolidation is chosen for each commodity on each arc in its path. Constraints (3) ensure that sufficient capacity is paid for on each arc to support the consolidations chosen for that arc. Constraints (4) ensure that all commodities in the consolidation chosen for an arc dispatch at the same time. Constraints (5) ensure that the dispatch times for arcs on the path of a commodity agree with their travel times. Constraints (6) ensure the dispatch time decision variables occur within the corresponding time windows. Constraints ( 7), (8), and (9) define the decision variables and their domains. While the big-M value in constraints (4) may lead to a weak formulation, the formulation can be tightened without rendering any solutions infeasible by setting

M kk i = β k i -α k i .

Results

To compare the two formulations computationally, we generated a set of instances based on a portion of the network of a United States-based LTL carrier. Specifically, we considered a portion of the network that consists of 25 terminals (e.g. |N | = 25) and 530 physical moves between terminals (e.g. |A| = 530). Cost, capacity, and travel time data were provided by the carrier. In addition, the carrier provided a load plan that prescribed paths through the terminal network for pairs of terminals (o, d) ∈ N based on their customer base at the time. We randomly generated instances that vary in the number of commodities (100,150,200,250,300,350,400,450, and 500) and other instance parameters.

We compare solving the consolidation-based formulation Cons-SNSP(C) with solving three variants of the time-expanded network formulation. The first, labeled "T-E Network" consists of solving the time-expanded network formulation over a complete time-expanded network. The second, labeled "Reduced T-E Network" consists of solving the formulation over a reduced timeexpanded network. The third, labeled "Reduced T-E Network + Valid inequalities" consists of solving the formulation over a reduced time-expanded network but strengthened with valid inequalities. However, a downside of the consolidation-based formulation is its enumerative nature. Thus, we have also developed a column generation-based algorithm to solve this formulation wherein consolidations are generated dynamically during the course of its execution. We will present this algorithm in detail during the talk. However, in Table 2 we report results from executing this algorithm on the same instances discussed above. We see that the algorithm dramatially

  n ij , denote all sets of commodities that can dispatch on arc (i, j) at the same time. Recalling that[α k i , β k i ]denotes the time window during which commodity k must depart from node i, we have that C g = {k g 1 , . . . , k g mg } ∈ C ij ∀g = 1, . . . , n ij if and only if ∩

	mg q=1 [α i , β kq

Table 1 -

 1 % solved, time to termination, averaged over all instances

	Method	% Solved Time (sec.)
	T-E Network	44.44%	4,094.54
	Reduced T-E network	58.33%	3,080.16
	Reduced T-E Network + Valid inequalities	66.67%	2,531.47
	Consolidation	77.78%	1,770.94

Table 2 -

 2 Comparison of IP -ColGen with Cons-SNSP(C) by number of commodities

	# Commodities	100	150	200	250	300	350	400	450	500	Average
	IP-ColGen
	% Solved	100%	100%	100%	100%	100%	100%	100%	100%	100%	100%
	Time (sec.)	0.52	1.10	1.58	2.79	6.52	8.16	28.56	75.16	131.18	28.40
	Cons-SNSP(C)
	% Solved	100.00%	100.00%	100.00%	100.00%	75.00%	75.00%	100.00%	25.00%	25.00%	77.78%
	Time (sec.)	0.52	1.41	3.00	5.43	1,807.23	1,857.88	606.80	5,462.01	6,194.21	1,770.94

outperforms solving the static, a priori formulation.