M Binois
email: binois@emse.fr

D Ginsbourger
email: ginsbourger@stat.unibe.ch

O Roustant
email: roustant@emse.fr

Quantifying uncertainty on Pareto fronts with Gaussian Process conditional simulations

Keywords: Multi-objective optimization, Attainment function, Vorob'ev Expectation, Expected Hypervolume Improvement, Kriging

Multi-objective optimization algorithms aim at finding Pareto-optimal solutions. Recovering Pareto fronts or Pareto sets from a limited number of function evaluations are challenging problems. A popular approach in the case of expensive-to-evaluate functions is to appeal to metamodels. Kriging has been shown efficient as a base for sequential multi-objective optimization, notably through infill sampling criteria balancing exploitation and exploration such as the Expected Hypervolume Improvement. Here we consider Kriging metamodels not only for selecting new points, but as a tool for estimating the whole Pareto Front and quantifying how much uncertainty remains on it at any stage of Kriging-based multi-objective optimization algorithms. Our approach relies on the Gaussian Process interpretation of Kriging, and bases upon conditional simulations. Using concepts from random set theory, we propose to adapt the Vorob'ev expectation and deviation to capture the variability of the set of non-dominated points. Numerical experiments illustrate the potential of the proposed workflow, and it is shown on examples how Gaussian process simulations and the estimated Vorob'ev deviation can be used to monitor the ability of Kriging-based multi-objective optimization algorithms to accurately learn the Pareto front.

Introduction

The interest in Multi-Objective Optimization (MOO) has been growing over the last decades, resulting in the development of numerous dedicated methods, especially in evolutionary MOO [START_REF] Deb | Introduction to evolutionary multiobjective optimization[END_REF]. These methods are able to cope with challenging problems occurring when few information about the properties of the objective functions is available (black-box optimization). In some situations, as for example in car crash safety design [START_REF] Liao | Multiobjective optimization for crash safety design of vehicles using stepwise regression model[END_REF], another difficulty comes from a limited budget of evaluations, because of expensive experiments or high fidelity simulations.

In this context, see e.g. [START_REF] Santana-Quintero | A review of techniques for handling expensive functions in evolutionary multiobjective optimization[END_REF] for a review, a common approach is to rely on a surrogate model or metamodel to alleviate the computational costs of the optimization process. In particular, Kriging metamodels have proven to be efficient because they not only give a response surface but also a quantification of prediction uncertainty. In mono-objective optimization, this property has been extensively used following the principles of the Efficient Global Optimization (EGO) algorithm [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] where the Expected Improvement criterion is used to balance between exploitation and exploration. Extensions to MOO have been developed, from scalarization approaches [START_REF] Knowles | ParEGO: a hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems[END_REF][START_REF] Zhang | Expensive multiobjective optimization by MOEA/D with Gaussian process model[END_REF] to the use of multi-objective expected improvement criteria such as the Expected Hypervolume Improvement [START_REF] Emmerich | Hypervolume-based expected improvement: Monotonicity properties and exact computation[END_REF].

While results about the optimality of solutions from aggregation approaches have been reported (see e.g. [START_REF] Miettinen | Nonlinear multiobjective optimization[END_REF]), things are more difficult to analyze for MOO and even worse in metamodel based MOO, where an additional source of uncertainty due to surrogate modeling must be taken into account. Besides, the study of the convergence in evolutionary MOO is an ongoing subject of research [START_REF] Wagner | A taxonomy of online stopping criteria for multi-objective evolutionary algorithms[END_REF].

Inspired by what has been proposed for Kriging-based excursion sets estimation in Chevalier et al. [START_REF] Chevalier | Estimating and quantifying uncertainties on level sets using the Vorobev expectation and deviation with Gaussian process models[END_REF][START_REF] Chevalier | Fast uncertainty reduction strategies relying on Gaussian process models[END_REF], we propose here to use notions from the theory of random sets [START_REF] Molchanov | Theory of random sets[END_REF] for quantifying uncertainty on Pareto fronts, through conditional simulations. The latter are used to estimate the probability that any given point in the objective space is dominated, which is known in performance assessment of multi-objective optimizers as the attainment function [START_REF] Grunert Da Fonseca | The attainment-function approach to stochastic multiobjective optimizer assessment and comparison[END_REF]. From this we obtain a metamodel-based estimation of the Pareto front using the Vorob'ev expectation [START_REF] Molchanov | Theory of random sets[END_REF], with a value of the associated uncertainty: the Vorob'ev deviation. At each stage of the sequential optimization process, an insight is provided to the practitioner about convergence and possibilities of further improvements. Furthermore, with two or three objectives the proposed approach makes it possible to visualize the variability around the estimation of the Pareto front in the objective space.

The paper is organized as follows: Section 2 details notions in MOO and in Gaussian Process Regression upon which the proposed approach is based. In particular Section 2.4 is dedicated to conditional simulations. In Section 3, we propose an original definition of uncertainty using the Vorob'ev expectation and deviation. Finally, Section 4 is dedicated to applications of the proposed methodology to three different test cases, where the potential of the approach to quantify uncertainty and monitor convergence within a sequential MOO algorithm is illustrated.

Multi-objective optimization and Gaussian Process Regression

Notions in MOO

Multi-objective optimizers aim at optimizing (say minimizing) several objectives at once: f 1 (x), . . . , f m (x) with x = (x 1 , . . . , x d ) T a vector of decision variables in E (usually E ⊂ R d ) and f : E → R m the vector valued function whose coordinates are f 1 , . . . , f m . As the objectives are usually in competition, the existence of an optimal solution minimizing all objectives simultaneously cannot generally be taken for granted. This leads to the definition of compromise solutions following the Pareto dominance: a vector is said to be dominated if there exists another vector which is not worse in any objective and better for at least one of them. If a vector is not dominated by any other vector, it is said to be optimal in the Pareto sense.

The set of optimal (or non-dominated) points in E is called Pareto set and the corresponding image by f , composed of non-dominated vectors, is called Pareto front. Multi-objective optimization algorithms aim at finding non-dominated objective vectors as close as possible to the true underlying Pareto front, creating a discrete approximation sometimes called a Pareto front approximation [START_REF] Zitzler | Quality assessment of Pareto set approximations[END_REF].

Kriging / Gaussian Process Regression

