S Mouton
email: serge.mouton@u-bordeaux1.fr

D Teissandier
email: denis.teissandier@u-bordeaux1.fr

P Sébastian

J P Nadeau

MANUFACTURING REQUIREMENTS IN DESIGN: THE RTM PROCESS IN AERONAUTICS

Keywords: Resin Transfer Molding (RTM), assembly, residual stresses, industrialisation, strain, geometric specification for non-rigid parts, Constraint Satisfaction Problem (CSP) solver

A sub-unit of an aeronautical structure (fuselage, fin, wing, etc.) consists of a set of components fixed rigidly together. One of today's major industrial challenges is to produce these sub-units out of composite materials in order to increase the level of integration and reduce mass and cost. This article describes a procedure to assist in the industrialisation of aeronautical components produced from composite materials in a design for manufacturing context. In a multi-expertise approach, the problem of optimising integration is combined with the feasibility of injection for the Resin Transfer Molding process. This approach then takes into account admissible manufacturing deviations, defined from a classification of the structure parts. The limits set for admissible deviations guarantee the mechanical behaviour of the assembled component and the requirements of the assembly as a whole. Finally, an industrialisation solutions space is defined. A constraint satisfaction problem solver is used to carry out this research with a spar from a horizontal plane in an aircraft used to illustrate the procedure.

INTRODUCTION

When taking industrialisation into account in the design cycle [START_REF] Herrmann | New Directions in Design for Manufacturing[END_REF][START_REF] Kim | Weight Minimization of Composite Laminated Plates with Multiple Constraints[END_REF][START_REF] Nardari | Simultaneous Engineering in Design and Manufacture using RTM Process[END_REF] two distinct and complementary aspects must be considered: technical feasibility and economic feasibility. The technical feasibility of the RTM process imposes specific shapes and dimensions. Moreover, the RTM process can cause strain as parts are removed from the moulds, caused by the architecture of the laminate and the great complexity of physical phenomena around the parameters governing the RTM process. These strains are characterised by manufacturing deviations from which the true geometry of the part can be defined in relation to its design geometry (CAD model). In order FPCM-9 (2008) The 9 th International Conference on Flow Processes in Composite Materials Montréal (Québec), Canada to integrate a part into a sub-unit, it is often necessary to distort it. This is due to the highly overconstrained architecture of sub-units for aeronautical structures. If a part undergoes strain when it is integrated into a sub-unit this gives rise to mechanical stresses which could mean that the functional requirement of the sub-unit is not respected and could affect the mechanical behaviour of the part itself. With a classification of the structure components it is possible to identify those components for which significant strains are acceptable. These must then be quantified in order to define the induced stresses in the component and the stresses added to the sub-unit during the assembly phase. The aim of our study is to integrate these stresses generated in the part and the sub-unit by the assembly operation into the choices in the flow-process grid. The problems involved in integrating industrialisation into the design cycle are already taken into consideration in the expert rules provided by many aircraft manufacturers. These rules are the result of work by experts within the company. They have a limited area of validity and restrict the area available for design investigation. In this study our aim is to develop knowledge models in order to be more effective in finding the level of integration of the parts into a sub-unit, thus reducing mass and controlling costs.

CLASSIFICATION OF COMPOSITE STRUCTURAL COMPONENTS

A classification of the structural components identifies two types of component: type_1 and type_2. Each manufactured part is in contact with other parts in a sub-unit, which are called the adjacent parts. The adjacent parts will exert efforts via the contact surfaces, generating mechanical assembly stresses which are then distributed throughout the manufactured part and the sub-unit. These stresses can compensate for manufacturing deviations to a greater or lesser extent, in correlation with the strain caused by the shaping process. For type_1 components, mechanical assembly stresses have a great influence on the manufacturing deviations. Thus a fairly large defect can be specified in the general shape of the part in its free state. Stresses are characterised by geometric specifications associated with a stress condition [START_REF]Dimensioning and Tolerancing of Non Rigid Parts[END_REF], and the unit is shown on the definition of the part. However, mechanical assembly stresses should not exceed a certain maximum limit, determined by two criteria:

one criterion intrinsic to the manufactured part: mechanical behaviour, -one criterion intrinsic to the sub-unit: stress condition beyond which the sub-unit will be too distorted, thus jeopardising respect for certain functional requirements.

For type_2 components, assembly stresses have very little influence on strain generated by the process. As a result, manufacturing deviations cannot be compensated by contacts from adjacent parts. This results in a specification for a much more serious default in the general shape of the part in its free state than for the type_1 components. When the expert rules are applied to type_2 components, these components have symmetry from a geometrical point of view and also from the point of view of laminate architecture. This type of design limits the risk of strain but does not allow for the optimisation of mass and function integration. The use of specifications for each stress condition often proves to be unhelpful. In the design phase this distinction between groups of components will influence the laminate definition.

LAMINATE DESIGN AND RTM PROCESS

For parts with large dimensions, matter needs to be put in the places where stresses will occur hence the laminate is organised by zone. This is indispensable in order to reduce component mass. Each zone is a stack of fibers, with the juxtaposition of the zones ensuring variations in thickness. Each zone is defined by its shape, the number of fibers it contains and the orientation of each fiber in relation to a common frame of reference. Designing by zone in this way is linked with expert rules. In one of his papers, Gay describes an example of design by zone [START_REF] Gay | Matériaux Composites", 5e édition révisée[END_REF]. This organisation defines the optimised architecture for the laminate. To limit strain due to shaping, the unidirectional fibers are positioned using ply symmetry. Ply symmetry is the symmetrical arrangement of the stacks according to a design where the laminate is divided into two halves. In the RTM process, the stiffening pieces used are generally woven. This woven stiffener has a number of unidirectional continuous fibers or associated and oriented materials. In aeronautics, cost and delays in certifying composite material constituents can limit the types of stiffener that can be used. During the laminate design phase, an optimised architecture is defined for the laminate, however, using certified woven stiffeners may be the source of deviations between the optimised laminate architecture and the manufactured laminate architecture, and Fig. 1 illustrates this using an example showing the differences. To simplify the diagram, this example does not follow the rules of technological minimum [START_REF] Gay | Matériaux Composites", 5e édition révisée[END_REF]. This difference between the architectures impacts on the thickness of the laminate, and hence on the mass of the component and may give rise to strains. Organisation of the stacking is defined in a second phase in order to reach a compromise between reduction in mass and symmetry of stacking. We plan to use asymmetrical stacking for type_1 components. This asymmetry produces strains that are acceptable for type_1 components although not for type_2 components. Moreover, strains associated with using the RTM process are not caused solely by the architecture of the laminate. The different heat expansions of the constituent parts and the thermal cycles of polymerisation can also produce strain as parts are removed from the mould [START_REF] Ruiz | Numerical Analysis of Cure Temperature and Internal Stresses in Thin and Thick RTM Parts[END_REF].

Laminate architecture optimized

DIFFERENT GEOMETRIES FOR A TYPE_1 COMPONENT

During the industrialisation phase for a type_1 component, the aim is to estimate acceptable limits for manufacturing deviations in the RTM process according to the maximal limit of mechanical assembly stresses. As well as specifying the acceptable limits for manufacturing deviations, industrial practices and ISO standards for geometric specifications also lay down a stress condition that derives either from the criterion defined for manufactured parts or the criterion for the sub-unit, as defined in the paragraph above, "Classification of structural composite components". To quantify the strains associated with the shaping process we will first define several component states. These definitions will be based on behaviours of the product in the design cycle. Each state is associated with specific stress conditions (Fig. 2). Free state [START_REF]Dimensioning and Tolerancing of Non Rigid Parts[END_REF]: the manufactured component is subjected only to the action of gravity and to the residual stresses of the shaping process. The Sub-unit state corresponds to the mechanical stresses of assembly: the manufactured part is integrated into the sub-unit to which it belongs. Lastly, we define an Airplane state where the aircraft is in normal operational conditions. The Primary stresses state corresponds to the inspection stage, where manufacturing deviations in the part must respect the two criteria: assembly requirements and mechanical behaviour.

