Yassine Abidi

Moncef Mahjoub

N Éjib Zemzemi

IONIC PARAMETERS ESTIMATION IN MULTI-SCALE CARDIAC ELECTROPHYSIOLOGY MODELLING

Keywords: Optimal control with PDE constraints, bidomain model, physiological ionic model, maximal conductance parameters

In this work, we present an optimal control formulation for the bidomain model in order to estimate maximal conductance parameters in cardiac electrophysiology multiscale modelling. We consider a general Hodgkin-Huxley formalism to describe the ionic exchanges at the microscopic level. We treat the desired parameters as control variables in a cost function minimizing the gap between the measured and the computed transmembrane potentials. First, we establish the existence of an optimal control solution and we formally derive the optimality system. Second, we propose a strategy for solving the estimation problem for both single and multiple parameters cases. Our algorithm is based on a gradient descent method, where the gradient is obtained by solving an adjoint problem. Both the state and the adjoint problems are solved using the finite element method. Numerical simulations for single and multiple conductances estimations show the capability of this approach to identify the values of sodium, calcium and potassium ion channels conductances of the Luo Rudy phase I model.

INTRODUCTION

In the last decades, the mathematical modelling of cardiac electrical activity has been recognized as one of the potential approaches capable of revealing diagnostic information about the heart. The electrical behavior of the cardiac tissue is described by a system consisting of partial differential equations (PDEs) coupled to a system of ordinary differential equations (ODEs) modelling the transmbrane ionic exchanges at the microscopic level. The model of the cell membrane dynamics is expressed using the Hodgkin-Huxley (HH) formalism [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. This model was adapted by Denis Nobel for Purkinje heart cells [START_REF] Noble | Cardiac action and pacemaker potentials based on the hodgkin-huxley equations[END_REF]. Since that, several other models have been introduced to describe the electrical activity of the cell membrane in the myocardium. Beeler and Reuter [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF] introduced the model for ventricular cells. Later, Di Francesco and Noble [START_REF] Francesco | A model of cardiac electrical activity incorporating ionic pumps and concentration changes[END_REF] proposed a model that considers ion pumps, which allows different chemical species such as potassium, sodium and calcium to return to their stable state. Moreover, a family of more complex models based always on the HH formalism taking into account the physiological behavior of the ion channels [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] Hund | Rate dependence and regulation of action potential and calcium transient in a canine cardiac ventricular cell model[END_REF][START_REF] Tusscher | Cell model for efficient simulation of wave propagation in human ventricular tissue under normal and pathological conditions[END_REF][START_REF] Grandi | A novel computational model of the human ventricular action potential and Ca transient[END_REF] have been proposed in the litterature. In all these models, maximal conductances play an important role in the generation of the electrical potential, the determination of some pathological conditions or in the study of the effect of drugs. Several authors investigated the optimization of parameters or source terms in the mathematical model of cardiac electrophysiology. Brandao et al. [START_REF] Brandao | Theoretical analysis and control results for the fitzhugh-nagumo equation[END_REF] studied the theoretical analysis and the controllability of the optimization subject to the FitzHugh-Nagumo model. Later, systematic analysis and numerical studies for the optimal control of monodomain and bidomain model are presented in [START_REF] Chamakuri | Numerical solution for optimal control of the reaction-diffusion equations in cardiac electrophysiology[END_REF][START_REF] Ainseba | Analysis of an optimal control problem for the tridomain model in cardiac electrophysiology[END_REF][START_REF] Casas | Sparse optimal control of the schlögl and fitzhugh-nagumo systems[END_REF][START_REF] Kunisch | Optimal control of the bidomain system (iii): Existence of minimizers and first-order optimality conditions[END_REF][START_REF] Chamakuri | Optimal control approach to termination of re-entry waves in cardiac electrophysiology[END_REF][START_REF] Andreianov | Solvability analysis and numerical approximation of linearized cardiac electromechanics[END_REF][START_REF] Chamakuri | Application of optimal control to the cardiac defibrillation problem using a physiological model of cellular dynamics[END_REF], and more recently for the optimal control of bidomain-bath model using Mitchell-Shaeffer model in 3D geometries [START_REF] Chamakuri | Pde constrained optimization of electrical defibrillation in a 3d ventricular slice geometry[END_REF][START_REF] Bendahmane | A 3d boundary optimal control for the bidomain-bath system modeling the thoracic shock therapy for cardiac defibrillation[END_REF][START_REF] Chamakuri | Primal-dual active set strategy for large scale optimization of cardiac defibrillation[END_REF]. In [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] authors established an experimental estimation of the intracellular and extracellular conductivities. Yan and Veneziani [START_REF] Yang | Estimation of cardiac conductivities in ventricular tissue by a variational approach[END_REF] use a variational procedure for the estimation of cardiac conductivities from measures of the transmembrane and extracellular potentials available at some sites of the tissue. Beretta et al. [START_REF] Beretta | An inverse problem for a semilinear parabolic equation arising from cardiac electrophysiology[END_REF] define a numerical reconstruction procedure for the inverse problem of detecting a spherical inhomogeneity from boundary measurements of the electric potential. Moreover, the identification from measurements of surface potentials has been tackled in an optimization framework for numerical purposes [START_REF] Álvarez | Shape reconstruction of cardiac ischemia from non-contact intracardiac recordings: A model study[END_REF][START_REF] Chávez | Inverse problem of electrocardiography: Estimating the location of cardiac ischemia in a 3d realistic geometry[END_REF][START_REF] Nielsen | On the use of the resting potential and level set methods for identifying ischemic heart disease: An inverse problem[END_REF]. In an other work [START_REF] Ngoma | Parameter identification for a non-differentiable ionic model used in cardiac electrophysiology[END_REF], authors propose a strategy to optimize a non differentiable cost function obtained from a fit of activation times map. Recently, the identification of the maximal conductances has been subject of theoretical studies in the monodomain [START_REF] Abidi | On the identification of multiple space dependent ionic parameters in cardiac electrophysiology modelling[END_REF][START_REF] Boulakia | Stability estimates for some parameters of a reaction-diffusion equation coupled with an ode[END_REF][START_REF] Ainseba | Stability of conductivities in an inverse problem in the reaction-diffusion system in electrocardiology[END_REF][START_REF] Lassoued | Stability results for the parameter identification inverse problem in cardiac electrophysiology[END_REF] and bidomain [START_REF] Wu | Carleman estimate for a linearized bidomain model in electrocardiology and its applications[END_REF][START_REF] Abidi | Ionic parameters identification of an inverse problem of strongly coupled pdes system in cardiac electrophysiology using carleman estimates[END_REF] models. In an other work [START_REF] Bouyssier | Parameters estimation approach for the mea/hipsc-cm asaays[END_REF], authors propose a numerical approach for the analysis of drugs effects on the electrical activity of hiPSC-CMs based on multi-electrode array experiments where the drug acts directly on the maximal conductance of the targeted ion channel. In our work, we still use a variational procedure for the estimation of a set of ionic maximal conductances by an optimal control approach. This approach is based on the minimization of an appropriate cost functional that depends on the maximal conductances and a measures of the transmembrane potential.

The paper is organized as follows: In Section 2, we briefly recall the bidomain model and the general structure of cardiac cellular membrane models describing the electrical wave propagation and the ionic exchange at the cell membrane respectively. We also present some preliminary material, including relevant notations, assumptions and regularities. Section 3 is devoted to the optimal control problem. We prove the existence of the control and we provide a formal derivation of the adjoint equations and the first order optimality condition, which are the basis for numerical resolution. The description of the numerical approach to solve the primal and the adjoint state equations and the optimization problem are explained in section 4. Finally, in Section 5, numerical results are presented for the different test cases. A summary and concluding remarks are given then in the last section.

MATHEMATICAL MODEL

Let Ω Ă R d pd ě 1q be a bounded connected open set whose boundary Γ " BΩ is regular enough, (Ω Ă R 3 being the natural domain of the hearth). Let T ą 0 be a fixed time horizon. We will use the notation Q " Ω ˆp0, T q, and Σ " Γ ˆp0, T q.

We introduce a parabolic-elliptic system called bidomain model, coupled to a system of ODEs. This model was proposed in the late 1970s by Tung [START_REF] Tung | A bidomain model for describing ischemic myocardial d-c potentials[END_REF] and is now the generally accepted model of electrical behaviour of cardiac tissue (see Henriquez [START_REF] Henriquez | Simulating the electrical behavior of cardiac tissue using the bidomain model[END_REF], Keener and Sneyd [START_REF] Keener | Mathematical Physiology, I: Cellular Physiology[END_REF]), can be written as:

(2.1) $ ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' % A m `Cm B t v `Iion p ¯ , v, w, zq ˘´divpσ i ∇vq " divpσ i ∇u e q `Am I app in Q, ´divpσ i ∇v `pσ i `σe q∇u e q " 0 in Q, B t w " F pv, wq in Q, B t z " Gp ¯ , v, w, zq in Q, σ i ∇v.ν `σi ∇u e .ν " 0 on Σ, σ i ∇v.ν `pσ i `σe q∇u e .ν " 0 on Σ, vpx, 0q " v 0 pxq, wpx, 0q " w 0 pxq, zpx, 0q " z 0 pxq in Ω,
where v : Q Ñ R is the transmembrane potential, u e : Q Ñ R is the extracellular electric potential, and σ i , σ e : Ω Ñ R dˆd are respectively the intra-and extracellular conductivity tensors. w : Q Ñ R k represent the gating variables and z : Q Ñ R m are the ionic intracellular concentration variables. A m is the surface to volume ratio of the cardiac cells, and C m ą 0 is the membrane capacitance per unit area. I app : Q Ñ R is the applied current source and ¯ :" t¯ i u 1ďiďN represent a set of maximal conductance parameters. The I ion is the current density flowing through the ionic channels and the functions F and G determine the evolution of the gating variables and intracellular concentrations, respectively, they are determined by an electrophysiological cell model. In an isolated heart conditions, no current flows out of the heart as expressed by the homogeneous Neumann boundary conditions. We denote by ν the unit normal to Γ outward of Ω.

2.1. Membrane models and ionic currents. Following work by Hodgkin and Huxley [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF], many HHbased models have later been developed for the cardiac action potential. In these models, the ionic current I ion through channels of the membrane, has the following general structure [START_REF] Franzone | Mathematical Cardiac Electrophysiology[END_REF]:

(2.2)

I ion p ¯ , v, w, zq " N ÿ i"1 ¯ i y i pvq k ź j"1 w p j,i j pv ´Ei pzqq,
where N is the number of ionic currents, ¯ i :" ¯ i pxq is the maximal conductance associated with the i th current, y i is a gating function depending only on the membrane potentiel v, p j,i are positive integers exponents and E i is the reversal potential for the i th current, which is given by Godman-Hodgkin-Katz (GHK) equation:

(2.3) E i " RT F ln ¨řj P R X j rX j s out `řk P R X ḱ rX ḱ s in ř j P R X j rX j s in `řk P R X ḱ rX ḱ s out ',
where R, T and F designate, respectively, the perfect gas constant, the temperature and the Faraday constant. P R X j,k represents the permeability for the ion X j,k , and rX j,k s in (respectively rX j,k s out ) is the intracellular (respectively extracellular) concentration of the ion X j,k . We recall that GHK equation is a generalized form of the Nernst equation.

For the ODEs, the dynamics of the gating variable w is described in the HH formalism by a system of ordinary differential equations which when w j is a gating variable (0 ď w j ď 1) are governed by the following equation:

(2.4) B t w j " F j pv, w j q :" α j pvqp1 ´wj q ´βj pvqw j , j " 1, . . . , k, where α j and β j two positive rational functions of exponentials in v. A general expression for both α j and β j is given by (2.5) C 1 exprC 2 pv `C3 qs `C4 pv `C5 q exprC 5 pv `C3 qs `C7 ,

where C 1 , C 2 , . . . , C 7 are constants. For the dynamics of the ionic concentration variables z, we have the system of ODE's:

(2.6)

B t z i " G i p ¯ , v, w, zq :" ´Ji p ¯ , v, w, log z i q `Hi p ¯ , v, w, zq, i " 1, . . . , m.
where the functions J i and H i are described in the Assumptions 2.1.

Conductivity tensors.

The anisotropic properties of the two media are modeled by an intracellular and extracellular conductivity tensors σ i pxq and σ e pxq. Generally, the conductivities σ i and σ e are two matrices given by σ j pxq " σ t j pxqI `pσ l j pxq ´σt j pxqqa l pxqa T l pxq, where σ l j and σ t j , j P ti, eu denote the intra-and extracellular conductivities along and transversal to the direction of the fiber (parallel to a l pxq), respectively, and I is the identity matrix.

Assumptions 2.1. (Stronger assumptions on the data in 2.1) We assume that:

(1) Ω Ă R 3 is a bounded domain with C 1,1 -boundary.

(2) σ i , σ e : Ω Ñ R 3ˆ3 are symmetric, positive definite matrix functions with W 1,8 pΩq-coefficients, obeying the uniform ellipticity conditions:

(2.7) 0 ď µ 1 }ξ} 2 ď ξ T σ i,e pxqξ ď µ 2 }ξ} 2 , @ξ P R 3 , @x P Ω, with µ 1 , µ 2 ą 0.

