Noureddine Lebaz
email: lebaz@insa-toulouse.fr

Arnaud Cockx

Mathieu Spérandio

Jérôme Morchain

Reconstruction of a distribution from a finite number of its moments: a comparative study in the case of depolymerization process

Keywords: Moment problem, population balance, reconstruction methods, maximum entropy

The resolution of the Population Balance Equation (PBE) using moment-based methods offers a high computational efficiency however, information on the time evolution of the probability density function (PDF) is out of reach. For this, several PDF reconstruction methods using a finite number of moments are proposed in the literature. In this contribution, three different methods (i.e. Beta Kernel Density Function based method, Spline based technique and the Maximum Entropy based approach) are tested and compared to the analytical solution of a depolymerization process. The Maximum Entropy method gives the most accurate approximations using only a set of six moments. This method is combined with the Quadrature Method of Moments (QMOM) for a simultaneous reconstruction during the PBE resolution. A three nodes and a four nodes quadrature are tested. The results show that the quality of the reconstruction is highly dependent on the accuracy of the computed moments.

Introduction

The recovery of a probability density function (PDF) knowing only a finite number of its moments is known as the finite-moment problem in mathematical analysis and arises in different scientific applications (e.g. physics, chemical engineering, economics) [START_REF] Gavriliadis | The truncated Stieltjes moment problem solved by using kernel density functions[END_REF][START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF]. This problem is generally declined in three different problem categories for the mono-variate case [START_REF] Abramov | An improved algorithm for the multidimensional moment-constrained maximum entropy problem[END_REF]):

• The Hausdorff moment problem: the PDF is supported on the closed interval [a, b] • The Stieltjes moment problem: the PDF is supported on [0, +∞)

• The Hamburger moment problem: the PDF is supported on the real line (-∞, +∞)

In chemical engineering, especially for particulate/dispersed systems (e.g. crystallization, polymerization/depolymerization, liquid-liquid extraction, multiphase systems), population balance models (PBM) are widely used for the description of the time evolution of the variable-based distribution (e.g. size, volume) undergoing elementary processes. Among these processes, one commonly finds breakage, aggregation/coalescence, nucleation, growth/dissolution . . . etc [START_REF] Ramkrishna | Population balance modeling. Promise for the future[END_REF]. The resolution of the population balance equation (PBE) is computationally intensive when using classical approaches (e.g. Monte Carlo methods [START_REF] Lin | Solution of the population balance equation using constantnumber Monte Carlo[END_REF], discretization methods [START_REF] Kumar | On the solution of population balance equations by discretizationI. A fixed pivot technique[END_REF]). This drawback is limiting when population balance modelling is coupled with Computational Fluid Dynamics (CFD) which is the case for multiphase systems. In order to overcome this issue, moment based methods (MOM:

Standard Method of Moments [START_REF] Hulburt | Some problems in particle technology: A statistical mechanical formulation[END_REF], QMOM: Quadrature Method of Moments [START_REF] Mcgraw | Description of Aerosol Dynamics by the Quadrature Method of Moments[END_REF]Marchisio et al., 2003a), DQMOM: Direct Quadrature Method of Moments [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF]) have been developed. The problem is reduced to the time-tracking of a finite number of the PDF moments offering by the fact a high computational efficiency. In QMOM the transported moments are calculated by reducing the PDF to an n-point distribution (sum of n weighted Dirac delta functions), the corresponding weights and abscissas are computed using specific algorithms [START_REF] John | On the efficiency and robustness of the core routine of the quadrature method of moments (QMOM)[END_REF]. In DQMOM, the weights and the abscissas of the initial n-point distribution are directly tracked instead of the moments.

Even though the moment-based methods are computationally efficient, information regarding either the shape or pointwise values of the PDF are identically out of reach. However, this level of knowledge is important for different applications as in the case of evaporation process where one needs to accurately evaluate the PDF at zero size [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] or in the enzymatic hydrolysis of particulate substrates where the loss of small solubilized particles produced by the enzymatic attacks has to be taken into account [START_REF] Lebaz | Population balance approach for the modelling of enzymatic hydrolysis of cellulose[END_REF]. To address this issue, different numerical techniques have been proposed for the reconstruction of the PDF from its moments reviewed by [START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF], mainly for the Hausdorff moment problem since in chemical engineering applications, the support of the PDF is known in most cases.

From a mathematical point of view, the moment problem has been extensively studied focusing on the conditions of existence of a unique or infinite solution(s) [START_REF] Shohat | The Problem of Moments[END_REF][START_REF] Akhiezer | The Classical Moment Problem: N.I. Akhiezer[END_REF][START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF]. Theoretically, a perfect reconstruction can be obtained using an infinity of PDF moments with an a priori restriction of the class of basis functions [START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF]. Numerically, this problem is known as a difficult inverse problem since the finite number of moments define a high ill-posed system of equations [START_REF] Athanassoulis | The truncated Hausdorff moment problem solved by using kernel density functions[END_REF]. Thus, there is no absolute method for reconstructing accurately the PDF from a finite number of its moments. From an experimental point of view the number of moments that can be obtained with a reasonable accuracy is limited.

The most intuitive reconstruction technique is to approximate the target PDF by a sum of elementary distributions (e.g. Gaussian, log-normal) [START_REF] Lee | Change of particle size distribution during Brownian coagulation[END_REF][START_REF] Diemer | A moment methodology for coagulation and breakage problems: Part 2moment models and distribution reconstruction[END_REF] when an information on the shape of the PDF is available a priori. Thus, the problem is reduced to a simple parameter fitting. This technique offers some advantages like the reduced number of moments required for the reconstruction, its simplicity and fastness but the shape of the target PDF has to be stipulated. This constitutes its principal limitation when the initial shape of the PDF changes dramatically during the course of the process. Inspired by this approach, [START_REF] Athanassoulis | The truncated Hausdorff moment problem solved by using kernel density functions[END_REF] proposed a robust technique based on the PDF approximation by a finite superposition of kernel density functions (KDF) for the Hausdorff moment problem, tested and validated successfully against both monomodal and bimodal PDFs. This approach has been extended for the Stieltjes moment problem using a generalized Gamma function as the kernel density function [START_REF] Gavriliadis | The truncated Stieltjes moment problem solved by using kernel density functions[END_REF][START_REF] Gavriliadis | Moment information for probability distributions, without solving the moment problem, II: Main-mass, tails and shape approximation[END_REF] but unfortunately it was not confronted to realistic cases. [START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF] proposed an innovative technique using a piecewise polynomial function without any a priori assumption about the shape of the target PDF. Linear, quadratic or cubic splines can be used, their coefficients are computed by solving an ill-conditioned linear system of equations. De [START_REF] De Souza | Reconstruction of a distribution from a finite number of moments with an adaptive spline-based algorithm[END_REF] improved the spline-based method by developing an adaptive algorithm in order to optimize the distribution of the grid nodes and capture more accurately the PDF in critical domains. This technique has been coupled recently with a PBE describing aggregation of urea particles [START_REF] Hackbusch | A numerical method for the simulation of an aggregation-driven population balance system[END_REF], droplet coalescence [START_REF] Bordás | Numerical methods for the simulation of a coalescence-driven droplet size distribution[END_REF] and pharmaceutical drying process (granules) [START_REF] Mortier | Comparison of techniques for reconstruction of a distribution from moments in the context of a pharmaceutical drying process[END_REF].

The third class of methods that will be considered in this paper is the Maximum entropy (ME) method which received increasing interest in the last two decades [START_REF] Mead | Maximum entropy in the problem of moments[END_REF][START_REF] Tagliani | Hausdorff moment problem and maximum entropy: a unified approach[END_REF][START_REF] Tagliani | Numerical aspects of finite Hausdorff moment problem by maximum entropy approach[END_REF][START_REF] Biswas | Function reconstruction as a classical moment problem: a maximum entropy approach[END_REF]. This technique is based on the maximization of the Shannon entropy from information theory, by solving a constraint optimization problem. [START_REF] Abramov | A practical computational framework for the multidimensional moment-constrained maximum entropy principle[END_REF] extended it to multidimensional problems. [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] coupled this technique with DQMOM in the case of droplet evaporation process.

In this paper we propose firstly an assessment of these methods in terms of accuracy, rapidity and minimum number of moments required for the reconstruction. The issue concerns the identification of the most efficient method to perform the reconstruction throughout the process (not only for a particular distribution). The performances of the three methods will be evaluated through a comparison of the reconstructed PDFs with an analytical solution of the PBE accounting for a breakage process. Finally, the implementation of the relevant method in the QMOM code for a simultaneous resolution of the PBE and reconstruction of the PDF is investigated. The outlines of this contribution are schematically summerized in figure 1.

Reconstruction methods

The purpose of the reconstrution techniques is to recover a function f (x), given its integer moments sequence µ n defined as:

µ n = x n f (x)dx, n = 0, 1, . . . (1)
We give hereafter a succinct description of the reconstruction techniques used in this contribution. For more detail, one can refer to the original works cited and/or to Appendix A. As common point, the bounded support [a,b] is rescaled to [0,1] for a general formulation. In the KDF based technique, the target PDF is approximated by a sum of weighted kernel density functions (equation 2) [START_REF] Athanassoulis | The truncated Hausdorff moment problem solved by using kernel density functions[END_REF]:

80 f (x) = I i=1 p i K(x; x i , h) (2)
where K(x; x i , h) are the KDFs, centered at x i , with bandwidth h, I is the total number of kernel density functions used for the reconstruction. The coefficients p i satisfy:

p i ≥ 0, i = 1, 2, . . . , I; I i=1 p i = 1 (3)
Since the target PDF is univariate, the KDF is chosen as a Beta Kernel. Thus, the objective is the determination of the optimal coefficients p i . For this, the finite-moment problem is reformulated as a constrained optimization problem aiming to find the coefficients p i which minimize the error between the set of initial 85 moments and those estimated via the sum of Beta KDFs. Numerically, Nonnegative Least Square algorithm (NNLS) can be used in this case. Optimization techniques are to be supplemented in order to improve the efficiency of this procedure. The properties of the Beta Kernel Function, the formulation of the constrained optimization problem, the optimization technique based on the shifted moments and the general algorithm of this method are described in Appendix A.1.

Spline-based method

This method is extensively described by [START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF]. The support of the target PDF [a, b] is

subdivided into n subintervals such as a = x 1 < x 2 < • • • < x n+1 = b. In each subinterval [x i , x i+1 ],
the target PDF is approximated by a piecewise polynomial s (l) (x) of degree l. The system of equations is detailed for cubic splines (l = 3). The unknowns are the four coefficients of the n splines. A smooth transition at the boundaries of the interval is assumed meaning that the PDF, its first and second derivatives are null at the boundaries, this gives 2 × 3 equations. The continuity of the splines, their first and second derivatives at the nodes provides 3(n -1) equations. The remaining (n -3) equations are supplemented by the moments. This leads to solve a 4n × 4n ill-conditioned linear system.

Since the set of the known moments limits the number of splines, one has to compute the reconstruction in the optimal support of f (x) thus, the reconstruction is computed iteratively starting from an initial

reconstruction f (0) (x) in an initial interval [x (0) 1 , x (0)
n+1 ]. To ensure the positivity of the reconstructed PDF a tolerance tol neg is introduced. Furthermore, the optimization of the distribution of the nodes improves the efficiency of the method especially for multimodal functions. More details on these techniques are given in Appendix A.2

Maximum Entropy method

The Maximum Entropy method is based on the maximization of the Shannon entropy H[f ] from information theory given by:

H[f ] = - ∞ 0 f (x)lnf (x)dx (4)
Under the moments constraints defined by equation 1.

The explicit representation of the ME approximation f M (x) of the target PDF takes the form:

f M (x) = exp   - N j=0 ξ j x j   (5)
To be supplemented by the (N + 1) constraints:

µ i = ∞ 0 x i f M (x)dx i = 0, 1, . . . , N (6)
The (N + 1) Lagrange's multipliers ξ 0 , . . . , ξ N are obtained through the resolution of the following set of (N + 1) nonlinear equations:

∞ 0 x i exp   - N j=0 ξ j x j   dx = µ i i = n, 1, . . . , N (7)
Numerically, iterative methods are used to solve equation 7 [START_REF] Mead | Maximum entropy in the problem of moments[END_REF]. In our case, we use the standard Newton method starting from an initial choice of the Lagrange's multipliers ξ (0) with a given tolerance for stopping the iterative procedure. Commonly used parameters are given in the results section.

Population balance modelling (PBM)

The mathematical modelling of polymer chain fragmentation processes (thermal, thermochemical and biological) is a fundamental issue in polymer science and engineering [START_REF] Mccoy | Discrete and continuous models for polymerization and depolymerization[END_REF][START_REF] Bose | Mathematical modelling and computer simulation of linear polymer degradation: Simple scissions[END_REF].

Population balance modelling (PBM) is the classical approach since it describes the dynamical evolution of the Chain Length Distributions (CLD) (or MWD: Molecular Weight Distributions) during depolymerization reactions [START_REF] Ziff | The kinetics of cluster fragmentation and depolymerisation[END_REF][START_REF] Madras | Time evolution to similarity solutions for polymer degradation[END_REF].

For a homogeneous system, the time evolution of the CLD undergoing fragmentation processes is described by the population balance equation (PBE) below, written in its continuous form:

∂n(L, t) ∂t = ∞ L β(L, L )γ(L )n(L , t)dL -γ(L)n(L, t) (8)
where n(L, t) is the length-based number density function (CLD), γ(L) the breakage frequency for a chain of length L, β(L, L ) is the breakage kernel giving the probability of obtaining a chain of length L from the breakup of a chain of length L . The first term on the RHS accounts for the formation (birth) of chains with length L resulting from the breakage of longer chains L . The last term is the death term due to the loss of chains of length L because of their breakup.

Except for some cases where analytical solutions exist, equation 8 is solved numerically. Since the resolution of this equation is not the objective of this contribution, an analytical solution is briefly described and used as reference for the moment problem.

Analytical solution of the PBE accounting for breakage process

Ziff [START_REF] Ziff | The kinetics of cluster fragmentation and depolymerisation[END_REF] considered linear polymers of length L and gave the analytical solution of the equation 8 with the breakage kernel β = 2/L and the breakage frequency γ = αL 2 (α is a factor of proportionality) which can be regarded as the number of attacks per time unit. Under these considerations, equation 8 can be rewritten as:

∂n(L, t) ∂t = 2α ∞ L L n(L , t)dL -αL 2 n(L, t) (9)
The general solution of the equation 9 is given as [START_REF] Ziff | The kinetics of cluster fragmentation and depolymerisation[END_REF]:

n(L, t) = e -αtL 2 n(L, 0) + 2αt ∞ L L n(L , 0)dL (10)
where n(L, 0) is the initial CLD.

In order to explicit this solution, we assume that the initial chain length distribution (CLD) follows a normal law (equation 11) with (m, σ) its mean and standard deviation respectively.

n(L, 0)) = 1 σ √ 2π e -(L-m) 2 2σ 2 (11)
The explicit solution in this specific case is given as (details are given in Appendix B) :

n(L, t) = 1 σ √ 2π e -αtL 2 (2αtσ 2 + 1)e -1 2 ( L-m σ ) 2 + mσαt √ 2π 1 -erf 1 √ 2 L -m σ ( 12
)
where erf refers to the error function.

Thus, starting with a CLD following a normal law with specific fragmentation process, one can obtain the exact CLD and compute its moments at any given time during the degradation reaction.

The Monovariate Quadrature Method of Moment (QMOM)

The PBE's resolution using Moment Methods is based on the time tracking of a finite set of the CLD moments instead of the CLD itself. This class of methods is widely used for its low computational cost thus, it can be implemented in CFD codes or extended to multidimensional problems. A succinct description of the well-known Quadrature Method of Moments (QMOM) proposed by [START_REF] Mcgraw | Description of Aerosol Dynamics by the Quadrature Method of Moments[END_REF] and validated by Marchisio et al. (2003b) is given in the specific case of breakage processes.

By a moment transformation of equation 8, the transport equation for the k th moment is:

∂µ k (L, t) ∂t = ∞ 0 L k ∞ 0 γ(L )β(L, L )n(L , t)dL dL - ∞ 0 L k γ(L)n(L, t)dL (13)
The Quadrature Method of Moments (QMOM) is based on the Gaussian quadrature of the continuous CLD [START_REF] Mcgraw | Description of Aerosol Dynamics by the Quadrature Method of Moments[END_REF]:

n(L, t) ≈ M i=1 w i (t)δ(L -L i (t)) ( 14
)
where M is the number of nodes i, L i is the property of the node (length), w i its weight, and δ is the Dirac function. Thus, the k th moment can be expressed as :

µ k (t) = ∞ 0 n(L)L k dL ≈ M i=1 w i (t)L i (t) (15)
where weights (w i ) and abscissas (L i ) are determined through the Product-Difference algorithm [START_REF] Gordon | Error Bounds in Equilibrium Statistical Mechanics[END_REF]. By substituting equations 14 and 15 in 13, we obtain (Marchisio et al., 2003b) :

∂µ k ∂t = M i=1 γ i b (k) i w i - M i=1 L i γ i w i (16) with b (k) i = ∞ 0 L k β(L, L i )dL (17)
Generaly, the QMOM requires at least a three nodes quadrature (M = 3) for an accurate time-tracking of the moments (Marchisio et al., 2003b).

MATLAB is used under a work-station comprising of an Intel R Core TM i7-3740QM CPU with a clock speed of 2.7 GHz and 16 GB of installed memory (RAM).

Results and discussion

The reasons for searching efficient reconstruction techniques when one deal with solving PBM are multiple. First of all, Methods of Moments are rather powerful in terms of numerical performance but since these are integral methods the connection with the physical distribution (Chain Length Distribution or Molecular Weight Distribution) is lost. Moreover, high order moments of an experimental distribution are often hard to evaluate and the direct comparison of numerical and experimental moments remains difficult to analyze.