INDUSTRIALISATION MODEL AND SOLUTIONS SPACE

The level of functional integration of an aeronautical sub-unit is determined by the ratio of the number of components assembled to the number of functions the sub-unit must carry out. This level is optimal if the sub-unit consists of a single component. The architecture of a component can give rise to several industrialisation solutions all of which would enable the component to be produced and to carry out the functions for which it has been designed. Our model must therefore be able to define a group of solutions which are compatible with the multi-expert stresses used.

The stresses are formalised according to expert knowledge derived from the engineering and design department and the industrialisation departments and take the form of expert rules, arithmetic and logic relations. If expert software tools are used, then experimental designs can be set up from which response surfaces models can be generated. Each solution is defined by a free state. This set of solutions forms the industrialisation solutions space. The modelling variables (Fig. 3) are divided into four groups: design variables: engineering department and material (in the aeronautics industry this group is constrained by certification); industrialisation variables: process and laminate, in our study the sub-group process is limited to RTM [START_REF] Durai Prabhakaran | Design for 'X'-Abilities of RTM Products -A Graph Theoretic Approach[END_REF];

A0 Define a space of solutions of industrialization

Performance variables

-injection feasibility, -manufacturing deviations, -level of functional integration.

Manufacturing variables

Solveur of constraint satisfaction problem

Design variables

Manufacturing constraints

Engineering department

Material

-material characteristics. laminate -laminate zone, -laminate architecture.

-permeability.

-Expert rules : deviations characterization RTM.

-Expert rules : manufacturing architecture definition of laminate.

-Expert rules : characterization of the cuts which represent areas of different thickness.

-Expert rules : process cycles RTM.

-Surfaces of responses.

-Arithmetic and logic relations.

Geometry free state

Manufacturing solution

Process chosen RTM

performance variables: functional integration, feasibility injection and manufacturing deviations; auxiliary variables: set of variables needed to define constraints [START_REF] Rais-Rohani | Analysis and Optimisation of Primary Aircraft Structures Based on Strength, Manufacturing, and Cost Requirements[END_REF].

The resolution strategy based on exploring a solutions space produces value domains {D} associated with the variables {V}, these may be design variables, industrialisation variables and performance variables. Each solution thus defined will satisfy all the constraints of the problem {C} and will be formalised by the free geometric state of the component. This set of variables and constraints represents a constraint satisfaction problem. Values associated with the performance variables (performance indicators) will guide choice in the solutions space. Due to the heterogeneity of the variables being manipulated, the use of multi-expert constraints and the exploration of a solutions space, we considered using a constraint satisfaction problem solver [START_REF] Scaravetti | Structuring of Embodiment Design Problem Based on the Product Lifecycle[END_REF]. Fig. 4 shows the processing of variables and their values, resulting in an industrialisation solution. The variables below are generally sets of elementary variables which define a parameter that is relevant to the design cycle for the component:

nominal geometry: defined by a set of variables associated with the initial shapes and dimensions of the component. A second set of variables defines the shapes, dimensions and positions of additional geometries to the component's initial geometry, -contact surface: defined by a set of variables describing the shapes, dimensions and positions of the contact surfaces, the number of variables may change with the level of functional integration, -component type: this variable has two values (type_1 or type_2), -material characteristics: defined by a set of variables describing the type of stiffener and the type of resin, -laminate zone: defined by a set of variables describing the different shapes that are possible at the ends of the zones, -laminate architecture: defined by a set of variables describing the stacking of the fibers, independently of the type of stiffener and the type of component, permeability: defined by a set of variables describing the types of resin flow through the preform in the tools, injection feasibility: defined for each solution by a variable whose value is a function of a porosity rate, manufacturing deviations: defined for each solution by a variable whose value is a function of the stress condition, functional integration: defined for each solution by a variable whose value is assessed based on the number of functions carried out by the component. A constraint satisfaction problem solver is an effective procedure to apply to an industrialisation model, taking into account causality when processing the variables and then prioritising these variables. In terms of digital processing, our work consisted of proposing a modelling framework (meta-model). Using meta-models, the performance of constraint satisfaction problem solvers can be improved in terms of calculation time and solution quality [START_REF] Chenouard | Résolution par Satisfaction de Contraintes Appliquée à l'Aide à la Décision en Conception Architecturale[END_REF]. In our study, three stages were necessary to finally define a solution: -definition of two characteristics relating to the laminate and one relating to the geometry of the component: definition of the stacking of the laminate according to the nature of the stiffener and the component type and definition of the geometry of the end of the fibers that make up the zones, definition of the nominal geometry of the component adapted to the process (joint radius, maximum dimensions, etc.). -estimation of solution injectability using a software tool. For this stage, permeability values must be defined, these values will be obtained experimentally or semi-analytically. Injectability can be estimated by measuring porosity rates from injection simulations. -definition of admissible manufacturing deviations which depend on the condition imposed on the constraint when defined industrially. To obtain this value, solution injectability must be greater than a fixed percentage of porosity rates; if this is not the case the solution is abandoned. Results from the first stage combined with the characteristics of the material will produce a definition of manufacturing deviations. This definition guarantees the mechanical behaviour of the component and of the structure into which it will be inserted.

One of the major difficulties found with this type of constraint satisfaction problem is that the set of variables, the set of constraints and the set of domains may evolve during the solving process. This is traditionally identified as a dynamic constraint satisfaction problem. Moreover, the value domains associated with the variables may change across continuous and discrete spaces. Our problem falls into the category of dynamic mixed constraint satisfaction problem [START_REF] Rossi | Handbook of Constraint Programming[END_REF].

INDUSTRIAL APPLICATION

This study is part of a collaborative work with the Dassault-Aviation company. We chose to work on a component currently produced by RTM, i.e., the spars that make up the horizontal plane of the Falcon [START_REF] Alix | RTM Technology Development: New Methodologies for Process Control Improvement[END_REF]. Assembly of the sub-unit that includes the spars is done using a rigid assembly tool, the spars are type_1 components. The spars are made of carbon fibre and epoxy resin. The solutions space defined by our industrialisation model will be compared with industrial solutions in order to estimate the validity of the modelling.

CONCLUSIONS AND FUTURE PERSPECTIVES

The aim of this aid to industrialisation is to explore a solutions space obtained using a constraint satisfaction problem solver. Constraints are a combination of design and industrialisation variables where each solution is associated with performance indicators. The final choice depends on industrial strategies. This aid to industrialisation can be used during several phases in the life cycle of a product: -in detailed design, during exchanges between the design department and the industrialisation department. Using this tool it is possible to understand the effect of integrating new functions into the free state geometry of the component and hence estimate repercussions on the RTM shaping phase and the assembly phase, -in industrialisation, by using specifications for non-rigid parts, to take into account strains due to the process of defining and implementing compliance criteria for manufactured parts, -in manufacture, to reduce the requirements in terms of geometry for the RTM process.

This work could be associated with post-RTM strain predictions [START_REF] Dong | Dimension Variation Prediction and Control for Composites[END_REF]. If strains could be predicted, then a tool could be produced to guarantee that the component's free state is respected.