(3) We use the regularized form of the variable y i pvq in hyperbolic functions introduced in [START_REF] Endresen | A theory for the membrane potential of living cells[END_REF]. In this case, y i pvq is a C 8 function with respect to the variable v for i " 1, . . . , N , and then is locally Lipschitz since v is bounded, similarly for the function v Þ Ñ y i pvqv. (4) J i and H i are locally Lipschitz continuous functions where:

(2.8)

J i P C 2 pR ˚ˆRˆR k ˆRq, 0 ă g ˚pwq ď BJ i Bτ p¯ , v, w, τ q ď g ˚pwq, ˇˇˇB J i Bv p¯ , v, w, 0q ˇˇˇď L v pwq, g ˚, g ˚, L v belong to C 1 pR k , R `q, and
(2.9) H i P C 2 pR ˚ˆR ˆRk ˆp0, `8q m q X LippR ˚ˆR ˆr0, 1s k ˆp0, `8q m q.

(5) The intial values belong to the following spaces: v 0 P H 2 pΩq, w 0 P L 2 pΩq k , and z 0 P L 2 pΩq m , with log z 0 :" plog z 0,1 , . . . , log z 0,m q P L 2 pΩq m . Moreover, v 0 satisfies the compatibility condition described in [START_REF] Veneroni | Reaction-diffusion systems for the macroscopic bidomain model of the cardiac electric field[END_REF].

(6) The extracellular potential u e has zero average on Ω, i.e ż Ω u e dx " 0.

(7) I app verifies the compatibility condition:

ż Ω I app px, tqdx " 0, @t P p0, T q.

Under the assumptions (2.1), the existence and uniqueness of the bidomain model 2.1 for more general physiological ionic models based on the classical Hodgkin-Huxley formalism can be found in [START_REF] Veneroni | Reaction-diffusion systems for the macroscopic bidomain model of the cardiac electric field[END_REF].

In our numerical computations, we consider the physiological Luo-Rudy phase I model (LR1) [START_REF] Luo | A model of the ventricular cardiac action potential: Depolarization, repolarization, and their interaction[END_REF], which consists of 6 ionic currents:

(2.10) I ion " I N a `Isi `IK `IK1 `IKp `Ib , which are fast sodium current (I N a ), slow inward calcium current (I si ), time dependent potassium current (I K ), time independent potassium current (I K1 ), plateau potassium current (I Kp ) and background current (I b ). The time dependent currents I N a , I si and I K , depend on six activation and inactivation gates m, h, j, d, f, x, and one intracellular concentration variable of Calcium rCa 2`s i , which are governed by ODEs of the form:

(2.11)

dw dt " α w pvqp1 ´wq ´βw pvqw, for w " m, h, j, d, f, x, d dt rCa 2`s i " ´10 ´4I si `0.07p10 ´4 ´rCa 2`s i q,
For details on formulation of those functions and the parameters used in our computations, we refer to the original paper of LR1 model [START_REF] Luo | A model of the ventricular cardiac action potential: Depolarization, repolarization, and their interaction[END_REF].

OPTIMAL CONTROL PROBLEM

In this section, we set the optimal control problem, for which the numerical experiments were carried out. Suppose that v meas is the desired state solution at the cardiac domain, we look for the set of parameters ¯ :" t¯ i u 1ďiďN that solves the following minimization problem:

(3.1) pPq $ ' & ' % min ¯ PC ad Ip ¯ q " 1 2 ` 1 ż Q |vp ¯ q ´vmeas | 2 dxt ` 2 ż Ω | ¯ | 2 dx ˘,
subject to the coupled PDE system (2.1), and ¯ P C ad , where I is the quantity of interest, 1 , 2 are the regularization parameters and v is the state variable. In addition, a Tikhonov-like regularization term ż Ω | ¯ | 2 dx used to weigh the impact of the regularization in the minimize procedure. C ad is the admissible domain for control given by (3.2) C ad " t ¯ P L 8 pΩq N : ¯ pxq P rm, M s N , @x P Ωu.

3.1. Optimal conditions and dual problem. In this subsection, we formally derive the optimality system associated to (3.1). Let's denote by J the cost function

J pv, ¯ q " 1 2 ` 1 ż Q |v ´vmeas | 2 dxt ` 2 ż Ω | ¯ | 2 dx ˘.
If vp ¯ q is a solution of (2.1), then we immediately have J pvp ¯ q, ¯ q " Ip ¯ q. We follow a Lagrangian approach and introduce the Lagrange functional:

(3.3)
Lpv, u e , w, z, ¯ , λ ˚q " J pv, ¯ q ´λ˚S pv, u e , w, z, ¯ q, where λ ˚:" pp, q, r, sqpx, tq denote the Lagrange multipliers, and Spv, u e , w, z, ¯ q is the state equation. Then,

Lpv, u e , w, z, ¯ , λ ˚q " J pv, ¯ q ´żQ pA m `Cm B t v `Iion p ¯ , v, w, zq ´Iapp ˘dxdt ´żQ p `´divpσ i ∇vq ´divpσ i ∇u e q ˘dxdt ´żQ q `´divpσ i ∇v `pσ i `σe q∇u e q ˘dxdt ´żQ r. `Bt w ´F pv, wq ˘dxdt ´żQ s. `Bt z ´Gp ¯ , v, w, zq ˘dxdt,

The first order optimality system is given by the Karusch-Kuhn-Tucker (KKT) conditions which result from equating the partial derivatives of L with respect to v, u e , w and z equal to zero. First, invoking integration by parts with respect to the temporal variable and Green's formula with respect to the spatial variable, we observe that Herein, we impose the terminal conditions (3.9) ppx, T q " 0, rpx, T q " 0, and spx, T q " 0 in Ω,

(3.5) ă BL Bv , δv ą "ă BJ Bv , δv ą ´λ˚ă BS Bv , δv ą "ă BJ Bv , δv ą `żQ tA m pC m B t p ´pB v I ion q `r.B v F `s.B v Guδvdxdt `żQ tdivpσ i ∇pq `divpσ i ∇qquδvdxdt ´żΣ tσ i ∇p.ν `σi ∇q.νuδvdSdt, (3.6)
Collecting the previous results with the derivatives of J , we get the following adjoint problem:

(3.10) $ ' ' ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' ' ' % ´Am pC m B t p ´pB v I ion q ´divpσ i ∇pq ´divpσ i ∇qq ´pB v F q T r ´pB v Gq T s " 1 pv ´vmeas q, in Q, ´divpσ i ∇p `pσ i `σe q∇qq " 0 in Q, ´Bt r `Am pB w I ion ´pB w F q T r ´pB w Gq T s " 0 in Q, ´Bt s `Am pB z I ion ´pB z Gq T s " 0 in Q,
ppx, T q " 0, rpx, T q " 0, and spx, T q " 0, in Ω, with the boundary conditions for the adjoint states

(3.11) $ & % ´σi ∇p.ν ´σi ∇q.ν " 0, on Σ, ´σi ∇p.ν ´pσ i `σe q∇q.ν " 0 on Σ,
and in addition the zero mean condition:

(3.12)

ż Ω qptqdx " 0 holds for a.e. t P p0, T q.

Based on the adjoint equations, to find the optimal conditions, we need to totally differentiate this functional J pv, ¯ q with respect to the control variable ¯ . In addition, the state equation Spv, u e , w, z, ¯ q " 0 yields the sensitivity equation:

(3.13) BS B ¯ `BS Bu du d ¯ " 0,
where u is the state variable. Eventually, using the Gâteaux derivatives of J as follows:

(3.14) DJ D ¯ " BJ B ¯ `BJ Bu du d ¯ .
Using (3.13) and the adjoint equation

BL Bu " BJ Bu ´λ˚B S Bu " 0, we get (3.15) DJ D ¯ " BJ B ¯ `λ˚B S Bu du d ¯ " BJ B ¯ ´λ˚B S B ¯ " BL B ¯ .
Then, we have

(3.16) ă DJ D ¯ , δ ¯ ą"ă BL B ¯ , δ ¯ ą" 2 ż Ω ¯ .δ ¯ dx ´żQ A m p B B ¯ I ion .δ ¯ dxdt `żQ p BG B ¯ q T s.δ ¯ dxdt,
where p BG B ¯ q T denotes the transpose of the Jacobian matrix of G P R m in point ¯ P R N . We deduce the optimality condition of the optimal control problem (3.1):

(3.17) DJ D ¯ " 0 R N ñ 2 ¯ ´Am ż T 0 p B B ¯ I ion dt `ż T 0 p BG B ¯ q T sdt " 0 R N , @x P Ω.
3.2. Existence of the control. In the following theorem, we show that our optimal control problem (3.1) has a unique solution.

Theorem 3.1. Under the previous Assumptions 2.1, there exists at least one solution ¯ ˚of the optimal control problem (3.1).

Proof. The goal is to prove that there exists ¯ ˚such that J p ¯ ˚q " inf

¯ PC ad J p ¯ q.
Since the cost functional J : C ad Ñ R is bounded from below, with the admissible domain C ad (3.2) is a closed and convex non-empty subset of the Banach space L 8 pΩq N , and J is coercive verifying the property:

(3.18) @p ¯ pkq q kPN P C N ad , lim kÑ8 ˇˇ¯ pkq ˇˇ" `8 ñ lim kÑ8 J p ¯ pkq q " `8,
Then, there exists a bounded minimining sequence

¯ pkq P C N ad such that (3.19) lim kÑ8 J p ¯ pkq q " inf ¯ PC ad J p ¯ q.
Since L 8 pΩq N is a non-reflexive space, it follows from the sequential Banach-Alaoglu theorem that there exists a subsequence of ¯ pkq , which we denote by the same symbol, such that

(3.20) ¯ pkq á ¯ ˚weakly-* in L 8 pΩq N .
Let pv pkq , u pkq e , w pkq , z pkq q be the associated solution of the bidomain model (2.1) with ¯ pkq , by the weak compactness property, see for instance Theorem 5.1 p.58 in [START_REF] Lions | Quelques méthodes de résolution des problemes aux limites non linéaires[END_REF], there exists a subsequence of pv pkq , u pkq e , w pkq , z pkq , ¯ pkq q, denoted by the same indices, such that (3.21)

v pkq á v in L 2 p0, T ; H 1 pΩqq, and d dt v pkq á d dt v in L 2 p0, T ; pH 1 pΩqq ˚q, u pkq e á u e in L 2 p0, T ; H 1 pΩqq, w pkq á w in L 2 p0, T ; L 2 pΩqq k , and B Bt w pkq á B Bt w in L 2 pQq k , z pkq á z in L 2 p0, T ; L 2 pΩqq m , and B Bt z pkq á B Bt z in L 2 pQq m .
By the Rellich-Kondrachov compactness theorem (see [START_REF] Evans | Partial differential equations[END_REF], p. 286), since v pkq , u pkq e are bounded in L 2 p0, T ; H 1 pΩqq, we can assume that v pkq Ñ v and u pkq e Ñ u e strong in L 2 pQq. Since v pkq , w pkq j , for j " 1, . . . , k, and z pkq i for i " 1, . . . , m, are weakly convergent and ψφ 1 ptq P L 2 pQq X L 2 p0, T ; H 1 pΩqq, for any fixed φ P Dp0, T q and ψ P C 8 0 pR d q, we naturally have

(3.22) ´ż T 0 ż Ω d dt v pkq ptqψφdxt " ż T 0 ż Ω v pkq ptqψφ 1 ptqdxt ÝÑ kÑ`8 ż T 0 ż Ω vptqψφ 1 ptqdxt, (3.23) ´ż T 0 ż Ω d dt w pkq j ptqψφdxt " ż T 0 ż Ω w pkq j ptqψφ 1 ptqdxt ÝÑ kÑ`8 ż T 0 ż Ω w j ptqψφ 1 ptqdxt, (3.24) ´ż T 0 ż Ω d dt z pkq i ptqψφdxt " ż T 0 ż Ω z pkq i ptqψφ 1 ptqdxt ÝÑ kÑ`8 ż T 0 ż Ω z i ptqψφ 1 ptqdxt.
Let the corresponding bilinear and continuous forms on H 1 pΩq ˆH1 pΩq associated with σ i,e be a i,e as follows:

(3.25) a i,e pu, ϕq " ż Ω σ i,e ∇u.∇ϕdx, and ψφ P L 2 pQq X L 2 p0, T ; H 1 pΩqq, for any fixed φ P Dp0, T q and ψ P C 8 0 pR d q, we have F j pv pkq , w pkq j q " α j pv pkq q `1 ´wpkq j ˘´β j pv pkq qw pkq j , " α j pv pkq q ´`α j pv pkq q ´αj pvq ˘wpkq j ´`β j pv pkq q ´βj pvq ˘wpkq j ´`α j pvq `βj pvq