In some cases, the distribution itself is necessary to evaluate some terms in the equations. To finish with, experimental and numerical distributions are not always known over the same range: the range of experimental methods is limited by lower and upper detection bounds. For all these reasons, an accurate reconstruction of a PDF from its moments is an essential tool for those involved in the analysis of experimental results with the help of mathematical models such as PBM solved with a method of moments.

Reconstruction of the analytical solution

As described before, the three reconstruction methods will be compared to the analytical solution given in section 3-1. We assume that the initial CLD is given by a normal law centered at a Degree of Polymerization (DP) of 60 (the mean number of monomers constituting a given polymer) with a standard deviation (σ) of 7. Equation 12gives the time evolution of the CLD for (β = 2/L and γ = αL 2 ). This solution is illustrated in figure 2 with α = 10 -3 and t

∈ [t i , t f ] = [0, 1].
Figure 2 shows a rapid degradation of the long polymers since the breakage frequency is proportional to L 2 , the concentration of the small polymer chains increases all along the depolymerization process.

For sake of clarity, the ability of the different methods to reconstruct the time evolution of the CLD is described separately and illustrated for five different instants:

t 0 = t i , t 1 = 1 4 t f , t 2 = 1 2 t f , t 3 = 3 4 t f and
t 4 = t f before addressing a comparative analysis. At any given instant, the reconstruction is initiated by computing the set of the moments of the analytical CLD given in equation ( 12). The integrations leading to the associated moments, equation (1), are performed using the Matlab function quad with a tolerance set to 10 -6 . Since the number of the moments required for the reconstruction is critical, the methods are tested with the same set of moments transformed systematically to the normalized interval [0,1] before rescaling the reconstruction result (figure 1). In order to compare the different reconstruction methods, we introduce the mean absolute error between the reconstructed CLD (f M (x)) and the analytical solution used as a reference (f (x)):

Er(t) = 1 q q i=1 |f (x, t) -f M (x, t)| ( 18
)
where q is the length of the discretized support.

Kernel Density Function-based method

The algorithm of this method was implemented in Matlab using the set of moments of each target CLD with their associated shifted moments by setting x = 0.8 in order to increase the number of the constraints and thus improve the efficiency of this technique. The number of Beta Kernel Density Functions (BKDFs) I is set to 15 thus, the support [0,1] is divided into 15 equal intervals, each BKDF is centered on the midpoint of its corresponding interval. [START_REF] Athanassoulis | The truncated Hausdorff moment problem solved by using kernel density functions[END_REF] showed that the optimal bandwidth parameter h varies in the interval [0.04,0.09]. In order to identify the optimal values for h, the interval is scanned and the reconstruction is conducted at different parameter values before selecting the most accurate. Once the BKDFs parameters are fixed, the associated Beta PDF parameters are calculated and the coefficients B n,i

and Bn,i (equations A.3 and A.8) are computed. The constraint minimization problem (equation A.10) is solved using LSQLIN subroutine.

Figure 3 shows the results for the five different instants using the first six moments of each target CLD.

Although the initial CLD is reconstructed with great accuracy most probably because of its simple form, the other CLDs are not recovered. Their singular shape mainly around the origin makes their reconstruction challenging. The reconstructions are oscillating around the target CLDs and take the value zero at the origin. To be more accurate, this technique requires a large number of constraints which means a higher number 210 of moments. In figure 4, the same problem is treated using the first eleven moments of each target CLD. The figure 4 shows that the use of a high number of moments can improve the reconstruction quality as illustrated in Table 1 where the mean absolute errors (equation 18) are calculated in both cases. But, since the technique is treating a highly ill-posed problem, the solution of the constraint system of equations may diverge. This is for example the case when the reconstruction of the initial CLD using eleven moments is attempted (figure 4). To sum up, this method requires a finetuning of the numerical parameters and an increase in the number of moments does not always guaranties an easier reconstruction. However if the reconstruction is successful, the accuracy increases with the number of moments.

Spline-based method

The adaptive algorithm with cubic splines proposed by De Souza et al. ( 2010) is implemented and used in this contribution. The Matlab script for the original algorithm was freely downloaded from (http://www.unimagdeburg.de/isut/LSS). The default tolerances proposed by the authors (tol red and tol sing ) are kept unchanged. This procedure is adopted deliberately because when the reconstruction is conducted simultaneously with the PBE resolution, one has to set the parameters of the algorithm a priori. The authors recommended to set tol neg at a very small negative value instead of 0. In our case, tol neg = -10 -1 is used initially. The results of the reconstruction using the spline based method with six moments are shown in figure 5.

Only the initial CLD has been recovered using the parameters described before. For the other instants, the algorithm has been stopped after 10 3 iterations. The modification of the negative tolerance to tol neg = -5.10 -3 leads to the convergence of the algorithm for t 1 but the reconstruction is highly oscillating. As mentioned by [START_REF] Mortier | Comparison of techniques for reconstruction of a distribution from moments in the context of a pharmaceutical drying process[END_REF], the algorithm is highly sensitive to the negative tolerance leading to different reconstructions for the same target PDF as illustrated in figure 6. This is a main drawback since tol neg may affect the algorithm and even its rate of convergence. In fact, for the example given in figure 6, the convergence is reached after 27 iterations for the case where tol neg = -10 -1 when the algorithm needs 42 iterations for the second case (tol neg = -10 -2 ).

Furthermore, a smooth transition at the boundaries of the recontruction interval was assumed. This is the case only for the initial distribution thus, one has to change the boundary conditions during the process. This is not reachable without an a priori information on the shape of the target PDF. The non-linear system in equation 7 is solved using Newton method starting from an initial choice of the Lagrange's multipliers as ξ (0) = (-ln(µ 0 )/(x max -x min ), 0, . . . , 0) with a tolerance of 10

). The moments are scaled to [0,1] and normalized before such as µ 0 = 1. The results are given in figure 7. Unlike the two previous methods, the Maximum Entropy technique reconstructs the analytical CLDs with high accuracy even when they are singular using only the first six moments. The Newton method used in this case converges in a few iterations (< 10) making the coupling of the PBE resolution with the PDF reconstruction possible. In addition, except for the tolerance of the Newton method, there is no additional parameter to set which simplifies its use. Table 2 gives the comparison between the three methods according to the quality of the reconstructed CLDs (mean absolute error) and the computational cost using the first six moments. The ME method gives the best results for this case of study not only for the reconstruction efficiency but also for the number of moments required and the computational cost among other criterions as resumed in Table 3. For these reasons, this method is selected to be implemented in the QMOM code for coupling the PBE resolution in the case of depolymerization process with the CLD reconstruction.

Criterion

QMOM coupled with ME based method

The main objective of the development of reconstruction methods from a finite set of moments, especially in the chemical engineering field, is to have access to the shape of the PSD and/or some relevant pointwise values during the process evolution. This kind of information is necessary when some aspects of the problem cannot be expressed from the moments only. This offers a valuable addition when the PBE is solved using moment methods coupling by the fact the computational efficiency of such methods with a simultaneous information on the PDF shape.

In this section, the ME based method tested successfully against an analytical solution giving the time evolution of the CLD undergoing breakage process is implemented in the QMOM code resolving the same breakage problem. Thus, on the one hand, the analytical CLDs are available for comparison, on the other hand, the PBE accounting for the same breakage process is solved using QMOM giving the time evolution of the moments which are used for the reconstruction (see figure 1) .

For this numerical approach, a three nodes quadrature is used (M=3), the abscissas (L i ) and the weights (w i ) are computed using the Product-Difference algorithm [START_REF] Gordon | Error Bounds in Equilibrium Statistical Mechanics[END_REF]. The system of ordinary differential equations (equation 16) is integrated using ode45 with a tolerance fixed to 10 -6 . First, the QMOM implementation is validated by confronting the computed moments against those predicted by the analytical solution. The results for the six first moments are given in figure 8. The two different methods lead to the same time evolution of the moments thus, the QMOM implementation is validated.

Since a three nodes quadrature is used, six moments are tracked. This set of moments is used as input for the ME based method. The reconstruction is conducted simultaneously with the PBE resolution. The result at different process time is shown in figure 9.

The reconstructions based on the numerical moments computed via QMOM (figure 9) are less accurate than those obtained previously with the exact analytical moments (figure 7). Note that the same number of moments is used in the two cases. This can only comes from the numerical error introduced by QMOM.

In order to improve the accuracy of the method, one can increase the number of the quadrature nodes. For quantifying the gain in terms of accuracy, we introduce the global error function below: We give in figure 10 the error E(k) induced by QMOM for both a three and four nodes quadrature for the six first moments. Note that the use of a four nodes quadrature implies a time tracking of the first eight moments of the CLD.

E(k) = |µ k -µ * k | µ * k (19)
The use of a four nodes quadrature in QMOM improves the accuracy of the computed moments, which consitutes a trivial and expected result. Note however that the use of a higher number of nodes induces supplementary moments meaning more equations in the system (equation 16) and may lead to near-singular matrix limiting the efficiency of PD algorithm. Now, if the first six moments computed using a four nodes quadrature are used for the reconstruction, the result is completely different as shown in figure 11 andin Table 4.