Fig. 1

 1 Fig. 1 Optimized and manufactured laminate architectures.

Fig. 2

 2 Fig. 2 Different states of a type_1 component.

Fig. 3

 3 Fig. 3 Modelling to determine industrially feasible solutions (IDEF0 formalism).

Fig. 4

 4 Fig. 4 Meta-model to determine industrialisation solutions (IDEF0 formalism).

ACKNOWLEDGEMENTS

I would like to thank Cyrille Clermont from Dassault-Aviation for his contribution to this study.


Nathalie Villa-Vialaneix
email: nathalie.villa@univ-paris1.fr.

Matthieu Vignes

Nathalie Viguerie

Magali San

Inferring Networks from Multiple Samples with Consensus LASSO

Keywords: network inference, Gaussian graphical model, regulation network, LASSO

solution whilst a second penalty is used to make condition-specific networks consistent with a consensual network. This "consensual network" is introduced to represent the dependency structure between genes, which is shared by all conditions. We show that different "consensus" penalty can be used, some integrating prior (e.g., bibliographic) knowledge and others that are adapted along the optimization scheme. In all situations, the proposed double penalty can be expressed in terms of a LASSO problem and hence, solved using standard approaches which address quadratic problems with L 1 -regularization. This approach is combined with a bootstrap approach and is made available in the R package therese 1 . Our proposal is illustrated on simulated datasets and compared with independent estimations and alternative methods.

It is also applied to a real dataset to emphasize the differences in regulatory networks before and after a low-calorie diet.

Abstract: Networks are very useful tools to decipher complex regulatory relationships between genes in an organism. Most work address this issue in the context of i.i.d., treated vs. control or time-series samples. However, many data sets include expression obtained for the same cell type of an organism, but in several conditions. We introduce a novel method for inferring networks from samples obtained in various but related experimental conditions. This approach is based on a double penalization: a first penalty aims at controlling the global sparsity of the

Introduction

The recent development of high-throughput techniques produces huge datasets where thousand of gene expressions are simultaneously measured. However, the number of observations is comparatively very small, and those are often measured in a variety of experimental conditions. One of the big challenges of modern Systems Biology is to understand the influence of controlled experimental conditions on the functioning of living organisms. This question is usually addressed by searching for the differences between gene expressions pertaining to the conditions (hence for "differentially expressed genes"). A more comprehensive look at the roles of the genes of an organism can be obtained by deciphering the interactions of these genes with each other; finding which regulation pathways are modified by a given experimental condition gives an interesting insight on the influence of the condition on the living system as a whole.

One of the most popular approach to understand the complex relationships existing between the expression of a large set of genes is to infer a co-expression network from a transcriptomic dataset. In such a model, nodes of the network represent the genes and an edge is meant to stand for a regulatory link between the two nodes it connects. A large number of different methods have been proposed to infer such networks: using correlations ("relevance network", [START_REF] Butte | Mutual information relevance networks: functional genomic clustering using pairwise entropy measurements[END_REF]), Bayesian networks [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF][START_REF] Pearl | Bayesian Networks[END_REF],

Gaussian Graphical Model [START_REF] Edwards | Introduction to Graphical Modelling[END_REF][START_REF] Schäfer | An empirical bayes approach to inferring largescale gene association networks[END_REF]... When observations are collected in different conditions, a naive approach would be to independently infer a network for each condition and to compare them. However, this method is not suited to highlight specific differences and shared motifs of regulation phenomenons. Moreover, since the number of observations is often too small, inferring networks independently (assuming that a common functioning exists in most scenarii) leads to emphasize irrelevant differences. Several proposals have already been made to overcome this issue: [START_REF] Chiquet | Inferring multiple graphical structures[END_REF][START_REF] Danaher | The joint graphical lasso for inverse covariance estimation accross multiple classes[END_REF][START_REF] Mohan | Structured learning of Gaussian graphical models[END_REF] use a modified Gaussian graphical model and [START_REF] Jung | Inférence conjointe de réseaux de gènes dans de multiples états[END_REF] proceeds with two steps with a clustering prior to the inference. The proposal developed in this paper is close to that of [START_REF] Chiquet | Inferring multiple graphical structures[END_REF][START_REF] Danaher | The joint graphical lasso for inverse covariance estimation accross multiple classes[END_REF][START_REF] Mohan | Structured learning of Gaussian graphical models[END_REF]: a Gaussian graphical model is used and two interpretable penalties are added to the likelihood. The first penalty aims at inferring sparse solutions; the second penalty is used to make networks obtained in different conditions consistent with a consensual network. The "consensual network" is introduced to represent high-level dependencies between genes, i.e., a common functioning of the living organism under study, in most situations. It can either include prior (e.g., bibliographic) knowledge or be expressed from the condition-specific networks. Finally, the estimation is made more robust by using a bootstrap approach.

The paper is organized as follows: Section 2 describes the double penalty approach. Section 3 explains our proposal for estimating the networks with a bootstrap strategy. Finally, Section 4 provides experimental results on simulations.

cLasso

In the Gaussian graphical model (GGM) framework, the classical objective is to estimate the graph of conditional dependencies between p variables (usually modeling gene expressions), (X j ) j=1,...,p , from n i.i.d. observations of the variables, namely (X ij ) i=1,...,n , ∀ j ∈ {1, . . . , p}. Each p-dimensional vector X i. is assumed to be the realization of a Gaussian random variable N (0, Σ). In this framework, non-zero entries of the concentration matrix K = Σ -1 exactly encode actual edges (between genes) in the conditional dependency graph. In the present section, we describe how this framework can be extended to the case where observations are obtained from different samples, each sample being measured in a given (but related) experimental condition.

Inferring multiple networks with GGM

