Cardiac memory phenomenon, time-fractional order nonlinear system and bidomain-torso type model in electrocardiology

Abstract: Cardiac memory, also known as the Chatterjee phenomenon, refers to the persistent but reversible T-wave changes on ECG caused by an abnormal electrical activation pattern. After a period of abnormal ventricular activation in which the myocardial repolarization is altered and delayed (such as with artificial pacemakers, tachyarrhythmias with wide QRS complexes or ventricular pre-excitation), the heart remembers and mirrors its repolarization in the direction of the vector of “abnormally” activated QRS complexes. This phenomenon alters patterns of gap junction distribution and generates changes in repolarization seen at the level of ionic-channel, ionic concentrations, ioniccurrent gating and action potential. In this work, we propose a mathematical model of cardiac electrophysiology which takes into account cardiac memory phenomena. The electrical activity in heart through torso, which is dependent on the prior history of accrued heartbeats, is mathematically modeled by a modified bidomain system with time fractional-order dynamics (which are used to describe processes that exhibit memory). This new bidomain system, that I name “ memory bidomain system”, is a degenerate nonlinear coupled system of reaction-diffusion equations in shape of a fractionalorder differential equation coupled with a set of time fractional-order partial differential equations. Cardiac memory is represented via fractional-order capacitor (associate to capacitive current) and fractional-order cellular membrane dynamics. First, mathematical model is introduced. Afterward, results on generalized Gronwall inequality within the framework of coupled fractional differential equations are developed. Next, the existence and uniqueness of solution of state system are proved as well as stability result. Further, some preliminary numerical applications are performed to show that memory reproduced by fractional-order derivatives can play a significant role on key dependent electrical properties including APD, action potential morphology and spontaneous activity.


Modeling motivation
The heart is an electrically controlled mechanical pump which drives blood flow through the circulatory system vessels (through deformation of its walls), where electrical impulses trigger mechanical contraction (of various chambers of heart) and whose dysfunction is incompatible with life. The coordinated contraction of heart and the maintenance of heartbeat are controlled by a complex network of interconnected cardiac cells electrically coupled by gap junctions and voltage-gated ion channels. The evaluation of bioelectrical activity in heart is then a very complex process which uses different phenomenological mechanism and subject to various perturbations, and physiological and pathophysiological variations.
The electrical system of a normal heart is highly organized in a steady rhythmic pattern. This normal heartbeat is called sinus rhythm. Irregular or abnormal heartbeats, called arrhythmias, are caused by a change in propagation and/or formation of electrical impulses, that regulate a steady heartbeat, causing a heartbeat that is too fast or too slow, that can remain stable or become chaotic (irregular and disorganized). Many times, arrhythmias are harmless and can occur in healthy people without heart disease; however, some of these rhythms can be serious and require special and efficiency treatments. Fibrillation is one type of arrhythmia and is considered the most serious cardiac rhythm disturbance. It occurs when the heart beats with rapid, erratic electrical impulses (highly disorganized almost chaotic activation). This leads to quivering (or fibrillation) of heart chambers rather than normal contraction, which then leads the heart to lose its ability to pump enough blood through circulatory system. The treatment therapy of these diseases, when it becomes troublesome or when it can present a danger, often uses electrical impulses to stabilize cardiac function and restore the sinus rhythm, by implanting the patients with active cardiac devices (electrotherapy). For example, in case of cardiac rhythms that are too slow, the devices transmit electronic impulses and ensure that periodic contractions of heart are maintained at a hemodynamically sufficient rate; and in the case of a fast or irregular heart rate, the devices monitor heart rate and, if needed, treat episodes of tachyarrhythmia (including tachycardia and/or fibrillation) by transmitting automatically impulses to either give defibrillation shocks or cause overstimulation (via an ICDs * ) or synchronize contraction of left and right ventricles.
After cessation of a transient period of abnormal ventricular activation (arrhythmia or pacing) in which the myocardial repolarization is altered and delayed (such as with artificial pacemakers [82], tachyarrhythmias with wide QRS complexes, intermittent left bundle branch block or ventricular pre-excitation observed in Wolff-Parkinson-White syndrome [45]), the heart remembers and mirrors its repolarization in the direction of vector of "abnormally" activated QRS complexes [66]. This remodeling of electrical properties of myocardium is characterized by persistent but reversible T-wave changes on the surface electrocardiogram (ECG). The scope, significance and direction of T-wave deviation depend on duration and direction of abnormal electrical activation. Moreover, these changes are often confused with pathological conditions manifesting with T-wave deviations, such as (acute) myocardial ischemia or infarction. This cumulative and complex phenomenon is named cardiac memory and can persist up for several weeks after normal ventricular conduction is restored. Heart is considered as network of cardiac oscillators communicating via gap junctions between neighboring cells and through voltage gated ion channels (that are activated by changes in electrical membrane potential near channel): this phenomenon alters patterns of gap junction distribution and generates changes in repolarization seen at the level of ionic-channel, ionic concentrations, ionic-current gating and action potential.
The existence of cardiac memory has been known for many years and resulted in a large number of publications, see for example [23, 37, 44, 58-60, 67, 72, 73, 75] and the references therein. Yet despite all this, this phenomemon is under-recognized and there is still limited information regarding its physiological significance and practical clinical implications because this phenomenon is considered to be a relatively benign pathophysiologic finding. Unfortunately, the work of past few years has shown that despite the mostly benign nature of cardiac memory in healthy individuals under some conditions, it can be the trigger of more complex arrhythmias and then requires emergency care of the patient. Other works estimate that clinically administered antiarrhythmic drugs alter expression of cardiac memory and that generate changes in repolarization could in turn alter the effects of these drugs (see e.g. [4,63]). Moreover, cardiac memory may also lead to unnecessary and invasive diagnostic investigation that put patients under unnecessary risks (see e.g., [71]).
Then, recognition of cardiac memory as a serious potential cause of T-wave changes and the analysis of T-wave morphology (throughout the ECG during narrow-and wide-QRS rhythms) are critical to help differentiate T-wave changes due to myocardial ischemia from those induced by cardiac memory, and consequently sustainably establish the clinical relevance of this phenomenon, facilitates diagnosis and increases efficiency of cardiac disorders treatment.
Consequently, this has greatly emphasized the need for model and methodologies capable of predicting and understanding the dynamic mechanisms of different sources of electrical instability in heart like cardiac action potential (AP) repolarization alternans which is influenced by memory. At cellular level, alternans is generally manifests as cyclic, beat-to-beat alternations between long and short action potential and/or intracellular Ca transient, and is frequently associated with the development of ventricular tachycardia and fibrillation. It is generally agreed that disturbances of bi-directional (mutual) relationship between transmembrane potential (or membrane voltage) φ and Ca-sensitive ionic currents play a key role in generation of alternans. Because membrane voltage φ is strongly affected by Ca-sensitive ionic currents and intracellular Ca loading in turn is strongly influenced by φ-dependent ionic currents. This complex bidirectional coupling influences and controls the amplitude and duration of action potentials (APD) through various time-and voltage-dependent ionic currents. The classical bidomain system is commonly used for modeling propagation of electrophysiological waves in cardiac tissue. Motivated by above discussions, to take into account the critical effect of memory in propagation of electrophysiological waves, together with other critical cardiac material parameters, we propose and analyze a new bidomain model, that I name "memory bidomain system" by incorporating memory effects. In next section, we shall present derivation of this memory bidomain model.