Improving the accuracy of the tracked moments has considerably improved the reconstruction quality without reaching the same performance as with the analytical solution (figure 7). This shows that, if the ME based method is coupled with QMOM, one has to compute the moments with a high accuracy. For the reconstruction purpose, it is wise to track some additional moments in order to improve the accuracy of the set of moments that will be used.

To avoid the use of the PD algorithm for computing the weights (w i ) and the abscissas (L i ) all along the process which may introduce a subsequent error for the computed moments, the Direct Quadrature Method of Moments (DQMOM) presents an alternative [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF]. In this method, the PD algorithm is only called once, for initializing the quadrature nodes. This method has been investigated and comparative results as with QMOM are obtained. The error on computing the initial (L i , w i ) is transported throughout the process time. In addition, it has been shown that more accurate and robuste algorithms can be used instead of the PD algorithm for computing the nodes [START_REF] John | On the efficiency and robustness of the core routine of the quadrature method of moments (QMOM)[END_REF]. The use of the Wheeler algorithm (called also LQMDA for Long Quotient-Modified Difference Algorithm or Chebyshev algorithm) [START_REF] Wheeler | Modified moments and Gaussian quadratures[END_REF], recommended for computing optimal weights and abscissas [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF][START_REF] John | On the efficiency and robustness of the core routine of the quadrature method of moments (QMOM)[END_REF], leads to the same results as those shown before. So in the case studied here, combining DQMOM with the Wheeler algorithm does not seem to improve the accuracy of the reconstruction. The advantage of this algorithm is to avoid the limitation of the PD algorithm when near-singular matrix arise.

Conclusion

The reconstruction of a full PDF knowing only a finite/small number of its moments is of critical importance especially for particulate systems where information on the PDF shape is needed for process follow-up.

Different reconstruction methods are available in literature but none of them is able to recover all the possible PDF shapes perfectly.

In this contribution, three different methods are investigated and their results confronted firstly against an analytical solution in the specific case of depolymerization process. The BKDF based method approximates the target PDFs by a finite sum of elementary Beta Kernel Density Functions and gives relative accurate results but requires a high number of moments for complex PDF shapes. The spline based method is highly dependent on the tolerance parameters. Finetuning of these parameters affects the rate of convergence of the algorithm and the final reconstructed PDF and depends on each single target PDF. Furthermore, the transition at the boundaries of the reconstruction interval has to be assumed which is not reachable when no information about the shape of the target PDF is available. As consequence, this method is not adapted for process monitoring when the distribution is needed at different time steps. The ME based method gives the best results and reconstructs the analytical target PDFs with high accuracy using only a set of six moments.

The method converges within a few iterations which makes it relevant to be coupled with the PBE resolution methods. This has been investigated in the second part.

Once the Quadrature Method of Moments accounting for breakage processes is validated against the analytical solution, the computed moments are used as inputs for the reconstruction method. The results

In order to optimize the reconstruction support and thus increase the accuracy of the method, the values of the reconstructed PDF f (k) (x) are checked at the boundaries of the support [x

(k) 1 , x (k) n+1 ] (i.e. in the subintervals [x (k) 1 , x (k) 2 ] and [x (k) n , x (k) n+1 ]
). For this, these two subintervals are divided into 10 equidistant small subintervalls, the values of f (k) (x) at the nodes of these small subintervals are checked if they are sufficiently small against the maximal value of f (k) in order to reduce eventually the size of the reconstruction support. This test is given for the right boundary as [START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF]: It is further checked that any value of f (k) (x) at the nodes x i (i = 1, . . . , n + 1) and at the midpoints of the subintervals are greater than a given tolerance tol neg (equation A.14). This test is implemented to ensure the positivity of the PDFs.

  10 j=1 f (k) (x nj ) 2 + f (k) (x (k) n+1 ) 2   1/2 ≤ tol red f (k) max (A.
f min,max := min j=1,...,k f (j) max , f (j) (x i ) f min,max ≥ tol neg , f (j) (x i6 ) f min,max ≥ tol neg (A.14)
with tol neg ≤ 0. The problem of the ill-conditioning of the linear system of equations is treated by regulizing the system. This is achieved by neglecting small singular values according to a given tolerance tol sing using a pseudo inverse routine of Matlab.

Note that in [START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF], an equidistant grid is used for the distribution of the nodes x i . This is not adapted for non-smooth and/or multimodal distributions. In the adaptive algorithm by [START_REF] De Souza | Reconstruction of a distribution from a finite number of moments with an adaptive spline-based algorithm[END_REF], this problem is solved by introducing a non-equidistant grid where the nodes are repositionned appropriately in order to capture accurately the PDF's critical domains.

Appendix B: Analytical solution of the PBE in the case of breakage process

We give in this appendix the details for the analytical solution in its explicit form (equation 12).

By substituting equation 11 in 10, we obtain :

n(L, t) = 1 σ √ 2π e -tL 2 e -1 2 ( L-µ σ ) 2 + 2t ∞ L L e -1 2 L -µ σ 2 dL (B.1)
We note I the integral term and proceed to the change of variable :

L -µ = L -→ L = L + µ -→ dL = d L (B.2)
The integral term is written as:

I = L e -1

Figure 1 :

 1 Figure 1: Schematic representation of the organigram of the contribution

Figure 2 :

 2 Figure2: Time evolution of the CLD given by the analytical solution of[START_REF] Ziff | The kinetics of cluster fragmentation and depolymerisation[END_REF]

Figure 3 :

 3 Figure 3: Comparison between the target CLDs (continuous line) and their reconstructions (dashed line) using the BKDF based technique with a sequence of six analytical moments

Figure 4 :

 4 Figure 4: Comparison between the target CLDs (continuous line) and their reconstructions (dashed line) using the BKDF based technique with a sequence of eleven analytical moments

Figure 5 :

 5 Figure 5: Comparison between the target CLDs (continuous line) and their reconstructions (dashed line) using the spline based technique with a sequence of six analytical moments

Figure 7 :

 7 Figure 7: Comparison between the target CLDs (continuous line) and their reconstructions (dashed line) using the ME based technique with a sequence of six analytical moments

  General comparison between the three different reconstruction methods

Figure 8 :

 8 Figure 8: The time evolution of the six first moments: comparison between the analytical solution (continuous line) and QMOM with three nodes (dashed lines)

Figure 9 :

 9 Figure 9: Comparison between the target CLDs (continuous line) and their reconstructions (dashed line) using the ME based technique with a sequence of six moments computed via QMOM with a three nodes quadrature

Figure 10 :

 10 Figure 10: Comparison of the error induced by QMOM for the six first moments when using a three nodes quadrature (blue line) and a four nodes quadrature (red line)

Figure 11 :

 11 Figure 11: Comparison between the target CLDs (continuous line) and their reconstructions (dashed line) using the ME based technique with a sequence of six moments computed via QMOM with a four nodes quadrature

  13)If equation A.13 is satisfied, the reconstruction support is reduced by setting x

Table 1 :

 1 Reconstruction accuracy using the BKDF based technique with 6 and 11 moments

		Mean absolute error
	Time
		BKDF: 6 moments BKDF: 11 moments
	t 0	0.0002	/
	t 1	0.0061	0.0014
	t 2	0.0026	0.0033
	t 3	0.0028	0.0024
	t 4	0.0068	0.0023

Table 2 :

 2 Comparison of the three different reconstruction methods (using the first six analytical moments) in terms of mean absolute error and computational cost

		BKDF based method Spline based method ME based method
	Time
		Error	Cost (s)	Error	Cost (s)	Error	Cost (s)
	t 0	0.0002	6.48	0.0011	6.48	0.0000	0.27
	t 1	0.0061	6.47	0.0105	5.86	0.0008	0.26
	t 2	0.0026	6.39	/	/	0.0005	0.23
	t 3	0.0028	6.31	/	/	0.0004	0.20
	t 4	0.0068	6.30	/	/	0.0004	0.20

