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

1.   Introduction and mathematical setting of the problem

1.1.   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.

^* The so-called implantable cardioverter defibrillators

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,59,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) $\phi$ and Ca-sensitive ionic currents play a key role in generation of alternans. Because membrane voltage $\phi$ is strongly affected by Ca-sensitive ionic currents and intracellular Ca loading in turn is strongly influenced by $\phi$-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.

1.2.   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). Tissue-level 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), $\phi_{i}$ and $\phi_{e}$, (or, equivalently, extracellular potential $\phi_{e}$ and the transmembrane voltage $\phi = \varphi_{i}-\varphi_{e}$) coupled with cellular state variables $u$ describing cellular membrane dynamics and torso potential state variable $\phi_{s}$. This is a system of non-linear partial differential equations (PDEs) coupled with ordinary differential equations (ODEs), in the physical region $\Omega$ (occupied by excitable cardiac tissue, which is an open, bounded, and connected subset of $d$-dimensional Euclidean space $\mathbf{R}^d$, $d\leq 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 $\Omega$ (Figure 1) which is a bounded open subset with a sufficiently regular boundary $\partial \Omega$, and during the final fixed time horizon $T > 0$, as follows

or equivalently

Here, $\partial^{\alpha}_{0^+}$ and $\partial^{\beta}_{0^+}$ denote the forward Caputo fractional derivatives with $\alpha$ and $\beta$ be real values in $]0, 1]$ and the unknowns are the potentials $\varphi_i$, $\varphi_{e}$, $\varphi_s$ and a single ionic variable $u$ (e.g. gating variable, concentration, etc.).

The heart's spatial domain is represented by $\Omega_{\rm H}$ which is a bounded open subset, and by $\Gamma_{\rm H} = \partial \Omega_{\rm 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 $\Omega_{\rm B}$ with a piecewise smooth boundary $\Gamma_{\rm B} = \Gamma_{\rm H} \cup \Gamma_T$ where $\Gamma_T$ is the thorax surface. The whole domain is denoted by the $\Omega = \overline{\Omega}_{\rm H}\cup {\Omega}_{\rm B}$. In $\Omega_{\rm H}$, the transmembrane potential is $\phi = \varphi_{i}-\varphi_{e}$ where $\varphi_{e}$ and $\varphi_{i}$ are the transmembrane, extracellular and intracellular potentials, respectively, and in $\Omega_{\rm B}$, $\varphi_s$ is thoracic medium electric potential. The parameter $\mathfrak{c}_{\alpha}$ is $\mathfrak{c}_{\alpha} = \kappa C_{\alpha} > 0$, where $C_{\alpha}$ is the membrane capacitance per unit area and $\kappa$ is the surface area-to-volume ratio (homogenization parameter). The membrane is assumed to be passive, so the capacitance $C_{\alpha}$ 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^{\alpha}_c = \mathfrak{c}_{\alpha}\partial_{0^{+}}^{\alpha} \phi$, in parallel, with $\alpha$ 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 $\partial_{0^{+}}^{\beta} u$.

^† which traditionally describes the recovery of APD as a function of the interbeat interval

The coupling between equations in $\Omega_{\rm H}$ and equation in $\Omega_{\rm B}$ (in systems (1.1) and (1.2)) is operate at the heart-torso interface $\Gamma_{\rm H}$. The continuity conditions at the interface can be given by

where ${\mathbf n}$ is the outward normal to $\Gamma_{\rm H}$.

The tensors $\mathcal{K}_i$ and $\mathcal{K}_e$ are the conductivity tensors describing anisotropic intracellular and extracellular conductive media, and tensor $\mathcal{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 ${\cal I}$ is equal to $\kappa {I}_{ion}$, where the nonlinear operator ${I}_{ion}$ describes the sum of transmembrane ionic currents across cell membrane with $u$. The nonlinear operator $\mathcal{G}$ is representing the ionic activity in myocardium. Functional forms for ${\cal I}$ and ${\cal G}$ are determined by an electrophysiological cell model. The source terms are $I_{i} = \kappa f_{i}$, $I_{e} = \kappa 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

where ${\mathbf n}_T$ is the outward normal to $\Gamma_T$.

In absence of a grounded electrode, the bidomain equations are a naturally singular problem since $\varphi_{e}$ and $\varphi_{s}$, in system (1.2), only appear in the equations and boundary conditions through their gradients. Moreover, the states $\varphi_{e}$ and $\varphi_{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)$)

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_{\rm 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 $\varphi_{e}$ and $\varphi_{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))

Remark 1.1. 1.Condition (1.6) is used for pressure in oceanography (see e.g. ^[15]).

2. The functions $\mathcal{K}_{i}$, $\mathcal{K}_{e}$, $\mathcal{H}$ and $\mathcal{G}$ depend on the fiber extension ratio.

3. We can suppose, for example, that $f_{\rm H}$ is only a time dependent source and is of the form

where $\chi_{\Omega^{(i)}_{\rm H}}$ is the characteristic function of set $\Omega^{(i)}_{\rm H}$, $i = 1, 2$. The support regions $\Omega^{(1)}_{\rm H}$ and $\Omega^{(2)}_{\rm 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,29,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.

2.   Assumptions, notations and some fundamental inequalities

Let $\mho\subset \mbox{IR}^{{m}}$, $m\geq 1$, be an open and bounded set with a smooth boundary and $\mho_T = \mho\times (0, T)$. We use the standard notation for Sobolev spaces (see ^[1]), denoting the norm of $W^{q, p}(\mho)$ ($q \in \mbox{IN}^{{}}$, $p\in [1, \infty]$) by $||.|| _{W^{q, p}}$. In the special case $p = 2$, we use $H^{q}(\mho)$ instead of $W^{q, 2}(\mho)$. The duality pairing of a Banach space $X$ with its dual space $X'$ is given by $\langle., .\rangle_{X', X}$. For a Hilbert space $Y$ the inner product is denoted by $(., .)_{Y}$ and the inner product in $L^2(\Omega)$ is denoted by $(., .)$. For any pair of real numbers $r, s\geq 0$, we introduce the Sobolev space $H^{r, s}(\mho_T)$ defined by $H^{r, s}(\mho_T)\! = \!L^{2}(0, T;H^{r}(\mho))\cap H^{s}(0, T;L^{2}(\mho)),$ which is a Hilbert space normed by $\left(|| v|| _{L^{2}(0, T;H^{r}(\mho))}^{2}+|| \!v|| _{H^{s}(0, T;L^{2}(\mho))}^{2}\right)^{1/2},$ where $H^{s}(0, T;L^{2}(\mho))$ denotes the Sobolev space of order $s$ of functions defined on $(0, T)$ and taking values in $L^{2}(\mho)$, and defined (see ^[52]) by $H^{s}(0, T;L^{2}(\mho))\! = \![H^{q}(0, T, L^{2}(\mho)), L^{2}(\mho_T)]_{\theta}$, for $s = (1-\theta)q$ with $\theta\in (0, 1)$ and $q \in \mbox{IN}^{{}}$, and $H^{q}(0, T;L^{2}(\mho))\! = \! \left\{v \in L^{2}(\mho_T) | \; \frac{\partial ^{j}v}{\partial t^{j}} \in L^{2}(\mho_T), \text{for } 1\leq j\leq q \right\}$. For a given Banach space $\mathbb{X}$, with norm $\Vert. \Vert _{\mathbb{X}}$, of functions integrable on $\mho$, we define its subspace $\mathbb{X}|_{\mbox{IR}^{{}}} = \Big\lbrace u\in\mathbb{X}, \int\limits_\mho u = 0 \Big\rbrace$ that is a Banach space with norm $\Vert. \Vert _{\mathbb{X}}$, and we denote by $[u]$ the projection of $u\in \mathbb{X}$ on $\mathbb{X}|_{\mbox{IR}^{{}}}$ such that $[u] = u-\frac{1}{mes(\mho)}\int\limits_\mho u d\mathbf{x}$ (with $mes(\mho)$ standing for Lebesgue measure of $\mho$).

Lemma 2.1. (Poincaré-Wirtinger inequality) Assume that $1 \leq p \leq \infty$ and that $\mho$ is a bounded connected open subset of $\mathit{\mbox{IR}}^{{d}}$ with a sufficiently regular boundary $\partial \Omega$ (e.g., a Lipschitz boundary). Then there exists a Poincaré constant $C$, depending only on $\Omega$ and $p$, such that for every function $u$ in Sobolev space $W^{1, p}(\mho)$, we have

