Mohamed Amine Bouhlel
email: mohamed.bouhlel@onera.fr

Nathalie Bartoli
email: nathalie.bartoli@onera.fr

Abdelkader Otsmane
email: abdelkader.otsmane@snecma.fr

Joseph Morlier
email: joseph.morlier@isae.fr

Improving kriging surrogates of high-dimensional design models by Partial Least Squares dimension reduction

Keywords: Kriging, Partial Least Squares, Experiment design, Metamodels

Engineering computer codes are often computationally expensive. To lighten this load, we exploit new covariance kernels to replace computationally expensive codes with surrogate models. For input spaces with large dimensions, using the kriging model in the standard way is computationally expensive because a large covariance matrix must be inverted several times to estimate the parameters of the model. We address this issue herein by constructing a covariance kernel that depends on only a few parameters. The new kernel is constructed based on information obtained from the Partial Least Squares method. Promising results are obtained for numerical examples with up to 100 dimensions, and significant computational gain is obtained while maintaining sufficient accuracy.

Number of principal components retained x 1 × d vector x j j th element of a vector x X n × d matrix containing sampling points y n × 1 vector containing simulation of X x (i) i th training point for i = 1, . . . , n (a 1 × d vector) w (l) d × 1 vector containing X weights given by the l th PLS iteration for l = 1, . . . , h X (0) X X (l-1) Matrix containing residual of inner regression of (l -1) st PLS iteration for l = 1, .

. . , h k(•, •)

Covariance function N (0, k(•, •)) Distribution of a Gaussian process with mean function 0 and covariance function k(•, •) x t Superscript t denotes the transpose operation of the vector x

Introduction and main contribution

In recent decades, because simulation models have striven to more accurately represent the true physics of phenomena, computational tools in engineering have become ever more complex and computationally expensive. To address this new challenge, a large number of input design variables, such as geometric representation, are often considered. Thus, to analyze the sensitivity of input design variables or to search for the best point of a physical objective under certain physical constraints (i.e., global optimization), a large number of computing iterations are required, which is impractical when using simulations in real time. This is the main reason that surrogate modeling techniques have been growing in popularity in recent years. Surrogate models, also called metamodels, are vital in this context and are widely used as substitutes for time-consuming highfidelity models. They are mathematical tools that approximate coded simulations of a few well-chosen experiments that serve as models for the design of experiments. The main role of surrogate models is to describe the underlying physics of the phenomena in question. Different types of surrogate models can be found in the literature, such as regression, smoothing spline [START_REF] Wahba | Smoothing noisy data with spline functions. Estimating the correct degree of smoothing by the method of generalized cross-validation[END_REF][START_REF] Wahba | Spline models for observational data[END_REF]), neural networks [START_REF] Haykin | Neural networks: a comprehensive foundation[END_REF], radial basis functions [START_REF] Buhmann | Radial basis functions: theory and implementations[END_REF] and Gaussian-process modeling [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF].

In this article, we focus on the kriging model because it estimates the prediction error. This model is also referred to as the Gaussian-process model [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] and was presented first in geostatistics (see, e.g., [START_REF] Cressie | Spatial prediction and ordinary kriging[END_REF] or [START_REF] Goovaerts | Geostatistics for natural resources evaluation (applied geostatistics)[END_REF]) before being extended to computer experiments and machine learning [START_REF] Schonlau | Computer experiments and global optimization[END_REF][START_REF] Sasena | Flexibility and efficiency enhancements for constrained global design optimization with Kriging approximations[END_REF][START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF][START_REF] Picheny | Adaptive designs of experiments for accurate approximation of a target region[END_REF]). The kriging model has become increasingly popular due to its flexibility in accurately imitating the dynamics of computationally expensive simulations and its ability to estimate the error of the predictor. However, it suffers from some well-known drawbacks in high dimension, which may be due to multiple causes. For starters, the size of the covariance matrix of the kriging model may increase dramatically if the model requires a large number of sample points. As a result, inverting the covariance matrix is computationally expensive. The second drawback is the optimization of the subproblem, which involves estimating the hyperparameters for the covariance matrix. This is a complex problem that requires inverting the covariance matrix several times. Some recent works have addressed the drawbacks of high-dimensional Gaussian processes [START_REF] Hensman | Gaussian processes for big data[END_REF][START_REF] Damianou | Deep gaussian processes[END_REF][START_REF] Durrande | Covariance kernels for simplified and interpretable modeling. A functional and probabilistic approach[END_REF] or the large-scale sampling of data [START_REF] Sakata | An efficient algorithm for Kriging approximation and optimization with large-scale sampling data[END_REF]. One way to reduce CPU time when constructing a kriging model is to reduce the number of hyper-parameters, but this approach assumes that the kriging model exhibits the same dynamics in all directions [START_REF] Mera | Efficient optimization processes using kriging approximation models in electrical impedance tomography[END_REF].

Thus, because estimating the kriging parameters can be time consuming, especially with dimensions as large as 100, we present herein a new method that combines the kriging model with the Partial Least Squares (PLS) technique to obtain a fast predictor. Like the method of principle components analysis (PCA), the PLS technique reduces dimension and reveals how inputs depend on outputs. PLS is used in this work because PCA only exposes dependencies between inputs. Information given by PLS is integrated in the covariance structure of the kriging model to reduce the number of hyper-parameters. The combination of kriging and PLS is abbreviated KPLS and allows us to build a fast kriging model because it requires fewer hyper-parameters in its covariance function; all without eliminating any input variables from the original problem. In general, the number of kriging parameters is equal to the number of dimension which is reduced to at maximum 4 parameters with our approach.

The KPLS methods is used for many academic and industrial verifications, and promising results have been obtained for problems with up to 100 dimensions. The cases used in this paper do not exceed 100 input variables, which should be quite sufficient for most engineering problems. Problems with more than 100 inputs may lead to memory difficulties with the toolbox Scikitlearn (version 0.14), on which the KPLS method is based.

This paper is organized as follows: Section 2 summarizes the theoretical basis of the universal kriging model, recalling the key equations. The proposed KPLS model is then described in detail in Section 3 by using the kriging equations. Section 4 compares and analyzes the results of the KPLS model with those of the kriging model when applied to classic analytical examples and some complex engineering examples. Finally, Section 5 concludes and gives some perspectives.

Universal kriging model

To understand the mathematics of the proposed methods, we first review the kriging equations. The objective is to introduce the notation and to briefly describe the theory behind the kriging model. Assume that we have evaluated a cost deterministic function at n points x (i)

(i = 1, . . . , n) with x (i) = x (i) 1 , . . . , x (i) d ∈ B ⊂ R d ,
and we denote by X the matrix [x (1)t , . . . , x (n)t ] t . For simplicity, B is considered to be a hypercube expressed by the product between intervals of each direction space, i.e., B = d j =1 [a j , b j ], where a j , b j ∈ R with a j ≤ b j for j = 1, . . . , d. Simulating these n inputs gives the outputs y = [y (1) , . . . , y (n) ] t with y (i) = y(x (i) ) for i = 1, . . . , n. We use ŷ(x) to denote the prediction of the true function y(x) which is considered as a realization of a stochastic process Y (x) for all x ∈ B. For the universal kriging model [START_REF] Roustant | DiceKriging, DiceOptim: two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF][START_REF] Picheny | Adaptive designs of experiments for accurate approximation of a target region[END_REF], Y is written as

Y (x) = m j =1 β j f j (x) + Z(x), (1)
where, for j = 1, . . . , m, f j is a known independent basis function, β j ∈ R is an unknown parameter, and Z is a random variable defined by Z(x) ∼ N (0, k), with k being a stationary covariance function, also called a covariance kernel. The kernel function k can be written as

k(x, x ) = σ 2 r(x, x ) = σ 2 r xx ∀ x, x ∈ B, ( 2
)
where σ 2 is the process variance and r xx is the correlation function between x and x . However, the correlation function r depends on some hyper-parameters θ and it is considered to be known. We also denote the n × 1 vector as

r xX = [r xx (1) , . . . , r xx (n) ] t and the n × n covariance matrix as R = [r x (1) X , . . . , r x (n) X ].

Derivation of prediction formula

Under the hypothesis above, the best linear unbiased predictor for y(x), given the observations y, is

ŷ(x) = f(x) t β + r t xX R -1 y -F β , ( 3
)
where f(x) = [f 1 (x), . . . , f m (x)] t is the m × 1 vector of basis functions, F = f(x (1) ), . . . , f(x (n) ) t is the n × m matrix, and β is the vector of generalized least-square estimates of β = [β 1 , . . . , β m ] t , which is given by

β = ⎡ ⎢ ⎣ β1 . . . βm ⎤ ⎥ ⎦ = F t R -1 F -1 F t R -1 y. (4)
Moreover, the universal kriging model provides an estimate of the variance of the prediction, which is given by

s 2 (x) = σ 2 1 -r t xX R -1 r xX , (5)
with

σ 2 = 1 n y -F β t R -1 y -F β . ( 6
)
For more details of the derivation of the prediction formula, see, for instance, [START_REF] Sasena | Flexibility and efficiency enhancements for constrained global design optimization with Kriging approximations[END_REF] or [START_REF] Schonlau | Computer experiments and global optimization[END_REF]. The theory of the proposed method has been expressed in the same way as for the universal kriging model. The numerical examples in Section 4 use the ordinary kriging model, which is a special case of the universal model, but with f(x) = {1} (and m = 1). For the ordinary kriging model, (3), (4), and (6) are then replaced by the equations given in Appendix A.

