Alessandro Barone
email: alebarone23@gmail.com

Michele Giuliano Carlino

Alessio Gizzi

Simona Perotto

Alessandro Veneziani

Efficient Estimation of Cardiac Conductivities: a Proper Generalized Decomposition Approach

Keywords: Computational Electrophysiology, Model Order Reduction, Data Assimilation, Proper Generalized Decomposition, Parameter Identification

While the potential groundbreaking role of mathematical modeling in electrophysiology has been demonstrated for therapies like cardiac resynchronization or catheter ablation, its extensive use in clinics is prevented by the need of an accurate customized conductivity identification. Data assimilation techniques are, in general, used to identify parameters that cannot be measured directly, especially in patient-specific settings. Yet, they may be computationally demanding. This conflicts with the clinical timelines and volumes of patients to analyze. In this paper, we adopt a model reduction technique, developed by F. Chinesta and his collaborators in the last 15 years, called Proper Generalized Decomposition (PGD), to accelerate the estimation of the cardiac conductivities required in the modeling of the cardiac electrical dynamics. Specifically, we resort to the Monodomain Inverse Conductivity Problem (MICP) deeply investigated in the literature in the last five years. We provide a significant proof of concept that PGD is a breakthrough in solving the MICP within reasonable timelines. As PGD relies on the offline/online paradigm and does not need any preliminary knowledge of the high-fidelity solution, we show that the PGD online phase estimates the conductivities in real-time for both two-dimensional and three-dimensional cases, including a patient-specific ventricle.

Introduction

The use of mathematical computational models in medicine is a consolidated approach (sometimes called in silico working, alongside the traditional in vivo and in vitro analyses) to have a deep understanding of physio-pathological dynamics. However, the complete translation of these models to clinical practice is prevented by several factors, including the need for their accurate customization to reflect the specific patient's condition [START_REF] Winslow | Computational medicine: translating models to clinical care[END_REF]. This requires accurate image processing procedures for the morphology retrieval as well as the personalized identification of model parameters. The parameters often cannot be measured and their fine tuning can be obtained only after data assimilation procedures (see, e.g., [START_REF] Bertagna | Data assimilation in cardiovascular fluid-structure interaction problems: an introduction[END_REF][START_REF] Bertoglio | Sequential parameter estimation for fluid-structure problems: Application to hemodynamics[END_REF][START_REF] D'elia | Applications of variational data assimilation in computational hemodynamics[END_REF][START_REF] Veneziani | Inverse problems in cardiovascular mathematics: toward patient-specific data assimilation and optimization[END_REF]), i.e., specific methodologies to combine available measures and numerical solutions for an enhanced modeling.

Different approaches can be used to include a measurable quantity into the mathematical modeling process.

In sequential stochastic procedures, the parameters are included as unknowns (or state variables) of the problem and are subject to an estimation procedure generally aimed to minimize the variance of the estimate or to maximize a probability density function. In variational approaches, the parameters to identify are used as control variables in the minimization of the mismatch between the results of the numerical model and the observed dynamics. In this way, the knowledge of observable quantities is converted into an estimate of the non-observable ones.

Data assimilation strategies [START_REF] Law | Data assimilation[END_REF][START_REF] Asch | Data assimilation: methods, algorithms, and applications[END_REF] are generally very powerful and able of an accurate personalization. However, they have intrinsic limitations and, in general, they may entail significant computational costs. Thus, as clinical problems often require relatively short timelines, empirical approaches are generally preferred to privilege efficiency over accuracy. In this way, the core knowledge of physiological and pathological dynamics expressed by physic-based models is discarded. This may be detrimental for patients whose features do not fall into the set of data used for the empirical methods.