A common solution to perform optimization under a tight evaluation budget is to appeal to a mathematical surrogate of the objective function. Here we focus on a class of probabilistic metamodels relying on Gaussian random fields. Originating from geostatistics with a technique named Kriging [START_REF] Matheron | Principles of geostatistics[END_REF], predicting with such a metamodel is known in the machine learning community as Gaussian Process Regression (GPR) [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. These Kriging/GPR metamodels have the property of interpolating the observations when noiseless data is considered (deterministic case). Furthermore, due to the probabilistic nature of these metamodels, they also provide a quantification of the prediction uncertainty.

Without loss of generality, here we do not assume a priori any stochastic dependency between the responses f 1 , . . . , f m and we treat them separately since the use of dependent models is significantly more cumbersome and has not been shown to perform better in state of the art studies [START_REF] Kleijnen | Multivariate versus univariate kriging metamodels for multi-response simulation models[END_REF][START_REF] Svenson | Multiobjective optimization of expensive black-box functions via expected maximin improvement[END_REF]. Following the settings of Gaussian Process Regression [START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF][START_REF] Sacks | Design and analysis of computer experiments[END_REF], each of the objective functions f i is supposed to be a sample path of a random field Y i :

Y i (.) = g T i (.)β i + Z i (.) (Universal Kriging)
where g i (.) T is a vector of known basis functions, β i a vector of unknown coefficient and Z i (.) is a zero mean Gaussian process (GP) with given covariance function, or kernel, k i . With n evaluations at the same locations for the objectives {Y i (x 1 ) = y i,1 , . . . , Y i (x n ) = y i,n , 1 ≤ i ≤ m} denoted A n , the predictor (or Kriging mean) and the prediction covariance (also referred to as Kriging covariance) of Universal Kriging (UK) are expressed as:

m i,n (x) = g i (x) T βi + k i,n (x) T K -1 i,n y i,n -G i,n βi c i,n (x, x') = k i (x, x') -k i,n (x) T K -1 i,n k i,n (x') + g i (x) T -k i,n (x) T K -1 i,n G i,n G T i,n K -1 i,n G i,n -1 g i (x') T -k i,n (x') T K -1 i,n G i,n T where y i,n = (y i,1 , . . . , y i,n ), K i,n = (k i (x s , x t )) 1≤s,t≤n , k i,n (x) = (k i (x, x 1 ), . . . , k i (x, x n )) T , G i,n = g i (x 1 ) T , . . . , g i (x n ) T T and βi = G T i,n K -1 i,n G i,n -1 G T i,n K -1 i,n y i,n .
Note that from a Bayesian point of view, assuming that the Y i are Gaussian conditionally on β i and putting an improper uniform prior on β i , it is known [START_REF] Handcock | A bayesian analysis of kriging[END_REF] that the Universal Kriging mean and covariance coincide with the conditional expectation and covariance of Y i knowing A n , respectively:

m i,n (x) = E(Y i (x)|A n ) and c i,n (x, x') = cov (Y i (x), Y i (x')|A n ).
The covariance functions are chosen according to prior hypothesis about the unknown functions, such as regularity, sparsity, possible symmetries, etc. [START_REF] Ginsbourger | Invariances of random fields paths, with applications in Gaussian process regression[END_REF]. While there exists a variety of admissible covariance functions, the most commonly used are the stationary "Gaussian" and "Matérn" kernels [START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF]. Maximum likelihood estimation or cross validation techniques [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of gaussian processes with model misspecification[END_REF] are typically employed to estimate values for the kernel hyperparameters [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. An example of a Kriging model with constant unknown trend and Matérn (ν = 5/2) kernel is depicted in Figure 1a.

Multi-objective expected improvement

Sequential approaches in MOO aim at adding new observations with a balance between exploration and exploitation. Similar to [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF], several extensions of the EGO algorithm have been proposed for MOO. The main idea is to derive criteria in the vein of the Expected Improvement by defining a generalization of the notion of improvement for multiple objectives. Popular methods include scalarization approaches like ParEGO [START_REF] Knowles | ParEGO: a hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems[END_REF] or MOEAD-EGO [START_REF] Zhang | Expensive multiobjective optimization by MOEA/D with Gaussian process model[END_REF] or truly multi-objective methods based on the definition of improvement functions over the current Pareto front P n defined by the current observations. Considered improvement functions are respectively based on Euclidean distance [START_REF] Keane | Statistical improvement criteria for use in multiobjective design optimization[END_REF], Hypervolume [START_REF] Emmerich | Single-and multiobjective evolutionary optimization assisted by Gaussian random field metamodels[END_REF][START_REF] Wagner | On expected-improvement criteria for model-based multi-objective optimization[END_REF] or Maximin distance [START_REF] Bautista | A sequential design for approximating the pareto front using the expected pareto improvement function[END_REF][START_REF] Svenson | Multiobjective optimization of expensive black-box functions via expected maximin improvement[END_REF] i.e. respectively the distance to the closest point of P n , the volume added over P n and an axis-wise distance to P n .

In the applications of Section 4, we use the Expected Hypervolume Improvement to sequentially add points. This criterion has been successfully applied to problems with limited budget [START_REF] Emmerich | Hypervolume-based expected improvement: Monotonicity properties and exact computation[END_REF], enjoys some beneficial properties [START_REF] Emmerich | Hypervolume-based expected improvement: Monotonicity properties and exact computation[END_REF] and furthermore is related to the concept of attainment function which is of particular importance in what follows.

Conditional simulations

From the Universal Kriging metamodels presented in Section 2.2, it is possible to generate samples interpolating the available evaluation results, called conditional simulations. They can be generated using a variety of methods, from matrix decomposition to spectral methods, as presented in [START_REF] Journel | Geostatistics for conditional simulation of ore bodies[END_REF][START_REF] Hoshiya | Kriging and conditional simulation of Gaussian field[END_REF][START_REF] Diggle | Model-based geostatistics[END_REF][START_REF] Roustant | DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF]. Examples of such conditional simulations are displayed in Figure 1b.

They have been applied in mono-objective optimization in [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF] as a tool to estimate an information gain when no analytical formula is available, as opposed to the Expected Improvement. Until now, the computation of multi-objective Expected Improvement relies either on analytical formulas or on Monte Carlo estimation with draws of the posterior distribution at a single location x. In contrast, conditional simulations consist in drawing posterior realizations of the unknown function at multiple points, say E p : {e 1 , . . . , e p } ⊂ E. Since exact methods essentially depend on p × p covariance matrices, the number of simulation points is typically limited by storage and computational cost. Despite this limitation, conditional simulations (based on matrix decomposition) prove useful for Pareto front estimation, as presented in the next sections.

Quantification of uncertainty

In this section, we assume that a Gaussian process model (see Section 2.2) has been estimated for each objective function from a set of n observations A n . These models allow us to generate conditional Pareto front realizations and further estimate the uncertainty on the Pareto front with concepts from random sets theory.

Conditional simulations for MOO: generation of conditional Pareto fronts and corresponding attained sets

Here we use conditional simulations to generate so-called conditional Pareto fronts (CPF). The first step is to simulate a finite number (say N ) of vector-valued GP conditional simulations

Y (1) 1 , . . . , Y (1) m , . . . , Y (N ) 1 , . . . , Y (N )
m at some simulation points in the design space. Selecting the corresponding non-dominated simulation points and simulated responses then provides conditional Pareto sets and fronts as summarized in Algorithm 1 and illustrated in Figure 2.

Note that what we denote by CPF are actually approximations of conditional Pareto fronts, just like conditional simulations of Gaussian random fields are often approximated realizations relying on a finite number of points. Simulation points can be fixed for all the simulations or changed from one simulation to the other. Fixed simulation points accelerate the simulation generation but they may introduce a bias and lead to missing worthwhile areas. On the other hand, modifying the simulation points increases the computational burden but is more exploratory, which might be an asset in high dimensions. Besides, the procedure used to choose the location of simulation points may impact the results. Accordingly, two sampling strategies are investigated in Section 4.2.

Y (i) j = Y (i) j (e 1 ), . . . , Y (i)
j (e p ) . end for Determine the Pareto set and front of Y

(i) 1 , . . . , Y (i) m .

end for

From now on we focus on the use of CPFs since the decision maker is mostly interested in visualizing results in the objective space. Each CPF is composed of non-dominated points in the objective space. They have been considered to assess the performances of MO optimizers [START_REF] Fonseca | Exploring the performance of stochastic multiobjective optimisers with the second-order attainment function[END_REF][START_REF] Zitzler | Quality assessment of Pareto set approximations[END_REF] under the term Random Non-dominated Point (RNP) sets: sets of random vectors in R m , non-dominated with respect to each other and with random finite cardinality (see e.g. [START_REF] Grunert Da Fonseca | The attainment-function approach to stochastic multiobjective optimizer assessment and comparison[END_REF]). An alternative view is to consider the set of all objective vectors dominated by a realization of an RNP set, called an attained set. Realizations of RNP sets and the corresponding attained sets are presented in Figure 3.

In the mono-objective case, GP models provide analytical expressions of the expectation (Kriging mean) and uncertainty (variance of the pointwise prediction). It would be interesting to get their counterpart for the attained sets of simulated CPFs. Nevertheless, defining an expectation and/or an index of variability for sets is not straightforward and requires concepts from random sets theory [START_REF] Molchanov | Theory of random sets[END_REF].

Basics from random sets theory: quantifying uncertainty with the Vorob'ev deviation

Before introducing related notions for CPFs, let us recall some general definitions. Set-valued random elements, in particular random closed sets [START_REF] Molchanov | Theory of random sets[END_REF] received attention in the probability literature over the last decades. There exist several candidate notions to define the mean of a random closed set, see [START_REF] Molchanov | Theory of random sets[END_REF] (Chapter 2). We choose a rather intuitive one, based on the notion of coverage function: This definition has been applied in the Kriging framework to estimate sets of critical input values [START_REF] Chevalier | Estimating and quantifying uncertainties on level sets using the Vorobev expectation and deviation with Gaussian process models[END_REF][START_REF] Chevalier | Fast uncertainty reduction strategies relying on Gaussian process models[END_REF]. It uses the Vorob'ev expectation, based on the upper level sets

Q β = {z ∈ D, p Y (z) ≥ β}, called β-quantiles.
Definition 2 (Vorob'ev expectation and deviation). Denoting by µ the Lebesgue measure on R m and assuming that E(µ(Y)) < +∞, the Vorob'ev expectation is the β * -quantile

Q β * such that E(µ(Y)) = µ(Q β * )
if this equation has a solution, and otherwise it is defined from the condition µ(Q

β ) ≤ E(µ(Y)) ≤ µ(Q β * ), ∀β > β * . The associated Vorob'ev deviation is the quantity E(µ(Q β * ∆Y))
, where ∆ denotes the symmetric difference between sets, i.e.

Q β * ∆Y = (Q β * ∪ Y) \ (Q β * ∩ Y).
The Vorob'ev expectation is a global minimizer of the deviation among all deterministic closed sets with volume equal to the average volume of Y (see [START_REF] Molchanov | Theory of random sets[END_REF] for a proof): for any set M with µ(M

) = E (µ (Y)) , E (µ (Q β * ∆Y)) ≤ E (µ (M ∆Y)).

Quantification of uncertainty on Pareto fronts using random set theory

In the MOO literature, the study of distribution location and spread of an attained set X rely on the attainment function α X [START_REF] Grunert Da Fonseca | The attainment-function approach to stochastic multiobjective optimizer assessment and comparison[END_REF]: the probability for a given point in the objective space to be dominated by an RNP set, α X = P(x ∈ Y).

Attained sets are also closed1 and unbounded subsets in R m . Hence, the attained sets obtained with the simulated CPFs can be considered as realizations of a random closed set and are denoted by Y i , (i = 1, . . . , N ). Looking again at Definition 1, one can see that the attainment function is in fact a coverage function. For proofs about the equivalence of the distribution of an RNP set and the corresponding attained set as well as for a definition of the attainment function in terms of coverage function, we refer to [START_REF] Grunert Da Fonseca | The attainment-function approach to stochastic multiobjective optimizer assessment and comparison[END_REF]. This reference establishes the link between optimization results and random closed sets, in a case where the attainment function is computed from several runs of optimizers. Their comparison is then performed based on statistical hypothesis testing procedures.

In practice the attainment function is estimated by taking the proportion of RNP sets that dominates a given vector in the objective space: Definition 3 (Empirical attainment function). Given a sample of attained sets Y 1 , . . . , Y N , the empirical attainment function is defined as: αN :

R m → [0, 1], αN (z) = 1 N N i=1 1 {z∈Yi} where 1 {z∈Yi} = 1 if z ∈ Y i , 0 otherwise.
An example of an empirical attainment function is presented in Figure 4, showing where in the objective space there is a high probability to improve on the current Pareto front.

Definition 2 requires that Y is bounded for its Vorob'ev expectation to exist. Hence it is necessary to define a reference point R to bound the integration domain. Since the Lebesgue measure of an attained set with respect to the reference point is the hypervolume indicator of the corresponding RNP set, denoted by I H (., R), the choice of R has a similar influence (see e.g. [START_REF] Auger | Hypervolume-based multiobjective optimization: Theoretical foundations and practical implications[END_REF]). Unless there is previous knowledge about the range of the objectives, we choose R as the maximum of each objective reached by the conditional simulations.

The Pareto frontier of the Vorob'ev expectation provides us with an estimate of the Pareto front, as illustrated in Figure 4. The value of the Vorob'ev deviation gives an idea about the variability of the simulated CPF and can be monitored as observations are added, as will be shown in Section 4. The procedure to determine the value β * corresponding to the Vorob'ev expectation (Vorob'ev threshold) as well as of the Vorob'ev deviation is described in Algorithm 2.

Remark 1. The last equality in Algorithm 2 comes from the following: Ω 1 z∈A∆B µ(dz) = I H2 (A, B, R) + I H2 (B, A, R) where I H2 (A, B, R) is the binary hypervolume indicator, defined for instance in [START_REF] Zitzler | Performance assessment of multiobjective optimizers: An analysis and review[END_REF]: the volume dominated by A and not by B, i.e.

I H2 (A, B, R) = I H (A ∪ B, R) -I H (B, R).
The determination of the Vorob'ev threshold β * requires the volumes of β-quantiles:

µ(Q β ) = Ω 1 αN (z)≥β µ(dz)
They can be estimated by numerical quadrature, i.e. by computing the values of αN on a grid when there are few objectives or relying on Monte Carlo schemes. The other integrations are performed using hypervolume computation procedures and their complexity is related to the difficulty of computing the hypervolume in general. Hence measuring the uncertainty with the Vorob'ev deviation would be possible with any number of objectives but would require approximating integrals to keep the computations affordable. Note that since we take the objectives separately, simulating with more objectives simply requires computing the additional conditional simulations corresponding to those objectives.

Y (i) 1 , R 1 × • • • × min i∈(1,...,N ) Y (i) m , R m 4:
Determine the average volume of the attained sets Y i :

E(µ(Y)) ≈ 1 N N i=1 Ω 1 {z∈Yi} µ(dz) = 1 N N i=1 I H (Y i , R) 5:
Find the value of the Vorob'ev threshold β * by dichotomy: set a = 0, b = 1: Remark 2. [START_REF] Grunert Da Fonseca | The attainment-function approach to stochastic multiobjective optimizer assessment and comparison[END_REF] proposes the use of the Vorob'ev median (Q 0.5 ) if no compact set is chosen for integration. While removing the problem of fixing the reference point, no equivalent of the Vorob'ev deviation seems available in this case.

while b -a < ǫ do if µ Q a+b 2 < E(µ(Y)) then b = a+b 2 else a = a+b 2 end if end while, β * = a+b
E(µ(Q β * ∆Y)) ≈ 1 N N i=1 Ω 1 (z∈Qβ * ∆Yi) µ(dz) = 1 N N i=1 (2I H (Q β * ∪ Y i , R) -I H (Q β * , R) -I H (Y i , R))
From a practical point of view, it is also useful for visualization purpose with few objectives to display the superposition of all the symmetric differences by defining an analogue of the attainment function:

Definition 4 (Symmetric-deviation function). The function δ Y : z ∈ R m → P (z ∈ Q β * ∆Y) is called the symmetric-deviation function of Y.
δ Y is the coverage function of Q β * ∆Y. It is estimated with the empirical symmetric-deviation function:

δN (z) = 1 N N i=1 1 {z∈Q β * ∆Yi} .
Deviation of the simulations

Application

Two-dimensional bi-objective test problems

In this section, we illustrate the benefits of the proposed methodology for estimating Pareto fronts. We consider the following two variable, bi-objective optimization problems from the literature: (P1) The problem presented in [START_REF] Parr | Improvement criteria for constraint handling and multiobjective optimization[END_REF], which has a convex Pareto front:

f 1 (x) = b 2 - 5.1 4π 2 b 2 1 + 5 π b 1 -6 2 + 10 1 - 1 8π cos(b 1 ) + 1 f 2 (x) = -(10.5 -b 1 )(b 1 + 5.5)(b 2 + 0.5) - 1 30 b 2 - 5.1 4π 2 b 2 1 -6 2 - 1 3 1 - 1 8π cos(b 1 ) + 1 where b 1 = 15x 1 -5, b 2 = 15x 2 and x 1 , x 2 ∈ [0, 1].
(P2) The ZDT3 problem [START_REF] Zitzler | Comparison of multiobjective evolutionary algorithms: Empirical results[END_REF] which has a disconnected Pareto front.

For each example, we start with a set of few observations that allow fitting initial Gaussian process models for the two objective functions. Then we add new points sequentially by maximizing the Expected Hypervolume Improvement, based on the formula detailed in [START_REF] Emmerich | Hypervolume-based expected improvement: Monotonicity properties and exact computation[END_REF]. At each step, the Gaussian process models are updated and their hyperparameters re-estimated. These models are then used to simulate CPFs, from which we compute the estimates of the Vorob'ev mean and the measures of uncertainty: Vorob'ev deviation and symmetric-deviation function. Since the integration domain varies as points are added, the values are displayed divided by the volume of this integration domain. The following test problems are fast to compute, so it is possible to compare the outcome of the proposed workflow to a reference Pareto front by using the volume of the symmetric difference.

The results are presented in Figure 6 and Figure 7, showing the evolution of the estimated Pareto fronts with the corresponding uncertainty around it. For the problem (P1) the sequence is detailed, demonstrating the strength of the proposed approach for giving insights about the uncertainty on the Pareto front. In particular, the uncertainty measures are helpful for choosing a minimal number of observations for approximating the Pareto front: while 10 initial observations may not be enough (Figure 6a) regarding the large symmetricdeviation, adding 10 more observations dramatically reduces the uncertainty (Figure 6c).

The conclusions are similar for problem (P2) Figure 7, where the Pareto front is disconnected, starting this time with 20 observations and sequentially adding ten more observations by Expected Hypervolume Improvement maximization. This example makes clear that the Vorob'ev expectation refers to the area dominated by the Pareto front. When the latter is disconnected, the dominated area's frontier is horizontal in the corresponding parts. The position of the cuts in the Pareto front depends on the model and simulations, resulting in an higher variation around cuts: a small change in the extent of a peak impacts the beginning of the next one (where the symmetric difference volume depends on the size of the disconnection).

One can note that the approximations obtained are dependent on the model accuracy. In Figure 6a the approximation is clearly too optimistic, due to an underestimation of the range parameters of the used Matérn covariance kernel (ν = 5/2). Similarly, for one objective, the expected value of the minimum would be misleading at the beginning.

Finally, we propose the use of the Vorob'ev deviation as a stopping criterion when the Pareto front location is known. An empirical rule could be a threshold on the Vorob'ev deviation (e.g. expected volume of the symmetric difference less than 1% of the integration volume) and detection of stagnation (e.g. under the threshold for several evaluations). On the examples Figures 6d and7b, by considering the two last evaluations, the result would have been to stop for problem (P1) and continue for problem (P2). Note that this simple criterion would fail if the estimation of hyperparameters is misleading. A more robust version would be to integrate the uncertainty on the hyperparameters estimation in a full Bayesian framework [START_REF] Diggle | Model-based geostatistics[END_REF]. More sophisticated values could also be derived inspired from [START_REF] Wagner | A taxonomy of online stopping criteria for multi-objective evolutionary algorithms[END_REF].

Additional experiments on conditional simulations

The aforementioned methodology depends on the number and location of simulation points used to obtain the CPFs from the Gaussian process models (cf. Algorithm 1 and first step of Algorithm 2). As a first study, we have compared two sampling strategies: uniform sampling and space-filling sampling relying on a Sobol sequence, again with 2 variables. The objective functions are taken as sample paths of centered Gaussian processes with Matérn covariance kernel (ν = 5/2), with range parameters equal to 0.3/ √ 3 for f 1 , and 0.5/ √ 3 for f 2 .

We compute the approximation error over the set of non-dominated points obtained from the two sampling strategies. To compare the results with a reference set obtained with an NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] with archiving, three er- The figure description is the same as in Figure 6.