Remark 2.1. From the Poincaré-Wirtinger inequality, we can deduce that the $H^{1}$ semi-norm and the $H^{1}$ norm are equivalent in the space $H^{1}(\mho)|_{\mbox{IR}^{{}}}$.

Remark 2.2. Let $q$ be a nonnegative integer. We have the following results (see e.g. ^[1])

(i) $H^{q}(\mho)\subset L^{p}(\mho)$, $\forall p\in[1, \frac{2m}{m-2q}]$, with continuous embedding (with the exception that if $2q = m$, then $p \in[1, +\infty[$ and if $2q > m$, then $p \in[1, +\infty]$).

(ii) (Gagliardo-Nirenberg inequalities) There exists $C > 0$ such that

where $0\leq \theta < 1$ and $p = \frac{2m}{m-2\theta q}$ (with the exception that if $q-m/2$ is a nonnegative integer, then $\theta$ is restricted to 0).

Finally, we introduce the spaces:

● $\mathbb{H}_{\rm H} = L^2(\Omega_{\rm H})$, $\mathbb{V}_{\rm H} = H^1(\Omega_{\rm H})$, $\mathbb{V} = H^1(\Omega)$ (endowed with their usual norms) and $\mathbb{U}_{\rm H} = \mathbb{V}_{\rm H}|_{\mbox{IR}^{}}$,

● $\mathbb{U}_{\rm HB} = \Big\lbrace \psi\in H^1(\Omega)\; | \; \int\limits_{\Omega_{\rm H}} \psi d\mathbf{x} = 0 \Big\rbrace$.

We will denote by $\mathbb{V}_{\rm H}'$ (resp. $\mathbb{U}_{\rm H}'$) the dual of $\mathbb{V}_{\rm H}$ (resp. of $\mathbb{U}_{\rm H}$). We have the following continuous embeddings, where $p\geq2$ if $d = 2$ and $2\leq p \leq 6$ if $d = 3$, $p^\prime$ is such that $\frac{1}{p^\prime}+\frac{1}{p} = 1$

and the injections $\mathbb{V}_{\rm H} \subset \mathbb{H}_{\rm H}$ and $\mathbb{U}_{\rm H} \subset \mathbb{H}_{\rm H}|_{\mbox{IR}^{{}}}$ are compact. We can now introduce the following spaces (where $q > 1$, $p\geq 2$ and $\frac{1}{p}+\frac{1}{p'} = 1$)

Remark 2.3. Space $\mathbb{D}_{p}(0, T)$ is equipped with norm: $|| u|| _{\mathbb{D}_{p}} = \max(|| u|| _{L^{p}({\cal Q}_{\rm H})}, || u || _{L^{2}(0, T; \mathbb{V}_{\rm H})})$ and its dual with norm: $|| v|| _{\mathbb{D}'_{p}}\! = \!\inf_{v = v_1+v_2}(|| v_1|| _{L^{p'}({\cal Q}_{\rm H})}+ || v_2 || _{L^{2}(0, T; \mathbb{V}_{\rm H}')})$.

Definition 2.1. A real valued function ${\cal H}$ defined on $D\times \mbox{IR}^{{q}}$, $q\geq 1$, is a Carathéodory function iff ${\cal H}(.; \mathbf{v})$ is measurable for all $\mathbf{v} \in\mbox{IR}^{{q}}$ and ${\cal H}(y; .)$ is continuous for almost all $y\in D$.

Our study involves the following fundamental inequalities, which are repeated here for review:

(i) Hölder's inequality: $\int_{D} \Pi_{i = 1, k}f_{i}dx\leq \Pi_{i = 1, k}|| f_{i}|| _{L^{q_{i}}(D)}$,

where $|| f_{i}|| _{L^{q_{i}}(D)} = \left(\int_{D} \mid f_{i}\mid^{q_{i}} dx\right)^{1/q_{i}} \; \; \text{and}\; \; \sum_{1\leq i \leq k}\frac{1}{q_{i}} = 1.$

(ii) Young's inequality ($\forall a, b > 0$ and $\epsilon > 0$): $ab\leq \frac{\epsilon}{p}a^{p}+\frac{\epsilon^{-q/p}}{q}b^{q}, \; for\; p, q\in]1, +\infty[\; and \; \frac{1}{p}+\frac{1}{q} = 1.$

(iii) Minkowski's integral inequality:

Finally, we denote by ${\mathfrak L}(A; B)$ the set of linear and continuous operators from a vectorial space $A$ into a vectorial space $B$, and by ${\cal R}^{*}$ the adjoint operator to a linear operator ${\cal 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.

(H1) We assume that the conductivity tensor functions $\mathcal{K}_{\theta}\in W^{1, \infty}(\overline{\Omega}_{\rm H})$, $\theta \in \{i, e\}$ and $\mathcal{K}_{s}\in W^{1, \infty}(\overline{\Omega}_{\rm B})$ are symmetric, positive definite matrix functions and that they are uniformly elliptic, i.e., there exist constants $0 < K_1 < K_2$ such that ($\forall \psi \in\mbox{IR}^{{d}}$)

Remark 2.4. We can emphasize a specificity of tensors $\mathcal{K}_{e}$ and $\mathcal{K}_{i}$ (see e.g., ^[25]).

1. The tensors $\mathcal{K}_{e}(\mathbf{x})$ and $\mathcal{K}_{i}(\mathbf{x})$ have the same basis of eigenvectors $Q(\mathbf{x}) = (q_k(\mathbf{x}))_{1\leq k\leq d}$ in $\mbox{IR}^{{d}}$, which reflect the organization of muscle in fibers, and consequently $\mathcal{K}_{i}(\mathbf{x}) = Q(\mathbf{x})\Lambda_{i}(\mathbf{x})Q(\mathbf{x})^T$ and $\mathcal{K}_{e}(\mathbf{x}) = Q(\mathbf{x})\Lambda_{e}(\mathbf{x})Q(\mathbf{x})^T$, where $\Lambda_{i}(\mathbf{x}) = diag((\lambda_{i, k})_{1\leq k\leq d})$ and $\Lambda_{e}(\mathbf{x}) = diag((\lambda_{e, k})_{1\leq k\leq d})$.

2. The muscle fibers are tangent to $\Gamma$ so that (for $\theta \in \{i, e\}$) : $\mathcal{K}_{\theta}{\mathbf n} = \lambda_{\theta, d}{\mathbf n}, \; \; \text{a.e., in}\; \; \Gamma,$ with $\lambda_{\theta, d}(\mathbf{x})\geq \lambda > 0$, $\lambda$ a constant.

The operators $\mathcal{I}$ and $\mathcal{G}$ which describe electrophysiological behavior of the system can be taken as follows (affine functions with respect to $u$)

where $\hbar$ is a sufficiently regular function and $\mathcal{I}_{0}(\mathbf{x}, t;\phi) = g_1(\mathbf{x}; \phi)+g_2(\mathbf{x}, t;\phi)$ with $g_1$ is an increasing function on $\phi$. Moreover, the operators $\mathcal{I}_{0}$, $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ appearing in $\mathcal{I}$ and $\mathcal{G}$, are supposed to satisfy the following assumptions.

${\rm {\mathbf (H2)}_{p}}$ The operators $\mathcal{I}_{0}$, $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ are Carathéodory functions from $(\Omega \times \mbox{IR}^{{}}) \times \mbox{IR}^{{}}$ into $\mbox{IR}^{{}}$ and continuous on $\phi$ (as in ^[12]). Furthermore, for some $p \geq 2$ if $d = 2$ and $p \in [2,6]$ if $d = 3$, the following requirements hold

(i) there exist constants $\beta_i \geq 0\; (i = 1, \dots, 10)$ such that for any $v\in \mbox{IR}^{}$

where $E_g$ is the primitive of $g_1$.

(ii) there exist constants $\mu_1 > 0, \mu_2 > 0, \mu_i\geq 0$ (for $i = 3, 7$) such that for any $(v, w)\in \mbox{IR}^{{2}}$:

In order to assure the uniqueness of solution we assume that

(H3) The Nemytskii operators $\mathcal{I}$ and $\mathcal{G}$ satisfy Carathéodory conditions and there exists some $\mu > 0$ such the operator $F_\mu:\mbox{IR}^{{2}}\rightarrow \mbox{IR}^{{2}}$ defined by

satisfies a one-sided Lipschitz condition (see e.g., ^[14,70]): there exists a constant $C_{L} > 0$ such that ($\forall \mathbf{v}_i = (v_{i}, w_{i})\in \mbox{IR}^{{2}}, \; i = 1, 2$)

Lemma 2.2. (^[9]) Assume that $F_{\mu}$ is differentiable with respect to $(\phi, u)$ and denote by $\lambda_1(\phi, u) \leq \lambda_2(\phi, u)$ the eigenvalues of symmetrical part of Jacobian matrix $\nabla F_\mu(\phi, u)$:

If there exists a constant $C_F$ independent of $\phi$ and $u$ such as:

then $F_{\mu}$ satisfies the hypothesis (H3).

Lemma 2.3. (^[9]) Let assumptions (2.3), (H1) and ${\rm {\mathbf (H2)}_{p}}$ be fulfilled. For $(\phi, u)\in L^p(\Omega) \times \mathbb{H}$ and $a.e.$, $t$, there exist constants $C_{i} > 0$ ($i = 1, 6$) such that

where $p'$ is such that $\frac{1}{p}+\frac{1}{p'} = 1$.

In this work we assume that $p = 4$. The considered functions $\mathcal{I}_{i}$, in this paper, include the three classical type models in which assumptions ${\mathbf(H1)}, {\rm {\mathbf (H2)}_{4}}$ and (H3) are satisfied namely the Rogers-McCulloch ^[57] (RM), Fitz-Hugh-Nagumo ^[39] (FHN) and Aliev-Panfilov ^[61](LAP) models as follows. The function $\mathcal{I}_{0}$ is defined by a cubic reaction term of the form $\mathcal{I}_{0}(.; v) = b_{1}(.)v(v-r)(v-1),$ and the functions $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ are given by

where $b_{i}\in W^{1, \infty}({\cal Q})$, $i = 1, 3$, are sufficiently regular functions from ${\cal Q}$ into $\mbox{IR}^{{+, *}}$ and $r\in[0, 1]$. We obtain easily the following Lemma.

Lemma 2.4. The following properties hold. For all $v_{1}$, $v_{2}$ in $\mathit{\mbox{IR}}^{{}}$ we have

For the sake of simplicity, we shall write $\mathcal{I}_{i}(\psi)$, ${\cal I}(\psi, v)$ and ${\cal G}(\psi, v)$ in place of $\mathcal{I}_{i}(\mathbf{x}, t;\psi)$, ${\cal I}(\mathbf{x}, t;\psi, v)$ and ${\cal G}(\mathbf{x}, t;\psi, v)$, respectively (for $i = 0, 2$).

3.   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.

3.1.   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 $\gamma\in]0, 1]$ and $X$ be a Banach space, we start from a formal level and assume the given functions $f: t\in (a, T) \to f(t)\in X$ and $g: t\in (a, T) \to f(t)\in X$ are sufficiently smooth (with $-\infty < a < T < \infty$).

Definition 3.1. The forward and backward Riemann-Liouville fractional integrals of fractional order $\gamma$ on $(a, T)$ are defined, respectively, by ($t\in (a, T)$)

where ${\mathbf \Gamma}(z) = \int_{0}^{\infty} e^{\tau}\tau^{z-1} d\tau$ is the Euler ${\mathbf \Gamma}$-function.

The basic equality for the fractional integral is (from Fubini's Theorem and the relationship between ${\mathbf \Gamma}$-function and ${\mathbf \beta}$-function)

and holds for a $L^{p}$-function $f$ ($1\leq p\leq \infty$).

Definition 3.2. The forward Riemann-Liouville and Caputo derivatives of fractional order $\gamma$ on $(a, T)$ are defined, respectively, by ($t\in (a, T)$)

From (3.2) we can deduce the following relation between fractional integral and Caputo derivative

Definition 3.3. The backward Riemann-Liouville and Caputo derivatives of fractional order $\gamma$, on $(a, T)$ are defined, respectively, by ($t\in (a, T)$)

By substitution $a\to -\infty$, 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\in (-\infty, T)$)

Remark 3.1. 1. For $\gamma \to 1$ the forward (resp. backward) Riemann-Liouville and Caputo derivatives of fractional order $\gamma$ of $f$ converge to the classical derivative $\frac{df}{dt}$ (resp. to $-\frac{df}{dt}$). Moreover, Riemann-Liouville fractional derivative of fractional order $\gamma$ of constant function $t \to f(t) = k$ is not 0 since $D^{\gamma}_{a^{+}}f(t) = \frac{k}{{\mathbf \Gamma}(1-\gamma)}\frac{d}{dt}(\int_{a}^{t}(t-\tau)^{-\gamma}d\tau) = \frac{k(t-a)^{-\gamma}}{{\mathbf \Gamma}(1-\gamma)}.$

2. 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 $f\in C^1([a, T], X)$ we have

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^{\gamma}_{a^{+}}$ is a continuous operator from

(i) $L^p(a, T)$ into $L^p(a, T)$, for any $p \geq 1$,

(ii) $L^p(a, T)$ into $L^r(0, T)$, for any $p\in (1, 1/\gamma)$ and $r\in [1, p/(1-\gamma p)]$,

(iii) $L^p(a, T)$ into $C^{0, \gamma-1/p}([a, T])$, for any $p\in]1/\gamma, +\infty[$,

(iii) $L^{1/\gamma}(a, T)$ into $L^p(a, T)$, for any $p\in [1, +\infty)$

(iv) $L^{\infty}(a, T)$ into $C^{0, \gamma}([a, T])$.

From Lemma 3.1 and (3.4) we can deduce the following corollary.

Lemma 3.2. Let $X$ be a Banach space and $\gamma\in]0, 1]$. Suppose the Caputo derivative $\partial^{\gamma}_{a^{+}}f\in L^{p}(a, T; X)$ and $p > \frac{1}{\gamma}$, then $f\in C^{0, \gamma-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.3. Let $0 < \gamma\leq 1$ and $p, q\geq 1$ with $1/p+1/q\leq 1+\gamma$. Then

(i) if $g$ is a $L^p$-function on $(a, T)$ and $g$ is a $L^q$-function on $(a, T)$, then

(ii) if $f\in I^{\gamma}_{T^{-}}(L^{p})$ and $g\in I^{\gamma}_{a^{+}}(L^{q})$, then

Lemma 3.4. Let $0 < \gamma\leq 1$, and $g$ be $L^p$-function on $(a, T)$ (for $p\geq 1$) and $f$ be absolutely continuous function on $[a, T]$. Then {

(i) $\int_{a}^{T} (\partial^{\gamma}_{a^{+}}f(\tau), g(\tau))_{X}d\tau = -\int_{a}^{T} (D^{\gamma}_{T^{-}}g(\tau), f(\tau))_{X}d\tau+ |(I^{1-\gamma}_{T^{-}}\left[g\right](\tau), f(\tau))_{X}|_{a}^{T},$

(ii) $\int_{a}^{T} (D^{\gamma}_{a^{+}}f(\tau), g(\tau))_{X}d\tau = -\int_{a}^{T} (D^{\gamma}_{T^{-}}g(\tau), f(\tau))_{X}d\tau+ (I^{1-\gamma}_{T^{-}}\left[g\right](T^{-}), f(T))_{X},$

(iii) $\int_{a}^{T} (D^{\gamma}_{T^{-}}f(\tau), g(\tau))_{X}d\tau = -\int_{a}^{T} (D^{\gamma}_{a^{+}}g(\tau), f(\tau))_{X}d\tau- (I^{1-\gamma}_{a^{+}}\left[g\right](a^{+}), f(a))_{X}.$

From ^[85], we can deduce the following Lemma.

Lemma 3.5. (A generalized Gronwall's inequality)

Assume $\gamma > 0$, $h$ is a nonnegative function locally integrable on $(0, T)$ and $b$ is a nonnegative, bounded, nondecreasing continuous function defined on $[0, T)$. Let $f$ be a nonnegative and locally integrable function on $(0, T)$ with

Then (for $t\in(0, T)$)

If in addition $h$ is a nondecreasing function on $(0, T)$, then

The used function $E_{\theta_1, \theta_2}$ is the classical two-parametric Mittag-Leffler function (usually denoted by $E_{\theta_1}$ if $\theta_2 = 1$) which is defined by

The function $E_{\theta_1, \theta_2}$ is an entire function of the variable $z$ for any $\theta_1, \theta_2 \in \mbox{lC}^{{}}$, $Re(\theta_1) > 0$.

Next we declare a compactness theorem in Hilbert spaces. Assume that $\mathcal{X}_{0}$, $\mathcal{X}_{1}$ and $\mathcal{X}$ are Hilbert spaces with

The Fourier transform of $f: \mbox{IR}^{{}} \to \mathcal{X}_{1}$ is defined by $\widehat{f}(\tau) = \int_{-\infty}^{\infty}\exp(-2i\pi s\tau)f(s)ds$ and we have Lemma (see e.g., ^[68])

Lemma 3.6. Let $\gamma\in (0, 1)$ and $f$ be a $L^1$-function in $\mathit{\mbox{IR}}^{{}}$ with compact support. Then

(i) $\widehat{I^{\gamma}_{-\infty}{f}}(\tau) = (2i\pi \tau)^{-\gamma} \widehat{f}(\tau),$

(ii) $\widehat{\partial^{\gamma}_{-\infty}{f}}(\tau) = (2i\pi \tau)^{\gamma} \widehat{f}(\tau).$

We define the Hilbert space $\mathcal{W}^{\gamma}(\mbox{IR}^{{}}; \mathcal{X}_{0}, \mathcal{X}_{1})$, for a given $\gamma > 0$, by

endowed with the norm

For any subset $K$ of $\mbox{IR}^{{}}$, we define the subspace $\mathbb{W}_{K}^{\gamma}$ of $\mathbb{W}^{\gamma}$ 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 $\mathcal{X}_{0}$, $\mathcal{X}_{1}$ and $\mathcal{X}$ be Hilbert spaces with the injection (3.9). Then for any bounded set $K$ and any $\gamma > 0$, the injection of $\mathbb{W}_{K}^{\gamma}(\mathit{\mbox{IR}}^{{}}; \mathcal{X}_{0}, \mathcal{X}_{1})$ into $L^{2}(\mathit{\mbox{IR}}^{{}}; \mathcal{X})$ is compact.

3.2.   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.

Theorem 3.2. Assume that $\gamma_1$ and $\gamma_2$ be in]0, 1] with $\gamma_2\leq \gamma_1$ and $\Psi$, $h_{i}$, for $i = 1, 2$, are nonnegative functions locally integrable on $(0, T)$. If $f$ and $g$ are nonnegative functions locally integrable, with $f(0) = f_0$, $g(0) = g_0$, and satisfy the inequalities (for $t\in (0, T)$)

where $c_i > 0$, $i = 1, 2$, $\eta\geq 0$ and $\xi\geq 0$ are four constants with $\eta\leq \frac{\xi{\mathbf \Gamma}(\gamma_1)}{{\mathbf \Gamma}(\gamma_2)T^{(\gamma_1-\gamma_2)}}$. Then we have the following estimate

where $F = f+g$, $F_0 = f_0+g_0$, $d_{0} = \frac{T^{(\gamma_1-\gamma_2)}{\mathbf \Gamma}(\gamma_2)}{{\mathbf \Gamma}(\gamma_1)}$ and $d_1 = c_1 d_0+c_2.$

The function $E_{., .}$ is the Mittag-Leffler function which is defined by (3.8).

Proof. From (3.4) we can deduce that

Put $F = f+g$, then

We can deduce that (since $t\in (0, T)$ and $\gamma_1-\gamma_2\geq 0$)

Since $\eta{\mathbf \Gamma}(\gamma_2)T^{(\gamma_1-\gamma_2)}-\xi{\mathbf \Gamma}(\gamma_1)\leq 0$, then

Thus

where $d_{0} = \frac{T^{(\gamma_1-\gamma_2)}{\mathbf \Gamma}(\gamma_2)}{{\mathbf \Gamma}(\gamma_1)}$ and $d_1 = \frac{c_1{\mathbf \Gamma}(\gamma_2)T^{(\gamma_1-\gamma_2)}+c_2{\mathbf \Gamma}(\gamma_1)}{{\mathbf \Gamma}(\gamma_1)} = c_1 d_0+c_2$.

Finally, from Lemma 3.5 and relation (3.2), we get (since $I^{k\gamma_2}_{0^{+}}\left[F\right](0) = F_0\; I^{k\gamma_2}_{0^{+}}\left[.1\right] = \frac{F_0 t^{k\gamma_2}}{{\mathbf \Gamma}(1+k\gamma_2)}$)

This completes the proof.

Corollary 3.1. Assume that assumptions of Theorem 3.2 hold, and that $h_i$, for $i = 1, 2$, are nondecreasing functions on $(0, T)$. Then, if $f$ and $g$ are nonnegative functions locally integrable, with $f(0) = f_0$, $g(0) = g_0$, and satisfy the inequalities (3.10), we have the following estimate

where $F = f+g$, $F_0 = f_0+g_0$, $d_{0} = \frac{T^{(\gamma_1-\gamma_2)}{\mathbf \Gamma}(\gamma_2)}{{\mathbf \Gamma}(\gamma_1)}$ and $d_1 = \frac{c_1{\mathbf \Gamma}(\gamma_2)T^{(\gamma_1-\gamma_2)}+c_2{\mathbf \Gamma}(\gamma_1)}{{\mathbf \Gamma}(\gamma_1)} = c_1 d_0+c_2$.

Proof. Since $h_i$, for $i = 1, 2$, are nondecreasing functions then

Consequently, from (3.11) we can deduce the result. This completes the proof.

Corollary 3.2. Assume that assumptions of Theorem 3.2 hold, and that $h_i$ are $L^{\frac{1}{\alpha_i}}$ nondecreasing functions on $(0, T)$, with $\alpha_{i}\in [0, \gamma_i[$, for $i = 1, 2$. Then, if $f$ and $g$ are nonnegative functions locally integrable, and satisfy the inequalities (3.10), we have the following estimate

where $r_i = (\frac{\gamma_i-\alpha_i}{1-\alpha_i})^{{1-\alpha_i}}$, for $i = 1, 2$, $F = f+g$, $F_0 = f_0+g_0$, ${d_0} = \frac{{{T^{({\gamma _1} - {\gamma _2})\;\;}}{\bf{\Gamma }}({\gamma _2})}}{{{\bf{\Gamma }}({\gamma _1})}}$ and ${d_1} = \frac{{{c_1}{\bf{\Gamma }}({\gamma _2}){T^{({\gamma _1} - {\gamma _2})}}\;\; + {c_2}\;{\bf{\Gamma }}({\gamma _1})}}{{{\bf{\Gamma }}({\gamma _1})}} = {c_1}{d_0} + {c_2}$.

Proof. From Hölder's inequality, we obtain (for $i = 1, 2$)

where $r_i = (\frac{\gamma_i-\alpha_i}{1-\alpha_i})^{{1-\alpha_i}}$.

According to (3.11), we have

Thus

This completes the proof.

In the sequel we will always denote $C$ some positive constant which may be different at each occurrence.

4.   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.

4.1.   Variational formulation and preliminary results

First we define the state $\psi$, as the extracellular cardiac potential in $\Omega_{\rm H}$, and the thoracic potential in $\Omega_{\rm B}$, by

Then the potential $\phi$ can be written as $\phi = \varphi_i-\psi|_{\Omega_{\rm H}}$. From the interface coupling condition (1.3), it follows that if $\varphi_e\in H^{1}(\Omega_{\rm H})$ and $\varphi_s\in H^{1}(\Omega_{\rm B})$ then $\psi\in H^{1}(\Omega)$. Similarly, we introduce the corresponding conductivity tensor $\mathcal{K}_g\in W^{1, \infty}(\overline{\Omega})$ as follows :

We now define the following forms

Proposition 4.1. (i) The forms $A_k$ (for $k = i, e, s$) and ${\mathbf A_g}$, are symmetric bilinear continuous forms on $H^{1}$.

(ii) The forms $A_k$, (for $k = i, e, s$) and ${\mathbf A_g}$, are coercive on $H^1$ (we denote by $\nu_k$, for $k = i, e, s, g$, their coercivity coefficients).

Proof. (i) and (ii) are easily obtained providing that properties of tensors $\mathcal{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).

Definition 4.1. A quadruplet of state functions $(\varphi_i, \phi, \psi, u)$ such that : $\varphi_i\in L^2(0, T; \mathbb{V}_{\rm H})$, $\phi\in \mathbb{D}_{4}(0, T) \cap L^{\infty}(0, T;L^2(\Omega_{\rm H}))$, $\partial^{\alpha}_t \phi\in \mathbb{D}'_{4}(0, T)$, $\psi\in L^2(0, T; \mathbb{V})$, $u\in L^{\infty}(0, T; L^3(\Omega_{\rm H}))$ and $\partial^{\beta}_t u\in L^2(0, T; L^2(\Omega_{\rm H}))$ is said to be a weak solution to system (1.1) if it satisfies the following weak formulation (since $I_i = -I_e = f_{\rm H}$)

for all $v_i\in\mathbb{V}_{\rm H}, \; v\in\mathbb{V}$, $\rho\in\mathbb{V}_{\rm H}$.

The results of next section concern the existence, uniqueness and regularity of solution of problem (1.1).

4.2.   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).

4.2.1.   Non degenerated problem study

Let $\epsilon > 0$ be a small positive number. Hence, our approximation system read as (since $I_i = -I_e = f_{\rm H}$)

The auxiliary initial conditions for $\varphi^\epsilon_i$ and $\varphi^\epsilon_s$, needed by problem (4.4), are defined by introducing two arbitrary functions $\varphi_{e, 0}$ and $\varphi_{s, 0}$. The weak formulation of (4.4) is given by (for all $v_i\in\mathbb{V}_{\rm H}, \; v\in\mathbb{V}$, $\rho\in\mathbb{V}_{\rm H}$ and $a.e$. $t\in (0, T)$)

where $\psi^\epsilon$ and the initial condition $\psi_{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)).

Theorem 4.1. Let assumptions (H1) and ${\rm {\mathbf (H2)}_{4}}$ be fulfilled, $\beta\geq \alpha$ and $0 < \epsilon\leq 1$ be given. Then for $(\varphi_{i, 0}, \phi_{0}, \psi_{0}, u_0)$ and $f_{\rm H}$ given such that $(\varphi_{i, 0}, \phi_{0}, \psi_{0}, u_0)\in \mathbb{V}_{\rm H}^2\times \mathbb{V}\times L^{3}(\Omega_{\rm H})$ and $f_{\rm H}\in L^{\frac{2}{\alpha_1}}(0, T;L^{2}(\Omega_{\rm H}))$ with $0 < \alpha_1 < \alpha$, there exists a solution $(\varphi^{\epsilon}_{i}, \phi^{\epsilon}, \psi^{\epsilon}, u^{\epsilon})$ of problem (4.5) verifying

Moreover, if $\alpha > 1/2$ and $\beta > 1/2$ we have the continuity of the solution at $t = 0$.

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\geq 1}$ be a Hilbert basis and orthogonal in $L^{2}$ of $\mathbb{V}_{\rm H}$, $(f_{k})_{k\geq 1}$ be a Hilbert basis of $\mathbb{U}_{\rm H}$, $(g_{k})_{k\geq 1}$ be a Hilbert basis of $\mathbb{W}_{\rm B} = \{v\in H^{1}(\Omega_{\rm B})\; |\; v|_{\Gamma_{T}} = 0\}$. We denote by $\tilde{f}_{k}$ an extension of $f_{k}$ in $H^{1}(\Omega)$, $\tilde{g}_{k}$ an extension of $g_{k}$ in $H^{1}(\Omega)$ with $\tilde{g}_{k}|_{\Omega_{\rm H}} = 0$ and by $(e_{k})_{k\geq 1}$ a Galerkin basis of $\mathbb{U}_{\rm HB}$ which is defined as $e_{2k} = \tilde{g}_{k}$ and $e_{2k-1} = \tilde{f}_{k}$ (for all $k > 1$).