Now assuming that the p gene expressions are measured from k samples, each corresponding to a specific experimental condition, the following model can be set: (X c j ) j=1,...,p, c=1,...,k are k Gaussian p-dimensional vector, N (0, Σ c ). A total of n observations are available: (X In the following, our goal is twofold:

• inferring k sparse graphs that model gene regulations in the k conditions;

• finding one consensual graph that models a "shared" functioning between conditions.

The GGM framework is used for the inference. As previously explained, the concentration matrices K c = (Σ c ) -1 need be estimated and the entries of these matrices exactly measure conditional dependencies between variables (X c j ) j through partial correlation coefficients,

s c jj ′ = Cor X c j , X c j ′ |(X c l ) l =j, j ′ because of the rela- tion s c jj ′ = - K c jj ′ √ K c jj K c j ′ j ′ [26].
These quantities can be estimated by considering the following (k × p) linear regression problems [START_REF] Meinshausen | High dimensional graphs and variable selection with the lasso[END_REF]:

∀ c = 1, . . . , k, ∀ j = 1 . . . , p, X c j = X c \j β c j + ǫ c j , (1)
where X c \j is the matrix X c = X c ij i=1,...,nc, j=1,...,p deprived from its j-th column

X c j , β c j = (β c jj ′ ) j ′ =j is a (p -1)
-dimensional vector and ǫ c j is a Gaussian centered error. In the Gaussian framework, it can be shown that the coefficients of the linear model are related to the previous quantities by

β c jj ′ = - K c jj ′ K c jj .
The k × p linear models of Equation ( 1) can be jointly estimated by maximizing a pseudo-likelihood:

L(K|X) = k c=1 p j=1 nc i=1 log P X c ij |X c i,\j , K c j .
(2) [START_REF] Chiquet | Inferring multiple graphical structures[END_REF] proved that maximizing the pseudo-likelihood of Equation ( 2) over matrices (K c ) c is equivalent to minimizing the following p quantities simultaneously:

∀ j = 1, . . . , p, 1 2 β T j Σ \j\j β j + β T j Σ j\j . ( 3
)
The p problem of Equation ( 3) are (p -1) × k-dimensional quadratic optimization problems and we specify:

• ∀ c = 1, . . . , k, β c j = (β c jj ′ ) j ′ =j ∈ R p-1 , where β c jj ′ = (K c ) -1 jj K c jj ′ ; • β j = β 1 j , . . . , β k j T ∈ R k(p-1) ; • Σ \j\j is the block diagonal matrix Σ \j\j = Diag Σ 1 \j\j , . . . , Σ k \j\j , having dimensions k(p -1) × k(p -1); • Σ j\j is the k(p -1)-dimensional vector, Σ 1 j\j , . . . , Σ k j\j .
However, this approach leads to matrices without non-zero entries. Moreover, when (n c ) c are not larger than p, the estimation of (β j ) j becomes trickier and pseudo-inverse methods lead to highly unstable results. Using the additional assumption that conditional dependency graphs are indeed sparse, a standard approach is to add a L 1 -penalization to the likelihood of Equation (2) ("Graphical LASSO", see [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]) or, alternatively, to consider p independent L 1 -penalized problems derived from those of Equation (3), see [START_REF] Chiquet | Inferring multiple graphical structures[END_REF][START_REF] Meinshausen | High dimensional graphs and variable selection with the lasso[END_REF]. The latter, more direct approach, has been reported to be more accurate in terms of edge detection in [START_REF] Villers | Assessing the validity domains of graphical Gaussian models in order to infer relationships among compo-nents of complex biological systems[END_REF].

Using a "consensus" penalty

In the previous section, the conditional dependency graphs are obtained from each sample independently. The assumption that the graphs issued from the different experimental conditions should be somehow alike, is not integrated into the model.

Especially in the case where the sample sizes are low, such an assumption should help to predict edges more accurately. Various techniques exist to address this issue: [START_REF] Chiquet | Inferring multiple graphical structures[END_REF] proposed to replace the covariance matrices Σ c by mixing it with the covariance matrices corresponding to the other conditions. Alternatively, some authors suggest to penalize the pseudo-likelihood by a penalty that can explicitly deal with the similarity between condition-specific graphs via different strategies:

• [5]
proposed two kinds of Group-LASSO type penalties:

P ((K c ) c ) = ij c (K c ij ) 2 (Group-LASSO) and P ((K c ) c ) = ij c (K c ij ) 2 + + c (K c ij ) 2 -
(sign-coherent Group-LASSO or "Cooperative LASSO"). The group-LASSO penalty globally controls sparsity and inferred edges are common to all conditions. The sign coherent option of their penalization scheme offers the possibility to enforce an edge to encode either an activating or repressing process but not both: it provides strongly similar networks between conditions and has been proven to be efficient in case of experimental conditions leading to small changes in the regulations.

However, for some particular applications (e.g., certain forms of cancer that lead to a complete re-organization of the living system), the assumption that the relations between two genes is always a repressing/enhancing relation is not biologically desirable;

• [6] used the penalty

P ((K c ) c ) = c =c ′ K c -K c ′ 1 ,
where . 1 is the standard L 1 -norm, which commands a strong similarity across conditions. This approach would lead to very similar condition-specific network, allowing only a few differences. Unlike the Cooperative Lasso approach described above, no special sign-coherent assumption is required but this method is more suited when condition-specific networks are not supposed to be very different;

• [START_REF] Mohan | Structured learning of Gaussian graphical models[END_REF] introduced the penalty

P ((K c ) c ) = c =c ′ p j=1 j K c j -K c ′ j 2 , where
. 2 is the standard L 2 -norm and K j is the j-th column of K. Hence, this approach encourages the support of K c -K c ′ to be the union of a given set of columns. Hence, this penalty only provides some flexibility to a few nodes to differ among conditions while all the other nodes have the same pattern of interactions.

The main idea of our proposal, we coined cLasso, is similar to the latter approaches, but using a softer penalization scheme than group-Lasso type penalties.

This choice aims at better estimating the edges that are not similar and also do not need to assume a particular origin for the differences between conditions. The k inferred graphs, G c , are forced towards a "consensual" graph: the resulting graphs are different from each others, but these differences can be controlled. This idea is tackled by using a penalized ML framework in which two penalties are introduced:

• the first one is a sparse penalization which controls the number of edges in every graph G c ;

• the second one is a L 2 penalization that aims at limiting the differences between the (K c ) c=1,...,k .

More precisely, ∀ j = 1, . . . , p, a consensual regression coefficient, β cons j , is introduced, that can be defined from the sample-dependent coefficients β c j or can be fixed by a user, including, in particular, prior biological knowledge. This coefficient represents a kind of "global" solution, that is condition-independent. It is used by replacing the minimization problems described in Equation ( 3) by the following double-penalized minimization problems:

1 2 β T j Σ \j\j β j + β T j Σ j\j + λ β j 1 + µ k c=1 β cons j -β c j 2 2 . ( 4
)
In Equation ( 4), β cons j is used to model the "consensus". In the following section, different types of consensus are described, and the practical computation of the solution is derived from the different cases. All described solutions lead to the optimization of quadratic problems penalized by the L 1 -norm.

Contrary to the other approaches presented above, the second penalty of Equation ( 4) is a soft one, that does not control drastically the number of different edges between conditions but rather limits them. It is thus advisable in the case where the number of differences is not too low and where the user really wants to see the differences across the conditions. Also, contrary to [START_REF] Danaher | The joint graphical lasso for inverse covariance estimation accross multiple classes[END_REF][START_REF] Mohan | Structured learning of Gaussian graphical models[END_REF], our proposal does not rely on a penalty which complexity increases quadratically with the number of conditions (this might be a problem if the number of conditions is high). Finally, as explained in Section 2.3.1, the definition of a consensus network can integrate prior knowledge that can help estimating the network with an increased accuracy.

Remark 1 As shown in Section 3.1, any choice for β cons j that leads to obtain a minimization problem that can be expressed as:

convex part + λ β j 1
is a valid consensus choice that can be solved using a common framework. In particular, this includes any consensus that is expressed as a linear combination of the estimated coefficients β c j (Section 2.3.2) or (fixed) a priori consensus (Section 2.3.1).

Consensus choices 2.3.1 A fixed consensus

When a prior information is known on the network (e.g., a bibliographic network), a natural choice is to use it for β cons . In this case, β cons is fixed in advance and does not depend on (β c j ) j : it does not need to be estimated. However, if no prior information is available, the network estimated from all the samples considered as a whole or any combination of networks obtained with independent estimations can be used for consensus and considered as a (fixed) a priori information network.

Proposition 1 Using a fixed β cons j , Equation ( 4) is equivalent to minimizing the following standard quadratic problem with L 1 -penalty:

1 2 β T j B 1 (µ)β j + β T j B 2 (µ) + λ β j 1 , ( 5
)
where

B 1 (µ) = Σ \j\j + 2µI k(p-1)
, with I k(p-1) the k(p -1)-identity matrix and

B 2 (µ) = Σ j\j -2µI k(p-1) β cons j with β cons j a k(p -1)
vector that only depends on the prior β cons j .

Proof (and exact values for B 1 and B 2 ): The L 2 -penalty of Equation ( 4) can be re-written as:

k c=1 β cons j -β c j 2 2 = k c=1 (β c j ) T β c j -2(β c j ) T β cons j + β cons j 2 2
Noticing that β cons

j 2
2 is a fixed value that does not depend on the estimated coefficients β c j , it follows that minimizing Equation ( 4) is equivalent to minimizing:

1 2 β T j Σ \j\j + 2µI k(p-1) β j + β T j Σ j\j -2µI k(p-1) β cons j + λ β j 1 ,
where β cons j is the vector (β cons j ) T , . . . , (β cons j ) T T .

An averaged consensus

When no prior information is given, an intuitive and convenient choice for the consensus is to simply average the estimators over the different samples:

β cons j = k c=1 nc n β c j .
In this case, β cons j is a linear combination of the (β c j ) c , which is an interesting feature, as explained in Proposition 2. Notice that the choice of averaging the coefficients β c j is almost equivalent in terms of networks (i.e., in terms of nonzero entries) as having a consensus which is the union of the condition-dependent networks. 4) can be re-written as the following standard quadratic problem with L 1 -penalty:

Proposition 2 Using β cons

j = k c=1 nc n β c j , Equation (
1 2 β T j S j (µ)β j + β T j Σ j\j + λ β j 1 ( 6
)
where S j (µ) = Σ \j\j + 2µA T A where A is a k(p -1) × k(p -1)-matrix that does not depend on j.

Proof (and exact value for

A): If, ∀ c = 1, . . . , k, U c = nc n I p-1 (with I p-1 the unit matrix having dimension p -1) and V c = 1 -nc n I p-1 , then β c j -β cons j = A c β j ,
where

A c is the (p -1) × k(p -1)-matrix [-U 1 , . . . , -U c-1 , V c , -U c+1 , ..., -U p ].
Then,

k c=1 β cons j -β c j 2 2 = k c=1 β T j A T c A c β j
and thus, setting

A =         A 1 . . . A k         , implies that k c=1 β cons j -β c j 2 2 = β T j A T Aβ j
which concludes the proof.

Remark 2 Because the term β T j A T Aβ j is a quadratic term in β, the formulation of the minimization problem given in Equation ( 4) is not a direct penalization of the ML optimization. More specifically, minimizing Equation ( 4) is equivalent to minimizing the following penalized ML:

L(K|X) -λ K 1 - µ n k c=1 (D c ) -1/2 K cons,c -K c 2 2
where

• K 1 = k c=1 K c 1 = k c=1 p j,j ′ =1 |K c jj ′ |; • D c = Diag(K c 11 , K c 22 , . . . , K c pp ); • K cons,c j\j = k t=1 nt n K c jj K t jj K t j\j .
Note that, as explained in [START_REF] Chiquet | Inferring multiple graphical structures[END_REF], estimating (K c jj ) j is not relevant to unveil the graph structure so, in practice, these values are set equal to Σ -1 jj . Hence, from the ML point of view, there is no definition of a consensual concentration matrix since this quantity depends on the sample (the average is weighted differently depending on the sample).

In practice, in every task, the variables are previously scaled and K c jj are all set equal to one, which leads to the following equivalent formulation of the optimization problem

L(K|X) -λ K 1 - µ n k c=1 K cons -K c 2 2 ,
where

K cons j\j = k t=1 nt n K t j\j .
Remark 3 The penalty of [START_REF] Mohan | Structured learning of Gaussian graphical models[END_REF] can be re-written as:

c =c ′ K c j -K cons j + K cons j -K c ′ j 2 = c =c ′ K c j -K cons j 2 + K c ′ j -K cons j 2 + 2 K c j -K cons j , K cons j -K c ′ j = (k -1) k c=1 K c j -K cons j 2 + 2 c =c ′ K c j -K cons j , K c ′ j -K cons j
Then, in the case of the averaged consensus, an edge (j, j ′ ) is in the consensus network if and only if it is in at least one of the condition-specific networks. In

particular, for k = 2, K 1 jj ′ -K cons jj ′
= 0 means that (j, j ′ ) is an edge in the consensual network and is not an edge in the network specific to condition 1. It is thus also an edge in the network specific to condition 2 (as there is only two conditions). In conclusion, when k = 2, K 2 j -K cons j , K 1 j -K cons j = 0 and thus the consensus penalty is very similar to the penalty proposed in [START_REF] Mohan | Structured learning of Gaussian graphical models[END_REF]. However, for k > 2, the situation might be more complicated: a condition-specific edge can be specific to more than one condition and thus the equality

K c j -K cons j , K c ′ j - K cons j
= 0 is no more guaranteed for the averaged consensus. Conversely, nullity of the scalar product (and thus the equivalence between the consensus penalty and the penalty proposed by [START_REF] Mohan | Structured learning of Gaussian graphical models[END_REF]) would be obtained, for instance, if an edge that is specific to a condition is present in only one of the condition-specific networks.

This property does not seem to be desirable on a biological point of view.

Computational aspects

This section will provide computational details on the cLasso methods. First, the method used to solve the optimization problems introduced above is described and then, a bootstrap approach is introduced to help decreasing the false positive rate and to help increasing the prediction accuracy when dealing with small sample size problem.

cLasso optimization

The cLasso problem is solved by minimizing the p sub-problems of Equations ( 5)

and [START_REF] Danaher | The joint graphical lasso for inverse covariance estimation accross multiple classes[END_REF]. The objective function of all the problems that can be decomposed into:

a convex part C(β j ) = 1 2 β T j Q 1 j (µ)β j + β T j Q 2 j (µ)
, convex in (β j ) j and that does not depend on λ; a norm penalty P(β j ) = c 1 nc β j 1 that is non differentiable at 0, with respect to β j .

The non differentiability of P shrinks the LASSO estimate toward 0 and potentially sets β jl = 0 for several indexes l, as explained in [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. In [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], the LASSO optimization problem is solved by a quadratic programming method, which is used to perform the estimation of the (β jl ) l together with a variable selection. Since then, several authors have proposed more efficient approaches to solve the LASSO optimization problem: [START_REF] Fu | Penalized regressions: the bridge versus the lasso[END_REF] developed the so-called "shooting algorithm" that starts from an unconstrained least-square solution and uses a coordinate descent. Unfortunately, this algorithm is not applicable in the case of sparse problems as soon as n < p. Others proposed to use differentiable approximations of P, such as [START_REF] Lee | Efficient l 1 regularized logistic regression[END_REF] that takes advantage of the approximation β 1 ≃ j β 2 j + ǫ. Finally, [START_REF] Osborne | On the LASSO and its dual[END_REF] uses a method that is efficient for medium-size problems and suited to the case n < p. We used a similar strategy, which is close to the one described in [START_REF] Chiquet | Inferring multiple graphical structures[END_REF]: it is based on a greedy update of an "active set" that progressively gathers together all non-zeros coefficients of the different sub-problems. At each step of the algorithm, the coefficients are estimated only for the variables that are included in the active set.

More precisely, for a fixed value of λ, starting from a vector β A j of non-zero coefficients on the active set A, the method first solves the so-called "master problem" given by Equations ( 5) or [START_REF] Danaher | The joint graphical lasso for inverse covariance estimation accross multiple classes[END_REF], which is differentiable, because, by definition, the coefficients of β A j are not null. This is done by using the sub-gradient

∂ β j [C(β j ) + λP(β j )].
Then, the set of active variables is updated by adding the variables that violate the most the first-order optimality condition. The algorithm stops when

• for all l ∈ A, β jl = 0 and

∂ β j [C(β j ) + λP(β j )] l = Q 1 j (µ)β j + Q 2 j (µ) l + λ sign(β jl ) = 0
• for all l ∈ A, β jl = 0 and

Q 1 j (µ)β j + Q 2 j (µ) l ≤ λ
Further details on the method can be found in [START_REF] Chiquet | Inferring multiple graphical structures[END_REF], that uses the same optimization scheme for the so-called "intertwined LASSO" method.

Finally, the method is applied to a whole set of λ values, starting from largest (i.e., from the one that yields to the strongest constrain) and using the optimal β j as a prior for solving the problem with the next smaller λ. This method is implemented in the R package therese, that can be downloaded from http:// therese-pkg.r-forge.r-project.org/.

Bootstrapped cLasso

As demonstrated in [START_REF] Bach | Bolasso: model consistent lasso estimation through the bootstrap[END_REF], the LASSO converges to the selection of all the variables included in the true model (true positives) with probability one but asymptotically selects all other variables (false positives), with a strictly positive probability. In practice, this means that using the LASSO algorithm yields to a rather high number of false positive edges in the network estimation. To overcome this difficulty, [START_REF] Bach | Bolasso: model consistent lasso estimation through the bootstrap[END_REF] proposes the so-called "Bolasso" method, that combines LASSO performed on bootstrap samples. Bootstrapping [START_REF] Efron | Nonparametric estimates of standard error: the jackknife, the bootstrap and other methods[END_REF] is a resampling technique that consists in creating new samples of the same size as the original by sampling randomly with replacement from the original dataset. Its aim is to estimate the sampling distribution of almost any statistics and thus to estimate the accuracy for these statistics.

In Bolasso, LASSO is run on a large number of bootstrap samples and the intersection of the variables selected in every bootstrap sample are finally kept. It is proved that this approach is a consistent model selection method.

Hence, in order to improve the false positive rate of the approach described above, we use a similar methodology, only taking into account the fact that the typical sample size in transcriptomic experiments is far from being close to the asymptotic case. More precisely, instead of intersecting the edges selected in every bootstrap sample, the number of times an edge is selected by all computations run on each bootstrap sample is used as a quality measure of the edge. Only the most frequently selected edges, those that are selected more than a given number denoted by T 2 ∈ {1, . . . , N boot }, are finally included in the estimated network.

In practice, for every bootstrap sample, Equation ( 5) or Equation ( 6) is solved for a list of several values of λ and a fixed value for µ, using the method described in Section 3.1. A given value of λ, depending on the bootstrap sample, is retained which corresponds to the first time in the path (i.e., to the largest λ) for which the number of estimated edges is larger than a target value, T 1 . T 1 is fixed to a rather high value to avoid missing relevant edges. The complete procedure is described in Algorithm 1.

The impact of T 1 and T 2 is discussed further in the simulations of Section 4.1.2. N c (j, j ′ ) ← N c (j, j ′ ) + 1 end for 12: end for 13: return List of edges for sample c: {(j, j ′ ) : N c (j, j ′ ) > T 2 }

Application

The simulations described in this section have been performed using R version 3.0 and the packages glasso2 , SIMoNe3 and JGL4 . Bootstrap was performed using a parallel implementation with the package doMC.

Simulated data

The method is first illustrated on simulated data. These experiments use one of the graphs provided at http://www.comp-sys-bio.org/AGN/data.html and created by Pedro Mendes (Virginia Bioinformatics Institute and State University; see [START_REF] Hoops | COPASI -a COmplex PAthway SImulator[END_REF]). More precisely, the graph "scale-free Century 007"5 was used to test the method. This network has 100 nodes (corresponding genes) and 200 edges (corresponding to gene interactions): the density of the (undirected) network is thus approximately equal to 4%. The term "scale free" indicates that the network has been generated from a preferential attachment model, as described in [START_REF] Barabási | Emergence of scaling in random networks[END_REF].

Additionally, the edges of the network are colored: half are "red" and half are "blue", which will differentiate a positive from a negative correlation between two variables.

Data generation

Several artificial expression datasets were generated from the graph described above. More precisely,

• k child networks were created by randomly rewiring a given ratio r of the edges of the original network. Hence, two child networks have approximately 100(1 -2r)% of shared edges. Loops and multiple edges were forbidden during the rewiring process but the color of the edges was preserved. Each of these k networks is used to model one experimental condition;

• n c expression data were then generated from a Gaussian multivariate variable with a covariance matrix Σ for which the conditional dependency structure corresponded to one of the child network. In addition, the edge colors were used to define the sign of the partial correlation: blue edges corresponded to a negative partial correlation (mimicking inhibition) and red edges (mimicking activation) to a positive one.

Several experiments were designed with various values for k (2, 4 or 5), r (varied between 5%, 10% and 20%) and the respective sample sizes n 1 , . . . , n k : 2 × 20, 2 × 30, 2 × 50, 5 × 20 and 4 × 30. Only small sample sizes (no less than 50 observations) were used to fit realistic experimental conditions in which only a few observations per condition are generally available. The resulting child networks had no more than 40% of different edges.

Figure 1 illustrates the generation process on an example: the "scale free Century 007" graph is displayed as well as two of its children, obtained by rewiring 5% of the edges.

Bootstrap analysis

In this section, we investigate the effect of T 1 and T 2 on the performance of the algorithm. This analysis is made using the results obtained from the expression data generated with 5% of rewired edges, 2 conditions, each containing 20 observations. For this network, 100 bootstrap samples were extracted: this number is low compared to standard recommendations but, for one hand, the approach is computationally expensive and, for the other hand, a previous work [START_REF] Allouche | A panel of learning methods for the reconstruction of gene regulatory networks in a systems genetics context[END_REF] showed that the benefit of bootstrapping was achieved with the combination of 30 to 40 bootstrap samples. Also, several values of T 1 and µ have been tested: T 1 ∈ {250, 300, 500} and µ ∈ {0.1, 1}. The performance of the different parameters are compared by Figure 1: The "scale free Century 007" graph and two resulting child networks, obtained with 5% of rewired edges. Green dotted edges are shared edges whereas red solid ones are condition specific edges. The vertex positions result from a force-directed placement algorithm as in [START_REF] Fruchterman | Graph drawing by force-directed placement[END_REF] and are common to all three networks so that the edges can easily be compared.

means of the F statistics:

F = 2 × precision × recall precision + recall ,
where the precision is the ratio of retrieved edges that are in the true network (true positive edges among positive edges) and the recall is the ratio of true edges that are retrieved by the method (true positive edges among the edges in the original network). F is the harmonic mean of the precision and of the recall and computes a trade-off between the two quantities.

For each condition and each pair of parameters {T 1 , µ}, the F statistics were calculated along the precision/recall values obtained for different values of values of T 2 (bootstrap estimation). Then, the pairs of parameters {T 1 , µ} were compared based on the averaged F over the conditions: the "best pair" of parameters is the one that maximizes the maximum averaged F along the path of T 2 values, the maximum F being used as a way to find the best compromise between precision/recall. According to this method, the best pair {T 1 , µ} for the expression data described above was {500, 1}, as shown in the level plot of Figure 2 When the value of T 1 is set to a rather high value, 500 (which is much larger than the true number of edges), and when µ is equal to 1, Figure 3 over the 100 bootstrap samples and correspond to a precision about 25% and a recall about 30%. These points give inferred networks with a resulting density slightly lower that the true network density (2.5-3.5% instead of 4%). This illus- trates the fact that, if there is a prior knowledge on a targeted density, a good strategy could be to choose T 2 so that resulting networks fit this targeted density.

When the value of T 1 is equal to a smaller value (250 which is larger than the true number of edges), the evolution of the density versus T 2 and the precision/recall curve are given in Figure 5. The conclusions are very similar except, of course, that for a given value of T 2 , the densities of the resulting networks are lower. Otherwise, the distributions of the number of times a given edges is selected by the algorithm in the bootstrap samples are quite similar and the best F value is also obtained for networks that have densities slightly slower than the true density. However, as shown in Figure 2, T 2 is a less important parameter for the method performance, as compared to µ. Optimal parameters, according to the maximum F statistics, are given in Table 1 for all simulations. As expected, µ needs be smaller in the case where the two conditions correspond to more different networks (i.e., when the number of rewired edges is larger) but generally, using a rather high value for T 1 is the strategy that provides the best results. The effect of the bootstrap on the performance is shown in the last column of this table, which contains the percentage of increase of the corresponding maximum F compared to the direct approach. Bootstrap only improves the performances when the percentage of rewired edges is moderate (lower than 10%) or when there are many different conditions. The counter-performance of bootstrapping could be explained by the fact that it enforces the joint effect and thus fails to estimate edges specific to the condition, that are less numerous in those cases. Additionally, this might be a very high-dimension issue [START_REF] Verzelen | Minimax risks for sparse regressions: ultra-high-dimensional phenomenons[END_REF]: when the n p ratio allows us to draw model inference but is at the limit of producing reasonably accurate estimates, the use of a bootstrap procedure produces a set of highly unstable estimates, which lead to fewer robust estimated edges. As a consequence, the model estimate might focus on those edges which are supported by many conditions and does not detect finer pattern of dependencies in the data. 1: Best parameters of the bootstrap cLasso for each simulation according to the maximum F along the path of T 2 values, and percentage of increase of the best F value compared to the direct (i.e., un-bootstrapped) approach.

µ T 1 %

Performance comparison

In this section, cLasso is compared to alternative methods for inferring graphs from expression data. More precisely, for each expression dataset described in Section 4.1.1, the following methods are applied to infer the k conditional dependency networks corresponding to the k different conditions:

• the graphical Lasso method, as described in [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] and hereafter denoted by gLasso: the k networks corresponding to the k different conditions are inferred independently. Hence, the comparison with this method aims at showing the effect of jointly inferring the networks instead of independently;

• the intertwined Lasso, the cooperative Lasso and the group Lasso methods, as described in [START_REF] Chiquet | Inferring multiple graphical structures[END_REF] and hereafter denoted by iLasso, coopLasso and grou-pLasso, respectively. These methods are used to provide a comparison with other joint inference methods. Also notice that the data generation provides sign-coherent networks (i.e., the different child networks are very likely to have the same sign for partial correlations corresponding to shared edges), which should favor the cooperative Lasso method;

• the fused graphical Lasso, as described in [START_REF] Danaher | The joint graphical lasso for inverse covariance estimation accross multiple classes[END_REF], denoted by fgLasso. After a few tests, the second regularization parameter, which controls the similarity accross conditions, was set to the value 0.1 for all simulations;

• the consensus Lasso method, as described in Section 2. The two choices of consensus described in Section 2.3 are tested with, for a priori network:

the mother network (i.e., the true network used to generate the child networks, which is never known in practice but is the closest thing we have in this simulation from a bibliographic network),

or, for comparing a naive two-step approach with the averaged consensus described in Section 2.3.2, a network which is the mean over the conditions of independent estimations (i.e., estimations obtained with µ = 0).

These methods are denoted by cLasso (m) (for the averaged consensus method described in Section 2.3.2), cLasso (p) (for the method using the true prior) and cLasso (2) (for the naive two-step approach that uses a mean over conditions of independant estimations). Notice that the method using as a prior the mother network is clearly favored in this comparison, since even if the child networks are not identical to the mother network, they are very related to it. The comparison with this method should be used to understand what is the effect of integrating true prior knowledge in the estimation. µ was set equal to 1.

For each method, the inference is performed for a whole path of λ values and the corresponding precision and recall are calculated for each value of λ. A bootstrap version with 100 bootstrap samples of each of these methods is also implemented with T 1 = 500. The number of times a given edge in a given condition is selected is then used to calculate precision/recall values for different values of T 2 .

Precisions are recalls are calculated by comparing the estimated conditionspecific networks with the children networks they are generated from. We do not compare directly the consensus network with the mother network because we are interested in testing the ability of the method to estimate the common edges as well as the condition-specific edges. The F statistics is used as a way to compare the different methods, as in Section 4.1.2. First, averaged F , over the different conditions, are calculated along the precision/recall values obtained for different values of λ (direct estimation) or for different values of T 2 (bootstrap estimation).

Then, the maximum of these values (for recall and precision values larger than 0.05, to avoid extremely bad values of the precision or of the recall) is used as a way to compare the performance of the different methods. The results are given in Tables 2 (direct estimation) and 3 (bootstrap estimation), for each of the 6 methods described above (the best-case scenario for the bootstrapped cLasso corresponds to the parameters already listed in Table 1).

Method

gLasso iLasso groupLasso coopLasso fgLasso cLasso (m) cLasso (p) cLasso ( 2 Several conclusions can be drawn from these results. For a moderate ratio of rewired edges (smaller than 10%), bootstrapping improves the performances of all methods, except for iLasso (also, the increase is very limited for the coopLasso method). The increase is particularly interesting when the sample size is small and/or the number of samples is high. On the contrary, when the ratio of rewired edges is equal to 20%, bootstrapping only improves the performances of cLasso with prior, and, only for 4-5 samples having the smallest sizes, of gLasso, grou-pLasso and cLasso (m).

As expected, the overall performance is strongly increased when a relevant prior is added (the best F is often 3 times larger), which shows that this strategy should probably be used when such an information is available. When this is not the case, fgLasso often obtains the best results. Otherwise, coopLasso, bootstrap cLasso or iLasso also have good performances. Bootstrap cLasso (m) seems to be useful in the case of a moderate number of rewired edges and when the sample size is smaller. The naive two-step approach, which requires two estimations instead of one, often leads to deteriorated performances as compared to cLasso (m) which is twice faster. Finally, direct iLasso is advised for the largest number of rewired edges and coopLasso is to be preferred when the number of rewired edges is small but the sample size larger.

Computational times needed 6 for the different estimations are very different:

• the time needed to estimate one of the condition-specific networks with glasso (independent estimations with graphical LASSO as described in [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]) is approximately equal to 1 second for 25 values of λ (to be multiplied by the by the will to emphasize only a few different edges between the two conditions and hence to focus on the most relevant differences between the two conditions.

Notice that the possibility to monitor µ allows to infer networks that are more or less consensual, depending on what your prior is. The choice for T 1 was directed by the fact that we wanted to obtain very sparse networks, easily readable for the biologist, which, in the case of approximately 200 nodes, requires to have a very low density. 100 bootstrap samples where used to estimate the edges by the cLasso approach. The distribution of the number of times an edge is selected for expression data after the diet was very similar to what was found before the diet. Hence, in order to favor a high precision (at the cost of maybe a low recall), only the edges that appeared in at least 80 bootstrap samples were selected. This yielded to networks having respectively 316 and 315 edges (with a density about 1.3%). These networks had 292 edges in common (i.e., approximately 90% of the total number of inferred edges).

The histogram of the number of times a given edge is selected over the 100 bootstrap sample is given in Figure 6 (for the network corresponding to gene expression data before the diet). It has to be noted that pairs of variables that were never selected over the 100 bootstrap samples have been removed from the histogram (it corresponded to approximately 70% of the 221×220 potential edges). Other shared edges are highly probable, such as the one between AZGP1 and GPD1L which are two known biomarkers of the metabolic syndrome [START_REF] Viguerie | Determinants of human adipose tissue gene expression: impact of diet, sex, metabolic status and cis genetic regulation[END_REF]. However, quite interestingly, at least one condition specific edge is also expected: the genes PCK2 and CIDEA are the best biomarkers, among this set of genes, for the weight loss, and the main difference between the two conditions is indeed the weight loss [START_REF] Viguerie | Determinants of human adipose tissue gene expression: impact of diet, sex, metabolic status and cis genetic regulation[END_REF].

The resulting networks are displayed in Figure 7. They have approximately 92% of their edges in common. The full biological validation of such a network is unrealistic (because of the very limited knowledge available in this area) but some of the interactions make sense. For instance, some regulatory relationships shared between the two conditions are already known, such as the relation between FADS1 and FADS2, which encode two desaturase enzymes from the same cluster gene with similar regulation by dietary composition [START_REF] Mangravite | Fatty acid desaturase regulation in adipose tissue by dietary composition is independent of weight loss and is correlated with the plasma triacylglycerol response[END_REF].

Conclusion

We have proposed the cLasso method, which is used for jointly inferring networks in the case of multiple and dependent expression data. This method relies on the definition of a consensual solution, which in our case, is simply the mean between the different conditions. The different networks are forced toward this consensus by a L 2 -penalty whilst the sparsity of the solution is handled by an additional L 1 penalty. The solution proposed in this paper can be reformulated as a LASSO problem similar to the ones described in [START_REF] Meinshausen | High dimensional graphs and variable selection with the lasso[END_REF][START_REF] Chiquet | Inferring multiple graphical structures[END_REF] and the method is implemented in the R package therese. Experiments were conducted, using a bootstrapped approach based on the cLasso method and showed that this method is reliable.

Future work should address the issue of unbalanced sample sizes between conditions, and of the choice of µ: a naive selection based on out-of-bag MSE has been proven inefficient so far for selecting the best value. However, this parameter can also be useful for the biologist to include prior knowledge about how similar the condition-specific networks should be: using different µ provide a family of solutions with different fractions of common edges, among which the biologist is free to choose.

Algorithm 1 3 : 4 : 5 : 6 :) j,c 7 : 8 :

 1345678 Bootstrap cLasso 1: require: list of genes {1, . . . , p}; list of individuals {1, . . . , n} individuals' sample number c 1 , . . . , c n with c i ∈ {1, . . . , k} gene expressions X (dimension n × p) parameters µ ∈ R (L 2 -regularization parameter) and T 1 , T 2 ∈ N (number of edges selected) 2: initialize: ∀ c ∈ {1, . . . , k}, ∀ j, j ′ ∈ {1, . . . , p}, N c (j, j ′ ) ← 0 for b = 1 → P do Sample at random with replacement in {1, . . . , n} return bootstrap sample B b Use B b to solve Equation (6) or Equation (5) for a full set of λ values return (β c,λ,b j ) j,c,λ Find λ max := arg max λ j,j ′ ,c I β c,λ,b j =0 > T 1 return (β c,b j ) j,c := (β c,λmax,b j for all j, j ′ ∈ {1, . . . , p} do if β c,b j,j ′ = 0 then 9:

Figure 2 :

 2 Figure 2: Maximum F along the path of T 2 values for different parameters µ and T 1

  gives an indication on the influence of T 2 on the density of the inferred networks. The histogram displays the distribution of the number of times a given edge is chosen by the algorithm over the 100 bootstrap samples. Notice that only a few edges are very frequently selected by the bootstrap method, whereas the targeted density of 4% is obtained by keeping edges that are selected only about 40 times (i.e., 40% of the bootstrap samples).

Figure 3 : 1 Figure 4

 314 Figure 3: Distribution of the number of times an edge is selected over 100 bootstrap samples for the first condition (left) and Evolution of the density versus T 2 for the 2 conditions (right). T 1 = 500, µ = 1

Figure 4 :

 4 Figure 4: Precision/Recall curve (with varying T 2 ) for the total number of edges of the true child networks compared to the corresponding inferred network. T 1 = 500, µ = 1 (each curve correspond to one of the k = 2 child networks.

Figure 5 :

 5 Figure 5: Evolution of the density versus T 2 for the 2 conditions (left) and Precision/Recall curve (with varying T 2 ) for the total number of edges of the true child networks compared to the corresponding inferred network. T 1 = 250, µ = 1

Figure 6 :

 6 Figure 6: Distribution of the number of times an edge is selected over 100 bootstrap samples for the expression data before the diet.

Figure 7 :

 7 Figure 7: Networks inferred by bootstrapped cLasso: before (left) and after lowcalorie diet (right). Blue edges are shared edges and pink edges are condition specific edges. Some gene names are given, that are commented in the text.

  1 ij ) i=1,...,n 1 , j=1,...,p . . . (X k ij ) i=1,...,n k , j=1,...,p , with

c n c = n and, for all c and all i, (X c ij ) j=1,...,p are i.i.d. observations of N (0, Σ c ).

Table 2 :

 2 Summary of the performance for the different method in terms of the maximum value of the F statistics. The best method for each couple of percentage of rewired edges and network sizes, is emphasized with bold face.

	)

Table 3 :

 3 Summary of the performance for the different method (bootstrap version) in terms of the maximum value of the F statistics. The best method for each couple of percentage of rewired edges and network sizes, is emphasized with bold face.

therese can be downloaded on R-Forge, from http://therese-pkg.r-forge.r-project. org/.

http://cran.r-project.org/web/packages/glasso/index.html

http://stat.genopole.cnrs.fr/logiciels/simone

http://cran.r-project.org/web/packages/JGL/index.html

http://www.comp-sys-bio.org/AGN/Century/index.html

Computational times are reported for a 4-core laptop, Intel(R) Core(TM) i5-3360M CPU @ 2.80GHz, RAM 8Go DDR3, with OS Kubuntu Linux 12.04.

Acknowledgements

We thank the two anonymous reviewers for their useful comments and suggestions.

We also thank Jean-Michel Poggi for thorough remarks that gave useful hints to use the method for integrating prior knowledge and Nicolas A. Edwards who helped testing several joint inference methods during his internship.

number of conditions);

• the time needed to estimate 5 joined networks with simone (implementing the methods described in [START_REF] Chiquet | Inferring multiple graphical structures[END_REF]; the time is reported for "cooperative LASSO") is approximately equal to 1 minute 30 seconds for 100 values of λ;

• the time needed to estimate 2 joined network with JGL (implementing the methods described in [START_REF] Danaher | The joint graphical lasso for inverse covariance estimation accross multiple classes[END_REF]; the time is reported for "fused LASSO") is approximately equal to for 25 values of λ and the time needed to estimate 5 joined networks with this method is approximately equal to 2 hours 30 minutes. Notice that the path of λ has been performed manually as this package is the only one that does not propose a regularization path for the sparse parameter;

• finally, the time needed to estimate 5 joined networks with therese is approximately equal to 2 minutes 30 seconds for 100 values of λ (and a little bit less than twice this value for the naive two-step approach).

Real data: effect of a diet on regulatory network

As an application to a real biological data set, we analyzed gene expression data described in [START_REF] Viguerie | Determinants of human adipose tissue gene expression: impact of diet, sex, metabolic status and cis genetic regulation[END_REF]. More precisely, the expression of 221 genes are used. These were obtained for 204 obese women before and after a 8-week low-calorie die (LCD) with the objective of more than 8% weight loss. Considering the two time steps of the analysis as two different conditions, we used the bootstrapped cLasso approach with µ = 1 and T 1 = 1000 (that corresponded to a targeted density of approximately 4%). The choice of a rather high regularization parameter µ was directed

Authors' biography