ror indicators are used: hypervolume difference, epsilon and R2 quality indicators [START_REF] Zitzler | Performance assessment of multiobjective optimizers: An analysis and review[END_REF]. The tests are repeated one hundred times. The results presented in Figure 8 show that space-filling sampling slightly outperforms uniform sampling to get an accurate estimation of the Pareto front. Additional tests performed showed a similar behavior with 3 variables but no difference for 10 variables, where the number of points considered was too small.

Here, only the impact of the error introduced by discretizing conditional simulations is studied. Given a few hundred to a few thousands of simulation points, it should remain negligible considering a relatively low number of variables, as is usually the case in application examples up to six variables and six objectives [START_REF] Svenson | Computer experiments: Multiobjective optimization and sensitivity analysis[END_REF]. Alternatively, resorting to approximate spectral simulation methods such as the truncated Karhunen-Loève expansion could be considered.

Conclusion and perspectives

We presented an original methodology to estimate and visualize the uncertainty of Pareto front approximations, based on Gaussian process conditional simulations. More precisely, the attainment function provides an estimation of the probability of dominating a given point in the objective domain. Then, a global uncertainty measure was defined relying on the theory of random sets through the concept of Vorob'ev deviation. It indicates the confidence of the model on the approximation of the attained set. The last tool is a visualization of the region of confidence for two or three objectives. Application on an higher number of objectives would be feasible, requiring the use of Monte Carlo methods for the computation of the various integrals. As illustrated on two bi-objective problems with convex or disconnected Pareto fronts, these measures can also be used as a basis to define stopping criteria in a sequential framework.