Modeling and formulation of the problem
Mathematical and computational cardiac electrophysiological modeling is now an important field in applied mathematics. Indeed, nowadays, heart and cardiovascular diseases are still the leading cause of death and disability all over the world. That is why we need to improve our knowledge about heart behavior, and more particularly about its electrical behavior. Consequently we want strong methods to compute electrical fluctuations in the myocardium to prevent cardiac disorders (as arrhythmia). Tissuelevel cardiac electrophysiology, which can provide a bridge between electrophysiological cell models at smaller scales, and tissue mechanics, metabolism and blood flow at larger scales, is usually modeled using the coupled bidomain equations, originally derived in [78], which represent a homogenization of the intracellular and extracellular medium, where electrical currents are governed by Ohm's law (see also e.g. [46] for a review and an introduction to this field). The model was modified and extended to include heart tissue surrounded by a conductive bath or a conductive body (see e.g. [64] and [74]). From mathematical viewpoint, the classical bidomain system is commonly formulated in terms of intracellular and extracellular electrical potentials of anisotropic cardiac tissue (macroscale), φ i and φ e , (or, equivalently, extracellular potential φ e and the transmembrane voltage φ = ϕ i − ϕ e ) coupled with cellular state variables u describing cellular membrane dynamics and torso potential state variable φ s . This is a system of non-linear partial differential equations (PDEs) coupled with ordinary differential equations (ODEs), in the physical region Ω (occupied by excitable cardiac tissue, which is an open, bounded, and connected subset of d-dimensional Euclidean space R d , d ≤ 3). The PDEs describe the propagation of electrical potentials and ODEs describe the electrochemical processes.
Cardiac memory can affect considerably the resulting electrical activity in heart and thus the cardiac disorders therapeutic treatment. It is then necessary to introduce the impact of memory on dynamical behaviors of such a system. Memory terms can cause dynamical instabilities (as Alternans) and give rise to highly complex behavior including oscillations and chaos. In order to take into account the influence of cardiac memory and inward movement of u into the cell which prolongs depolarization phase of action potential, we propose a new bidomain model. In this new model, classical bidomain model has been modified via time fractional-order dynamical system, arising due to cellular heterogeneity, which are used to describe processes that exhibit memory. More precisely, the derived system with memory (or history), is a nonlinear coupled reaction-diffusion model in shape of a set of time fractional-order differential equations (FDE) coupled with a set of time fractional-order partial differential equations, in the torso-heart's spatial domain Ω (Figure 1) which is a bounded open subset with a sufficiently regular boundary ∂Ω, and during the final fixed time horizon T > 0, as follows  Here, ∂ α 0 + and ∂ β 0 + denote the forward Caputo fractional derivatives with α and β be real values in ]0, 1] and the unknowns are the potentials ϕ i , ϕ e , ϕ s and a single ionic variable u (e.g. gating variable, concentration, etc.).
The heart's spatial domain is represented by Ω H which is a bounded open subset, and by Γ H = ∂Ω H we denote its piecewise smooth boundary. A distinction is made between the intracellular and extracellular tissues which are separated by the cardiac cellular membrane. The surrounding tissue within the thorax is modeled by a volume conduction Ω B with a piecewise smooth boundary Γ B = Γ H ∪ Γ T where Γ T is the thorax surface. The whole domain is denoted by the Ω = Ω H ∪ Ω B . In Ω H , the transmembrane potential is φ = ϕ i − ϕ e where ϕ e and ϕ i are the transmembrane, extracellular and intracellular potentials, respectively, and in Ω B , ϕ s is thoracic medium electric potential. The parameter c α is c α = κC α > 0, where C α is the membrane capacitance per unit area and κ is the surface area-tovolume ratio (homogenization parameter). The membrane is assumed to be passive, so the capacitance C α can be assumed to be not a function of state variables. Classically this membrane is assimilate to a simple parallel resistor-capacitor circuit. However various studies showed clearly that a passive membrane may be more appropriately modeled with a non-ideal capacitor, in which the current-voltage relationship is given by a fractional-order derivative (see e.g [84]). So, according to the passivity of tissue we can assimilate this membrane to electrical circuit with a resistor associate to the ionic current (I ion ) and a capacitor associate to the capacitive current, I α c = c α ∂ α 0 + φ, in parallel, with α usually ranging from 0.5 to 1 ( Figure 2). Moreover, since the electrical restitution curve (ERC) † is affected by the action potential history through ionic memory, we have represented the memory via ionic variable u by a time fractional-order dynamic term ∂ β 0 + u. Figure 2. Modeling of the membrane as resistor and non-ideal capacitor coupled in parallel (where 0 < α ≤ 1 is the fractional order of capacitor).
The coupling between equations in Ω H and equation in Ω B (in systems (1.1) and (1.2)) is operate at the heart-torso interface Γ H . The continuity conditions at the interface can be given by where n is the outward normal to Γ H . The tensors K i and K e are the conductivity tensors describing anisotropic intracellular and extracellular conductive media, and tensor K s represents the conductivity tensor of thoracic medium. The electrophysiological ionic state u describes a cumulative way of effects of ion transport through cell membranes (which describes e.g., the dynamics of ion-channel and ion concentrations in different cellular compartments). The operator I is equal to κI ion , where the nonlinear operator I ion describes the sum of transmembrane ionic currents across cell membrane with u. The nonlinear operator G is representing the ionic activity in myocardium. Functional forms for I and G are determined by an electrophysiological cell model. The source terms are I i = κ f i , I e = κ f e and I = I i + I e , where f i and f e describe intracellular and extracellular stimulation currents, respectively. To close the system, we impose the following initial and boundary conditions initial conditions where n T is the outward normal to Γ T . In absence of a grounded electrode, the bidomain equations are a naturally singular problem since ϕ e and ϕ s , in system (1.2), only appear in the equations and boundary conditions through their gradients. Moreover, the states ϕ e and ϕ s are only defined up to the same constant. Such problems have compatibility conditions determining whether there are any solution to the PDEs. This is easily found by integrating the second and third equations of (1.2) over the domain and using the divergence theorem with boundary conditions. Then the following conservation of the total current is derived (a.e. in (0, T )) (1.5) Consequently, we must choose current I such that the compatibility condition (1.5) is satisfied. For the choose of intracellular stimulus, it is natural and usual (as is common) to take I i = −I e (denoted in the sequel by f H ) i.e., setting the total current I to zero. It is clear that this condition does not correspond to zero extracellular stimulus, but that the extracellular stimulus current and the intracellular stimulus current are equal in magnitude and opposite in direction. Moreover, the functions ϕ e and ϕ s are defined within a class of equivalence, regardless of the same time-dependent function. This function can be fixed, for example, by setting the condition, (a.e. in (0, T)) where χ Ω (i) H is the characteristic function of set Ω (i) H , i = 1, 2. The support regions Ω (1) H and Ω (2) H can be considered to represent an anode (positive electrode) and a cathode (negative electrode) respectively.
In recent years, various problems concerning biological rhythmic phenomena and memory processes which can be included in a cardiac model in many ways (via delay model or time fractional dynamical model), have been studied. Via delayed system we can cited e.g., [11,17,20,40,41,43,69] and the references therein. For problems associated with bidomain models with time-delay, the literature is limited, to our knowledge, to [8,9,27]. In these references, in order to take into account the influence of disturbance in data and the time-varying delays on propagation of electrophysiological waves in heart, the authors have developed a new mathematical model and have considered the theoretical analysis as well as numerical simulations (with real data) and validation of developed model. Via time fractional dynamical model, the literature is also limited, we can cite e.g. [26,31]. In [26], the authors study, using a minimal cardiac cell model, the effects of a fractional-order time diffusion for the voltage and for the ionic current gating. They have shown numerically the interest of modeling memory through fractional-order and that with the model it is able to analyze the influence of memory on some electrical properties as spontaneous activity and alternans. Concerning problems associated with classical bidomain models various methods and technique, as evolution variational inequalities approach, semi-group theory, Faedo-Galerkin method and others, the studies of well-posedness of solutions have been derived in the literature (see e.g., [12,16,18,19,24,80] and the references therein); for development of multiscale mathematical and computational modeling of bioelectrical activity in myocardial tissue and the formation of cardiac disorders (as arrhythmias), and their numerical simulations, which are based on methods as finite difference method, finite element method or lattice Boltzmann method, have been receiving a significant amount of attention (see e.g., [5, 6, 21, 28-30, 35, 38, 46, 51, 65, 79, 81] and the references therein).
The rest of paper is organized as follows. In next section, we give some preliminaries results useful in the sequel. In Section 3, some preliminaries results concerning fractional calculus are given and results on generalized Gronwall inequality to a coupled fractional differential equations are developed. In Section 4 we shall prove the existence, stability and uniqueness of weak solutions of derived model, under some hypotheses for data and some regularity of nonlinear operators. Numerical experiments, using a modified Lattice Boltzmann Method for numerical simulations of the derived memory bidomain systems, are described in Section 5. In Section 6, conclusions are discussed.

