Otman Maalal

Marc Prat
email: mprat@imft.fr

René Peinador

Didier Lasseux

Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms

The pore size distribution of a porous medium is often estimated from the retention curve or the invading fluid flow rate curve using simple relationships more or less explicitly based on the consideration that the porous medium is made of a bundle of cylindrical parallel tubes. This type of determination is tested using pore network simulations. Starting from two-or three-dimensional networks, the characteristics of which are known apriori, the estimation of the throat size distribution (TSD) is performed using the standard methods in the case of drainage. Results show a significant discrepancy with the input data. The disagreement is more pronounced when the fluid flow rate curve is employed together with the parallel tubes assumption. The physical origins of these shortcomings are identified. A method, based on pore network simulations combined with a genetic algorithm and the hill climbing algorithm, is then designed, which makes simultaneous use of the nonwetting fluid flow rate curve and the retention curve of the medium. Very significant improvement is achieved in the estimation of the TSD using this procedure.

I. INTRODUCTION

The pore size distribution (PSD) is a key structure characteristic of a porous medium of major importance to model flow and transport processes. A simple method for determining the PSD is to use the retention curve, i.e., the fluid potential to saturation relationship [1]. As discussed for instance in [2,3], this is often performed assuming more or less explicitly that the pore system can be represented by a bundle of parallel cylindrical capillary tubes of circular cross section with a random distribution of radii. However, except for a few particular systems, for example the Anopore membrane discussed in [4], the microstructure of a porous medium is generally significantly different from a simple bundle of parallel cylindrical tubes. Most porous materials are rather characterized by a system of highly interconnected pores. As a result, the information on the PSD obtained from simple interpretation based on this representation is questionable and must be cautiously assessed [3]. Another issue is that the PSD obtained in this way makes no difference between the throats, i.e., the constrictions in the pore space, and the pore bodies, defined here as the local larger cavities. Making such a distinction between pore throats and pore bodies is a classical feature of pore network models (PNMs) [5]. This is justified by the fact that constrictions are largely controlling the capillary effects since menisci tend to get preferentially pined in constrictions. They are also decisive for transport phenomena since, for example, the pressure drop between two adjacent pore bodies can be essentially controlled by the constriction connecting them when this constriction is sufficiently narrow compared to the pore body size [6,7]. This distinguishes the PSD, which corresponds to the pore bodies in the PNM partition of the pore space, from the throat size distribution (TSD), i.e., the distribution of constriction sizes. In this context, the question arises as to whether it is possible to identify both the PSD and TSD, or at least the TSD, from relatively basic information such as the retention curve or the flow rate curve, the latter denoting the flow rate to fluid potential relationship. With the development of imaging techniques, e.g., micro-x-ray tomography, this approach may be itself questioned since the TSD and the PSD can be determined by means of various algorithms applied to segmented images of the microstructure [8,9]. However, imaging devices are expensive and the segmentation of the pore space in pore bodies and pore throats is rarely straightforward, requiring the treatment of a huge amount of data. Moreover, obtaining reliable images for materials with pore bodies and throats in the submicronic range is not an easy task. As a result, cheaper and more routine characterization techniques are still quite desirable.

In what follows, the focus is laid on the TSD. The intrinsic permeability and porosity, together with the retention and nonwetting fluid flow rate curves, are supposed to be known, i.e., are available from measurement for a pair of fluids for which the wetting fluid is perfectly wetting. In addition to the assessment of two classical TSD extraction methods, the objective of the paper is to explore whether this information is sufficient to characterize the TSD. Taking into account the nonwetting fluid flow rate is motivated by the fact that the evaluation of a porosimetry technique, referred to as the fluid-fluid displacement porosimetry (FFDP) [10,11], is also of interest. In this technique, a perfectly wetting liquid is displaced by a nonwetting one through a series of pressure steps and the nonwetting flow rate through the sample is recorded for each pressure step. The flow rate vs pressure curve is then harnessed to extract the TSD [10]. In the context of FFDP, the situation under study in the present paper can actually be seen as a quasistatic version, referred to as QSFFDP, where QS stands for quasistatic. Contrary to the classical FFDP technique where the viscous pressure drop across the sample in the nonwetting fluid is on the order of the capillary pressure, the viscous pressure drop remains small compared to the capillary pressure. Actually, this corresponds to the situation considered in practice for determining the nonwetting fluid relative permeability. Considering the flow rate curve so as to use the same convention as in the FFDP or the nonwetting fluid relative permeability curve is strictly equivalent in the quasistatic limit under consideration here. The classical FFDP (as opposed to the QSFFDP) will be the object of a future study.

On this basis, the idea is then to rely on pore network simulations. The approach is similar to the one presented in [3]. Starting from a perfectly known microstructure with interconnected pore bodies, referred to as the reference network, an immiscible displacement is performed numerically for conditions similar to the ones employed for the classical retention curve measurement. The TSD is extracted from the retention curve and flow rate data using classical extraction methods based on the bundle of capillary tubes model. The TSD obtained with this procedure is compared to the actual TSD. The decisive feature here is that the TSD is perfectly known a priori for the reference medium. This allows illustrating and assessing the limitations of the classical methods for the TSD determination. Then, the question arises as to how the use of the immiscible displacement data can be improved when the microstructure is clearly not a simple system of parallel tubes. This amounts to considering an ill-posed inverse problem in which the objective is to determine the TSD. Regarding the method used to compute the displacement in the network, main possibilities can be divided into direct simulations like lattice Boltzmann methods [12] or volume of fluid models [13], etc., and PNMs [5]. As will be explained below, the TSD identification method requires computing the immiscible displacement many times. In this respect, PNMs are well adapted as this technique requires much shorter computational time compared to direct simulations.

The present paper is organized as follows. In Sec. II, two classical methods for determining the PSD (more exactly the TSD) from the nonwetting fluid flow rate curve or the retention curve are recalled. The pore network models are presented in Sec. III. The algorithms used to simulate the displacement and compute the flow rate and retention curves are presented in Sec. IV. The assessment of the TSD determination from the classical methods is presented in Sec. V. In Sec. VI, a method to determine the TSD from both the nonwetting fluid normalized flow rate and retention curves using PNM simulations combined with an optimization method is presented. A short discussion is presented in Sec. VII. Conclusions are drawn in Sec. VIII.

II. FLOW RATE POROSIMETRY AND RETENTION CURVE POROSIMETRY

Suppose the retention curve, P c (S), flow rate curve, Q NW (P c ), and porosity, ε, are known, P c (S) being the capillary pressure, Q NW the nonwetting fluid flow rate for the quasistatic fluid distribution corresponding to P c (S), and S the wetting fluid saturation. Considering the generalized Darcy's law implies

Q NW = A kk r (P c ) μ P L ( 1
)
where μ is the nonwetting fluid viscosity, A is the porous medium cross-sectional surface area, L is the sample length, k is the porous medium intrinsic permeability, and k r is the nonwetting fluid relative permeability. In addition, P is the pressure difference between the inlet and outlet in the nonwetting fluid. In the quasistatic limit under consideration here, P is small compared to the capillary pressure. Throughout the paper, a constant P is considered independently of the capillary pressure.

To extract the TSD from the retention curve, a pore size, r, is associated to a value of P c using the Young-Laplace equation

P c = 2γ cosθ r (2)
where the contact angle is θ = 0 in our case (perfectly wetting displaced fluid) and γ is the interfacial tension. For simplicity, trapping phenomena [14] are assumed negligible so that S varies in the interval [0,1]. The wetting fluid saturation, S, corresponding to r is then simply defined as [15]

S = ∫ r 0 V p (r) f (r)dr ∫ +∞ 0 V p (r) f (r)dr (3)
where V p (r) is the volume of a pore of size r and f (r) is the pore size distribution, i.e., the pore size probability density function (p.d.f.). Deriving Eq. ( 3) with respect to r yields

dS dr = V p (r) f (r) ∫ +∞ 0 V p (r) f (r)dr . ( 4
)
Since, by definition,

∫ +∞ 0 V p (r) f (r)dr = ε V PM
, where V PM is the porous medium volume, this leads to

f (r) = εV PM V p (r) dS dr . ( 5
)
In discrete form, Eq. ( 5) can be expressed as

f k (r AB ) = εV PM V p (r AB ) S(P cB ) -S(P cA ) r B -r A ( 6
)
where P cA and P cB are two close values of P c and r AB = (r A + r B )/2, r A and r B being determined from Eq. ( 2). Within the framework of the model of a bundle of parallel cylindrical tubes, the volume of a pore is expressed as

V p (r) = π r 2 L
where L is the tube length assuming straight tubes. Under these circumstances, Eq. ( 6) is expressed as

f k (r AB ) = εV PM π L r -2 AB S(P cB ) -S(P cA ) r B -r A . ( 7
)
The value of f k is then normalized according to

d k (r AB ) = f k n i=1 f i . (8)
Here, n is the number of classes considered to compute the discrete p.d.f. Hence, d k (r AB ) represents the fraction of pores (number of pores divided by total number of pores) the sizes of which are in class k.

The method for determining the TSD from the flow rate Q NW (P c ) is similar and also explicitly relies on the model of a bundle of parallel cylindrical tubes. It is also similar to the one used in the FFDP technique [10]. In the bundle of tubes model, the flow rate is computed applying Poiseuille's law in each tube containing the nonwetting fluid:

Q NW = Pπ 8μL +∞ r min r 4 f Q (r)dr (9)
where the tube length is again L assuming straight tubes and

r min = 2γ cosθ P c , ( 10
)
f Q being again the TSD. In principle, f Q = f . However, a different notation is used since it is not obvious that both procedures lead to the same TSD when the porous medium microstructure is different from a bundle of parallel tubes.

Combining Eqs. ( 1) and ( 10), assuming P = const P c independently of the capillary pressure, and deriving the ensuing equation with respect to r lead to

dQ NW dr = - Pπ 8 μL r 4 f Q (r) (11)
and thus

f Q (r) = - 8 μL Pπ r -4 dQ NW dr . ( 12
)
Equation ( 12) can be expressed in discrete form as [10] f Qk (r AB ) = -

8 μL Pπ r -4 AB Q NWB -Q NWA r B -r A . ( 13
)
The value of f Qk is then normalized as

d Qk (r AB ) = f Qk n i=1 f Qi ( 14
)
In what follows, water is considered as the wetting fluid and air as the displacing fluid. Thus γ = 0.072 N/m and μ = μ nw = 10 -5 Pa s.

III. INTERCONNECTED MODEL POROUS MEDIUM