Further work is needed to analyze the different kinds of uncertainty and biases that may occur when applying the proposed methodology. In addition, techniques for simulating efficiently over more points and updating simulations with new observations [START_REF] Chevalier | Fast update of conditional simulation ensembles[END_REF] should be considered, as well as optimization of simulation points locations with re-interpolation or re-simulation [START_REF] Oakley | Bayesian uncertainty analysis for complex computer codes[END_REF]. Another direction for future research includes the integration of the proposed uncertainty estimate in a Stepwise Uncertainty Reduction (SUR) strategy [START_REF] Chevalier | Fast parallel kriging-based stepwise uncertainty reduction with application to the identification of an excursion set[END_REF] as an infill criterion.

Figure 1 :

 1 Figure 1: (a) Example of Kriging model based on observations at n = 7 locations (black points), with Kriging mean (bold line) and Kriging pointwise 95% prediction intervals (shaded area). The Gaussian predictive distribution for x = 0.7 is represented by the vertical dashed line. (b) Conditional simulations (colored thin lines) from the fitted UK metamodel.

2 Figure 2 :

 22 Figure 2: Conditional Pareto sets and fronts corresponding to the GP models Y 1 , Y 2 , based on the observations An represented by blue triangles. Left and center: examples of two conditional simulations of Pareto sets (top) and fronts (bottom), where the simulations are performed on a regular 100 × 100 grid. The simulation points and simulated responses are plotted with dots. The corresponding non-dominated points are represented by red squares. Right: contour plot of the probability density of optimal points in the decision space estimated from 30 conditional simulations (top) and superposition of simulated conditional Pareto fronts (bottom).