For all $n\in \mbox{IN}^{{*}}$, we denote by $W_{i, n} = span(w_1, \cdots, w_n)$, $W_{e, n} = span(f_1, \cdots, f_n)$, $W_{s, n} = span(g_1, \cdots, g_n)$ and $W_{g, n} = span(e_1, \cdots, e_{2n})$ the spaces generated, respectively, by $(w_{k})_{n\geq k\geq 1}$, $(f_{k})_{n\geq k\geq 1}$, $(g_{k})_{n\geq k\geq 1}$ and $(e_{k})_{2n\geq k\geq 1}$, and we introduce the orthogonal projector $L_{j, n}$ on the spaces $W_{j, n}$ (for $j = i, e, s, g$).

For each $n$, we would like to define the approximate solution $(\phi^\epsilon_{n}, \varphi^\epsilon_{i, n}, \psi^\epsilon_{n}, u^\epsilon_{n})$ of the problem (4.4). Setting

where $\boldsymbol{\omega}^\epsilon_i = (\varpi^\epsilon_{i, n, l})_{l = 1, n}$, $\boldsymbol{\omega}^\epsilon = (\varpi^\epsilon_{n, l})_{l = 1, n}$, $\vartheta^\epsilon = (\upsilon^\epsilon_{n, l})_{l = 1, n}$ and $\Pi^\epsilon = (\pi^\epsilon_{n, m})_{m = 1, 2n}$ are unknown functions, and replacing $(\phi^\epsilon, \varphi^\epsilon_{i}, \psi^\epsilon, u^\epsilon)$ by $(\phi^\epsilon_{n}, \varphi^\epsilon_{i, n}, \psi^\epsilon_{n}, u^\epsilon_{n})$ in (4.4), we obtain $\forall l = 1, n$, $m = 1, 2n$ and $a.e$. $t\in (0, T)$, the system of Galerkin equations (since in $\Omega_{\rm H}$, $\phi^\epsilon_{n} = \varphi^\epsilon_{i, n}-\psi^\epsilon_{n}$)