(3.26) ż T 0 a i,
˘wpkq j . For i " 1, . . . , m, (3.28)
G i p ¯ pkq , v pkq , w pkq , z pkq q " `Gi p ¯ pkq , v pkq , w pkq , z pkq q ´Gi p ¯ ˚, v pkq , w pkq , z pkq q Gi p ¯ ˚, v pkq , w pkq , z pkq q ´Gi p ¯ ˚, v, w, zq Gi p ¯ ˚, v, w, zq.

and

(3.29)

I ion p ¯ pkq , v pkq , w pkq , z pkq q " N ÿ i"1 ¯ pkq i y i pv pkq q k ź j"1 w p j,i ,pkq j pv pkq ´C `γ i logpz pkq i qq.
Starting by (3.27), we have v pkq Ñ v a.e in Q and α j is continuous (of class C 8 ), so that α j pv pkq q Ñ α j pvq a.e in Q, and α j pv pkq q is bounded in L 2 pQq. It follows from a classical result, see Lemma 1.3 in [START_REF] Lions | Quelques méthodes de résolution des problemes aux limites non linéaires[END_REF], that α j pv pkq q á α j pvq weak in L 2 pQq:

(3.30) @φ P Dp0, T q, ż T 0 ă α j pv pkq ptqq, φptqψ ą dt ÝÑ kÑ`8 ż T 0 ă α j pvptqq, φptqψ ą dt.
As `αj pvq `βj pvq ˘φptqψ P L 2 pQq, the weak convergence of w pkq j in L 2 pQq implies that: (3.31)

@φ P Dp0, T q, ż T 0 ă `αj pvq`β j pvq ˘wpkq j ptq, φptqψ ą dt ÝÑ kÑ`8 ż T 0 ă `αj pvq`β j pvq ˘wj ptq, φptqψ ą dt.
The remaining term in α j is such that @φ P Dp0, T q, } `αj pv pkq q ´αj pvq ˘φptqψ} L 2 pQq ÝÑ kÑ`8 0, } `βj pv pkq q ´βj pvq ˘φptqψ} L 2 pQq ÝÑ kÑ`8 0.

Since }w pkq j } L 2 pQq is bounded, we finally have (3.35) @φ P Dp0, T q, ż T 0 ă `αj pv pkq q ´αj pvq ˘wpkq j ptq, φptqψ ą dt ÝÑ kÑ`8 0,

(3.36) @φ P Dp0, T q, ż T 0 ă `βj pv pkq q ´βj pvq ˘wpkq j ptq, φptqψ ą dt ÝÑ kÑ`8 0.

Let the non-linear term (3.28) associated with t ¯ pkq , ¯ ˚u be t Ḡk , Ḡ˚u , for i " 1, . . . , m:

(3.37)

Ḡk "ă G i p ¯ pkq , v, w, zq, ϕ ą, and Ḡ˚" ă G i p ¯ ˚, v, w, zq, ϕ ą .

We have v pkq Ñ v, w pkq Ñ w and z pkq Ñ z a.e in Q, and G i is a Lipschitz-continuous function, then for any fixed φ P Dp0, T q and ψ P C 8 0 pR d q, we have

(3.38) ż T 0 ă `Gi p ¯ ˚, v pkq ptq, w pkq ptq, z pkq ptqq ´Gi p ¯ ˚, vptq, wptq, zptqq ˘, φptqψ ą dt ÝÑ kÑ`8 0.
In the same way, since ¯ pkq á ¯ ˚weak in L 8 pΩq N and G i is a Lipschitz-continuous function, for any fixed φ P Dp0, T q and ψ P C 8 0 pR d q, we have (3.39) ż T 0 ă `Gi p ¯ pkq , v pkq ptq, w pkq ptq, z pkq ptqq ´Gi p ¯ ˚, v pkq ptq, w pkq ptq, z pkq ptqq ˘, φptqψ ą dt ÝÑ kÑ`8 0.

Similary, for the non-linear ionic term (3.29), Let

I ion p ¯ pkq , v pkq , w pkq , z pkq q " I ion p ¯ pkq , v pkq , w pkq , z pkq q ´Iion p ¯ ˚, v pkq , w pkq , z pkq q `Iion p ¯ ˚, v pkq , w pkq , z pkq q ´Iion p ¯ ˚, v, w, zq 0, @φ P Dp0, T q, ψ P C 8 0 pR d q.

`Iion p ¯ ˚, v, w, zq. Since ¯ pkq á ¯ ˚weak in L 8 pΩq N , v pkq Ñ v,
Also, since I ion is linear with respect to ¯ :" t¯ i u 1ďiďN , for any fixed φ P Dp0, T q and ψ P C 8 0 pR d q, we have (3.42) ż T 0 ă `Iion p ¯ pkq , v pkq , w pkq , z pkq q ´Iion p ¯ ˚, v pkq , w pkq , z pkq q ˘, φptqψ ą dt

" ż T 0 ă `N ÿ i"1 `¯ pkq i ´¯ i ˘yi pv pkq q k ź j"1 w p j,i ,pkq j pv pkq ´C `γ i logpz pkq i qq, φptqψ ą dt ÝÑ kÑ`8 0
Since Ω has Lipschitz boundary, it then satisfies the segment condition (see [START_REF] Adams | Sobolev spaces[END_REF], p.68, definition. 3.21), We say C 8 0 pR d q is dense in H 1 pΩq, (see [START_REF] Adams | Sobolev spaces[END_REF] theorem. 3.22). Therefore from (3.42), we have (3.43) ż T 0 ă `Iion p ¯ pkq , v pkq , w pkq , z pkq q´I ion p ¯ ˚, v pkq , w pkq , z pkq q ˘, φptqϕ ą dt ÝÑ kÑ`8 0, @φ P Dp0, T q, ϕ P H 1 pΩq.

Combined (3.43) with limit (3.41), it implies (3.44) ż T 0 ă `Iion p ¯ pkq , v pkq , w pkq , z pkq q´I ion p ¯ ˚, v, w, zq ˘˘, φptqϕ ą dt ÝÑ kÑ`8 0, @φ P Dp0, T q, ϕ P H 1 pΩq.

Similarly, from (3.39)-(3.38), we can show (3.45) ż T 0 ă `Gi p ¯ pkq , v pkq , w pkq , z pkq q ´Gi p ¯ ˚, v, w, zq ˘, φptqϕ ą dt ÝÑ kÑ`8 0, @φ P Dp0, T q, ϕ P H 1 pΩq.

By passing to the limit in the equations satisfied by pv pkq , u pkq e , w pkq , z pkq , ¯ pkq q, we obtain that pv, u e , w, z, ¯ ˚q is a solution to the bidomain system (2.1).

We conclude that ¯ ˚realizes the minimum of J . In fact, from the weak-* convergence (3.20), it follows that ¯ pkq is also weakly convergent in L 2 pΩq N .

By the weak lower semi-continuity of the regularisation term

ż Ω | ¯ ˚|2 dx in L 2 pΩq N , we deduce that (3.46) J pv, ¯ ˚q ď lim kÑ`8 inf J pv pkq , ¯ pkq q, " inf ¯ PC ad J pv, ¯ q,
and the existence of a minimizer is proved. Now, we demonstrate the numerical procedure to solve the optimality system (3.1).

NUMERICAL APPROXIMATION

In this section, we give a brief overview of the space and time discretization techniques to solve the primal (2.1) and adjoint (3.10) equations numerically. We use a finite element method (FEM) for the spatial discretization and a semi-implicit Euler scheme for the temporal discretization. We solve the optimal control problem (3.1) using the gradient descent method.

' ż Ω A m `Cm B t v `Iion p ¯ , v, w, zq ´Iapp ˘ψdx
`żΩ σ i p∇v `∇u e q.∇ψdx " 0, @ψ P H 1 pΩq, for a.a. t P p0, T q, ' ż Ω `σi ∇v `pσ i `σe q∇u e ˘.∇ψdx " 0, @ψ P H 1 pΩq, with ż Ω ψdx " 0, for a.a. t P p0, T q, ' ż Ω pB t w ´F pv, wqq.ψdx " 0, @ψ P L 2 pΩq k , for a.a. t P p0, T q, ' ż Ω pB t z ´Gp ¯ , v, w, zqq.ψdx " 0, @ψ P L 2 pΩq m , for a.a. t P p0, T q, with vpx, 0q " v 0 pxq, wpx, 0q " w 0 pxq, and zpx, 0q " z 0 pxq, for almost all x P Ω.