f 2 f 1 Figure 3 :

 13 Figure 3: Example of 3 realizations of RNP sets (points, triangles and squares) and the corresponding attained sets (shaded areas).

Definition 1 (

 1 Coverage function). Let Y be a random closed set on a topological space D (here D ⊂ R m equipped with the topology induced by the Euclidean distance). The coverage function p Y is defined by p Y : x ∈ D → P(x ∈ Y).

Algorithm 2 3 :

 23 General procedure for estimating the Vorob'ev expectation and deviation 1: Generate N CPFs (see Algorithm 1). 2: if R is unknown then find the extremal values for the objectives over the RNP sets realizations: Define the integration domain: Ω = min i∈(1,...,N )

2 6 :

 6 Estimate the Vorob'ev deviation:

Figure 4 :

 4 Figure 4: Example of empirical attainment function (level sets, left) and corresponding Vorob'ev expectation (shaded area, right), based on 200 conditional simulations. The estimate of the underlying Pareto front (dashed line) is the Pareto frontier of the Vorob'ev expectation. The ten observations are marked by blue triangles.

Figure 5

 5 Figure 5 presents an example of a symmetric difference between two sets and an empirical symmetricdeviation function. This shows the variability around the estimated Pareto front: dark areas indicate regions where the estimation of the Pareto front is not known precisely.

Figure 5 :

 5 Figure 5: Left: symmetric difference between the Vorob'ev expectation (the level line of the Vorob'ev threshold is represented by the red dashed line) and a simulated CPF's attained set (blue dotted line). Right: illustration of the deviation around the estimated Pareto front corresponding to Figure 4 with an example of empirical symmetric-deviation function (level plot).

  0.05 0.10 0.15 0.20 0.25 0.30 Number of observations Vorob'ev deviation Distance to the true Pareto front (d)

Figure 6 :

 6 Figure 6: Evolution of the deviation with new observations added using Expected Hypervolume Improvement for Problem (P1). The shaded area represents the image of E by f with a thicker border for the Pareto front. Observations are marked with blue triangles and the blue solid line represents the current Pareto front. The dashed line is the estimated Pareto front, with the corresponding values of the symmetric-deviation in level plot. Bottom right: evolution of the Vorob'ev deviation scaled by the current integration volume (black solid line with circles) and evolution of the distance to the real Pareto front measured with the volume of the symmetric difference with the estimation (red dotted line with crosses).

Figure 7 :

 7 Figure 7: Evolution of the deviation with new observations added using Expected Hypervolume Improvement for Problem (P2). The figure description is the same as in Figure 6.

Figure 8 :

 8 Figure8: Hypervolume difference, epsilon and R2 quality indicators between random reference Pareto fronts and approximations of them. A hundred conditional simulations for each of the two objectives considered are generated on a bi-dimensional grid with 225 points. Each of them is re-interpolated by Kriging and a reference Pareto front is obtained by applying an NSGA-II[START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] with archiving. This reference set is used to compare the quality of approximations of the Pareto front obtained by uniform sampling (red crosses) or with a Sobol sequence (blue points) with respect to the number of simulation points. The error bars indicate the quantiles of level 5% and 95% for the hundred repetitions.

  Algorithm 1 Simulation of N conditional Pareto sets and fronts for i = 1, 2, . . . , N do Choose p simulation points E p = e 1 , . . . , e p in E (fixed or different at each iteration). for j = 1, 2, . . . , m do Generate a conditional simulation at e 1 , . . . , e p for the j th objective:

As a finite union of closed sets (hyper quadrants).

Acknowledgments

This work has been conducted within the frame of the ReDice Consortium, gathering industrial (CEA, EDF, IFPEN, IRSN, Renault) and academic ( École des Mines de Saint-Étienne, INRIA, and the University of Bern) partners around advanced methods for Computer Experiments. The authors also thanks two anonymous reviewers for their helpful suggestions.


Lassoued Jamila
email: lassoued.jamila@gmail.com

Zemzemi Nejib
email: nejib.zemzemi@inria.frrapportde

Numerical Investigation in the usefulness of Proper Orthogonal Decomposition method in cardiac electrophysiology

Keywords: monodomain model, reduced order method, POD, mitchell schaeffer model, ionic parameters modèle monodomaine, méthode de reduction de modèle, POD, modèle de mitchell schaeffer, paramètres ionique

Numerical simulation of cardiac electrophysiology is very time consuming. Reduced order method have been recently used in different fields including cardiac electrophysiology. In this paper we use a reduced order method based on the proper orthogonal decomposition (POD), and we propose to evaluate the accuracy of this method while changing different parameters in the model. To describe the propagation of the action potential in the myocardium, we use the monodomain model which is a reaction diffusion PDE system coupled to a dynamic system of ODEs representing the time evolution of the electrophysiology in the cell membrane. We build the reduced order model using a set of parameters, afterwards, we evaluate the accuracy of the reduced model while changing the parameters of the model. We numerically analyze the sensitivity of the reduced order method to the model parameters including the change of the whole ionic model.

Introduction

In order to make the numerical simulation of electric activity of the heart useable in clinical applications, we need to develop a model that takes into account the physiological properties of the heart and provides realistic propagation of the electrical wave and especially not too expensive in terms of computation cost.

The electrical wave in the heart is governed by a system of reaction diffusion equations called bidomain model, it is coupled to an ODE system representative the cellular activity.