where $(L_{i, n}\phi_{0}, L_{i, n}\varphi_{i, 0}, L_{g, n}\psi_{0}, L_{i, n}u_{0})$ satisfies (by construction)

The function $F^{\epsilon, i}_{n}$, $G^{\epsilon}_{n}$ and $H^{\epsilon}_{n}$ are defined by (for $(v, w, \rho)\in \mathbb{V}_{\rm H}\times \mathbb{V}\times \mathbb{V}_{\rm H}$)

Step 1. We show first that, for every $n$, the system (4.8) admits a unique local solution. The system (4.8) is equivalent to an initial-value for a system of nonlinear fractional differential equations for functions $(\boldsymbol{\omega}^\epsilon_i, \Pi^\epsilon, \vartheta^\epsilon)$

where ${\mathcal F}_{n} :[0, T]\times \mbox{IR}^{{4n}}\to \mbox{IR}^{{4n}}$ is a known nonlinear vector function and is given by

Thanks to assumptions ${\rm {\mathbf (H2)}_{4}}$ the function ${\mathcal F}_{n}$ is Carathéodory function.

Since, by construction the matrix $\mathcal{M} \in \mbox{IR}^{{4n\times 4n}}$ is symmetric and invertible, then by applying $\mathcal{M}_{n}^{-1}$ to the equation of problem (4.11), we obviously get

