Mathematical analysis of cardiac electromechanics with physiological ionic model

This paper is concerned with the mathematical analysis of a coupled elliptic-parabolic system modeling the interaction between the propagation of electric potential coupled with general physiological ionic models and subsequent deformation of the cardiac tissue. A prototype system belonging to this class is provided by the electromechanical bidomain model, which is frequently used to study and simulate electrophysiological waves in cardiac tissue. The coupling between muscle contraction, biochemical reactions and electric activity is introduced with a so-called active strain decomposition framework, where the material gradient of deformation is split into an active (electrophysiology-dependent) part and an elastic (passive) one. We prove existence of weak solutions to the underlying coupled electromechanical bidomain model under the assumption of linearized elastic behavior and a truncation of the updated nonlinear diffu-sivities. The proof of the existence result, which constitutes the main thrust of this paper, is proved by means of a non-degenerate approximation system, the Faedo-Galerkin method, and the compactness method.


Introduction 4
The heart is the muscular organ that contracts to pump blood throughout the body. Failure 5 in its contraction leads to sudden cardiac death which is classified as the main cause of mortality 6 in the world. The contraction of the heart is initiated by an electrical signal called action poten-7 tial starting in the sinoatrial node. The electrical signal then travels through the atria and the used nondegenerate approximation systems to which they applied the Faedo-Galerkin scheme [3], 1 Bourgault et al. introduced a "Bidomain" operator and used a semigroup approach [4], Matano 2 and Mori derived global classical solutions [5] and Veneroni proved the existence and uniqueness 3 of a strong solution with more involved ionic models using a fixed point approach with strong 4 assumptions on the initial data [6]. 5 Still at the macroscopic level, cardiac deformation can be modeled by the equations of motion 6 for a hyperelastic material, written in the reference configuration. However, like any living tissue, 7 there is a difficulty in applying the principles of force balance to cardiac tissue due to its ability to 8 actively deform. In other words, its contraction is influenced by intrinsic mechanisms taking place 9 at the microscopic level. This ability is taken into account in the literature following different 10 approaches. One common option is to assume that stresses are additively decomposed into active 11 and passive parts and it is called the active stress formulation ([7, 8, 9, 10, 11])). In this paper, 12 we follow the active strain formulation, [12,13], where the deformation gradient is factorized into 13 active and passive factors, and fiber contraction rewrites in the mechanical balance of forces as 14 a prescribed active deformation. Furthermore, this decomposition incorporates the micro-level 15 information on fiber contraction and fiber directions in the kinematics [14]. These mechanisms 16 essentially translate into a dependence of the strain energy function on auxiliary internal state 17 variables, which represent the level of mechanical tissue activation passed across scales [15]. For 18 comparisons between the two approaches in terms of numerical implementation, constitutive is- 19 sues, and stability, we refer the reader to [16,17]. 20 Mathematical analysis of general nonlinear elasticity can be found in [18,19], whereas applica-  27 and derived constraints on the initial data. Andreianov et al. also assumed linearized elasticity 28 equations but they adopted the active strain formulation and employed the bidomain model cou-29 pled with FitzHugh-Nagumo ionic model. This is the setting we employ in the present work, but 30 we use a general physiological ionic model which kinetics overlap with Beeler-Reuter model [28] or 31 Luo-Rudy model [29]. The electrical to mechanical coupling is obtained by considering that the 32 active part of deformation incorporates the effect of calcium dynamics. We also consider that the 33 evolution of electrical potential, governed by the bidomain equations, depends on the displacement 34 which enters into the equations upon a change of coordinates from Eulerian to Lagrangian. 35 Putting our contributions into perspective, we first note that up to the author's knowledge, exis-36 tence of solution of an electromechanical model coupled with physiological ionic model has never 37 been rigourously mathematically anaylzed. Moreover our paper admits a rigorous mathematical 38 treatment, yielding the existence of weak solutions of our model. We point out that our model is 39 degenerate, strongly nonlinear and so no maximum principle applies. We want to mention that 40 we have not been able to prove uniqueness of weak solutions because of the presence of nonlinear 41 lower-order terms in our model. Furthermore, comparing to the work [6] (where the author proves of such solutions to fabricate weak solutions of the original systems. Again convergence is achieved 1 by a priori estimates and compactness arguments. On the technical side, we point out that the 2 passage to the limit in the pressure term is not straightforward due to the artificial compressibility 3 assumption along with the use of "Navier-type" boundary conditions. 4 The contents of this paper are organized as follows. Section 2 describes the cardiac electrome-5 chanical model we adopt, presenting the equations of passive nonlinear mechanics, the bidomain 6 system, and the active-strain-based coupling strategy. We also list the basic assumptions of the 7 model and provide a definition of weak solution. In Section 3 we state and prove the solvability 8 of the continuous problem employing Faedo-Galerkin approximations and compactness theory to 9 obtain the existence of solution of a regularized problem in the first place. Then the existence of 10 weak solutions for the original problem is given in Section 4 by using (one more time) a priori 11 estimates and compactness arguments. In Section 5, we close our contribution with some remarks 12 and discussion of future directions.  2.1. A general nonlinear elasticity problem. From the mechanical view point, we consider the heart as a homogeneous continuous material occupying in the initial undeformed configuration a bounded domain Ω R ⊂ R d (d = 3) with Lipschitz continuous boundary ∂Ω R . Its deformation is described by the equations of motion written in the reference configuration Ω R . The current configuration is the deformed configuration denoted by Ω. We look for the deformation field φ : Ω R → R d that maps a material particle occupying initially the position X to its current position x = φ(X). We denote by F := ∇ X φ, the deformation gradient tensor where ∇ X is the gradient operator with respect to the material coordinates X, noting that det(F) > 0. The cardiac tissue is also assumed to be a hyperelastic incompressible material. In other words, there exists a strain stored energy function W = W(X, F), differentiable with respect to F, from which constitutive relations between strain and stresses are obtained. In addition, the first Piola stress tensor P, which represents force per unit undeformed surface is given by: where Cof (·) is the cofactor matrix, and p is the Lagrange multiplier associated to the incompressibility constraint: det(F) = 1 and interpreted as " hydrostatic pressure". The balance equations in the reference configuration for deformations and pressure read as: Find φ, p such that ∇ · P(F, p) = g in Ω R , 2.2. The bidomain equations. The electrophysiological aspect of the heart is incorporated in the model through the widely used bidomain equations [1]. The unknowns are the intracellular (i) and extracellular (e) electric potentials v the gating or recovery variables w = w(t, x) = (w 1 , · · · , w k ), and the concentration variable z = z(t, x) = (z 1 , · · · , z m ) at (t, x) ∈ Ω T := (0, T ) × Ω, where T is the final time instant. Cardiac electrical conductivity is represented in the global coordinate system by the orthotropic tensors n}, are the intra-and extracellular conductivities along, transversal, and normal to the fibers' direction, respectively. The direction of the fibers is a local quantity used to determine the principal directions of propagation, thus we have d s = d s (x), s ∈ {l, t, n}. The externally applied stimulation currents corresponding to the intra-and extracellular spaces are represented by the functions I i s and I e s , respectively. The bidomain equations are given by: Here c m is the capacitance and χ is the membrane surface area per unit volume. For simplicity, we shall suppose that χ = 1 and c m = 1. Problem (2.3) is provided with homogeneous Neumann boundary conditions on the intra-and extracellular potentials. In the physiological membrane model, the ionic current I ion has the following general form Herein, d i is the maximal conductance associated with the i t h current, f i is a gating function depending on the transmembrane potentiel v, n i,j is positive integer and E := r i log ze zi is equilibrium (Nernst) potential (r i is a constant and z e is an extracellular concentration). Moreover, the dynamics of the gating variable w is described in the Hodgkin-Huxley formalism by a system of ODEs governed by the following equation The functions α j and β j are positive with the following form where ρ 1,κ , ρ 3,κ , ρ 4,κ , v ≥ 0 and ρ 2,κ , ρ 5,κ > 0 are constants.