Table 4 :

 4 Reconstruction accuracy using the ME based technique with six moments: comparison between the use of the moments computed with a three nodes quadrature and those obtained with a four nodes quadrature

		Mean absolute error
	Time
		QMOM-ME: M=3 QMOM-ME: M=4
	t 0	0.0000	0.0000
	t 1	0.0027	0.0010
	t 2	0.0037	0.0007
	t 3	0.0054	0.0019
	t 4	0.0080	0.0031

  By substituting equation B.6 in B.1, we obtain the final explicit solution given in equation 12.

	Nomenclature Abreviations
	a, b BKDF B n , Bn CFD	2σ 2 Lower and upper bounds of the reconstruction interval L 2 d L + µ e -1 2σ 2 Beta Kernel Density Function L 2 d L Coefficients of the BKDF system of equations Computational Fluid Dynamics	(B.3)
	Recall that : E(k) CLD	Global error function Chain Length Distribution
	e -ax 2 Mean absolute error for the reconstruction dx = Degree of Polymerization √ π 2 √ a erf (x Probability Density Function DQMOM Direct Quadrature Method of Moments E r (t) DP f (x)	√	a)	(B.4)
	and f M (x) KDF		The reconstructed Probability Density function Kernel Density Function
	g(x, κ, λ) LQDMA	Beta Probability Density Function Long Quotient-Modified Difference Algorithm
	h LSQLIN H[f ] ME	xe -ax 2 Constrained Linear Least Square Bandwidth of the BKDF Shannon entropy Maximum Entropy	dx = -	1 2a	e -ax 2	(B.5)
	Thus : I MWD		Number of BKDFs used for the reconstruction Molecular Weight Distribution
	I = σ 2 e -1 2 ( L-µ σ ) 2 moment's order Non-Negative Least Square K(x, x * , h) Beta kernel density function + µσ k, n NNLS PBE Population Balance Equation	π 2	1 -erf	1 √ 2	L -µ σ	(B.6)
	L, L PBM		Polymer chains length Polupation Balance Modelling
	L i PD		Abscissae of the Gaussian quadrature Product Difference Algorithm
	m PDF		Mean of the normal law Probability Density Function
	M QMOM	Number of nodes of the Gaussian quadrature Quadrature Method of Moments
	n(L, t)	Length based number density function
	N		Length of the set of moments
	p i		Weight coefficient, ∈ [0, 1]
	q		Length of the discretized interval of reconstruction
	t, t i , t f		Time, intial time, final time
	w i		Weights of the Gaussian quadrature
	x		Variable of the PDF
	x		Parameter of the shifted moments, ∈ [0, 1]
	x * , x i		Centers of the BKDFs
	Greek Symbols
	α	Factor of proportionality for the breakage frequency
	β	Breakage kernel
	γ	Breakage frequency
	Γ	Gamma function
	κ, λ Parameters of the Beta Probability Density Function
	µ	The set of moments
	σ	Standard deviation for the normal law
	ξ	Lagrange's multipliers

show that the ME based technique is highly dependent on the moments accuracy. When a three nodes quadrature is used for tracking the abscissas and the weights via the PD algorithm, the sets of associated moments lead to unsatisfactory reconstructions. Increasing the number of nodes improves significantly the moments accuracy and by the way the quality of the reconstruction. No doubt that the issue is common to all reconstruction methods: accurate values of the moments are required.

In summary, the ME based method is preferred than the BKDF and spline based methods because it requires a smaller number of moments and offers the best computational cost/reconstruction accuracy. The cases where a frequent reconstruction is needed along with the PBE resolution, the moments have to be computed accurately to acheive a good reconstruction quality using the ME based method. This is the main drawback when dealing with experimental moments. Further confrontation especially when different processes are coupled (e.g. breakage, aggregation, growth) is needed for assessing the relevance of this method.

Aknowledgments

The authors would like to thank Toulouse White Biotechnology and Région Midi-Pyrénées for their financial support. The Beta PDF is given, in its normalized form with κ and λ as parameters, by:

where Γ is the Gamma function. The Beta kernel is then defined by :

The Beta KDF parameters (x * and h) are related to the Beta PDF parameters (κ and λ) (see Appendix

A [START_REF] Athanassoulis | The truncated Hausdorff moment problem solved by using kernel density functions[END_REF]).

The moments of the Beta kernel can be analytically calculated thanks to the equation A.3 :

By substituting equation A.3 in 2, one obtains the moments µ * n of the reconstructed distribution :

For determining the unknown coefficients p i , the finite-moment problem is reformulated as a constrained optimization problem (equation A.5 ):

Equation A.5 can be solved numerically using a Nonnegative Least Square algorithm (NNLS). For improving the efficiency of this method and reducing the number of moments required for the reconstruction, different optimization techniques are available in literature such as the use, in addition of the given moments, of the shifted moments as described in the original work by [START_REF] Gavriliadis | Moment data can be analytically completed[END_REF].

The shifted moments of f (x) are defined by:

Since, in the finite-moment problem, f (x) is the target, the shifted moments μn are derived from the integer moments µ n as follows:

The shifted moments of the Beta KDF are calculated by:

Thus, the shifted moments μ * n of the reconstructed PDF are expressed as:

The minimization problem (equation A.5) is reformulated as :

This new constraint minimization problem is numerically solved using NNLS or LSQLIN (Matlab subroutines).

The algorithm of this method can be reduced to the following steps: In order to optimize the reconstruction support and thus increase the accuracy of the spline-based method, the values of the reconstructed PDF f (k) (x) are checked at the boundaries of the support [x

). For this, these two subintervals are divided into 10 equidistant small subintervalls, the values of f (k) (x) at the nodes of these small subintervals are checked if they are sufficiently small against the maximal value of f (k) in order to reduce eventually the size of the reconstruction support. This test is given for the right boundary as [START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF]:

n+1 ], the nodes are redistributed equidistantly.

It is further checked that any value of f (k) (x) at the nodes x i (i = 1, . . . , n + 1) and at the midpoints of the subintervals are greater than a given tolerance tol neg (equation A.12). This test is implemented to ensure the positivity of the PDFs.

with tol neg ≤ 0. The problem of the ill-conditioning of the linear system of equations is treated by regulizing the system. This is achieved by neglecting small singular values according to a given tolerance tol sing using a pseudo inverse routine of Matlab.

Note that in [START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF], an equidistant grid is used for the distribution of the nodes x i . This is not adapted for non-smooth and/or multimodal distributions. In the adaptive algorithm by [START_REF] De Souza | Reconstruction of a distribution from a finite number of moments with an adaptive spline-based algorithm[END_REF], this problem is solved by introducing a non-equidistant grid where the nodes are repositionned appropriately in order to capture accurately the PDF's critical domains.

This method is extensively described by [START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF]. The support of the target PDF

the target PDF is approximated by a piecewise polynomial s (l) (x) of degree l. The system of equations is detailed for cubic splines (l = 3). The unknowns are the four coefficients of the n splines. A smooth transition at the boundaries of the interval is assumed meaning that the PDF, its first and second derivatives are null at the boundaries, this gives 2 × 3 equations. The continuity of the splines, their first and second derivatives at the nodes provides 3(n -1) equations. The remaining (n -3) equations are supplemented by the moments. This leads to solve a 4n × 4n ill-conditioned linear system.

Since the set of the known moments limits the number of splines, one has to compute the reconstruction in the optimal support of f (x) thus, the reconstruction is computed iteratively starting from an initial reconstruction f (0) (x) in an initial interval [x

n+1 ].


S Hoang

A Georget

J Asselineau

Anne-Gaëlle Venier

Camille Leroyer

Anne-Marie M Rogues

Rodolphe Thiébaut
email: rodolphe.thiebaut@u-bordeaux.fr

A-G Venier

R Thie ´baut

Risk factors for colonization and infection by Pseudomonas aeruginosa in patients hospitalized in intensive care units in France

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

Introduction

Despite persistent and constant efforts made in Healthcare associated infections prevention, cross infections remain an issue especially in intensive care units (ICU). In 2012, 23% of patients acquired a healthcare associated infection during their hospitalization in an ICU that is a risk four times higher than in other hospital units [START_REF] Thiolet | Enquête nationale de pre ´valence des infections nosocomiales et des traitements anti-infectieux en e ´tablissement de sante ´, France[END_REF].

Cross infections are associated with high mortality rates and higher costs by complicating patients' cares and lengthening hospital stay [START_REF] Pittet | Nosocomial Bloodstream Infection in Critically III Patients: Excess Length of Stay, Extra Costs, and Attributable Mortality[END_REF][START_REF] Laupland | Cost of intensive care unit-acquired bloodstream infections[END_REF][START_REF] Sa ´nchez-Vela ´zquez | The Burden of Nosocomial Infection in the Intensive Care Unit: Effects on Organ Failure, Mortality and Costs. A Nested Case-Control Study[END_REF][START_REF] Tacconelli | Epidemiology, medical outcomes and costs of catheter-related bloodstream infections in intensive care units of four European countries: literatureand registry-based estimates[END_REF].

Pathogens responsible for those infections are often multi-drug resistant (MDR) bacteria. Among those pathogens, Pseudomonas aeruginosa is often encountered [START_REF] Ortona | A study on the incidence of nosocomial infections in a large university hospital[END_REF] and is responsible for severe infections, difficult to manage, such as ventilator associated pneumonia, bacteremia or skin infections, mainly in immunocompromised patients [START_REF] Manfredi | Pseudomonas spp. complications in patients with HIV disease: an eight-year clinical and microbiological survey[END_REF] with already a poor baseline prognosis. Indeed this bacteria, once it has colonized the patient digestive mucosa or skin, can lead to infections by immune deficiency and skin or mucosa breach by indwelling invasive device for example. When focusing on the bacteria strains it seems that P.aeruginosa related infections are more likely to occur to P.aeruginosa colonized patients, which highlights the link between the two conditions [START_REF] Go ´mez-Zorrilla | Prospective Observational Study of Prior Rectal Colonization Status as a Predictor for Subsequent Development of Pseudomonas aeruginosa Clinical Infections[END_REF][START_REF] Cohen | A prospective survey of Pseudomonas aeruginosa colonization and infection in the intensive care unit[END_REF].