Assumptions, notations and some fundamental inequalities
Let ⊂ IR m , m ≥ 1, be an open and bounded set with a smooth boundary and T = × (0, T ). We use the standard notation for Sobolev spaces (see [1]), denoting the norm of W q,p ( ) (q ∈ IN, p ∈ [1, ∞]) by . W q,p . In the special case p = 2, we use H q ( ) instead of W q,2 ( ). The duality pairing of a Banach space X with its dual space X is given by ., . X ,X . For a Hilbert space Y the inner product is denoted by (., .) Y and the inner product in L 2 (Ω) is denoted by (., .). For any pair of real numbers r, s ≥ 0, we introduce the Sobolev space H r,s ( T ) defined by H r,s ( T ) = L 2 (0, T ; H r ( )) ∩ H s (0, T ; L 2 ( )), which is a Hilbert space normed by v 2 L 2 (0,T ;H r ( )) + v 2 H s (0,T ;L 2 ( )) 1/2 , where H s (0, T ; L 2 ( )) denotes the Sobolev space of order s of functions defined on (0, T ) and taking values in L 2 ( ), and defined (see [52]) by H s (0, T ; L 2 ( )) = [H q (0, T, L 2 ( )), L 2 ( T )] θ , for s = (1 − θ)q with θ ∈ (0, 1) and q ∈ IN, and H q (0, T ; L 2 ( )) = v ∈ L 2 ( T )| ∂ j v ∂t j ∈ L 2 ( T ), for 1 ≤ j ≤ q . For a given Banach space X, with norm . X , of functions integrable on , we define its subspace X| IR = u ∈ X, u = 0 that is a Banach space with norm . X , and we denote by [u] the projection of u ∈ X on X| IR such that udx (with mes( ) standing for Lebesgue measure of ). (ii) (Gagliardo-Nirenberg inequalities) There exists C > 0 such that where 0 ≤ θ < 1 and p = 2m m−2θq (with the exception that if q − m/2 is a nonnegative integer, then θ is restricted to 0).
Finally, we introduce the spaces: • H H = L 2 (Ω H ), V H = H 1 (Ω H ), V = H 1 (Ω) (endowed with their usual norms) and U H = V H | IR , We will denote by V H (resp. U H ) the dual of V H (resp. of U H ). We have the following continuous  Our study involves the following fundamental inequalities, which are repeated here for review: (i) Hölder's inequality: (ii) Young's inequality (∀a, b > 0 and > 0): ab ≤ p a p + −q/p q b q , f or p, q ∈]1, +∞[ and 1 Finally, we denote by L(A; B) the set of linear and continuous operators from a vectorial space A into a vectorial space B, and by R * the adjoint operator to a linear operator R between Banach spaces.
From now on, we assume that the following assumptions hold for the nonlinear operators and tensor functions appearing in our model.
The operators I and G which describe electrophysiological behavior of the system can be taken as follows (affine functions with respect to u) where is a sufficiently regular function and I 0 (x, t; φ) = g 1 (x; φ) + g 2 (x, t; φ) with g 1 is an increasing function on φ. Moreover, the operators I 0 , I 1 and I 2 appearing in I and G, are supposed to satisfy the following assumptions.
(H2) p The operators I 0 , I 1 and I 2 are Carathéodory functions from (Ω × IR) × IR into IR and continuous on φ (as in [12]). Furthermore, for some p ≥ 2 if d = 2 and p ∈ [2, 6] if d = 3, the following requirements hold (i) there exist constants β i ≥ 0 (i = 1, . . . , 10) such that for any v ∈ IR where E g is the primitive of g 1 .
such that for any (v, w) ∈ IR 2 : In order to assure the uniqueness of solution we assume that (H3) The Nemytskii operators I and G satisfy Carathéodory conditions and there exists some µ > 0 such the operator F µ : IR 2 → IR 2 defined by satisfies a one-sided Lipschitz condition (see e.g., [14,70]): there exists a constant C L > 0 such that Lemma 2.2. ( [9]) Assume that F µ is differentiable with respect to (φ, u) and denote by λ 1 (φ, u) ≤ λ 2 (φ, u) the eigenvalues of symmetrical part of Jacobian matrix ∇F µ (φ, u): If there exists a constant C F independent of φ and u such as: then F µ satisfies the hypothesis (H3).
The function I 0 is defined by a cubic reaction term of the form , and the functions I 1 and I 2 are given by (a) for RM type model : where b i ∈ W 1,∞ (Q), i = 1, 3, are sufficiently regular functions from Q into IR +, * and r ∈ [0, 1]. We obtain easily the following Lemma.

Fractional calculus and a generalized Gronwall's inequality
Fractional differential equations have been studied by many investigations in recent years and the idea of defining a derivative of fractional order (non-integer order) dates back to Leibnitz [49]. Since 19th century, different authors have considered this problem e.g., Riemann, Liouville, Hadmard, Hardy, Littlewood and Caputo among others. Fractional integrals and derivatives have proved to be useful in real applications, since they arise naturally in many biological phenomena such as viscoelasticity, neurobiology and chaotic systems, see for instance [3,34,36,55,62,76,83].
The classical and the most used form of fractional calculus is given by the Riemann-Liouville and Caputo derivatives. In contrast to this nonlocal Riemann-Liouville derivative operator, when solving differential equations, is the use of Caputo fractional derivative [22] in which it is not necessary to define the fractional order initial conditions.

Definitions and basic results
The object of this section is to give a brief introduction to some definitions and basic results in fractional calculus in the Riemann-Liouville sense and Caputo sense. Let γ ∈]0, 1] and X be a Banach space, we start from a formal level and assume the given functions Definition 3.1. The forward and backward Riemann-Liouville fractional integrals of fractional order γ on (a, T ) are defined, respectively, by (t ∈ (a, T )) The basic equality for the fractional integral is (from Fubini's Theorem and the relationship between Γ-function and β-function) and holds for a L p -function f (1 ≤ p ≤ ∞).
Definition 3.2. The forward Riemann-Liouville and Caputo derivatives of fractional order γ on (a, T ) are defined, respectively, by (t ∈ (a, T )) From (3.2) we can deduce the following relation between fractional integral and Caputo derivative The backward Riemann-Liouville and Caputo derivatives of fractional order γ, on (a, T ) are defined, respectively, by (t ∈ (a, T )) By substitution a → −∞, the following definition is obtained.
Definition 3.4. The Liouville-Weyl fractional integral and the Caputo fractional derivative on the real axis are defined, respectively, by (t ∈ (−∞, T ))

It is possible to show that the difference between Riemann-Liouville and Caputo fractional
derivatives depends only on the values of f on endpoints a and T . More precisely, for From [42], we have the following results Lemma 3.1. (Continuity properties of fractional integral in L p spaces on (a, T )) The fractional integral I γ a + is a continuous operator from (i) L p (a, T ) into L p (a, T ), for any p ≥ 1, From Lemma 3.1 and (3.4) we can deduce the following corollary.
Lemma 3.2. Let X be a Banach space and γ ∈]0, 1]. Suppose the Caputo derivative ∂ γ a + f ∈ L p (a, T ; X) and p > 1 γ , then f ∈ C 0,γ−1/p ([a, T ]; X). We also need for our purposes the fractional integration by parts in the formulas (see for instance [2,47,86]) Lemma 3.4. Let 0 < γ ≤ 1, and g be L p -function on (a, T ) (for p ≥ 1) and f be absolutely continuous function on [a, T ]. Then From [85], we can deduce the following Lemma.
If in addition h is a nondecreasing function on (0, T ), then The used function E θ 1 ,θ 2 is the classical two-parametric Mittag-Leffler function (usually denoted by The function E θ 1 ,θ 2 is an entire function of the variable z for any θ 1 , θ 2 ∈ l C, Re(θ 1 ) > 0. Next we declare a compactness theorem in Hilbert spaces. Assume that X 0 , X 1 and X are Hilbert spaces with X 0 → X → X 1 being continuous and X 0 → X is compact.
The Fourier transform of f :

s)ds and we have
Lemma (see e.g., [68]) Lemma 3.6. Let γ ∈ (0, 1) and f be a L 1 -function in IR with compact support. Then We define the Hilbert space W γ (IR; X 0 , X 1 ), for a given γ > 0, by From Lemma 3.6 and similar arguments to drive Theorem 2.2 in [77], we have the following compactness result.
Theorem 3.1. Let X 0 , X 1 and X be Hilbert spaces with the injection (3.9). Then for any bounded set K and any γ > 0, the injection of W γ K (IR; X 0 , X 1 ) into L 2 (IR; X) is compact.