4
The choice of the membrane model to be used is reflected in the functions I ion (v, w, z), R(v, w) and 5 G(v, w, z). For a physiological description of the action potential, we will consider a fairly general The electrical to mechanical coupling is done through the "active strain model" [13] where the deformation gradient F is factorized into a passive component F p and an active component F a , F = F p F a . The tensor F p acts at the tissue level and accounts for both deformation of the material needed to insure compatibility and possible tension due to external loads. The tensor F a represents the distortion that dictates deformation at the fiber level and depends on the electrophysiology through the relation, [14]: where γ s , s ∈ {l, t, n} are quantities that depend on the electrophysiology equations. Such a factorization of the deformation tensor F assumes the existence of an intermediate configuration between the reference and the current frames. In that configuration, the strain energy function depends solely on the deformation at the macroscale F p , [30]: and the Piola stress tensor is given by : (see also Refs. [14,30] ourselves to a phenomenological description of local activation in terms of the gating variables. 10 The scalar fields γ l , γ t and γ n can be written as functions of a parameter γ: 11 γ l,t,n = γ l,t,n (γ), where γ l,t,n : R → [−Γ l,t,n , 0] are Lipschitz continuous monotone functions. The values Γ l,t,n should be small enough, in order to ensure that det(F a ) stays uniformly far from zero, for γ ∈ R.
The mechanical-to-electrical coupling is achieved by a change of variables in the bidomain equations from the current configuration (Eulerian coordinates) to the reference configuration (Lagrangian coordinates), which leads to a conduction term depending on the deformation gradient F. Summarizing, the active strain formulation for the electromechanical activity in the heart is written as follows [30]: where Q T := (0, T ) × Ω R . Here, according to the above discussion, we should take (2.7) The system of equations (2.5) has to be completed with suitable initial conditions for v, w, γ, z 3 and with boundary conditions on v i,e and on the elastic flux a(·, ·, ·, ·). But det(F) = 1, so one can use the approximation tr(F − I) 0, hence, ∇ · φ = tr(F) tr(I) = n. Now, when u denotes the displacement i.e. u = φ(X)−X, the above condition becomes ∇·u = 0, which is the linearized incompressibility condition. We also linearize the flux in (2.6) with respect to F using Taylor series' expansion of Cof(F) about I, given by: and we obtain 8 a(x, γ, F, p) := µFC −1 a (x, γ) − pI.
(2.9) 2.5. The problem to be solved and its weak formulation. For simplicity of notation, we will use Ω and Ω T to denote Ω R and Q T respectively in all what follows, unless otherwise specified. Let us consider the following class of problems: in Ω, for a.e. t ∈ (0, T ), (2.10) in Ω T , (2.14) in Ω T , (2.16) in Ω T . 2 results in the compatibility constraint (2.31) below). The initial data are: For simplicity we take m = 1 in the concentration variable z. The following properties of the 4 model (2.10)-(2.17) and (2.18)-(2.20) are instrumental for the subsequent analysis: is a family of symmetric tensors, uniformly bounded and positive definite: is a family of symmetric matrices, uniformly bounded and positive definite: (A.6) the functions R, G and I ion are given by the kinetics of a general physiological ionic model and it can be verified that the assumptions, stated below, are satisfied by several gating and ionic concentration variables in Beeler-Reuter or Luo-Rudy ionic models. We assume that the function R(v, w) : where α j and β j , j = 1, · · · , k are positive rational functions of exponentials in v such 10 that: dα j dv and dβ j dv are uniformly bounded, (2.21) for some constant C α,β > 0. The function I ion : R × R k × (0, +∞) → R has the general 12 form: where I j ion ∈ C 0 (R × R k ) and satisfies the condition: and I z ion is such that: where Θ,Θ, L belong to C 0 (R, R + ) and C 1,I , . . . , C 5,I are positive constants. Finally the 1 function G is given by: where a 1 , a 2 , a 3 are positive physiological constants that vary from one ion to another. In 3 our case, we only consider z to correspond to the intracellular calcium concentration.