P.aeruginosa transmission is known to be partly endogenous and exogenous but the precise roles of one and another remain uncertain. Human digestive flora has been described as the main source of endogenous transmission that can be enhanced by antibiotic pressure [START_REF] Go ´mez-Zorrilla | Prospective Observational Study of Prior Rectal Colonization Status as a Predictor for Subsequent Development of Pseudomonas aeruginosa Clinical Infections[END_REF]. Hydric environment such as faucets, invasive devices exposure and other colonized patients are responsible for exogenous transmission. Medical staff has been identified as a key carrier from patients to patients (cross-transmission) [START_REF] Obadia | Detailed Contact Data and the Dissemination of Staphylococcus aureus in Hospitals. Salathe ´M, e ´diteur[END_REF]. Although genotyping methods allow tracking the spread of bacteria, the respective contribution of the various factors remains unclear.

Hence, the factors associated with P.aeruginosa colonization or infection are variable from one study to another [START_REF] Falagas | Risk factors for the isolation of multi-drug-resistant Acinetobacter baumannii and Pseudomonas aeruginosa: a systematic review of the literature[END_REF] including: length of hospitalization [START_REF] Bergmans | Cross-colonisation with Pseudomonas aeruginosa of patients in an intensive care unit[END_REF] severity of underlying disease [START_REF] Bergmans | Cross-colonisation with Pseudomonas aeruginosa of patients in an intensive care unit[END_REF][START_REF] Talon | Risks and Routes for Ventilator-Associated Pneumonia with Pseudomonas aeruginosa[END_REF][START_REF] Cobos-Trigueros | Acquisition of Pseudomonas aeruginosa and its resistance phenotypes in critically ill medical patients: role of colonization pressure and antibiotic exposure[END_REF] exposure to invasive procedures [START_REF] Talon | Risks and Routes for Ventilator-Associated Pneumonia with Pseudomonas aeruginosa[END_REF][START_REF] Cobos-Trigueros | Acquisition of Pseudomonas aeruginosa and its resistance phenotypes in critically ill medical patients: role of colonization pressure and antibiotic exposure[END_REF][START_REF] Halwani | Cross-transmission of nosocomial pathogens in an adult intensive care unit: incidence and risk factors[END_REF], contamination of water tap [START_REF] Bergmans | Cross-colonisation with Pseudomonas aeruginosa of patients in an intensive care unit[END_REF][START_REF] Thuong | Epidemiology of Pseudomonas aeruginosa and risk factors for carriage acquisition in an intensive care unit[END_REF][START_REF] Berthelot | Prospective study of nosocomial colonization and infection due to Pseudomonas aeruginosa in mechanically ventilated patients[END_REF][START_REF] Agodi | Pseudomonas aeruginosa carriage, colonization, and infection in ICU patients[END_REF], close contact with contaminated patients [START_REF] Cobos-Trigueros | Acquisition of Pseudomonas aeruginosa and its resistance phenotypes in critically ill medical patients: role of colonization pressure and antibiotic exposure[END_REF], antibiotic selective pressure [START_REF] Thuong | Epidemiology of Pseudomonas aeruginosa and risk factors for carriage acquisition in an intensive care unit[END_REF][START_REF] Berthelot | Prospective study of nosocomial colonization and infection due to Pseudomonas aeruginosa in mechanically ventilated patients[END_REF]. The heterogeneity of the factors reported is due partly to the complexity of the measurement of exposures, the definition of the outcome and of the methods used for the analysis of the association with P.aeruginosa colonization or infection. In a previous report, we identified individual risk factors and environmental factors associated to P.aeruginosa acquisition in a cohort of 1314 patients hospitalized in ICU in France [START_REF] Venier | Risk factors for Pseudomonas aeruginosa acquisition in intensive care units: a prospective multicentre study[END_REF]. However, we did not distinguish between colonization and infection, both were merge in the same outcome although factors could be different in the two situations with different clinical consequences. The aim of our study was to assess the distinct role of environment, antibiotic and patient condition toward both colonization and infection, and to identify an interaction between environment and antibiotic selective pressure.

Materials and methods

Study setting

In 2009, an observational prospective multicentric study, entitled DYNAPYO, was performed, involving ten French ICUs for a five months study period. There were two surgical, four medical and four mixed ICUs (both medical and surgical units). Six university hospitals (Besanc ¸on, Bordeaux, Garches, Lyon, Montpellier and Paris) and two general hospitals (Lens and Tourcoing) took part in the study and included each one ICU, excepted for Besanc ¸on and Lyon which included two ICUs. These ICUs amounted between 9 and 20 beds and between 10 and 47 water taps. Data collection regarding environment and patient population was performed by trained healthcare professionals and was registered into a secure online electronic form.

Sample size

In 2006 a pilot study was performed. In that study, the patients exposed to hydric contamination had a P.aeruginosa acquisition risk of 1.7. Within 6 months, the incidence of P.aeruginosa acquisition was 5.7%. To confirm such size effect with a statistical power of 80% and type I error of 5%, at least 112 events should be recorded. According to the expected incidence, it required to include around 2000 patients. However, the present study was stopped prematurely because the incidence of acquisition was higher than expected. Every newly hospitalized adult patient in ICUs with no P.aeruginosa carriage within the 48 first hours after admission was included. No P.aeruginosa carriage means that samples (rectal oropharyngeal or sputum swabs) were found negative at admission and within the 48 first hours. During the follow-up period, they were weekly screened (rectal oropharyngeal or sputum swabs) until they died or were discharged. Data were collected through medical record review: those about clinical and medical conditions and prior antibiotic use were collected at admission, and follow-up data such as hydric environment pressure, indwelling invasive device exposure and antibiotic support were regularly recorded.

To estimate the hydric environment pressure, faucets were weekly sampled in the morning before use and disinfection. Specific sample protocol has been described in the previous report [START_REF] Venier | Risk factors for Pseudomonas aeruginosa acquisition in intensive care units: a prospective multicentre study[END_REF]. No hygiene protocol modification could happen within the study period because of these results as the units were blinded.

Samples were then analyzed by the bacteriology laboratories where aerobic cultures were performed.

Data were made anonymous. All patients were informed of the survey and their rights. Signed consents were not required. Local ethical committee [Comite ´Consultatif sur le Traitement de l'Information en matière de Recherche dans le domaine de la Sante ´(CCTIRS) et Commission de l'Informatique et des Liberte ´s (CNIL)] approved the study.

Determination of colonization and infection

The first patient colonization and the first patient infection by P.aeruginosa were the two outcomes of interest. Patient colonization was defined as the presence of P.aeruginosa among at least one of the screening sample site from the 48 th hour until discharge. If any colonization was identified within the 48 th hour, patients were considered as imported cases and were not included in the analysis but still considered in exposure factor calculation. Colonization status was unknown for healthcare workers so that no change in usual isolation procedure occurred during the study.

Patient infection was defined according to current recommendations from REA-RAISIN network [START_REF]Surveillance des infections nosocomiales en re ´animation adulte-Protocole 2014[END_REF], as clinical and biological infection features associated with the isolation of the bacteria in any samples (blood culture, urine culture. . .).

Main exposure factors

We focus on two exposure factors of main interest. First, hydric environment pressure was measured by sampling performed at the faucets located in the patient room: Tap water contamination at the entry in the patient room was the first main exposure factor. Antibiotic treatments prescribed since the day before, were the second main exposure factor and were daily recorded. We distinguished exclusive inactive antibiotics from active antibiotics against P.aeruginosa which includes at least one of the following: ureido and carboxypenicillins, antipseudomonal cephalosporins, carbapenems, colimycin fosfomycin, fluoroquinolones, and aminoglycosides. The antibiotic activity was determined by the theoretical sensivity of antibiotics on a wild-type P.aeruginosa.

Other risk factors

Other individual risk factors were: age, sex, immucompromized status, IGSII score (i.e patient disease severity at the admission in ICUs predicting mortality), Charlson comorbidity index for chronic disease, cause of admission, prior antibiotics prescription before admission, prior hospitalization within the year prior to admission, history of prior carriage of P.aeruginosa (i.e prior colonization or diagnostic of P.aeruginosa infection). Environment risk factors were: any specialization of the ICU (medical, surgical or mixed), center and prior occupation by a colonized patient before entry in the room. Risk factors associated with medical care were: invasive mechanical ventilation exposure and broncho-endoscopy exposure.

Statistical analysis

Analysis was performed with a multi-state model that included five states and eight transitions (Fig 1). Baseline state 1-Hospitalization was the baseline status of patients hospitalized without colonization nor infection. States 2-Colonization and 3-Infection reflected the occurrence of the first episode of colonization and infection. We assumed that the state 3-Infection was systematically preceded by the state 2-Colonization. Arrows represent transitions from a state to another. The model assumed that transitions were markovian, which means that history of past states (i.e. duration in the state) was not taken into account in the estimation [START_REF] Putter | Tutorial in biostatistics: competing risks and multi-state models[END_REF]. A discrete non-homogeneous model was assumed: discrete because time unit is the day; nonhomogeneous because time-dependent variables were used as adjustment variables [START_REF] De Uña-A ´lvarez | Multi-state models for the analysis of time-to-event data[END_REF]. We estimated the associations between risk factors and events of interest by modelling transition intensities. The transition intensities represent the instantaneous risk of going from one state to another (e.g. acquiring infection from carrying P.aeruginosa).