where ${\mathcal S}_{n} = \mathcal{M}_{n}^{-1}\mathcal{F}_{n}$.

The existence of a local absolutely continuous function $(\omega^\epsilon_i, \Pi^\epsilon, \vartheta^\epsilon)$ on an interval $[0, \tilde{T}]$, with $\tilde{T}\in]0, T]$ is insured by the standard FODE theory see ^[33,32] (see also e.g. ^[48]). Thus, we have a local solution $(\phi^\epsilon, \varphi_i^\epsilon, \psi^\epsilon, u^\epsilon)$ of problem (4.8) on $[0, \tilde{T}]$.

N.B. The constants which will be introduced in the sequel, are independent of $n$ and $\epsilon$ unless otherwise specified.

Step 2. We next derive a priori estimates for functions $\phi^\epsilon, \varphi_i^\epsilon, \psi^\epsilon, u^\epsilon$, which entail that $\tilde{T} = T$, by applying iteratively the step 1. For simplicity, in next step we omit the " $\tilde{}\;$" on $T$. Now, we set

where $\vartheta_{n, k}$, $\pi_{n, p}$ and $\zeta_{n, l}$ are absolutely continuous coefficients. Then, from (4.8), the approximation solution satisfies the following weak formulation

Taking ${\mathbf h}_{n} = \varphi_{i, n}^\epsilon$, ${\mathbf e}_{n} = \psi^{\epsilon}_{n}$ and ${\mathbf v}_{n} = u^\epsilon_{n}$ and using the uniform coercivity of the forms $A_i$ and ${\mathbf A_g}$, we can deduce (since $\psi^\epsilon_{n} = \varphi^\epsilon_{e, n}$ and $\phi^\epsilon_{n} = \varphi^\epsilon_{i, n}-\varphi^\epsilon_{e, n}$ in $\Omega_{\rm H}$)

Then the following inequality