26
Our main result in this paper is the following theorem: uniformly bounded and the functions u(t, ·), v i,e (t, ·) are in H 1 (Ω) 3 and H 1 (Ω) respectively. 6 We also note that passage to the limit in the pressure term p is not straightforward because it is 7 not possible to establish an a priori uniform estimate in L 2 (Ω T ) due to the use of the "artificial 8 compressibility" which utility becomes clearer in the following section.

19
In order to impose the compatibility condition (2.31), we let We observe that span{µ l } l∈N is dense in the space H 1,0 (Ω), given as in Definition 2.1. Fur-20 thermore, we orthonormalize the basis (µ l ) l∈N by the Gram-Schmidt process, and we denote the 21 new basis by (µ l ) l∈N that is orthonormal in L 2 (Ω). For m ≥ 0, we introduce the finite di-

25
We are looking for a discrete solution (3.1) Upon discretization, we obtain a system of ODEs coupled to a system of algebraic equations to . We obtain the for l = 0, · · · , m. Having no initial conditions on the functions u, p, v i and v e in the original 2 problem, we need to supplement our system with initial conditions. We define the functions: and Ω v e,0 dx = 0. We further select u 0 = 0 and an arbitrary p 0 . The 4 initial data of the ODE system are then given by for j = 1, · · · , k. Using the orthonormality of the bases, we can write (3.2) as a system of ordinary differential equations in the coefficients: To be concise, we detail in the following paragraph how the bidomain equations can be treated to 1 obtain the ODE system. We first note that using v m = v i,m − v e,m , we have: Replacing v i,ε,m and v e,ε,m by their expressions as in (3.1), we obtain for l = 0, · · · , m: By the L 2 -orthonormality of the bases, the above equations can be rewritten in the form: where F k , k = i, e assemble all the terms not containing time derivatives. The latter system is equivalent to a system written as: and A = (a lr ) with a lr = Ω ω l µ r . In order to write: v i,m v e,m = M −1 b, we need to prove that the matrix M is invertible. For this sake, we expand it as: It is enough to prove that the matrix N : and ξ e = (ξ e,0 , · · · , ξ e,m ) T ∈ R m+1 . Then Thus the matrix M is positive definite, hence invertible. Consequently, the whole system (3.2) 1 can be written as a system of ordinary differential equations in the form y (t) = f (t, y(t)).  satisfy the following weak formulations for each fixed t, which will be the starting point for deriving 17 a series of a priori estimates: for j = 1, · · · , k. To simplify the notation, we perform the derivations in the following three 1 Lemmata while omitting the subscript ε, m. We start first by obtaining estimates on the gating 2 and concentration variables (w ε,m and z ε,m ) that are needed to prove the uniform bounds. In the 3 following Lemma, we show that the gating variables w j , j = 1, . . . , k satisfy the universal bounds 4 0 ≤ w j ≤ 1.
Proof of (3.14): We have by (2.16) and (2.30) We rewrite it as: After squaring both sides, we obtain: Then we integrate over (0, t),to get: Therefore, by (3.11), (3.20) and (2.24) we find for some constant c 12 > 0. 1 Using the above estimates on z and w, we shall control the L 2 norm of I ion by the L 2 norm of 2 v and this result will be later used to reach a uniform in ε and m estimate on v ε,m .
Using (3.11) and (3.12), one obtains Finally, integrate (3.22) over (0, t) × Ω and use (3.14) along with the condition that z 0 is in L 2 (Ω), 7 to get (3.21). 8 We recall that in order to establish the passage to the limit as m → ∞, we need to bound the 9 solutions of the discrete regularized problem in various Banach spaces, making use of the preceding 10 estimates.