A generalized Gronwall inequality to a coupled fractional differential equation
In this section we present a new generalized Gronwall inequality with singularity in the context of a coupled fractional differential equations.

11)
where The function E .,. is the Mittag-Leffler function which is defined by (3.8).
Then, if f and g are nonnegative functions locally integrable, and satisfy the inequalities (3.10), we have the following estimate Proof. From Hölder's inequality, we obtain (for i = 1, 2) where (3.20) Thus (3.21) This completes the proof. In the sequel we will always denote C some positive constant which may be different at each occurrence.

Well-Posedness of the memory bidomain-torso system
In this section, we prove the existence and uniqueness of weak solution to problem (1.1), under Lipschitz and boundedness assumptions on the non-linear operators.

Variational formulation and preliminary results
First we define the state ψ, as the extracellular cardiac potential in Ω H , and the thoracic potential in Ω B , by ψ = ϕ e in Ω H and ψ = ϕ s in Ω B . (4.1) Then the potential φ can be written as φ = ϕ i − ψ| Ω H . From the interface coupling condition (1.3), it follows that if ϕ e ∈ H 1 (Ω H ) and ϕ s ∈ H 1 (Ω B ) then ψ ∈ H 1 (Ω). Similarly, we introduce the corresponding conductivity tensor K g ∈ W 1,∞ (Ω) as follows : We now define the following forms Proof. (i) and (ii) are easily obtained providing that properties of tensors K k , for k = i, e, s, g, and (2.2) are satisfied.
We can now define the weak solution to bidomain-thorax coupled model problem (1.1).