Analyses were performed with the package "MSM" developed by C. Jackson in the R 1 software [START_REF] Jackson | Multi-state models for panel data: the msm package for R[END_REF]. Weekly screenings for both patients and water taps led to interval censoring which was taken into account by the model. Hydric environment exposure was assumed to be constant between two screenings. In case of punctual missing values regarding follow-up data, the package allows missing at random (MAR) imputation; in case of baseline or all along missing data, patients were removed from analyses. To minimize loss of patients, conditional imputation procedure was implemented to replace missing data for IGSII score, history of prior carriage of P.aeruginosa, prior hospitalization within the year prior to admission and prior antibiotics prescription before admission. To answer the aims of the study, we ran the model transition by transition with the following order of interest: "Hospitalization-Colonization", "Colonization-Infection", "Colonization-Death", "Infection-Death" and "Hospitalization-Death". We did not focus on the transition to discharge, as the event was only considered as a competing risk.

We first performed an univariable analysis including ten baseline factors and six binary time-dependent factors in each transition. Our three main exposure factors were forced in the model whatever their association, apart from hydric environment pressure which was specifically studied in the transition to colonization for the sake of clinical plausibility. We performed a multivariable analysis including significant variables with a p-value<25% using a step-bystep forward approach for each transition (p<5%). Adjustment for site was performed and we evaluated the interactions between hydric environment and both active and inactive antibiotics selective pressure in the final model in the transition leading to colonization. Hazard ratios and their 95% confidence intervals were provided as well as likelihood ratio test p-values. Loglinearity of the risk for quantitative variables has been checked.

Results

Study population

In the DYNAPYO cohort, 1700 patients were newly hospitalized, including 1314 P.aeruginosafree patients who were enrolled in the study. Patients characteristics are shown in Tables 1 and2.

The median age was 59 years and 11% of patients were immunocompromised. The median IGSII score was 42 points (Interquartile Range IQR: 30; 56). The median (IQR) length of stay was 6 days (3; 13) and 232 patients (17.7%) died in a median delay of 7.5 days (3; 15). First colonization and first infection were respectively estimated to occur at 10 days (7; 16) and 12 days (8; 19). During the follow-up, 70% were mechanically ventilated and among them, almost 90% were ventilated at the admission.

A total of 855 positive hydric environment samples among 4999 screened were positive for P. aeruginosa (17.1%).

Colonization and infection incidences

Among the 1314 patients, 201 patients were colonized during a total of 11 915 days follow-up, corresponding to an incidence of 16.9 colonizations per 1000 patient-days ( 95% CI [14.6-19.2]). An infection was observed in 87 patients during a 12 608 days of follow-up, representing 43.3% of P.aeruginosa carriers and an incidence of 6.9 infections per 1000 patient-days ( 95% CI [5.5-8.3]).

Risk factors analysis

1305 patients were included in the analysis in the state 1-Hospitalization whereas 9 patients were excluded from this analysis because of missing data for one or more variables. Among the 1305 patients, 180+200+914 moved to another state, meaning that 11 stayed in the state 1-Hospitalization at the end of the study (right-censoring). Also two and four patients stayed Table 3 shows the main results after adjustment for the site. The risk of colonization was significantly increased when inactive antibiotic were prescribed (HR = 1.60; 95% CI = [1.15-2.21]), whereas it was decreased when active antibiotic were prescribed (HR = 0.67; 95% CI = [0.48-0.93]). Tap water contamination at the entry in the patient room enhanced the risk of colonization by +66% (HR = 1.66; 95% CI = [1.01-2.75]). History of P. aeruginosa carriage and invasive ventilation were also significant independent risk factors of P. aeruginosa colonization: HR = 4.03; 95% CI = [1.96-8.27] and HR = 4.70; 95% CI = [2.66-8.31], respectively. Tests for interactions between the hydric environment and inactive or active antibiotics were not statistically significant (p-values 0.28 and 0.78, respectively).

The factors associated with P.aeruginosa infection were limited. Among colonized patients, the prescription of active antibiotic tended to be associated with a lower risk of infection occurrence (HR = 0.64; 95% CI = [0.41-1.01]). There was no significant association with inactive antibiotic support (HR = 0.77; 95% CI = [0.47-1.26]).

Older age (HR = 1.12; 95% CI = [1.01-1.23]) increased the risk of death during hospitalization. Whatever the colonized status, a higher IGSII score (HR = 1.46; 95% CI = [1.34-1.60] and HR = 1.33; 95% CI = [1.06-1.66]) and the diagnosis of coma (HR = 1.79; 95% CI = [1.13-2.84] and HR = 3.69; 95% CI = [1.20-11.33]) were associated with an increased risk of death.

Discussion

In an endemic situation, this study showed the constant presence and role played by P aeruginosa with incidence density rates of 16.9 colonizations and 6.9 infections per 1000 patient-days. Hydric environment was identified as a consistent and independent risk factor, but also individual and health risk factors such as prior carriage of the bacteria, mechanical ventilation and antibiotic use were found to contribute to the transmission of the bacteria. No interaction between antibiotic use and hydric environment contamination was detected. Therefore, an inappropriate prescription of antibiotics does not seem to increase the risk of contamination by the hydric environment.

These results highlighted the consistent role of the hydric environment in P.aeruginosa transmission and the independent effect of both active and inactive antibiotics support in colonization occurrence. The prevalence of colonized tap water was of 17.1%, quite similarly in others studies [START_REF] Zhou | Sources of sporadic Pseudomonas aeruginosa colonizations/infections in surgical ICUs: Association with contaminated sink trap[END_REF]. Our results on risk factors are consistent with those from previously published studies [START_REF] Talon | Risks and Routes for Ventilator-Associated Pneumonia with Pseudomonas aeruginosa[END_REF][START_REF] Thuong | Epidemiology of Pseudomonas aeruginosa and risk factors for carriage acquisition in an intensive care unit[END_REF][START_REF] Agodi | Pseudomonas aeruginosa carriage, colonization, and infection in ICU patients[END_REF][START_REF] Bertrand | Endemicity, molecular diversity and colonisation routes of Pseudomonas aeruginosa in intensive care units[END_REF][START_REF] Rogues | Contribution of tap water to patient colonisation with Pseudomonas aeruginosa in a medical intensive care unit[END_REF]. Environment has been identified as a major risk factor in the exogenous transmission. Petignat et al. have shown the impact of infection control measures targeting faucets toward P.aeruginosa colonization [START_REF] Petignat | Exogenous Sources of Pseudomonas aeruginosa in Intensive Care Unit Patients: Implementation of Infection Control Measures and Follow-Up With Molecular Typing[END_REF]. Prior carriage of P.aeruginosa and mechanical ventilation were also found as important risk factors [START_REF] Falagas | Risk factors for the isolation of multi-drug-resistant Acinetobacter baumannii and Pseudomonas aeruginosa: a systematic review of the literature[END_REF]. Active antibiotics are described as protective factors [START_REF] Cobos-Trigueros | Acquisition of Pseudomonas aeruginosa and its resistance phenotypes in critically ill medical patients: role of colonization pressure and antibiotic exposure[END_REF][START_REF] Berthelot | Prospective study of nosocomial colonization and infection due to Pseudomonas aeruginosa in mechanically ventilated patients[END_REF][START_REF] Martı ´nez | Influence of antipseudomonal agents on Pseudomonas aeruginosa colonization and acquisition of resistance in critically ill medical patients[END_REF], whereas inactive antibiotics have already been identified as risk factors [START_REF] Go ´mez-Zorrilla | Prospective Observational Study of Prior Rectal Colonization Status as a Predictor for Subsequent Development of Pseudomonas aeruginosa Clinical Infections[END_REF][START_REF] Thuong | Epidemiology of Pseudomonas aeruginosa and risk factors for carriage acquisition in an intensive care unit[END_REF]. In the present study, the antibiotics activity was defined according to the theoretical knowledge of the effect of each molecule and not by the antibiogram for any participant. This could represent an important limit in our study. However, the observed opposite effects of active and inactive antibiotics is expected: inactive antibiotics, involving an unbalance of the digestive microflora can indeed allow the bacteria to proliferate, in contrary to active antibiotics which remove the bacteria.

Although we did not find any evidence for an interaction between the environment and antibiotics, it does not mean that there was none. Either there is no interactive role between antibiotics and environment, or there is one which could have been missed by a lack of statistical power or by an ignoring confounding factor.