This model, which consists of a reaction diffusion non linear system, is coupled to dynamic system modelling the cellular ionic currents. It is alsoknown to be very expensive from a computation point of view [START_REF] Scacchi | Multilevel schwarz and multigrid preconditioners for the bidomain system[END_REF][START_REF] Franzone | A parallel solver for reaction-diffusion systems in computational electrocardiology[END_REF]. In this work, we choose the monodomain model [START_REF] Clements | Activation dynamics in anisotropic cardiac tissue via decoupling[END_REF][START_REF] Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models[END_REF], a simplified version of the bidomain model, we reduce the complexity of its discretization using a Galerkin basis that apriori contains sufficient informations about the expected solution. There are different methods to construct this basis. Following, [START_REF] Boulakia | Reduced-order modeling for cardiac electrophysiology, application to parameter identification[END_REF][START_REF] Boulakia | Numerical simulation of electrocardiograms[END_REF], we choose to use the Proper Orthogonal Decomposition (POD) method. This approach has been used in many fields of science and engineering, including in parameter estimation for cardiac electrophysiology [START_REF] Boulakia | Reduced-order modeling for cardiac electrophysiology, application to parameter identification[END_REF]. But to the best of our knowledge, no work has been done in analyzing the accuracy of this method with respect to changing the parameters of the model. In this work we propose to address the problem of the sensitivity of POD with respect to parameter changing.

The POD is a method to derive reduced models by projecting the system onto subspaces spanned by a basis of elements that contains the main features of the expected solution. This method essentially provides an orthonormal basis for representing the given data in a certain least squares optimal sense. Truncation of the optimal basis provides a way to find optimal lower dimensional approximations of the given data. To be able to use POD in parameter estimation like in [START_REF] Boulakia | Reduced-order modeling for cardiac electrophysiology, application to parameter identification[END_REF], we think that a prior study should tackle the sensitivity of the POD solution with respect to the parameters that we want to estimate. Of course, we have a rigorous quantification of the error when the POD method is used to solve the problem with the same parameter with which the POD basis has been constructed. But we don't have any a priori estimation of the error while changing parameters, especially the parameters of the ionic model.

The numerical study that we will address in this work compares the POD solution to the full finite element solution of the monodomain problem when we change the parameters of the ionic model. Moreover, we also check the sensitivity of the method while changing the whole ionic model. In section 2, we will present the mathematical models and the reduced-order method we will be using in this work. In section 3 we will carry out numerical simulation studying the sensitivity of the POD method. The main conclusions of the study are then summarized in section 4.

Modelling and numerical methods

Electric model

The bidomain model is based on the assumption that, at the cell scale, the cardiac tissue can be viewed as partitioned into two ohmic conducting media, separated by the cell membrane: intracellular, made of the cardiac cells, and extracellular which represents the space between them (see [START_REF] Neu | Homogenization of syncytial tissues a comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart[END_REF][START_REF] Pennacchio | Multiscale modeling for the bioelectric activity of the heart[END_REF]13,[START_REF] Sundnes | Computing the electrical activity in the heart[END_REF]). After an homogenization process (see [START_REF] Pennacchio | Multiscale modeling for the bioelectric activity of the heart[END_REF]), the intra-and extracellular domains can be supposed to occupy the whole heart volume Ω H . The heart is assumed to be isolated, the propagation of the electric wave in the myocardium is then described by the following system of equations:

                   χ m ∂V m ∂t + I ion (V m , w) -div(σ i ∇V m ) -div(σ i ∇u e ) = I app in Ω H × (0, T ) -div((σ i + σ e )∇u e ) -div(σ i ∇V m ) = 0 in Ω H × (0, T ) ∂w ∂t + G(V m , w) = 0 in Ω H × (0, T ) σ i ∇V m .n = 0 on Σ. σ e ∇u e .n = 0 on Σ, (1)
where Ω H and Σ denote respectively the heart domain and its boundary. The time domain is given by [0, T ]. And χ m the membrane capacitance per area unit. The variables V m and u e denote the action potential and the extracellular potential, σ i and σ e the intra-and extracellular conductivity tensors. The term I app is a given source current, used in particular to initiate the activation, and I ion (V m , w) represents the ionic current across the membrane. In this study we use the monodomain model [START_REF] Potse | A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart[END_REF], where we assume that the ratios of the anisotropy of the electrical conductivity tensor are the same in the intra and extracellular media. Following [START_REF] Clements | Activation dynamics in anisotropic cardiac tissue via decoupling[END_REF], the reaction-diffusion equation governing the action potential reads as follow:

                   χ m ∂V m ∂t + I ion (V m , w) -div(σ∇V m ) = I app in Ω H × (0, T ) ∂w ∂t + G(V m , w) = 0 in Ω H × (0, T ) σ∇V m .n = 0 on Σ, (2)
where

σ def = σ i (σ i + σ e ) -1
σ e is the bulk conductivity.

Numerical methods

In order to solve numerically the monodomain model, we use finite element method. After discretizing (2) in time, we start by writing the variational formulation of the problem

         w n+1 = w n -δtG(V n m , w n+1 ) in Ω H χ m Ω H V n+1 m ψ + δt Ω H σ∇V n+1 m ∇ψ = χ m Ω H V n m ψ +δt Ω H (I n+1 app -I ion (V n m , w n+1 ))ψ ∀ ψ ∈ H n (Ω H ) (3)
Projecting the variational formulation on the discretized finite element space, we obtain the following linear system

     w n+1 = w n -δtG(V n m , w n+1 ) in Ω H χ m M V n+1 m + δtSV n+1 m = χ m M V n m + δtM (I n+1 app -I ion (V n m , w n+1 )) (4)
where M = (

Ω H e i e j ) i,j=1,••• ,n , S = ( Ω H σ∇e i ∇e j ) i,j=1,••• ,n
and (e i ) 1≤i≤n is the Galerkin finite element basis. The matrix M is the mass matrix and S is the stiffness matrix. These functions I ion (V m , w), G(V m , w) depends on the used ionic model. In this study, the dynamics of w and I ion are described by different ionic models including a two-variables model phenomenological by Mitchell and Schaeffer [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF], a four-variables model by Bueno [START_REF] Bueno-Orovio | Minimal model for human ventricular action potentials in tissue[END_REF] and a complete physiological cell membrane ionic model by TenTusscher [START_REF] Tusscher | Cell model for efficient simulation of wave propagation in human ventricular tissue under normal and pathological conditions[END_REF].

Reduced order method

For the sake of completeness, we briefly recall the principle of the proper orthogonal decomposition method. The POD is a linear process that includes determining an optimal orthogonal basis of eigen modes in the sense of energy. This means that no base is also capable of capturing a higher amount of energy with the same number of modes. There are two phases in this method, the first is the generation of a reduced order basis and the second is solving the reduced order problem.

Generation of the reduced order basis

The problem to be solved can be stated as follows: how to approximate a function u(x, t) as follows

u(x, t) ≃ K k=1 a k (t)Φ k (x) (5)
where (Φ k ) K k=1 is a new basis and (a k (t)) K k=1 is a temporal coefficients. Therefore, we need to find the new basis (Φ k ) K k=1 and the temporal coefficients (a k (t)) K k=1 . Assume that we know the values of the function u in m locations

x 1 , x 2 , • • • , x m in different p moments.
Let's denote A, the snapshot matrix containing the values of the function u, where each colonne correspond to a given time:

A =          u(x 1 , t 1 ) u(x 1 , t 2 ) . . u(x 1 , t p ) u(x 2 , t 1 ) u(x 2 , t 2 ) . . u(x 2 , t p ) . . . . . . . . . . . . u(x m , t 1 ) u(x m , t 2 ) . . u(x m , t p )          ∈ R m×p
Solving the problem of approximation ( 5) is then equivalent to determining the orthonormal basis Φ k (x) K k=1 with K ≤ p solution of the problem min