Existence, uniqueness and regularity results
The purpose of this section is to prove the well-posedness of the degenerate bidomain-thorax coupled model (1.1). To overcome this issue, we introduce a specific non-degenerate approximate to system (1.1).

Non degenerated problem study
Let > 0 be a small positive number. Hence, our approximation system read as (since with the initial condition where ψ and the initial condition ψ 0 are defined as in (4.1).
The following results concern the existence and regularity of the weak solution to problem (4.4) (i.e. of (4.5)).
Proof. To establish the existence result of weak solution to system (4.4) we apply the Faedo-Galerkin method, derive a priori estimates, and then pass to the limit in the approximate solutions using compactness arguments. We approximate the system equations by projecting them onto finite n dimensional subspaces, then we take the limit in n. For this, let (w k ) k≥1 be a Hilbert basis and orthogonal in L 2 of V H , ( f k ) k≥1 be a Hilbert basis of U H , (g k ) k≥1 be a Hilbert basis of We denote byf k an extension of f k in H 1 (Ω),g k an extension of g k in H 1 (Ω) withg k | Ω H = 0 and by (e k ) k≥1 a Galerkin basis of U HB which is defined as e 2k =g k and e 2k−1 =f k (for all k > 1 ). For all n ∈ IN * , we denote by W i,n = span(w 1 , · · · , w n ), W e,n = span( f 1 , · · · , f n ), W s,n = span(g 1 , · · · , g n ) and W g,n = span(e 1 , · · · , e 2n ) the spaces generated, respectively, by (w k ) n≥k≥1 , ( f k ) n≥k≥1 , (g k ) n≥k≥1 and (e k ) 2n≥k≥1 , and we introduce the orthogonal projector L j,n on the spaces W j,n (for j = i, e, s, g).
N.B. The constants which will be introduced in the sequel, are independent of n and unless otherwise specified.

(4.35)
Finally for the external forcing we have (using Hölder inequality) According to (4.35), (4.36), (4.19), (4.21), and the fact, from (4.9) and ∈ (0, 1), the quantities ψ n (0), φ n (0) and ϕ i,n (0) are uniformly bounded in H 1 , we can derive from (4.25) the following estimate (for all t ∈ (0, T ))  We can deduce that (from (4.27) and (4.37)) In order to show that the local solution can be extended to the whole time interval (0; T ), we assume that we have already defined a solution of (4.12) on [0, T k ] and we shall define the local solution on [T k , T k+1 ] (where 0 < T k+1 − T k is small enough) by using the obtained a priori estimates and the fractional derivative ∂ γ T + k (for γ = α or β) with beginning point T k . Consequently, by iteration process, we thus obtain that Faedo-Galerkin solutions are well-defined on (0, T ). So, we omit the details.
We can now show the well-posedness of memory bidomain-torso (degenerate) problem (1.1).
(i) If α = β, we have the following relation  (ii) If α < β and if the following assumption holds the relation (4.46) is also obtained.
Remark 4.1. In the previous theoretical work, our main results investigate the well-posedness of system (1.1) with an abstract class of ionic models, including some classical models as Rogers-McCulloch, Fitz-Hugh-Nagumo and Aliev-Panfilov. We can consider other ionic model type including Mitchell-Schaeffer model,see [56] (which has a slightly different structure). This two-variable model can be defined with operators I and G as (in which the state variables are dimensionless and scaled) where 0 < φ gate < 1, τ in , τ out , τ open and τ close are given positive constants. The cubic function φ 2 (1 − φ) describes the voltage dependence of the inward current. This model is well-known to be valid under the assumption 0 < τ in τ out min(τ open , τ close ). In order to guarantee the well-posedness of system (1.1) with Mitchell-Schaeffer ionical model, we can use the following regularized version of ionic operator G where the differentiable function 0 ≤ h ζ ≤ 1 is given by with ζ a positive parameter. The operator G ζ (φ, u) can be written as G ζ (φ, u) = I 2,ζ (φ) + ζ (φ)u where (4.60) According to the definition of tanh, we can deduce that lim u). The regularized Mitchell-Schaeffer model has a slightly different structure compared to models in (2.3) because in this model, ζ depend on φ through the function h ζ . Since h ζ is sufficiently regular, the arguments of this paper can be adapted with some slight necessary modifications to analyze the well-posedness of this regularized Mitchell-Schaeffer ionical model. In more general, the study developed in this paper remains valid (with some necessary modifications) if we consider the operator G in the form of G(.; φ, u) = I 2 (.; φ) + (.; φ)u (i.e. a general form of Hodgkin-Huxley model including Beeler-Reuter and Luo-Rudy ionic models described by continuous or regularized discontinuous functions, see [7,53,54]) with G Carathéodory function from (Ω × IR) × IR 2 into IR and locally Lipschitz continuous function on (φ, u) and, I 2 and sufficiently regulars.