Note that the assumption of known covariance with known hyper-parameters θ is unrealistic in reality. For this reason, the covariance function is typically chosen from among a parametric family of kernels. Table 12 in Appendix B gives some examples of typical stationary kernels. The number of hyper-parameters required for the estimate is typically greater than (or equal to) the number of input variables. In this work, we use in the following a Gaussian exponential kernel:

k(x, x ) = σ 2 d i=1 exp -θ i x i -x i 2 ∀ θ i ∈ R + .
By applying some elementary operations to existing kernels, we can construct new kernels. In this work, we use the property that the tensor product of covariances is a covariance kernel in the product space. More details are available in [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], [START_REF] Durrande | Covariance kernels for simplified and interpretable modeling. A functional and probabilistic approach[END_REF], [START_REF] Bishop | Pattern recognition and machine learning (information science and statistics[END_REF], [START_REF] Liem | Surrogate models and mixtures of experts in aerodynamic performance prediction for mission analysis[END_REF].

Estimation of hyper-parameters θ

The key point of the kriging approximation is how it estimates the hyper-parameters θ, so its main steps are recalled here, along with some mathematical details.

One of the major challenges when building a kriging model is the complexity and difficulty of estimating the hyper-parameters θ, in particular when dealing with problems with many dimensions or with a large number of sampling points. In fact, using (3) to make a kriging prediction requires inverting an n × n matrix, which typically has cost of O n 3 , where n is the number of sampling points [START_REF] Braham | 12th international symposium on modeling and optimization in mobile, ad[END_REF]. The hyper-parameters are estimated using maximum likelihood (ML) or cross validation (CV), which are based on observations. Bachoc compared the ML and CV techniques [START_REF] Bachoc | Cross Validation and Maximum Likelihood estimation of hyper-parameters of Gaussian processes with model misspecification[END_REF]) and concluded that, in most cases studied, the CV variance is larger. The ML method is widely used to estimate the hyper-parameters θ; it is also used in this paper. In practice, the following log-ML estimate is often used:

log -ML(θ) = - 1 2 n ln(2πσ 2 ) + ln(det R(θ)) +(y -Fβ) t R(θ) -1 (y -Fβ)/σ 2 . ( 7
)
Inserting β and σ 2 given by ( 4) and ( 6), respectively, into the expression (7), we get the following so-called concentrated likelihood function, which depends only on the hyper-parameters θ:

log-ML(θ) = - 1 2 [n ln σ 2 + ln det R(θ)] = - 1 2 n ln 1 n (y -F(F t R -1 F) -1 F t R -1 y) t ×R -1 (y -F(F t R -1 F) -1 F t R -1 y) + ln det R . ( 8
)
To facilitate reading, R(θ) has been replaced by R in the last line of (8).

Maximizing ( 8) is very computationally expensive for high dimensions and when using a large number of sample points because the (n × n) matrix R(θ) in ( 8) must be inverted. The maximization problem is often solved using genetic algorithms (see Forrester and Sobester 2008 for more details). In this work, we use the derivative-free optimization algorithm COBYLA that was developed by [START_REF] Powell | A direct search optimization method that models the objective and constraint functions by linear interpolation[END_REF]. COBYLA is a sequential trust-region algorithm that uses linear approximations for the objective and constraint functions.

Figure 1 recalls the principal stages of building a kriging model, and each step is briefly outlined below:

1. The user must provide the initial design of experiments (X, y) and the parametric family of the covariance function k. 2. To derive the prediction formula, the kriging algorithm assumes that all parameters of k are known.

Fig. 1 The main steps for building an ordinary kriging model 3. Under the hypothesis of the kriging algorithm, we estimate hyper-parameters θ from the concentrated likelihood function given by ( 8) and by using the COBYLA algorithm. 4. Finally, we calculate the prediction (3) and the associated estimation error ( 5) after estimating all hyperparameters of the kriging model.

Kriging model combined with Partial Least Squares

As explained above, maximizing the concentrated likelihood ( 8) can be time consuming when the number of covariance parameters is large, which typically occurs in large dimension. Solving this problem can be accelerated by combing the PLS method and the kriging model. The θ parameters from the kriging model represent the range in any spatial direction. Assuming, for instance, that certain values are less significant for the response, then the corresponding θ i (i = 1, . . . , d) will be very small compared to the other θ parameters. The PLS method is a well-known tool for high-dimensional problems and consists of maximizing the variance by projecting onto smaller dimensions while monitoring the correlation between input variables and the output variable. In this way, the PLS method reveals the contribution of all variablesthe idea being to use this information to scale the θ parameters.

In this section we propose a new method that can be used to build an efficient kriging model by using the information extracted from the PLS stage. The main steps for this construction are as follows:

1. Use PLS to define weight parameters. 2. To reduce the number of hyper-parameters, define a new covariance kernel by using the PLS weights. 3. Optimize the parameters.

The key mathematical details of this construction are explained in the following.

Linear transformation of covariance kernels

Let x be a vector space over the hypercube B. We define a linear map given by

F : B -→ B , x -→ [α 1 x 1 , . . . , α d x d ] , (9)
where

α 1 , . . . , α d ∈ R and B is a hypercube included in R d (B can be different from B). Let k be an isotropic covariance kernel with k : B × B → R. Since k is isotropic, the covariance kernel k(F (•), F (•))
depends on a single parameter, which must be estimated. However, if α 1 , . . . , α d are well chosen, then we can use k(F (•), F (•)). In this case, the linear transformation F allows us to approach the isotropic case [START_REF] Zimmerman | A network design criterion for estimating selected attributes of the semivariogram[END_REF]. In the present work, we choose α 1 , . . . , α d based on information extracted from the PLS technique.

Partial Least Squares

The PLS method is a statistical method that finds a linear relationship between input variables and the output variable by projecting input variables onto a new space formed by newly chosen variables, called principal components (or latent variables), which are linear combinations of the input variables. This approach is particularly useful when the original data are characterized by a large number of highly collinear variables measured on a small number of samples. Below, we briefly describe how the method works. For now, suffice it to say that only the weighting coefficients are central to understanding the new KPLS approach. For more details on the PLS method, please see [START_REF] Helland | On structure of partial least squares regression[END_REF], [START_REF] Frank | A statistical view of some chemometrics regression tools[END_REF], [START_REF] Alberto | Partial Least Squares regression on symmetric positive-definite matrices[END_REF].

The PLS method is designed to search out the best multidimensional direction in X space that explains the characteristics of the output y. After centering and scaling the (n × d)-sample matrix X and the response vector y, the first principal component t (1) is computed by seeking the best direction w (1) that maximizes the squared covariance between t (1) = Xw (1) and y:

w (1) = arg max w w t X t yy t Xw such that w t w = 1. ( 10
)
The optimization problem ( 10) is maximized when w (1) is the eigenvector of the matrix X t yy t X corresponding to the eigenvalue with the largest absolute value; the vector w (1) contains the X weights of the first component. The largest eigenvalue of problem ( 10) can be estimated by the power iteration method introduced by [START_REF] Lanczos | An iteration method for the solution of the eigenvalue problem of linear differential and integral operators[END_REF].

Next, the residual matrix from X = X (0) space and from y = y (0) are calculated; these are denoted X (1) and y (1) , respectively:

X (1) = X (0) -t (1) p (1) , y (1) = y (0) -c 1 t (1) , ( 11
)
where p (1) (a 1 × d vector) contains the regression coefficients of the local regression of X onto the first principal component t (1) , and c 1 is the regression coefficient of the local regression of y onto the first principal component t (1) .

Fig. 2 Upper left shows construction of two principal components in X space. Upper right shows prediction of y (0) . Bottom left shows prediction of y (1) . Bottom right shows final prediction of y

The system (11) is the local regression of X and y onto the first principal component.

Next, the second principal component, which is orthogonal to the first, can be sequentially computed by replacing X by X (1) and y by y (1) to solve the system (10). The same approach is used to iteratively compute the other principal components. To illustrate this process, a simple three-dimensional (3D) example with two principal components is given in Fig. 2. In the following, we use h to denote the number of principal components retained.

The principal components represent the new coordinate system obtained upon rotating the original system with axes, x 1 , . . . , x d [START_REF] Alberto | Partial Least Squares regression on symmetric positive-definite matrices[END_REF]. For l = 1, . . . , h, t (l) can be written as

t (l) = X (l-1) w (l) = Xw (l) * . ( 12
)
This important relationship is used for coding the method.

The following matrix W * = [w (1) * , . . . , .w (h) * ] is obtained by [START_REF] Manne | Analysis of two Partial-Least-Squares algorithms for multivariate calibration[END_REF]:

W * = W P t W -1 ,
where W = [w (1) , . . . , w (h) ] and P = [p (1) t , . . . , p (h) t ].

The vector w (l) corresponds to the principal direction in X space that maximizes the covariance of X (l-1) t y (l-1) y (l-1) t X (l-1) . If h = d, the matrix W * = [w (1) * , . . . , w (d) * ] rotates the coordinate space (x 1 , . . . , x d ) to the new coordinate space (t (1) , . . . , t (d) ), which follow the principal directions w (1) , . . . , w (d) .

As mentioned in the introduction, the PLS method is chosen instead of the PCA method because the PLS method shows how the output variable depends on the input variables, whereas the PCA method focuses only on how the input variables depend on each other. In fact, the hyper-parameters θ for the kriging model depend on how each input variable affects the output variable.