(Φ k ) K k=1 ∈mathbbR m p i=1 u(x, t i ) - K k=1 (u(x, t i ), Φ k (x))Φ k (x) 2 . ( 6
)
The solution of the minimization problem ( 6) is given by the truncated Singular Value Decomposition (SVD) in the order K of the snapshots matrix A. We refer the reader interested in more details to [START_REF] Rathinam | A new look at proper orthogonal decomposition[END_REF][START_REF] Volkwein | Model reduction using proper orthogonal decomposition[END_REF] for a detailed explanation of the link between SVD and POD. Using SVD we obtain,

A = U ΣV T
where U and V correspond respectively to the the left singular vectors of A of dimension m × m and the right singular vectors of A of dimension p × p. and Σ is a diagonal matrix containing the σ 1 , σ 2 , • • • , σ r , called the singular values of A (and A T ) such as

σ 1 σ 1 , • • • , σ r with r = min(m, p).
The solution of ( 6) is given by the left singular vectors of A truncated to order K: we find the matrix Φ as the K first columns of the matrix U

Φ = U K ∈ R m×K (7)
The accuracy of the POD method is related to the following energy quantification [START_REF] Volkwein | Model reduction using proper orthogonal decomposition[END_REF]:

p j=1 A j - K i=1 (A T i Φ i ) Φ i 2 = d i=K+1 λ i
where A j is the jth colonne of the matrix A and d the rank.

The number K of eigen modes is choosen such that d i=K+1 λ i is less than some tolerance.

Solving the reduced order problem

The POD basis is used by applying a projection procedure, Instead of being projected on the full finite element basis, the variational formulation (3) would now be projected on the POD basis Φ = (φ k ) K k=1 . Knowing the reduced order solution V n pod at time t n , we look for

V n+1 pod = K k=1 v n+1 k φ k . The unknown vector v n+1 = (v n+1 k ) K k=1
is obtained by solving the following reduced order problem,

Φ T (χ m M + δtS)Φv n+1 = Φ T (χ m M V n pod + δtM (I n+1 app -I ion (V n pod , w n+1 ))) (8)
The left hand side matrix Φ T (χ m M + δtS)Φ is the reduced order matrix and could be easily inverted once for all in our problem. The right hand side is the finite element right hand-side vector projected on the POD basis.

For both equations, the finite element system (4) and the reduced order (8), we use the LU method to determine the full finite element and the POD solutions. In the first case, the order of the matrix χ m M +δtS is 1681, then it requires 3.166710 9 floating point operations, while the second requires only 3.376910 4 . This improvement in the number of operations will be illustrated in the following subsection 2.4, when we will present a comparison of the solutions obtained.

Numerical comparison between the full finite element and the POD solutions

In this section, we present numerical simulations about the POD approximation of the monodomain model. For the sake of simplicity and reproducibility, the heart domain is the square Ω H = [0, 1]×[0, 1] and the unit is cm. It is discretized on the x and y-axis with a space step 0.025 cm. The time domain is [0, 500], its unit is ms and the time step is 0.5 ms. The heart is stimulated in a region at the left bottom corner [0, 0.025] × [0, 0.025] during 2 ms. The parameter of the ionic model are given in Table 1. As mentioned in the previous section we start by generating a solution of the monodomain problem. This solution is stored in the snapshot matrix used to construct the POD basis. In this example we have 1681 degrees of freedom. Capturing 99.99 % the finite element solution energy requires only 37 modes. We use these 37 modes as a basis and we solve the POD problem.

τ
Using a matlab code, the full finite element solution costs 10.34 s, where the ODE system costs 0.31 s and the linear system costs 10.03 s. The solution of the reduced order problem is computed in 0.5 s, where 0.31 for the ODE system and 0.19 s for the linear system. In Figure 1, we show three snapshots of the full finite element solution (top) and the POD solution (bottom) at the depolarization phase, and in the Figure 2 at the repolarization phase. The pattern of the solution is the same and the wave front is accurately captured.

Numerical results

For a given set of parameters and using the standard finite element method, we can compute the action potential and then construct the snapshots matrix A. Using SVD method described above we obtain a reduced order basis that is able to generate the solution up to the considered tolerence. The problem comes when we change some parameters and we try to solve the monodomain problem using the same reduced order basis. Our concern is how the solution of the monodomain problem will behave when using the POD method? What is the accuracy of the method? To answer to these question, we propose to change some parameters in the ionic model and assess the error between the POD and the full finite element solutions. In the paragraph 3.1, we propose to perform this assessment on a phenomenological ionic model with 5 parameters [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF] and on a physiological model with a large number of parameters [START_REF] Tusscher | Cell model for efficient simulation of wave propagation in human ventricular tissue under normal and pathological conditions[END_REF]. In paragraph 3.3, we propose to vary the whole ionic model and test if the POD method with a basis computed using Micheall Shaeffer model is able reproduce the solution of Ten Tusscher model.

Sensitivity of the POD solution to Micheall Shaeffer model parameters

In this paragraph, we use the Micheall and Shaeffer model to describe the electrical activity of the cell membrane. This model is used to describe the phenomenon of excitability while keeping a low complexity, and with only two state variables: the action potential V m and a recovery variable w. The dynamics of w and I ion are given by:

I ion (V m , w) = w τ in V 2 m (V m -1) - V m τ out G(V m , w) =            w -1 τ open if V m v gate w τ close if V m >v gate (9)
The time constants τ in , τ out are respectively related to the length of the depolarization and repolarization periods, τ open and τ close are the characteristic times of gate opening and closing respectively and v gate corresponds to the change-over voltage. Using the method of constant variation, we can easily determine the analytical expression of the solution in the interval [t n , t n+1 ]. For a given time t ∈ [t n , t n+1 ] we have

• if V m (t n ) v gate , w(t n+1 ) = 1 -(1 -w(t n ))e (t-tn)/τopen . • if V m (t n )>v gate , w(t n+1 ) = w(t n )e (t-tn)/τopen .
The original values of these parameters are given in Table 1 ( [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF]).

RR n°7003

In order to test the accuracy of the POD solution, we first compute the POD basis using the original parameters of the model. Then we propose to perform for each parameter some simulation, where the considered parameter vary between half and three halfs its original value given in Table 1. In table 2, we provide the different considered values for each parameter. We evaluate the error between the reduced order solution and the full finite element solution. All of these simulations are performed in the time domain T = 500 ms with a discretization step dt = 0.5 ms. In Figure 3, we show the evolution of the L 2 relative errors with respect to parameters τ close (top left), τ open (top right), τ out (bottom left) and τ in (bottom right). For parameters τ close , τ open and τ out the relative error is less than 1%. However, for parameter τ in the error significantly increases when the parameter is far from its original value. In particular for τ in less than 0.175 the error is higher than 1% and it reaches 10% for τ in = 0.15. This means that POD basis constructed with the original parameters is able to approximate the solution with a good accuracy for different values of parameters τ close , τ open and τ out . But for τ in , the accuracy is acceptable (less than 1%) when τ in ∈ [0.175, 0.45]. In Figure 4, we show the time course of the L 2 relative error between the finite element solution to a POD solution computed for τ in = 0, 45. The magnitude of the error is more important in the depolarization phase. In the plateau phase the error is very low but is increases with repolarization phase.

τ

Physiological model