Numerical applications
In this section, we shall present some preliminary numerical simulations of the problem (1.1) to illustrate our result with showing two-dimensional simulation. Two of the most used ionical models, FitzHugh-Nagumo model [39] and Mitchell-Shaeffer model [56], are considered. This numerical study investigates the impact of fractional-order φ dynamics and fractional-order u dynamics (i.e. the influence of the capacitor and the ionic variable memory) on key rate dependent electrical properties including APD, the amplitude of the action potential and spontaneous activity.
Although part of the theoretical analysis was limited to the case β ≥ α, in order to show the impact of β on the evolution of system, we will also present some results for α fixed to 1 and β variable. We will explore the behavior of the two models when first they are stimulated by a brief current and second when they are stimulated twice in sequence.
In these two situations we fix first β at 1 and we vary α, and in a second time we fix α at 1 and we vary β to analyze the impact of each time-fractional order derivative on system behavior, one after the other.
In order to solve numerically system (1.1) we have analyzed different approximation methods of Caputo fractional derivative ∂ γ 0 + , γ ∈]0, 1]. In this paper, we consider the following approximation (for sufficiently regular function ψ at time t ∈]t n−1 , t n ] (n = 1, 2, . . .), with t k = k∆t for k = 0, 1, . . . and given time step ∆t) where It is clear that I 1 memory (ψ history ), for γ = 1, is a null function. The state function ψ history is corresponding to the prior history of ψ.
Note that the previous fractional-order dynamic is the sum of first-order dynamics and an operator I γ memory representing hypothetical memory effects.
The term I γ memory depends on the prior history of accrued heartbeats and characterize the influence of capacitive/ionic memory on the dynamical of system. This memory operator I γ memory can be interpreted as a new force which is added over time to the external applied forcing. So, for example, this adding force can invert at some time points the sign of the external applied stimulus ‡ and then generates changes in myocardial polarization.
Here, we perform various numerical experiments to investigate the impact of different values of φ fractional-order α and u fractional-order β across the time in (0, T ) and we present some numerical ‡ Thus the depolarizing effect of the stimulus becomes repolarizing and vice versa results corresponding to the middle point of heart domain Ω H for φ, ϕ e and u, and in the middle point of torso domain Ω B for ϕ s . We assume that the domain heart-torso is in a square region Ω = [−L T /2, −L T /2] × [−L T /2, −L T /2] and the domain Ω H corresponds to an approximate shape of heart, where L T is the torso length and L H is the characteristic length of Ω H (Figure 3).
The boundary of domain Ω H is defined by (for (x, y) ∈ Γ H ) :  We use some biophysical parameters which are presented in Table 1 and we fix L H = 15cm and In those both applications, we impose the following initial conditions Finally, we consider T = 1, K i = K e = 0.003s.cm −1 I d (with I d identity matrix) and we define (as in [28]) the external stimulus f H in short time a i < t < a i + T act , i = 0, 1 with a 0 = 0 (electroshock for example), as where I app is the amplitude of external applied stimulus with I i,app = b i 10 4 , and the functions χ H , χ [a i ,T act +a i ] and χ prop are defined by The last characteristic is associated to wave propagation as a diagonal which evolves from to leftbottom side to right-top side of Ω H . The function Φ, which corresponds to the shape of electrical wave, is defined by To perform the numerical simulations, we have developed a numerical scheme reliable, efficient, stable and easy to implement in context of such "memory bidomain systems" by generalizing the modified Lattice Boltzmann Method introduced in previous works [9,27,28]. Nota Bene: In all figures first-order dynamic (corresponding to the case α = β = 1) is shown in black and the last curve drawn is in dotted line.