Let V h Ă H 1 pΩq be the finite dimensional subspace of piecewise linear basis functions with respect to the spatial grid for the approximation of electrical potentials, gating variables and concentration variables. The approximate solutions of the vectors V P R M and U P R M and matrices W P R kˆM and Z P R mˆM are expressed in the form: (

V ptq "

M ÿ i"1 v i ptqω i , U ptq " M ÿ i"1
u i e ptqω i , W ptq " W pjq ptq

( k j"1 " M ÿ i"1 w pjq i ptqω i ( k j"1 , Zptq " Z pj 1 q ptq ( m j 1 "1 " M ÿ i"1 z pj 1 q i ptqω i ( m j 1 "1 ,
where tω i u M i"1 denote the basis functions, and M is the number of nodal points at the tissue domain. We recall that the variables k and m correspond to the number of gating variables w and intracellular concentration variables z P R m , respectively, verifying the ODEs of the bidomain system (2.1).

The semi-discretization of the primal equations in space results in the differential algebraic system is given as follows:

A m C m M B Bt V " ´Ai V ´Ai U `Am M `Iapp ´Iion p ¯ , V , W pjq , Z pj 1 q q ˘, (4.3) A ie U " ´Ai V , (4.4) M B Bt W pjq " F pjq pV , W pjq q, forj " 1, ..., k, (4.5) M B Bt Z pj 1 q " G pj 1 q p ¯ , V , W , Zq, forj 1 " 1, ..., m, (4.6)
along with initial conditions for V , W pjq and Z pj 1 q , where A ie " tă pσ i `σe q∇ω i , ∇ω j 2 ąu M i,j 2 "1 and A i " tă σ i ∇ω i , ∇ω j 2 ąu M i,j 2 "1 are the stiffness matrices and M " tă ω i , ω j 2 ąu M i,j 2 "1 is the mass matrix. The vector I app is defined by comuting at each component I app i ptq " I app px i , tq. The i th component of the vector

I ion `¯ , V , W , Z ˘is equal to I ion `¯ , v i ptq, w pjq i ptq ( k j"1 , z pj 1 q i ptq ( m j 1 "1
The resulting system is then discretized in time using a semi implicit scheme (Backward Euler formula). Let N P N be a given integer and consider a uniform partition rt n , t n`1 s for 0 ď n ď N ´1, with t n " n∆t, of the time interval of interest r0, T s, with a time step ∆t " T {N .

Denote by pw n , z n , v n , u n e q the approximated solution obtained at time t n . Then pw n`1 , z n`1 , v n`1 , u n`1 e q is computed as follows: For 0 ď n ď N ´1 1. Solve the ionic model, computing w n`1 , z n`1 in each node of Ω.

Using the EDO formulation (2.4), we get an explicit numerical solution for each time step, being able to recover variables w j , for j " 1, . . . , k, as The explicit formula is employed for all the gating variables w of the ionic model (which are the majority of the unknowns of the ionic model). Therefore for the intracellular concentrations z, we apply the backward Euler method: for i " 1, . . . , m,

(4.8) z n`1 i ´zn i ∆t " G i p ¯ , v n , w n`1 , z n`1 q, in Ω.
2. Evaluate in each node the ionic current I ion p ¯ , v n , w n`1 , z n`1 q.

3. With the computed ionic current, solve the first equation of the bidomain model (2.1) to compute the transmembrane potential v n`1 P V h using a first order semi-implicit scheme: (4.9) .∇φ e dx `żΩ σ i ∇v n`1 .∇φ e dx " 0.

ż Ω A m C m v n`1 ´vn ∆t φdx`ż Ω σ i ∇v n`1 .∇φdx " ż Ω A m `In`1 app ´Iion p ¯ , v n , w n`1 , z n`1 q ˘φdx´ż Ω σ i ∇u n e .∇φdx.
Therefore, the equations (4.3) and (4.4) become the following matrix representation: (4.11)

p A m C m ∆t M `Ai qV n`1 " A m C m ∆t M V n `Am M `Iapp n`1 ´Iion p ¯ , V n , W n`1 , Z n`1 q ˘´A i U n , A ie U n`1 " ´Ai V n`1 .

4.2.

Space and time Discretization of the adjoint problem. We use an analogous numerical algorithm as for the primal problem (2.1) to discretize the dual equations (3.10) in the tissue domain. The approximate solutions P , Q, R and S can be expressed in the form: (4.12)

P ptq " M ÿ i"1 p i ptqω i , Qptq " M ÿ i"1 q i ptqω i , Rptq " R pjq ptq ( k j"1 " M ÿ i"1 r pjq i ptqω i ( k j"1 , Sptq " S pj 1 q ptq ( m j 1 "1 " M ÿ i"1 s pj 1 q i ptqω i ( m j 1 "1 .
The following semi-discrete form of the dual problem is obtained:

(4.13) ´Am C m M B Bt P `Am M pB V I ion q T P " ´Ai P ´Ai Q `M pB V F q T R `M pB V Gq T S ` 1 M pV ´V meas q, A ie P " ´Ai Q, ´B Bt R `Am pB W I ion q T P " pB W F q T R `pB Z Gq T S,
´B Bt S `Am pB Z I ion q T P " pB Z Gq T S, with terminal conditions P pT q " R pjq pT q " S pj 1 q pT q " 0, @j " 1, . . . , k, @j 1 " 1, . . . , m. The derivative expression terms are defined as follows (4.14)

B V I ion p ¯ , V , W , Zq :" B Bv I ion `¯ , v i ptq, w pjq i ptq ( k j"1 , z pj 1 q i ptq ( m j 1 "1 ˘(M i"1 P R M , B W I ion p ¯ , V , W , Zq :" B Bw I ion `¯ , v i ptq, w pjq i ptq ( k j"1 , z pj 1 q i ptq ( m j 1 "1 ˘(M i"1 P R kˆM , B Z I ion p ¯ , V , W , Zq :" B Bz I ion `¯ , v i ptq, w pjq i ptq ( k j"1 , z pj 1 q i ptq ( m j 1 "1 ˘(M i"1 P R mˆM , B V F pV , W q :" B Bv F `vi ptq, w pjq i ptq ( k j"1 ( M i"1 P R kˆM , B W F pV , W q :" B Bw F `vi ptq, w pjq i ptq ( k j"1 ( M i"1 P R kˆkˆM , B V Gp ¯ , V , W , Zq :" B Bv G `¯ , v i ptq, w pjq i ptq ( k j"1 , z pj 1 q i ( m j 1 "1 ˘(M i"1 P R mˆM . B W Gp ¯ , V , W , Zq :" B Bw G `¯ , v i ptq, w pjq i ptq ( k j"1 , z pj 1 q i ( m j 1 "1 ˘(M i"1 P R kˆmˆM . B Z Gp ¯ , V , W , Zq :" B Bz G `¯ , v i ptq, w pjq i ptq ( k j"1 , z pj 1 q i ( m j 1 "1 ˘(M i"1 P R mˆmˆM .
For i " 1, ...M , let's denote by R n i " tR pjq,n i u k j"1 P R k , the value of the adjoint state with respect to w and Z n i " tZ pj 1 q,n i u m j 1 "1 P R m , the value of the adjoint state with respect to z at the node x i and at time t n , the computation of these two variables is performed point-wise as follows.

Step 1: For i=1,. . . , M, solve (4.15)

$ ' & ' % `1 ∆t I ´`pB W F q n i ˘T ˘Rn i " 1 ∆t R n`1 i ´Am `pB W I ion q n`1 i ˘T P n`1 i ``pB W Gq n`1 i ˘T S n`1 i , `1 ∆t I ´ppB Z Gq n i q T ˘Sn i " 1 ∆t S n`1 i ´Am ppB Z I ion q n`1 i q T P n`1 i .
Step 2:

(4.16) $ ' ' ' ' ' ' & ' ' ' ' ' ' % `Am C m ∆t M `Ai ˘P n " A m M "`C m ∆t 1 ´pB V I ion q n`1 ˘˚P n`1 ‰ ` 1 M pV n ´V n meas q `M m ÿ j 1 "1 pB V G pj 1 q,n q ˚Spj 1 q,n ´Ai Q n`1 `M k ÿ j"1 pB V F pjq q n ˚Rpjq,n , A ie Q n " ´Ai P n .
Here, 1 P R M is a vector filled with the value 1 and the multiplication "*" is performed point-wise.

Remark 4.1. The retrograde problem (3.10) is fully linear, we use a first order semi-implicit scheme to solve it. The reason is that we separate the ODE system variables R and S from the PDE variables P and Q.

We also solve the bidomain problem sequentially, we first compute P and then we compute Q. This follows the same scheme developed for the primal problem in [START_REF] Fernández | Decoupled time-marching schemes in computational cardiac electrophysiology and ecg numerical simulation[END_REF].

Optimization algorithm.

Given an initial guess of maximal conductance parameters ¯ guess , we solve the optimization problem using the following algorithm based on a gradient descent method. Here, F unc and Grad are positive constants defining the desired tolerance on the cost function and its gradient respectively. The coefficient α is positive and could be fixed or updated at each iteration and M axIterN umber stands for the maximal number of iterations in the optimization procedure.

NUMERICAL RESULTS

In this section, numerical results on the basis of two different test cases are presented. In all tests, the computational domain Ω " r0, 1s ˆr0, 1s Ă R 2 of size 1 cm 2 is fixed and a triangular discretization is used which consists of 5718 elements and 2960 nodes. During the simulations, we fix the time step length ∆t " 0.1ms. The termination of the optimization algorithm is based on the following condition: F unc " 10 ´8 and Grad " 10 ´8. For all the following tests, the desired transmembrane potential v meas is simulated with the physiological Luo Rudy phase I model with its original control parameters using a prescribed time course of a stimulation current which its magnitude is I app ptq " 80µA{cm 2 and its duration is 1ms.

Single maximal conductance parameter estimation:

In the first case, we consider to optimize separately three of the maximal conductance parameters ¯ N a , ¯ si , ¯ K1 representing three different ion channels: sodium, calcium and potassium, respectively. Moreover, if the stopp criteria conditions are not satisfied, the algorithm terminates within a prescribed number of iterations. Here the maximum number of iteration parameter M axIterN umber " 20.

Maximal conductance parameter of the fast inward sodium current ¯ N a :

In this test, we present a numerical results of the estimation of the parameter ¯ N a . Since this parameter is mainly important in the depolarization phase, we consider the cost functional in the time window r0ms, 20mss of the simulation. The exact value N a is equal to 23 mS{cm 2 . We generate the measurement v meas by solving the forward problem using the exact value of N a and we start our optimization procedure using a guess value N a,guess " 1 2 N a " 11.5 mS{cm 2 . Since the cost function depends on the parameters As shown in Fig 1 for both cases, the optimization algorithm converges to the desired control value. But the accuracy is better with 2 " 0.001 than 2 " 0.05 as shown in Table1. From now, we fix 1 " 1 and 2 " 0.001. In order to test the robustness of the algorithm, we add different levels of gaussian noise to the measured data v meas , and we solve the optimization problem following Algorithm 1 for each value of noise. As shown in Fig 1, the algorithm converges for all levels of noise. Table 1 shows that the accuracy is altered with the noise. But for 15% of noise, the relative error on the estimated value of ¯ N a is under 2%. Figure 2 shows the evolution of the cost function Ip¯ N a q and the norm of its gradient with respect to the optimization iterations for different regularization parameter values 2 and noise levels on the measured potential. The table 2 shows the relative error of obtained solution with respect to the 0%, 5% and 10% noise levels. We can see that it converge from the fourth iteration and accuracy the obtained optimal solution of ¯ si seems to be less sensitive to noise compared to optimal solution of N a . The minimum value of the cost functional J(ρ si ) The results in the table 3 show that the optimal solution of ¯ K1 is more sensitive to the noise than ¯ si and less sensitive to noise than ¯ N a . The minimum value of the cost functional J(ρ K1 ) The gradient value of the cost functional J(ρ K1 ) ǫ 2 =0.001, and 0% noise of v meas ǫ 2 =0.001, and 5% noise of v meas ǫ 2 =0.001, and 10% noise of v meas FIGURE 6. Left: Log scale plot of the cost function Ip¯ K1 q. Right: Log scale plot of the norm of its gradient during the optimization procedure.

ǫ 2 =0
ǫ 2 =0

Simultaneous estimate of the maximal conductances parameters:

In this section, we aim to optimize all the maximal conductance parameters of LR1 model simultaneously.

The parameters are: ¯ N a , ¯ si , ¯ K ,¯ K1 ,¯ Kp and ¯ b . The presented numerical results are based on fixing the regularisation parameters 1 " 1, 2 " 0, and adding different levels of gaussian noise to the measured data. First, we have to say that all our attempts to estimate multiple parameters from a single observation v meas failed. Based on the theoretical work [START_REF] Abidi | Ionic parameters identification of an inverse problem of strongly coupled pdes system in cardiac electrophysiology using carleman estimates[END_REF], the stability and the uniqueness of the parameter estimation inverse problem could be guaranteed using as much measurements as the number of the parameters to be optimized. In this work, we choose to modify the position of the stimuli in space in order to obtain the desired number of observations. Suppose that v l meas is the observation associated to the l th applied current I l app , we look for the set of parameters ¯ :" t¯ i u 1ďiďN that solves the following minimization problem:

(5.1) pPq

$ ' ' & ' ' % min ¯ PC ad Ip ¯ q " 1 2 N S ÿ l"1 ż Q ˇˇv l p ¯ q ´vl meas ˇˇ2 dxt,
subject to: v l is solution of the coupled PDE system (2.1) where I app " I l app .

Here N S is the number of stimulus currents with N S ě N . Without loss of generality, using the same strategy presented in paragraph 3.1, the computation of the cost function gradient for multiple observations is given by

(5.2) DJ D ¯ " N S ÿ l"1 ´żQ p BG B ¯ q T s l dxdt ´żQ A m p l B B ¯ I ion dxdt ¯,
where pp l , q l , r l , s l q is the adjoint state of pv l , u l e , w l , z l q solution of the system (2.1) where I app " I l app . Figure 7 shows the evolution of the cost function with respect to the optimization iteration. We remark that in the noise free case (red line) the cost function achieves 10 ´8, whereas for the noisy observation cases the cost function almost stagnates from iteration 30 achieving 2.76e-3 (respectively 1.09e-2) for 5% (respectively 10%) of noise level. Table 4 provides the error percentage of each of the estimated parameters with respect of the noise level. We remark that the accuracy of the estimated parameters is higher than 98% even for 5 and 10 % of noise level. We, particularly, notice that the error on the ¯ N a parameter is less than 0.056% for the three cases which might reveal the sensitivity of the model with respect to the ¯ N a .

Noise levels (%)

¯ 4. Error percentage of the optimal control solution for the different noise levels on the measured potentials pv l meas q 1ďlď6 .

N a ¯ si ¯ K ¯ K1 ¯ Kb ¯ b 0 % 0.
Figure 8 shows snapshots of the action potential distribution on the computational domain for the exact parameters (left), the initial guess parameters (middle) and the estimated parameters (right). First, we remark that for all the parameters the optimization algorithm converges to the desired value. We also remark that during the 15 first iteration the evolution of the parameters in the case of noise-free and noisy data was almost the same. Later in the iterations the noise has an effect on the parameters values. We also remarked that with 5 and 10% of noise on the measured data, the optimization algorithm stops at the 35 th iteration because it couldn't improve anymore the cost function, while for the noise free case the cost function achieves 10 ´8. In this paper, we presented an original approach for the estimation of maximal conductance parameters in the multi-scale cardiac electrophysiology modelling. We used the bidomain system coupled to the physiological Luo Rudy phase I ionic model. We formulated the parameter estimation problem as an optimization procedure in an optimal control problem where the cost function represents the misfit between the measured action potential and the solution of the bidomain model. There are two results in this work. The first part is related to the mathematical analysis of the maximal conductance parameters estimation. Our result is a proof of the existence of a minimizer of the cost function. The proof is inspired from the work of Yan and Veneziani [START_REF] Yang | Estimation of cardiac conductivities in ventricular tissue by a variational approach[END_REF] who proved the existence of the conductivity parameters estimation solution. The second part is related to the numerical estimation of the six conductance parameters introduced in the Luo Rudy phase I ionic model. In order to minimize the cost function, we used a gradient descent method where the gradient is computed by solving an adjoint problem. We detailed the numerical schemes used for the computation of the primal and adjoint problems. Then we conducted different numerical simulation to solve the problem. First, our numerical results show the capability of this approach to estimate separately each of the maximal conductance parameters with a single observation of the action potential in the computational domain. Second, our attempts to solve the optimization problem for multiple parameters estimation using a single measurement failed. Third, using a strategy of multiple observations generated by variating the position of the stimuli in the computational domain, the optimization algorithm were able to converge to the desired parameters with a very satisfactory accuracy. This study shows also that the optimization procedure is robust with respect of noise. Results show that by adding different levels of noise (0%, 5% and 10%), the error on the estimated conductances was less than 2% in the three cases and for all parameters. The challenge is to explore the capability of this method to estimate these physiological parameters when dealing with real life measurement. Finally, we have to say that this study is preliminary and that we didn't explore all of the potential of the optimal control approach. The method here presented allows the estimation of space dependent parameters. This would be subject of our future research.

  pv pkq q ´αj pvq ˘wpkq j ptqφptqψdxdt ˇˇˇď } `αj pv pkq q ´αj pvq ˘φptqψ} L 2 pQq }wpkq j } L 2 pQq ,and similarly for β j , (3.33) ˇˇˇż T 0 ż Ω `βj pv pkq q ´βj pvq ˘wpkq j ptqφptqψdxdt ˇˇˇď } `βj pv pkq q ´βj pvq ˘φptqψ} L 2 pQq }w pkq j } L 2 pQq .Since v pkq Ñ v a.e in Q and α j , β j are continuous (of class C 8 ), we have(3.34)

4. 1 .

 1 Space and time Discretization of the bidomain model. Here we give a description of the spatial discretization of the primal problem by a finite element method based on the weak formulation. A weak solution triple pv, w, zq satisfies: (4.1)

  0, . . . , n, w j p0q " w j,0 .

4 . 1 eP

 41 With the computed transmembrane potential v n`1 , solve the second equation of the bidomain model (2.1) to compute the extracellular potential u n`

1 and 2

 2 used to make a balance between the function of interest `żQ |vp N a q ´vmeas | 2 dxt ˘and the regularization term `żΩ | N a | 2 dx ˘, we first run the optimization procedure with 1 " 1 and we vary 2 from 0.05 to 0.001.

FIGURE 2 .

 2 FIGURE 2. Left: Log scale plot of the cost function Ip¯ N a q. Right: Log scale plot of the norm of its gradient during the optimization procedure.

1 FIGURE 7 .

 17 FIGURE 7. Log scale plot of the cost function Ip¯ q during the optimization procedure.

FIGURE 8 .

 8 FIGURE 8. 2D visualization of solution v at time t " 8msec with exact, guess and optimal parameter values on the left, middle and right respectively.

Figures 9 ,

 9 Figures 9, (respectively, 10 and 11) show the evolution of the parameters ¯ N a and ¯ si (respectively, ¯ K , ¯ K1 and ¯ Kp , ¯ b ) during the optimization iterations. Results are shown for different levels of noise (0%, 5% and 10%) on the measured data. In each figure, we show the target value of the parameter ¯ exact (constant value) and the evolution of the optimized parameters during the iteration.

14 FIGURE 9 .

 149 FIGURE 9. The evolution of the optimal control solution ¯ N a (left) and ¯ si (right) during the optimization iterations.

1 FIGURE 10 .

 110 FIGURE 10. The evolution of the optimal control solution ¯ K (left) and ¯ K1 (right) during the optimization iterations.

FIGURE 11 .

 11 FIGURE 11. The evolution of the optimal control solution ¯ Kp (left) and ¯ b (right) during the optimization iterations.

FIGURE 12 .

 12 FIGURE 12. Action potential v with different value of ¯

  e pu pkq ptq, φptqψqdt ÝÑ

		ż T
	kÑ`8	0	a i,e puptq, φptqψqdt.
	Concerning the non-linear terms in (4.1), we use the equations (2.2)-(2.4)-(2.6) and Assumptions 2.1 to
	write: For j " 1, . . . , k,
	(3.27)

  w pkq Ñ w and z pkq Ñ z a.e in Q, and I ion is a bounded continuous function in L 2 pQq, we have (3.41) ż T

0 ă `Iion p ¯ ˚, v pkq , w pkq , z pkq q´I ion p ¯ ˚, v, w, zq ˘, φptqψ ą dt ÝÑ kÑ`8

  Algorithm 1 Optimization of the maximal conductance parameters ¯ ¯ " ¯ guess , Solve state problem, Solve adjoint problem, while Ip ¯ q ą F unc & } DI D ¯ } ą Grad & iter ď M axIterN umber do ¯ " ¯ ´α ˆDI

	D ¯ .
	Solve state problem,
	Solve adjoint problem,
	Compute the cost function and its gradient,
	end while
	¯ optim " ¯ .

TABLE 1 .

 1 Relative errors of the optimal control solution ¯ N a for all cases.

		25				Optimal control values of ρ Na
		23											2	Noise on v meas (%)	} N a,exact ´ N a } } N a,exact }	ˆ10 2
	N a ρ	20						ρ Na exact ǫ 2 =0.05, and 0% noise of v meas	0.05 0.001	0 % 0 %	1.117 % 0.195 %
		15						ǫ 2 =0.001, and 0% noise of v meas ǫ 2 =0.001, and 5% noise of v meas	0.001	5 %	0.196 %
		11.5						ǫ 2 =0.001, and 10% noise of v meas ǫ 2 =0.001, and 15% noise of v meas	0.001	10 %	0.22 %
		10	0	2	4	6	8	10	12	14	16	18	0.001	15 %	1.58 %
							Optimization iterations
				FIGURE 1. The optimal control so-
				lution for the optimization of N a
				for different values of 2 and differ-
				ent levels of noise.

TABLE 2 .

 2 5.1.2. Maximal conductance parameter of the slow inward-calcium related current ¯ si :In this test, we present a numerical results for the optimization of the parameter ¯ si . Since this parameter acts on the plateau phase, we performed the optimization on a time window [0ms, 400ms]. We consider the initial guess value ¯ si,guess " 3 2 ¯ si,exact " 0.135 mS{cm 2 . The Fig3shows the evolution of the parameter ¯ si during the optimization procedure. Relative errors the optimal control solution ¯ si for all cases.

	0.15
	0.1									2	Noise on v meas (%)	} si,exact ´ si } } si,exact }	ˆ10 2
										0.001	0 %	3.26e-7 %
										0.001	5 %	6.8e-4 %
	0.05
										0.001	10	2.37e-3
	0
	0	1	2	3	4	5	6	7	8	9
		FIGURE 3. The optimal control so-
		lution for the optimization of si .

  .001, and 0% noise of v meas ǫ 2 =0.001, and 5% noise of v meas ǫ 2 =0.001, and 10% noise of v meas Log scale plot of the cost function Ip¯ si q. Right: Log scale plot of the norm of its gradient during the optimization procedure.5.1.3. Maximal time-independent potassium maximal conductance parameter ¯ K1 :In this test, we present a numerical results for the optimization of the parameter ¯ K1 . Since this parameter acts on the repolarization phase, we performed the optimization on a time window [0ms, 400ms]. The initial guess considered is ¯ K1,guess " 3 2 ¯ K1,exact " 0.90705 mS{cm 2 . Fig5shows the evolution of the parameter ¯ K1 during the optimization procedure. The table 3 shows the relative error of the obtained solution with respect to the noise level.

	1
	0.95
	0.8 0.85 0.9			2		Noise on v meas (%)	} K1,exact ´ K1 } } K1,exact }	ˆ10 2
	0.75			0.001	0 %	0.0019 %
	0.65 0.7			0.001	5 %	0.0042 %
	0.6			0.001	10 %	0.083 %
	0.55 0.5	1 FIGURE 5. The evolution of the op-2 3 4 5 6 7 timal control solution ¯ K1 during the optimization iterations. FIGURE 4. Left: 0	8	log 10 ||∇J (ρ si )|| 10 -10 10 -5 10 0 10 5	1	Optimization iterations 4 5 6 The gradient value of the cost functional J(ρ si ) 3 7 ǫ 2 =0.001, and 0% noise of v meas 8 ǫ 2 =0.001, and 5% noise of v meas ǫ 2 =0.001, and 10% noise of v meas TABLE 3. Relative errors of the op-2 timal control solution ¯ K1 for differ-ent noise levels.	9

  .001, and 0% noise of v meas ǫ 2 =0.001, and 5% noise of v meas ǫ 2 =0.001, and 10% noise of v meas

		10	0
	1 )||	10 -2
	10 ||∇J (ρ K	10	-4
	log	10	-6
		10	-8	1	2	3	4	5	6	7	8
								Optimization iterations

  0098 % 0.1967 % 1.0509 % 0.0764 % 0.0578 % 0.1181 %

	5 %	0.0552 % 0.4470 % 1.9545 % 0.5021 % 0.3921 % 0.0043 %
	10 %	0.0498 % 0.23 % 0.9467 % 0.3302 % 0.4888 % 0.0356 %
	TABLE

TUNIS EL MANAR UNIVERSITY, ENIT-LAMSIN, BP 37, LE BELV ÉD ÈRE 1002 TUNIS, TUNISIA. E-mail address: yassine.abidi@enit.rnu.tn E-mail address: moncef.mahjoub@lamsin.rnu.tn

INRIA, BORDEAUX SUD-OUEST, 200 AVENUE DE LA VIELLE TOUR

TALENCE CEDEX FRANCE. E-mail address: nejib.zemzemi@inria.fr

Acknowledgements. This work was supported by the French National Research Agency, grant references ANR-10-IAHU04-LIRYC. This work has also been supported by EPICARD cooperative research program, funded by INRIA international laboratory LIRIMA. The LAMSIN researcher's work is supported on a regular basis by the Tunisian Ministry of Higher Education, Scientific Research and Technology.


Samuel Marrazzo
email: samuel.marrazzo@polytechnique.edu

Thierry Gonçalvès-Novo

François Millet

Jean- Christophe Chanteloup

Low temperature diode pumped active mirror Yb 3+ :YAG disk laser amplifier studies

An experimental study of a static helium gas gap heat switch concept for laser amplification is presented. High single pass gains with large co-sintered ceramic Yb:YAG disks are recorded in the 80-200K temperature range on a diode pumped active mirror amplifier.

Introduction

High average power and energy Diode-Pumped Solid-State Lasers (DPSSL) have been driving strong and continuous interest among the laser physicist community over the past decade. Whereas early stage work on DPSSL was confined on room temperature operated programs [START_REF] Hein | Diode pumped chirped pulse amplification to the joule level and beyond[END_REF][START_REF] Bayramian | The Mercury project: a high average power, gas-cooled laser for inertial fusion energy development[END_REF][START_REF] Yasuhara | 213 W average power of 2.4 GW pulsed thermally controlled Nd:glass zigzag slab laser with a stimulated Brillouin scattering mirror[END_REF][START_REF] Chanteloup | 6.6 J / 2 Hz Yb:YAG diode-pumped laser chain activation[END_REF] with a still un-surpassed record average power set at 550W (55J/10Hz) with the Mercury project [START_REF] Bayramian | The Mercury project: a high average power, gas-cooled laser for inertial fusion energy development[END_REF], a crucial shift took place in the recent years when almost all teams moved towards DPSSL chains relying on at least one low temperature (100 to 200K) operated amplifier [START_REF] Koerner | High efficiency nanosecond pulse amplification based on diode-pumped cryogenic-cooled Yb:YAG[END_REF][START_REF] Divoky | Design of high-energyclass cryogenically cooled Yb3+:YAG multislab laser system with low wavefront distortion[END_REF][START_REF] Gonçalvès-Novo | 14 J/2 Hz Yb 3+ :YAG diode pumped solid state laser chain[END_REF][START_REF] Siebold | Highenergy diode-pumped D 2 O-cooled multislab Yb:YAG and Yb:QX-glass lasers[END_REF][START_REF] Banerjee | DiPOLE: a 10 J, 10 Hz cryogenic gas cooled multi-slab nanosecond Yb:YAG laser[END_REF][START_REF] Divoky | 1-J operation of monolithic composite ceramics with Yb:YAG thin layers: multi-TRAM at 10-Hz repetition rate and prospects for 100-Hz operation[END_REF].

Tradeoffs between thermal management, extraction efficiency, Laser damage threshold and Amplified Spontaneous Emission (ASE) [START_REF] Peters | Amplified spontaneous emission 0.1.threshold condition[END_REF][START_REF] Allen | Amplified spontaneous emission and external signal amplification in an inverted medium[END_REF] management are at the core of DPSSL amplifier engineering and design [START_REF] Chanteloup | Current status on high average power and energy diode pumped solid state lasers[END_REF]. The Laboratoire pour l'Utilisation des Lasers Intenses (LULI) has been engaged in the development of a DPSSL chain named Lucia where amplifying stages rely on Yb 3+ :YAG active mirrors cooled from the High Reflectivity (HR) coated back surface [START_REF] Gonçalvès-Novo | 14 J/2 Hz Yb 3+ :YAG diode pumped solid state laser chain[END_REF][START_REF] Albach | Beyond the current 10J energy level for the Lucia laser system with a cryogenically cooled second amplifier head[END_REF].

Thermal and ASE issues specific to this architecture are detailed in [START_REF] Albach | Deformation of partially pumped active mirrors for high average-power diode-pumped solid-state lasers[END_REF][START_REF] Albach | Influence of ASE on the gain distribution in large size, high gain Yb3+:YAG slabs[END_REF]. The mitigation of ASE deleterious effect (like parasitic lasing) can be performed through several well-known approaches [START_REF] Guch | Parasitic suppression in large aperture disk lasers employing liquid edge claddings[END_REF][START_REF] Ertel | ASE suppression in a high energy Titanium sapphire amplifier[END_REF][START_REF] Gonçalvès-Novo | From 10 to 30 joules with the Lucia laser system: update on current performance and cryogenic amplifier development[END_REF] when operation takes place at room temperature. On the other hand, operating in the 100-200K range makes it almost inevitable to rely on co-sintered Cr 4+ /Yb 3+ :YAG ceramics [START_REF] Yagi | Y 3 Al 5 O 12 ceramic absorbers for the suppression of parasitic oscillation in high-power Nd:YAG lasers[END_REF][START_REF] Albach | Experimental cross evaluation of large size ceramic and crystalline Yb 3+ :YAG laser gain media performance at high average power[END_REF]. Transparent laser ceramics [START_REF] Ikesue | Fabrication and optical properties of high-performance polycrystalline Nd:YAG ceramics for solid-state lasers[END_REF] have been made available in large size only recently [START_REF] Ueda | Scalable ceramic lasers[END_REF]. Although it offers interesting perspectives in terms of laser gain media engineering [START_REF] Akiyama | Laser ceramics with rare-earth-doped anisotropic materials[END_REF], a recent cross evaluation study [START_REF] Albach | Large size crystalline vs. co-sintered ceramic Yb 3+ :YAG disk performance in diode pumped amplifiers[END_REF] of crystalline versus ceramic Yb:YAG amplifier disks pointed out how sensitive to stress (mechanical or thermal) ceramics could be. A careful thermal management minimizing thermal stress is indeed requested to make sure a ceramic based Yb:YAG medium will not totally depolarize the beam it aims at amplifying.

Within the Lucia project, a cooling approach based on a static helium gas gap acting as the heat transport medium between the hot YAG disk and a heatsink was proposed [START_REF] Lucianetti | Active-mirror-laser-amplifier thermal management with tunable helium pressure at cryogenic temperatures[END_REF]. This thermal management concept is known as a "heat switch" or "thermal switch" [START_REF] Duband | A thermal switch for use at liquid helium temperature in space-born cryogenic systems[END_REF][START_REF] Jones | A thermal switch for use at low temperature[END_REF] which allows controlling the heat-conduction paths in the enclosed helium thin layer and consequently the Yb:YAG medium temperature. Based on this architecture, we have built a prototype laser head minimizing the radial temperature gradient generating thermal lens or depolarization through birefringence. This paper describes this active mirror amplifier prototype based on this concept and its qualification in terms of thermal management and gain performances.

Low temperature active mirror amplifier cooled with a static helium gas gap

Thermal management principle

While Lucia DPSSL chain first active mirror amplifier is cooled with high pressure water jets directed towards the reflecting surface [START_REF] Gonçalvès-Novo | 14 J/2 Hz Yb 3+ :YAG diode pumped solid state laser chain[END_REF], a fundamentally different solution has been devised [START_REF] Lucianetti | Active-mirror-laser-amplifier thermal management with tunable helium pressure at cryogenic temperatures[END_REF] for this low-temperature amplifier. The idea is to rely on static helium gas enclosed into a thin (<100 µm) cell while the pressure is adjusted at a level below atmospheric pressure. This is very different from the cooling technique originally developed for Mercury [START_REF] Bayramian | The Mercury project: a high average power, gas-cooled laser for inertial fusion energy development[END_REF] and later implemented in DiPOLE [START_REF] Banerjee | DiPOLE: a 10 J, 10 Hz cryogenic gas cooled multi-slab nanosecond Yb:YAG laser[END_REF] DPSSLs for instance. In these cases, the amplifier heads are indeed cooled with a high velocity helium flow circulating in between several Yb doped slabs; and special care is requested in the design to minimize vibrations and helium pollution. Lucia low temperature amplifier head does not face such considerations and the gain medium is preserved from being altered by impurities and vibrations potentially carried by the helium forced flow.

The principle is based on the dependence of helium thermal conductivity k with the temperature T, the pressure p and the gas gap thickness e as illustrated in the formula:

1 1 2 ( ) 8 1 1 ( ) ( ) 1 1 3 3 bulk bulk k T T k T k T e p R T α α -     ⋅ = ⋅ + ⋅ ⋅ + -       ⋅ ⋅ ⋅ ⋅    
where k bulk (T) is the helium bulk conductivity, R is the specific gas constant for Helium and α 1 and α 2 are the thermal accommodation factors [START_REF] Lucianetti | Active-mirror-laser-amplifier thermal management with tunable helium pressure at cryogenic temperatures[END_REF].

Within the helium layer located between the cold source (copper at the temperature of liquid nitrogen (LN 2 ), 77 K) and the hot source (the Yb:YAG subjected to the pump laser flux), the thermal exchanges occur through two processes: He molecule-wall and moleculemolecule interactions (Fig. 1). The prevalence of one of these two processes (or its simultaneous occurrence) is related to p and T for a given gas gap thickness e. Calculating the molecules mean free path Λ as a function of these variables allows to define three different heat conduction regimes: Viscous (Λ << e), transition (e/10< Λ <10e) and free molecular (Λ >> e) regimes as illustrated in Fig.

2.

Fig. 2. Thermal conductivity of Helium as a function of pressure for temperature of 160 K (plain lines) and 300 K (dashed lines) and various Helium gap thicknesses.

Co-sintered ceramics as gain medium disk

Taking advantage of Yb reabsorption at 1030 nm was proven to be quite efficient to mitigate ASE-triggered effect like parasitic oscillations when operating a Yb:YAG amplifier in the vicinity of room temperature [START_REF] Albach | Current status of the Lucia laser system[END_REF][START_REF] Albach | Amplified spontaneous emission in large size, high gain Yb 3+ :YAG amplifiers: numerical modeling and experimental test bench for foreseen kJ-range diode pumped solid state laser facilities[END_REF]. However, when operation below 150K is requested, 1030 nm reabsorption almost vanishes (cross section < 5.10 -22 cm 2 ) as illustrated in [START_REF] Chanteloup | Defining the optimal gradient doped:YAG disk for room and low temperature diode pumped solid-state laser operations[END_REF]. Absorbing spontaneously emitted radiations at this temperature level therefore requires a YAG disk with a Cr 4+ doped periphery. In order to investigate that approach, a first generation of Cr 4+ /Yb 3+ :YAG cosintered ceramic was manufactured by Konoshima, Japan (Fig. 3, left) where the Cr 4+ doped radial thickness is 5 mm whereas doping concentration is estimated to be 0.25at%. Although, very efficient to prevent parasitic lasing, it proved to carry a too high Cr 4+ doping concentration due to a 6 cm -1 linear absorption coefficient (Table 1). The absorption power was indeed too strong and the resulting heat load located at the Cr 4+ /Yb 3+ interface generated internal stress [START_REF] Albach | Experimental cross evaluation of large size ceramic and crystalline Yb 3+ :YAG laser gain media performance at high average power[END_REF][START_REF] Albach | Large size crystalline vs. co-sintered ceramic Yb 3+ :YAG disk performance in diode pumped amplifiers[END_REF]. A second generation of 3 disks (Fig. 3, right) was then processed with a 4x larger Cr 4+ doped radial thickness of 20 mm coupled to ~4x reduction in absorption coefficient (Table 1). The actual doping level, very likely below 0.1at%, was not made available by the manufacturer. Whereas one side is AR coated, the other (in contact with Helium) is HR coated, both for pump (940 nm) and extraction (1030 nm) wavelengths.

Laser head description

The Lucia water cooled first laser head used to reach the 14 J milestone [START_REF] Gonçalvès-Novo | 14 J/2 Hz Yb 3+ :YAG diode pumped solid state laser chain[END_REF] is pumped at a rather high intensity level (15 to 20 kW/cm 2 ) [START_REF] Chanteloup | Wavelength tunable, 264 J Laser diode array for 10 Hz/1ms Yb:YAG pumping[END_REF]. The high emission cross section of Yb:YAG at low temperature [START_REF] Fan | Cryogenic Yb3+doped solid-state lasers[END_REF][START_REF] Brown | Yb: YAG absorption at ambient and LF cryogenic temperatures[END_REF] allows a significant reduction in pump budget and this second laser head requires therefore only 5.5 kW/cm 2 (~1 ms pulse duration) within a 6 cm 2 homogeneously filled elliptical footprint (43 kW Lastronics pump engines with 76% efficiency). Figure 4 details the head where a first copper block, as heatsink, is brazed on a LN 2 filled Dewar. For operation/maintenance purpose, a second removable copper element is screwed onto the first one; a ~100 µm layer of Indium ensures an efficient thermal contact. The other face of that copper element is carefully rectified and acts as one wall of the low pressure helium gas gap. The second wall is the HR coated side of the YAG disk. The pressure in the helium gas gap is adjusted to the required one with a vacuum pump and a source of helium gas. The thickness of the helium layer is fixed by the thickness of a thermal insulating element made of polyimide (VESPEL ® ) and the torque applied to the disk mount screws applying pressure on a pair of metallic seals. These later seals (HELICOFLEX ® ) are specifically designed to prevent any helium leaks even at low temperature. The removable copper element hosts also a set of heating cartridges giving an extra degree of freedom for temperature control. For instance these cartridges can be used to set the initial temperature at the desired value (i.e. above the ~90K at which the gain medium is usually thermalized before pumping starts).

Experimental qualifications

Thermal

In order to validate the thermal management concept exposed in section 2.1 and [START_REF] Lucianetti | Active-mirror-laser-amplifier thermal management with tunable helium pressure at cryogenic temperatures[END_REF], it was necessary to have access to the gain medium surface and inner temperatures. For this purpose, two 15 mm deep holes were drilled at 3 and 7 mm depth on the periphery of a 77 mm diameter 2 at% doped Yb:YAG crystal of 10 mm thickness. Platinum PT100 probes were inserted in both holes while a thermocouple was strapped onto a thin annular copper disk placed on the top (altitude 0 mm i.e. AR coating surface) of the YAG crystal as illustrated on Fig. 5. Fig. 5. Left, picture of a Yb:YAG crystal in its mount. The three thermal probes allow temperature measurement at 0, 3 and 7 mm depths. While the amplifier is operated at 1, 2 and 5 Hz, a below 5K axial gradient is observed. Helium pressure is set at 700 mbar, gap thickness at 400 µm while pump intensity reaches 5.5 kW/cm 2 homogeneously distributed over the pump area illustrated in red on left picture.

An axial temperature gradient below 5K was observed. However it must be kept in mind that these measurements are performed a couple of centimeters away from the central illuminated part. It therefore does not give information about the gradient experienced by the crystal where the heat load is at its maximum, i.e. at its center. But the main purpose of this experiment was to identify the optimum range of pressure sensitivity as illustrated on Fig. 6.

We confirm a sensitive pressure range between 1 and 10 mbar. Above that level, the helium layer thermal conductivity is not evolving anymore as illustrated in Fig. 2 and therefore allows cooling control around 130K for a 2 Hz pump rate. Similar experiments were performed with a much thinner gas gap of 55 µm and an identical pump intensity of 5.5 kW/cm 2 for which a minimal cooling temperatures of 110K (5Hz, 27.5W/cm 2 ) and 150K (10 Hz, 55W/cm 2 ) were achieved.

It should also be pointed out that this experimental configuration allowed us to record Yb:YAG absorption cross sections between 80 and 300 K as illustrated on Fig. 2 of [START_REF] Chanteloup | Defining the optimal gradient doped:YAG disk for room and low temperature diode pumped solid-state laser operations[END_REF]. Fig. 6. Yb:YAG ceramic temperature versus helium pressure in a static 400 µm thick cell. Isomean free path (MFP) curves are displayed for 40, 4 and 0.4 µm. The first one defines the limits between transition (blue) and viscous (green) regimes as defined in section 2.1. Data points were collected while the 5.5 kW/cm 2 pump intensity was delivered at 2 Hz with a 1 ms pulse duration.

Gain

[25] was dedicated at cross evaluation of 45 mm diameter co-sintered ceramics (Generation 1) and crystals at room temperature. This new prototype allows us to extend that study by comparing 77 mm disks operated at low temperature. A narrow (<0.1nm) 1030.3 nm CW 60 mW laser source was first used as a probe beam. Experimental conditions were identical for both disks.

The Yb:YAG crystal performances are clearly affected by the lack of ASE control scheme (no 1030 nm absorption of Yb at these temperatures and absence of Cr ions playing a similar role). After a few hundreds of microseconds, strong gain saturation takes place. Even at relatively high (200K) temperature, the gain never exceed a value of 2 (left, Fig. 7). These ASE plateaus are not observed on the Cr 4+ /Yb 3+ :YAG co-sintered ceramic gain curves and a 50% increase in single pass gain can be attained (right, Fig. 7). For these graphs and the followings, the reported temperatures are recorded through a peripheral copper ring in contact with the AR coated surface of the YAG disk. Maximum gain values recorded with the ceramics are reported on Fig. 8 for different adjusted pressures in the gas gap (i.e. temperature control). The blue square are the maxima of Fig. 7 right graph curves (2 Hz, 150 µm) whereas the red circles have been obtained similarly while operating the amplifier at 5 Hz with a 55 µm thickness Helium gap). The single pass gains increase with decreasing temperature until a temperature in the 160-190 K range where it somehow saturates and then decreases. In order to explain this behavior, one could invoke ASE losses but, although likely to play a negative role, ASE does not appear to play a significant one when carefully looking at build-up phase at Fig. 7 (right image). The reason has in fact to be linked to the nature of the source used for this experiment, i.e. a 1030.3 nm narrow one. Figure 9 shows that, in order to reach much higher gains, a source tunable down to 1029.5 nm is requested. Fig. 9. Maximum emission cross section (red, left scale) and wavelength allowing to reach that maximum (green, right scale) versus temperature. Reaching high cross section values around 1.10 -19 cm 2 requires a 1029.5 nm source. This graph was derived from data kindly shared by [35]. Gain measurements were then performed while adjusting the wavelengths to better match emission cross section maxima. Figure 10 displays single pass gain obtained with a 1029.7 nm narrow (100 pm) 35 mW CW source. Gain values vary between 12 and 20 with YAG temperatures across the 142 to 84 K range.

Similar experiment where performed with Lucia front end nominal configuration with a 20 mm diameter beam of 10 ns pulses at 2 Hz at 1029.7 nm; maximum energy per pulse being half a Joule. Figure 11 shows that single pass gain stays above 6 in the 100 to 460 mJ explored energy range. The YAG temperature was close to 140K. Any Helium pressure adjustment imposes a temperature stabilization period which depends of the operating pressure range (1, 10 or 100 mbar). A ~1% stabilization can require up to half an hour. As soon as the point of operation is defined (helium gap thickness, pressure, pump load…), it should not be necessary to adjust the helium pressure.

Conclusions

We have experimentally validated at low temperature the static helium gas gap heat switch concept for active mirror laser amplifier. A prototype was built and integrated into the Lucia laser chain. It can host YAG disks up to 77 mm in diameter in an active mirror configuration. A simple helium pressure adjustment between 1 and 10 mbars allows us to explore the amplifier behavior over a 100K wide temperature range around 140K. Colder or warmer operating temperatures could be achieved easily by adjustment of the VESPEL ring thickness allowing the helium layer to be tuned from tens of µm to hundreds of µm.

Thanks to Konoshima optimization work in terms of ultralow Cr 4+ doping concentration, a new generation of co-sintered Cr 4+ /Yb 3+ :YAG disks were used allowing a 2x single pass gain improvement from 3 to around 7. Foreseen evolution of this low temperature amplifier will be devoted to a 2nd generation amplifier able to operate in the same temperature range without the need of LN 2 .

Fig. 1 .

 1 Fig. 1. Schematic illustration of heat transfer processes at work between faces of a helium filled cavity.

Fig. 3 .

 3 Fig. 3. Illustration of both generations of Cr 4+ /Yb 3+ :YAG co-sintered ceramics.

Fig. 4 .

 4 Fig. 4. CAD section view of the low temperature laser head.

Fig. 7 .

 7 Fig. 7. Single pass gain temporal behavior for two 1 cm thick, 77 mm in diameter Yb:YAG disks operated in the low temperature Lucia amplifier. The 5.5 kW/cm 2 pump is limited at 1 ms duration (blue area) with a 2 Hz repetition rate. The left graph pictures gain curves obtained in the crystal case while the right graph is dedicated to the co-sintered ceramic. The Helium gas gap thickness is 150+/-5 µm.

Fig. 8 .

 8 Fig. 8. Maximum gain observed versus Cr 4+ /Yb 3+ :YAG co-sintered ceramic temperature. Extraction wavelength is 1030.3 nm. Pump intensity is 5.5 kW/cm 2 with 1 ms pulse duration. The blue squares are the maximum values observed on Fig. 7 (right), i.e. obtained at 2 Hz with a 150+/-5 µm Helium gas gap thickness. The red disks are gain maxima obtained at 5 Hz with a 55+/-5 µm gap.

Fig. 10 .

 10 Fig. 10. Small signal gain recorded at 1029.7 nm while pumping the Cr 4+ /Yb 3+ :YAG cosintered ceramic at 5.5kW/cm 2 intensity and single shot, 0.5, 1 and 2 Hz repetition rate. Respective temperatures are 84, 98, 117, 142 K. Helium pressure is 1 bar.

Fig. 11 .

 11 Fig. 11. Left, blue: single pass gain recorded in pulse mode (10 ns / 2 Hz) with the Cr 4+ /Yb 3+ :YAG co-sintered ceramic. Right, Red: Amplified energy in a single pass.

Table 1 . Co-sintered ceramics specifications.

 1

	Ceramics	Thickness	Yb 3+	Cr 4+	Cr 4+	Yb 3+	Cr 4+	Absorption	Transmission
			doped	doped	doped	doping	doping	coefficient	over t
			Ø	Ø	radial
					thickness
					t
	Unit	mm	mm	mm	mm	at%	at%	cm -1	%
	Generation 1	7	35	45	5	2	0.25	6	4.98
	Generation 2-A							2.92	0.29
	Generation 2-B	10	37	77	20	1	not known	1.65	3.71
	Generation 2-C							1.58	4.24

Acknowledgments

The research leading to these results has received funding from LASERLAB-EUROPE (grant agreements no. 228334 and n° 284464, EC's Seventh Framework Programme). The authors thank J. Hein and J. Körner for sharing emission cross section data used for Figs. 8 and9.


Abdeldjalil Aissa

El Bey

Karim Abed-Meraim

Yves Grenier

A Aïssa-El-Bey

UNDERDETERMINED BLIND SEPARATION OF AUDIO SOURCES IN TIME-FREQUENCY DOMAIN

come

Underdetermined blind source separation of audio sources in time-frequency domain

INTRODUCTION

Blind source separation (BSS) considers the recovery of unobserved original sources from several mixtures observed at the output of a set of sensors. Each mixture contains a combination of the sources that results from the mixing medium between the sources and the sensors. The term "blind" indicates that no a priori knowledge of the sources and the medium structure is available. To compensate for this lack of information, the sources are usually assumed to be statistically independent. Blind source separation has application in different areas, such as communications, speech processing, image processing and biomedical engineering [START_REF]Blind estimation using higher-order statistics[END_REF]. A challenging problem of BSS occurs when there are more sources than sensors, and this is referred to as underdetermined blind source separation (UBSS). A time-frequency based UBSS algorithm has been recently proposed in [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequancy masking[END_REF][START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF] to successfully separate speech sources using time-frequency distributions (TFDs). This algorithm provides good separation performance when the sources are disjoint in the TF plane. It also provides the separation of TF quasi-disjoint sources, that is the sources are allowed to have a small degree of overlapping in the TF plane. However, the intersection points in the TF plane are not directly treated. More precisely, a point at the intersection of two sources is clustered "randomly" to belong to one of the sources. As a result, the source that picks up this point now contains some information from the other source while the later source loses some information of its own. However, for the other source, there is an interference at this point, hence the separation performance may degrade if no treatment is provided for this. An increasing in the number of intersection points degrades the separation quality. In this paper, we propose another algorithm, combining the TF-UBSS with subspace projection, that allows an explicit treatment of the intersection points. The main assumption used in this algorithm is that the number of sources simultaneously present at an intersection point in the TF plane cannot exceed the number of sensors.

PROBLEM FORMULATION

Data model

Let s(t) = [s 1 (t), . . . , s N (t)]
T represent the N nonstationary source signals. The source signals are transmitted through a medium so that an array of M linear sensors picks up a set of mixed signals represented by an M -dimensional vector

x(t) = [x 1 (t), . . . , x M (t)]
T . We consider the instantaneous mixing medium that is modeled as follows

x(t) = As(t) + w(t), (1)
where A = [a 1 , . . . , a N ] is the mixing matrix and w(t) = [w 1 (t), . . . , w M (t)] T is the observation noise vector. We assume that any M columns of A are linearly independent. The goal of BSS is to recover s(t) from x(t). When M < N , the problem becomes UBSS. Let Ω 1 and Ω 2 be the TF supports (i.e. the points of TF plane where the local energy of the considered sources is non-zero) of two sources s 1 (t) and s 2 (t).

If Ω 1 ∩ Ω 2 = ∅, the sources are said to be non-disjoint in the TF plane. The second assumption is that the sources are not necessarily disjoint, and in particular, there exist, at most, simultaneously (M -1) sources at the same TF point. However, we still assume that there exists for each source signal a region R i in the TF plane where it exists alone, i.e. the energy of the other sources are negligible at the TF points within the considered region. For this reason, a TF representation is commonly referred to as a time-frequency distribution (TFD). TFDs have been applied to a wide variety of engineering problems. Specifically, they have been successfully used for signal recovery at low signal-to-noise ratio (SNR), accurate estimation of the instantaneous frequency (IF), signal detection in communication, radar processing and for the design of time-varying filter. For more details on TFDs and related methods, see for example the recent comprehensive reference [START_REF] Boashash | Time Frequency Signal Analysis and Processing: Method and Applications[END_REF]. The method presented in this paper, uses the Short-Time Fourier Transform (STFT) that is defined as:

Time-frequency representation

S x (t, f ) = m=+∞ m=-∞ h(t -m)x(m)e -j2πf m (2)
where h(t) is the Hamming window.

CLUSTER-BASED TF-UBSS APPROACH FOR DISJOINT SOURCES

In this section, we briefly review the STFT method in [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequancy masking[END_REF], and propose a cluster-based linear TF-UBSS algorithm using STFT to avoid some of the drawbacks in [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequancy masking[END_REF]. First, under the transformation into the TF domain using the STFT, the model in (1) leads to:

S x (t, f ) = AS s (t, f ), (3)
where

S xi (t, f ) = m=+∞ m=-∞ h(t -m)x i (m)e -j2πf m (4a) S x (t, f ) = [S x1 (t, f ), . . . , S xM (t, f )] T . ( 4b
)
and S s (t, f ) is the N × 1 source STFT vector. To avoid processing all TF points (and hence to reduce computational cost), we apply first a noise thresholding as that for each timeslice (t, f ):

If S x (t, f 0 ) max f { S x (t, f ) } > ǫ, then keep (t, f 0 ), ( 5
)
where ǫ is a small threshold (typically, ǫ = 0.05). Then, the set of all selected points, Ω, is expressed by Ω = N i=1 Ω i , where Ω i is the TF support of the source s i (t). Under the assumption that all sources are disjoint in the TF domain, (3) is reduced to

S x (t, f ) = a i S si (t, f ), ∀(t, f ) ∈ Ω i , ∀i = 1, • • • , N. (6)
where the source STFT vector has been reduced to only the STFT of the source s i (t). Now, in [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequancy masking[END_REF], the structure of the mixing matrix is particular as such it has only 2 rows (i.e. the method uses only 2 sensors) and the first row of the mixing matrix contains all 1. Then, ( 3) is expanded to

S x1 (t, f ) S x2 (t, f ) = 1 . . . 1 a 2,1 . . . a 2,N    S s1 (t, f ) . . . S sN (t, f )    , (7)
and ( 6) to

S x1 (t, f ) S x2 (t, f ) = 1 a 2,i S si (t, f ),
which results in

a 2,i = S x2 (t, f ) S x1 (t, f ) . (8)
Therefore, all the points for which the ratios on the right-hand side of ( 8) have the same value form the TF support Ω i of a single source, say s i (t). Then, the STFT estimate of s i (t) is computed by:

Ŝsi (t, f ) = S x1 (t, f ), ∀(t, f ) ∈ Ω i , 0, otherwise.
Finally, the source estimate ŝi (t) is obtained by converting Ŝsi (t, f ) to the time domain using inverse STFT [START_REF] Griffin | Signal estimation from modified short-time Fourier transform[END_REF]. For more details, refer to this paper. It is observed that the structure of the mixing matrix, as expressed in (7) has some limiting factors. First, the extension of the UBSS method in [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequancy masking[END_REF] to more than two sensors is not obvious. Second, the division on the right-hand side of ( 8) is prone to error if the denominator is close to zero. To avoid the above mentioned problems, we propose here a modified version of the method valid for any number of sensors. This method is now referred to as the cluster-based linear TF-UBSS algorithm. The clustering method proceeds as follows: first compute the spatial direction vectors by:

v(t, f ) = S x (t, f ) S x (t, f ) , (t, f ) ∈ Ω, (9)
and force them, without loss of generality, to have the first entry real and positive. Next, we cluster these vectors into N classes {C i | i = 1, • • • , N }, using the k-mean clustering algorithm [START_REF] Frank | The data analysis handbook[END_REF]. The collection of all points, whose vectors belong to the class C i , now forms 2. Vector clustering by (9) and [START_REF] Frank | The data analysis handbook[END_REF].

3. Mixing matrix and source STFT estimation by ( 10) and (11).

4. Source TF synthesis by [START_REF] Griffin | Signal estimation from modified short-time Fourier transform[END_REF].

the TF support Ω i of the source s i (t). Then, the column vector a i of A is estimated as the centroid of this set of vectors:

âi = 1 |C i | (t,f )∈Ωi v(t, f ), ( 10
)
where |C i | is the number of vectors in this class. Therefore, we can estimate the STFT of each source s i (t) (up to scalar constant) by:

Ŝsi (t, f ) = âH i S x (t, f ), ∀ (t, f ) ∈ Ω i , (11)
since, from (6), we have

âH i S x (t, f ) = âH i a i S si (t, f ) ∝ S si (t, f ), ∀ (t, f ) ∈ Ω i .
This algorithm is summarized in Table 1.

SUBSPACE-BASED TF-UBSS APPROACH FOR NON-DISJOINT SOURCES

We propose here to use an appropriate subspace projection to estimate the TFDs of the individual sources, under the previously stated data assumptions. Under the TF non-disjoint condition, consider a source point (t, f ) ∈ Ω at which there are K contributing sources s α1 (t), . . . , s αK (t), with K < M . Then, (3) is reduced to the following

S x (t, f ) = ÃS s(t, f ), ∀(t, f ) ∈ Ω (12)
where à and s are defined by:

s = [s α1 (t), . . . , s αK (t)] T , (13a) Ã = [a α1 , . . . , a αK ]. ( 13b
)
Let Q Ã be the orthogonal projection matrix onto the noise subspace of Ã. Then, Q Ã can be computed by:

Q Ã = I -Ã ÃH Ã -1 ÃH . (14)
We have the following observation:

Q Ãa i = 0, i ∈ {α 1 , . . . , α K } Q Ãa i = 0, otherwise . (15)
Table 2. Subspace-based TF-UBSS algorithm 1. STFT computation.

2. Single-source point selection; mixing matrix estimation by, k-mean algorithm.

3. For all source points, perform subspace-based TFD estimation of sources by ( 14), ( 16) and (17).

4. Source TF synthesis by [START_REF] Griffin | Signal estimation from modified short-time Fourier transform[END_REF].

If A is known or a priori estimated, then this observation gives us the criterion to detect the indices α 1 , . . . , α K ; and hence, the contributing sources at the considered TF point (t, f ). In practice, to take into account noise, one detects the column vectors of à minimizing:

{α 1 , . . . , α K } = arg min β1,...,βK Q ÃS x (t, f ) | Ãβ (16)
where Ãβ = [a β1 , . . . , a βK ].

Next, TFD values of the K sources at TF point (t, f ) are estimated by:

Ŝs (t, f ) ≈ Ã# S x (t, f ). (17)
where Ã# represents the pseudo-inverse of Ã. Now, to apply the above procedure, we need to estimate A first. This is performed here by clustering all the spatial direction vectors in (9) as for the preview TF-UBSS algorithm. Then within each class C i we estimate the far-located vectors from the centroid (in the simulation we estimate vectors v(t, f ) such that: v(t, f )-âi > 0.8 max

v(t,f )∈Ωi v(t, f )-âi
leading to a reduced size class Ci . This is to essentially keep the vectors corresponding to the TF region R i (which are ideally equal to the spatial direction a i of the considered source signal). Finally, the i th column vector of A is estimated as the centroid of Ci . Table 2 provides a summary of the subspace projection based TF-UBSS algorithm.

SIMULATIONS AND RESULTS

Simulation results are illustrated in the figures below. In this simulation, we have used uniform linear array of M = 3 sensors. It receives signals from N = 4 independent speech sources, lasting 8192 samples. In figure 2, the upper line represents the original source signals, the second line represents the M mixtures and the bottom one represents the sources estimates by our algorithm. In figure 3, we compare the performance of our method with the TF-UBSS method of

CONCLUSION

This paper introduces a new approach for blind separation of non-disjoint and nonstationary sources using TFDs. The proposed method can separate more sources than sensors and provides, in the case of non-disjoint sources, a better separation quality than the method proposed in [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequancy masking[END_REF]. This method is based on a vector clustering procedure that estimates the mixing matrix A, and subspace projection to separate the sources at the intersection points in the TF plane.

Fig. 1 .

 1 Fig. 1. (a)-TF disjoint, (b) TF non-disjoint

Fig. 2 .

 2 Fig. 2. Blind source separation example for 4 speech sources and 3 sensors in instantaneous mixture case: the upper line represents the original source signals, the second line represents the M mixtures and the bottom one represents estimates of sources by our algorithm.

Fig. 3 .

 3 Fig. 3. NMSE versus SNR for 4 speech sources and 3 sensors: comparison of the performance of our algorithm with the modified TF-UBSS

Table 1 .

 1 Cluster-based TF-UBSS algorithm 1. Mixture STFT computation by (4).

Table 1

 1 (i.e. modified method of[START_REF] Yilmaz | Blind separation of speech mixtures via time-frequancy masking[END_REF]). The plots represents the average normalized MSE (NMSE) of the estimated sources versus the SNR in dB. For the subspace-based method we have used K = 2 for all TF points. As can be observed, a significatif gain is obtained, thanks to our subspace projection.

						0.8																	1
						0.6																	0.8
																							0.6
						0.4
																							0.4
						0.2																	0.2
						0																	0
						-0.2																	-0.2
																							-0.4
						-0.4
																							-0.6
						-0.6																	-0.8
						-0.8	0	1000	2000	3000	4000	5000	6000	7000	8000								-1	0	1000	2000	3000	4000	5000	6000	7000	8000
									(a) s 1 (t)													(b) s 2 (t)
						0.8																	1
						0.6																	0.8
																							0.6
						0.4
																							0.4
						0.2																	0.2
						0																	0
						-0.2																	-0.2
																							-0.4
						-0.4
																							-0.6
						-0.6																	-0.8
						-0.8	0	1000	2000	3000	4000	5000	6000	7000	8000								-1	0	1000	2000	3000	4000	5000	6000	7000	8000
									(c) s 3 (t)													(d) s 4 (t)
	0.8														0.6													0.6
	0.6
															0.4													0.4
	0.4
															0.2													0.2
	0.2
	0														0													0
	-0.2
															-0.2													-0.2
	-0.4
															-0.4													-0.4
	-0.6
	-0.8	0	1000	2000	3000	4000	5000	6000	7000	8000				0	1000	2000	3000	4000	5000	6000	7000	8000				0	1000	2000	3000	4000	5000	6000	7000	8000
				(e) x 1 (t)								(f) x 2 (t)								(g) x 3 (t)
						0.8
						0.6
						0.4
						0.2
						0
						-0.2
						-0.4
						-0.6
						-0.8	0	1000	2000	3000	4000	5000	6000	7000	8000
									(h) s 1 (t)