In this paragraph we perform a same study as for Micheall-Shaeffer model but with physiological model Ten Tusscher [START_REF] Tusscher | Cell model for efficient simulation of wave propagation in human ventricular tissue under normal and pathological conditions[END_REF]. The Ten Tusscher Model represents a physiological model of cardiac cell human ventricular. It is established from experimental measurements on the most major ionic currents, and it includes a dynamic basis for calcium. it shows that the differences in the morphologies can be explained by certain characteristics of two currents, potassium current transient outward and slow delayed rectifier, which differ from one type to another cell. The model is described by a set of differential equations (see [START_REF] Ten Tusscher | A model for human ventricular tissue[END_REF]). In particular we are interested in evaluating the POD solution to these parameters [START_REF] Tusscher | Cell model for efficient simulation of wave propagation in human ventricular tissue under normal and pathological conditions[END_REF]:

• g N a : the maximal Fast sodium (N a + ) current I N a conductance (the original parameter is 14.838).

• g kr : the maximal rapid delayed rectifier potassium current I Kr conductance (the original parameter is 0.153).

• g to : the maximal transient outward current I to conductance (the original parameter is 0.294).

We compute the space and time L 2 relative error between the full finite element and the POD solutions for the parameters given in Table 3. In Figure 5 we see that for g Kr and g to , the solution is suffitiently accurate. However, the parameter g N a the error is low small when the parameter is close to the original value but it dramatticalt increases when the parameter is so far.

g N

Sensitivity of the POD solution to different models

In this paragraph, we will test the capability of a POD basis constructed using a phenological model (Micheall-Shaeffer) to capture the propagation of the electrical wave governed by a physiological ionic model. Firstly, we generate a finite element solution using the Micheall-Shaeffer model, from this solution we construct a POD basis. After that, we solve the reduced order problem of the monodomain equation coupled to the Ten Tusscher ionic model. We then compute the L 2 relative error between this solution and the full finite element solution of mondomain coupled to Ten Tusscher model. In Figure 6 (left), we plot the time course of this error. We observe that the error is high, mainly at the depolarisation phase where it reaches 10%. We performed the same test for the Bueno et al. model [START_REF] Bueno-Orovio | Minimal model for human ventricular action potentials in tissue[END_REF][START_REF] Bueno-Orovio | Mathematical modeling and spectral simulation of genetic diseases in the human heart[END_REF]. In Figure 6 (right), we remark that the error is about 1% for this case. These results show that the scale of the action potential values plays an important role. In fact, the action potential values in the Micheall-Shaeffer model are between 0 mV and 1 mV, whereas they are between -86 mV and 25 mV in the Ten Tusscher model and between 0 mV and 1.4 mV in the Bueno et al. model. For this reason we propose to rescale the Micheall-Shaeffer action potential to the new model before generating the POD basis. In Figure 7 (left (respectively right)), we plot the relative error between the reduced order solution using the new POD basis constructed from data rescaled to Ten Tusscher (respectively, Bueno) model and the full finite element solution. We see in both plots that changing the scale before generating the POD basis improves the accuracy of the solution.

and conclusions

We have presented in this work a reduced order approach based on POD method for the computation of the electrical activity of the heart. This approach not expensive in terms of computational cost and is very accurate when we reproduce the solution on which we have build the POD basis. In order to evaluate the usefulness of this approach in parameter estimation problem, we build a POD basis using the original parameters of the ionic model and we computed the L 2 relative error between the finite elements solution and the reduced order solution for different parameters. For M-S ionic model we conclude that we obtain a good accuracy when we vary the values of τ open , τ close and τ out . Whereas for the parameter τ in , we have seen that there is a trust interval, out of it the solution is not sufficiently accurate. For the physiological Ten Tusscher and Panfilov model we have seen that the error between POD solution and the full finite elements solution is acceptable for parameters g kr and g to and important for parameter g N a out of a trust interval. Both parameters τ in for M-S and g N a for the Ten Tusscher and Panfilov model play an important role in the depolarization phase and the wave front velocity. This explains the sensitivity of the POD solution to these parameters. In order to obtain a good accuracy out of the trust interval, one should update the POD basis.

Since the propagation of the electrical wave has the same pattern if we only change the ionic model, we have tested the case when we use a basis generated from M-S model to solve the Monodomain problem coupled to other ionic models (Ten Tusscher and Pafilov, Bueno). We observed that the POD solution is computed with a good accuracy for the Bueno model. On the contrary, the error is very important for the Ten Tusscher and Panfilov model. We found that this unaccuracy is due to the difference in scale between the transmembrane potential range in different models. Rescaling the transmembrane potential before generating the POD basis comes with an improve of the POD solution mainly for the Ten Tusscher and Panfilov model where the range of the action potential is considerably higher than for M-S model. In future works, we aim to use the POD approach in the parameters estimation problem in cardiac electrophysiology. In particular, we will use physiological optical mapping measurements in order to personalize the electrical model on wedge preparation of animal tissue.

TimeFigure 1 :

 1 Figure 1: Top (respectively, bottom):Snapshots of the full finite element (respectively POD) solution at times 10, 30 and 40 ms (from left to right) at the depolarization phase

Figure 2 :

 2 Figure 2: Top (respectively, bottom):Snapshots of the full finite element (respectively POD) solution at times 100, 200 and 300 ms (from left to right) at the repolarization phase

Figure 3 :

 3 Figure 3: The error between the finite element solution and the POD solution with respect to τ close (top left), τ open (top right), τ out (bottom left) and τ in (bottom right)In Figure4, we show the time course of the L 2 relative error between the finite element solution to a POD solution computed for τ in = 0, 45. The magnitude of the error is more important in the depolarization phase. In the plateau phase the error is very low but is increases with repolarization phase.

Figure 4 :

 4 Figure 4: The time course of the spacial L 2 relative error between the solution to a POD solution computed for τ in = 0, 45.

Figure 5 :

 5 Figure 5: Left (respectively right, bottom) the error between the finite elements solution and the POD solution with respect to g N a (respectively g Kr , g to )

Figure 6 :

 6 Figure 6: Left (respectively, right), the error between the Micheall Shaeffer full finite element solution and the POD solution of the monodomain equation coupled to the Ten Tusscher (respectively, Bueno) ionic modelFor this reason we propose to rescale the Micheall-Shaeffer action potential to the new model before generating the POD basis. In Figure7(left (respectively right)), we plot the relative error between the reduced order solution using the new POD basis constructed from data rescaled to Ten Tusscher (respectively, Bueno) model and the full finite element solution. We see in both plots that changing the scale before generating the POD basis improves the accuracy of the solution.

7 :

 7 (respectively, right), the error between the full finite element solution and the POD solution with scaleted Ten Tusscher (respectively, Bueno) POD basis

Table 1 :

 1 Original cell membrane parameters.

	close	τ open	τ out	τ in	V gate
	150	120	6	0.3	0.13

Table 2 :

 2 Parameters used for POD simulations.

	close	τ open	τ out	τ in
	75	60	3	0.15
	100	80	4	0.2
	125	100	5	0.25
	150	120	6	0.3
	175	140	7	0.35
	200	160	8	0.4
	225	180	9	0.45

Table 3 :

 3 Cell membrane parameters.

	a	g Kr	g to
	7.4190	0.0765	0.1470
	9.8920	0.1020	0.1960
	12.3650	0.1275	0.2450
	14.8380	0.1530	0.2940
	17.3110	0.1785	0.3430
	19.7840	0.2040	0.3920
	22.2570	0.2295	0.4410

RR n°7003

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399