FitzHugh-Nagumo ionical model
The first ionical model we choose is the simplest FitzHugh-Nagumo phenomenological model (see [39]). This model can be defined with the operators I and G as: and For this model, we investigate only the impact of membrane capacitive memory on action potential properties in absence of ionic variable memory (i.e., we fix β to 1 and we vary α). In Figures 4-6 (response to single stimulus) and in Figures 7-9 (response to two stimuli in sequence) we plot the time evolution of φ, ϕ e and u, at midpoint of heart, for different values of α with β = 1. We observe first that we have always the same kind of behavior. Second, we find that the amplitude varies depending on α values. For state ϕ e , the effect fractional-order dynamic is negligible except for the extremum of ϕ e (as the value of α becomes smaller, the value of the minimum point becomes smaller too).

Mitchell-Schaeffer ionical model
The second ionic model corresponds to Mitchell-Shaeffer biophysical ionic model (see [56]). This two-variable model can be defined with operators I and G as: These operators depend on the change-over voltage φ gate , the resting potential φ min , the maximum potential φ max , and on times constants τ in , τ out , τ open and τ close , the two times τ open and τ close , respectively controlling the durations of the action potential and of the recovery phase, and the two times τ in and τ out , respectively controlling the length of depolarization and repolarization phases. These constants are such that τ in < τ out < min(τ open , τ close ). First, we investigate the influence of membrane capacitive memory on action potential properties in absence of ionic variable memory (i.e., we fix β to 1 and we vary α). In Figures 10-13 (response to a single stimulus) and in Figures 14-17 (response to two stimuli in sequence), we plot the time evolution of φ, ϕ e and u, at midpoint of heart and ϕ s at at midpoint of torso, for different values of α with β = 1. We find first that following the stimuli cessation, spontaneous action potentials are triggered and persist for the duration of simulation. Second, the spontaneous action potential cycle length decreases as fractional-order α decreases, and we notice therefore an increase in the rate of spontaneous activity (when α decreases). We observe the same kind of behaviors for the other state variables than for φ, with an acceleration of spontaneous activity when α decreases.  Figure 17. Response to a multiple stimuli (two stimuli in sequence). Time evolution of the torso potential ϕ s at midpoint of torso with β = 1 and α = 1, 0.925, 0.9, 0.85.
We next investigate the influence of ionic variable memory on action potential properties in absence of membrane capacitive memory. In Figures 18-21 (single stimulus) and in Figures 22-25 (two stimuli in sequence), we plot the time evolution of φ, ϕ e and u, at midpoint of heart and ϕ s at at midpoint of torso, for different values of β with α = 1. We observe that the smaller the value of β becomes, the more APD is shortened. Therefore the ionic current memory alters the action potential duration. Moreover from some value of APD we observe the triggering of spontaneous chain activities.  Figure 25. Response to a multiple stimuli (two stimuli in sequence). Time evolution of the torso potential ϕ s at midpoint of torso with α = 1 and β = 1, 0.925, 0.9, 0.85, 0.8.