Construction of new kernels for KPLS models

Let B be a hypercube included in R d . As seen in the previous section, the vector w

(1) * is used to build the first principal component t (1) = Xw

(1) * , where covariance between t (1) and y is maximized. The scalars w (1) * 1 , . . . , w (1) * d can then be interpreted as measuring the importance of x 1 , . . . , x d , respectively, for constructing the first principal component where its correlation with the output y is maximized. However, we know that the hyper-parameters θ 1 , . . . , θ d (see Table 12 in Appendix B) can be interpreted as measuring how strongly the variables x 1 , . . . , x d , respectively, affect the output y. Thus, we define a new kernel k kpls1 : B ×B → R given by k 1 (F 1 (•), F 1 (•)) with k 1 : B × B → R being an isotropic stationary kernel and

F 1 : B -→ B, x -→ w (1) * 1 x 1 , . . . , w (1) * d x d . ( 13
)
F 1 goes from B to B because it only works for the new coordinate system obtained by rotating the original coordinate axes, x 1 , . . . , x d . Through the first component t (1) , the elements of the vector w (1) * reflect how x depends on y. However, such information is generally insufficient, so the elements of the vector w

(1) * are supplemented by the information given by the other principal components t (2) , . . . , t (h) . Thus, we build a new kernel k kpls1:h sequentially by using the tensor product of all kernels k kplsl , which accounts for all this information in only a single covariance kernel:

k kpls1:h (x, x ) = h l=1 k l (F l (x) , F l (x )), ( 14
)
with k l : B × B → R and

F l : B -→ B x -→ w (l) * 1 x 1 , . . . , w (l) * d x d .
(

) 15
If we consider the Gaussian kernel applied with this proposed approach, we get (i) We build KPLS models based on h = 1, then h = 2, . . . principal components. (ii) We choose the number of components that minimizes the leave-one-out cross-validation error.

k(x, x ) = σ 2 h l=1 d i=1 exp -θ l w (l) * i x i -w (l) * i x i 2 , ∀ θ l ∈ [0, +∞[.
The flowchart given in Fig. 3 shows how the information flows through the algorithm, from sample data, PLS algorithm, kriging hyper-parameters, to final predictor. With the same definitions and equations, almost all the steps for constructing the KPLS model are similar to the original steps for constructing the ordinary kriging model. The exception is the third step, which is highlighted in the solid-red box in Fig. 3. This step uses the PLS algorithm to define the new parametric kernel k kpls1:h as follows:

a. initialize the PLS algorithm with l = 1; b. if l = 1, compute the residual of X (l-1) and y (l-1) by using system (11); ) is what differentiates this approach from the ordinary kriging approach c. compute X weights for iteration l; d. define a new kernel k kpls1:h by using ( 14); e. if the number of iterations is reached, return to step 3, otherwise continue; f. update data considering l = l + 1.

Note that, if kernels k l are separable at this point, the new kernel given by ( 14) is also separable. In particular, if all kernels k l are of the exponential type (e.g., all Gaussian exponentials), the new kernel given by ( 14) will be the same type as k l . The proof is given in Appendix C.

Numerical examples

We now present a few analytical and engineering examples to verify the proper functioning of the proposed method. The ordinary kriging model with a Gaussian kernel provides the benchmark against which the results of the proposed combined approach are compared. The Python toolbox Scikit-learn v.014 [START_REF] Pedregosa | Scikit-learn: machine learning in python[END_REF]) is used to implement these numerical tests. This toolbox provides hyper-parameters for the ordinary kriging. The computations were done on an Intel ® Celeron ® CPU 900 2.20 GHz desktop PC. For the proposed method, we combined an ordinary kriging model with a Gaussian kernel with the PLS method with one to three principal components.

Analytical examples

We use two academic functions and vary the characteristics of these test problems to cover most of the difficulties faced in the field of substitution models. The first function is g07 [START_REF] Michalewicz | Evolutionary algorithms for constrained parameter optimization problems[END_REF] with 10 dimensions,

y g07 (x) = x 2 1 + x 2 2 + x 1 x 2 -14x 1 -16x 2 +(x 3 -10) 2 + 4(x 4 -5) 2 + (x 5 -3) 2 +2(x 6 -1) 2 + 5x 2 7 + 7(x 8 -11) 2 +2(x 9 -10) 2 + (x 10 -7) 2 + 45, -10 ≤ x i ≤ 10, fori = 1, . . . , 10.
For this function, we use experiments based on a latin hypercube design with 100 data points to fit models.

The second function is the Griewank function (Regis and Shoemaker 2013), which is used because of its complexity, as illustrated in Fig. 4 for the two-dimensional (2D) case. The function is

y Griewank (x) = d i=1 x 2 i 4000 - d i=1 cos x i √ i + 1, -600 ≤ x i ≤ 600, for i = 1, . . . , d.
Two types of experiments are done with this function. The first is defined over the interval [-600, 600] and has varying dimensions (2,5,7,10,20,60). This experiment serves to verify the effectiveness of the proposed approach in both low and high dimensions. It is based on the latin hypercube design and uses n data points to fit models, as mentioned in Table 1. The second type of experiment is defined over the interval [-5, 5], where the Griewank function is more complex than for the first type of experiment (cf. Figs. 4 and5). Over this reduced interval, experiments are done with 20 and 40 dimensions (20D and 40D) and with 50, 100, 200, and 300 sampling points. To analyse the robustness of the method, ten experiments, each with a different latin hypercube design, are used for this case.

To compare the three approaches (i.e., the g07 function and the Griewank function of the intervals [-600, 600] and [-5, 5]), 5000 random points are computed and the results are stored in a database. The following relative error is used to compare the performance of the ordinary kriging model with the KPLS model:

Error = || Ŷ -Y|| 2 ||Y|| 2 100, (16)
where || • || 2 represents the usual L 2 norm, and Ŷ and Y are the vectors containing the prediction and the real values of random points, respectively. The CPU time required to fit models is also noted ("h" refers to hours, "min" refers to minutes, and "s" refers to seconds).

Comparison with g07 function

The results listed in Table 2 show that the proposed KPLS surrogate model is more accurate than the ordinary kriging model when more than one component is used.

Using just one component gives almost the same accuracy as the ordinary kriging model. In this case, only a single θ hyper-parameter from the space correlation needs be estimated compared to ten θ hyper-parameters for the [-600, 600] Table 3 compares the ordinary kriging model and the KPLS model in two dimensions.

Comparison with complex Griewank function over interval

If two components are used for the KPLS, we expect to obtain the same accuracy and time cost for the two approaches because the difference between the two models consists only of a transformation of the search-space coordinates when a Gaussian kernel is used (the space in which the θ hyper-parameters exist). In this case, the KPLS model with only one component degrades the accuracy of the solution.

Tables 4,5, and 6, show the results for 5, 7, and 10 dimensions, respectively.

Varying the number of principal components does not significantly affect the accuracy of the model. The gain in computation time does not appear upon increasing the number of principal components: the computation time is reduced when we use the KPLS model. Upon increasing the number of principal components, the CPU time for constructing the KPLS model increases but In the 20-dimension (20D) example (Table 7), using KPLS with only one principal component leads to a poor relative error (10.15 %) compared with other models. In this case, two principal components are required to build the combined model. The CPU time remains less than that for the ordinary kriging model (11.7 s vs 107 s).

The results in Table 8 for the KPLS model with 60 dimensions (60D) show that this model is faster than the ordinary kriging model. Compared with the kriging model, the CPU time is reduced 42-fold when one principal component is used and over 17-fold when three principal components are used.

Thus, for the Griewank function over the interval [-600, 600] and at the highest dimensions, the majority of the results obtained for the analytical examples are better when using the KPLS model than when using the ordinary kriging model. The proposed method thus appears interesting, particularly in terms of saving CPU time while maintaining good accuracy. [-5, 5] As shown in Fig. 4, the Griewank function looks like a parabolic function. This is because, over the interval [-600, 600], the cosine part of the Griewank function does not contribute significantly compared with the sum of x 2 i /4000. This is especially true given that the cosine part is a product of factors each of which is less than unity. If we reduce the interval from [-600, 600] to [-5, 5], we can see why the Griewank function is widely used as a multimodal test function with a very rugged landscape and a large number of local optima (see Fig. 5). Compared with the interval [-600, 600], the opposite happens for the interval [-5, 5]: the "cosine part" dominates; at least for moderate dimensions where the product contains few factors. For this case, which seems very difficult, we consider 20 and 60 input variables. For each problem, ten experiments based on the latin hypercube design are built with 50, 100, 200, and 300 sampling points. The mean and the standard error are given in Tables 14 and15 For 20 input variables and 50 sampling points, the KPLS model gives a more accurate solution than the ordinary kriging model, as shown in Fig. 7a. The rate of improvement with respect to the number of sampling points is less for the KPLS model than for the kriging model (cf. Fig. 7b-d). Nevertheless, the results shown in Fig. 8 indicate that the KPLS model leads to an important reduction in CPU time for the various number of sampling points.

Comparison with complex Griewank function over interval