In order to improve the efficiency of rigorous identification procedures, model reduction techniques can be used to replace high-fidelity models with educated surrogates, yet including the core knowledge of the problem at hand. Model reduction is an important topic developed over the years in the engineering and mathematical literature as the increased computational power was not able to cover the even more progressively rising complexity of problems addressed by mathematical modeling (see, e.g., [START_REF] Amsallem | Nonlinear model order reduction based on local reduced-order bases[END_REF][START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Quarteroni | Reduced basis methods for partial differential equations: an introduction[END_REF][START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF][START_REF] Luo | Proper orthogonal decomposition methods for partial differential equations[END_REF][START_REF] Costabal | Multi-fidelity classification using gaussian processes: Accelerating the prediction of large-scale computational models[END_REF]). Life sciences and medicine certainly fall in this picture.

In this work, we address a data assimilation problem of cardiac electrophysiology. The knowledge of the differential equations for the potential propagation in the cardiac tissue is quite consolidated, as witnessed by the specific literature (see, e.g., [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF][START_REF] Glass | Theory of heart: biomechanics, biophysics, and nonlinear dynamics of cardiac function[END_REF][START_REF] Pullan | Mathematically modelling the electrical activity of the heart: from cell to body surface and back again[END_REF]). Modeling improvements are mainly devoted to the micro-, meso-and macroscopic (i.e., cell, tissue and organ levels) description of the ions dynamics at cellular and subcellular level [START_REF] Fenton | Models of cardiac cell[END_REF], to their behavior at the cell-cell interface [START_REF] Lenarda | A modeling framework for electro-mechanical interaction between excitable deformable cells[END_REF][START_REF] Weinberg | Ephaptic coupling rescues conduction failure in weakly coupled cardiac tissue with voltagegated gap junctions[END_REF], and to the spatio-temporal coupling among the different cardiac components resulting in synchronized emerging phenomena [START_REF] Loppini | Spatiotemporal correlation uncovers characteristic lengths in cardiac tissue[END_REF]. These models feature parameters that are quite hard to measure in vivo and data assimilation procedures have been recognized as a viable approach [START_REF] Graham | Estimation of the bidomain conductivity parameters of cardiac tissue from extracellular potential distributions initiated by point stimulation[END_REF][START_REF] Abidi | Ionic parameters estimation in multi-scale cardiac electrophysiology modelling[END_REF]. In particular, variational techniques for the estimation of the tensor of the cardiac conductivities of the classical Bidomain and Monodomain models have been addressed and analyzed in [START_REF] Yang | Estimation of cardiac conductivities in ventricular tissue by a variational approach[END_REF]. The reliability of this approach for solving the Bidomain/Monodomain Inverse Conductivity Problem has been demonstrated on several benchmarks, covering synthetic as well as in vitro cases [START_REF] Barone | Numerical sensitivity analysis of a variational data assimilation procedure for cardiac conductivities[END_REF][START_REF] Barone | Experimental validation of a variational data assimilation procedure for estimating space-dependent cardiac conductivities[END_REF].

However, the efficiency of such procedures needs to be properly addressed, as the computational cost of the iterative mismatch minimization is generally high, especially when dealing with real geometries. This problem has been promptly recognized as a bottleneck, and several Reduced-Order Models (ROMs) have been investigated [START_REF] Boulakia | Reduced-order modeling for cardiac electrophysiology. application to parameter identification[END_REF][START_REF] Gerbeau | Reduced order model in cardiac electrophysiology with approximated lax pairs[END_REF][START_REF] Pagani | Numerical approximation of parametrized problems in cardiac electrophysiology by a local reduced basis method[END_REF][START_REF] Yang | Efficient estimation of cardiac conductivities via pod-deim model order reduction[END_REF][START_REF] Costabal | Multi-fidelity classification using gaussian processes: Accelerating the prediction of large-scale computational models[END_REF]. The proposed approaches rely on the construction of a surrogate, as a combination of basis functions generally built moving from previous solutions (called "snapshots") for a predetermined set of values for the parameters. For instance, in the Proper Orthogonal Decomposition (POD) considered in [START_REF] Yang | Efficient estimation of cardiac conductivities via pod-deim model order reduction[END_REF], the snapshots are smartly selected based on the concept of domain of influence in the space of the conductivities, and the surrogate is constructed after a Singular Value Decomposition (SVD) of the snapshot matrix.

In this paper, we explore a ROM procedure that does not require any a priori knowledge of the solution, even though it still relies on an offline/online paradigm. The Proper Generalized Decomposition (PGD) is a model reduction technique introduced by F. Chinesta et al. [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations: a Primer[END_REF], specifically devised to efficiently evaluate a parametrized differential problem when varying the parameter values. The idea is to treat the problem parameters as additional independent variables and to compute the solution over an extended domain, inclusive of the range the parameters are expected to belong to. As the actual numerical computation is performed in a highly dimensional domain during the offline stage, a special representation is introduced, where the surrogate solution is factorized with respect to each independent variable (including the parameters) or low-dimensional groups of variables. By advocating an iterative approach to compute this factorization, PGD proved to be an effective technique to reliably approximate several parametric problems [START_REF] Ammar | Reduction of the chemical master equation for gene regulatory networks using proper generalized decompositions[END_REF][START_REF] Niroomandi | Real-time simulation of biological soft tissues: a PGD approach[END_REF][START_REF] Signorini | Proper generalized decomposition solution of the parameterized Helmholtz problem: application to inverse geophysical problems[END_REF].

Actually, in the online phase, the solution is ready to be promptly evaluated for any value of the parameters as well as of the independent variables. This is particularly effective in a variational parameter identification procedure, when the solution for different guesses of the parameter(s) is tested in an iterative minimization process.

Based on these general properties, in this paper we explore the use of PGD to solve the Monodomain Inverse Conductivity Problem (MICP). After a short introduction to PGD in Section 2, we present the Monodomain model, the MICP and its basic features in Section 3. Successively, we introduce the specific technical aspects of the application of PGD to the Monodomain problem in Section 4. We detail the PGD approximation to the Monodomain problem in Section 4.1 and address implementation details in Section 4.2.

In Section 4.3 we assess the accuracy of the PGD model over a number of test cases. We start with a twodimensional (2D) test problem, yet in a morphology based on a previous study on canine tissues. Then, we extend the procedure to three-dimensional (3D) problems, including a patient-specific left ventricle. Finally, we apply PGD to the solution of MICP in Section 5. All the results pinpoint that, in spite of an offline stage that may be quite demanding in terms of CPU time, the PGD dramatically reduces the parameter identification to nearly real-time computations on common computational architectures. These preliminary results encourage further work for the extension of the methodology to the clinical practice (Section 6).

The PGD in a Nutshell

For the sake of completeness, we provide a short introduction to the PGD, to recall the basic ideas and provide some notation and terminology. For a more complete introduction, the reader is referred to [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations: a Primer[END_REF][START_REF] Ladevèze | The LATIN multiscale computational method and the proper generalized decomposition[END_REF][START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the proper generalized decomposition[END_REF].

In what follows, Ω ⊂ R d is the physical space domain, with d = 2, 3. We denote by V a Hilbert space, generally a subspace of the Sobolev space H 1 (Ω) [START_REF] Lions | Non Homogeneous Boundary Value Problems and Applications[END_REF]. With q ∈ R δ we denote a vector of parameters the problem we consider depends on. Let us focus on a generic parametrized elliptic problem in the weak form:

find u ∈ V s.t. a q (u, v) = F q (v) ∀v ∈ V, (1)
where, for any given admissible value of the parameters in q, a q : V × V → R and F q : V → R are a bilinear, continuous and coercive and a linear continuous form, respectively. More precisely, we assume that the parameter q i , for i = 1, . . . δ, ranges in an admissible interval S i , so that the admissible set for q is the cuboid S ≡ S 1 ×S 2 . . .×S δ . According to a PGD procedure, we regard the solution u as a function of x ∈ R d and q ∈ R δ . Specifically, we assume that u ∈ W =V ⊗ δ l=1 L 2 (S l ). Thus, the extended weak formulation, including the parameters as independent variables, is obtained by integrating (1) over S and coincides with the (d + δ)-dimensional problem: find u ∈ W s.t.

A(u, v) = F(v), ∀v ∈ W, (2)
with

A : W × W → R and F : W → R given by A(u, v) ≡ S a q (u, v)dq, F(v) ≡ S F q (v)dq, (3)
respectively. With a little abuse of notation, we keep denoting the unknown and the test functions with the same symbols, u, v, although their dependence on the parameters in q.

As the dimension of the problem is now increased, the numerical solution (e.g., with a generic Galerkin approach) may be problematic. We therefore proceed under the assumption of separability, that is typi-cally postulated for unsteady problems to separate the space to the time dependence of the solution (semidiscretization). This means that the reduced solution is regarded as the linear combination of factors breaking up the dependence on the different independent variables into the product of low-dimensional separated functions. Consequently, the PGD approximation of the problem takes the form

u m (x; q) = m k=1 u x k (x) δ l=1 u l k (q l ), (4)
for m ∈ N + , where functions u x k ∈ V h ⊂ V account for the x-dependence of the solution, while functions

u l k ∈ Q l ⊂ L 2 (S l
), for l = 1, 2, . . . δ, carry the dependence on the parameters, V h and Q l being discrete spaces with dim(V h ) = N x h and dim(Q l ) = N l , respectively. The separability assumption is applied also to the extended bilinear and linear forms in [START_REF] Bertagna | Data assimilation in cardiovascular fluid-structure interaction problems: an introduction[END_REF], so that we have

A(u, v) = Na α=1 a x α (u x ψ , v x β ) δ l=1 a l α (u l ψ , v l β ), F(v) = N F ϕ=1 F x ϕ (v x β ) δ l=1 F l ϕ (v l β ), (5)
with in a separate form as well [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations: a Primer[END_REF]. Now, we define the space of the functions in (4), i.e. the PGD space,

a x α : V × V → R and a l α : L 2 (S l ) × L 2 (S l ) → R bilinear forms, for α = 1, . . . , N a , F x ϕ : V → R and F l ϕ : L 2 (S l ) → R,
W m = w m (x; q) = m k=1 w x k (x) δ l=1 w l k (q l ), with w x k ∈ V h , w l k ∈ Q l , l = 1, . . . , δ, x ∈ Ω, q ∈ S . (6)
Differently from a standard Galerkin approach, the separated functions, u x k , u l k in (4), are not selected apriori (e.g., as piecewise or Gaussian polynomials). Conversely, the computation of these functions is the result of a progressive construction customized on the problem to solve, performed in an iterative way that takes advantage of the factorization in (4). In more details, the computation of u m is based on two steps: (i) a "greedy" weighted residual step to enrich the approximation u m-1 ∈ W m-1 to u m ∈ W m , by adding the new contribution (or mode)

δ m u ≡ u x m (x) δ l=1
u l m (q l ); (ii) a factorized computation of the m-th term, δ m u, based on a fixed-point iteration method called Alternating-Direction Strategy (ADS) [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations: a Primer[END_REF][START_REF] Ladevèze | The LATIN multiscale computational method and the proper generalized decomposition[END_REF][START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the proper generalized decomposition[END_REF].

As for (i), we aim to find the new mode, δ m u, by solving the weighted residual problem

A(u m-1 + δ m u, v m ) = F(v m ), ∀v m ∈ W m , (7)
where u m-1 ∈ W m-1 is the available approximation (for m = 1, typically, one sets u 0 = 0). The progressive addition of new terms stops when the PGD solution no longer changes significantly. Technically, this leads to the (external) convergence criterion

u m u 1 ≤ tol e , (8)
where the tolerance tol e is user-defined and the norm can be selected in different ways (e.g., in

H 1 (Ω) ⊗ δ l=1 L 2 (S l )
). The rationale is that the new modes add a progressively less relevant information to the solution (as it happens, for instance, with a Sturm-Liouville Eigenvalue expansion [START_REF] Pinchover | An Introduction to Partial Differential Equations[END_REF]), so that the lefthand side of (8) does actually reduce when m increases.

As for (ii), we perform the computation of δ m u in [START_REF] Asch | Data assimilation: methods, algorithms, and applications[END_REF] by solving iteratively the component dependence on each variable (or set of variables) independently, in an alternating direction framework. Introducing the fixed-point iterative index, j ≥ 1, let us assume to have a guess δ m,j-1 u(x; q) = u x m,j-1 (x)

δ l=1
u l m,j-1 (q l ) for δ m u (for j = 1, δ m,0 u(x; q) is chosen according to the boundary conditions of the problem at hand). We start computing the contribution u x m,j by solving the d-dimensional problem

A(u x m,j δ l=1 u l m,j-1 , v m ) = F(v m ) -A(u m-1 , v m ), (9)
where v m (x; q) = v x m (x) δ l=1 u l m,j-1 (q l ), for any v x m ∈ V h . Successively, we tackle the dependence on the parameter q r , for r = 1, . . . , δ, by solving the 1D problem

A(u x m,j u r m,j r-1 l=1 u l m,j δ l=r+1 u l m,j-1 , v m ) = F(v m ) -A(u m-1 , v m ), (10)
where u r m,j is the unknown and we choose now v m (x; q) = u x m,j (x)v r m (q r ) r-1 l=1 u l m,j (q l ) δ l=r+1 u l m,j-1 (q l ), for any v r m ∈ Q r . By solving [START_REF] Quarteroni | Reduced basis methods for partial differential equations: an introduction[END_REF] for r = 1, . . . , δ, we complete the computation of δ m u in [START_REF] Asch | Data assimilation: methods, algorithms, and applications[END_REF], by setting

δ m u(x; q) = u x m (x) δ l=1 u l m (q l ), with u x m (x) = u x m, j (x), u l m (q l ) = u l m,
j (q l ), j denoting the fixed-point iteration index ensuring the convergence to problems ( 9)- [START_REF] Quarteroni | Reduced basis methods for partial differential equations: an introduction[END_REF]. A standard stopping criterion for the (internal) iterative loops in [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF] and [START_REF] Quarteroni | Reduced basis methods for partial differential equations: an introduction[END_REF] reads

δ m,j u -δ m,j-1 u δ m,j-1 u ≤ tol i , (11)
with tol i a user-defined tolerance and δ m,τ u(x; q) = u x m,τ (x) δ l=1 u l m,τ (q l ), with τ = j -1, j. Once this criterion is fulfilled, the PGD approximation u m (x; q) = u m-1 (x; q) + δ m u(x; q) is available. All the steps involved in ( 9)-( 10) can be written in an algebraic form that we skip here for the sake of brevity (see, e.g., [START_REF] Perotto | Model reduction by separation of variables: a comparison between hierarchical model reduction and proper generalized decomposition[END_REF] for the explicit computations).

Algorithm 1 PGD algorithm (Offline phase) Set: m ← m + 1;

5:

Set: j ← 0;

6:

Set δ m,0 u;

7:
while ( δ m,j u -δ m,j-1 u / δ m,j-1 u ≤ tol i & j ≤ j max ) do end while 14:

Set:

u x m (x) ← u x m,j (x); 15:
Set: u r m (q r ) ← u r m,j (q r );

16:

Set: u m (x; q) ← u m-1 (x; q) + u x m (x) δ r=1 u r m (q r ); 17: end while

The PGD procedure is summarized in Algorithm 1. Implementation details related specifically to the monodomain problem are deferred to Section 4.2. In particular, among the input parameters, m max and j max introduce an upper bound to the maximum number of enrichment and fixed-point iterations, respectively. Algorithm 1 is the so-called offline phase of the PGD method. After this phase, the reduced solution is available for any value of the independent variables, including the admissible set S of the parameters. Thus, the cost of the online phase reduces to the cost of assembling the solution u m (x; q) = m k=1

u x k (x) δ l=1 u l k (q l ).
This is a strategic property in view of a data assimilation procedure, where we need to iteratively estimate the solution of the problem of interest for different values of the parameters, to minimize the mismatch from the observations. At this time, we do not have theoretical results on the convergence of the loops involved in the offline phase. Actually, the ultimate goal of this work is assessing the reliability of the PGD model reduction, and verifying its computational competitiveness when solving the cardiac conductivities estimation problem.

The Monodomain Model

The Forward Problem

Propagation of the electrical signal in the cardiac tissue is the result of complex multiscale dynamics occurring over the whole heart, yet based on cellular kinetics. Well-established experimental evidences, originating from Hodgkin-Huxley famous studies, supported the description of the kinetics of ionic currents by means of local nonlinear reaction terms within a reaction-diffusion model [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. These reaction functions require constitutive laws generally denoted by ionic models. They consist of a set of ordinary differential equations capturing the dynamics of the ionic currents flowing throughout the cell membrane (generalization including intra-cellular reaction-diffusion processes, involving, e.g., Ca 2+ ions, are also present in the literature [START_REF] Keener | Mathematical Physiology[END_REF]). In particular, they describe the time evolution of the so-called gating variables, which control the ionic fluxes and ensure the fulfillment of ion balance laws over the whole cell.

The mathematical modeling of the electrical activity of the heart has been the subject of many works. One of the most popular descriptions is the Bidomain model [START_REF] Pullan | Mathematically modelling the electrical activity of the heart: from cell to body surface and back again[END_REF][START_REF] Sundnes | Computing the electrical activity in the heart[END_REF]. This stems from a homogenized description of both the intra-and the extra-cellular potentials, that are defined on the same spatial domain, dropping a fine microscopic description of the intra-and extra-cellular regions. The model has been demonstrated to be accurate and reliable (specifically for defibrillation applications [START_REF] Luther | Low-energy control of electrical turbulence in the heart[END_REF]), yet its degenerate nature of parabolic partial differential equation system (where the matrix multiplying the vector of time derivatives is actually singular) makes the efficient numerical solution quite challenging. This has justified an abundant literature on the efficient solution of the Bidomain model, in terms of algorithms, preconditioners and high-performance computing implementations (see, e.g., [START_REF] Colli Franzone | A parallel solver for reaction-diffusion systems in computational electrocardiology[END_REF][START_REF] Gerardo-Giorda | A model-based block-triangular preconditioner for the bidomain system in electrocardiology[END_REF][START_REF] Gerardo-Giorda | Optimized schwarz coupling of bidomain and monodomain models in electrocardiology[END_REF][START_REF] Linge | Numerical solution of the bidomain equations[END_REF][START_REF] Mirabella | An a posteriori error estimator for model adaptivity in electrocardiology[END_REF][START_REF] Pavarino | Multilevel additive schwarz preconditioners for the bidomain reactiondiffusion system[END_REF][START_REF] Pavarino | Parallel multilevel schwarz and block preconditioners for the bidomain parabolic-parabolic and parabolic-elliptic formulations[END_REF][START_REF] Perego | An efficient generalization of the rush-larsen method for solving electrophysiology membrane equations[END_REF][START_REF] Plank | Algebraic multigrid preconditioner for the cardiac bidomain model[END_REF][START_REF] Vigmond | Computational techniques for solving the bidomain equations in three dimensions[END_REF][START_REF] Vigmond | Solvers for the cardiac bidomain equations[END_REF][START_REF] Hurtado | Non-conforming finite-element formulation for cardiac electrophysiology: an effective approach to reduce the computation time of heart simulations without compromising accuracy[END_REF] to mention a few).

A more classical approach in computational electrocardiology, also able to reduce the computational costs, relies on the Monodomain model. This is actually a model reduction of the Bidomain equations formulated in terms of the transmembrane potential (the difference between intra-and extra-cellular voltages), under the assumption that the anisotropy of the intra-and extra-cellular spaces is the same, i.e., the conductivity in the extra-cellular region is proportional to the intra-cellular conductivity. Even if this assumption is generally not completely justified, this model has been recognized as a possible trade-off between accuracy and efficiency.

Besides, it gains the same level of information in the case of physiological spatio-temporal propagation (i.e., the object of the present study). For this reason, in this paper, we focus on the Monodomain system, coupled, in particular, with the simplified two-variable phenomenological Rogers-McCulloch (RM) model [START_REF] Clayton | Models of cardiac tissue electrophysiology: progress, challenges and open questions[END_REF], featuring a single gating variable. More complex models, both physiological and phenomenological, are available in the literature [START_REF] Fenton | Models of cardiac cell[END_REF][START_REF] Luo | A model of the ventricular cardiac action potential. depolarization, repolarization, and their interaction[END_REF][START_REF] Luo | A dynamic model of the cardiac ventricular action potential. i. simulations of ionic currents and concentration changes[END_REF][START_REF] Ten Tusscher | A model for human ventricular tissue[END_REF]. However, the purpose of this work is to assess a methodological procedure, so we defer the extension to more complex models to future works.

The Monodomain equations we refer to read

     ∂ t u = ∇ • (σ∇u) -I ion (u, w) + I app in Q, d t w = g(u, w) in Q, (12)
with initial and boundary conditions According to the standard notation, symbols ∂ t , d t , ∇ and ∇• represent the partial and total derivatives in time, the spatial gradient and divergence operator respectively. The Neumann-type boundary condition (13) 2 , usually adopted in this context, mimics an insulated tissue [START_REF] Pullan | Mathematically modelling the electrical activity of the heart: from cell to body surface and back again[END_REF]. In mathematical terms, Neumann-type conditions try to minimize the sensitivity of arbitrary data on the simulation results [START_REF] Cherry | Effects of boundaries and geometry on the spatial distribution of action potential duration in cardiac tissue[END_REF].

     u(x, 0) = u 0 (x), w(x, 0) = w 0 (x) in Ω, σ∇u • n = 0 on ∂Q, (13)
We assume the computational domain Ω to be homogeneous and anisotropic, meaning that the conductivity tensor σ is uniform in space. Also this assumption will be removed in future works [START_REF] Barone | Experimental validation of a variational data assimilation procedure for estimating space-dependent cardiac conductivities[END_REF]. We define with {a l , a t , a n } the orthonormal fields related to the structure of the myocardium which define the longitudinal, transverse and normal direction of the fibers, respectively [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF][START_REF] Pullan | Mathematically modelling the electrical activity of the heart: from cell to body surface and back again[END_REF]. Accordingly, we identify with (σ l , σ t , σ n ) the longitudinal, transverse and normal conductivity parameters. In the most general form, the conductivity tensor can be decomposed as

σ = σ l a l a T l + σ t a t a T t + σ n a n a T n . (14)
Moreover, for d = 3 and under the hypothesis of axial isotropy, i.e., σ t = σ n and a n ≡ e 3 with {e 1 , e 2 , e 3 } the canonical basis in R 3 , the conductivity tensor reduces to

σ = σ t I + (σ l -σ t )a l a T l , (15)
with I∈ R 3×3 the identity tensor. The conductivity values belong to separate admissible intervals, namely,

σ l ∈ [m l , M l ] and σ t ∈ [m t , M t ], ∆ l = M l -m l and ∆ t = M t -m t denoting the corresponding length, with M t < m l .
In fact, it is observed that electrical conduction along the length of myocytes is faster than along the transverse direction [START_REF] Spach | Anisotropy of cardiac tissue: a major determinant of conduction[END_REF][START_REF] Spach | Cell size and communication: role in structural and electrical development and remodeling of the heart[END_REF], so that σ l should be greater than σ t .

Let θ : Ω → [0, 2π] be the angle between the longitudinal fiber direction, a l , and the x-axis direction, e 1 , at any point of the domain Ω. It holds that a l = (cos θ)e 1 + (sin θ)e 2 . Therefore, the conductivity tensor

RM parameters G V th V p η 1 η 2 η 3 1.5 Ω -1 cm -2 13 mV 100 mV 4.4 Ω -1 cm -2 0.012 ms -1 1 ms -1 Table 1: Rogers-McCulloch model parameters.
Here Ω denotes the Ohm, the unit of electrical resistance.

explicitly reads σ =       σ l cos 2 θ + σ t sin 2 θ (σ l -σ t ) cos θ sin θ 0 (σ l -σ t ) cos θ sin θ σ l sin 2 θ + σ t cos 2 θ 0 0 0 σ t       for d = 3, (16)
and

σ =    σ l cos 2 θ + σ t sin 2 θ (σ l -σ t ) cos θ sin θ (σ l -σ t ) cos θ sin θ σ l sin 2 θ + σ t cos 2 θ    for d = 2. ( 17
)
The RM model is based on a cubic polynomial formulation for I ion and one gating variable which allow to describe the most relevant feature of the cardiac action potential at a limited computational cost. The model reads

I ion (u, w) = Gu 1 - u V th 1 - u V p + η 1 uw, g(u, w) = η 2 u V p -η 3 w . (18)
The parameters G, V th , V p , η 1 , η 2 , η 3 are provided in Table 1.

The Monodomain Inverse Conductivity Problem

The ultimate goal of our research stems from the clinical need of identifying the cardiac conductivities, σ l , σ t , in a patient-specific setting. This problem presents several practical challenges, for the limited access to data in vivo. Also, the available literature in the field shows that the range of these parameters may be quite large [START_REF] Clerc | Directional differences of impulse spread in trabecular muscle from mammalian heart[END_REF][START_REF] Roberts | Influence of cardiac fiber orientation on wavefront voltage, conduction velocity, and tissue resistivity in the dog[END_REF][START_REF] Roberts | Effect of tissue anisotropy on extracellular potential fields in canine myocardium in situ[END_REF], while the numerical sensitivity of the Monodomain solution to cardiac conductivities is generally high [START_REF] Johnston | A sensitivity study of conductivity values in the passive bidomain equation[END_REF].

A variational data assimilation approach was proposed and analyzed in [START_REF] Graham | Estimation of the bidomain conductivity parameters of cardiac tissue from extracellular potential distributions initiated by point stimulation[END_REF][START_REF] Yang | Estimation of cardiac conductivities in ventricular tissue by a variational approach[END_REF]. The basic idea is to infer the conductivities from available measures of the transmembrane potential retrieved from the surface of the tissue at certain instants of the propagation. The variational procedure achieves this assimilation by finding the conductivities that minimize the mismatch between the available data and the results of the Monodomain solution. Formally, this leads to the so-called Monodomain Inverse Conductivity Problem (MICP): find σ l and σ t in the tensor σ to minimize the functional

J (σ) = 1 2 T 0 Ω obs (u(σ) -u meas ) 2 dxdt, (19)
subject to ( 12)-( 13). Specifically, u meas denotes the experimental data measured on the observation domain, Ω obs ⊂ Ω, and u is function of σ through the [START_REF] Luo | Proper orthogonal decomposition methods for partial differential equations[END_REF][START_REF] Costabal | Multi-fidelity classification using gaussian processes: Accelerating the prediction of large-scale computational models[END_REF][START_REF] Lenarda | A modeling framework for electro-mechanical interaction between excitable deformable cells[END_REF]..

An existence analysis for this problem is reported in [START_REF] Yang | Estimation of cardiac conductivities in ventricular tissue by a variational approach[END_REF], whereas numerical and experimental validations are extensively discussed in [START_REF] Barone | Numerical sensitivity analysis of a variational data assimilation procedure for cardiac conductivities[END_REF][START_REF] Barone | Experimental validation of a variational data assimilation procedure for estimating space-dependent cardiac conductivities[END_REF]. As it is promptly realized, the iterative minimization procedure based on the introduction of Lagrange multipliers and the solution of the Monodomain adjoint problem is computationally demanding. As a matter of fact, when following an optimize-then-discretize approach, we need to evaluate the Monodomain system and its adjoint, forward and backward in time, at each minimization iteration. This led to the introduction of model reduction techniques, based either on a Proper Orthogonal Decomposition (POD) paradigm [START_REF] Boulakia | Reduced-order modeling for cardiac electrophysiology. application to parameter identification[END_REF][START_REF] Pagani | Numerical approximation of parametrized problems in cardiac electrophysiology by a local reduced basis method[END_REF][START_REF] Yang | Efficient estimation of cardiac conductivities via pod-deim model order reduction[END_REF] or the Lax-pairs [START_REF] Gerbeau | Reduced order model in cardiac electrophysiology with approximated lax pairs[END_REF]. The POD paradigm requires the offline computation of snapshots for different values of the parameters. Successively, these snapshots are reduced via SVD and combined to form a basis for a rapid evaluation of the solution for new values of the parameters.

As demonstrated in [START_REF] Yang | Efficient estimation of cardiac conductivities via pod-deim model order reduction[END_REF], the selection of the snapshots is critical for the successful achievement of an efficient reduced solution, and specific techniques are required [START_REF] Pagani | Numerical approximation of parametrized problems in cardiac electrophysiology by a local reduced basis method[END_REF].

On the contrary, as already pointed out, the PGD approach does not need the computation of any snapshot in the offline phase, and directly computes the solution in the incremented independent-variable space (in our case, the extended independent variables are the space, the time and the two conductivities).

The actual performance of PGD will be illustrated in the next section.

PGD Model Reduction of the Monodomain Problem

Formulation of the Reduced Model

To solve the Monodomain system with the PGD approach, we perform first a semi-discretization in time of the problem. This is, a priori, not necessary, as the time could be considered as an independent variable of the problem (like the space and the conductivities). However, we prefer to eliminate the time dependence by resorting to a traditional finite difference discretization. Also, we decouple the PDE from the ODE ionic model. By defining a time step ∆t > 0, the time interval [0, T ] is uniformly divided into N subintervals, (t n , t n+1 ), with t 0 = 0, t N = T , being t n+1 = t n + ∆t = (n + 1)∆t, for any n = 0, . . . , N -1. A popular strategy for the semi-discretization in time is a semi-implicit approach that automatically linearizes the problem at each time step. With this approach and using a first order discretization in time, at each time t n+1 , for (u n , w n ) available, we solve the problems

     w n+1 = ∆t g(u n , w n+1 ) + w n in Ω, u n+1 -∆t ∇ • (σ∇u n+1 ) = ∆t I n+1 app -∆t I ion (u n , w n+1 ) + u n in Ω, (20)
completed by the boundary conditions (13) 2 . For each n = 0, . . . , N -1, (18,20) lead to

w n+1 = η 2 ∆t (1 + ∆t η 2 η 3 )V p u n + w n 1 + ∆t η 2 η 3 . (21)
Consequently, the ionic term reads

I ion (u n , w n+1 ) = Gu n -G 1 V p + 1 V th (u n ) 2 + G V p V th (u n ) 3 + η 1 u n w n+1 . ( 22
)
At each step, we apply the PGD model reduction to (20) 2 , whose weak formulation, for V = H 1 (Ω), reads

a σ (u n+1 , v) = F σ (v) ∀v ∈ V, (23)
with

a σ (u n+1 , v) = Ω ∆t σ∇u n+1 • ∇v dΩ + Ω u n+1 v dΩ, (24)
F σ (v) = Ω (∆t I n+1 app -∆t I ion (u n , w n+1 ) + u n )v dΩ. ( 25
)
Notice that, for the sake of generality, we retain the index σ for the functional F in ( 23) to be consistent with our notation in Section 2. In fact, in our specific problem, the functional at the right hand side is independent of the parameters.

To apply the PGD approach, we introduce the space W ≡ H 1 (Ω) ⊗ L 2 (S l ) ⊗ L 2 (S t ), where S l and S t are the admissible ranges for σ l and σ t , respectively. Then, the extended variational formulation (2) features

A(u, v) = S l St a σ (u, v) dσ l dσ t , F(v) = S l St F σ (v) dσ l dσ t . (26)
Again, we preserve the same notation, u and v, despite the dependence on σ l and σ t . To perform the model reduction, we introduce the PGD space

W m = w m (x; σ l , σ t ) = m k=1 w x k (x)w σ l k (σ l )w σt k (σ t ), with w x k ∈ V h , w σ l k ∈ Q σ l h , w σt k ∈ Q σt h , x ∈ Ω, σ l ∈ S l , σ t ∈ S t , (27)
with V h , Q σ l h and Q σt h discrete subspaces of V , L 2 (S l ) and L 2 (S t ). The PGD approximation for the solution at time t = t n+1 is defined as

u n+1 (x; σ l , σ t ) ≈ m n+1 u k=1 u n+1,x k (x)u n+1,σ l k (σ l )u n+1,σt k (σ t ). ( 28
)
We observe that the modal index m n+1 u may vary, a priori, at each time-step. The different contributions of the sum are computed up to the fulfillment of criterion (8), while each term

u n+1,x k , u n+1,σ l k , u n+1,σt k , is
computed by the ADS strategy, breaking the solution into an iterative fixed-point solver, and alternatively solving for x, σ l and σ t , as described in Section 2.

The factorization of the forms A and F in (5) leads to

A(u, v) = 10 α=1 a x α (u x ψ , v x β )a l α (u l ψ , v l β )a t α (u t ψ , v t β ), F(v) = N F α=1 F x α (v x β )F l α (v l β )F t α (v t β ), (29)
with ψ and β varying in the corresponding modal index range, and where the bilinear and linear factors are explicitly provided in Table 2 and 3, respectively. The notation (•, •) stands for the standard L 2 -product in Ω. For all the details concerning the derivation of the quantities in Tables 2 and3, we refer to [START_REF] Carlino | Model reduction by proper generalized decomposition in electrocardiology[END_REF]. In particular, the factorization of the linear functional F(v) demands a separated representation also for w,

given by

w n+1 (x; σ l , σ t ) ≈ m n+1 w k=1 w n+1,x k (x)w n+1,σ l k (σ l )w n+1,σt k (σ t ), (30)
although, in practice, the gating variable is computed via equation [START_REF] Graham | Estimation of the bidomain conductivity parameters of cardiac tissue from extracellular potential distributions initiated by point stimulation[END_REF].

By a direct computation for N F in [START_REF] Yang | Efficient estimation of cardiac conductivities via pod-deim model order reduction[END_REF], one verifies that

N F = 1 + 2m u + 2 m u 1 + 3 m u 2 + m u 3 + m u m w = 1 + 4m u + 3 2 m u (m u -1) + 1 6 m u (m u -1)(m u -2) + m u m w ,
the time index being omitted to simplify notation. Actually, in [START_REF] Barone | Experimental validation of a variational data assimilation procedure for estimating space-dependent cardiac conductivities[END_REF], I app contributes with a unique function, while the term u n requires m u functions. For the ionic term [START_REF] Abidi | Ionic parameters estimation in multi-scale cardiac electrophysiology modelling[END_REF], the linear term requires m u functions, the quadratic contribution demands mu Exploiting factorizations [START_REF] Pagani | Numerical approximation of parametrized problems in cardiac electrophysiology by a local reduced basis method[END_REF] and [START_REF] Yang | Efficient estimation of cardiac conductivities via pod-deim model order reduction[END_REF] in the inner loop ( 9)-( 10), the generic fixed-point iteration, j, of the [START_REF] Boulakia | Reduced-order modeling for cardiac electrophysiology. application to parameter identification[END_REF]. Lines for α = 1, . . . , 9 refer to the second order term in [START_REF] Barone | Numerical sensitivity analysis of a variational data assimilation procedure for cardiac conductivities[END_REF], while the last line (α = 10), refers to the reactive contribution.

α a x α (u x ψ , v x β ) a l α (u l ψ , v l β ) a t α (u t ψ , v t β ) 1 ∆t(cos 2 (θ)∂ x u x ψ , ∂ x v x β ) (σ l u l ψ , v l β ) (u t ψ , v t ) 2 ∆t(sin 2 (θ)∂ x u x ψ , ∂ x v x β ) (u l ψ , v l β ) (σ t u t ψ , v t β ) 3 ∆t(cos (θ) sin (θ)∂ y u x ψ , ∂ x v x β ) (σ l u l ψ , v l β ) (u t ψ , v t β ) 4 -∆t(cos (θ) sin (θ)∂ y u x ψ , ∂ x v x β ) (u l ψ , v l β ) (σ t u t ψ , v t β ) 5 ∆t(cos (θ) sin (θ)∂ x u x ψ , ∂ y v x β ) (σ l u l ψ , v l β ) (u t ψ , v t β ) 6 -∆t(cos (θ) sin (θ)∂ x u x ψ , ∂ y v x β ) (u l ψ , v l β ) (σ t u t ψ , v t β ) 7 ∆t(sin 2 (θ)∂ y u x ψ , ∂ y v x β ) (σ l u l ψ , v l β ) (u t ψ , v t β ) 8 ∆t(cos 2 (θ)∂ y u x ψ , ∂ y v x β ) (u l ψ , v l β ) (σ t u t ψ , v t β ) 9 ∆t(∂ z u x ψ , ∂ z v x β ) (σ l u l ψ , v l β ) (u t ψ , v t β ) 10 (u x ψ , v x β ) (u l ψ , v l β ) (u t ψ , v t β )
Table 3: Factorization of the PGD extended linear form in [START_REF] Boulakia | Reduced-order modeling for cardiac electrophysiology. application to parameter identification[END_REF]. The first column provides the number of components for any row. As for the indices, we have k, γ, ζ = 1, . . . , m n u and ξ = 1, . . . , m n+1 w .

F x α (v x β ) F l α (v l β ) F t α (v t β ) 1 ∆t(I n+1 app , v x β ) (1, v l β ) (1, v t β ) m u (u n,x k , v x β ) (u n,σ l k , v l β ) (u n,σt k , v t β ) m u -∆tG(u n,x k , v x β ) (u n,σ l k , v l β ) (u n,σt k , v t β ) m u ∆tG( 1 Vp + 1 V th )((u n,x k ) 2 , v x β ) ((u n,σ l k ) 2 , v l β ) ((u n,σt k ) 2 , v t β ) mu 2 2∆tG( 1 Vp + 1 V th )(u n,x k u n,x γ , v x β ) (u n,σ l k u n,σ l γ , v l β ) (u n,σt k u n,σt γ , v t β ) with γ > k m u -∆t G VpV th ((u n,x k ) 3 , v x β ) ((u n,σ l k ) 3 , v l β ) ((u n,σt k ) 3 , v t β ) 2 mu 2 -3∆t G VpV th ((u n,x k ) 2 u n,x ζ , v x β ) ((u n,σ l k ) 2 u n,σ l ζ , v l β ) ((u n,σt k ) 2 u n,σt ζ , v t β ) with ζ = k mu 3 -6∆t G VpV th (u n,x k u n,x ζ u n,x γ , v x β ) (u n,σ l k u n,σ l ζ u n,σ l γ , v l β ) (u n,σt k u n,σt ζ u n,σt γ , v t β ) with γ > ζ > k m u m w ∆tη 1 (u n,x k w n+1,x ξ , v x β ) (u n,σ l k w n+1,σ l ξ , v l β ) (u n,σt k w n+1,σt ξ , v t β ) three-step ADS reads                                                  10 α=1 a x α (u n+1,x m,j , v x m )a l α (u n+1,σ l m,j-1 , u n+1,σ l m,j-1 )a t α (u n+1,σt m,j-1 , u n+1,σt m,j-1 ) = N F α=1 F x α (v x m )F l α (u n+1,σ l m,j-1 )F t α (u n+1,σt m,j-1 ) - 10 α=1 m-1 k=1 a x α (u n+1,x k , v x m )a l α (u n+1,σ l k , u n+1,σ l m,j-1 )a t α (u n+1,σt k , u n+1,σt m,j-1 ) ∀v x m ∈ V h , 10 α=1 a x α (u n+1,x m,j , u n+1,x m,j )a l α (u n+1,σ l m,j , v l m )a t α (u n+1,σt m,j-1 , u n+1,σt m,j-1 ) = N F α=1 F x α (u n+1,x m,j )F l α (v l m )F t α (u n+1,σt m,j-1 ) - 10 α=1 m-1 k=1 a x α (u n+1,x k , u n+1,x m,j )a l α (u n+1,σ l k , v l m )a t α (u n+1,σt k , u n+1,σt m,j-1 ) ∀v l m ∈ Q σ l h , 10 α=1 a x α (u n+1,x m,j , u n+1,x m,j )a l α (u n+1,σ l m,j , u n+1,σ l m,j )a t α (u n+1,σt m,j , v t m ) = N F α=1 F x α (u n+1,x m,j )F l α (u n+1,σ l m,j )F t α (v t m ) - 10 α=1 m-1 k=1 a x α (u n+1,x k , u n+1,x m,j )a l α (u n+1,σ l k , u n+1,σ l m,j )a t α (u n+1,σt k , v t m ) ∀v t m ∈ Q σt h , (31)
where, to simplify the notation, we have set m n+1 u = m. Moreover, we have highlighted in bold the term each equation has to be solved for.

Implementation Details

For the discrete spaces V h , Q σ l h and Q σt h in [START_REF] Gerbeau | Reduced order model in cardiac electrophysiology with approximated lax pairs[END_REF] we choose piecewise linear finite elements (FE) [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]. At the time t = t n+1 , the initial guess for the external loop is set to u n+1,x The linear systems following the discretization of the three ADS steps in [START_REF] Ammar | Reduction of the chemical master equation for gene regulatory networks using proper generalized decompositions[END_REF] are solved using the conjugate gradient method with a standard incomplete LU (ILU) right preconditioner [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. While this strategy is appropriate to solve the first ADS step, (31) 1 , that may lead to large linear systems associated with the discretization of the physical space, the remaining steps concern the parametric space and, in general, they require the solution of smaller linear systems. Direct solvers may be more efficient in such cases, depending on the linear algebra package used. The optimization of this part of the implementation will be object of future works. Simulations were carried out in serial, on a workstation equipped with Intel 6-Core i7-7800X CPU 3.50GHz and 64 GB of RAM. The code was implemented in LifeV [START_REF] Bertagna | The lifev library: engineering mathematics beyond the proof of concept[END_REF][START_REF] Passerini | Validation of an open source framework for the simulation of blood flow in rigid and deformable vessels[END_REF], an object oriented C++ parallel finite element library based on the Trilinos project [START_REF] Heroux | An Overview of Trilinos[END_REF], developed by different groups worldwide.

Numerical Results for the Reduction of the Forward Problem

Here we investigate the performance of the PGD technique to reduce the Monodomain model [START_REF] Luo | Proper orthogonal decomposition methods for partial differential equations[END_REF] in terms of accuracy of the solution and computational efficiency. We focus on realistic geometries, both in 2D and 3D.

A 2D Canine Tissue Geometry

We consider a realistic 2D geometry of a portion of a canine ventricular tissue [START_REF] Gizzi | Effects of pacing site and stimulation history on alternans dynamics and the development of complex spatiotemporal patterns in cardiac tissue[END_REF]. The computational domain features 22747 degrees of freedom (DOFs). Fig. 1(a) shows the mesh and the cardiac fiber structure that was roughly approximated by looking at the anatomy of the tissue. As previously mentioned, the cardiac tissue is assumed homogeneous, so that the conductivity fields are constant and identified by the 2D vector σ = (σ l , σ t ) with the local preferred orientation represented in the right panel of Fig. 1(a). The parameters, σ l and σ t , belong to the admissible space S = S l × S t = [m l , M l ] × [m t , M t ] = [0.06, 0.09] cm 2 /ms × [0.01, 0.04] cm 2 /ms, (∆ l = ∆ t = 0.03 cm 2 /ms). These bounds were manually tuned so to reproduce the realistic wavefront propagation velocities observed in vitro experiments [START_REF] Gizzi | Effects of pacing site and stimulation history on alternans dynamics and the development of complex spatiotemporal patterns in cardiac tissue[END_REF][START_REF] Fenton | Termination of atrial fibrillation using pulsed low-energy far-field stimulation[END_REF]. Domain S has been discretized with 250 × 250 DOFs, while the simulation time step is set to ∆t = 0.2 ms until T = 30 ms. One stimulus of I app = 250 mV/ms is applied at the top of the domain for a duration of 2 ms. For the sake of the computational costs, we focus on the wavefront propagation, whereas insights into the PGD approximation for the whole action potential are provided in Appendix A. One of the primary goals of the following numerical experiments is to investigate the interplay between the tolerance tol e (associated to the number of modes) and the reliability of the solution. In fact, the lower the tolerance, the more accurate is expected to be the PGD solution. However, this entails higher computational costs in the offline phase as more modes are required to converge. We compare the full FE Monodomain solution, assumed as the reference solution, with the PGD approximation, varying tol e =10 -4 , 10 -5 ,10 -6 , and the conductivities. The tolerance of the ADS fixed point iterations, tol i , is set to 10 -2 .

Fig. 2(a) shows the number of modes as a function of the time. We notice that more modes are needed as the excitable wave travels through the tissue and the dynamics become more involved. Then, this number suddenly drops at around t=22 ms, when the wavefront propagation terminates. As expected, the number of modes needed for convergence increases with lower tolerance values. This results in a heavier computational burden of the offline phase, as highlighted in Table 4. However, the computation of the PGD solution in the online phase is remarkably inexpensive. Evaluating the PGD approximation at a new pair of parameters, (σ l , σ t ), only takes at most 0.2 s in contrast with 15 s of the corresponding FE solution, thus reducing the computation time of two orders of magnitude. The easy evaluation of the PGD solution is extremely helpful to set up a fast solver for the inverse problem, as we will verify in the next section.

In terms of accuracy, the quality of the PGD solution depends on the conductivity values which, in this particular case where we do not perform any change in the model reaction parameters, are mostly determined by the parameter σ t . Fig. 2(b) displays the percentage L 2 -relative error between the FE and the PGD transmembrane potential obtained with tol e =10 -4 , σ l = 0.09 cm 2 /ms, and for σ t = 0.01, 0.02, 0.03, 0.035, 0.04 cm 2 /ms, corresponding to the conduction velocities (CVs) = 0.15, 0.21, 0.25, 0.27, 0.29 cm/ms, ordered from the lowest to the highest value. Regardless of the wavefront velocity, the discrepancy between the FE and the PGD approximation is minimal after the excitation ends. On the other hand, during the potential propagation, the error increases at faster CVs, reaching almost 40% for σ t = 0.04 cm 2 /ms.

The accuracy of the PGD solution improves for slow CVs. For instance, for σ t = 0.01, 0.02 cm 2 /ms, the error is always below 10%. Therefore, we postulate that the PGD basis is informative enough to reproduce slow excitation waves, yet it needs further enrichment to accurately capture faster wavefront propagation. This is confirmed in Fig. 2(c) comparing the percentage L 2 -relative error on the potential when reducing the tolerance. Also, let us introduce the anisotropy ratio R a defined as R a ≡ σ l /σ t . The solution of the Monodomain equation in known to be sensitive to R a (see e.g. [START_REF] Yang | Efficient estimation of cardiac conductivities via pod-deim model order reduction[END_REF]), so we analyze three different values typical for the cardiac tissue [START_REF] Pullan | Mathematically modelling the electrical activity of the heart: from cell to body surface and back again[END_REF], R a = 2, 4, 6, corresponding to σ = (0.07, 0.035) cm 2 /ms, σ =(0.08, 0.02) cm 2 /ms and σ =(0.084, 0.014) cm 2 /ms, respectively. The case R a = 2 features the highest CVs as σ t takes the greatest value (0.035 cm 2 /ms), whereas the cases R a = 4, 6 result in a slower propagation.

The approximation is overall inaccurate at the initial stages of the propagation because of the lack of regularity of the stimulus function I app that abruptly goes to zero after 2 ms. Then, for the high CV case, R a = 2, the maximum of the error decreases from roughly 20% for tol e =10 -4 to approximately 10% and 3%

for tol e =10 -5 , 10 -6 , respectively. Therefore, enriching the PGD basis is necessary to ensure a reasonable accuracy in the case of high CVs. For low CV cases, R a = 4, 6, the approximation slightly improves in case of a richer PGD basis. However, setting tol e =10 -4 already guarantees a good approximation as the error is below 10%. A qualitative comparison between FE and PGD transmembrane potentials varying tolerance tol e and anisotropy ratio is provided in Fig. 3. The snapshots shown are related to the times featuring the highest error. As for the PGD solution obtained with tol e =10 -4 , the CVs are accurately captured for the slow CV cases, R a = 4, 6 in (b) and (c), whereas, for the high CV case, R a = 2 in (a), the PGD wavefront propagation is slower than the FE one, so that CV is underestimated. As expected, a smaller tolerance improves the accuracy, with an associated increment of the computational cost (compare the PGD solution for tol e =10 -6 with the FE approximation). In general, setting tol e =10 -5 seems to be the most convenient trade-off between accuracy and efficiency. The absolute value of the error between the FE and the PGD potential for the high CV case, R a = 2, is shown in Fig. 4 for two different times. The error increases as the wavefront propagation approaches the boundary independently of the selected (external) tolerance, confirming the influence of the boundary conditions on the simulations [START_REF] Cherry | Effects of boundaries and geometry on the spatial distribution of action potential duration in cardiac tissue[END_REF].

An additional investigation aims at understanding the dependence of the accuracy on other parameters than the tolterance tol e . Since tol e =10 -4 gives a good computational efficiency, we keep this value and explore possible ways to improve the quality of the resulting PGD approximation in case of high CVs by changing other numerical parameters. In particular, three strategies are investigated:

(a) refining the discretization of the admissible space S;

(b) lowering the tolerance, tol i , for the inner fixed point iterations;

(c) narrowing the admissible interval for σ t since the CVs mostly depends on it.

The first two strategies attain at the selection of discretization parameters, while the third one refers to our a priori knowledge on the parameter to estimate.

For strategy (a), S is discretized with 500 × 500 as opposed to 250 × 250 DOFs used for the previous tests. Regarding (b), we test the accuracy of the PGD solution decreasing the tolerance, tol i , from 10 -2 to 10 -4 . The percentage L 2 -relative error on the potential yielded by strategies (a) and (b) is shown in Fig. 5(a) for σ = (0.07, 0.035) cm 2 /ms. Both the approaches do not result in a lower approximation error.

In fact, the refinement of the partition of S do not give any improvement, while a lower tolerance tol i is slightly beneficial only in the initial stage of the wave propagation. We also notice that these two strategies entail higher computational costs, either for the size of each problem (strategy (a)) or the number of nternal iterations (strategy (b)).

On the other hand, approach (c) successfully reduces the error. In Fig. 5(b), we analyze the performance of PGD in terms of percentage L 2 -relative error on the potential and number of modes, for ∆ t = 0.01, 0.02 cm 2 /ms, and compare this trend with the results obtained for ∆ t = 0.03 cm 2 /ms, ∆ l being set to 0.03 cm 2 /ms. Narrower ∆ t yields an improvement of the PGD approximation, with ∆ t = 0.01 cm 2 /ms leading to the most accurate results (the error reduces from roughly 20% to 3%). This is reasonable since restricting 20 the admissible interval for σ t implies less values to explore. Therefore, the more a priori knowledge we have on the parameter values, the better the approximation becomes. Strategy (c) is also beneficial to efficiency, as fewer modes are needed for convergence thus containing the computational demand of the offline phase.

A 3D Canine Tissue Geometry

To assess the impact of the geometry size on the quality and efficiency of the PGD, we consider a 3D canine tissue geometry obtained by extruding, along the z-axis, the 2D domain used in the previous section.

The resulting mesh is 0.5 cm thick and discretized with 136482 DOFs. The cardiac fiber structure on the top surface of the mesh is the same as in the 2D case and then it is extruded along the z-axis. We consider the most general admissible space for the conductivity parameters S = [0.06, 0.09] cm 2 /ms × [0.01, 0.04] cm 2 /ms, discretized with 250 × 250 DOFs. As done in the 2D test, we analyze the anisotropy ratios R a = 2, 4, 6 corresponding to σ = (0.07, 0.035) cm 2 /ms, σ = (0.08, 0.02) cm 2 /ms and σ = (0.084, 0.014) cm 2 /ms, respectively. Having more degrees of freedom in the mesh of the domain when compared to the 2D case clearly results in an increment of the CPU time. Table 5 shows that the offline phase actually takes 8 and 41 hours for tol e =10 -4 and 10 -5 , respectively. The case tol e =10 -6 considered in the 2D setting was not explored here because computationally unaffordable. As for the online phase, adopting the PGD technique is extremely convenient since the evaluation of the reduced solution requires at most 0.75 s as opposed to 330 s for the FE approximation.

As shown in Fig. 6, the trend of the error is rather similar to the 2D case, when varying the tolerance and the anisotropy ratio (compare with Fig. 2(c)). For tol e =10 -4 , the case R a = 2, corresponding to high CVs, features higher error than for the other two values. Decreasing the tolerance to 10 -5 qualitatively improves the solution. A qualitative comparison between the FE and the PGD transmembrane potential for different tolerances and conductivities is carries out in Fig. 7. As already noted in the 2D experiment, in the high CV case, R a = 2, the PGD solution obtained with tol e =10 -4 underestimates the CVs, yielding a poor approximation of the FE solution. On the contrary, in the cases R a = 4, 6, although the CVs are slightly overestimated, the PGD solution is sufficiently close to the FE approximation. Overall, the PGD approximation improves for tol e =10 -5 providing a better matching of the wavefront propagation of the FE solution. Finally, we remark that the anisotropy ratio R a does not have a significant impact on the accuracy of the PGD approach.

A Real Ventricle

In view of clinical applications, we run simulations on a real left ventricular geometry reconstructed from SPECT images [START_REF] Gerardo-Giorda | A model-based block-triangular preconditioner for the bidomain system in electrocardiology[END_REF]. The excitation wave was simulated on a mesh with 22470 DOFs. Fig. 1(b) displays the 3D mesh and a realistic representation of the fiber structure used in the simulation. The fiber orientation was first obtained on an ellipsoidal domain and then adapted to the real domain, following the strategy proposed in [START_REF] Gerardo-Giorda | A model-based block-triangular preconditioner for the bidomain system in electrocardiology[END_REF]. The conductivity values are the same as in the previous tests. In this case, the resulting CVs are similar, whereas the influence of the anisotropy ratio R a on the wavefront curvature is more evident.

One stimulus of I app = 250 mV/ms is applied at the ventricular endocardium, for a duration of 2 ms. The electrical propagation simulated with the FE method for σ = (0.08, 0.02) cm 2 /ms is shown in Fig. 8. Since the adopted spatial discretization is coarser than the one used for the 3D canine geometry, we also included tol e =10 -6 in our analysis. As shown in Table 6, the offline phase takes 0.8, 4.1 and 32 hours for tol e =10 -4 , 10 -5 , 10 -6 , respectively. On the other hand, we emphasize that the online phase is extremely inexpensive as it requires at most 0.2 s as opposed to 60 s of the FE simulation. The trend of the error on the potential between the FE and the PGD solution when varying the tolerance, tol e , and the anisotropy ratio is displayed in Fig. 9. Unlike the previous tests, the three anisotropy ratios show a similar error for tol e =10 -4 , with a loss of accuracy for R a = 4, 6 in comparison to the canine tissue tests. However, the approximation is overall reliable as the error is around 10%. The negative impact of the discontinuity of I app on the accuracy of the PGD technique is more evident in these experiments as the error at the initial stage of the simulation is around 30%. The reliability of the approximation does not significantly benefit from lower tolerance values of tol e . Only the case R a = 2 presents a slight reduction of the error, although less significant than the tests with the canine geometry. For R a = 4 and 6 we do not have substantial improvement. We argue that, in these cases, the accuracy is determined by the numerical errors of the space discretization and the linear systems solution. The results are confirmed in Fig. 10 with a visual comparison of the different solutions. Notice that a higher error is observed at the ventricular apex, in particular in the cases R a = 4, 6, because of the complex swirling pattern of the cardiac fibers and the geometric curvature of the ventricle that affect the accuracy of the simulations.

Solving the PGD-MICP

We verify now that PGD may significantly accelerate the minimization of the misfit functional [START_REF] Weinberg | Ephaptic coupling rescues conduction failure in weakly coupled cardiac tissue with voltagegated gap junctions[END_REF] when solving the MICP. Once the PGD solution is available, the value of the functional J as a function of the conductivities is readily computed. This enables the use of nonlinear constrained optimization algorithms, not necessarily developed for differential problems. Specifically, we tested our method with a generic solver like the sequential quadratic programming (SQP) method [START_REF] Nocedal | Numerical optimization[END_REF], implemented by the fmincon routine in MATLAB R . The possibility of using this kind of solvers is clearly a practical advantage of the online phase.

We present the results about the conductivity evaluation in a synthetic setting, where the data used as observations are generated on a spatial discretization more refined in comparison with the one used for the model reduction. At each time step, we also add a Gaussian noise with zero mean and standard deviation equal to p max x,t |u|, where p is the noise-to-signal ratio. Synthetic measurements were recorded every dt snap = 2 ms for a global duration of T = 30 ms, so that 15 voltage recordings are used to calculate J . The observation domain, Ω obs , consists of 8000 points equally distributed in the domain, which is comparable with the number of observation points characterizing standard optical mapping recordings [START_REF] Gizzi | Effects of pacing site and stimulation history on alternans dynamics and the development of complex spatiotemporal patterns in cardiac tissue[END_REF].

The search is constrained to the admissible parametric space S = S l × S t with S l = [0.06, 0.09] cm 2 /ms, S t = [0.01, 0.04] cm 2 /ms (∆ l = ∆ t = 0.03 cm 2 /ms). We set σ 0 = (0.06, 0.025) cm 2 /ms as initial guess. We evaluate the performance of the PGD method to solve the MICP varying the tolerance tol e and the anisotropic ratio, R a = 2, 4, 6 corresponding to σ = (0.07, 0.035) cm 2 /ms, σ =(0.08, 0.02) cm 2 /ms and σ =(0.084, 0.014) cm 2 /ms, respectively. The computational cost for solving the MICP with PGD is negligible and cannot be directly compared with the variational technique used in [START_REF] Barone | Numerical sensitivity analysis of a variational data assimilation procedure for cardiac conductivities[END_REF][START_REF] Barone | Experimental validation of a variational data assimilation procedure for estimating space-dependent cardiac conductivities[END_REF], as PGD takes advantage of the offline phase. However, we will demonstrate that the PGD-based solution of the MICP is advantageous even when including offline costs.

A 2D Canine Tissue Geometry

Considering the 2D canine mesh used in Section 4.3, we expect to get a more precise estimation for R a = 4 since the PGD approximation showed to be more accurate than for the other choices of the conductivity pairs (see Fig. 2(c)).

The reliability of the estimation is assessed by looking specifically at the value of σ t since it controls the CVs.

Table 7 gathers the results varying tol e , the anisotropy ratio and p. Regardless of the tolerance selected and the value of p, the best estimates are obtained for σ exact = (0.08, 0.02) cm 2 /ms, the PGD approximation being in such a case closer to the FE solution. The recovery of the parameter σ t is particularly precise and exhibits a low sensitivity to p, meaning that we can reconstruct the true propagation dynamics with a reliable prediction for the CVs. On the contrary, more variability affects σ l estimates. This is in accordance with the results in Section 4.3, where we noticed that the PGD solution depends more on the parameter σ t , which controls the CVs, while it is less sensitive to the value of σ l . As for the other conductivity pairs, the estimation obtained with tol e =10 -4 is less precise, especially for the high CV case, R a = 2, consistently with the fact that the PGD error is higher in this case, as shown in Fig. 2(c). An overall improvement of the results is evident for lower tolerance values. For tol e =10 -6 , the estimates are accurate and robust with respect to the presence of noise in the data.

In view of a better estimation in the case of high CVs, we restrict the admissible range for σ t , by choosing a length ∆ t = 0.01 and 0.02 cm 2 /ms, respectively. The same strategy has been used in Section 4.3 to reduce the error of the PGD approximation (see Fig. 5). Table 8 shows a more precise estimation when using ∆ t = 0.01 cm 2 /ms both for tol e =10 -4 and 10 -5 . This agrees with the results in Fig. 5 where, for ∆ t = 0.01 cm 2 /ms, the PGD method is able to capture fast propagating excitable waves. Regarding ∆ t = 0.02 cm 2 /ms, an improvement is visible only for tol e =10 -5 . Results are robust with respect to the percentage, p, of noise.

The extremely fast evaluation of the reduced solution in the PGD online phase makes inversion remarkably inexpensive. In fact, solving the MICP using the PGD approximation of the transmembrane potential requires, at most, only 30 s (see Table 9).

The results presented so far suggest some practical ideas on using the PGD for solving the MICP, even

with tol e = 10 -4 . In a sort of Predictor-Corrector approach, the PGD with a large tolerance can be used to inform a second-level computation. The latter can be performed either with PGD again but with a much narrowed exploration interval for the parameters (as suggested by strategy (c) in the previous Section) or with a fill-order solution, i.e. using the classical FE solver. The latter idea was used here in the case of σ exact = (0.07, 0.035) cm 2 /ms and tol e =10 -4 . Solving the MICP with a FE approximation and starting from the PGD estimate σ 0 = (0.0657, 0.0328) cm 2 /ms takes 527 s to reach convergence with the final estimation σ = (0.07, 0.0349) cm 2 /ms, to compare with 3293 s needed starting from σ 0 = (0.06, 0.025) cm 2 /ms to reach a similar estimate, precisely, σ = (0.0702, 0.0360) cm 2 /ms.

A 3D Canine Tissue Geometry

The MICP with PGD technique has been tested with the 3D canine geometry in Section 4.3. The case tol e =10 -6 has not been considered because the offline phase would be too expensive. As the PGD approximation for the propagation of the transmembrane potential in 3D yields results similar to the 2D case (see Section 4.3), we do expect results qualitatively comparable to the ones obtained in the 2D framework also for the MICP-PGD procedure. This is confirmed in Table 10. We obtain better results for σ exact = (0.08, 0.02) cm 2 /ms with high precision for the parameter σ t and a general improvement of the estimation when decreasing the tolerance value to 10 -5 . We highlight also for this test case the huge computational saving provided by the combination of MICP with PGD. In fact, despite the large number of DOFs of the 3D mesh (roughly 135K), the solution of the inverse problem only requires about 180 s when tol e = 10 -5 , as highlighted in Table 9.

p = 1%, [σ] = mm 2 /ms R a = 2 R a = 4 R a = 6 σ exact =(7.
p = 1%, [σ] = mm 2 /ms R a = 2 R a = 4 R a = 6 σ exact = (7.00, 3.50) σ exact = (8.00, 2.00) σ exact = (8. [START_REF] Keener | Mathematical Physiology[END_REF] Table 10: 3D canine tissue: conductivity estimation provided by the PGD-MICP procedure varying the external tolerance, the anisotropy ratio and the percentage of noise. For the sake of readability, we display the conductivities in mm 2 /ms.

A Real Ventricle

Finally, we analyze the estimation obtained in the real left ventricular test case. As previously discussed, the PGD approximation applied to this geometry is more sensitive to the anisotropy ratio, R a , and features higher errors (see Fig. 9). Therefore, the precision of the inversion may be affected. From Table 11, we note that using the PGD basis obtained with tol e =10 -4 leads to acceptable results only for σ exact = (0.07, 0.035) cm 2 /ms, whereas in the other tests the estimates it is fairly inaccurate, especially for σ l . The estimates become more accurate as the tolerance decreases and reasonably match the exact conductivities in the case σ exact = (0.07, 0.035) cm 2 /ms and σ exact = (0.08, 0.02) cm 2 /ms, whereas the estimation of σ l for σ exact = (0.084, 0.014) cm 2 /ms still lacks accuracy even for tol e =10 -6 . This agrees with the error pattern shown in Fig. 9, where, for these particular parameter values, the discrepancy between FE and PGD solutions increases when reducing the tolerance. A possible strategy, in this case, is to refine the space discretization for the PGD solver. Table 9 confirms the computational advantages led by solving the MICP with PGD, estimations being obtained in, at most ,40 s.

Conclusions and Perspectives

The patient-specific customization of mathematical models is a crucial step to bring numerical modeling into the clinical routine. Unfortunately, the high computational costs of standard data assimilation procedures conflict with clinical time frames. Ideally, one would like to have a nearly real-time estimate of the (3) the cost of the offline phase could be dumped if recyclable on different geometries. This calls for the construction of a PGD library of offline solutions on a reference geometry to be successively mapped onto a real patient-specific morphology. While this approach might slow down the online phase, the overall benefit for a large pool of patients and, eventually, for clinical applications could be potentially high. Finally, there is the need of a deep theoretical analysis of the convergence of the iterative solvers of the PGD, to identify rigorously the conditions that guarantee the reliability of the approach and, possibly, appropriate acceleration techniques.

Encouraged by the extremely positive results reached in this work, we plan to pursue these developments in the next future.

  linear forms, for ϕ = 1, . . . , N F , where N a , N F ∈ N + and indices ψ, β = 1, . . . , m strictly depend on the differential operators involved in the definition of A and F, the test function v being rewritten

  with Q = Ω × (0, T ) the space-time domain with boundary ∂Q = ∂Ω × [0, T ]; [0, T ] a fixed time interval; u the transmembrane potential ([mV]); w the gating variable characterizing the RM model; σ the conductivity tensor ([cm 2 /ms]), in general function of the local coordinates, x, and of the local fiber direction, a; I ion the total ionic current flowing through the cell membrane; g the kinetic dynamics of the ionic quantities; I app the external stimulation current (e.g., electrical pacing or synaptic input); n the outward unit normal vector to the boundary ∂Ω; U 0 and w 0 the initial value for the potential and for the gating variable, respectively.

1 + mu 2 functions, the cubic term leads to mu 1 + 2 mu 2 + mu 3 functions

 12123 and the mixed term, in u n and w n+1 , m u m w functions. Finally, we remark that while N a (equal to 10 in (29)) only depends on the differential terms identifying the bilinear form A and remains unchanged during the PGD iterations and in time, N F changes during the enrichment iterations and, additionally, it changes in time, due to the time dependence in F σ .

  the contributions determined at the m-th ADS iteration are initially set to u n+1,Neumann-type boundary condition (13) 2 . As for the tolerances tol e and tol i in Algorithm 1, we will select different values to investigate the impact of this choice on the performance of the PGD.

Figure 1 :

 1 Figure1: (a) Realistic geometry of a 2D portion of a canine ventricular tissue: computational mesh (left) and approximation of the fiber structure (right). The segments represent the local cardiac fiber direction, while the colorbar indicates the value of the local fiber angle with respect to the x-axis (unit of measurement is in radians). (b) 3D mesh (left) and myocardial fiber orientation from two different viewpoints (center-right) for the real ventricle simulation (image from[START_REF] Yang | Estimation of cardiac conductivities in ventricular tissue by a variational approach[END_REF]).

Figure 2 :

 2 Figure 2: 2D canine tissue: (a) Trend of the PGD modes for different tolerances. (b) Percentage L 2 -relative error between the FE and the PGD potential, for different conductivity values and tole = 10 -4 . (c) Percentage L 2 -relative error between the FE and the PGD potential for different tolerances and anisotropy ratios.

  e =10 -4 PGD tol e =10 -5PGD tol e =10 -6

Figure 3 : 6 Figure 4 :

 364 Figure 3: 2D canine tissue: FE and PGD solutions for different anisotropy ratios and tolerances. The green arrow in (a) points out the stimulation site. (a) Ra = 2, σ = (0.07, 0.035) cm 2 /ms, t = 17 ms. (b) Ra = 4, σ = (0.08, 0.02) cm 2 /ms, t = 22 ms. (c) Ra = 6, σ = (0.084, 0.014) cm 2 /ms, t = 25 ms.

Figure 5 :

 5 Figure 5: 2D canine tissue: possible improvements of PGD approximation with tole=10 -4 in case of high CVs (σ = (0.07, 0.035) cm 2 /ms). (a) Percentage L 2 relative error on the potential when increasing the number of DOFs for the discretization of S versus when lowering tolerance tol i . (b) Percentage L 2 -relative error on the potential (left panel) and trend of the number of modes (right panel) when narrowing the admissible interval for σt.

R a = 2 ,Figure 6 :

 26 Figure 6: 3D canine tissue: percentage L 2 -relative error between the FE and the PGD potential for different tolerances and anisotropy ratios.

  e = 10 -4 PGD tol e = 10 -5Abs. error tol e = 10 -5

Figure 7 :

 7 Figure 7: 3D canine tissue: FE and PGD solutions for different anisotropy ratios and tolerances (three left panels). The green arrow in (a) points out the stimulation site. (a) Ra = 2, σ = (0.07, 0.035) cm 2 /ms, t = 17 ms. (b) Ra = 4, σ = (0.08, 0.02) cm 2 /ms, t = 22 ms. (c) Ra = 6, σ = (0.084, 0.014) cm 2 /ms, t = 25 ms. The right panel shows the distribution of the absolute value of the error on the potential [mV] between FE and PGD solutions.

Figure 8 :

 8 Figure 8: Real ventricle: snapshots of the excitation wave propagation (visualized from two different perspectives) simulated using the FE method. The white arrow in snapshot t=5 ms, (top) highlights the stimulation site triggering the propagation.

R a = 2 ,

 2 σ = (0.07, 0.035) R a = 4, σ = (0.08, 0.02) R a = 6, σ = (0.084, 0.014)

Figure 9 :

 9 Figure 9: Real ventricle: percentage L 2 -relative error between the FE and the PGD potential for different tolerances and anisotropy ratios.

FE

  PGD tol e = 10 -4 PGD tol e = 10 -5 PGD tol e = 10 -6 Abs. error tol e = 10 -6

Figure 10 :

 10 Figure 10: Real ventricle: FE and PGD solutions for different anisotropy ratios and tolerances (four left panels) at t = 20 ms. The stimulation point is located at the ventricular endocardium, thus it is not visible from this visual perspective. (a) Ra = 2, σ = (0.07, 0.035) cm 2 /ms. (b) Ra = 4, σ = (0.08, 0.02) cm 2 /ms. (c) Ra = 6, σ = (0.084, 0.014) cm 2 /ms. The right panel shows the distribution of the absolute value of the error on the potential [mV] between FE and PGD solutions.

Figure A. 11 :

 11 Figure A.11: 2D canine tissue geometry: (a) comparison between the FE and the PGD solutions at the mesh nodes featuring the highest approximation error. (b) L 2 -norm of the absolute error of the PGD approximation for the potential with respect to the FE solution, for σ = (0.08, 0.02) cm 2 /ms and tole=10 -4 .

Figure A. 12 :

 12 Figure A.12: 2D canine tissue geometry: comparison between the FE and the PGD approximations for the wavefront (a) and waveback (b) propagation [mV], and associated L 2 -norm of the absolute error for σ = (0.08, 0.02) cm 2 /ms and tole=10 -4 .

1 :

 1 Input: u 0 , S, tol e , tol i , m max , j max 2: Set: m ← 0;3: while ( u m / u 1 ≤ tol e & m ≤ m max ) do

	4:

Table 2 :

 2 Factorization of the PGD extended bilinear form in

Table 4 :

 4 10 -4 10 -5 10 -6 2D canine tissue: CPU times of the PGD approach for different values of tole, and time demanded by the FE solver.

	FE

Table 5 :

 5 3D canine tissue: CPU times of the PGD approach for different values of tole, and time demanded by the FE solver.

		10 -4 10 -5	FE
	Offline [h]	8	41	/
	Online [s]	0.6	0.75	330

Table 6 :

 6 Real ventricle: CPU times of the PGD approach for different values of tole, and time demanded by the FE solver.

		10 -4 10 -5 10 -6	FE
	Offline [h]	0.8	4.1	32	/
	Online [s]	0.1	0.15	0.2	60

  R a = 6 σ exact = (7.00, 3.50) σ exact = (8.00, 2.00) σ exact = (8.40, 1.40) tol e =10 -4

		00, 3.50)	σ exact =(8.00, 2.00)	σ exact =(8.40, 1.40)
	tol e =10 -4	(6.57,3.28)	(8.11,1.93)	(8.98,1.29)
	tol e =10 -5	(6.93,3.96)	(7.71,1.96)	(8.57,1.34)
	tol e =10 -6	(6.90,3.56)	(7.80,1.96)	(8.22,1.35)
		p = 5%, [σ] = mm 2 /ms
		R a = 2	R a = 4	R a = 6
		σ exact = (7.00, 3.50) σ exact = (8.00, 2.00) σ exact = (8.40, 1.40)
	tol e =10 -4	(6.56,3.27)	(7.91,1.93)	(8.97,1.31)
	tol e =10 -5	(6.95,3.96)	(7.76,1.95)	(8.78,1.35)
	tol e =10 -6	(6.83,3.56)	(7.84,1.97)	(8.31,1.35)
		p = 10%, [σ] = mm 2 /ms
		R a = 2	R a = 4
		(6.55,3.24)	(8.05,1.93)	(8.98,1.30)
	tol e =10 -5	(6.96,3.97)	(7.73,1.95)	(8.44,1.36)
	tol e =10 -6	(6.82,3.56)	(7.91,1.96)	(8.26,1.35)

Table 7 :

 7 2D canine tissue: conductivity estimation provided by the PGD-MICP procedure varying the external tolerance, the anisotropy ratio and the percentage of noise. For the sake of readability, we display the conductivities in mm 2 /ms. =0.03 (6.57,3.28) (6.56,3.27) (6.[START_REF] Luo | A model of the ventricular cardiac action potential. depolarization, repolarization, and their interaction[END_REF],3.24) ∆ t =0.02 (7.85,3.71) (7.86,3.73) (7.84,3.71) ∆ t =0.01 (6.91,3.33) (6.90,3.32) (6.91,3.35) tol e = 10 -5 , [σ] = mm 2 /ms p = 1% p = 5% p = 10% ∆ t =0.03 (6.93,3.96) (6.95,3.96) (6.96,3.97) ∆ t =0.02 (6.40,3.66) (6.39,3.66) (6.40,3.65) ∆ t =0.01 (6.80,3.35) (6.75,3.35) (6.82,3.34)

	tol e = 10 -4 , [σ] = mm 2 /ms
	p = 1%	p = 5%	p = 10%
	∆ t

Table 8 :

 8 2D canine tissue: estimation of the conductivity σexact = (0.07, 0.035) provided by the MICP-PGD procedure varying the external tolerance, the length of the admissible domain ∆t for σt and the percentage of noise. For the sake of readability, we display the conductivities in mm 2 /ms.

		2D canine tissue 3D canine tissue Real ventricle
	tol e =10 -4	22	154	32
	tol e =10 -5	25	180	36
	tol e =10 -6	30	/	40

Table 9 :

 9 Execution time [s] associated with the MICP-PGD approach varying the external tolerance. The times are average values since they may vary depending on the realization of the noise random variable and on the initial guess of the optimization procedure.

  R a = 6 σ exact = (7.00, 3.50) σ exact = (8.00, 2.00) σ exact = (8.40, 1.40) tol e =10 -4 R a = 6 σ exact = (7.00, 3.50) σ exact = (8.00, 2.00) σ exact = (8.40, 1.40) tol e =10 -4

				, 1.40)
	tol e =10 -4	(6.78,3.21)	(8.07,1.98)	(8.98,1.32)
	tol e =10 -5	(6.97,3.98)	(7.82,1.99)	(8.38,1.37)
			p = 5%, [σ] = mm 2 /ms
		R a = 2	R a = 4
		(6.90,3.20)	(8.11,1.97)	(8.96,1.31)
	tol e =10 -5	(6.98,4.00)	(8.56,1.99)	(8.38,1.39)
		p = 10%, [σ] = mm 2 /ms
		R a = 2	R a = 4
		(6.81,3.25)	(8.13,1.96)	(9.00,1.32)
	tol e =10 -5	(6.97,4.00)	(7.92,1.99)	(8.45,1.39)

Acknowledgements

AB and AV acknowledge the support of the NSF Project DMS 1412963 and the XSEDE Consortium, Grant TG-ASC160069. MGC wants to thank the Department of Mathematics at Emory University for hosting him during part of the research and the support of the Politecnico di Milano for the travel expenses. SP acknowledges the INdAM-GNCS Project 2019 Advanced intrusive and non-intrusive model order reduction techniques and applications. AG acknowledges the support of the Italian National Group for Mathematical Physics (GNFM-INdAM) Project Constitutive modeling of active fiber-distributed media.

R a = 4 R a = 6 σ exact = (7.00, 3.50) σ exact = (8.00, 2.00) σ exact = (8. [START_REF] Keener | Mathematical Physiology[END_REF] patient-specific parameters. Specific model reduction methods that may retain the clinical accuracy of the full mode (i.e., the accuracy needed for diagnosing or decision-making), with a significant improvement of the computational efficiency are crucial, for instance, in uncertainty quantification and parameter estimation (see, e.g., [START_REF] Brandes | Isogeometric hierarchical model reduction: from analysis to patientspecific simulations[END_REF][START_REF] Alvarez | Transversally enriched pipe element method (tepem): An effective numerical approach for blood flow modeling[END_REF][START_REF] Guzzetti | Propagating uncertainties in large-scale hemodynamics models via network uncertainty quantification and reduced-order modeling[END_REF][START_REF] Fleeter | Multilevel and multifidelity uncertainty quantification for cardiovascular hemodynamics[END_REF][START_REF] Seo | Multi-fidelity estimators for coronary circulation models under clinically-informed data uncertainty[END_REF]).

In the specific field of electrophysiology, the need of an accurate estimation of the conductivities in the Monodomain problem is associated with many possible applications, for instance when dealing with the optimal placement of a pace-maker or the identification of ablation sites [START_REF] Kunisch | Optimal control of the bidomain system (i): The monodomain approximation with the rogers-mcculloch model[END_REF][START_REF] Kunisch | Optimal control of the bidomain system (ii): Uniqueness and regularity theorems for weak solutions[END_REF][START_REF] Kunisch | Optimal control of the bidomain system (iii): Existence of minimizers and firstorder optimality conditions[END_REF][START_REF] Nagaiah | Optimal control approach to termination of re-entry waves in cardiac electrophysiology[END_REF][START_REF] Petras | Tissue drives lesion: computational evidence of interspecies variability in cardiac radiofrequency ablation[END_REF]. Parameter estimation calls for accurate and efficient model reduction techniques. However, consolidated model reduction techniques, like the POD, may suffer from a non optimal selection of the snapshots. In this paper, a snapshotfree model order reduction technique like PGD, thus circumventing one of the main drawbacks of POD. In spite of a significant offline cost, the minimal cost of the online phase makes PGD very competitive with respect to other model reduction techniques, in particular, when involved in multi-query contexts, such as the resolution of inverse problems.

What presented here is just a preliminary step of a more complex series of possible developments. Among these: (1) the extension to more complex ionic models (and, eventually, to the Bidomain model) is not trivial, as the modeling of the ionic terms requires a specific factorization of the increased linear functional F; (2) the rapid solution of the online phase makes affordable the introduction of Uncertainty Quantification techniques, like the Bayesian ones. This is a critical step for the self-assessment of the quality of the parameter estimation;

Appendix A. PGD Approximation of the Whole Action Potential

The simulations considered in the body of this work only focused on the wavefront propagation to limit computational costs related to the offline phase. Here we investigate the accuracy of the PGD method to reconstruct the waveback propagation. For this purpose, we set T=150 ms to capture the whole Action Potential (AP). We refer to the 2D canine tissue geometry with σ = (0.08, 0.02) cm 2 /ms and tol e =10 -4 .

Fig. A.11 (a)

compares the FE and the PGD transmembrane potentials at the mesh nodes featuring the highest approximation error. Error is negligible in the resting phase (in this case, roughly in the interval [0,20] ms and after 120 ms), when u = 0 mV, and in the plateau phase ( [START_REF] Abidi | Ionic parameters estimation in multi-scale cardiac electrophysiology modelling[END_REF]115] ms), when the transmembrane potential is slowly decreasing. On the contrary, the PGD approximation is less accurate when the transmembrane potential sharply changes. In particular, the error is more significant in the interval [START_REF] Loppini | Spatiotemporal correlation uncovers characteristic lengths in cardiac tissue[END_REF][START_REF] Abidi | Ionic parameters estimation in multi-scale cardiac electrophysiology modelling[END_REF] ms characterizing the wavefront propagation than in the interval [115,120] Extending the PGD technique to more complex phenomenological models, such as the Mitchell-Schaeffer [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF], the Fenton-Karma [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation[END_REF] and the Minimal [START_REF] Bueno-Orovio | Minimal model for human ventricular action potentials in tissue[END_REF] ones, that provide a more accurate reproduction of the cardiac AP, may improve the quality of the waveback approximation. As expected, a larger simulation time results in a heavier computational demand during the offline stage. In fact, Algorithm 1 requires 12 hours as opposed to 0.9 hours for the shorter simulation with T = 30.5 ms (see Table 4).


Yannick De Oliveira
email: ydeoliveira@moulon.inra.fr

Julien Cornouiller

Guy-Ross Assoumou-Ella

Johann Joets

Alain Charcosset

Thaliadb A database dedicated to association genetics in plants

Introduction

Diversity and association genetics studies lead to manipulate a large number of individual, lines, clones and/or populations. Moreover, emergence of highthroughput technologies for both genotyping and phenotyping generates a large amount of data. These need to be stored and managed in order to perform requests and organize datasets to conduct association genetics studies. The Thaliadb database has been developed with this aim. It manages genetic resources, phenotyping and genotyping data, and also population structure information. Thaliadb enables data extraction in formats used by genetic association softwares. Seed lot history management project This project aim to develop a database to store and manage seed lot history data (relations between seed lots, stock etc.)

Data Structure

The web service we are developing will allow request on phenotyping and genotyping data.

GnpAsso project (PI D. Steinbach)

The aim of GnpAsso project [5] is to develop generic tools to manage and exploit results from association genetics studies using genotyping and phenotyping data.

The goal is to store significant associations in the GnpAsso database.

Thaliadb is integrated with SniPlay [6] pipeline to process datasets.

Thaliadb enables data extraction and visualization at multiple level :

-A view of a single entity (population, loci,...; picture on the left) with all data related to this entity -A table view (on the top) describing several entities. Thaliadb allows user to request and extract data in format compatible with association genetic software.

-On the right, a map view (Google map) with the geographic repartition of accession genotypes. This example shows the frequencies of the allele B of VGT1-MITE locus for accessions located in Europe.

[5] http://www.agence-nationale-recherche.fr/projet-anr/?tx_lwmsuivibilan_pi2[CODE]=ANR-10-GENM-0006

[6] A. Dereeper, S. Nicolas, L. Le Cunff, R. Bacilieri, A. Doligez, J-P Peros, M Ruiz, P This, SNiPlay: a webbased tool for detection, management and analysis of SNPs. Application to grapevine diversity projects. BMC Bioinformatics 12:134, 2011.

Funding : Thaliadb is developed with funding from ANR project ANR-10-GENM-0006 (PI D. Steinbach)