Conclusions
Modeling and control of electrical cardiac activity represent nowadays a very valuable tool to facilitate diagnosis and maximize the efficiency and safety of treatment for cardiac disorders. In this work, we have developed a new bidomain model of the cardiac electrophysiology which takes into account cardiac memory phenomena, named "memory bidomain model". Cardiac memory is a non-anecdotal phenomena, the impact of which in clinical practice is significant, particularly at the level of diagnostic and decision-making challenges. The derived model is a degenerate nonlinear coupled system of reaction-diffusion equations in shape of a fractional-order ODE coupled with a set of time fractional-order PDEs. Cardiac memory is represented via fractional-order capacitor and fractional-order cellular membrane dynamics. The existence of a weak solution as well as regularity and Lipschitz continuity of map solution results are established, with an abstract class of ionic models, including some classical models as Rogers-McCulloch, Fitz-Hugh-Nagumo and Aliev-Panfilov. This theoretical analysis is completed by numerical results that validate the interest of developed model. For this, some preliminary numerical simulations are performed to illustrate our results by comparing dynamic restitution curves obtained by numerically solving the first-order model and time fractional-order models (which reproduce critical memory effects). FitzHugh-Nagumo model and Mitchell-Shaeffer model are considered in these numerical applications. We observe that memory can alter the key dependent electrical properties including APD, action potential morphology and spontaneous activity.
These first observations demonstrate that memory effects can play a significant role in cardiac electrical dynamics (whether normal or disease states) and then motivates the continuation and deepening of these studies.
For predicting and acting on phenomena and processes occurring inside and surrounding cardiac medium, it would be interesting to study control and regulation problems in order to determine § the best optimal prognostic values of sources in presence of disturbances and pacing history, using the approach developed in [8,10,13]. The resolution of this type of problem makes it possible in practice to improve the stimuli applied during e.g., defibrillation sessions. This stimulus must be adjusted according to the patient's metabolism and type of cardiac problem (cardiac arrest rhythms, bradycardia, tachycardia, etc.).