Since, from ${\rm {\mathbf (H2)}_{4}}$, $\mu_1 \mathcal{I}(\phi^\epsilon_{n}, u^\epsilon_{n})\phi^\epsilon_{n}+\mathcal{G}(\phi^\epsilon_{n}, u^\epsilon_{n})u^\epsilon_{n}\geq \mu_2\mid \phi^\epsilon_{n}\mid^{4}-\mu_3(\mu_1\mid \phi^\epsilon_{n}\mid^{2}+\mid u^\epsilon_{n}\mid^{2})-\mu_4$ and $\mid \mathcal{G}(\phi^\epsilon_{n}, u^\epsilon_{n})u^\epsilon_{n}\mid \leq c_1(1+\mid \phi^\epsilon_{n}\mid^{2}\mid u^\epsilon_{n}\mid+\mid u^\epsilon_{n}\mid^2)$, we obtain (from Hölder's inequality)

with the constants $\eta$ chosen appropriately. In particular we have

where $\Psi^\epsilon_{n} = \mathfrak{c}_{\alpha}|| \phi^\epsilon_{n}|| ^{2}_{L^{2}(\Omega_{\rm H})}+(|| \sqrt{\epsilon}\psi^\epsilon_{n}|| ^{2}_{L^{2}(\Omega)}+|| \sqrt{\epsilon}\varphi^\epsilon_{i, n}|| ^{2}_{L^{2}(\Omega_{\rm H})})$ and $U^\epsilon_{n} = || u^\epsilon_{n}|| ^{2}_{L^{2}(\Omega_{\rm H})}$.

According to regularity $f_{\rm H}\in L^{2/\alpha_1}(0, T, L^{2}(\Omega_{\rm H}))$ and by choosing $\eta = \frac{\Gamma(\beta)\mu_2/\mu_1}{T^{\beta-\alpha}\Gamma(\alpha)}$), we can get by Corollary 2 and the fact that quantity $\Psi^\epsilon_{n}(0)+U^\epsilon_{n}(0)$ is uniformly bounded (from (4.9) and $0 < \epsilon < 1$)

for all $t\in (0, T)$ and then $\Psi^\epsilon_{n}(t)$ and $U^\epsilon_{n}(t)$ is uniformly bounded with respect to $\epsilon$ and to $n$. This ensures that

Moreover, from the first inequality of (4.16), relations (3.4) and (4.9), we have that for $t\in (0, T)$ (since $\int_{0}^{t}(t-\tau)^{\alpha-1}d\tau = \frac{1}{\alpha}t^{\alpha}\leq \frac{1}{\alpha}T^{\alpha}$)

where $C_0 = (1+|| f_{\rm H}|| ^{2}_{L^{2/\alpha_1}(0, T, L^{2}(\Omega_{\rm H}))})$.

We can deduce that (since $\phi^\epsilon_{n} = \varphi^\epsilon_{i, n}-\varphi^\epsilon_{e, n}$ and $\psi^\epsilon_{n} = \varphi^\epsilon_{e, n}$ in $\Omega_{\rm H}$)

Hence (since $\alpha-1 \leq 0$ and $(t-\tau)\in (0, T)$, $\forall \tau \in (0, t)$ then $(t-\tau)^{\alpha-1}\geq T^{\alpha-1}$)

Using Lemma 3 and results (4.19) and (4.22), we get

Now we estimate the fractional derivative $\partial^{\alpha}_{0^+}(\phi^\epsilon_{n}, \varphi^\epsilon_{i, n}, \psi^\epsilon_{n})$ and $\partial^{\beta}_{0^+}u^\epsilon_{n}$. Taking ${\mathbf h}_{n} = \partial^{\alpha}_{0^+}\varphi_{i, n}^\epsilon$, ${\mathbf e}_{n} = \partial^{\alpha}_{0^+}\psi^{\epsilon}_{n}$ and ${\mathbf v}_{n} = \partial^{\beta}_{0^+}u^\epsilon_{n}$ and using uniform coercivity of forms $A_i$ and ${\mathbf A_g}$, we can deduce (since $\psi^\epsilon_{n} = \varphi^\epsilon_{e, n}$ and $\phi^\epsilon_{n} = \varphi^\epsilon_{i, n}-\varphi^\epsilon_{e, n}$ in $\Omega_{\rm H}$)

Then

Let us now estimate right-hand side of equations of system (4.25).

According to ${\rm {\mathbf (H2)}_{4}}$ and regularity of $\hbar$, we can deduce that

and then from second equation of (4.25) (according to (4.19) and (4.21))

Before to estimate the term $\int_{\Omega_{\rm H}} I_{0^{+}}^{\alpha}\left[\mathcal{I}(\phi^\epsilon_{n}, u^\epsilon_{n})\partial^{\alpha}_{0^+}\phi^\epsilon_{n}\right](t)d\mathbf{x}$, we prove that $u^\epsilon_{n}$ is uniformly bounded in $L^{\infty}(0, T, L^{3}(\Omega_{\rm H}))$. Since $u^\epsilon_{n}$ satisfies

then, according to regularity of $\hbar$, we have

Consequently (since from assumption (2.4) we have $\mid \mathcal{I}_{2}(\phi^\epsilon_{n}) \mid \leq \beta_5 + \beta_6\mid \phi^\epsilon_{n} \mid^{2}$)

This implies (since $\int_{0}^{t}(t-\tau)^{\beta-1}d\tau\leq T^{\beta}/\beta$)

and then (using Minkowski inequality and the continuous embedding of $H^{1}$ in $L^{6}$)

Using (4.19), (4.21) and the fact, from (4.9) the quantities $u^\epsilon_{n}(0)$ is uniformly bounded in $L^3$, we obtain

and then (from Lemma 3.5)

Consequently

We can now estimate the term $\int_{\Omega_{\rm H}} I_{0^{+}}^{\alpha}\left[\mathcal{I}(\phi^\epsilon_{n}, u^\epsilon_{n})\partial^{\alpha}_{0^+}\phi^\epsilon_{n}\right](t)d\mathbf{x}$.

Since $\mathcal{I}(\phi^\epsilon_{n}, u^\epsilon_{n}) = \mathcal{I}_0(\phi^\epsilon_{n})+\mathcal{I}_1(\phi^\epsilon_{n}) u^\epsilon_{n}$ and $\mathcal{I}_0(\phi^\epsilon_{n}) = g_1(\phi^\epsilon_{n})+g_2(\phi^\epsilon_{n})$ with $g_1$ is an increasing function, we have first (from ${\rm {\mathbf (H2)}_{4}}$)

Moreover, since $g_1$ is an increasing function then the primitive $E_g$ of $g_1$ is convex. Hence, from ^[50], we can deduce (since $g_1$ is and independent on time) $\partial^{\alpha}_{0^+} E_g(\phi^\epsilon_{n})\leq g_1(\phi^\epsilon_{n})\partial^{\alpha}_{0^+}\phi^\epsilon_{n}$ and then

Now, we estimate the term $-I_{0^{+}}^{\alpha}\left[\partial^{\alpha}_{0^+} E(\phi^\epsilon_{n})\right]$.

Since, from ${\rm {\mathbf (H2)}_{4}}$, $\mu_6 \mid \phi^\epsilon_{n}\mid^2 -\mu_7\leq E_g(\phi^\epsilon_{n}) \leq c_2(1+\mid \phi^\epsilon_{n}\mid^4)$ then