The QS immiscible displacement is simulated considering simple model porous media referred to as PNMs. In a PNM, the pore space is represented as a network of pore bodies connected by narrower channels. The pore bodies correspond to locally larger volumes whereas the channels correspond to the throats of the pore space, i.e., the narrower passages, also referred to as the constrictions, between pore bodies. Although it is possible to construct the network from digital images of real microstructures (see, e.g., [16]), simpler networks are sufficient for testing the QSFFDP. In what follows, square [two-dimensional (2D)] and cubic [three-dimensional (3D)] networks are considered. As illustrated in Fig. 1, pore bodies are spheres of radius r p located at the nodes of a regular grid. The distance between two adjacent nodes is the lattice spacing, a, with a = 350 μm in the simulations presented below. The size of the network is the number of nodes (pore bodies) in each direction [for example, Fig. 1(a) shows a 8 × 5 square network]. Throats are cylindrical channels of radius r t . The throat size is randomly distributed according to a given p.d.f. Similarly, the pore body size can be distributed according to a given p.d.f. with the constraint that the pore body size is equal to or greater than the largest throat to which it is connected. In what follows, several types of p.d.f. will be considered, namely, uniform, Gaussian, log normal, and bimodal. Also, it can be noted that the network coordination number, i.e., the number of throats connected to a pore body, is 4 for the 2D square network and 6 for the 3D cubic network.

IV. FLUID-FLUID IMMISCIBLE QUASISTATIC DISPLACEMENT SIMULATION