Similar results occur for the 60D Griewank function (Fig. 9). The mean RE obtained with the ordinary kriging model improves from approximately 1.39 to 0.65 % upon increasing the number of sampling points from 50 to 300 (cf. Fig. 9a,d). However, a very important reduction in CPU time is obtained, as shown in Fig. 10. The CPU time required for the KPLS model is hardly visible because it is much, much less than that required by the ordinary kriging model. We thus show in Fig. 11 the CPU time required by the KPLS model alone to show the different CPU times required for the various configurations (KPLS1, KPLS2, and KPLS3). For Griewank function over the interval [-5, 5], the KPLS method seems to perform well when the number of observations is small compared to the dimension d. In this case, the standard separable covariance function for the ordinary kriging model is almost impossible to use because the number of parameters to be estimated is too large compared with the number of observations. Thus, the KPLS method seems more efficient in this case.

Industrial examples

The following engineering examples are based on results of numerical experiments done at SNECMA on multidisciplinary optimization. The results are stored in tables.

Aerospace turbomachinery consists of numerous blades that transfer energy between air and the rotor. The disks with compressor blades are particularly important because they must satisfy the dual criteria of aerodynamic performance and mechanical stress. Blades are mechanically and aerodynamically optimized by searching parameter space for an aerodynamic shape that ensures the best compromise that satisfies a set of constraints. The blade, which is a 3D entity, is first divided into a number of radial 2D profiles whose thickness is a given percentage of the distance from the hub to the shroud (see Fig. 6).

A new 3D blade is constructed by starting with the 2D profiles and then exporting them to various meshing tools before analyzing them in any specific way. The calculation is integrated into the Optimus platform (Noesis Solutions 2009), which makes it possible to integrate multiple engineering software tools (computational structural mechanics, computational fluid dynamics, . . . ) into a single automated work flow. Optimus, which is an industrial software package for screening variables, optimizing design, and analyzing the sensitivity and robustness, explores and analyzes the results of the work-flow to optimize product design. Toward this end, it uses high-fidelity codes or a reduced model of these codes. It also exploits a wide range of approximation models, including the ordinary kriging model.

Input variables designate geometric hyper-parameters at different percent height and outputs are related to aerodynamic efficiency, vibration criteria, mechanical stress, geometric constraints, and aerodynamic stress. Three numerical experiments are considered: 16) and the CPU time required to fit the model are reported in Tables 91011.

The relative errors for the four models are very similar: the KPLS model results in a slightly improved accuracy for the solutions y 1 , y 2 , y 4 from tab 1 , y 1 from tab 2 , and y 1 from tab 3 but degrades slightly the solution y 3 from tab 1 . The main improvement offered by the proposed model relates to the time required to fit the model, particularly for a large number of training points. Table 10 shows that, with only one principal component, the CPU time is drastically reduced compared with the Optimus model. More precisely, for tab 2 , the ordinary kriging model requires 1 h 30 min whereas the KPLS1 model requires only 11 s and provides better accuracy. In addition, the results for KPLS2 and 11). One other point of major interest for the proposed method is its natural compatibility with sequential enrichment techniques such as the efficient global optimization strategy (see [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF].

Dimensional limits

This project is financed by SNECMA and most of their design problems do not exceed 100 input variables. In addition, the toolbox Scikit-learn (version 0.14) may have memory problems when a very large number of input variables is considered. Thus, problems with more than 100 input variables are not investigated in this work. However, by optimizing memory access and storage, this limit could easily be increased.

Conclusion and future work

Engineering problems that require integrating surrogate models into an optimization process are receiving increasing interest within the multidisciplinary optimization community. Computationally expensive design problems can be solved efficiently by using, for example, a kriging model, which is an interesting method for approximating and replacing high-fidelity codes, largely because these models give estimation errors, which is an interesting way to solve optimization problems. The major drawback involves the construction of the kriging model and in particular the large number of hyper-parameters that must be estimated in high dimensions. In this work, we develop a new covariance kernel for handling this type of higher-dimensional problem (up to 100 dimensions). Although the PLS method requires a very short computation time to estimate θ, the estimate is often difficult to execute and computationally expensive when the number of input variables is greater than 10. The proposed KPLS model was tested by applying it to two analytic functions and by comparing its results to those tabulated in three industrial databases. The comparison highlights the efficiency of this model for up to 99 dimensions. The advantage of the KPLS models is not only the reduced CPU time, but also in that it reverts to the kriging model when the number of observations is small relative to the dimensions of the problem. Before using the KPLS model, however, the number of principal components should be tested to ensure a good balance between accuracy and CPU time.

An interesting direction for future work is to study how the design of the experiment (e.g., factorial) affects the KPLS model. Furthermore, other verification functions and other types of kernels can be used. In all cases studied herein, the first results with this proposed method reveal significant gains in terms of computation time while still ensuring good accuracy for design problems with up to 100 dimensions. The implementation of the proposed KPLS method requires minimal modifications of the classic kriging algorithm and offers further interesting advantages that can be exploited by methods of optimization by enrichment. The covariance functions are written as functions of the ith component m i = |x ix i | with θ i ≥ 0 and p i ∈ [0, 2] for i = 1, . . . , d The covariance functions are written as functions of the ith component m

Fig. 3

 3 Fig. 3 Main steps for constructing KPLS model. The solid-red box (step 3 relative to PLS) is what differentiates this approach from the ordinary kriging approach

Fig. 4

 4 Fig.4A two-dimensional Griewank function over the interval[-600, 600]

  to what is required by industry in terms of dimensions,

Fig. 5

 5 Fig.5A two-dimensional Griewank function over the interval[-5, 5]

  in Appendix D. To better visualize the results, boxplots are used to show CPU time and the relative error RE given by Figs. 7, 8, 9, 10 and 11 in Appendix D.

Fig. 6

 6 Fig. 6 Example of 2D cut of blade (c is chord; CG is gravity center; β 1 is angle for BA; β 2 is angle for BF; Ep is maximum thickness)

  (i) The first experiment is denoted tab 1 and contains 24 input variables and 4 output variables. It has 99 training points and 52 validation points. The outputs are denoted y 1 , y 2 , y 3 , and y 4 . (ii) The second experiment is denoted tab 2 and contains 10 input variables and only 1 output variable. It has 1295 training points and 500 validation points. (iii) The third experiment is denoted tab 3 and contains 99 input variables and 1 output variable. It has 341 training points and 23 validation points.

  * i (x ix i )| with θ l ≥ 0 and p l ∈ [0, 2] for l = 1, . . . , h where F l (x) i is the ith coordinate of F l (x). If we insert (20) in (14) we getk kpls1:h (x, x ) F l (x) i , F l x i ),corresponding to an one-dimensional kernel. Hence, k kpls1:h is a separable kernel. In particular, if we consider a generalized exponential kernel with p 1 = • • • = p h = p ∈ [0, 2], we obtainψ i x i , x i = σ

Fig. 7

 7 Fig. 7 RE for Griewank function in 20D over interval [-5, 5]. Experiments are based on the 10 latin hypercube design

Table 13

 13

	in Appendix B presents new KPLS kernels based
	on examples from Table 12 (also in Appendix B) that con-
	tain fewer hyper-parameters because h	d. The number of
	principal components is fixed by the following leave-one-
	out cross-validation method:

Table 1

 1 Number of data points used for latin hypercube design for the Griewank test function

Table 2

 2 Results for g07 function in ten dimensions with 100-point latin hypercube

	Surrogate	Error (%)	CPU time
	Ordinary kriging	0.013	5.14 s
	KPLS (1 component)	0.014	0.11 s
	KPLS (2 components)	0.0015	0.43 s
	KPLS (3 components)	0.0008	0.44 s

Table 3

 3 Griewank function in two dimensions with 70-point latin hypercube over the interval[-600, 600]

	Surrogate	Error (%)	CPU time
	Ordinary kriging	5.50	0.09 s
	KPLS (1 component)	7.23	0.04 s
	KPLS (2 components)	5.50	0.10 s

Table 4

 4 Griewank function in five dimensions with 100-point latin hypercube over the interval[-600, 600]

	Surrogate	Error (%)	CPU time
	Ordinary kriging	0.605	0.55 s
	KPLS (1 component)	0.635	0.12 s
	KPLS (2 components)	0.621	0.31 s
	KPLS (3 components)	0.623	0.51 s
	ordinary kriging model. Increasing the number of com-
	ponents improves the accuracy of the KPLS surrogate
	model. Whereas the PLS method treats only linearly related
	input and output variables, this example shows that the
	KPLS model can treat nonlinear problems. This result is
	not contradictory because (23) shows that the KPLS model
	is equivalent to the kriging model with specific hyper-
	parameters.

Table 5

 5 Griewank function in seven dimensions with 200-point latin hypercube over the interval[-600, 600]

	Surrogate	Error (%)	CPU time
	Ordinary kriging	0.138	3.09 s
	KPLS (1 component)	0.141	0.25 s
	KPLS (2 components)	0.138	0.52 s
	KPLS (3 components)	0.141	0.94 s

Table 6

 6 Griewank function in ten dimensions with 300-point latin hypercube over the interval[-600, 600]

	Surrogate	Error (%)	CPU time
	Ordinary kriging	0.052	21 s
	KPLS (1 component)	0.033	0.6 s
	KPLS (2 components)	0.035	2.41 s
	KPLS (3 components)	0.034	3.58 s
	still remains lower than for ordinary kriging. For these
	three examples, the combined approach with only one
	PLS component offers sufficient accuracy with a CPU
	time reduced 35-fold for 10 dimensions (i.e., 21 s for
	the ordinary kriging model and 0.6 s for the combined
	model).