According to (4.31), (4.32) and (4.34), we can deduce that

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 $\epsilon\in (0, 1)$, the quantities $\psi^\epsilon_{n}(0)$, $\phi^\epsilon_{n}(0)$ and $\varphi^\epsilon_{i, n}(0)$ are uniformly bounded in $H^1$, we can derive from (4.25) the following estimate (for all $t\in (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 $\partial^{\gamma}_{T_{k}^{+}}$ (for $\gamma = \alpha$ or $\beta$) 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.

Step 3. We are now ready to prove existence of solutions to system (4.5). From results (4.19), (4.22), (4.29) and (4.38), Theorem 3.1 and compactness argumentand, it follows that there exist $(X^{\epsilon}; u^{\epsilon}) = (\phi^\epsilon, \varphi_i^\epsilon, \psi^\epsilon; u^\epsilon)$ and $(\tilde{X}^{\epsilon}; \tilde{u}^\epsilon) = (\tilde{\phi}^\epsilon, \tilde{\varphi}_i^\epsilon, \tilde{\psi}^\epsilon; \tilde{u}^\epsilon)$ such that there exists a subsequence of $(X^{\epsilon}_{n}; u^{\epsilon}_{n}) = (\phi^\epsilon_{n}, \varphi_{i, n}^\epsilon, \psi^\epsilon_{n}; u^\epsilon_{n})$ also denoted by $(X^\epsilon_{n}; u_{n}^\epsilon)$, such that

First we show that $(\partial^{\alpha}_{0^+}X^\epsilon_{n}; \partial^{\beta}_{0^+} u_n^\epsilon)$ exists in the weak sense and that $(\partial^{\alpha}_{0^+}X^{\epsilon}; \partial^{\beta}_{0^+}u^\epsilon) = (\tilde{X}^\epsilon; \tilde{u}^{\epsilon}).$ Indeed, we take $\omega\in C_{0}^{\infty}(0, T)$ and $v\in \mathbb{V}_{\rm H}$ (then $\omega v\in \mathbb{D}_{4}(0, T)$), by the weak convergence and Lebesgue's dominated convergence arguments we have

Consequently, $\tilde{\phi}^\epsilon = \partial^{\alpha}_{0^+}\phi^\epsilon$ in the weak sense. In the same way we prove, in the weak sense, that $(\partial^{\alpha}_{0^+}\varphi^\epsilon_{i}, \partial^{\alpha}_{0^+}\psi^\epsilon, \partial^{\beta}_{0^+}u^\epsilon) = (\tilde{\varphi}^\epsilon_{i}, \tilde{\psi}^\epsilon, \tilde{u}^{\epsilon})$.

Consider $\kappa\in \mathcal{D}(]0, T[)$ and $({\mathbf h}_n, {\mathbf k}_n, {\mathbf \rho}_n)\in W_{i, n}\times W_{g, n}\times W_{i, n}$. According to (4.8), we can deduce that

According to (4.21), (4.39) and to density properties of spaces spanned by $(w_l)$ and $(e_m)$, and using similar arguments to derive (4.40), it is easy to pass to limit ($n\to \infty$) in linear terms. For the nonlinear terms, we take into account assumption ${\rm {\mathbf (H2)}_{4}}$ and we use the classical technique based on taking the difference between the sequence and its limit in the form of sum of two terms such that the first uses weak convergence result and the other uses strong convergence result. So we omit the details. The limit $(\phi^\epsilon, \varphi_i^\epsilon, \psi^\epsilon, u^\epsilon)$ then satisfies the following system (for all $({\mathbf h}, {\mathbf k}, \rho)\in \mathbb{V}_{\rm H}\times \mathbb{V}\times \mathbb{V}_{\rm H}$)

In the case of $\alpha > 1/2$ and $\beta > 1/2$, the continuity of solution at $t = 0$ and the equalities $\varphi_i^{\epsilon}(0^{+}) = \varphi_{i, 0}$, $\phi^{\epsilon}(0^{+}) = \phi_0$, $\psi^{\epsilon}(0^{+}) = \psi_0$ and $u^{\epsilon}(0^{+}) = u_0$, is a consequence of Lemma 3.2. This completes the proof.

We can now show the well-posedness of memory bidomain-torso (degenerate) problem (1.1).

4.3.   Degenerate problem study

Similar to derive results (4.19), (4.22), (4.29) and (4.38), we have the following a priori estimates for problem (4.4).

Lemma 4.1. Assume the conditions (H1) and ${\rm {\mathbf (H2)}_{4}}$ hold, $f_{\rm H}\in L^{2/\alpha_1}(0, T;L^{2}(\Omega))$ and $u_{0}\in L^{3}(\Omega_{\rm H})$.

(i) If $\phi_{0}, \varphi_{i, 0}$ and $\psi_{0}$ are $L^{2}$-functions, then there exists a constant $C > 0$ (independent on $\epsilon$) such that (for all $t\in (0, T)$)

(ii) If $\phi_{0}, \varphi_{i, 0}$ and $\psi_{0}$ are $H^{1}$-functions, then there exists a constant $C > 0$ (independent on $\epsilon$) such that (for all $t\in (0, T)$)

We are now ready to prove existence of weak solutions to system (1.1).

Theorem 4.2. Let assumptions (H1) and ${\rm {\mathbf (H2)}_{4}}$ be fulfilled and $\beta\geq \alpha$. Then for $(\phi_{0}, u_0)$ and $f_{\rm H}$ given such that $(\phi_{0}, u_0)\in \mathbb{V}_{\rm H}\times L^{3}(\Omega_{\rm H})$ and $f_{\rm H}\in L^{\frac{2}{\alpha_1}}(0, T;L^{2}(\Omega_{\rm H}))$ with $0 < \alpha_1 < \alpha$, there exists a weak solution $(\varphi_{i}, \phi, \psi, u)$ to problem (1.1) verifying

Moreover, if $\beta > 1/2$ then $u\in C([0, T];L^{2}(\Omega_{\rm H}))$ and $u(0^{+}) = u_0$ and if $\alpha > 1/2$ then $\phi\in C([0, T];L^{2}(\Omega_{\rm H}))$ and $\phi(0^{+}) = \phi_0$.

Proof. From Lemma 4.1, similar arguments to derive (4.40), Theorem 3.1 and compactness argument, it follows that there exist $X = (\phi, \varphi_i, \psi)$ and $u$ such that there exists a subsequence of $(X^{\epsilon}; u^{\epsilon}) = (\phi^\epsilon, \varphi_{i}^\epsilon, \psi^\epsilon; u^\epsilon)$ also denoted by $(X^{\epsilon}; u^{\epsilon})$ such that

Thus, using (4.45) and similar argument to derive (4.41), we can deduce the existence of weak solution $(\phi, \varphi_i, \psi, u)$ satisfying regularity (4.44).

Finally, in case of $\alpha > 1/2$ and $\beta > 1/2$, the continuity of solution at $t = 0$ and the equalities $\phi^{\epsilon}(0^{+}) = \phi_0$, and $u^{\epsilon}(0^{+}) = u_0$, is a consequence of Lemma 3.2. This completes the proof..

Theorem 4.3. Let assumptions (H1), ${\rm {\mathbf (H2)}_{4}}$ and (H3) be fulfilled and $\beta\geq \alpha$. Then the solution $(\varphi_i, \phi, \psi, u)$ to system (1.1) is unique. Moreover, let $(\phi^{(i)}_{0}, u^{(i)}_0, f^{(i)}_{\rm H})$ be given such that $(\phi_{0}^{(i)}, u_0^{(i)})\in \mathbb{V}_{\rm H}\times L^{3}(\Omega_{\rm H})$ and $f^{(i)}_{\rm H}\in L^{\frac{2}{\alpha_1}}(0, T;L^{2}(\Omega_{\rm H}))$, for $i = 1, 2$, with $0 < \alpha_1 < \alpha$. For $(\varphi^{(j)}_i, \phi^{(j)}, \psi^{(j)}, u^{(j)})$ a weak solutions to system (1.1), which corresponds to data $(\phi^{(j)}_0, u^{(j)}_0, f^{(j)}_{\rm H})$ (for $j = 1, 2$), the following Lipschitz continuity relationship is satisfied.

(i) If $\alpha = \beta$, we have the following relation

(ii) If $\alpha < \beta$ and if the following assumption holds

the relation (4.46) is also obtained.

Proof. As uniqueness result is a consequence of relation (4.46) (see further), we start by showing this relation.

Since $(\varphi^{(j)}_i, \phi^{(j)}, \psi^{(j)}, u^{(j)})$ is a weak solution to system (1.1), for $j = 1, 2$, then $(\varphi_{i}, \phi, \psi, u) = (\varphi_{i}^{(1)}-\varphi_{i}^{(2)}, \phi^{(1)}-\phi^{(2)}, \psi^{(1)}-\psi^{(2)}, u^{(1)}-u^{(2)})$ satisfies

where $(\phi_0, u_0, f_{\rm H}) = (\phi^{(1)}_0-\phi^{(2)}_0, u^{(1)}_0-u^{(2)}_0, f_{\rm H}^{(1)}-f_{\rm H}^{(2)})$.

Then (since $\psi = \varphi_e$ and $\phi = \varphi_i-\varphi_e$ in $\Omega_{\rm H}$),

(i) If $\alpha = \beta$, according to (H3), we have

and then

So (using Hölder inequality),

By Lemma 3.5, we can deduce that first

and then, from (4.51)

(ii) Assume now that $\alpha < \beta$, then, according to (H3), we get

On account of assumption (4.47), Hölder's inequality, continuous embedding of $H^{1}$ in $L^{4}$ and the fact that $\phi_i\in L^{\infty}(0, T, H^{1}(\Omega_{\rm H}))$, for $i = 1, 2$, we have

So

with the constants $\eta_1$ and $\eta_2$ chosen appropriately. According to estimates (4.55) and by choosing $\eta_1 = \frac{1}{\mathfrak{c}_{\alpha}}\min(\nu_i, \nu_g)$, estimate (4.54) becomes

with $\Psi(t) = || \nabla \varphi_i || ^{2}_{L^{2}(\Omega_{\rm H})} +|| \nabla \psi|| ^{2}_{L^{2}(\Omega)}$ and $\zeta = \frac{1}{\mathfrak{c}_{\alpha}}\min(\nu_i, \nu_g)$.

Consequently, by choosing $\eta_2 < \frac{\Gamma(\beta)\zeta}{T^{\beta-\alpha}\Gamma(\alpha)}$, we can deduce first from Corollary 2 that

and then, from (4.56)

Finally, if the difference between data is null i.e., $\phi_0 = 0$, $u_0 = 0$ and $f_{\rm H} = 0$, we can conclude from (4.53) (or (4.57)), (1.4) and (1.6) that $(\phi, \varphi_i, \psi, u) = (0, 0, 0, 0)$ and then the uniqueness of solution. This completes the proof.

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 $\mathcal{I}$ and $\mathcal{G}$ as (in which the state variables are dimensionless and scaled)

where $0 < \phi_{gate} < 1$, $\tau_{in}$, $\tau_{out}$, $\tau_{open}$ and $\tau_{close}$ are given positive constants. The cubic function $\phi^2(1-\phi)$ describes the voltage dependence of the inward current. This model is well-known to be valid under the assumption $0 < \tau_{in}\ll\tau_{out}\ll\min (\tau_{open}, \tau_{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 $\mathcal{G}$

where the differentiable function $0\leq h_\zeta \leq 1$ is given by

with $\zeta$ a positive parameter.

The operator $\mathcal{G}_{\zeta}(\phi, u)$ can be written as $\mathcal{G}_{\zeta}(\phi, u) = \mathcal{I}_{2, \zeta}(\phi)+\hbar_{\zeta}(\phi) u$ where

According to the definition of $\tanh$, we can deduce that $\lim\limits_{\zeta\rightarrow 0} h_{\zeta}(\phi) = \left\{ \begin{array}{ll} 1 { if }~~ \phi < \phi_{gate}, \\ 0 { if }~~ \phi > \phi_{gate} \end{array} \right.$ and then $\lim\limits_{\zeta\rightarrow 0} \mathcal{G}_\zeta(\phi, u) = \mathcal{G}(\phi, u)$. The regularized Mitchell-Schaeffer model has a slightly different structure compared to models in (2.3) because in this model, $\hbar_{\zeta}$ depend on $\phi$ through the function $h_\zeta$. Since $h_\zeta$ 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 $\mathcal{G}$ in the form of $\mathcal{G}(.; \phi, u) = \mathcal{I}_{2}(.; \phi)+\hbar(.; \phi) 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 $\mathcal{G}$ Carathéodory function from $(\Omega \times \mbox{IR}^{{}}) \times \mbox{IR}^{{2}}$ into $\mbox{IR}^{{}}$ and locally Lipschitz continuous function on $(\phi, u)$ and, $\mathcal{I}_{2}$ and $\hbar$ sufficiently regulars.

5.   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 $\phi$ 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 $\beta\geq \alpha$, in order to show the impact of $\beta$ on the evolution of system, we will also present some results for $\alpha$ fixed to 1 and $\beta$ 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 $\beta$ at 1 and we vary $\alpha$, and in a second time we fix $\alpha$ at 1 and we vary $\beta$ 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 $\partial_{0^+}^{\gamma}$, $\gamma\in]0, 1]$. In this paper, we consider the following approximation (for sufficiently regular function $\psi$ at time $t\in]t_{n-1}, t_{n}]$ ($n = 1, 2, \ldots$), with $t_k = k\Delta t$ for $k = 0, 1, \ldots$ and given time step $\Delta t$)

where

It is clear that $I^{1}_{memory}\;(\psi_{history})$, for $\gamma = 1$, is a null function. The state function $\psi_{history}$ is corresponding to the prior history of $\psi$.

Note that the previous fractional-order dynamic is the sum of first-order dynamics and an operator $I^{\gamma}_{memory}$ representing hypothetical memory effects.

The term $I^{\gamma}_{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^{\gamma}_{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.

^‡ Thus the depolarizing effect of the stimulus becomes repolarizing and vice versa

Here, we perform various numerical experiments to investigate the impact of different values of $\phi$ fractional-order $\alpha$ and $u$ fractional-order $\beta$ across the time in $(0, T)$ and we present some numerical results corresponding to the middle point of heart domain $\Omega_H$ for $\phi$, $\varphi_e$ and $u$, and in the middle point of torso domain $\Omega_B$ for $\varphi_s$. We assume that the domain heart-torso is in a square region $\Omega = [-L_T/2, -L_T/2]\times [-L_T/2, -L_T/2]$ and the domain $\Omega_H$ corresponds to an approximate shape of heart, where $L_T$ is the torso length and $L_H$ is the characteristic length of $\Omega_H$ (Figure 3).

The boundary of domain $\Omega_H$ is defined by (for $(x, y)\in \Gamma_H$) : $\vert x-a\vert^3 + (y-b -c|x-d|^2)^2-e = 0,$ where $a, b, c, d, e$ are given positive constants. Here we choose $a = 0.5$, $b = 3.5$, $c = 0.5$, $c = 1$, $d = 0.5$ and $e = 1/15$.

We use some biophysical parameters which are presented in Table 1 and we fix $L_H = 15cm$ and $L_T = 45cm$. In those both applications, we impose the following initial conditions

Finally, we consider $T = 1$, $\mathcal{K}_{i} = \mathcal{K}_{e} = 0.003 s.cm^{-1}{\mathbf I}_d$ (with ${\mathbf 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 $\chi_{H}$, $\chi_{[a_i, T_{act}+a_i]}$ and $\chi_{prop}$ are defined by

The last characteristic is associated to wave propagation as a diagonal which evolves from to left-bottom side to right-top side of $\Omega_H$. The function $\Phi$, 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 $\alpha = \beta = 1$) is shown in black and the last curve drawn is in dotted line.

5.1.   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 $\mathcal{I}$ and $\mathcal{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 $\beta$ to 1 and we vary $\alpha$). 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 $\phi$, $\varphi_e$ and $u$, at midpoint of heart, for different values of $\alpha$ with $\beta = 1$. We observe first that we have always the same kind of behavior. Second, we find that the amplitude varies depending on $\alpha$ values. For state $\varphi_e$, the effect fractional-order dynamic is negligible except for the extremum of $\varphi_e$ (as the value of $\alpha$ becomes smaller, the value of the minimum point becomes smaller too).

5.2.   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 $\mathcal{I}$ and $\mathcal{G}$ as:

These operators depend on the change-over voltage $\phi_{gate}$, the resting potential $\phi_{min}$, the maximum potential $\phi_{max}$, and on times constants $\tau_{in}$, $\tau_{out}$, $\tau_{open}$ and $\tau_{close}$, the two times $\tau_{open}$ and $\tau_{close}$, respectively controlling the durations of the action potential and of the recovery phase, and the two times $\tau_{in}$ and $\tau_{out}$, respectively controlling the length of depolarization and repolarization phases. These constants are such that $\tau_{in} < \tau_{out} < \min (\tau_{open}, \tau_{close}).$

First, we investigate the influence of membrane capacitive memory on action potential properties in absence of ionic variable memory (i.e., we fix $\beta$ to 1 and we vary $\alpha$). 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 $\phi$, $\varphi_e$ and $u$, at midpoint of heart and $\varphi_s$ at at midpoint of torso, for different values of $\alpha$ with $\beta = 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 $\alpha$ decreases, and we notice therefore an increase in the rate of spontaneous activity (when $\alpha$ decreases). We observe the same kind of behaviors for the other state variables than for $\phi$, with an acceleration of spontaneous activity when $\alpha$ decreases.

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 $\phi$, $\varphi_e$ and $u$, at midpoint of heart and $\varphi_s$ at at midpoint of torso, for different values of $\beta$ with $\alpha = 1$. We observe that the smaller the value of $\beta$ 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.

6.   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.).

^§ from some observations, cardiac MRI-Measurement or desired target

Conflict of interest

The author declares there is no conflicts of interest in this paper.