5
The remaining part of this subsection is devoted to proving Theorem 3.1.

15
• and p ε,m p ε weak star in L ∞ (0, T ; L 2 (Ω)) and weakly in L 2 (Ω T ). 16 We also observe that from the sequences U i,ε,m and U e,ε,m introduced in the proof of Lemma 17 3.5, we can extract subsequences such that: Moreover, knowing that ∂ t U i,ε,m and ∂ t U e,ε,m are uniformly bounded in L 2 (0, T ; (H 1 (Ω)) ), we 21 obtain, by compactness and uniqueness of the limit, the following strong convergence: and a.e. in Ω T , 23 U e,ε,m → U e,ε := v ε − εv e,ε in L 2 (Ω T ) and a.e. in Ω T .

17
• u ε u weakly in L 2 (0, T ; H 1 (Ω) 3 ) and ∇u ε ∇u in L 2 (Ω T ) 3×3 . 18 We briefly note that in the distribution sense ε ∂ t u ε , ψ → 0, in D (0, T ). Similarly, ε ∂ t p ε , ψ → 0, in D (0, T ). Due to the "artificial compressibility" used in the proof, we were not able to obtain a bound on 22 ∂ t p ε that is independent of ε except in L 1 (0, T ; L 2 (Ω)), (see remark 3.1), which is not a reflexive 23 space. So in order to pass to the limit in the term involving the pressure, we made a detour by
Consequently, we have, for all q ∈ L 2 (Ω), lim ε→0 Ω (p ε − p)q dx = 0. 1 Therefore, according to all of the preceding convergence results, and repeating some of the 2 arguments of the previous section, we have for all ψ ∈ H 1 (Ω) 3 , ρ in L 2 (Ω), µ ∈ H 1,0 (Ω) (given as  In summary, we consider that in our work, we have paved the way towards addressing the 10 solvability of cardiac electromechanics coupled with physiological ionic models. We used a mathe- insight is needed to mathematically analyze more realistic models. 19