Modeling of immiscible two-phase flow displacements on a network has been the subject of many works and the reader is referred to [5,17,[START_REF] Sahimi | Flow and Transport in Porous Media and Fractured Rock: From Classical Methods to Modern Approaches[END_REF] for reviews. Here, the interest is on the displacement of a wetting fluid by a nonwetting one in the network. This process is known as drainage. In the quasistatic limit considered for the present paper, it was shown that drainage can be simulated using a simple algorithm referred to as the invasion percolation (IP) algorithm [START_REF] Lenormand | [END_REF][20][21]. This algorithm can be summarized as follows. A pressure difference, P ck , between the two phases is applied by imposing a pressure increment in the nonwetting phase in contact with one face of the network, identified as the inlet face. The opposite face is the outlet, whereas all the other faces of the network are kept impervious. At each pressure increment, all the throats containing a meniscus for which the capillary pressure threshold, P cth , as expressed by Eq. ( 2), is such that P ck P cth are identified. The throat for which (P ck -P cth )/P ck is maximum, hence corresponding to the interfacial throat of maximum diameter among the interfacial throats such that (P ck -P cth ) P ck > 0, is invaded together with the adjacent pore body if the latter was not already invaded in a previous step. If this pore body was occupied by the wetting fluid, new menisci are then positioned at the entrance of the throats occupied by the wetting fluid and connected to the newly invaded pore body. In other words, the list of interfacial throats is updated so as to remove the invaded throat and add the throats connected to the newly invaded pore body and occupied by the wetting fluid. This procedure is repeated until there are no unstable menisci remaining in the network, i.e., until there exists no throat for which (P ck -P cth ) P ck > 0. The resulting wetting -nonwetting fluid distribution corresponds to the equilibrium distribution at the specified value of P ck . For a sufficiently low pressure difference between the two fluids, the wetting phase displacement is not sufficient for the invading phase to reach the network outlet. The first equilibrium for which the invading fluid reaches the outlet corresponds to the breakthrough (BT) and is characterized by P ck = P cBT . While incrementing the displacing fluid pressure, the fluid distribution at BT is the first for which the nonwetting fluid can flow through the network, which means that Q NW = 0 for P ck < P cBT while Q NW > 0 for P ck P cBT . Note that the invasion rules described above correspond to a quasistatic drainage process without trapping. This means that an element, pore body, or throat, occupied by the wetting fluid, is assumed to be always connected to the outlet face, possibly through liquid films [22] present in the pore body and throat surface roughness and/or in a subnetwork of small crevices at the throat and pore body walls. Trapping rules could be implemented (see for instance [14]), but, for simplicity, they are not considered in what follows. It can be argued that neglecting trapping is not consistent with the considered cylindrical throat shape since this shape favors trapping phenomena as opposed to throats of a polygonal cross section which favor corner flows [22] and thus the wetting fluid hydraulic connectivity throughout the network. However, what is really important within the scope of the present paper is to consider the same flow physics in the reference network and the networks in the optimization procedure (see Sec. VI). The idea is to keep the simulations as simple as possible in the present effort. For this reason, trapping phenomena have been neglected without explicitly considering corner flows and thus throats of a polygonal cross section.

Each pressure increment provides a point in the graph of the retention curve, the fluid potential being the capillary pressure (i.e., the pressure difference between the nonwetting and wetting fluids) and the saturation being computed from the equilibrium fluid distribution. Next, the computation of the nonwetting fluid flow rate, Q NW , after breakthrough is performed as follows [5]. At each saturation equilibrium, a pressure difference, P, much smaller than the capillary pressure characterizing this equilibrium, is applied in the nonwetting fluid between the network inlet and outlet. For the resulting flow, the mass conservation equation is expressed at each invaded pore body i:

Q NW i, j = 0 ( 15
)
where Q NW i, j is the nonwetting fluid volume flow rate between pore body i and adjacent pore body j. If the throat between pore bodies i and j is occupied by the wetting fluid then Q NW i, j = 0. If the throat is occupied by the nonwetting fluid, then

Q NW i, j = πr 4 t i, j 8μ NW l i, j P i, j (16)
where r t i, j and l i, j are the radius and length of the throat linking the two pore bodies, respectively, while P i, j is the pressure drop in the nonwetting phase between the two neighboring pore bodies. Taking into account the pressure boundary conditions at the inlet (P inlet = P + P outlet ) and outlet (P outlet = Cte) and the zero flux condition imposed on the lateral faces of the network, Eqs. ( 15) and ( 16) lead to a linear system for the pressure field in the nonwetting fluid. This system is numerically solved using the conjugate gradient method. Once the pressure field is obtained, Q NW can be computed either at the outlet or inlet using Poiseuille's law applied to all nonwetting fluid throats connected to the outlet or inlet, respectively. In summary, the fluid-fluid distribution is computed, in a first step, neglecting the viscous pressure drop in each fluid phase on the grounds that capillary effects are dominant. For this given fluid-fluid distribution, the flow is computed in a second step as indicated above. This is carried out by sampling the capillary pressure, P c , with successive increasing values, P ck , the former step allowing the determination of the retention curve and the latter allowing the determination of the flow rate curve. Examples of fluid distribution in a 2D network are reported in Fig. 2.

Figure 3 shows the retention curve and the nonwetting fluid normalized flow rate Q n = Q NW Q max obtained with the above algorithm for a 50 × 50 square network with Gaussian distributions for both the TSD and PSD having mean radii of 30 and 90 μm and standard deviations of 5 and 10 μm, respectively. Here, Q max is the flow rate that would be obtained when the network is fully invaded and the same P as the one used to compute Q NW is applied, which means that Q max = A k μ P L , and thus, from Eq. ( 1), that Q n = k r .

V. ASSESSMENT OF TSD EXTRACTION FROM THE PARALLEL TUBES MODELS

As illustrated in Fig. 3, the PNM allows us to simulate the flow rate and retention (or capillary) curves with the decisive advantage that the TSD (and the PSD) of the medium is perfectly known a priori; it is referred to as the reference TSD in the following. Then, starting from the computed curves, the standard determination procedures can be used, i.e., Eqs. ( 6) and (13), in order to extract the TSD (or PSD) and compare the result with the reference TSD (or PSD). Since the throats control the displacement process, it is reasonable to expect that only the TSD can be determined. The TSD obtained from the parallel tubes model using either Eq. ( 6) or Eq. ( 13) is referred to as the simple model TSD, i.e., the SMTSD. A 2D square 50 × 50 network and a 3D cubic 20 × 20 × 20 network were considered for the comparison.

For each network, the Q n (P c ) and P c (S) curves were computed and the SMTSD extracted using Eqs. ( 6) and (13). In order to get smoother results, the extracted data, for a given type of network, were averaged over 100 realizations of the network. This corresponds to 100 realizations of each network type for which throat and pore body sizes are randomly assigned according to the given p.d.f. This means that the fluid-fluid displacement was computed 100 times for each case under consideration. Then the SMTSD was extracted from each flow rate curve and each retention curve, yielding 2 × 100 TSD. Then, for either the flow rate curves or the retention curves, the average TSD is computed from the 100 TSD. This average is presented and discussed below. A "case" is defined by the dimensionality of the network, i.e., square or cubic; the size of the network; the selected p.d.f., i.e., uniform or Gaussian; and the parameters of the p.d.f., namely, the minimum and maximum radii, r min and r max , respectively, for the uniform distribution or the mean size and standard deviation, r moy and σ , respectively, for the pore bodies and the throats in the case of a Gaussian p.d.f. The comparison between the reference TSD and the SMTSD for various cases when the retention curve information [Eq. ( 6)] is used is presented in Fig. 4 and Table I. As can be seen, the SMTSD does not match the reference TSD, in particular for the uniform distribution for which the uniform nature is not retrieved. For the 3D networks, it can be seen that the SM procedure introduces skewness towards large throat size. In brief, the simple model (SM) does not contain the necessary information to predict the reference TSD.

As illustrated in Fig. 5 and Table I, the SMTSD is somewhat less well estimated from the flow rate curve compared to the case where the capillary curve information is used.

In particular, it can be noticed that the larger throats are not identified.

A major shortcoming is that some throat sizes are not retrieved at all. In the case of the 2D network with a uniform distribution for instance, the throat sizes in the upper half of the distribution are not retrieved, in particular while employing the flow rate curve. Thus, only half of the throat sizes is actually identified. The same problem can be noticed with the 3D network, but to a lesser extent since the range of nonidentified throat sizes is narrower (but here again the larger throats are not identified). Insights on this result can be gained from percolation theory [START_REF] Stauffer | Introduction to Percolation Theory[END_REF]. In the percolation theory, a certain fraction, p, of throats are "active." This fraction of active throats (or "bonds" in the language of percolation theory) is distributed randomly. Suppose p is progressively increased starting from p = 0. For a particular value, p c , of this fraction, a percolating path (a path of interconnected active throats connecting the inlet to the outlet) forms for the first time in the network. The quantity p c is the percolation threshold of the network the value of which depends on the network [START_REF] Stauffer | Introduction to Percolation Theory[END_REF]. For a square network in bond percolation, p c = 0.5, whereas p c ≈ 0.25 for a cubic network. In the present context, an active throat is a throat that can be invaded, i.e., a throat such that P ck P cth , and the breakthrough corresponds to the percolation threshold since there is no nonwetting fluid path flowing through the network before breakthrough. The percolation threshold can thus be related to the TSD:

p c = ∞ r c f (r)dr = 1 - r c 0 f (r)dr = 1 -F (r c ) ( 17
)
where r c is the percolation radius, i.e., the radius given by P cth = P cBT , f being the reference TSD p.d.f., and F the corresponding cumulative distribution function [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF]. Thus, the threshold or percolation radius, r c , represents the critical radius above which the flow rate is nonzero. This means that, when the bundle of parallel tubes model is considered, all the throats of the network having radii greater than the threshold radius, which are hence invaded by the nonwetting fluid prior to the BT, cannot be determined by the parallel tubes model.

From Eq. ( 17), r c is given by

r c = r max (1 -p c ) + p c r min ( 18
)
for a uniform p.d.f., where r max and r min are the maximum and the minimum radii of the (reference) TSD, respectively. For a Gaussian TSD, the threshold radius is given by

r c = r moy + √ 2 σ erf -1 (1 -2p c ) (19)
where erf -1 is the inverse error function.

For a square lattice (p c = 0.5), the above relationships give r c = r max +r min 2 and r c = r moy for the uniform and Gaussian p.d.f., respectively. For a cubic lattice (p c ≈ 0.25), one obtains r c = 3 4 r max + 1 4 r min and r c = r moy + √ 2 2 σ for the uniform and Gaussian p.d.f., respectively. The values are reported in Fig. 5. The analytical threshold radii are perfectly consistent with the results shown in the TSD comparisons in Fig. 5. In summary, the parallel tubes model is not the appropriate representation able to identify the throat sizes greater than the critical radius corresponding to the network percolation threshold. The greater the percolation threshold, the larger the range of missing throats in the TSD determination with the parallel tubes model.

VI. BEYOND THE TSD EXTRACTION FROM THE PARALLEL TUBES MODEL

The observation above calls upon more sophisticated procedures to extract the TSD of a porous material from data available from a primary drainage process. As shown in the previous section, the standard and commonly used models to extract the TSD do not allow determining the reference TSD from the Q n (P c ) or P c (S) data. This unsatisfactory situation is explained by two reasons. First, for a given pressure step, all the newly invaded throats, regardless of their sizes, are assigned to a single throat size value corresponding to the imposed pressure difference related to the Young-Laplace equation (2). Second, when the determination of the TSD is based on the flow rate curve, all the radius values larger than the percolation radius cannot be determined since the flow rate is zero until breakthrough is reached, hence providing no information on all the throats of radius larger than the critical value. Thus, the question arises as to whether a better use of the Q n (P c ) and P c (S) curves can be performed in order to predict the reference TSD. An attempt in this direction is presented in the following.

A. TSD determination as an inverse problem

In essence, the problem to be solved is an inverse problem relying on the following question: knowing the Q n (P c ) and P c (S) curves, can the TSD of the corresponding medium be reliably predicted? In this procedure, the direct problem is of course to predict the Q n (P c ) and P c (S) curves of the porous medium under consideration. In addition, the inverse problem is solved using again the PNM and the direct problem is solved as in the previous section using the IP algorithm in the absence of trapping, considering the pressure and no flux boundary conditions already mentioned in Sec. IV. Actually, a somewhat simplified version of this inverse problem is tackled, assuming that the pore network structure is known. This means that the minimization procedure to fit Q n (P c ) and P c (S) is carried out for a given type of network, i.e., square or cubic for instance. In other words, for a given type of network (a 50 × 50 2D square network for instance), the objective is to correctly identify the TSD. The assumption that the pore network structure is known a priori is further discussed in Sec. VII below. To build and test the procedure, a particular realization of the network, with specified p.d.f. parameters, is used as the reference network the TSD of which is hence known. It shall be viewed as the porous medium to be characterized for which the experimental measurement of Q n (P c ) and P c (S) would have been performed in practice and the TSD of which is to be determined.

The inverse problem is solved using an optimization method based on a searching algorithm which is a hybrid version of the genetic and hill climbing algorithms [START_REF] Su | [END_REF]. The genetic algorithm [26,27] is a searching algorithm based on the biological evolution laws including reproduction, mutation, and natural selection. A population of individuals is first defined. In our case, the individuals are realizations of the pore network. Every individual has genes. A gene in a pore network is either a throat or a pore body characterized by its size. The algorithm is relying on a fitness function the value of which is estimated from a direct solution for each individual. The fitness function, denoted fit(r t ) in its generic form, is an objective function comparing the curve obtained by direct simulation for the individual to the reference one, r t being a vector the components of which are the individual throat radii. The genetic algorithm aims at minimizing this fitness function [26].

The overall algorithm can be summarized by the following steps.

The first step consists in generating a series of N 0 parent realizations (N 0 = 16 in this paper) to which evolutionary laws and natural selection are applied. In the most general case, both the throat and pore body sizes vary from one element (pore body or throat) to the other in the network. Referring for instance to the computation of the normalized flow rate, it is clear from Eq. ( 16) that the local volume flow rate through a throat depends not only on the throat radius, but also on the adjacent pore bodies radii. Indeed, by construction, the length of a throat l i, j is expressed as l i, j = ar pi, jr pi+1, j (20) where a is the lattice spacing and r pi, j and r pi+1, j are the radii of the two adjacent pore bodies. As a result, the Q n (P c ) and P c (S) curves depend not only on the TSD but also on the PSD. Thus, the pore body size should be kept as a variable in the optimization process in order to keep generality. However, as already mentioned, the drainage process is essentially sensitive to the TSD and the impact of the PSD on the flow rate can be reasonably neglected. Thus, the pore body size is taken as constant (equal to 80 μm) in what follows. Consequently, generation of individuals is carried out by only varying the throat sizes. The minimum and maximum radii of the TSD are determined from the invasion percolation procedure carried out on the reference network. The minimum throat radius, r min 0 , is given by Eq. ( 2) at the corresponding pressure step when the saturation becomes equal to zero. Conversely, the maximum throat radius, r max 0 , cannot be as exactly determined as r min 0 . It is approximately specified by multiplying by a factor 1.3 the radius corresponding to the first pressure step at which the saturation is no longer equal to 1.

In order to introduce diversity in the initial population, individuals are generated within N 0b class intervals. The ith interval has also a uniform p.d.f. with minimum and maximum radii, r min i and r max i , respectively, given by r min i = (N 0b +i)r min 0 +(N 0b -i)r max 0 2N 0b and r max i = (N 0b -i)r min 0 +(N 0b +i)r max 0 2N 0b

. The number of individuals generated in the ith interval is N 0i so that N 0 = N 0b i=1 N 0i . Here, N 0b = 4 and N 0i = 4. The normalized nonwetting fluid flow rate and capillary pressure retention curves, together with the corresponding fitness function of interest, are computed for each individual parent. After this initialization procedure, the genetic algorithm is executed according to the following steps.

(1) The natural selection is applied following a tournament method [26]. This is achieved by randomly choosing a subset of N 1 individuals among the generation of individuals (here, N 1 = 5) and the two ones featuring the best fitness values are selected.

(2) Two evolutionary operators, namely, crossover and mutation, are applied to the two selected parent individuals. Crossover is the process by which the two selected parent individuals give birth to a child realization by randomly exchanging throat sizes (genes). Mutation is an operation applied separately to both selected parent realizations. It consists in altering genes (the sizes of randomly chosen throats). Mutations are generally supposed to happen rarely. In the present case, the mutation rate is taken equal to 0.01. In our case this means that in a 50 × 50 network, for instance, 12 throats are selected randomly in each direction. Then the sizes of these 24 throats are specified randomly between r min and r max according a uniform p.d.f.

(3) The Q n (P c ) and P c (S) curves and the fitness values for the two children are computed and, among the four individuals (the two parents and the two children), the two individuals having the best fitness values are selected. They are kept to form the new generation.

(4) Steps 1-3 are repeated until there are no more individuals from the current parent population and a new generation of N 0 individuals is formed.

(5) The process made of steps 1-4 is iterated until convergence is reached. This means that the generation of N 0 networks obtained at the end of steps 1-4 forms the new parent generation when convergence has not been reached. By convergence, it is meant that the fitness of all the individuals from one generation to the following one does not change by more than about 0.1%.

For all the cases under study in this paper, 200 generations (iterations) were enough to reach convergence at which the solution is estimated to be close to the global solution. Then a new step is performed using the hill climbing algorithm.

The solution obtained from the genetic algorithm is global in the sense that the population of all the individuals is optimized by adjusting their overall genes content. When considering a network with its given TSD, the genetic algorithm, which operates on randomly selected throats, keeps in fact this TSD. A network with a Gaussian TSD, for example, is such that throat radii close to the mean radius are more likely to be selected than others. This limits the population diversity. The objective of the hill climbing algorithm, which, like the genetic algorithm, is a searching algorithm, is to gradually modify the throat size p.d.f. within a single individual so as to converge toward the TSD of the reference network. To this purpose, the hill climbing algorithm is operated on the best fitted individual in the population once the genetic algorithm ends. It attempts to find a better solution by adjusting only a single variable among all the variables. By contrast to the genetic algorithm where the throat radii are adjusted, the hill climbing algorithm [START_REF] Su | [END_REF] aims at the adjustment of the throat densities. Once the lower and upper bounds of the TSD are specified at the starting point of the optimization process, the distribution can be divided into class intervals, or bins. Then the density of throat sizes belonging to each bin can be determined, i.e., the fraction of throats having a size within the radius interval corresponding to the bin. The hill climbing algorithm aims at determining the optimal density of each bin. To summarize, the genetic algorithm is a multidirectional search algorithm operating on radii belonging to different distribution bins, whereas the hill climbing algorithm is a unidirectional algorithm that operates only on radii belonging to a pair of defined bins. Actually, as the sum of all distribution bin densities must be equal to 1, the hill climbing algorithm operates on two bin densities by increasing one (the positive bin) and decreasing the other (the negative bin) with the same amount. Thus, two bins must be selected. In our case, the two bins of the best fitted individual resulting from the genetic algorithm are randomly chosen. A relatively small randomly chosen number of throats are interchanged between the negative and positive bins. If the fitness function is improved, the modified individual is retained and the hill climbing process is repeated until the fitness function does not improve. When it is not improved, another test is performed by randomly choosing a different couple of bins. The fitness function for the hill climbing algorithm is equivalent to the one considered for the genetic algorithm and is denoted by fit(d) in its generic form, where d is the distribution density vector. The difference between the genetic and hill climbing algorithms lies in the variables to control for minimizing the objective function. Fitness functions are detailed in the following section. The optimization algorithm is summarized in the chart depicted in Fig. 6.

B. Fitness functions

The direct problem solution provides the Q n (P c ) and P c (S) curves of the reference network, i.e., the one for which the TSD from either the normalized flow rate curve or the capillary pressure curve data, or both, is to be determined. As indicated above, the network is a 2D, 50 × 50 square network, with a lattice spacing of 350 μm and a chosen TSD as will be further presented in the results section. The initial N 0i individuals (also 2D, 50 × 50 square networks with a lattice spacing of 350 μm) are generated in the N 0b class intervals having a uniform TSD with their own minimum and maximum radii as mentioned above.

For the same pressure steps as the ones considered for the reference Q n (P c ) and P c (S) curves, the optimization process, following the genetic and hill climbing algorithms, is performed. Again, this implies to solve the direct problem for every individual considered in the optimization process. The TSD at each step of the optimization procedure is compared to the reference one by evaluating the fitness function, which is defined as the cumulative error function. If the optimization is based on the normalized flow rate (21) where n is the total number of pressure steps; Q n-ref j and Q n-ind j are, respectively, the reference and individual normalized flow rates corresponding to the jth pressure increment; and ω = r t , d depending on whether the algorithm is at the genetic or hill climbing stage. Again, flow rates are considered so as to use the same convention as for the FFDP [10]. However, it should be clear that Q n = k r in the quasistatic limit under consideration, where k r is the nonwetting fluid relative permeability.

Q n = Q NW Q max , this fitness function is given by fit Q n (ω) = n j=1 |Q n-ref j -Q n-ind j (ω)|, ω = r t , d
Similarly, if the optimization is based on the capillary pressure curve, the fitness function is defined as (22) where S ref j and S ind j are the reference and individual saturations at the jth pressure increment, respectively. As can be inferred from Eq. ( 22), where ε ref and ε ind are, respectively, the reference and individual porosities, the porosity is actually also considered as a variable to optimize. By definition, εS is computed as

fit εP c (ω) = n j=1 |ε ref S ref j -ε ind S ind j (ω)|, ω = r t , d
εS = V noninv V total ( 23
)
where V noninv is the volume of the pore space not invaded at the j th step and V total is the total volume of the porous medium. The wetting fluid volume fraction, εS, does not involve the total pore space volume, V pores , but only the volume of the noninvaded pore bodies and throats as a variable to adjust. This is preferable to the saturation as the latter depends on both V noninv and V pores that can have different impact on the fitness function. Actually, the variable that is adjusted is the porosity when the saturation is equal to 1; otherwise, it is the saturation function.

When both the flow rate and the capillary pressure curves are considered, then the fitness function is defined as

fit Q n -εP c (ω) = n j=1 |Q n-ref j -Q n-ind j (ω)| + f p n j=1 |ε ref S ref j -ε indk S ind j (ω)|, ω = r t , d (24)
where f p is a penalization factor that is taken as constant and computed as

f p = 1 ε ref . ( 25
)
The idea with the weighted definition in Eq. ( 24) is actually to give more importance to the capillary pressure curve in the optimization process.

C. Results

Results are first discussed for a reference network with a Gaussian TSD having a mean radius of 40 μm and a standard deviation of 10 μm. This reference TSD is shown in Fig. 7(a). As mentioned in Sec. VIB, the initial population that must be generated at the beginning of the optimization process is composed of realizations with uniform distributions, thus with a type of p.d.f. different from the Gaussian one of the reference network.

Results on optimized TSDs when the target data are only the capillary pressure curve [Eq. ( 22)] (left column), only the normalized flow rate curve [Eq. ( 21)] (middle column), and both the retention and normalized flow rate curves [Eq. ( 24)] (right column) are represented in Fig. 8. Corresponding quantitative results are reported in Table II.

As can be seen from Fig. 8 and Table II, using the capillary pressure leads to a better result than when only the normalized flow rate curve is employed. Not surprisingly, the best result is obtained when both the retention curve and the normalized flow rate curve are used in the optimization process. Indeed, comparing the reference TSD in Fig. 7(a) to the optimized one obtained with this procedure in Fig. 7(c) shows a very good agreement. This is also confirmed by the results in Table II. As illustrated in Fig. 8 (right column), the retention and normalized flow rate curves resulting from the optimization procedure reproduce very well the corresponding curves of the reference case. A similar test was performed when the reference TSD is uniform [Fig. 7(b)]. The results are represented in Fig. 8 and quantitative data in Table III. Here again, the TSD obtained considering the capillary pressure curve only is closer to the reference one than when the optimization is based on the normalized flow rate curve only. The best result is again obtained when both the capillary pressure and the normalized flow rate curves are used in the optimization procedure. However, as illustrated in Figs. 7(b) and 7(d), the optimization is less performant than in the case of a Gaussian TSD for the reference medium. This discrepancy can be explained by the uncertainty on the distribution upper bound associated to the difficulty of accurately determining the maximum throat radius, r max 0 , as mentioned in Sec. VIA. This inaccuracy leads us to introduce some throats greater than the maximum radius (50 μm) in the reference distribution during the optimization procedure. TABLE II. Comparison between the reference and optimized TSDs when the reference TSD is the Gaussian TSD depicted in Fig. 7(a). This effect is of course much less sensitive in the case of distributions of Gaussian or log-normal types for which the relative density of throats close to the distribution upper bound is much lower than in a uniform distribution.

P c (S)
As also shown in Fig. 9, the capillary and normalized flow rate curves of the optimized network almost perfectly match the ones of the reference network. However, the match for the capillary curve is not satisfactory for the lower values of the capillary pressure when only the normalized flow rate data are used in the optimization procedure. This also holds in the case of the Gaussian reference TSD (see Fig. 8). This is an additional illustration of the impact of the distribution upper bound specification and also of the percolation threshold issue discussed in Sec. V. Two additional tests were performed considering a reference log-normal distribution [as depicted in Fig. 10(a)] and the expected more challenging bimodal distribution depicted in Fig. 10(b). For these cases, only the full optimization proce-dure, i.e., the one based on both P c (S) and Q n (P c ), was used. As can be seen from Fig. 10 and Table IV, results are quite satisfactory. Here also, as shown in Fig. 11, the match between the reference capillary pressure and the normalized flow rate TABLE IV. Comparison between the reference and optimized TSDs using both the capillary pressure and normalized flow rate curves as target data, where the reference TSD is the log-normal one depicted in Fig. 10(a) or the bimodal one depicted in Fig. 10(b curves and the curves resulting from the optimization procedure is excellent. Finally, among the four reference TSDs considered here, results obtained for the uniform one are the less accurate, mainly because of the uncertainty associated with the distribution upper bound specification.

VII. DISCUSSION

To further illustrate the interest of the optimization approach compared to the conventional approaches based on the consideration of the parallel tubes model and either the normalized flow rate curve [Eq. (13)] or the Young-Laplace equation and the retention curve [Eq. (6)], the simple model (i.e., the parallel tubes model) was applied to the bimodal reference network [Fig. 10(b)]. The results are represented in Fig. 12 and reported in Table V.

As can be seen, the bimodal nature of the distribution cannot be identified using the simple model. The approach based on the normalized flow rate curve and the parallel tubes model, i.e., Eq. ( 13), only captures the throat sizes below the percolation threshold (as discussed in Sec. V). As a result, only throats in the range of the lower peak distribution are retrieved. When the capillary pressure curve combined with the Young-Laplace equation is used, i.e., Eq. ( 6), throats over a larger range of size are identified but the bimodal nature of the distribution is not at all retrieved. However, as indicated in Table V, the procedure based on the capillary pressure curve, i.e., Eq. ( 6), leads to a satisfactory estimate of the throat size mean value. Actually, as can be seen also from Table I, the use of the capillary pressure curve information in conjunction with the Young-Laplace equation for the parallel tubes model leads to the same result, i.e., the throat size mean value is reasonably well predicted. Note from Tables II and III that it is also the case with the optimization procedure: the throat size mean value is well predicted when the capillary pressure curve information is used [and also when both the P c (S) and Q n (P c ) information are used as can be seen from Tables II-V]. An explanation of this fact, i.e., the fact that the use of the retention curve in conjunction with the Young-Laplace equation and the parallel tubes model leads to a good estimate of the mean throat size, is proposed in the Appendix. 12. Reference bimodal TSD (a), SMTSD identified using Eq. ( 13), i.e., the normalized flow rate curve and the parallel tubes model (b), and SMTSD identified using Eq. ( 6), i.e., the capillary pressure curve and the parallel tubes model (c).

For perspective, we offer the following remarks on the above results. Although the TSD identification based on the optimization procedure clearly leads to much better results than the simple parallel tubes model, several points are worth being recalled or discussed. First, trapping phenomena were ignored. It would be interesting to investigate the possible impact of this phenomenon on the TSD determination. Second, the optimization procedure has been tested considering only the throat size variability. In other words, the pore body size variability (the PSD) was not considered. In this respect, it would be interesting to explore the ability of the optimization procedure to capture both the PSD and the TSD. Third, the optimization procedure has been performed for given pore network structures, i.e., square or cubic in our simulations. It can certainly be argued that, in general, the network structure is unknown a priori. In addition, the pore network that can represent a real porous medium is rarely structured but rather unstructured. This means that the distance between two adjacent pore bodies may vary, as well as the local coordination number, i.e., the number of pore bodies to which a given pore body is directly connected. Therefore, an interesting question to address will be to test the impact of the structure. In other words, if cubic networks are used, for instance, in the optimization procedure with reference data obtained from a medium the 3D pore network of which is obviously not cubic, will the procedure still be relevant, or at least provide significantly better results than those obtained using the much simpler procedures making use of Eq. ( 6) or Eq. ( 13)? Since the microstructures of porous media can be quite diverse, it can be surmised that the quality of the results will also be quite diverse depending on the medium pore network structure and how close it is from a cubic network. In principle, the optimization procedure can of course be used for pore networks different from cubic. However, this may be cumbersome, keeping in mind that many realizations of the network must be generated. Here, the use of stochastic network generators [28] might be an option. From a practical point of view, a simpler option may be to create a catalog of various pore network structures (cubic, Voronoi, etc.) and choose the one which seems the closest to the structure of the porous medium to be characterized. This choice could be guided by images from the microstructure or partial information on this structure. However, when detailed digital 3D images of the microstructure are available, it is probably simpler to try to extract the PSD and/or the TSD directly from the images [8,9]. Nevertheless, it would be interesting to test whether the use of the Q n (P c ) and P c (S) curves with the optimization procedure might help characterize the microstructure by testing different pore network structures available in the catalog. Another type of situation is when the throat or pore body sizes change due to some coupled phenomena (particle trapping, dissolution or precipitation, or biofilm development). Using the optimization procedure to characterize the pore network evolution can be much simpler than repeated 3D imaging. In this case, one could use a stochastic pore network generator based on data obtained on initial images of the microstructure. Nevertheless, this requires determining the retention curves and the nonwetting fluid normalized flow rate (as well as the porosity and the permeability) repeatedly. Finally, the network size considered in this paper remains quite small, and this was mainly motivated by the interest which was focused on the characterization of thin porous media with only few pores over their thickness [START_REF] Prat | Thin Porous Media: Handbook of Porous Media[END_REF]. One can wonder whether increasing the network size might affect the results. Although tests in this TABLE V. Comparison of the throat size mean values of the reference TSD with those obtained using the parallel tubes model and normalized flow rate curve, Eq. ( 13), or the capillary pressure curve combined with the Young-Laplace equation, Eq. ( 6), for the reference bimodal TSD depicted in Fig. 10(b direction could be interesting, it can be noted, for instance, that the problem associated with the network percolation threshold, when the flow rate information is used, will still be present whatever the network size. As a consequence, the simple models [based on Eq. ( 6) or Eq. ( 13)] are not expected to perform better whatever the network size. Another issue not addressed in the present paper is the possible impact of finite-size effects (see, e.g., [START_REF] Larson | [END_REF]), when considering relatively small networks. The finite-size effects are of course more important in small networks and might hamper the simulations. It is assumed that the network considered in the present paper is sufficiently large for the results to be valid for much larger pore networks.

It should also be noted that the optimization procedure is significantly more intense in terms of computational resources and computer programming compared to the conventional procedures based on the parallel tubes model, which actually correspond to negligible computational times. For a 50 × 50 network, the computational time of one simulation corresponding to the computation of the Q n (P c ) and P c (S) curves considering 60 pressure steps is on average 3.66 s (using a single CPU 2.5-GHz Xeon processor). The genetic algorithm is run considering 200 generations. Each generation includes 16 individuals, hence involving 16 Q n (P c ) and 16 P c (S) curves simulations. Thus the computational time associated with one generation is 16 3.66 ≈ 58.56 s. The computational time of the genetic algorithm step is therefore 200 × 58.56 = 11 712 s (3 h, 16 min). The hill climbing algorithm step is performed over a single network but typically requires between 500 and 700 Q n (P c ) and P c (S) curves simulations, which correspond to a computational time of about 35 min. The total computational time for a 50 × 50 network is therefore around 4 h. Since the latter is relatively important even for a small network, the question arises as to whether the method can be realistically used in practice. In this respect, one can first notice that the code is written in MATLAB. A C or FORTRAN version would be significantly more efficient. The IP algorithm used to compute the P c (S) curve could be made more effective, (see for instance [31,32]). Some algorithm parallelization [a genetic algorithm, a linear system solver for the computation of the Q n (P c ) curve, etc.] is also expected to considerably reduce the computational time. Also, as mentioned previously, the porosimetry method considered in the paper is primarily developed for the characterization of thin porous media, like fuel cell porous components, [16], filtration membranes, [10], etc. The number of pores over the thickness of such media can be quite small, less than 10 for the porous component considered in [16] for instance. Thus, the network need not necessarily be large. Also, network size for the computations of properties over a representative elementary volume is typically of a few tens of pores [5]. In summary, the possibility of seriously improving the computational performances, compared to our current code, are sufficiently significant for considering that the computational performances should not be an issue for using the method in practice.

VIII. CONCLUSION

Pore network simulations were developed to assess the performance of two conventional and simple methods aiming at the determination of the TSD of a porous medium from macroscopic data such as the retention curve and/or the normalized flow rate curve. It was shown that using the normalized flow rate curve alone, which is roughly an information similar to the one used in the FFPD technique [10], is not satisfactory since all the throats greater than the critical throat corresponding to the network percolation threshold cannot be identified. On the whole, if a parallel tubes model is considered, the procedure based on the capillary pressure curve combined with the Young-Laplace equation leads to better determination of the mean radius than the procedure based on the normalized flow rate curve. However, all the invaded pores at a given pressure step in both procedures are assigned to one pore size corresponding to the specified capillary pressure regardless of their actual sizes. Thus, these simple procedures do not make any distinction between the TSD and the PSD. In order to alleviate the limitations of these conventional TSD identification methods, a procedure combining pore network simulations and an optimization technique based on a genetic algorithm and the hill climbing algorithm was developed. This algorithm leads to much better results compared to the simpler methods. Referring again to the FFPD, the paper also indicates that the knowledge of the nonwetting fluid normalized flow rate is not sufficient to determine the TSD. Using the porosity and the retention curve as additional information is necessary to obtain satisfactory results.

However, the optimization procedure is significantly more intense in terms of computational resources and computer programming compared to the conventional procedures. Efficient algorithms, programming languages, and computational procedures (parallelization) are desirable for practical use of the optimization procedure so as to perform the optimization over a sufficiently large network in a reasonably short time. Also, a somewhat simple situation has been considered to test the optimization procedure. Only the TSD was taken into account and no variability in the pore body size was envisaged. The network structure and size were the same in the reference network and the various networks generated during the optimization procedure. Certainly, more work is needed to fully explore the capability of the approach proposed in the present paper with the aim of solving the ill-posed problem consisting in identifying the TSD (and the PSD) from macroscopic data such as the normalized flow rate and capillary pressure curves. Nevertheless, results obtained here are very encouraging and suggest to further investigate the way opened by the methodology developed in this paper.

FIG. 1 .

 1 FIG. 1. Square (2D) (a) and cubic (3D) (b) pore networks.

FIG. 2 .FIG. 3 .

 23 FIG. 2. Examples of fluid-fluid distribution computed in a 30 × 30 2D network with a uniform TSD characterized by r min = 20 μm and r max = 60 μm. Invading fluid in blue enters the network at the bottom face: (a) P c = 2900 Pa, S = 0.82, (b) P c = 3250 Pa, S = 0.71, (c) P c = 3650 Pa, S = 0.55, (d) P c = 4250 Pa, S = 0.12.

FIG. 4 .

 4 FIG. 4. Comparison of the TSD obtained from the simple model (SMTSD) using the retention curve P c (S), i.e., Eq. (6), and the reference (real) TSD for various cases: (a) 2D 50 × 50 network, uniform p.d.f., r max = 40 μm, r min = 20 μm, (b) 2D 50 × 50 network, Gaussian p.d.f., r moy = 20 μm, σ = 2 μm, (c) 3D 20 × 20 × 20 network, uniform p.d.f., r max = 40 μm, r min = 20 μm, (d) 3D 20 × 20 × 20 network, Gaussian p.d.f., r moy = 20 μm, σ = 2 μm.

FIG. 5 .

 5 FIG. 5. Comparison of the TSD obtained from the parallel tubes model (SMTSD) and the reference (real) TSD for various cases using the flow rate information, Eq. (13). r c is the critical radius corresponding to the network percolation threshold: (a) 2D 50 × 50 network, uniform p.d.f., r max = 40 μm, r min = 20 μm, r c = 30 μm, (b) 2D 50 × 50 network, Gaussian p.d.f., r moy = 20 μm, σ = 2 μm, r c = 20 μm, (c) 3D 20 × 20 × 20 network, uniform p.d.f., r max = 40 μm, r min = 20 μm, r c = 35 μm, (d) 3D 20 × 20 × 20 network, Gaussian p.d.f., r moy = 20 μm, σ = 2 μm, r c = 21.4 μm.

FIG. 6 .

 6 FIG. 6. Optimization algorithm chart.

FIG. 7 .

 7 FIG. 7. Top: Reference Gaussian TSD (a) and reference uniform TSD (b) for a 50 × 50 square network. Bottom: TSD obtained from the optimization procedure using both the normalized flow rate and capillary pressure curve as target data. The TSDs in (c) and (d) are to be compared with the TSDs in (a) and (b), respectively.

PFIG. 8 .

 8 FIG.8. Results of optimization procedure considering the capillary pressure curve only as target data (left column), the normalized flow rate curve only (middle column), and both the normalized flow rate (equivalent to k r ) and capillary pressure curves (right column), when the reference p.d.f. is Gaussian [as depicted in Fig.7(a)].

FIG. 9 .

 9 FIG.9. Optimization results considering only the capillary pressure curve as target data (left column), only the normalized flow rate curve (middle column), and both the normalized flow rate and capillary pressure curves (right column). The reference p.d.f. is the uniform one reported in Fig.7(b).

FIG. 10 .

 10 FIG. 10. Results of the optimization procedure using both the normalized flow rate and capillary pressure curves as target data for a 50 × 50 square network having a reference log-normal TSD (a) and bimodal TSD (b). Figures in the left and right columns are for the log-normal and bimodal distributions, respectively. The optimized TSDs in (c) and (d) are to be compared with the reference TSDs in (a) and (b), respectively.

FIG. 11 .

 11 FIG. 11. Results of the optimization procedure using both the normalized flow rate and capillary pressure curves as target data for a 50 × 50 square network having a reference log-normal TSD (a) and bimodal TSD (b). Figures in the left and right columns are for the log-normal and bimodal distributions, respectively. Comparison of the reference and optimized capillary (or retention) curves (a), (b) and normalized flow rate (or relative permeability) curves (c), (d).

  FIG.12. Reference bimodal TSD (a), SMTSD identified using Eq. (13), i.e., the normalized flow rate curve and the parallel tubes model (b), and SMTSD identified using Eq. (6), i.e., the capillary pressure curve and the parallel tubes model (c).

TABLE I .

 I Comparison between the reference TSD and the SMTSD (parallel tubes) obtained from the Young-Laplace equation and the capillary pressure (or retention) curve, Eq. (6), or from the normalized flow rate curve, Eq. (13), for various cases.Parallel tubes model and P c (S) Parallel tubes model and Q n (S)

		Statistical moments (μm) Reference TSD SMTSD Relative error (%) SMTSD	Relative error (%)
	Uniform p.d.f. 2D network	Mean	30	28.84	3.88	24.14	19.53
		Standard deviation	5.75	2.78	48.3	2.72	52.64
	Gaussian p.d.f. 2D network	Mean	20	19.84	0.79	18.08	9.58
		Standard deviation	2	0.92	54.04	1.4	30.1
	Uniform p.d.f. 3D network	Mean	30	30.94	3.14	25.38	15.4
		Standard deviation	5.75	4.79	16.77	3.77	34.44
	Gaussian p.d.f. 3D network	Mean	20	20.69	3.45	18.6	7
		Standard deviation	2	1.19	40.5	1.59	20.3

TABLE III .

 III Comparison between the reference TSD and the obtained TSD, where the reference TSD is the uniform TSD depicted in Fig.7(b).

			P c (S) only	Q n (P c ) only	Both P c (S) and Q n (P c )
	Statistical	Reference	Optimized	Relative	Optimized	Relative	Optimize TSD	Relative
	moments (μm)	TSD	TSD	error (%)	TSD	error (%)	TSD	error (%)
	Mean	40	39.88	0.29%	43.67	9.18	39.93	0.18
	Standard deviation	5.77	6.36	10.22	10.72	85.79	5.97	3.47
	Skewness	0.01	0.91		0.97		0.42
	Kurtosis	1.8	5.06	181	2.99	66.1	2.93	62.78
				only	Q n (P c ) only	Both P c (S) and Q n (P c )
	Statistical	Reference	Optimized	Relative	Optimized	Relative	Optimized	Relative
	moments (μm)	TSD	TSD	error (%)	TSD	error (%)	TSD	error (%)
	Mean	40	39.57	1.08%	43.53	8.83	40.24	0.6
	Standard deviation	9.9	11.37	14.85	15.18	53.33	9.77	1.31
	Skewness	-0.01	0.06		0.45		0.043
	Kurtosis	3	3.275	9.16	2.57	14.33	2.77	7.67

  ).

			Log-normal TSD			Bimodal TSD
	Statistical moments (μm) Reference TSD Optimized TSD Relative error (%) Reference TSD Optimized TSD Relative error (%)
	Mean	39.84	39.85	0.003	35.3	35.09	0.6
	Standard deviation	9.83	9.8	0.3	16.61	17.07	2.8
	Skewness	0.63	0.58	7.94	0.06	0.245
	Kurtosis	3.48	3.33	4.31	1.71	2.01	17.54

  ).

			Parallel tubes model
	Reference TSD	Q n (P c )		P c (S)
	Mean (μm)	Mean (μm)	Relative error (%)	Mean (μm)	Relative error (%)
	35.3	19.6	44.48	35.3	0.3

ACKNOWLEDGMENT

This research was supported by the project "PEMFC-SUDOE" (Project No. SOE1/P1/E0293), which is cofinanced by the European Regional Development Fund in the framework of the Interreg Sudoe program.

APPENDIX: MEAN THROAT SIZE FROM THE RETENTION CURVE

Some explanation of the fact that the use of the retention curve in conjunction with the Young-Laplace equation always

where f n (r) is the normalized p.d.f. According to Eq. ( 8),

Combining Eqs. (A1) and (A2) with Eq. ( 5) yields

To relate S and r the Brooks and Corey representation [33] for the retention curve may be invoked, which in our case takes the form

where P M and λ are two parameters that must be fitted to represent the retention curve. Introducing r M = 2γ P M , where r M actually corresponds to the radius of the largest throat at the network inlet, and using Eq. ( 2) for θ = 0 leads us to express Eq. (A5) as

Substituting Eq. (A6) into Eq. (A4) leads after some calculations and the assumption that

Fitting the values of P M and λ for the networks featuring a Gaussian and uniform TSD explored in this paper (see Sec. V) and applying Eq. (A7) lead to the results presented in Table VI.

These results justify why the capillary curve contains the required information that can restitute a relevant value of the mean throat radius.


A André

A K Touré

D Stien

Véronique Eparvier

2,5-diketopiperazines mitigate the amount of advanced glycation end products accumulated with age in human dermal fibroblasts

2,5-diketopiperazines mitigate the amount of advanced glycation end-products accumulated with age in human dermal fibroblasts

A. André *, † , A. K. Touré † , D. Stien ‡ and V. Eparvier † Cell culture ; chemical analysis ; skin physiology ; 2,5-diketopiperazines ; molecular networks ; AGEs/CML clearance.

INTRODUCTION

Louis-Camille Maillard described the glycation reaction for the first time in 1912 [START_REF] Maillard | Action des acides aminés sur les sucres, formation des mélanoïdines par voie méthodique[END_REF]. This naturally occurring non-enzymatic reaction between reducing sugars and aminoacid residues of proteins leads to the formation of a complex and heterogenous group of compounds named advanced glycation end-products (AGEs).

AGEs accumulate in the whole body during the lifetime, both in the extra-and intracellular environments [START_REF] Vlassara | Advanced Glycation in Health and Disease: Role of the Modern Environment[END_REF]. Long lifespan proteins of the body, like collagen, elastin and hemoglobin are well-known targets of irreversible modifications caused by the glycation reaction [START_REF] Kasper | Age-related changes in cells and tissues due to advanced glycation end products (AGEs)[END_REF].

At the extracellular level, proteins crosslinking induces inhibition of cell growth, impaired cell adhesion, and tissue dysfunction [START_REF] Paul | The effect of advanced glycation end-product formation upon cell-matrix interactions[END_REF]. At the intracellular level, AGEs cause cell damage directly by altering protein functions, and indirectly by producing reactive oxygen species (ROS) [START_REF] Valencia | Advanced glycation end product ligands for the receptor for advanced glycation end products: biochemical characterization and formation kinetics[END_REF]. Excessive levels of ROS can then alter cell detoxification systems such as the proteasome, as well as the enzymes involved in cell repair [START_REF] Ott | Role of advanced glycation end products in cellular signaling[END_REF]. AGEs can also bind to the Receptor of Advanced Glycation End products (RAGE), triggering the release of inflammatory mediators by the NF-B pathway [START_REF] Liu | AGEs increased migration and inflammatory responses of adventitial fibroblasts via RAGE, MAPK and NF-kappaB pathways[END_REF]. Finally the inflammation and oxidation reactions triggered by AGEs will result in the formation of more AGEs in cells. These phenomena also reduce the antioxidant and repair potential of cells, which will eventually lead to dysfunction of tissues and organs [START_REF] Badenhorst | Cross-linking influences the impact of quantitative changes in myocardial collagen on cardiac stiffness and remodelling in hypertension in rats[END_REF]. Consequently, AGEs are considered to be the main agents responsible for many age-related diseases (3) such as diabetic vascular complications [START_REF] Goldin | Advanced glycation end products: sparking the development of diabetic vascular injury[END_REF], atherosclerosis [START_REF] Kume | Immunohistochemical and ultrastructural detection of advanced glycation end products in atherosclerotic lesions of human aorta with a novel specific monoclonal antibody[END_REF], or Alzheimer's disease [START_REF] Baynes | The role of AGEs in aging: causation or correlation[END_REF].

The glycation reaction does not spare the skin. AGEs have been shown to accumulate in both the dermis and epidermis [START_REF] Gkogkolou | Advanced glycation end products: Key players in skin aging? Dermatoendocrinol[END_REF][START_REF] Kawabata | The presence of N(ε)-(Carboxymethyl) lysine in the human epidermis[END_REF]. Pentosidine, a fluorescent AGE, and carboxymethyl-lysine, an AGE that forms protein adducts, are the two main final glycation products found in the skin [START_REF] Gkogkolou | Advanced glycation end products: Key players in skin aging? Dermatoendocrinol[END_REF]. The accumulation of AGEs in the intra-and extracellular skin environment causes, among other things, a stiffening of collagen, elastin or fibronectin fibers, as well as an increase in melanogenesis [START_REF] Lee | Advanced glycation end products (AGEs) promote melanogenesis through receptor for AGEs[END_REF][START_REF] Gkogkolou | Advanced glycation end products: Key players in skin aging? Dermatoendocrinol[END_REF]. The skin is thus less elastic and more sensitive to the various mechanical stimuli it endures. Also, the complexion is more yellow and less homogeneous, which inevitably affects its youthful appearance overtime [START_REF] Gkogkolou | Advanced glycation end products: Key players in skin aging? Dermatoendocrinol[END_REF][START_REF] Ogura | Dermal carbonyl modification is related to the yellowish color change of photo-aged Japanese facial skin[END_REF][START_REF] Reihsner | Alterations of biochemical and twodimensional biomechanical properties of human skin in diabetes mellitus as compared to effects of in vitro non-enzymatic glycation[END_REF], as seen in people with diabetes [START_REF] Lee | Advanced glycation end products (AGEs) promote melanogenesis through receptor for AGEs[END_REF][START_REF] Dyer | Accumulation of Maillard reaction products in skin collagen in diabetes and aging[END_REF]. Most compounds on the cosmetic market act by preventing or reversing the first step of the glycation reaction, i.e., the binding between proteins and sugars leading to the formation early glycation products, also called Amadori products. Looking at their mode of action, none of these ingredients really affect the reduction in the amount of final AGEs present in elderly cells. Therefore, the dermatological, cosmetic and therapeutic interest in discovering compounds that reduce the amount of AGEs accumulated in skin is important.

In this context, we screened a bank of pure compounds and microbial extracts to find suitable candidates to reduce both the total amount of AGEs, and the amount of pentosidine in human dermal fibroblasts [START_REF] André | A method to quantify intracellular glycation in dermal fibroblasts using liquid chromatography coupled to fluorescence detection -Application to the selection of deglycation compounds of dermatological interest[END_REF]. During screening, molecules belonging to the family of 2,5-diketopiperazines, and an ethyl acetate extract of Sphingobacterium sp.

(SNB-CN13) showed a reduction in AGE content of fibroblasts by more than 20%. In this article, we present for the first time the activity of naturally occurring 2,5-DKP on the elimination of carboxymethyl lysine AGE (CML) in human dermal fibroblasts and evaluate their potential as active ingredients in dermatology and cosmetics.

MATERIAL AND METHODS

General experimental procedures. Optical rotations were measured on an Anton Paar

MCP 300 polarimeter in a 100 mm-long 350 µL cell. UV spectra were recorded using a Perkin-Elmer Lamba 5 spectrophotometer. CD spectra were measured at 25°C on a JASCO J-810 spectropolarimeter. Nuclear Magnetic Resonance (NMR) spectra were recorded on a Bruker 500 MHz spectrometer or a Bruker 600 MHz spectrometer equipped with a 1 mm inverse detection probe. Chemical shifts (δ) are reported in ppm based on TMS signal, and coupling constants (J) are reported in Hertz. HRESIMS measurements were performed using a Waters Acquity UHPLC system with column bypass coupled to a Waters Micromass LCT premier Time-of-Flight mass spectrometer (Milford, Ma, USA) equipped with an electrospray interface (ESI). X-ray diffraction data for compound 11 were obtained at room temperature on a Rigaku Rapid II diffractometer equipped with a rotating anode mm007 HF generator and Osmic mirrors (Cu K radiation, λ = 1.54187 Å) using -scans. Data were indexed, integrated and scaled using CrystalClear [START_REF] Otwinowski | Processing of X-ray diffraction data collected in oscillation mode[END_REF]. They were also corrected for polarization, Lorentz and absorption effects (FS_ABSCOR). X-ray diffraction data for compound 20 were obtained at room temperature on a Rigaku XtaLabPro diffractometer equipped with a microfocus source (MicroMax003_Mo) and multilayer confocal mirrors (Mo K radiation, λ = 0.71073 Å). Data were indexed, integrated and scaled using CrysalisPro. They were also corrected for polarization, Lorentz and absorption effects (SCALE3 ABSPACK). For each compound, the structure was solved with the ShelXT (23) structure solution program using Direct Methods and refined with the ShelXL (24) refinement package using Least Squares minimization. All non-hydrogen atoms were refined with anisotropic displacement parameters and H atoms have been added geometrically and treated as riding on their parent atoms. Molecular graphics were computed with Ortep 3 [START_REF] Farrugia | WinGX and ORTEP for Windows: an update[END_REF].

Flash chromatography was performed on a Grace Reveleris system with dual UV and ELSD detection equipped with a 120 g C18 column. For UV-based experiments, effluents were monitored at 254 and 280 nm. Analytical and preparative HPLCs were conducted on a Gilson system equipped with a 322 pumping device, a GX-271 fraction collector, a 171 diode array detector and a prepELSII detector electrospray nebulizer. Columns used for these experiments included: a Phenomenex Luna C18 or PFP (2) 5 µm 4.6 x 250 mm analytical column, and Phenomenex Luna C18 or PFP (2) 5 µm 21.2 x 250 mm preparative column. The flow rate was set to 1 or 21 mL/min respectively, using a linear gradient of water mixed with an increasing proportion of acetonitrile. All solvents were HPL grade and supplemented with 0.1% of formic acid. The Mueller-Hinton (MH) was purchased from Conda. TritonX-100 was from Euromedex (Euromedex France, Souffelweyersheim, France). DMSO was obtained from Carlo Erba (Carlo Erba Reagents S.A.S, Val de Reuil, France). Tween-20 was obtained from Acros Organics (Acros Organics, Geel, Belgium). Bovine serum albumin (BSA) was obtained from Dutscher (Dominique Dutscher SAS, Brumath, France). Paraformaldehyde was obtained from Electron Microscopy Science (EMS, Hartfield, PA, U.S.A). Mouse anti-CarboxyMethyl-Lysine (CML) antibody (KH011) was obtained from Cosmo Bio Co (Cosmo Bio Co LTD, Tokyo, Japan). Goat anti-mouse secondary antibody Alexa-488 conjugate (A11017) was obtained from Lifetechnologies (ThermoFisher Scientific, Waltham, MA, U.S.A).

Mounting medium with DAPI was obtained from Vector (Vector Laboratories Inc, Burlingame, CA, U.S.A).

Isolation and identification of Sphingobacterium spp. SNB-CN13. Pieces of termite

nests were collected and cultured as described in our previous work [START_REF] Nirma | Antifungal agents from Pseudallescheria boydii SNB-CN73 isolated from a Nasutitermes sp. termite[END_REF]. The procedure adopted for identification of the strain was also describe in the same article [START_REF] Nirma | Antifungal agents from Pseudallescheria boydii SNB-CN73 isolated from a Nasutitermes sp. termite[END_REF], and the strain was saved at -80 °C at the Métabolites de micro-organismes symbiotiques -Bioinspiration laboratory, ICSN/CNRS, France. The bacteria rDNA sequence obtained was aligned with sequences from GenBank, NCBI (http://www.ncbi.nlm.nih.gov/genbank/ consulted on December 15, 2019), indicating that SNB-CN13 was closely related to Sphigobacterium species. The sequence was deposited on the GenBank with the accession number FJ959369.1.

Culture of bacterium, extraction and isolation of compounds. Sphingobacterium sp.

was cultivated on solid Mueller-Hinton medium at 28°C for one week on 101 Petri dishes that were 14 cm in diameter (total 1.4 m 2 ). The contents of the Petri dishes were cut into small pieces, transferred into a large container and macerated with ethyl acetate for 24 hours. The organic solvent was collected via filtration, washed with water in a separating funnel, and evaporated to dryness under reduce pressure to yield 3.24 g (2.3 g/m 2 ).

A portion of the ethyl acetate extract (2.6 g) was fractionated using flash chromatography using a gradient of water mixed with increased proportions of acetonitrile (H2O:CH3CN 95:5 -80:20 -50:50 -20-80 -0:100) at 80 mL/min using a C18 120 g column. Five fractions were then collected based on the solvent steps described above. During fractionation, a precipitate was formed in a tube of fraction F2, corresponding to compound 8 (2 mg). Fractions F1 (18.7 mg) and F2 (191 mg) retained biological activity and were subjected to further fractionation. Fraction F1 was purified by analytical HPLC (Luna PFP (2), mobile phase H20 + 0.1% FA / CH3OH + 0.1% FA, isocratic elution 95:5 over 35 minutes, flow rate 0.6 mL/min) to yield compound 18 (1.5 mg, RT = 8.2 min) and 15 (2.6 mg, RT = 26.5 min). Fraction F2 was purified by preparative HPLC (Luna PFP (2), mobile phase H20 + 0.1% FA / CH3CN + 0.1% FA, isocratic elution 88:12 over 30 minutes, flow rate 21 mL/min) to afford compound 6 All chemical descriptions of the isolated compounds (NMR, HREIMS) are available in the additional material section of this publication.

Cell preparation method.

The skin explant of a 53-year-old donor comes from abdominal surgery and is prepared (hypodermis eliminated) and provided by the BioEC Laboratory (Longjumeaux, France).

The explant (epidermis + dermis) is then cut into discs of about 1 cm in diameter. The discs are placed in a T75 culture flange, the dermis in contact with the plastic of the flange. After a few minutes of rest the discs adhere to the culture flange. Culture medium is added (10 mL of DMEM, 20% FBS, 1% antibiotic/antifungal medium) so that the skin discs are not completely covered.

The flask is then placed in the oven at 37°C, 5% CO2, and the culture medium is changed every 3 days. About a week later, a heterogenous culture of fibroblasts and keratinocytes is observed in the mother box. When the cells have proliferated over the entire surface of the flask, the cells are detached with a trypsin solution and stored in liquid nitrogen. During the first culture starting again from the cryotubes, only the fibroblasts colonize the flask, in DMEM medium with 10% FBS.

Before and during each test, cells are stained with trypan blue and counted using an automatic cell counter TC-20 (Bio-Rad, Marnes-la-Coquette, France). The cell viability of the counted sample is automatically analysed by the device according to the trypan blue exclusion principle. The percentage of viability of treated cells has always been above 95%.

Biological activity -Immunofluorescence of CML. Dermal fibroblasts were plated at 5000 cells/cm 2 on glass coverslips set into a 24-wells-plate and cultured in DMEM supplemented with 10 % FBS for 24 hours. On day 2, treatment with the molecules was done for 24 hours, at the concentration of 1 µg/mL. On day 3, medium containing the treatment was removed, cells were washed with PBS, and fresh DMEM + 10 % FBS medium was added in wells. On day 4, the fibroblasts were fixed with 2 % paraformaldehyde solution for 12 minutes at room temperature, permeabilized with 0,2% TritonX-100 solution for 10 minutes and finally blocked with 5 % BSA, 0,2 % Tween-20 solution for further 10 minutes. Cells were then incubated with the primary antibody anti-CML (mouse) at a concentration of 7 µg/mL at 4°C overnight, followed by incubation with Alexa-488 anti-mouse secondary antibody diluted at 1:200 for 2 hours at room temperature. After washing steps, coverslips were mounted onto microscope slides using a DAPI containing mounting medium and observed under a confocal photonic microscope (Spinning-Disk microscope, Roper/Nikon). Five images were recorded by coverslips (ten pictures per treatment condition). All pictures were taken with the same exposure time. Fluorescence quantification of images recorded was done using ImageJ software.

UPLC-MS/MS analysis of SNB-CN13 extract

Extract was diluted at 2 mg/mL in MeOH. Metabolome profiling was conducted with a Thermo UHPLC-HRMS system. Analyses of extract (1.0 µL injected, 2 µg on column) was performed in electrospray positive ionisation mode in the 133.4-2000 Da range in centroid mode. The mass detector was an Orbitrap MS/MS FT Q-Exactive focus mass spectrometer. The analysis was conducted in FullMS data dependent MS 2 mode. In FullMS, resolution was set to 70,000 and AGC target was 3.10 6 . In MS2, resolution was 17,500, AGC target 10 5 , isolation window 0.4 Da, normalized collision energy 30, with 20 s dynamic exclusion. UHPLC column was a Phenomenex Luna Omega polar C18 150×2.1 mm, 1.6 µm. The column temperature was set to 42 °C, and the flow rate was 0.5 mL/min. The solvent system was a mixture of water (A) with increasing proportions of acetonitrile (B), both solvents modified with 0.1% formic acid. The gradient was as follows: 2 % B 3 min before injection, then from 1 min to 13 min, a linear increase of B up to 100 % followed by 100 % B for 5 min. The flow was discarded (not injected into the mass spectrometer) before injection and up to 1 min after injection.

MZmine 2 Data preprocessing Parameters

Raw files were directly processed using MZmine 2.32. The mass detection was realized keeping the noise at 20. The chromatogram building was achieved using a minimum time span of 0.01 min, minimum height of 1.0E7 and m/z tolerance of 0.0015 or 5 ppm.

The local minimum search deconvolution algorithm was used with the following settings: chromatographic threshold = 1 %, minimum retention time range = 0.05 min, minimum relative height = 30 %, minimum absolute height = 5.0E5, minimum ratio of peak top/edge = 1, and peak duration range = 0.01 -5.0 min. Chromatograms were deisotoped using the isotopic peaks grouper algorithm with an m/z tolerance of 0.0015 or 5 ppm, and an RT tolerance of 0.1 min. Peak alignment was performed using the RANSAC aligner method with the following parameters: m/z tolerance = 0.0015 or 5 ppm, absolute RT tolerance = 0.1 min, absolute RT tolerance after correction = 0.5 min, threshold value = 1). The peak list was gap filled with the peak finder module (intensity tolerance at 50%, m/z tolerance at 0.0015 or 5 ppm and absolute RT tolerance of 0.1 min).

Molecular Network analysis

A molecular network was created using the online workflow at Global Natural Products Social Networking (GNPS) (http://gnps.ucsd.edu) [START_REF] Wang | Sharing and community curation of mass spectrometry data with Global Natural Products Social Molecular Networking[END_REF]. The data was filtered by removing all MS/MS peaks within +/-17 Da of the precursor m/z. MS/MS spectra were window filtered by choosing only the top 6 peaks in the +/-50 Da window throughout the spectrum. Precursor ion mass tolerance and fragment ion mass tolerance were set at 0.01 Da. Further, consensus spectra that contained less than 1 spectrum were discarded.

A network was then created where edges were filtered to have a cosine score above 0.7 and more than 4 matched peaks. Further edges between two nodes were kept in the network if and only if each of the nodes appeared in each other's respective top 10 most similar nodes. The spectra in the network were then searched against GNPS' spectral libraries. The library spectra were filtered in the same manner as the input data. All matches kept between network spectra and library spectra were required to have a score above 0.7 and at least 4 matched peaks. Finally, molecular network was visualized using Cytoscape 3.4.0.

General procedure for the synthesis of symmetric 2,5-diketopiperazines. Thionyl chloride (SOCl2, 1.67 mL, 22.9 mmol) was added dropwise to a solution of amino acid (15.19 mmol) in methanol (MeOH, 25 mL) at 0 °C under argon. The mixture was stirred for 19 h at room temperature. The solution was evaporated under reduced pressure.

The residue was cooled at 0 °C and mixed with an aqueous solution of K2CO3 (50 %, 2.5 mL) for 15 minutes. The solution was then extracted with tert-butyl methyl ether (TBME), and the organic phase was dried with Na2SO4 before being evaporated under reduced pressure. The oily residue was then heated for 17 h between 80 and 105 °C depending on the amino acid used (Leu: 80 °C; Met: 80 °C; Pro: 105 °C). Then TBME (2 x 3 mL) was added to the reaction mixture and the precipitate was collected by filtration after cooling to afford the corresponding symmetric 2,5-diketopiperazine. The NRM, HRMS and crystallographic descriptions of the synthetized molecules can be found in the supplementary material.

RESULTS AND DISCUSSION

A preliminary examination of a bank of pure compounds and microbial extracts revealed that 2,5-diketopiperazines and an ethyl acetate extract of Sphingobacterium sp. SNB-CN13 reduced the AGE content of dermal fibroblasts by 28 and 23 %, respectively [START_REF] André | A method to quantify intracellular glycation in dermal fibroblasts using liquid chromatography coupled to fluorescence detection -Application to the selection of deglycation compounds of dermatological interest[END_REF].

The working hypothesis of this work was that the activity of the Sphingobacterium sp.

extract could be attributed to the presence of DKPs.

First, the SNB-CN13 extract was analyzed by LC-MS 2 . The collected data were used to build a molecular network to visualize the chemical diversity of the Sphingobacterium sp. extract (Fig. 1). Based on comparison with spectral databases and spectral similarities, we have identified two groups of related nodes containing DKPs in the network. By examining in detail the clustered MS 2 spectra in these nodes, we have identified successive characteristic losses of CO (-27.9 Da), NH3 (-45.0 Da), and CO again (-73.0 Da) from the parent ion, as well as characteristic fragment ions of amino acids residues formed after the cleavage of the amide bonds. A typical fragmentation scheme is provided in the supporting information (Scheme S1). Overall, this information confirmed the identification of the parent ions as pseudomolecular diketopiperazines ions [START_REF] Xing | Rapid screening for cyclo-dopa and diketopiperazine alkaloids in crude extracts of Portulaca oleracea L. using liquid chromatography/tandem mass spectrometry[END_REF]. Overall, 26 presumed DKPs were found in the crude extract of Sphingobacterium sp. (Table S1). All compounds were identified by comparison of their spectroscopic data and physicochemical properties with the reported in the literatures, and determined to be identical.

Following these results, the synthesis of diketopiperazines was studied. In order to develop a cosmetic ingredient for the organic market, the development strategy turned to a synthesis protocol based on the principles of green chemistry. In this process, amino acids were esterified with MeOH/HCl, and then condensed by heating after deprotonation of the amino group (Fig. 3) (44). Using this protocol, the three symmetrical DKPs cyclo(L-leucine-L-leucine) [START_REF] Lee | Advanced glycation end products (AGEs) promote melanogenesis through receptor for AGEs[END_REF], cyclo(L-methionine-L-methionine) [START_REF] Dyer | Accumulation of Maillard reaction products in skin collagen in diabetes and aging[END_REF], cyclo(L-proline-L-proline) [START_REF] André | A method to quantify intracellular glycation in dermal fibroblasts using liquid chromatography coupled to fluorescence detection -Application to the selection of deglycation compounds of dermatological interest[END_REF] were obtained without epimerization (Fig. 4). I. Tyr-L-Tyr) (7) to 28 % for cyclo(L-Leu-L-Leu) [START_REF] Lee | Advanced glycation end products (AGEs) promote melanogenesis through receptor for AGEs[END_REF]. Cyclo(L-Met-L-Met) [START_REF] Dyer | Accumulation of Maillard reaction products in skin collagen in diabetes and aging[END_REF].

Cyclo(Dha-L-Ile) (2), cyclo(L-Leu-L-Glu) (4) and cyclo(L-Pro-L-Met) (9) were inactive. It should also be noted that the symmetrical synthetic DKPs cyclo(L-Leu-L-Leu) [START_REF] Lee | Advanced glycation end products (AGEs) promote melanogenesis through receptor for AGEs[END_REF], cyclo(L-Met-L-Met) [START_REF] Dyer | Accumulation of Maillard reaction products in skin collagen in diabetes and aging[END_REF] and cyclo(L-Pro-L-Pro) [START_REF] André | A method to quantify intracellular glycation in dermal fibroblasts using liquid chromatography coupled to fluorescence detection -Application to the selection of deglycation compounds of dermatological interest[END_REF] showed the most important activity by significantly decreasing the intracellular amount of CML by 27.0-28.0 %.

Cell viability was assessed during cell counting using the trypan blue exclusion principle and revealed no decrease in fibroblast viability following treatment with diketopiperazines (viability >95%).

CML is a chemically very stable molecule. It is unlikely that the biological activity of DKP described in this paper results from a direct chemical action on CML. At the same time, it has been described that DKPs, due to their small steric size, are capable of binding to many cellular receptors [START_REF] Martins | Diketopiperazines: biological activity and synthesis[END_REF]. One of the hypotheses, which remains to be demonstrated, to explain the activity of diketopiperazines would be related to the AGER receptor. The AGER receptor has been described in many cells including fibroblasts [START_REF] Ott | Role of advanced glycation end products in cellular signaling[END_REF]. Its subunit 1 (AGER-1, OST 48), present on both the cell surface and the surface of the endoplasmic reticulum, specifically fixes terminal glycation products, without binding early glycation products such as Amadori products [START_REF] Lu | Advanced glycation endproduct (AGE) receptor 1 is a negative regulator of the inflammatory response to AGE in mesangial cells[END_REF][START_REF] Li | Molecular identity and cellular distribution of advanced glycation endproduct receptors: relationship of p60 to OST-48 and p90 to 80K-H membrane proteins[END_REF][START_REF] Cai | Advanced glycation end product (AGE) receptor 1 suppresses cell oxidant stress and activation signaling via EGF receptor[END_REF]. The binding of intra-and extracellular AGEs to this receptor triggers their degradation, and has a protective effect against oxidative stress [START_REF] Ott | Role of advanced glycation end products in cellular signaling[END_REF][START_REF] Lu | Advanced glycation endproduct (AGE) receptor 1 is a negative regulator of the inflammatory response to AGE in mesangial cells[END_REF]. With age or in the case of diabetes, the level of AGEs in the body is high. In both cases, there is also a decrease in the expression of the AGER-1 receptor [START_REF] Ott | Role of advanced glycation end products in cellular signaling[END_REF][START_REF] Lu | Advanced glycation endproduct (AGE) receptor 1 is a negative regulator of the inflammatory response to AGE in mesangial cells[END_REF]. This hypothesis that the activity of diketopiperazines on the elimination of AGEs in dermal fibroblasts is supported by an increase of a cellular detoxification process will have to be demonstrated in order to provide more details on the mechanism of action of diketopiperazines.

CONCLUSION

In this article we describe for the first time that some molecules belonging to the 2,5diketopiperazine family act as detoxifying agents, capable of eliminating CML, a final AGE accumulated with age in dermal fibroblasts. 2,5-DKPs are natural cyclic dipeptides, produced by 90% of Gram-negative bacteria, but also by Gram-positive bacteria, fungi, plants and mammals. These molecules are also found in the human body (in the brain and human fluids) [START_REF] Prasad | Food-Derived Neuroactive Cyclic Dipeptides[END_REF][START_REF] Prasad | Cyclo(His-Pro): its distribution, origin and function in the human[END_REF], and in industrially processed food (beer, bread, meat, chocolate…) [START_REF] Stark | Structures, sensory activity, and dose/response functions of 2,5-diketopiperazines in roasted cocoa nibs (Theobroma cacao)[END_REF][START_REF] Gautschi | Chemical Characterization of Diketopiperazines in Beer[END_REF].

We demonstrated that a de-replication strategy combining molecular network Although the mode of action of the diketopiperazines remains to be confirmed, and given the elimination activity of advanced glycation end-products and the relative ease of synthesis of symmetrical compounds, diketopiperazines can be considered as compounds of interest for dermatological and cosmetic applications with an anti-ageing aim.

( 2 .

 2 3 mg, RT = 14.8 min), 11 (4.2 mg, RT = 9.0 min) and sub-fractions F2-4, F2-5, F2-6, F2-7, F2-8, F2-9, F2-11, and F2-13 that were subjected to further purification. Subfraction F2-4 (4.3 mg) was purified by analytical HPLC (Luna PFP (2), mobile phase H20 + 0.1% FA / CH3OH + 0.1% FA, isocratic elution 85:15 over 40 minutes, flow rate 1 mL/min) to afford compounds 9 (0.4 mg, RT = 31.7 min), and 10 (1.2 mg, RT = 34.0 min). Sub-fraction F2-5 (1.5 mg) was purified by analytical HPLC (Luna PFP (2), mobile phase H20 + 0.1% FA / CH3OH + 0.1% FA, isocratic elution 77:23 over 20 minutes, flow rate 1 mL/min) to afford compound 4 (0.4 mg, RT = 17.8 min). Sub-fraction F2-6 (10 mg) was purified by analytical HPLC (Luna PFP (2), mobile phase H20 + 0.1% FA / CH3OH + 0.1% FA, isocratic elution 80:20 over 30 minutes, flow rate 1 mL/min) to afford compounds 17 (1.5 mg, RT = 15.5 min), and 5 (0.4 mg, RT = 25.1 min). Sub-fraction F2-7 (1.5 mg) was purified by analytical HPLC (Luna PFP (2), mobile phase H20 + 0.1% FA / CH3OH + 0.1% FA, isocratic elution 80:20 over 35 minutes, flow rate 1 mL/min) to yield compound 13 (0.2 mg, RT = 17.6 min). Sub-fraction F2-8 (6 mg) was purified by analytical HPLC (Luna PFP (2), mobile phase H20 + 0.1% FA / CH3OH + 0.1% FA, isocratic elution 80:20 over 25 minutes, flow rate 1 mL/min) to afford compounds 14 (2.1 mg, RT = 13.7 min), 7 (1.3 mg, RT = 17.0 min), and 2 (0.7 mg, RT = 19.3 min). Subfraction F2-9 (3.1 mg) was purified by analytical HPLC (Luna PFP (2), mobile phase H20 + 0.1% FA / CH3OH + 0.1% FA, isocratic elution 80:20 over 35 minutes, flow rate 1 mL/min) to yield compound 12 (1.4 mg, RT = 26.9 min). Sub-fraction F2-11 (2.4 mg) was purified by analytical HPLC (Luna PFP (2), mobile phase H20 + 0.1% FA / CH3OH + 0.1% FA, isocratic elution 80:20 over 30 minutes, flow rate 1 mL/min) to afford compounds 3 (0.4 mg, RT = 21.7 min), and 1 (0.5 mg, RT = 23.0 min). Finally, subfraction F2-13 (4.9 mg) was purified by analytical HPLC (Luna PFP (2), mobile phase H20 + 0.1% FA / CH3OH + 0.1% FA, isocratic elution 60:40 over 12 minutes, flow rate 1 mL/min) to afford compound 16 (0.7 mg, RT = 8.0 min).

Figure 1 .

 1 Figure 1. Molecular network of the crude extract of Sphingobacterium sp.. Clusters A and B contain DKP structures (yellow nodes). C. Example of MS 2 spectrum (m/z 197.1284; RT 4.55 min) showing ion fragments characteristic of amino acids and losses characteristic of the DKP skeleton.

Figure 2 .

 2 Figure 2. 2,5-Diketopiperazines 1-12 isolated from the crude extract of SNB-CN13 by bioguided fractionation.

Figure 3 .

 3 Figure 3. Synthesis process of symmetrical 2,5-DKP.

Figure 4 .

 4 Figure 4. Symmetrical 2,5-DKPs synthetized.

Figure 5 .

 5 Figure 5. Confocal microscopy images (x 40) of human dermal fibroblasts treated for 24 h with (a) DMSO as control, (b) compound 9, (c) compound 19, (d) compound 20, and (e) compound 21 at 1 µg/mL. The nucleus was labelled by DAPI (blue); The immunological labeling of CML was performed with Alexa 488 (green).

  the images recorded under the confocal photonic microscope was quantified. The expression of CML in fibroblasts treated with molecules 1-12 and 19-21 was compared to the control. The results were expressed as a percentage of the decrease in CML compared to the control group for each treatment and are presented in Table

a

  Significance levels when compared to negative control were calculated by one-way Anova followed by Dunett's multiple comparison test. *** p < 0.001, * p < 0.05, ns = not significant Among molecules 1-12 and 19-21, twelve of them significantly decreased the intracellular amount of CML in fibroblasts, with values ranging from 8.9 % for cyclo(L-

  technique and MS2 data analysis was an effective tool to detect the presence of 26 potential DKPs in the crude extract of Sphingobacterium sp. SNB-CN13 and to predict their amino acid composition. The chemical study led to the isolation of 12 of the 26 diketopiperazines identified by the molecular network technique. Evaluation of the biological activity of the isolated DKPs showed that these molecules can remove advanced glycation end products present in the skin with age. It is interesting to note that symmetrical DKPs synthesized following the principles of green chemistry are also active on the same target.

Table I .

 I Percentage of decrease in CML in fibroblasts treated with molecules 1-21.

	Molecules	Decrease of CML content (% compared to control)	p-value a