Table 7

 7 Griewank function in 20 dimensions with 400-point latin hypercube over the interval[-600, 600]

	Surrogate	Error (%)	CPU time
	Ordinary kriging	0.35	107 s
	KPLS (1 component)	10.15	1.16 s
	KPLS (2 components)	0.003	11.7 s
	KPLS (3 components)	0.002	16.23 s

Table 8

 8 Griewank function in 60 dimensions with 800-point latin hypercube over the interval[-600, 600]

	Surrogate	Error (%)	CPU time
	Ordinary kriging	11.47	293 s
	KPLS (1 component)	7.4	6.88 s
	KPLS (2 components)	6.04	12.57 s
	KPLS (3 components)	5.23	16.82 s

Table 10

 10 Results for tab 2 experiment data (10 input variables, 1 output variable y 1 ) obtained by using 1295 training points, 500 validation points, and error given by (16)Points used in tab 1 , tab 2 , and tab 3 come from previous computationally expensive computer experiments done at SNECMA, which means that this separation between training points and verification points was imposed by SNECMA. The goal is to compare the ordinary kriging model that is implemented in the Optimus platform with the proposed KPLS model. The relative error given by (

	10D	Surrogate	Error (%)	CPU time
	tab 2	Kriging	5.37	1 h 30 min
		KPLS1	5.07	11.69 s
		KPLS2	5.02	1 min 22 s
		KPLS3	5.34	7 min 34 s
	"Kriging" refers to the ordinary kriging optimus solution and
	"KPLSh" refers to the KPLS model with h principal components

Table 9

 9 Results for tab 1 experiment data (24 input variables, 4 output variables y 1 , y 2 , y 3 , y 4 ) obtained by using 99 training points, 52 validation points, and the error given by (16)

	24D	Surrogate	y 1		y 2		y 3		y 4
			Error (%)	CPU	Error (%)	CPU	Error (%)	CPU	Error (%)	CPU
				time		time		time		time
	tab 1	Kriging	0.082	8 s	4.45	8.4 s	8.97	8.17 s	6.27	8.12 s
		KPLS1	0.079	0.12 s	4.04	0.11 s	10.35	0.18 s	5.67	0.11 s
		KPLS2	0.079	0.43 s	4.06	0.69 s	10.33	0.42 s	5.67	0.19 s
		KPLS3	0.079	0.82 s	4.05	0.5 s	10.41	1.14 s	5.67	0.43 s

"Kriging" refers to the ordinary kriging Optimus solution and "KPLSh" refers to the KPLS model with h principal components

Table 11

 11 Results for tab 3 experiment data (99 input variables, 1 output variable y 1 ) obtained by using 341 training points, 23 validation points, and error given by (16)

	99D	Surrogate	Error (%)	CPU time
	tab 3	Kriging	0.021	20 min 02 s
		KPLS1	0.19	46.6 s
		KPLS2	0.03	2 min 15 s
		KPLS3	0.02	4 min 56 s
	"Kriging" refers to the ordinary kriging optimus solution and
	"KPLSh" refers to the KPLS model with h principal components
	KPLS3 models applied to a 99D problem are very promising
	(see Table

Table 12

 12 Examples of commonly used stationary covariance functions

	Covariance	Expression				Hyper-parameters θ	Number of
	functions								hyper-parameters
									to estimate
	Generalized exponential	σ 2	d i=1	exp(-θ i m	p i i )			(θ 1 , . . . , θ d , p 1 , . . . , p d )	2 d
	Gaussian exponential Matern 5 2 Matern 3 2	σ 2 σ 2 σ 2	d i=1 d i=1 i=1 d	exp(-θ i m 2 i ) 1 + √ 5θ i m i + 5 3 θ 2 i m 2 i exp(-1 + √ 3θ	√	5θ i m i )	(θ 1 , . . . , θ d ) (θ 1 , . . . , θ d )	d d

i m i exp(-√ 3θ i m i ) (θ 1 , . . . , θ d ) d

Table 13

 13 Examples of KPLS covariance functions

	Covariance functions	Expression										Hyper-parameters θ	Number of
															hyper-parameters
															to estimate
	Generalized exponential	σ 2	h	d	exp -θ l m (l) i	p l					(θ 1 , . . . , θ h , p 1 , . . . , p h )	2h	2d
			l=1	i=1
	Gaussian exponential	σ 2	h	d	exp -θ l m	(l) i	2					(θ 1 , . . . , θ h )	h	d
			l=1	i=1
	Matern 5 2	σ 2	h l=1	d i=1	1 +	√	5θ l m (l) i +	5 3	θ 2 l m (l) i	2	exp -√ 5θ l m (l) i	(θ 1 , . . . , θ h )	h	d
	Matern 3 2	σ 2	h	d	1 +	√	3θ l m	(l) i	exp -√	3θ l m	(l) i	(θ 1 , . . . , θ h )	h	d
			l=1	i=1

Table 15

 15 Results for Griewank function in 60D over interval[-5, 5] Ten trials are done for each test(50, 100, 200, and 300 training points)

	Surrogate	Statistic	50 points		100 points		200 points		300 points
			Error (%)	CPU time	Error (%)	CPU time	Error (%)	CPU time	Error (%)	CPU time
	Kriging	Mean	1.39	560.19 s	1.04	920.41 s	0.83	2015.39 s	0.65	2894.56 s
		std	0.15	200.27 s	0.05	231.34 s	0.04	239.11 s	0.03	728.48 s
	KPLS1	Mean	0.92	0.07 s	0.87	0.10 s	0.82	0.37 s	0.79	0.86 s
		std	0.02	0.02 s	0.02	0.007 s	0.02	0.02 s	0.03	0.04 s
	KPLS2	Mean	0.91	0.43 s	0.87	0.66 s	0.78	2.92 s	0.74	1.85 s
		std	0.03	0.54s	0.02	1.06 s	0.02	2.57 s	0.03	0.51 s
	KPLS3	Mean	0.92	1.57 s	0.86	3.87 s	0.78	6.73 s	0.70	20.01 s
		std	0.04	1.98 s	0.02	5.34 s	0.02	10.94 s	0.03	26.59 s

Acknowledgments

The authors thank the anonymous reviewers for their insightful and constructive comments. We also extend our grateful thanks to A. Chiplunkar from ISAE SUPAERO, Toulouse and R. G. Regis from Saint Joseph's University, Philadelphia for their careful correction of the manuscript and to SNECMA for providing the tables of experiment results. Finally, B. Kraabel is gratefully acknowledged for carefully reviewing the paper prior to publication.

Appendix A: Equations for ordinary kriging model

The expression (3) for the ordinary kriging model is transformed into (see [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF]

where 1 denotes an n-vector of ones and

In addition, ( 6) is written as 12.

Appendix B: Examples of kernels

Appendix C: Proof of equivalence kernel

For l = 1, . . . , h, k l are separable kernels (or a d-dimensional tensor product) of the same type, so ∃ φ l1 , . . . , φ ld such that

Appendix D: Results of Griewank function in 20D and 60D over interval [-5, 5]

In Tables 14 and15, the mean and standard deviation (std) of the numerical experiments with the Griewank function are given for 20 and 60 dimensions, respectively. To better visualize the results, boxplots are used in Figs. 7-11.


Elodie Mezzatesta

Annie Perraud

Cathy Vieillescazes

Carole Mathe
email: carole.mathe@univ-avignon.fr.

Analysis of balms taken from Egyptian human mummies using Solid Phase Extraction and GC-MS

Keywords: Balm, Bitumen, Beeswax, Resins, Mummies

The aim of this paper is to establish a protocol by Solid Phase Extraction-Gas Chromatography-Mass Spectrometry leading to a wide and fine qualitative chemical characterization of the several natural substances present in human mummies' balms, using a minimal quantity of samples. In this study, nine samples were analyzed from mummies dating back from the Third Intermediate Period to the Roman Period, and were provided by the Confluences Museum (Lyon, France).

Using Solid Phase Extraction, three fractions were examined in this protocol. The first one, eluted with hexane, concerned chemical families of hydrocarbons of bitumen.

The second, eluted with ethanol, enabled terpenic compounds to be characterized and beeswax. The last one, composed of diethyl ether with 2% of acetic acid, extracted carboxylic acids with a long aliphatic chain (fatty matter) and glycerides.

This study also allowed the characterization of non-saponified compounds from beeswax to be obtained while excluding the common saponification step. The analyzed mummification balms were shown to contain fatty matter, beeswax, bitumen and diterpenic resinous material. This one-pot Solid Phase Extraction-Gas Chromatography-Mass Spectrometry method was efficient in reducing both the number of analytical steps and the complexity of the archaeological balms subsequently analyzed by GC-MS.

Introduction

Mummification was a fundamental part of funerary practice in ancient Egypt. The Egyptians believed in immortality and rebirth in the afterlife and in order to be able to access this afterlife, it was necessary to preserve the deceased body. The application of balms was a key part of the mummification process in ensuring successful and lasting preservation of the body [START_REF] Dunand | Les Momies et La Mort En Egypte[END_REF]. Before applying these balms, however, several steps of the complex process of mummification were carried out.

Firstly, an excerebration was performed, followed by evisceration of the thoracoabdominal cavity, before dehydrating the entire body using natron. The complete process of mummification (including evisceration, dehydration and bandaging) took 70 days; within this period and after dehydration by natron, balms could be added inside the eviscerated cavities (thorax, abdomen) and the skull, and then all over the body which would subsequently be bandaged. However, the entire process could vary for each mummy, depending on various factors. All these stages were described by Herodotus (L'Enquête, II, 85-88) [START_REF] Barguet | Hérodote -Thucydide Oeuvres Complètes[END_REF] and Diodorus Siculus (Bibliothèque historique, I, 21) [START_REF] Vernière | Diodore de Sicile. Bibliothèque Historique. Les belles lettres[END_REF] and could differ according to the period and to the status of the deceased individuals.

Many varieties of natural substances were used in the formulation. These substances have certain pharmacological properties such as being hydrophobic, antibacterial, antifungal and odoriferous. The most frequently used substances were plant oils or animal fats, di-or triterpenic resins, beeswax, bitumen and spices [START_REF] Goyon | Rituels Funéraires de l'ancienne Egypte. Le Rituel de l'Embaumement[END_REF][START_REF] Nicholson | Ancient Egyptian Materials and Technology[END_REF][START_REF] Colombini | Characterization of the Balm of an Egyptian Mummy from the Seventh Century B[END_REF][START_REF] Abdel-Maksoud | A review on the materials used during the mummification processes in ancient Egypt[END_REF][START_REF] Marshall | Les Momies Egyptiennes[END_REF].The fatty substances, composed of esters of fatty acids with long chain between 12 and 24 carbon atoms, were used for their hydrophobic properties, the most abundant generally being palmitic (C16:0) and stearic (C18:0) acids. Sterols such as cholesterol are considered as a marker of animal fat, and campesterol as a marker of plant oils [START_REF] Clark | Tracing the Evolution of Organic Balm Use in Egyptian Mummification via Molecular and Isotopic Signatures[END_REF][START_REF] Bastien | Etude Chimique Des Substances Contenues Dans Une Collection Exceptionnelle de Poteries Provenant de Deir-El-Medineh (Egypte) : Une Population et Ses Produits[END_REF][START_REF] Gea | Characterization of ancient lipids in prehistoric organic residues: Chemical evidence of livestock-pens in rockshelters since early neolithic to bronze age[END_REF].

Plant resins were used in the formulation of balms for their antibacterial properties.

They are composed of terpenic compounds secreted by certain plants [START_REF] Nicholson | Ancient Egyptian Materials and Technology[END_REF][START_REF] Bruni | Identification of archaeological triterpenic resins by the non-separative techniques FTIR and 13C NMR: The case of Pistacia resin (mastic) in comparison with frankincense[END_REF].

Using beeswax limited rehydration of the embalmed body. Beeswax is composed of esters of palmitate with a very long carbon chain alcohols (from 22 to 34 carbon atoms) [START_REF] Régert | Chemical Alteration and Use of Beeswax Through Time: Accelerated Ageing Tests and Analysis of Archaeological Samples from Various Environmental Contexts[END_REF][START_REF] Régert | Characterisation of wax works of art by gas chromatographic procedures[END_REF]. Bitumen was used for its hydrophobic properties provided by three different apolar chemical families corresponding to linear alkanes and triterpenoids such as sterane and hopane compounds [START_REF] Bahn | The making of a mummy[END_REF][START_REF] Harrell | Sources of mummy bitumen in ancient Egypt and Palestine[END_REF][START_REF] Connan | The bituminous mixtures of Kavuşan Höyük (SE Turkey) from the end of the 3rd millennium (2000BC) to the Medieval period (AD 14th century): Composition and origin[END_REF].

Identifying all of the compounds present in each balm plays an important role in the understanding of archaeological data. Furthermore, such identification represents a significant challenge because of the various chemical reactions such as oxidation and hydrolysis, which can occur naturally over time (aging) or because of anthropogenic factors, such as the heating temperature.

Compounds of alteration can play an important part in identifying plant species. The molecular composition of samples can provide information on the nature of the constituents, the state of conservation of the sample and the treatment undergone by the material during the preparation and/or the formation of balms [START_REF] Peters | The Biomarker Guide: Volume 2, Biomarkers and Isotopes in Petroleum Systems and Earth History[END_REF][START_REF] Evershed | Organic residue analysis in archaeology: the archaeological biomarker revolution[END_REF][START_REF] Ménager | Analyse de baumes issus de tête momifiée (Thèbes)[END_REF]. To identify their composition, the most widespread analytical method described in the literature is Gas Chromatography coupled to Mass Spectrometry (GC-MS) [START_REF] Nicholson | Enlightening the past: Analytical proof for the use of Pistacia exudates in ancient Egyptian embalming resins: Gas Chromatography[END_REF][START_REF] Ribechini | An integrated analytical approach for characterizing an organic residue from an archaeological glass bottle recovered in Pompeii (Naples, Italy)[END_REF][START_REF] Bonaduce | Analytical approaches based on gas chromatography mass spectrometry (GC/MS) to study organic materials in artworks and archaeological objects[END_REF][START_REF] Charrié-Duhaut | The canopic jars of Rameses II: real use revealed by molecular study of organic residues[END_REF][START_REF] Bera | Dammar resin from the Eocene lignite of Bengal Basin, eastern India: Terpenoid composition and botanical origin[END_REF][START_REF] Jones | A prehistoric Egyptian mummy: Evidence for an 'embalming recipe' and the evolution of early formative funerary treatments[END_REF][START_REF] Tchapla | Characterisation of embalming materials of a mummy of the Ptolemaic era. Comparison with balms from mummies of different eras[END_REF][START_REF] Hamm | Headspace solid phase microextraction for screening for the presence of resins in Egyptian archaeological samples[END_REF] because it allows a large number of organic compounds from a complex matrix, even in traces, to be separated and identified. However, some chemical families, like bitumen, were difficult to accurately characterize because of their low concentration and of the high proportion of fatty acids present in the balms under study. Thus, because of the presence (i) of a mixture of natural substances, (ii) of alteration processes, (iii) of a very low proportion of some ingredients or, on the contrary, a high concentration of others, the use of GC-MS is often preceded by several sample pre-treatments. A pre-treatment sample was often necessary, such as saponification reactions for the detection of waxes [START_REF] Łucejko | Chemical analyses of Egyptian mummification balms and organic residues from storage jars dated from the Old Kingdom to the Copto-Byzantine period[END_REF][START_REF] Degano | Multi-analytical techniques for the study of pre-Columbian mummies and related funerary materials[END_REF], dichloromethane or diethyl ether extractions for the detection of free fatty acids and terpenic compounds, as well as fractionation for hydrocarbons [START_REF] Colombini | Characterization of the Balm of an Egyptian Mummy from the Seventh Century B[END_REF][START_REF] Jones | A prehistoric Egyptian mummy: Evidence for an 'embalming recipe' and the evolution of early formative funerary treatments[END_REF][START_REF] Łucejko | Chemical analyses of Egyptian mummification balms and organic residues from storage jars dated from the Old Kingdom to the Copto-Byzantine period[END_REF][START_REF] Łucejko | An analytical approach based on X-ray diffraction, Fourier transform infrared spectroscopy and gas chromatography/mass spectrometry to characterize Egyptian embalming materials[END_REF][START_REF] Ménager | Study of egyptian mummification balms by FT-IR spectroscopy[END_REF][33][START_REF] Goyon | A propos de la momie infantile du musée Calvet d'Avignon[END_REF][START_REF] Brettell | Organic residue analysis of Egyptian votive mummies and their research potential[END_REF]. All of these pre-treatments enable a fine chemical characterization of samples to be obtained, but pre-treatments also create the need for more archaeological material mass, when the quantities available for analysis are generally small.

The aim of this study therefore, was to develop a single and efficient protocol by SPE-GC-MS. The main objective of this paper is to determine the fine chemical composition of nine mummification balms and to reduce the number of experimental steps, thus reducing the amount of archaeological sample used thanks to a one-pot procedure.

Materials and Methods

Archaeological samples

Within the framework of a multidisciplinary research project entitled Human Egyptian Lyon Confluences Mummies (HELYCOM)-Mourir pour renaître, mummification balms of 9 human mummies were analyzed. In this research work, female and male mummies from Upper Egypt (around Thebes) dating from 600 BC to the Roman Period (300 AD) were studied. In this study, a total of 9 samples belonging to the Egyptian collection of the Confluences Museum in Lyon (France) were examined.

The samples referred to as 30000111 [START_REF] Goyon | A propos de la momie infantile du musée Calvet d'Avignon[END_REF], 30000139 [START_REF] Bahn | The making of a mummy[END_REF], 30000148 [START_REF] Peters | The Biomarker Guide: Volume 2, Biomarkers and Isotopes in Petroleum Systems and Earth History[END_REF], 30000286 [START_REF] Gea | Characterization of ancient lipids in prehistoric organic residues: Chemical evidence of livestock-pens in rockshelters since early neolithic to bronze age[END_REF] and 90001951A (50) were mummified heads, and 90001169 [START_REF] Scalarone | Ageing behaviour and pyrolytic characterisation of diterpenic resins used as art materials: colophony and Venice turpentine[END_REF], 90001258 [START_REF] Carpy | Molecular characterization of retene derivatives obtained by thermal treatment of abietane skeleton diterpenoids[END_REF], 90001259 (47) and 900015978 (66) were provided from complete mummies. The mummies studied were those of men and women dating from the Third Intermediate Period to the Roman Period and coming from different necropolises such as Thebes, Esna, Kom-Ombo and Deir-el-Medineh. The weight of the samples varied from 20 to 200 mg. They were collected at different location points on the mummies: on the mummified skin, on the textile covered by balms, or directly within the balm itself (Table 1).

Solvents and reagents

All solvents were of the highest purity grade. Hexane, tetrahydrofuran (THF) and N,O-Bis (trimethylsilyl)trifluoroacetamide/Trimethylchlorosilane (BSTFA/TMCS) were supplied by Sigma-Aldrich. Ethanol, dichloromethane (DCM) and diethylether (DEE) were supplied by Merck.

Saponification

10 mg of sample was extracted with 3 x 1 mL of THF aided by sonification (5 min) and then centrifuged at 6000 rpm (5 min). The supernatant was set aside to extract the solid pellet again. The solvent extracts were combined and then 2 mL of a solution of potassium hydroxide KOH 10% in MeOH/H2O (9/1, v/v) were added. The mixture was magnetically stirred and heated at 65°C during 1 hour. After evaporation 3 mL of pure water was added with 1 mL of HCl 5 M. The aqueous phase was washed with 3 x 5 mL of diethyl ether. The organic phases were combined and dried with anhydrous sodium sulphate, then filtered on filter paper. Excess reagent was evaporated to dryness under a stream of nitrogen. Trimethylsilylation was applied.

After evaporation, the derivatized sample was solubilized in hexane/DCM mixture (2/1, v/v) and filtered on PTFE cartridge before injection in GC-MS.

Solid Phase Extraction

Preliminary tests were realized to determine the best SPE conditions. In this study, according to the apolar properties to the wide chemical compounds present in human mummies' balms, three types of SPE columns were tested: C18, Cyano and Amino phases associated to different solvents of elution (hexane, methanol, ethanol, isopropanol, diethyl ether with 2% of acetic acid and dichloromethane). The NH 2 SPE column associated to hexane, ethanol, diethyl ether with 2% acetic acid exhibited the best results so this column was used in this study according to the following protocol.

10 mg of the sample were extracted with 1 mL of hexane/tetrahydrofuran (THF) (1/1) using ultrasound for five minutes and then centrifuged at 6000 rpm for five minutes.

The supernatant was set aside to extract the solid pellet a second and third time. The three fractions obtained were combined and evaporated to dryness under a stream of nitrogen and the mixture was then dissolved in 500 μL hexane/THF (1/1, v/v). This mixture was called the charge.

In parallel, an SPE cartridge strata NH 2 200 mg/3 mL (Phenomenex) was conditioned in an SPE vacuum manifold with 4 mL of hexane at a flow of 1.2 mL.min -1 . The charge was deposited on the cartridge. Then a first elution was carried out with 4 mL of hexane and collected (fraction 1). A second elution was performed with 4 mL of ethanol and was collected (fraction 2). Finally, a last elution was carried out with 3 mL of diethyl ether (DEE) including 2% acetic acid (fraction 3). Each fraction was evaporated to dryness under a stream of nitrogen and then trimethylsilylated with 200 µL of BSTFA/TMCS at 70°C during 30 min. Fractions 1 and 2 were solubilized with 60 μL hexane/DCM (2/1, v/v) and were injected into GC-MS at temperature gradient No.1. Fraction 3 was solubilized with 1.5 mL hexane/DCM (2/1) and injected into GC-MS by using gradient program No.2.

Each SPE parameter, namely the nature of the sorbent, the dimension of the cartridge, the sampling flow rate, the sample volume, the eluents and the desorption flow rate, was optimized.

GC-MS conditions

A Thermo Scientific Focus gas chromatographic system composed of a Thermo Scientific Al 3000 auto-sampler coupled with an ITQ 700 ion trap mass spectrometer was used (Thermo Fisher Scientific). The GC column was a fused silica capillary column Thermo trace GOLD TG-5MS (5% diphenyl / 95% dimethylpolysiloxane, 30 m length × 0.25 mm i.d. × 0.25 μm film thickness). Helium was the carrier gas with a constant flow of 1 mL.min -1 . 1 μL of each sample was injected with a splitless time of 1 min. The injector temperature was set at 250 °C. Mass spectra were recorded in electron impact mode with an electron ionization voltage of 70 eV, an ionization time of 25,000 μs and a mass range of 40-650 m/z. The transfer line, ion trap and manifold temperatures were respectively set at 300 °C, 200 °C and 50 °C. Detected compounds were identified from their retention time and interpretation of their mass spectra in comparison with standard compounds and using NIST database. The oven temperature was programmed as follows:

-Temperature gradient No.1: 160 °C, isothermal for one minute, then increased by 10 °C.min -1 up to 180 °C, constant for three minutes, followed by a second increase of 3 °C.min -1 up to 260°C, constant for five minutes, followed by a third increase of 2.5°C.min -1 up to 300°C and finally, a final increase of 10°C.min -1 up to 320°C, the temperature of which was maintained for three minutes. The total analysis time was fifty-nine minutes.

-Temperature gradient No.2: 115°C, constant for two minutes, then increased by 15°C.min -1 up to 220°C, followed by a second increase of 3°C.min -1 up to 250°C and finally, a last increase of 10°C.min -1 up to 320°C, the temperature of which was maintained for three minutes. The total analysis time was twentynine minutes.

Results and discussion

Owing to (i) the hypothetical occurrence of well-known natural substances, such as fatty matter or plant resins and (ii) to the presence of alcohol and/or carboxylic functional groups, the selected cartridge was the Strata NH 2 . It effectively allows the retention of compounds depending on the electrostatic attractions of the functional groups of the compounds of interest to the NH 2 grouping, bound at the surface of the SPE cartridge. The SPE protocol was carried out according to the schema described in Figure 1. All of the molecules detected in the nine archaeological samples are listed in Table 2.

Tetrahydrofuran solubilization

Solubilization of the sample is necessary beforehand, since depositing the solid matter directly on the SPE cartridge percolates the cartridge and inhibits elution. The solvent must facilitate the most exhaustive possible extraction of the compounds of interest, both qualitatively and quantitatively. Under these experimental conditions, the solvent should be moderately apolar in order to extract the compounds present in the balms. Several tests were carried out, comparing them with solvents usually described in the specialized literature, such as dichloromethane and diethyl ether.

Tetrahydrofuran/hexane (1/1) was chosen owing to a polarity close to that of chloroform or dichloromethane, but it is a water-miscible solvent and it should therefore have superior capacities for the solubilization of polar molecules [START_REF] Mezzatesta | Caractérisation Moléculaire de Baumes de Momies Humaines d'Egypte Ancienne[END_REF].

Fraction 1

The first eluted fraction was carried out using hexane, which is the most apolar solvent used in this study. This elution was start with an apolar solvent because the stationary phase is a polar phase, moreover, hexane is also the conditioning solvent of the cartridge. Hexane was able to elute the most apolar molecules present in the studied balms, such as aliphatic molecules, linear, branched, cyclized alkanes or polycyclic aromatic hydrocarbons (PAHs). An illustration of the obtained results in Sample No.18 is presented in Figure 2. The detection of the three chemical families of bitumen at m/z 57 (alkanes), m/z 191 (hopane) and m/z 217 (sterane) can be used for the characterization of bitumen [START_REF] Clark | Tracing the Evolution of Organic Balm Use in Egyptian Mummification via Molecular and Isotopic Signatures[END_REF]. The chemical families of hopanes and steranes existed in only a very low proportion in bitumen (5%), the characterization of these chemical families thus usually requiring special treatment [START_REF] Connan | La Momification Dans l'Egypte Ancienne: Le Bitume et Les Autres Ingrédients Organiques Des Baumes de Momies Ou Les Ingrédients Organiques Des Baumes de Momies Égyptiennes: Bitume, Cire d'abeillerésines[END_REF]. Identification of archaeological bitumen was generally performed by isolating and splitting the several compounds depending on their nature: saturated hydrocarbons, aromatic hydrocarbons, and other compounds using column chromatography [START_REF] Charrié-Duhaut | Molecular and isotopic archaeology: Top grade tools to investigate organic archaeological materials[END_REF][START_REF] Hauck | Molecular evidence of bitumen in the Mousterian lithic assemblage of Hummal (Central Syria)[END_REF].

The GC-MS analyses of Fraction 1 of the samples referred to as 34, 15, 18, 41 and 50 resulted in the detection of occurrence of the three chemical families of bitumen: alkanes, hopanes and steranes. These results confirmed the presence of bitumen in the balms under study.

In addition to these results, a PAH was detected and identified to retene. An anthropic degradation of Pinaceae resin such as thermal degradation can be deduced from the detection of retene. This latter molecule was formed mainly under very high temperatures [START_REF] Carpy | Molecular characterization of retene derivatives obtained by thermal treatment of abietane skeleton diterpenoids[END_REF] and its occurrence was observed in 5 balms and in a significant proportion in the balm of Sample 11. Contrary to the other diterpenic compounds eluted in the fraction 2, retene was identified in the fraction 1 because of its apolar property.

Fraction 2

Fraction 2, eluted with ethanol enabled most polar compounds such as diterpenoids, which are functionalized compounds, to be eluted, with the occurrence of alcohol, carboxylic acid, ketone, or aldehyde functions [START_REF] Clark | Tracing the Evolution of Organic Balm Use in Egyptian Mummification via Molecular and Isotopic Signatures[END_REF]. Diterpenic resin was detected in 8 out of 9 samples (34, 15, 18, 11, 41, 40, 47 and 66). For example, the sample referred to as 18 had many diterpenoids, which are characteristic of the use of Pinaceae resin, such as in the elution order: dehydroabietic acid (DHA, Rt= 23.13 min), 7-hydroxy-DHA (Rt=24.5 min) and 7-oxo-DHA (Rt= 26.3 min) (Figure 3). During the aging process, an initial predominant isomerization produces dehydroabietic acid. A natural degradation of this acid could occur the formation of products such as 7-oxo-DHA, 7-oxo-15-hydroxy-DHA [START_REF] Scalarone | Ageing behaviour and pyrolytic characterisation of diterpenic resins used as art materials: colophony and Venice turpentine[END_REF][START_REF] Romero-Noguera | Study of biodeterioration of diterpenic varnishes used in art painting: Colophony and Venetian turpentine[END_REF][START_REF] Mezzatesta | GC-MS and PCA analyses of diterpenoids degradation state in 21 human mummies of Ancient Egypt dating from New Kingdom to Graeco-Roman Period[END_REF].

Additionally, to the diterpenoids already identified, another compound, dehydroabietic methyl ester (DHAM, Rt = 20.5 min), was characterized in both of the samples referred to as 18 and 11. This molecule is produced by methanolysis. In fact, during heating at a very high temperature, resinous wood releases methanol, which reacts with the carboxylic acids of diterpenoids and this reaction leads to the formation of the corresponding methyl ester derivatives. DHAM allows wood tar to be characterized and this compound is totally absent when resin alone is heated [44].

Moreover, palmitate esters corresponding to unsaponified beeswax molecules were identified in this fraction (Rt= 55-58 min). Beeswax does not necessarily undergo hydrolysis with aging, and its molecules therefore can remain in their ester form [45].

In this eluted fraction, the presence of beeswax was directly characterized in samples 34, 15, 18, 66, 50 and 47. The characterization of beeswax is usually performed by initially carrying out a pre-treatment of the samples corresponding to a reaction of saponification [START_REF] Łucejko | Chemical analyses of Egyptian mummification balms and organic residues from storage jars dated from the Old Kingdom to the Copto-Byzantine period[END_REF]. With the purpose of confirming either presence or absence of beeswax, saponification was performed as a preliminary step. Results obtained through saponification are in accordance with aforementioned SPE-GC-MS data [START_REF] Mezzatesta | Caractérisation Moléculaire de Baumes de Momies Humaines d'Egypte Ancienne[END_REF].

In order to interpret the elution phenomenon of these apolar compounds in ethanol, it is necessary to consider the state of ionization of the several compounds and also to consider the occurrence of free fatty acids in the mixture. In this solution, the stationary amine phase was in the NH 3 + form and free fatty acids were in the COO - form, enabling ionic bonds therefore to be formed between the two forms and thus creating a type of semi-apolar grafted phase in the cartridge. Beeswax esters could form Van der Waals bonds (London interactions) with the carbon chain of free fatty acids, which would explain why they are not eluted with hexane in Fraction 1. The subsequent addition of ethanol, a protic solvent, capable thus of generating hydrogen bonds, could create an imbalance leading to the elution of esters.

Fraction 3

Fraction 3 was eluted with a mixture of diethyl ether with 2% acetic acid. In this experimental condition, carboxylic acids presented COOH functions, and the amine of the stationary phase was in NH 3 + form. No ionic bond was formed between these two chemical forms. For this reason, this fraction allowed fatty acids to be eluted. All samples contained monocarboxylic fatty acids (MC) and/or dicarboxylic acids (DC).

An illustration of the obtained results in sample 18 is presented in Figure 4. In the obtained GC-MS chromatograms, the relative proportion of MC represented on average, approximately 50% of the total composition of balms and this percentage could attain 90% depending on the archaeological sample studied. This relative percentage was calculated from the peak area of each detected compound from corresponding chromatograms. Monocarboxylic acids corresponded also to the main chemical compounds present in all of the analyzed balms.

Saturated monocarboxylic acids cannot provide much information in the archaeological context about the origin of the fat matter used. In fresh material, it is possible to examine the peak area ratio of MC 16:0/MC 18:0 (palmitic acid/stearic acid) to determine the nature of the fatty matter [START_REF] Clark | Tracing the Evolution of Organic Balm Use in Egyptian Mummification via Molecular and Isotopic Signatures[END_REF]. However, this method is not applicable in an archaeological context because the proportions of these molecules can change with the state of degradation of the material or can even be caused by a mixture of oils and fats. This method is therefore not suitable for the analysis of mummy balms. The same is true for the proportions of saturated fatty acids with a short carbon chain. They fail to provide much information because these molecules are certainly the result of fat/oil degradations.

In this study using GC-MS, therefore, the animal or plant nature of fats cannot be accurately defined because biomarkers of plant oils and/or animal fats were not

Conclusion

The studied mummification balms contained fatty matter, beeswax, bitumen and diterpenic resin. This efficient SPE-GC-MS protocol has thus enabled additional substances to be identified, substances such as hopanes and steranes of bitumen and palmitate esters of beeswax. It has also enabled markers of wood tar such as DHAM compound to be identified and has reduced the number of experiments such as the long pre-treatment of saponification or the splitting of the constituents of bitumen. The quantity of samples has also been reduced, which is very interesting from an archaeological point of view given the very small amount of material available, in accordance with cultural heritage directives concerning preservation.

This developed SPE protocol as applied to human balms of mummies has enabled each ingredient of the studied balms to be separated and concentrated. Moreover, this multiple extraction/partial elution, carried out during a single experimentation, has resulted in better pre-separation of the compounds deriving from a large panel of natural substances used in Ancient Egypt. The SPE protocol described in this study was also successfully applied and performed in the analyses of 61 samples taken from 43 human mummies from Upper Egypt [START_REF] Mezzatesta | Caractérisation Moléculaire de Baumes de Momies Humaines d'Egypte Ancienne[END_REF].

[44] Colombini, M. P., Modugno, F., Ribechini, E., GC/MS in the characterization of lipids, in: Colombini, Modugno (Eds.) Organic Mass Spectrometry in Art and Archaeology. 2009, pp 189-213.

[45] Čížová, K., Vizárová, K., Ház, A., Vykydalová, A., Cibulková, Z., Šimon, P., Study of the degradation of beeswax taken from a real artefact. J Cult Herit.

2019, 37, 103-112.

  detected. Owing to the use of GC-MS in archaeological context, the objective of this paper was not to determine the nature or the origin of the fatty matter (i.e. fatty acids and glycerides). Cholesterol was only detected in the samples referred to as 34, 11 and 40. The occurrence of this molecule did not allow animal origins to be identified with certainty because the occurrence could result from an external contamination of the deceased by capillarity or of a contemporary contamination during the transport or the storage of the objects. Samples 15, 18, 41, 47 66 and 50 contained numerous dicarboxylic acids. These types of acids were markers of degradation of fatty substances. Unsaturated fatty acids can undergo oxidation and form hydroperoxide intermediates to finally generate short monocarboxylic acids and dicarboxylic acids[START_REF] Bastien | Etude Chimique Des Substances Contenues Dans Une Collection Exceptionnelle de Poteries Provenant de Deir-El-Medineh (Egypte) : Une Population et Ses Produits[END_REF].

Figure 1 :

 1 Figure 1: Schema of the final Solid Phase Extraction protocol developed with the predominantly detected substances in each elution fraction.

Figure 2 :

 2 Figure 2: Chromatograms of sample 18, Fraction 1 and extracted signal at m/z 57, m/z 191 and m/z 217, respectively base peaks of alkanes, hopanes and steranes.

Figure 3 :

 3 Figure 3: TIC chromatogram of sample 18, Fraction 2.

Figure 4 :

 4 Figure 4: TIC chromatogram of sample 41, Fraction 3.

Figure 1 :

 1 Figure 1: Schema of the final Solid Phase Extraction protocol developed with the predominantly detected substances in each elution fraction.

Figure 2 :

 2 Figure 2: TIC partial chromatogram of sample 18, Fraction 1 and extracted signal at m/z 57, m/z 191 and m/z 217, respectively base peaks of alkanes, hopanes and steranes.

Figure 3 :

 3 Figure 3: TIC Partial chromatogram of sample 18, Fraction 2.

Figure 4 :

 4 Figure 4: TIC partial chromatogram of sample 41, Fraction 3

Table 1 :

 1 Sample presentation 528

	Inventory no.	Lab code	Sex and age	Site excavation of	Relative dating
	30000111		Elderly male	Gournah	Graeco-Roman Period
	30000139		Adult female	Gournah	Ptolemaic to Roman Period
	30000148		Elderly male	Thebes	Unknown
	30000286		Adult or elderly male	Esna	Third Intermediate Period to Late Period
	90001951A		Young male	adult	Thebes	Late Period to Ptolemaic Period
	90001169		Adult male		Gournah	Graeco-Roman Period
	90001258		Adult female	Kom-Ombo Ptolemaic to Roman Period
	90001259		Adult female	Unknown	Late Period to Ptolemaic Period
	90001597		Adult		Gournah	Unknown
	529

Acknowledgments

The authors acknowledge the anonymous reviewers for their careful reading of our manuscript and their many insightful comments and suggestions. The authors wish to thank the Museum of Confluences in Lyon and the LabEx ARCHIMEDE and Human Egyptian LYon COnfluences Mummies (HELYCOM) -Mourir pour renaître project. A part of this project was supported by LabEx ARCHIMEDE from the Investissement d'Avenir program ANR-11-LABX-0032-01. The research has been funded by a PhD grant given by The French Ministry of Higher Education and Research.

Conflict of interest

The authors declare that they have no conflict of interest.