Routine screening for P.aeruginosa carriage is not recommended [START_REF] Bernard | Pre ´vention de la transmission croise ´e: pre ´cautions comple ´mentaires contact[END_REF], but this could be questionned: first we identified prior P.aeruginosa carriage as a risk factor which could be linked to clinical predisposition to acquire the bacteria once again. Second, we found that 43% of colonized patients developed infectious conditions, and these results are consistent with Gomez-Zorilla et al. [START_REF] Go ´mez-Zorrilla | Prospective Observational Study of Prior Rectal Colonization Status as a Predictor for Subsequent Development of Pseudomonas aeruginosa Clinical Infections[END_REF] who ascertained that infections occurred much frequently among colonized than non-colonized patients (39% vs 3.4% p<0.001). Third, we highlighted the protective role of active antibiotics and featured that deaths occurred more likely with infected patients (22/86) than colonized patients (30/200). Thus, according to medical history and condition, for example for immunocompromised patients or among those that received common antibiotics (which are mainly inefficient against the bacteria), this screening associated with genotyping resistance test, could avoid delayed efficient antibiotics and so, avoid death [START_REF] Micek | Pseudomonas aeruginosa Bloodstream Infection: Importance of Appropriate Initial Antimicrobial Treatment[END_REF].

Thanks to the use of multistate model that allows to distinguish colonization from infection, we could demonstrate that indeed very few factors were associated with the risk of infection: the exposure to active antibiotic especially. Also, usual risk factors of death were isolated (age, IGSII, coma).

Furthermore as we carefully checked the link between anterior contaminated samples in faucets and patient incident colonization, these results strengthen what we found in the risk factors analysis.

Although a seven days interval screening seems to be a good compromise between feasibility and costs [START_REF] Zhou | Sources of sporadic Pseudomonas aeruginosa colonizations/infections in surgical ICUs: Association with contaminated sink trap[END_REF][START_REF] Rogues | Contribution of tap water to patient colonisation with Pseudomonas aeruginosa in a medical intensive care unit[END_REF], test performances with a moderate sensitivity [START_REF] Slekovec | Faut-il proce ´der à un de ´pistage du portage de Pseudomonas aeruginosa dans les services de re ´animation?[END_REF] and number of screening samples could have led to classification error. Colonization incidence could have been overestimated because of wrong inclusion of unknown P.aeruginosa carriers who had been then identified as incident cases. Hence this could explain the association between prior carriage of P.aeruginosa and colonization. Also, water tap sampling was also weekly collected and previous notes can also be applied: water tap status was considered constant until the next screening. This situation has probably led to a decrease of the association between water tap pressure and colonization or infection. However, a more intense screening is difficult to organize. Future studies may take advantage of new technologies such as the i-BIRD project [START_REF] Obadia | Detailed Contact Data and the Dissemination of Staphylococcus aureus in Hospitals. Salathe ´M, e ´diteur[END_REF] in which spatio-temporal data of contact between healthcare workers are measured by sensitive sensors.

In conclusion, hydric contamination and antibiotics pressure seem to remain key independent risk factors in P.aeruginosa colonization. These results advocate the need to carry on preventive and targeted interventions toward healthcare associated infections. This study encourages the settlement of cautious healthcare associated infections prevention measures such as hands hygiene, respect of aseptic rules and more regular faucets maintenance.

Fig 1 .

 1 Fig 1. Multistate model representation. https://doi.org/10.1371/journal.pone.0193300.g001

  /doi.org/10.1371/journal.pone.0193300.t001 respectively in the states 2-Colonization and 3-Infection at the end of the study and were right-censored. Fig 2 shows the evolution of those 1305 patients during hospitalization.

Fig 2 .

 2 Fig 2. Multistate model and number of subjects switching from one state to another. By definition, the 1305 patients started their hospitalization in state 1-hospitalization. n ij : number of patients moving from state i to state j. https://doi.org/10.1371/journal.pone.0193300.g002

Table 1 . Patients characteristics at admission.

 1 Dynapyo cohort, France 2009.

	Patient characteristics	Included patients (n = 1314)
		Total	Proportion in %
	Sex
	Male	772	58.8
	Female	542	41.2
	Immunodeficiency
	Yes	144	11.0
	No	1160	88.2
	Unknown	10	0.8
	Prior hospitalization in the year before admission
	Yes	598	45.5
	No	699	53.2
	Unknown	17	1.3
	Prior carriage of P. aeruginosa
	Yes	27	2.1
	No	1220	92.8
	Unknown	67	5.1
	Surgery within 30 days before admission
	Yes	283	21.5
	No	1026	78.1
	Unknown	5	0.4
	Prior antibiotic administration before admission
	Active antibiotics	51	3.9
	Inactive antibiotics	187	14.2
	Active and inactive antibiotics	136	10.4
	No antibiotics	866	65.9
	Unknown	74	5.6
	Cause of hospitalization
	Coma	388	29.5
	Acute respiratory insufficiency	367	27.9
	Septic shock	176	13.4
	Others	383	29.2
	Intensive care unit (ICU)
	Medical ICU	604	46.0
	Surgical ICU	179	13.6
	Mixed ICU	531	40.4
	Patient characteristics	Included patients (n = 1314)
		Total	Median (IQR)
	Age (in years)	1314	59 (43-73)
	IGS-II score at admission	1284	42 (30-56)

Table 2 . Exposure characteristics during the follow-up.

 2 Dynapyo cohort, France 2009.

	Exposure

characteristics Included patients (n = 1314) Total Proportion in % Mechanical ventilation exposure >24h

	Yes	903	68.7
	No	407	31.0
	Unknown	4	0.3
	Bronchoscopic endoscopy
	Yes	127	9.7
	No	1183	90.0
	Unknown	4	0.3
	Exposure to a room previously occupied by a colonized patient
	Yes	320	24.4
	No	990	75.3
	Unknown	4	0.3
	Tap water of the room contamined by P.aeruginosa
	Yes	315	24.0
	No	995	75.7
	Unknown	4	0.3
	Exposure to active antibiotics
	Yes	491	37.4
	No	818	62.2
	Unknown	5	0.4
	Exposure to inactive antibiotics
	Yes	601	45.7
	No	708	53.9
	Unknown	5	0.4
	https://doi.org/10.1371/journal.pone.0193300.t002

Table 3 . Multistate model after adjustment for the site: Risk factors for P. aeruginosa colonization, infection and risk factors for death

 3

	Variables	Hospitalization-Colonization	Colonization-Infection
		(n 1 = 1305) (n 2 = 200)		(n 2 = 200) (n 3 = 86)
		HR	95% CI	p-value	HR	95% CI	p-value
	Contaminated hydric environment exposure at the entry room	1.66	1.01-2.75	<0.05
	Active antibiotics exposure	0.67	0.48-0.93	0.02	0.64	0.41-1.01	<0.05
	Inactive antibiotics exposure	1.60	1.15-2.21	<0.01	0.77	0.47-1.26	0.30
	Prior P.aeruginosa carriage	4.03	1.96-8.27	<0.01
	Invasive mechanical ventilation	4.70	2.66-8.31	<0.0001
	Variables	Hospitalization-

. DYNAPYO Cohort-2009.

Death Colonization-Death Infection-Death (n 1 = 1305) (n 4 = 180) (n 2 = 200) (n 4 = 30) (n 3 = 86) (n 4 = 22)

		HR	95% CI	p-value	HR	95% CI	p-value	HR	95% CI	p-value
	Active antibiotics exposure	0.65	0.47-0.90	<0.01	0.41	0.19-0.87	0.02	1.67	0.35-8.07	0.51
	Inactive antibiotics exposure	0.83	0.60-1.14	0.24	0.86	0.36-2.04	0.77	0.66	0.25-1.71	0.41
	Age (% 10 years)	1.12	1.01-1.23	0.03
	IGS II score up to 24 hours (% 10 points)	1.46	1.34-1.60	<0.0001	1.33	1.06-1.66	0.01
	Admission diagnosis (ref.Others)			0.06			0.02
	Coma	1.79	1.13-2.84		3.69	1.20-11.33
	Acute respiratory distress	1.21	0.75-1.96		2.06	0.66-6.41
	Septic shock	1.16	0.70-1.92		0.51	0.10-2.68
	HR = Hazard ratio
	95% CI = 95% Confidence interval
	p-value: obtained with Likehood ratio test
	Interactions:
	Contaminated hydric environment à Active antibiotics exposure p = 0.78
	Contaminated hydric environment à Inactive antibiotics exposure p = 0.28
	https://doi.org/10.1371/journal.pone.0193300.t003

PLOS ONE | https://doi.org/10.1371/journal.pone.0193300March 9, 2018

Acknowledgments

The authors declare that they have no competing interests.

be written between Bordeaux University Hospital and the applicant for any request for data transfer outside the organization. The authors received no specific funding for this work.

The data underlying this study has been restricted by Commission Nationale Informatique et Liberte ´(CNIL: https:// www.cnil.fr/en/home). The data contains sensitive information (medical data) that could be, even anonymized, connected to patients through sufficient identification data in the data set (such as age,

Author Contributions

Conceptualization: A-G. Venier, A. M. Rogues, R. Thie ´baut.

Formal analysis: S. Hoang, A. Georget, J. Asselineau.