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

Aziz Belmiloudi

► To cite this version:

Aziz Belmiloudi. Cardiac memory phenomenon, time-fractional order nonlinear system and bidomain-torso type model in electrocardiology. AIMS Mathematics, 2021, 6 (1), pp.821-867. 10.3934/math.2021050 . hal-03093927

HAL Id: hal-03093927

https://hal.science/hal-03093927

Submitted on 25 May 2021

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research article

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

Aziz Belmiloudi*

Institut de Recherche Mathématique de Rennes (IRMAR), Université Européenne de Bretagne, 20 avenue Buttes de Coësmes, CS 70839, 35708 Rennes Cédex 7, France

* Correspondence: Email: aziz.belmiloudi@math.cnrs.fr.

**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, ionic-current 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 fractional-order 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.

**Keywords:** fractional-order dynamics; heart-torso coupling; cellular heterogeneity; integral inequality; weak solution; cardiac memory; memory induced T-wave; ionic models

**Mathematics Subject Classification:** 35R11, 34A08, 35K57, 35B30, 92C50, 35Q92
# 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.

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

<footnote>*The so-called implantable cardioverter defibrillators</footnote>
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 phenomenon 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 \le 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.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The derived “*Memory bidomain model*” is defined on heart domain $\Omega_H$, while $\Omega_B$ is the rest of body.</figcaption></figure>

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
$$ 
\begin{array}{l}
c_{\alpha} \partial_{0^{+}}^{\alpha} \phi + \mathcal{I}(.; \phi, u) - \operatorname{div}(\mathcal{K}_{i} \nabla \varphi_{i}) = I_{i}, \quad \text{in } Q_{H} = \Omega_{H} \times (0, T) \\
c_{\alpha} \partial_{0^{+}}^{\alpha} \phi + \mathcal{I}(.; \phi, u) + \operatorname{div}(\mathcal{K}_{e} \nabla \varphi_{e}) = -I_{e}, \quad \text{in } Q_{H} \\
-\operatorname{div}(\mathcal{K}_{s} \nabla \varphi_{s}) = 0, \quad \text{in } Q_{B} = \Omega_{B} \times (0, T)
\end{array}
\qquad (1.1) $$

$$ \partial_{0^{+}}^{\beta} u + \mathcal{G}(.; \phi, u) = 0, \text{ in } Q_{H} $$

subject to initial and boundary conditions (1.4),

or equivalently

$$ 
\begin{array}{l}
c_{\alpha} \partial_{0^{+}}^{\alpha} \phi + \mathcal{I}(.; \phi, u) - \operatorname{div}(\mathcal{K}_{i} \nabla \phi) = \operatorname{div}(\mathcal{K}_{i} \nabla \varphi_{e}) + I_{i}, \quad \text{in } Q_{H} \\
-\operatorname{div}((\mathcal{K}_{e} + \mathcal{K}_{i}) \nabla \varphi_{e}) = \operatorname{div}(\mathcal{K}_{i} \nabla \phi) + I, \quad \text{in } Q_{H} \\
-\operatorname{div}(\mathcal{K}_{s} \nabla \varphi_{s}) = 0, \quad \text{in } Q_{B}
\end{array}
\qquad (1.2) $$

$$ \partial_{0^{+}}^{\beta} u + \mathcal{G}(.; \phi, u) = 0, \text{ in } Q_{H} $$

subject to initial and boundary conditions (1.4).

Here, $\partial_{0^{+}}^{\alpha}$ and $\partial_{0^{+}}^{\beta}$ 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_H$ which is a bounded open subset, and by $\Gamma_H = \partial\Omega_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_B$ with a piecewise smooth boundary $\Gamma_B = \Gamma_H \cup \Gamma_T$ where $\Gamma_T$ is the thorax surface. The whole domain is denoted by the $\Omega = \bar{\Omega}_H \cup \Omega_B$. In $\Omega_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_B$, $\varphi_s$ is thoracic medium electric potential. The parameter $c_\alpha$ is $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 assimilated 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^\alpha = 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)$^\dagger$ 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$.

<footnote>$^\dagger$which traditionally describes the recovery of APD as a function of the interbeat interval</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Modeling of the membrane as resistor and non-ideal capacitor coupled in parallel (where $0 < \alpha \le 1$ is the fractional order of capacitor).</figcaption></figure>

The coupling between equations in $\Omega_H$ and equation in $\Omega_B$ (in systems (1.1) and (1.2)) is operate at the heart-torso interface $\Gamma_H$. The continuity conditions at the interface can be given by

$$ \varphi_e = \varphi_s, \quad \text{in } \Sigma_H = \Gamma_H \times (0, T) \\ (\mathcal{K}_e \nabla \varphi_e) \cdot \mathbf{n} = (\mathcal{K}_s \nabla \varphi_s) \cdot \mathbf{n}, \quad \text{on } \Sigma_H. \tag{1.3} $$

where $\mathbf{n}$ is the outward normal to $\Gamma_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 $\mathcal{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 $\mathcal{I}$ and $\mathcal{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

initial conditions

$$ \phi(., t = 0) = \phi_0, \quad u(., t = 0) = u_0, \quad \text{in } \Omega_H $$

and boundary conditions

$$ \varphi_e = \varphi_s, \quad \text{in } \Sigma_H \tag{1.4} $$

$$ (\mathcal{K}_e \nabla \varphi_e) \cdot \mathbf{n} = (\mathcal{K}_s \nabla \varphi_s) \cdot \mathbf{n}, \quad \text{on } \Sigma_H $$

$$ (\mathcal{K}_i \nabla \varphi_i) \cdot \mathbf{n} = 0, \quad \text{on } \Sigma_H $$

$$ (\mathcal{K}_s \nabla \varphi_s) \cdot \mathbf{n}_T = 0, \quad \text{on } \Sigma_T = \Gamma_T \times (0, T) $$

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

$$ \int_{\Omega_H} I d\mathbf{x} = 0. \qquad (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 $\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))

$$ \int_{\Gamma_H} \varphi_e d\mathbf{x} = \int_{\Gamma_H} \varphi_s d\mathbf{x} = p_H \quad \text{where } p_H \text{ is a fixed time-dependent function.} \qquad (1.6) $$

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

2. The functions $K_i, K_e, H$ and $G$ depend on the fiber extension ratio.

3. We can suppose, for example, that $f_H$ is only a time dependent source and is of the form

$$ f_H(\mathbf{x}, t) = \theta(t)(\chi_{\Omega_H^{(1)}}(\mathbf{x}) - \chi_{\Omega_H^{(2)}}(\mathbf{x})), \qquad (1.7) $$

where $\chi_{\Omega_H^{(i)}}$ is the characteristic function of set $\Omega_H^{(i)}$, $i = 1, 2$. The support regions $\Omega_H^{(1)}$ and $\Omega_H^{(2)}$ can be considered to represent an anode (positive electrode) and a cathode (negative electrode) respectively. $\square$

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.

**2. Assumptions, notations and some fundamental inequalities**

Let $\mathcal{U} \subset \mathbb{R}^m$, $m \ge 1$, be an open and bounded set with a smooth boundary and $\mathcal{U}_T = \mathcal{U} \times (0, T)$. We use the standard notation for Sobolev spaces (see [1]), denoting the norm of $W^{q,p}(\mathcal{U})$ ($q \in \mathbb{N}$, $p \in [1, \infty)$) by $\| \cdot \|_{W^{q,p}}$. In the special case $p = 2$, we use $H^q(\mathcal{U})$ instead of $W^{q,2}(\mathcal{U})$. The duality pairing of a Banach space $X$ with its dual space $X'$ is given by $\langle \cdot, \cdot \rangle_{X',X}$. For a Hilbert space $Y$ the inner product is denoted by $(\cdot, \cdot)_Y$ and the inner product in $L^2(\Omega)$ is denoted by $(\cdot, \cdot)$. For any pair of real numbers $r, s \ge 0$, we introduce the Sobolev space $H^{r,s}(\mathcal{U}_T)$ defined by $H^{r,s}(\mathcal{U}_T) = L^2(0, T; H^r(\mathcal{U})) \cap H^s(0, T; L^2(\mathcal{U}))$, which is a Hilbert space normed by $(\| v \|_{L^2(0,T;H^r(\mathcal{U}))}^2 + \| v \|_{H^s(0,T;L^2(\mathcal{U}))}^2)^{1/2}$, where $H^s(0, T; L^2(\mathcal{U}))$ denotes the Sobolev space of order $s$ of functions defined on $(0, T)$ and taking values in $L^2(\mathcal{U})$, and defined (see [52]) by $H^s(0, T; L^2(\mathcal{U})) = [H^q(0, T, L^2(\mathcal{U})), L^2(\mathcal{U}_T)]_\theta$, for $s = (1 - \theta)q$ with $\theta \in (0, 1)$ and $q \in \mathbb{N}$, and $H^q(0, T; L^2(\mathcal{U})) = \{v \in L^2(\mathcal{U}_T) | \frac{\partial^j v}{\partial t^j} \in L^2(\mathcal{U}_T), \text{ for } 1 \le j \le q\}$. For a given Banach space $X$, with norm $\| \cdot \|_X$, of functions integrable on $\mathcal{U}$, we define its subspace $X|_{\mathbb{R}} = \{u \in X, \int_{\mathcal{U}} u = 0\}$ that is a Banach space with norm $\| \cdot \|_X$, and we denote by $[u]$ the projection of $u \in X$ on $X|_{\mathbb{R}}$ such that $[u] = u - \frac{1}{\text{mes}(\mathcal{U})} \int_{\mathcal{U}} u dx$ (with $\text{mes}(\mathcal{U})$ standing for Lebesgue measure of $\mathcal{U}$).

**Lemma 2.1.** (Poincaré-Wirtinger inequality) Assume that $1 \le p \le \infty$ and that $\mathcal{U}$ is a bounded connected open subset of $\mathbb{R}^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}(\mathcal{U})$, we have

$$
\| [u] \|_{L^p(\mathcal{U})} \le C \| \nabla u \|_{L^p(\mathcal{U})}.
$$

**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(\mathcal{U})|_{\mathbb{R}}$. $\square$

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

(i) $H^q(\mathcal{U}) \subset L^p(\mathcal{U})$, $\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

$$ \| v \|_{L^p(\mathcal{V})} \le C \| v \|_{H^q(\mathcal{V})}^\theta \| v \|_{L^2(\mathcal{V})}^{1-\theta}, \forall v \in H^q(\mathcal{V}), $$

where $0 \le \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). $\square$

Finally, we introduce the spaces:

* $\mathbb{H}_H = L^2(\Omega_H)$, $\mathbb{V}_H = H^1(\Omega_H)$, $\mathbb{V} = H^1(\Omega)$ (endowed with their usual norms) and $\mathbb{U}_H = \mathbb{V}_H|_{\mathbb{R}}$,

* $\mathbb{U}_{HB} = \{\psi \in H^1(\Omega) \mid \int_{\Omega_H} \psi dx = 0\}$.

We will denote by $\mathbb{V}'_H$ (resp. $\mathbb{U}'_H$) the dual of $\mathbb{V}_H$ (resp. of $\mathbb{U}_H$). We have the following continuous embeddings, where $p \ge 2$ if $d = 2$ and $2 \le p \le 6$ if $d = 3$, $p'$ is such that $\frac{1}{p'} + \frac{1}{p} = 1$

$$ \begin{aligned} & \mathbb{V}_H \subset \mathbb{H}_H \subset \mathbb{V}'_H, \mathbb{U}_H \subset \mathbb{H}_H|_{\mathbb{R}} \subset \mathbb{U}'_H, \\ & \mathbb{V}_H \subset L^p(\Omega_H) \subset \mathbb{H}_H \equiv (\mathbb{H}_H)' \subset L^{p'}(\Omega_H) \subset \mathbb{V}'_H \end{aligned} \quad (2.1) $$

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

$$ \mathbb{D}_p(0, T) = L^p(\mathbb{Q}_H) \cap L^2(0, T; \mathbb{V}_H), \text{ and its dual } \mathbb{D}'_p(0, T) = L^{p'}(\mathbb{Q}_H) + L^2(0, T; \mathbb{V}'_H) \subset L^{p'}(0, T; \mathbb{V}'_H). $$

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

**Definition 2.1.** A real valued function $\mathcal{H}$ defined on $D \times \mathbb{R}^q$, $q \ge 1$, is a Carathéodory function iff $\mathcal{H}(.; v)$ is measurable for all $v \in \mathbb{R}^q$ and $\mathcal{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 \le \Pi_{i=1,k} \| f_i \|_{L^{q_i}(D)}$,

where $\| f_i \|_{L^{q_i}(D)} = (\int_D |f_i|^{q_i} dx)^{1/q_i}$ and $\sum_{1 \le i \le k} \frac{1}{q_i} = 1$.

(ii) Young's inequality ($\forall a, b > 0$ and $\epsilon > 0$): $ab \le \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:

$$ \left[ \int_{\Omega} \left( \int_{0}^{t} |f(x, s)| ds \right)^{p} dx \right]^{1/p} \le \int_{0}^{t} \left( \int_{\Omega} |f(x, s)|^{p} dx \right)^{1/p} ds, \text{ for } p \in ]1, +\infty[ \text{ and } t > 0. $$

Finally, we denote by $\mathcal{L}(A; B)$ the set of linear and continuous operators from a vectorial space A into a vectorial space B, and by $\mathcal{R}^*$ the adjoint operator to a linear operator $\mathcal{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 $K_{\theta} \in W^{1,\infty}(\bar{\Omega}_{\text{H}})$, $\theta \in \{i, e\}$ and $K_{s} \in W^{1,\infty}(\bar{\Omega}_{\text{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 \mathbb{R}^d)$

$$
\begin{align}
K_1 \|\psi\|^2 &\le \psi^T \mathcal{K}_{\theta} \psi \le K_2 \|\psi\|^2 \quad \text{in } \bar{\Omega}_{\text{H}}, \nonumber \\
K_1 \|\psi\|^2 &\le \psi^T \mathcal{K}_{s} \psi \le K_2 \|\psi\|^2 \quad \text{in } \bar{\Omega}_{\text{B}}. \tag{2.2}
\end{align}
$$

**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\le k\le d}$ in $\mathbb{R}^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}) = \text{diag}((\lambda_{i,k})_{1\le k\le d})$ and $\Lambda_e(\mathbf{x}) = \text{diag}((\lambda_{e,k})_{1\le k\le d})$.

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

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$)

$$
\begin{align}
\mathcal{I}(\mathbf{x}, t; \phi, u) &= \mathcal{I}_0(\mathbf{x}, t; \phi) + \mathcal{I}_1(\mathbf{x}, t; \phi)u, \tag{2.3} \\
\mathcal{G}(\mathbf{x}, t; \phi, u) &= \mathcal{I}_2(\mathbf{x}, t; \phi) + \hbar(\mathbf{x}, t)u, \nonumber
\end{align}
$$

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.

(H2)$_p$ The operators $\mathcal{I}_0$, $\mathcal{I}_1$ and $\mathcal{I}_2$ are Carathéodory functions from $(\Omega \times \mathbb{R}) \times \mathbb{R}$ into $\mathbb{R}$ and continuous on $\phi$ (as in [12]). Furthermore, for some $p \ge 2$ if $d = 2$ and $p \in [2, 6]$ if $d = 3$, the following requirements hold

(i) there exist constants $\beta_i \ge 0$ $(i = 1, \dots, 10)$ such that for any $v \in \mathbb{R}$

$$
\begin{align}
|\mathcal{I}_0(\cdot, v)| &\le \beta_1 + \beta_2|v|^{p-1}, \nonumber \\
|\mathcal{I}_1(\cdot, v)| &\le \beta_3 + \beta_4|v|^{p/2-1}, \nonumber \\
|\mathcal{I}_2(\cdot, v)| &\le \beta_5 + \beta_6|v|^{p/2}, \nonumber \\
|E_g(\cdot, v)| &\le \beta_7 + \beta_8|v|^p, |g_2(\cdot, v)| \le \beta_9 + \beta_{10}|v|^{p-2}, \tag{2.4}
\end{align}
$$

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

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

$$
\begin{align}
\mu_1 v \mathcal{I}(\cdot, v, w) + w \mathcal{G}(\cdot, v, w) &\ge \mu_2 |v|^p - \mu_3 (\mu_1 |v|^2 + |w|^2) - \mu_4, \tag{2.5} \\
E_g(\cdot, v) &\ge \mu_5 |v|^p - \mu_6 |v|^2 - \mu_7. \nonumber
\end{align}
$$

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}: \mathbb{R}^2 \to \mathbb{R}^2$ defined by

$$
F_{\mu}(\cdot; \mathbf{v}) = \begin{pmatrix} \mu(\mathcal{I}(\cdot; \mathbf{v})) \\ \mathcal{G}(\cdot; \mathbf{v}) \end{pmatrix}, \quad \forall \mathbf{v} = (v, w) \in \mathbb{R}^2, \qquad (2.6)
$$
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 \mathbb{R}^2$, $i = 1, 2$)

$$ (F_{\mu}(.; \mathbf{v}_1) - F_{\mu}(.; \mathbf{v}_2)) \cdot (\mathbf{v}_1 - \mathbf{v}_2) \geq -C_L \|\mathbf{v}_1 - \mathbf{v}_2\|^2. \quad (2.7) $$

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

$$ Q_{\mu}(\phi, u) = \frac{1}{2} (\nabla F_{\mu}(\phi, u)^{T} + \nabla F_{\mu}(\phi, u)). $$

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

$$ C_F \leq \lambda_1(\phi, u) \leq \lambda_2(\phi, u), \quad (2.8) $$

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

**Lemma 2.3.** ( [9]) Let assumptions (2.3), (H1) and (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

$$ \begin{aligned} \|\mathcal{I}(., t; \phi, u)\|_{L^{p'}(\Omega_H)} &\le C_1 + C_2 \|\phi\|_{L^{p'}(\Omega_H)}^{p/p'} + C_3 \|u\|_{\mathbb{H}}^{2/p'}, \\ \|\mathcal{G}(., t; \phi, u)\|_{L^2(\Omega_H)} &\le C_4 + C_5 \|\phi\|_{L^{p'}(\Omega_H)}^{p/2} + C_6 \|u\|_{\mathbb{H}}, \end{aligned} \quad (2.9) $$

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 (H1), (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

(a) for RM type model : $\mathcal{I}_1(.; v) = b_2(.)v, \mathcal{I}_2(.; v) = -b_3(.)v,$

(b) for FHN type model : $\mathcal{I}_1(.; v) = b_2(.), \mathcal{I}_2(.; v) = -b_3(.)v,$

(c) for LAP type model : $\mathcal{I}_1(.; v) = b_2(.)v, \mathcal{I}_2(.; v) = b_3(.)v(r + 1 - v),$

where $b_i \in W^{1,\infty}(Q)$, $i = 1, 3$, are sufficiently regular functions from $Q$ into $\mathbb{R}^{+,*}$ 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 $\mathbb{R}$ we have

$$ \mathcal{I}_0(.; v_1) - \mathcal{I}_0(.; v_2) = b_1(v_1 - v_2) (v_1^2 + v_2^2 + v_1v_2 - (r+1)(v_1 + v_2) + r) \text{ and} $$

(a) for RM type model : $\mathcal{I}_1(.; v_1) - \mathcal{I}_1(.; v_2) = b_2(v_1 - v_2), \mathcal{I}_2(.; v_1) - \mathcal{I}_2(.; v_2) = -b_3(v_1 - v_2),$

(b) for FHN type model : $\mathcal{I}_1(.; v_1) - \mathcal{I}_1(.; v_2) = 0, \mathcal{I}_2(.; v_1) - \mathcal{I}_2(.; v_2) = -b_3(v_1 - v_2),$

(c) for LAP type model : $\mathcal{I}_1(.; v_1) - \mathcal{I}_1(.; v_2) = b_2(v_1 - v_2),$

$$ \mathcal{I}_2(.; v_1) - \mathcal{I}_2(.; v_2) = b_3(v_1 - v_2)((r+1) - v_1 - v_2). $$

For the sake of simplicity, we shall write $\mathcal{I}_i(\psi)$, $\mathcal{I}(\psi, v)$ and $\mathcal{G}(\psi, v)$ in place of $\mathcal{I}_i(\mathbf{x}, t; \psi)$, $\mathcal{I}(\mathbf{x}, t; \psi, v)$ and $\mathcal{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))$

$$ 
\begin{aligned}
I_{a^+}^{\gamma} [f](t) &= \frac{1}{\Gamma(\gamma)} \int_{a}^{t} (t-\tau)^{\gamma-1} f(\tau) d\tau, \\
I_{T^{-}}^{\gamma} [f](t) &= \frac{1}{\Gamma(\gamma)} \int_{t}^{T} (\tau-t)^{\gamma-1} f(\tau) d\tau,
\end{aligned} 
\quad (3.1) $$

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

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

$$ I_{a^+}^{\gamma_1} [I_{a^+}^{\gamma_2} [f]] = I_{a^+}^{\gamma_1 + \gamma_2} [f] \quad (3.2) $$

and holds for a $L^p$-function $f$ ($1 \le p \le \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))$

$$ 
\begin{aligned}
D_{a^+}^{\gamma} [f](t) &= \frac{d}{dt} (I_{a^+}^{1-\gamma} [f](t)) = \frac{1}{\Gamma(1-\gamma)} \frac{d}{dt} \left( \int_a^t (t-\tau)^{-\gamma} f(\tau) d\tau \right), \\
\partial_{a^+}^{\gamma} [f](t) &= I_{a^+}^{1-\gamma} \left[ \frac{df(t)}{dt} \right] = \frac{1}{\Gamma(1-\gamma)} \int_a^t (t-\tau)^{-\gamma} \frac{df}{d\tau}(\tau) d\tau.
\end{aligned} 
\quad (3.3) $$

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

$$ f(t) = f(a) + I_{a^+}^{\gamma} [\partial_{a^+}^{\gamma} f](t) = f(a) + \frac{1}{\Gamma(\gamma)} \int_{a}^{t} (t-\tau)^{\gamma-1} \partial_{a^+}^{\gamma} f(\tau) d\tau. \quad (3.4) $$
**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))$

$$ 
\begin{aligned}
D_{T^{-}}^{\gamma} f(t) &= -\frac{d}{dt}(I_{T^{-}}^{1-\gamma}[f](t)) = -\frac{1}{\Gamma(1-\gamma)}\frac{d}{dt}\left(\int_{t}^{T}(\tau-t)^{-\gamma}f(\tau)d\tau\right), \\
\partial_{T^{-}}^{\gamma}[f](t) &= -I_{T^{-}}^{1-\gamma}\left[\frac{df(t)}{dt}\right] = -\frac{1}{\Gamma(1-\gamma)}\int_{t}^{T}(\tau-t)^{-\gamma}\frac{df}{dt}(\tau)d\tau.
\end{aligned}
\quad (3.5) $$

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))$

$$ 
\begin{aligned}
I_{-\infty}^{\gamma}[f](t) &= \frac{1}{\Gamma(\gamma)} \int_{-\infty}^{t} (t-\tau)^{\gamma-1} f(\tau) d\tau, \\
\partial_{-\infty}^{\gamma} f(t) &= I_{-\infty}^{1-\gamma} \left[ \frac{df}{dt}(t) \right] = \frac{1}{\Gamma(1-\gamma)} \int_{-\infty}^{t} (t-\tau)^{\gamma-1} \frac{df}{dt}(\tau) d\tau.
\end{aligned}
\quad (3.6) $$

**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_{a^+}^\gamma f(t) = \frac{k}{\Gamma(1-\gamma)} \frac{d}{dt} \left( \int_a^t (t-\tau)^{\gamma-1} d\tau \right) = \frac{k(t-a)^{\gamma-1}}{\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

$$ 
\begin{aligned}
D_{a^+}^{\gamma} f(t) &= \partial_{a^+}^{\gamma} f(t) + \frac{f(a)(t-a)^{-\gamma}}{\Gamma(1-\gamma)}, \\
D_{T^{-}}^{\gamma} f(t) &= \partial_{T^{-}}^{\gamma} f(t) + \frac{f(T)(T-t)^{-\gamma}}{\Gamma(1-\gamma)}.
\end{aligned}
\quad (3.7) $$

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

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

(ii) $L^p(a, T)$ into $L^p(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_{a^+}^\gamma 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 \le 1$ and $p, q \ge 1$ with $1/p + 1/q \le 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

$$ \int_a^T (f(\tau), I_{a^+}^\gamma [g](\tau))_X d\tau = \int_a^T (g(\tau), I_{T^-}^\gamma [f](\tau))_X d\tau, $$

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

$$ \int_a^T (f(\tau), D_{a^+}^\gamma g(\tau))_X d\tau = \int_a^T (g(\tau), D_{T^-}^\gamma f(\tau))_X d\tau. $$

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

$$ (i) \quad \int_a^T (\partial_{a^+}^\gamma f(\tau), g(\tau))_X d\tau = - \int_a^T (D_{T^-}^\gamma g(\tau), f(\tau))_X d\tau + |(I_{T^-}^{1-\gamma} [g](\tau), f(\tau))_X|_a^T, $$

$$ (ii) \quad \int_a^T (D_{a^+}^\gamma f(\tau), g(\tau))_X d\tau = - \int_a^T (D_{T^-}^\gamma g(\tau), f(\tau))_X d\tau + (I_{T^-}^{1-\gamma} [g](T^-), f(T))_X, $$

$$ (iii) \quad \int_a^T (D_{T^-}^\gamma f(\tau), g(\tau))_X d\tau = - \int_a^T (D_{a^+}^\gamma g(\tau), f(\tau))_X d\tau - (I_{a^+}^{1-\gamma} [g](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

$$ f(t) \le h(t) + b(t) I_{0^+}^\gamma [f](t). $$

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

$$ f(t) \le h(t) + \int_0^t \sum_{k=1}^\infty h(\tau)(t-\tau)^{k\gamma-1} \frac{(b(t))^k}{\Gamma(k\gamma)} d\tau. $$

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

$$ f(t) \le h(t) E_{\gamma,1}(b(t)t^\gamma). $$

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

$$ E_{\theta_1, \theta_2}(z) = \sum_{k=0}^{\infty} \frac{1}{\Gamma(k\theta_1 + \theta_2)} z^k. \qquad (3.8) $$

The function $E_{\theta_1, \theta_2}$ is an entire function of the variable $z$ for any $\theta_1, \theta_2 \in \mathbb{C}$, $\mathrm{Re}(\theta_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 \hookrightarrow X \hookrightarrow X_1 \text{ being continuous and } X_0 \hookrightarrow X \text{ is compact.} \qquad (3.9) $$

The Fourier transform of $f: \mathbb{R} \to X_1$ is defined by $\hat{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 $\mathbb{R}$ with compact support. Then

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

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

We define the Hilbert space $W^{\gamma}(\mathbb{R}; \mathcal{X}_0, \mathcal{X}_1)$, for a given $\gamma > 0$, by

$$ W^{\gamma}(\mathbb{R}; \mathcal{X}_0, \mathcal{X}_1) = \{v \in L^2(\mathbb{R}, \mathcal{X}_0) \mid \partial_{-\infty}^{\gamma} f \in L^2(\mathbb{R}, \mathcal{X}_1)\}, $$

endowed with the norm

$$ \| v \|_{W^{\gamma}} = \left( \| v \|_{L^2(\mathbb{R}, \mathcal{X}_0)} + \| \tau \rvert^{\gamma} \widehat{v} \|_{L^2(\mathbb{R}, \mathcal{X}_1)} \right)^{1/2}. $$

For any subset $K$ of $\mathbb{R}$, we define the subspace $W_K^{\gamma}$ of $W^{\gamma}$ by

$$ W_K^{\gamma}(\mathbb{R}; \mathcal{X}_0, \mathcal{X}_1) = \{v \in W^{\gamma}(\mathbb{R}; \mathcal{X}_0, \mathcal{X}_1) \mid \text{support of } v \subset K\}. $$

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 $W_K^{\gamma}(\mathbb{R}; \mathcal{X}_0, \mathcal{X}_1)$ into $L^2(\mathbb{R}; \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 \le \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)$)

$$ \begin{aligned} \partial_{0^+}^{\gamma_1} f(t) &\le h_1(t) + c_1(f(t) + g(t)) + \eta \Psi(t), \\ \partial_{0^+}^{\gamma_2} g(t) &\le h_2(t) + c_2(f(t) + g(t)), \end{aligned} \quad (3.10) $$

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

$$ F(t) \le F_0 E_{\gamma_2, 1}(d_1 t^{\gamma_2}) + \sum_{k=0}^{\infty} d_1^k I_{0^+}^{k\gamma_2+\gamma_1} [h_1](t) + \sum_{k=0}^{\infty} d_1^k I_{0^+}^{(k+1)\gamma_2} [h_2](t), \quad (3.11) $$

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

The function $E_{\gamma}$ is the Mittag-Leffler function which is defined by (3.8).
Proof. From (3.4) we can deduce that

$$
\begin{align}
\Gamma(\gamma_1)(f(t) - f_0) &\le \int_0^t (t-\tau)^{\gamma_1-1} h_1(\tau) d\tau + c_1 \int_0^t (t-\tau)^{\gamma_1-1} (f(\tau) + g(\tau)) d\tau \nonumber \\
&\quad + \eta \int_0^t (t-\tau)^{\gamma_1-1} \Psi(\tau) d\tau, \tag{3.12} \\
\Gamma(\gamma_2)(g(t) - g_0) &\le \int_0^t (t-\tau)^{\gamma_2-1} h_2(\tau) d\tau + c_2 \int_0^t (t-\tau)^{\gamma_2-1} (f(\tau) + g(\tau)) d\tau \nonumber \\
&\quad - \xi \int_0^t (t-\tau)^{\gamma_2-1} \Psi(\tau) d\tau. \nonumber
\end{align}
$$

Put $F = f + g$, then

$$
\begin{align}
(\Gamma(\gamma_1)\Gamma(\gamma_2))(F(t) - F_0) &\le \int_0^t (\Gamma(\gamma_2)(t-\tau)^{\gamma_1-1}h_1(\tau) + \Gamma(\gamma_1)(t-\tau)^{\gamma_2-1}h_2(\tau))d\tau \nonumber \\
&\quad + \int_0^t (\eta\Gamma(\gamma_2)(t-\tau)^{\gamma_1-1} - \xi\Gamma(\gamma_1)(t-\tau)^{\gamma_2-1})\Psi(\tau)d\tau \tag{3.13} \\
&\quad + \int_0^t (c_1\Gamma(\gamma_2)(t-\tau)^{\gamma_1-1} + c_2\Gamma(\gamma_1)(t-\tau)^{\gamma_2-1})F(\tau)d\tau. \nonumber
\end{align}
$$

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

$$
\begin{align}
(\Gamma(\gamma_1)\Gamma(\gamma_2))(F(t) - F_0) &\le \Gamma(\gamma_2) \int_0^t (t-\tau)^{\gamma_1-1} h_1(\tau) d\tau + \Gamma(\gamma_1) \int_0^t (t-\tau)^{\gamma_2-1} h_2(\tau) d\tau \nonumber \\
&\quad + (\eta \Gamma(\gamma_2) T^{(\gamma_1-\gamma_2)} - \xi \Gamma(\gamma_1)) \int_0^t (t-\tau)^{\gamma_2-1} \Psi(\tau) d\tau \tag{3.14} \\
&\quad + (c_1 \Gamma(\gamma_2) T^{(\gamma_1-\gamma_2)} + c_2 \Gamma(\gamma_1)) \int_0^t (t-\tau)^{\gamma_2-1} F(\tau) d\tau. \nonumber
\end{align}
$$

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

$$
\begin{align}
(\Gamma(\gamma_1)\Gamma(\gamma_2))(F(t) - F_0) &\le \Gamma(\gamma_2) \int_0^t (t-\tau)^{\gamma_1-1} h_1(\tau) d\tau + \Gamma(\gamma_1) \int_0^t (t-\tau)^{\gamma_2-1} h_2(\tau) d\tau \nonumber \\
&\quad + (c_1 \Gamma(\gamma_2) T^{(\gamma_1-\gamma_2)} + c_2 \Gamma(\gamma_1)) \int_0^t (t-\tau)^{\gamma_2-1} F(\tau) d\tau. \tag{3.15}
\end{align}
$$

Thus

$$
F(t) \le F_0 + I_{0+}^{\gamma_1} [h_1] (t) + I_{0+}^{\gamma_2} [h_2] (t) + d_1 I_{0+}^{\gamma_2} [F] (t), \quad (3.16)
$$

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

Finally, from Lemma 3.5 and relation (3.2), we get (since $I_{0+}^{k\gamma_2} [F] (0) = F_0 I_{0+}^{k\gamma_2} [.1] = \frac{F_0 I_{0+}^{k\gamma_2}}{\Gamma(1+k\gamma_2)}$)

$$
F(t) \le F_0 \sum_{k=0}^{\infty} \frac{1}{\Gamma(1+k\gamma_2)} (d_1 t^{\gamma_2})^k + \sum_{k=0}^{\infty} d_1^k I_{0+}^{k\gamma_2+\gamma_1} [h_1](t) + \sum_{k=0}^{\infty} d_1^k I_{0+}^{(k+1)\gamma_2} [h_2](t).
$$

This completes the proof.

<figure><img src="image_5.png" /></figure>
**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

$$F(t) \leq F_0 E_{\gamma_2,1}(d_1 t^{\gamma_2}) + t^{\gamma_1} h_1(t) E_{\gamma_2,\gamma_1+1}(d_1 t^{\gamma_2}) + t^{\gamma_2} h_2(t) E_{\gamma_2,\gamma_2+1}(d_1 t^{\gamma_2}), \quad (3.17)$$

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

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

$$I_{0^+}^{\gamma_2+\gamma_i}[h_i](t) \leq h_i(t) \frac{1}{\Gamma(k\gamma_2 + \gamma_i)} \int_0^t (t-\tau)^{k\gamma_2+\gamma_i-1} d\tau = h_i(t) \frac{1}{\Gamma(k\gamma_2 + \gamma_i + 1)} t^{k\gamma_2+\gamma_i}.$$

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

**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

$$F(t) \leq F(0)E_{\gamma_2,1}(d_1 t^{\gamma_2}) + \frac{t^{\gamma_1-\alpha_1}}{r_1} E_{\gamma_2,\gamma_1}(d_1 t^{\gamma_2}) \| h_1 \|_{L^{1/\alpha_1}(0,t)} \\ + \frac{t^{\gamma_2-\alpha_2}}{r_2} E_{\gamma_2,\gamma_2}(d_1 t^{\gamma_2}) \| h_2 \|_{L^{1/\alpha_2}(0,t)}, \quad (3.18)$$

$$\text{where } r_i = \left(\frac{\gamma_i - \alpha_i}{1 - \alpha_i}\right)^{1-\alpha_i}, \text{ for } i = 1, 2, \quad F = f + g, \quad F_0 = f_0 + g_0, \quad d_0 = \frac{T^{(\gamma_1 - \gamma_2)}\Gamma(\gamma_2)}{\Gamma(\gamma_1)} \quad \text{and} \\ d_1 = \frac{c_1\Gamma(\gamma_2)T^{(\gamma_1 - \gamma_2)} + c_2\Gamma(\gamma_1)}{\Gamma(\gamma_1)} = c_1d_0 + c_2.$$

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

$$\begin{aligned} \Gamma(k\gamma_2 + \gamma_i) I_{0^+}^{\gamma_2+\gamma_i}[h_i](t) &\leq \left(\int_0^t (t-\tau)^{\frac{\gamma_2+\gamma_i-1}{1-\alpha_i}} d\tau\right)^{1-\alpha_i} \left(\int_0^t h_i(\tau)^{\frac{1}{\alpha_i}} d\tau\right)^{\alpha_i} \\ &\leq t^{k\gamma_2} \left(\int_0^t (t-\tau)^{\frac{\gamma_i-1}{1-\alpha_i}} d\tau\right)^{1-\alpha_i} \left(\int_0^t h_i(\tau)^{\frac{1}{\alpha_i}} d\tau\right)^{\alpha_i} = \frac{t^{k\gamma_2+\gamma_i-\alpha_i}}{r_i} \left(\int_0^t h_i(\tau)^{\frac{1}{\alpha_i}} d\tau\right)^{\alpha_i}, \end{aligned} \quad (3.19)$$

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

According to (3.11), we have

$$F(t) \leq F_0 E_{\gamma_2,1}(d_1 t^{\gamma_2}) + \sum_{k=0}^{\infty} d_1^k \frac{t^{k\gamma_2+\gamma_1-\alpha_1}}{\Gamma(k\gamma_2+\gamma_1)r_1} \left(\int_0^t h_1(\tau)^{\frac{1}{\alpha_i}} d\tau\right)^{\alpha_1} \\ + \sum_{k=0}^{\infty} d_1^k \frac{t^{k\gamma_2+\gamma_2-\alpha_2}}{\Gamma(k\gamma_2+\gamma_2)r_2} \left(\int_0^t h_2(\tau)^{\frac{1}{\alpha_i}} d\tau\right)^{\alpha_2}. \quad (3.20)$$

Thus

$$F(t) \leq F_0 E_{\gamma_2,1}(d_1 t^{\gamma_2}) + \frac{t^{\gamma_1-\alpha_1}}{r_1} E_{\gamma_2,\gamma_1}(d_1 t^{\gamma_2}) \| h_1 \|_{L^{1/\alpha_1}(0,t)} \\ + \frac{t^{\gamma_2-\alpha_2}}{r_2} E_{\gamma_2,\gamma_2}(d_1 t^{\gamma_2}) \| h_2 \|_{L^{1/\alpha_2}(0,t)}. \quad (3.21)$$

This completes the proof. $\square$

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_H$, and the thoracic potential in $\Omega_B$, by

$$ \psi = \varphi_e \text{ in } \Omega_H \text{ and } \psi = \varphi_s \text{ in } \Omega_B. \quad (4.1) $$

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

$$ K_g = K_e \text{ in } \bar{\Omega}_H \text{ and } K_g = K_s \text{ in } \bar{\Omega}_B. $$

We now define the following forms

$$ \begin{aligned} A_i(w, v) &= \int_{\Omega_H} K_i \nabla w \cdot \nabla v d\mathbf{x}, & A_e(w, v) &= \int_{\Omega_H} K_e \nabla w \cdot \nabla v d\mathbf{x}, \\ A_s(w, v) &= \int_{\Omega_B} K_s \nabla w \cdot \nabla v d\mathbf{x}, & A_g(w, v) &= \int_{\Omega} K_g \nabla w \cdot \nabla v d\mathbf{x}. \end{aligned} \quad (4.2) $$

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

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

*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. $\square$

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}_H)$, $\phi \in D_4(0, T) \cap L^\infty(0, T; L^2(\Omega_H))$, $\partial_i^\alpha \phi \in D'_4(0, T)$, $\psi \in L^2(0, T; \mathbb{V})$, $u \in L^\infty(0, T; L^3(\Omega_H))$ and $\partial_i^\beta u \in L^2(0, T; L^2(\Omega_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_H$)

$$ \begin{aligned} \langle \epsilon_a \partial_{0+}^\alpha \phi, v_i \rangle_{V'_H, V_H} + \int_{\Omega_H} \mathcal{I}(.; \phi, u) v_i d\mathbf{x} + A_i(\varphi_i, v_i) &= \int_{\Omega_H} f_H v_i d\mathbf{x}, \\ \langle \epsilon_a \partial_{0+}^\alpha \phi, v \rangle_{V'_H, \mathbb{V}} + \int_{\Omega_H} \mathcal{I}(.; \phi, u) v d\mathbf{x} - A_g(\psi, v) &= \int_{\Omega_H} f_H v d\mathbf{x}, \quad (4.3) \\ (\partial_{0+}^\beta u, \rho)_{L^2(\Omega_H)} + \int_{\Omega_H} \mathcal{G}(.; \phi, u) \rho d\mathbf{x} &= 0, \end{aligned} $$

for all $v_i \in V_H$, $v \in \mathbb{V}$, $\rho \in V_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_H$)

$$ 
\begin{aligned}
& c_{\alpha} \partial_{0+}^{\alpha} \phi^{\epsilon} + \epsilon \partial_{0+}^{\alpha} \varphi_{i}^{\epsilon} + \mathcal{I}(\phi^{\epsilon}, u^{\epsilon}) - \operatorname{div}(\mathcal{K}_{i} \nabla \varphi_{i}^{\epsilon}) = f_{H}, \quad \text{in } Q_{H} \\
& c_{\alpha} \partial_{0+}^{\alpha} \phi^{\epsilon} - \epsilon \partial_{0+}^{\alpha} \varphi_{e}^{\epsilon} + \mathcal{I}(\phi^{\epsilon}, u^{\epsilon}) + \operatorname{div}(\mathcal{K}_{e} \nabla \varphi_{e}^{\epsilon}) = f_{H}, \quad \text{in } Q_{H} \\
& \epsilon \partial_{0+}^{\alpha} \varphi_{s}^{\epsilon} - \operatorname{div}(\mathcal{K}_{s} \nabla \varphi_{s}^{\epsilon}) = 0, \quad \text{in } Q_{B} \\
& \partial_{0+}^{\beta} u^{\epsilon} + \mathcal{G}(\phi^{\epsilon}, u^{\epsilon}) = 0, \quad \text{in } Q_{H}
\end{aligned}
\quad (4.4) $$

supplemented with the same initial and boundary conditions as in (1.4).

The auxiliary initial conditions for $\varphi_i^\epsilon$ and $\varphi_s^\epsilon$, 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}_H$, $v \in \mathbb{V}$, $\rho \in \mathbb{V}_H$ and a.e. $t \in (0, T)$)

$$ 
\begin{aligned}
& \langle c_{\alpha} \partial_{0+}^{\alpha} \phi^{\epsilon}, v_i \rangle_{\mathbb{V}'_{H}, \mathbb{V}_{H}} + \epsilon \langle \partial_{0+}^{\alpha} \varphi_{i}^{\epsilon}, v_i \rangle_{\mathbb{V}'_{H}, \mathbb{V}_{H}} + \int_{\Omega_H} \mathcal{I}(\phi^{\epsilon}, u^{\epsilon}) v_i d\mathbf{x} + A_i(\varphi_i^{\epsilon}, v_i) = \int_{\Omega_H} f_H v_i d\mathbf{x}, \\
& \langle c_{\alpha} \partial_{0+}^{\alpha} \phi^{\epsilon}, v \rangle_{\mathbb{V}'_{H}, \mathbb{V}_{H}} - \epsilon \langle \partial_{0+}^{\alpha} \psi^{\epsilon}, v \rangle_{\mathbb{V}'_{H}, \mathbb{V}} + \int_{\Omega_H} \mathcal{I}(\phi^{\epsilon}, u^{\epsilon}) v d\mathbf{x} - \mathbf{A}_g(\psi, v) = \int_{\Omega_H} f_H v d\mathbf{x}, \\
& (\partial_{0+}^{\beta} u^{\epsilon}, \rho)_{L^2(\Omega_H)} + \int_{\Omega_H} \mathcal{G}(\phi^{\epsilon}, u^{\epsilon}) \rho d\mathbf{x} = 0,
\end{aligned}
\quad (4.5) $$

with the initial condition

$$ (\phi^{\epsilon}, \varphi_{i}^{\epsilon}, \psi^{\epsilon}, u^{\epsilon})(t = 0) = (\phi_{0}, \varphi_{i,0}, \psi_{0}, u_{0}), $$

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 (H2)$_4$ be fulfilled, $\beta \ge \alpha$ and $0 < \epsilon \le 1$ be given. Then for $(\varphi_{i,0}, \phi_0, \psi_0, u_0)$ and $f_H$ given such that $(\varphi_{i,0}, \phi_0, \psi_0, u_0) \in \mathbb{V}_H^2 \times \mathbb{V} \times L^3(\Omega_H)$ and $f_H \in L^{\frac{2}{\alpha_1}}(0, T; L^2(\Omega_H))$ with $0 < \alpha_1 < \alpha$, there exists a solution $(\varphi_i^\epsilon, \phi^\epsilon, \psi^\epsilon, u^\epsilon)$ of problem (4.5) verifying

$$ 
\begin{aligned}
& \phi^{\epsilon} \in L^{\infty}(0, T; \mathbb{V}_H), \quad \partial_{0+}^{\alpha} \phi^{\epsilon} \in L^2(Q_H), \\
& \varphi_{i}^{\epsilon} \in L^{\infty}(0, T; \mathbb{V}_H), \quad \sqrt{\epsilon} \partial_{0+}^{\alpha} \varphi_{i}^{\epsilon} \in L^2(Q_H), \\
& \psi^{\epsilon} \in L^{\infty}(0, T; \mathbb{V}), \quad \sqrt{\epsilon} \partial_{0+}^{\alpha} \phi^{\epsilon} \in L^2(Q_H), \\
& u^{\epsilon} \in L^{\infty}(0, T; L^2(\Omega_H)), \quad \partial_{0+}^{\beta} u^{\epsilon} \in L^2(Q_H).
\end{aligned}
\quad (4.6) $$

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 \ge 1}$ be a Hilbert basis and orthogonal in $L^2$ of $V_H$, $(f_k)_{k \ge 1}$ be a Hilbert basis of $U_H$, $(g_k)_{k \ge 1}$ be a Hilbert basis of $W_B = \{v \in H^1(\Omega_B) \mid 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_H} = 0$ and by $(e_k)_{k \ge 1}$ a Galerkin basis of $U_{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 \mathbb{N}^*$, we denote by $W_{i,n} = \operatorname{span}(w_1, \dots, w_n)$, $W_{e,n} = \operatorname{span}(f_1, \dots, f_n)$, $W_{s,n} = \operatorname{span}(g_1, \dots, g_n)$ and $W_{g,n} = \operatorname{span}(e_1, \dots, e_{2n})$ the spaces generated, respectively, by $(w_k)_{n \ge k \ge 1}$, $(f_k)_{n \ge k \ge 1}$, $(g_k)_{n \ge k \ge 1}$ and $(e_k)_{2n \ge k \ge 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_n^\epsilon, \phi_{i,n}^\epsilon, \psi_n^\epsilon, u_n^\epsilon)$ of the problem (4.4). Setting

$$ 
\begin{align} 
\phi_n^\epsilon(\cdot, t) &= \sum_{l=1}^n \bar{\omega}_{n,l}^\epsilon(t)w_l, \quad u_n^\epsilon(\cdot, t) = \sum_{l=1}^n \nu_{n,l}^\epsilon(t)w_l, \tag{4.7} \\
\varphi_{i,n}^\epsilon(\cdot, t) &= \sum_{l=1}^n \bar{\omega}_{i,n,l}^\epsilon(t)w_l, \quad \psi_n^\epsilon(\cdot, t) = \sum_{l=1}^{2n} \pi_{n,l}^\epsilon(t)e_l 
\end{align} 
$$

where $\omega_i^\epsilon = (\bar{\omega}_{i,n,l}^\epsilon)_{l=1,n}$, $\omega^\epsilon = (\bar{\omega}_{n,l}^\epsilon)_{l=1,n}$, $\vartheta^\epsilon = (\nu_{n,l}^\epsilon)_{l=1,n}$ and $\Pi^\epsilon = (\pi_{n,m}^\epsilon)_{m=1,2n}$ are unknown functions, and replacing $(\phi^\epsilon, \phi_i^\epsilon, \psi^\epsilon, u^\epsilon)$ by $(\phi_n^\epsilon, \phi_{i,n}^\epsilon, \psi_n^\epsilon, u_n^\epsilon)$ 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_H$, $\phi_n^\epsilon = \phi_{i,n}^\epsilon - \psi_n^\epsilon$)

$$ 
\begin{align} 
(c_\alpha + \epsilon) \int_{\Omega_H} \partial_{0+}^\alpha \varphi_{i,n}^\epsilon w_l d\mathbf{x} - c_\alpha \int_{\Omega_H} \partial_{0+}^\alpha \psi_n^\epsilon w_l d\mathbf{x} 
&= - \int_{\Omega_H} \mathcal{I}(\varphi_{i,n}^\epsilon - \psi_n^\epsilon, u_n^\epsilon) w_l d\mathbf{x} - A_i(\varphi_{i,n}^\epsilon, w_l) + \int_{\Omega_H} f_H w_l d\mathbf{x} \nonumber \\
&:= \langle F_n^{\epsilon,i}, w_l \rangle \stackrel{\text{def}}{=} F_l^i(t; \omega_i^\epsilon, \Pi^\epsilon, \vartheta^\epsilon), 
\end{align} 
$$
$$ 
\begin{align} 
-c_\alpha \int_{\Omega_H} \partial_{0+}^\alpha \varphi_{i,n}^\epsilon e_m d\mathbf{x} + c_\alpha \int_{\Omega_H} \partial_{0+}^\alpha \psi_n^\epsilon e_m d\mathbf{x} + \epsilon \int_{\Omega} \partial_{0+}^\alpha \psi_n^\epsilon e_m d\mathbf{x} 
&= \int_{\Omega_H} \mathcal{I}(\varphi_{i,n}^\epsilon - \psi_n^\epsilon, u_n^\epsilon) e_m d\mathbf{x} - A_g(\psi_n^\epsilon, e_m) - \int_{\Omega_H} f_H e_m d\mathbf{x} \tag{4.8} \\
&:= \langle G_n^\epsilon, e_m \rangle \stackrel{\text{def}}{=} G_m(t; \omega_i^\epsilon, \Pi^\epsilon, \vartheta^\epsilon), 
\end{align} 
$$

$$ 
\begin{align} 
\int_{\Omega_H} \partial_{0+}^\beta u_n^\epsilon w_l d\mathbf{x} &= - \int_{\Omega_H} \mathcal{G}(\varphi_{i,n}^\epsilon - \psi_n^\epsilon, u_n^\epsilon) w_l d\mathbf{x} \\
&:= \langle H_n^\epsilon, w_l \rangle \stackrel{\text{def}}{=} H_l(t; \omega_i^\epsilon, \Pi^\epsilon, \vartheta^\epsilon), 
\end{align} 
$$

with the initial condition

$$ (\phi_n^\epsilon, \phi_{n,i}^\epsilon, \psi_n^\epsilon, u_n^\epsilon)(t=0) = (L_{i,n}\phi_0, L_{i,n}\varphi_{i,0}, L_{g,n}\psi_0, L_{i,n}u_0) $$

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)

$$ (L_{i,n}\phi_0, L_{i,n}\varphi_{i,0}, L_{g,n}\psi_0, L_{i,n}u_0) \xrightarrow{n \to \infty} (\phi_0, \varphi_{i,0}, \psi_0, u_0) \text{ strongly in } \mathbb{V}_H^2 \times \mathbb{V} \times L^3(\Omega_H). \quad (4.9) $$
The function $F_n^{\epsilon,i}$, $G_n^{\epsilon}$ and $H_n^{\epsilon}$ are defined by (for $(v, w, \rho) \in \mathbb{V}_H \times \mathbb{V} \times \mathbb{V}_H$)

$$ 
\begin{aligned}
\langle F_n^{\epsilon,i}, v \rangle &= - \int_{\Omega_H} \mathcal{I}(\phi_n^{\epsilon}, u_n^{\epsilon}) v d\mathbf{x} - A_i(\varphi_{i,n}^{\epsilon}, v) + \int_{\Omega_H} f_H v d\mathbf{x}, \\
\langle G_n^{\epsilon}, w \rangle &= \int_{\Omega_H} \mathcal{I}(\phi_n^{\epsilon}, u_n^{\epsilon}) w d\mathbf{x} - \mathbf{A}_g(\psi_n^{\epsilon}, w) - \int_{\Omega_H} f_H w d\mathbf{x}, \\
\langle H_n^{\epsilon}, \rho \rangle &= - \int_{\Omega_H} \mathcal{G}(\phi_n^{\epsilon}, u_n^{\epsilon}) \rho d\mathbf{x}.
\end{aligned} 
\quad (4.10) $$

**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 $(\omega_i^{\epsilon}, \Pi^{\epsilon}, \vartheta^{\epsilon})$

$$ 
\begin{aligned}
\mathcal{M}_n(\partial_{0+}^{\alpha} \omega_i^{\epsilon}, \partial_{0+}^{\alpha} \Pi^{\epsilon}, \partial_{0+}^{\beta} \vartheta^{\epsilon})^T &= \mathcal{F}_n(t; \omega_i^{\epsilon}, \Pi^{\epsilon}, \vartheta^{\epsilon}), \\
(\omega_i^{\epsilon}, \Pi^{\epsilon}, \vartheta^{\epsilon})(t=0) &= (\omega_{i,0}^{\epsilon}, \Pi_0^{\epsilon}, \vartheta_0^{\epsilon}),
\end{aligned} 
\quad (4.11) $$

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

$$ \mathcal{F}_n(t; \omega_i^{\epsilon}, \Pi^{\epsilon}, \vartheta^{\epsilon}) = ((F_l^i(t; \omega_i^{\epsilon}, \Pi^{\epsilon}, \vartheta^{\epsilon}))_{l=1,n}, (G_m(t; \omega_i^{\epsilon}, \Pi^{\epsilon}, \vartheta^{\epsilon}))_{m=1,2n}, (H_l(t; \omega_i^{\epsilon}, \Pi^{\epsilon}, \vartheta^{\epsilon}))_{l=1,n})^T. $$

Thanks to assumptions (H2)$_4$ the function $\mathcal{F}_n$ is Carathéodory function.

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

$$ 
\begin{aligned}
(\partial_{0+}^{\alpha} \omega_i^{\epsilon}, \partial_{0+}^{\alpha} \Pi^{\epsilon}, \partial_{0+}^{\beta} \vartheta^{\epsilon})^T &= \mathcal{S}_n(t; \omega_i^{\epsilon}, \Pi^{\epsilon}, \vartheta^{\epsilon}), \\
(\omega_i^{\epsilon}, \Pi^{\epsilon}, \vartheta^{\epsilon})(t=0) &= (\omega_{i,0}^{\epsilon}, \Pi_0^{\epsilon}, \vartheta_0^{\epsilon}),
\end{aligned} 
\quad (4.12) $$

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

The existence of a local absolutely continuous function $(\omega_i^{\epsilon}, \Pi^{\epsilon}, \vartheta^{\epsilon})$ on an interval $[0, \tilde{T}]$, with $\tilde{T} \in ]0, T]$ is insured by the standard FODE theory see [32, 33] (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 "~" on $T$. Now, we set

$$ \mathbf{h}_n(\cdot, t) = \sum_{k=1}^{n} \vartheta_{n,k}(t)w_k, \quad \mathbf{v}_n(\cdot, t) = \sum_{l=1}^{n} \zeta_{n,l}(t)w_l, \quad \mathbf{e}_n(\cdot, t) = \sum_{p=1}^{2n} \pi_{n,p}(t)e_p. $$

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

$$
\begin{align}
& (c_{\alpha} + \epsilon) \int_{\Omega_H} \partial_{0^+}^{\alpha} \varphi_{i,n}^{\epsilon} \mathbf{h}_n d\mathbf{x} - c_{\alpha} \int_{\Omega_H} \partial_{0^+}^{\alpha} \psi_n^{\epsilon} \mathbf{h}_n d\mathbf{x} \nonumber \\
& = - \int_{\Omega_H} \mathcal{I}(\varphi_{i,n}^{\epsilon} - \psi_n^{\epsilon}, u_n^{\epsilon}) \mathbf{h}_n d\mathbf{x} - A_i(\varphi_{i,n}^{\epsilon}, \mathbf{h}_n) + \int_{\Omega_H} f_H \mathbf{h}_n d\mathbf{x}, \nonumber \\
& - c_{\alpha} \int_{\Omega_H} \partial_{0^+}^{\alpha} \varphi_{i,n}^{\epsilon} \mathbf{e}_n d\mathbf{x} + c_{\alpha} \int_{\Omega_H} \partial_{0^+}^{\alpha} \psi_n^{\epsilon} \mathbf{e}_n d\mathbf{x} + \epsilon \int_{\Omega} \partial_{0^+}^{\alpha} \psi_n^{\epsilon} \mathbf{e}_n d\mathbf{x} \nonumber \\
& = \int_{\Omega_H} \mathcal{I}(\varphi_{i,n}^{\epsilon} - \psi_n^{\epsilon}, u_n^{\epsilon}) \mathbf{e}_n d\mathbf{x} - \mathbf{A}_g(\psi_n^{\epsilon}, \mathbf{e}_n) - \int_{\Omega_H} f_H \mathbf{e}_n d\mathbf{x}, \tag{4.13} \\
& \int_{\Omega_H} \partial_{0^+}^{\beta} u_n^{\epsilon} \mathbf{v}_n d\mathbf{x} = - \int_{\Omega_H} \mathcal{G}(\varphi_{i,n}^{\epsilon} - \psi_n^{\epsilon}, u_n^{\epsilon}) \mathbf{v}_n d\mathbf{x}, \nonumber
\end{align}
$$

with the initial condition

$$
(\phi_n^\epsilon, \phi_{n,i}^\epsilon, \psi_n^\epsilon, u_n^\epsilon)(t=0) = (L_{i,n}\phi_0, L_{i,n}\phi_{i,0}, L_{g,n}\psi_0, L_{i,n}u_0).
$$

Taking **h**<sub>n</sub> = φ<sub>i,n</sub><sup>ε</sup>, **e**<sub>n</sub> = ψ<sub>n</sub><sup>ε</sup> and **v**<sub>n</sub> = u<sub>n</sub><sup>ε</sup> and using the uniform coercivity of the forms A<sub>i</sub> and **A**<sub>g</sub>, we can deduce (since ψ<sub>n</sub><sup>ε</sup> = φ<sub>e,n</sub><sup>ε</sup> and φ<sub>n</sub><sup>ε</sup> = φ<sub>i,n</sub><sup>ε</sup> − φ<sub>e,n</sub><sup>ε</sup> in Ω<sub>H</sub>)

$$
\begin{align}
& \frac{1}{2} \left( c_{\alpha} \partial_{0^+}^{\alpha} \| \phi_n^{\epsilon} \|_{L^2(\Omega_H)}^2 + \epsilon (\partial_{0^+}^{\alpha} \| \psi_n^{\epsilon} \|_{L^2(\Omega)}^2 + \partial_{0^+}^{\alpha} \| \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2) \right) \nonumber \\
& \qquad + v_i \| \nabla \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2 + v_g \| \nabla \psi_n^{\epsilon} \|_{L^2(\Omega)}^2 + \int_{\Omega_H} \mathcal{I}(\phi_n^{\epsilon}, u_n^{\epsilon}) \phi_n^{\epsilon} d\mathbf{x} \nonumber \\
& \leq \int_{\Omega_H} f_H \phi_n^{\epsilon} d\mathbf{x}, \tag{4.14} \\
& \frac{1}{2} \partial_{0^+}^{\beta} \| u_n^{\epsilon} \|_{L^2(\Omega_H)}^2 \leq - \int_{\Omega_H} \mathcal{G}(\phi_n^{\epsilon}, u_n^{\epsilon}) u_n^{\epsilon} d\mathbf{x}. \nonumber
\end{align}
$$

Then the following inequality

$$
\begin{align}
& \frac{1}{2} \left( c_{\alpha} \partial_{0^+}^{\alpha} \| \phi_n^{\epsilon} \|_{L^2(\Omega_H)}^2 + \epsilon (\partial_{0^+}^{\alpha} \| \psi_n^{\epsilon} \|_{L^2(\Omega)}^2 + \partial_{0^+}^{\alpha} \| \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2) \right) \nonumber \\
& \quad + v_i \| \nabla \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2 + v_g \| \nabla \psi_n^{\epsilon} \|_{L^2(\Omega)}^2 + \frac{1}{\mu_1} \int_{\Omega_H} (\mu_1 \mathcal{I}(\phi_n^{\epsilon}, u_n^{\epsilon}) \phi_n^{\epsilon} + \mathcal{G}(\phi_n^{\epsilon}, u_n^{\epsilon}) u_n^{\epsilon}) d\mathbf{x} \nonumber \\
& \leq \int_{\Omega_H} f_H \phi_n^{\epsilon} d\mathbf{x} - \frac{1}{\mu_1} \int_{\Omega_H} \mathcal{G}(\phi_n^{\epsilon}, u_n^{\epsilon}) u_n^{\epsilon} d\mathbf{x}, \tag{4.15} \\
& \frac{1}{2} \partial_{0^+}^{\beta} \| u_n^{\epsilon} \|_{L^2(\Omega_H)}^2 \leq - \int_{\Omega_H} \mathcal{G}(\phi_n^{\epsilon}, u_n^{\epsilon}) u_n^{\epsilon} d\mathbf{x}. \nonumber
\end{align}
$$

Since, from (H2)$_4$, $\mu_1 \mathcal{I}(\phi_n^\epsilon, u_n^\epsilon)\phi_n^\epsilon + \mathcal{G}(\phi_n^\epsilon, u_n^\epsilon)u_n^\epsilon \ge \mu_2 |\phi_n^\epsilon|^4 - \mu_3(\mu_1 |\phi_n^\epsilon|^2 + |u_n^\epsilon|^2) - \mu_4$ and $|\mathcal{G}(\phi_n^\epsilon, u_n^\epsilon)u_n^\epsilon| \le c_1(1 + |\phi_n^\epsilon|^2|u_n^\epsilon| + |u_n^\epsilon|^2)$, we obtain (from Hölder's inequality)

$$
\begin{align}
& c_{\alpha} \partial_{0^+}^{\alpha} \| \phi_n^{\epsilon} \|_{L^2(\Omega_H)}^2 + \epsilon (\partial_{0^+}^{\alpha} \| \psi_n^{\epsilon} \|_{L^2(\Omega)}^2 + \partial_{0^+}^{\alpha} \| \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2) \nonumber \\
& \quad + 2v_i \| \nabla \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2 + 2v_g \| \nabla \psi_n^{\epsilon} \|_{L^2(\Omega)}^2 + \frac{\mu_2}{\mu_1} \int_{\Omega_H} |\phi_n^{\epsilon}|^4 d\mathbf{x} \nonumber \\
& \le c_2 + c_3 \| f_H \|_{L^2(\Omega_H)}^2 + c_4 (\| \phi_n^{\epsilon} \|_{L^2(\Omega_H)}^2 + \| u_n^{\epsilon} \|_{L^2(\Omega_H)}^2), \tag{4.16} \\
& \partial_{0^+}^{\beta} \| u_n^{\epsilon} \|_{L^2(\Omega_H)}^2 \le c_5 + c_6 (\| \phi_n^{\epsilon} \|_{L^2(\Omega_H)}^2 + \| u_n^{\epsilon} \|_{L^2(\Omega_H)}^2) + \eta \int_{\Omega_H} |\phi_n^{\epsilon}|^4 d\mathbf{x}, \nonumber
\end{align}
$$
with the constants $\eta$ chosen appropriately. In particular we have

$$ \partial_{0+}^{\alpha} \Psi_{n}^{\epsilon} + \frac{\mu_2}{\mu_1} \int_{\Omega_H} |\phi_n^{\epsilon}|^4 dx \leq c_7(1 + \| f_H \|_{L^2(\Omega_H)}^2) + c_8(\Psi_n^{\epsilon} + U_n^{\epsilon}), \\ \partial_{0+}^{\beta} U_n^{\epsilon} \leq c_9 + c_{10}(\Psi_n^{\epsilon} + U_n^{\epsilon}) + \eta \int_{\Omega_H} |\phi_n^{\epsilon}|^4 dx, \quad (4.17) $$

where $\Psi_n^\epsilon = c_\alpha \|\phi_n^\epsilon\|_{L^2(\Omega_H)}^2 + (\|\sqrt{\epsilon}\psi_n^\epsilon\|_{L^2(\Omega)}^2 + \|\sqrt{\epsilon}\varphi_{i,n}^\epsilon\|_{L^2(\Omega_H)}^2)$ and $U_n^\epsilon = \|u_n^\epsilon\|_{L^2(\Omega_H)}^2$.

According to regularity $f_H \in L^{2/\alpha_1}(0, T, L^2(\Omega_H))$ and by choosing $\eta = \frac{\Gamma(\beta)\mu_2/\mu_1}{T^{\beta-\alpha}\Gamma(\alpha)}$, we can get by Corollary 3.2 and the fact that quantity $\Psi_n^\epsilon(0) + U_n^\epsilon(0)$ is uniformly bounded (from (4.9) and $0 < \epsilon < 1$)

$$ \Psi_n^\epsilon(t) + U_n^\epsilon(t) \leq c_{11}(1 + \| f_H \|_{L^{2/\alpha_1}(0,T,L^2(\Omega_H))}^2), \quad (4.18) $$

for all $t \in (0, T)$ and then $\Psi_n^\epsilon(t)$ and $U_n^\epsilon(t)$ is uniformly bounded with respect to $\epsilon$ and to $n$. This ensures that

$$ \text{the sequences } (\phi_n^\epsilon, \sqrt{\epsilon}\varphi_{i,n}^\epsilon, u_n^\epsilon) \text{ and } (\sqrt{\epsilon}\psi_n^\epsilon) \text{ are bounded sets} \tag{4.19} \\ \text{of } L^\infty(0, T; L^2(\Omega_H)) \text{ and } L^\infty(0, T; L^2(\Omega)), \text{ respectively.} $$

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 $$

$$ \Gamma(\alpha)\Psi_n^\epsilon(t) + \int_0^t (t-\tau)^{\alpha-1}(2\nu_i \| \nabla\varphi_{i,n}^\epsilon \|_{L^2(\Omega_H)}^2 + 2\nu_g \| \nabla\psi_n^\epsilon \|_{L^2(\Omega)}^2 + \frac{\mu_2}{\mu_1} \| \phi_n^\epsilon \|_{L^4(\Omega_H)}^4) d\tau \quad (4.20) \\ \leq c_{11}C_0 $$

where $C_0 = (1 + \| f_H \|_{L^{2/\alpha_1}(0,T,L^2(\Omega_H))}^2)$.

We can deduce that (since $\phi_n^\epsilon = \varphi_{i,n}^\epsilon - \varphi_{e,n}^\epsilon$ and $\psi_n^\epsilon = \varphi_{e,n}^\epsilon$ in $\Omega_H$)

$$ \begin{aligned} I_{0+}^{\alpha} \left[ \| \nabla \phi_n^{\epsilon} \|_{L^2(\Omega_H)}^2 \right] (t) &= \frac{1}{\Gamma(\alpha)} \int_0^t (t-\tau)^{\alpha-1} \| \nabla \phi_n^{\epsilon} \|_{L^2(\Omega_H)}^2 d\tau \leq c_{12} C_0, \\ I_{0+}^{\alpha} \left[ \| \nabla \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2 \right] (t) &= \frac{1}{\Gamma(\alpha)} \int_0^t (t-\tau)^{\alpha-1} \| \nabla \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2 d\tau \leq c_{12} C_0, \\ I_{0+}^{\alpha} \left[ \| \nabla \psi_n^{\epsilon} \|_{L^2(\Omega)}^2 \right] (t) &= \frac{1}{\Gamma(\alpha)} \int_0^t (t-\tau)^{\alpha-1} \| \nabla \psi_n^{\epsilon} \|_{L^2(\Omega)}^2 d\tau \leq c_{12} C_0, \\ I_{0+}^{\alpha} \left[ \| \phi_n^{\epsilon} \|_{L^4(\Omega_H)}^4 \right] (t) &= \frac{1}{\Gamma(\alpha)} \int_0^t (t-\tau)^{\alpha-1} \| \phi_n^{\epsilon} \|_{L^4(\Omega_H)}^4 d\tau \leq c_{12} C_0. \end{aligned} \quad (4.21) $$

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}$)

$$ \text{the sequences } (\varphi_{i,n}^{\epsilon}) \text{ and } (\psi_{i,n}^{\epsilon}) \text{ are bounded sets} \\ \text{of } L^2(0, T; H^1(\Omega_H)) \text{ and } L^2(0, T; H^1(\Omega)), \text{ respectively} \quad (4.22) $$

and the sequence $(\phi_n^\epsilon)$ is bounded set of $L^4(0, T, L^4(\Omega_H)) \cap L^2(0, T; H^1(\Omega_H))$.

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

$$ \begin{aligned} & \left| \int_0^T \int_{\Omega_H} \mathcal{I}(\phi_n^\epsilon, u_n^\epsilon) \phi_n^\epsilon dx \right| \text{ is uniformly bounded with respect to } n, \\ & \left| \int_0^T \int_{\Omega_H} \mathcal{G}(\phi_n^\epsilon, u_n^\epsilon) u_n^\epsilon dx \right| \text{ is uniformly bounded with respect to } n. \end{aligned} \quad (4.23) $$
Now we estimate the fractional derivative $\partial_{0+}^{\alpha}(\phi_n^{\epsilon}, \varphi_{i,n}^{\epsilon}, \psi_n^{\epsilon})$ and $\partial_{0+}^{\beta} u_n^{\epsilon}$. Taking $\mathbf{h}_n = \partial_{0+}^{\alpha} \varphi_{i,n}^{\epsilon}, \mathbf{e}_n = \partial_{0+}^{\alpha} \psi_n^{\epsilon}$ and $\mathbf{v}_n = \partial_{0+}^{\beta} u_n^{\epsilon}$ and using uniform coercivity of forms $A_i$ and $\mathbf{A}_g$, we can deduce (since $\psi_n^{\epsilon} = \varphi_{e,n}^{\epsilon}$ and $\phi_n^{\epsilon} = \varphi_{i,n}^{\epsilon} - \varphi_{e,n}^{\epsilon}$ in $\Omega_H$)

$$
\begin{align}
& c_{\alpha} \| \partial_{0+}^{\alpha} \phi_n^{\epsilon} \|_{L^2(\Omega_H)}^2 + \epsilon (\| \partial_{0+}^{\alpha} \psi_n^{\epsilon} \|_{L^2(\Omega)}^2 + \| \partial_{0+}^{\alpha} \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2) \nonumber \\
& \qquad + \frac{\nu_i}{2} \partial_{0+}^{\alpha} \| \nabla \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2 + \frac{\nu_g}{2} \partial_{0+}^{\alpha} \| \nabla \psi_n^{\epsilon} \|_{L^2(\Omega)}^2 \nonumber \\
& \leq - \int_{\Omega_H} \mathcal{I}(\phi_n^{\epsilon}, u_n^{\epsilon}) \partial_{0+}^{\alpha} \phi_n^{\epsilon} d\mathbf{x} + \int_{\Omega_H} f_H \partial_{0+}^{\alpha} \phi_n^{\epsilon} d\mathbf{x}, \tag{4.24} \\
& \| \partial_{0+}^{\beta} u_n^{\epsilon} \|_{L^2(\Omega_H)}^2 = - \int_{\Omega_H} \mathcal{G}(\phi_n^{\epsilon}, u_n^{\epsilon}) \partial_{0+}^{\alpha} u_n^{\epsilon} d\mathbf{x}. \nonumber
\end{align}
$$

Then

$$
\begin{align}
& c_{\alpha} I_{0+}^{\alpha} \left[ \| \partial_{0+}^{\alpha} \phi_n^{\epsilon} \|_{L^2(\Omega_H)}^2 \right] (t) + \epsilon \left( I_{0+}^{\alpha} \left[ \| \partial_{0+}^{\alpha} \psi_n^{\epsilon} \|_{L^2(\Omega)}^2 \right] (t) + I_{0+}^{\alpha} \left[ \| \partial_{0+}^{\alpha} \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2 \right] (t) \right) \nonumber \\
& \qquad + \frac{\nu_i}{2} \| \nabla \varphi_{i,n}^{\epsilon} (t) \|_{L^2(\Omega_H)}^2 + \frac{\nu_g}{2} \| \nabla \psi_n^{\epsilon} (t) \|_{L^2(\Omega)}^2 \nonumber \\
& \leq \frac{\nu_i}{2} \| \nabla \varphi_{i,n}^{\epsilon} (0) \|_{L^2(\Omega_H)}^2 + \frac{\nu_g}{2} \| \nabla \psi_n^{\epsilon} (0) \|_{L^2(\Omega)}^2 \nonumber \\
& \qquad - \int_{\Omega_H} I_{0+}^{\alpha} \left[ \mathcal{I}(\phi_n^{\epsilon}, u_n^{\epsilon}) \partial_{0+}^{\alpha} \phi_n^{\epsilon} \right] (t) d\mathbf{x} + \int_{\Omega_H} I_{0+}^{\alpha} \left[ f_H \partial_{0+}^{\alpha} \phi_n^{\epsilon} \right] (t) d\mathbf{x}, \tag{4.25} \\
& \int_0^t \| \partial_{0+}^{\beta} u_n^{\epsilon} \|_{L^2(\Omega_H)}^2 d\tau = - \int_0^t \int_{\Omega_H} \mathcal{G}(\phi_n^{\epsilon}, u_n^{\epsilon}) \partial_{0+}^{\alpha} u_n^{\epsilon} d\mathbf{x} d\tau. \nonumber
\end{align}
$$

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

According to (H2)$_4$ and regularity of $\bar{\hbar}$, we can deduce that

$$
\begin{align}
& - \int_{\Omega_H} \mathcal{G}(\phi_n^{\epsilon}, u_n^{\epsilon}) \partial_{0+}^{\beta} u_n^{\epsilon} d\mathbf{x} \le c_0 \int_{\Omega_H} (1 + |\phi_n^{\epsilon}|^2 + |u_n^{\epsilon}|) |\partial_{0+}^{\beta} u_n^{\epsilon}| d\mathbf{x} \tag{4.26} \\
& \le \frac{1}{2} \| \partial_{0+}^{\beta} u_n^{\epsilon} \|_{L^2(\Omega_H)}^2 + c_1 (1 + \| \phi_n^{\epsilon} \|_{L^4(\Omega_H)}^4 + \| u_n^{\epsilon} \|_{L^2(\Omega_H)}^2) \nonumber
\end{align}
$$

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

$$
\| \partial_{0+}^{\beta} u_n^{\epsilon} \|_{L^2(0,t,L^2(\Omega_H))}^2 \le c_2(1 + \| \phi_n^{\epsilon} \|_{L^4(\Omega_H)}^4 + \| u_n^{\epsilon} \|_{L^{\infty}(0,T,L^2(\Omega_H))}^2) \le C. \quad (4.27)
$$

Before to estimate the term $\int_{\Omega_H} I_{0+}^{\alpha} [\mathcal{I}(\phi_n^{\epsilon}, u_n^{\epsilon}) \partial_{0+}^{\alpha} \phi_n^{\epsilon}] (t) d\mathbf{x}$, we prove that $u_n^{\epsilon}$ is uniformly bounded in $L^{\infty}(0, T, L^3(\Omega_H))$. Since $u_n^{\epsilon}$ satisfies

$$
u_n^{\epsilon}(\mathbf{x}, t) = u_n^{\epsilon}(\mathbf{x}, 0) - \frac{1}{\Gamma(\beta)} \int_0^t (t - \tau)^{\beta-1} \mathcal{G}(\phi_n^{\epsilon}, u_n^{\epsilon})(\mathbf{x}, \tau) d\tau,
$$

then, according to regularity of $\bar{\hbar}$, we have

$$
| u_n^{\epsilon}(\mathbf{x}, t) | \le c_3 \left( | u_n^{\epsilon}(\mathbf{x}, 0) | + \int_0^t (t-\tau)^{\beta-1} | \mathcal{I}_2(\phi_n^{\epsilon}) | d\tau + \int_0^t (t-\tau)^{\beta-1} | u_n^{\epsilon}(\mathbf{x}, \tau) | d\tau \right). \quad (4.28)
$$
Consequently (since from assumption (2.4) we have $|\mathcal{I}_2(\phi_n^\epsilon)| \le \beta_5 + \beta_6 |\phi_n^\epsilon|^2$)

$$|u_n^\epsilon(\mathbf{x}, t)| \le c_4 \left( |u_n^\epsilon(\mathbf{x}, 0)| + \int_0^t (t-\tau)^{\beta-1} d\tau \right) + c_5 \left( \int_0^t (t-\tau)^{\beta-1} |\phi_n^\epsilon|^2 d\tau + \int_0^t (t-\tau)^{\beta-1} |u_n^\epsilon(\mathbf{x}, \tau)| d\tau \right).$$

This implies (since $\int_0^t (t-\tau)^{\beta-1} d\tau \le T^\beta/\beta$)

$$\left( \int_{\Omega_H} |u_n^\epsilon(\mathbf{x}, t)|^3 d\mathbf{x} \right)^{1/3} \le c_6 \left( 1 + \| u_n^\epsilon(\mathbf{x}, 0) \|_{L^3(\Omega_H)} + \left[ \int_{\Omega_H} \left( \int_0^t (t-\tau)^{\beta-1} |u_n^\epsilon(\mathbf{x}, \tau)| d\tau \right)^3 d\mathbf{x} \right]^{1/3} \right. \\ \left. + \left[ \int_{\Omega_H} \left( \int_0^t (t-\tau)^{\beta-1} |\phi_n^\epsilon(\mathbf{x}, s)|^2 ds \right)^3 d\mathbf{x} \right]^{1/3} \right).$$

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

$$\| u_n^\epsilon(., t) \|_{L^3(\Omega_H)} \le c_7 \left( 1 + \| u_n^\epsilon(., 0) \|_{L^3(\Omega)} + \int_0^t (t-\tau)^{\beta-1} \| u_n^\epsilon(., \tau) \|_{L^3(\Omega_H)} d\tau \right. \\ \left. + \int_0^t (t-\tau)^{\beta-1} \| \phi_n^\epsilon(., \tau) \|_{H^1(\Omega_H)}^2 d\tau \right).$$

Using (4.19), (4.21) and the fact, from (4.9) the quantities $u_n^\epsilon(0)$ is uniformly bounded in $L^3$, we obtain

$$\| u_n^\epsilon(., t) \|_{L^3(\Omega_H)} \le c_8 + c_9 I_{0+}^\beta \left[ \| u_n^\epsilon \|_{L^3(\Omega_H)} \right] (t) \qquad (4.29)$$

Consequently

$$\text{the sequence } u_n^\epsilon \text{ is uniformly bounded in } L^\infty(0, T, L^3(\Omega_H)). \qquad (4.30)$$

We can now estimate the term $\int_{\Omega_H} I_{0+}^\alpha [\mathcal{I}(\phi_n^\epsilon, u_n^\epsilon) \partial_{0+}^\alpha \phi_n^\epsilon](t) d\mathbf{x}$.

Since $\mathcal{I}(\phi_n^\epsilon, u_n^\epsilon) = \mathcal{I}_0(\phi_n^\epsilon) + \mathcal{I}_1(\phi_n^\epsilon)u_n^\epsilon$ and $\mathcal{I}_0(\phi_n^\epsilon) = g_1(\phi_n^\epsilon) + g_2(\phi_n^\epsilon)$ with $g_1$ is an increasing function, we have first (from $(\mathbf{H}2)_4$)

$$\begin{aligned} & - \int_{\Omega_H} (g_2(\phi_n^\epsilon) + \mathcal{I}_1(\phi_n^\epsilon)u_n^\epsilon) \partial_{0+}^\alpha \phi_n^\epsilon d\mathbf{x} \le c_0 \int_{\Omega_H} (1 + |\phi_n^\epsilon|^2 + |u_n^\epsilon| \| \phi_n^\epsilon \| + |u_n^\epsilon|) \partial_{0+}^\alpha \phi_n^\epsilon d\mathbf{x} \\ & \le \frac{c_1}{4} \| \partial_{0+}^\alpha \phi_n^\epsilon \|_{L^2(\Omega_H)}^2 + c_2(1 + \| \phi_n^\epsilon \|_{L^4(\Omega_H)}^4 + \| u_n^\epsilon \|_{L^2(\Omega_H)}^2 + \| \phi_n^\epsilon \|_{L^6(\Omega_H)}^2) \| u_n^\epsilon \|_{L^3(\Omega_H)}^2. \end{aligned} \qquad (4.31)$$

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_{0+}^\alpha E_g(\phi_n^\epsilon) \le g_1(\phi_n^\epsilon) \partial_{0+}^\alpha \phi_n^\epsilon$ and then

$$-I_{0+}^{\alpha} [g_1(\phi_n^\epsilon) \partial_{0+}^{\alpha} \phi_n^\epsilon](t) \le -I_{0+}^{\alpha} [\partial_{0+}^{\alpha} E(\phi_n^\epsilon)](t). \qquad (4.32)$$
Now, we estimate the term $-I_{0+}^{\alpha} [\partial_{0+}^{\alpha} E(\phi_n^{\epsilon})]$.

$$ - \int_{\Omega_H} I_{0+}^{\alpha} [\partial_{0+}^{\alpha} E_g(\phi_n^{\epsilon})](t) d\mathbf{x} = - \int_{\Omega_H} E_g(\phi_n^{\epsilon})(t) d\mathbf{x} + \int_{\Omega_H} E_g(\phi_n^{\epsilon})(0) d\mathbf{x}. \quad (4.33) $$

Since, from (H2)$_4$, $\mu_6 | \phi_n^{\epsilon}|^2 - \mu_7 \le E_g(\phi_n^{\epsilon}) \le c_2(1 + |\phi_n^{\epsilon}|^4)$ then

$$ - \int_{\Omega_H} I_{0+}^{\alpha} [\partial_{0+}^{\alpha} E(\phi_n^{\epsilon})] d\mathbf{x} \le c_3(1 + \|\phi_n^{\epsilon}(0)\|_{H^1(\Omega_H)}^4 + \|\phi_n^{\epsilon}(t)\|_{L^2(\Omega_H)}^2). \quad (4.34) $$

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

$$ - \int_{\Omega_H} I_{0+}^{\alpha} [\mathcal{I}(\phi_n^{\epsilon}, u_n^{\epsilon}) \partial_{0+}^{\alpha} \phi_n^{\epsilon}] (t) d\mathbf{x} \le \frac{c_{\alpha}}{4} I_{0+}^{\alpha} \left[ \| \partial_{0+}^{\alpha} \phi_n^{\epsilon} \|_{L^2(\Omega_H)}^2 \right. \\ \left. \qquad + c_4 (\| u_n^{\epsilon} \|_{L^{\infty}(0,T,L^2(\Omega_H))}^2 + \| \phi_n^{\epsilon} \|_{L^{\infty}(0,T,L^2(\Omega_H))}^2) \right. \\ \left. \qquad + c_5 (1 + I_{0+}^{\alpha} [\| \phi_n^{\epsilon} \|_{L^4(\Omega_H)}^4] + \| \phi_n^{\epsilon}(0) \|_{H^1(\Omega_H)}^4), \right. \\ \left. \qquad + c_6 I_{0+}^{\alpha} [\| \phi_n^{\epsilon} \|_{H^1(\Omega_H)}^2] \| u_n^{\epsilon} \|_{L^{\infty}(0,T,L^3(\Omega_H))}^2. \right] \quad (4.35) $$

Finally for the external forcing we have (using Hölder inequality)

$$ I_{0+}^{\alpha} \left[ \int_{\Omega_H} f_H \partial_{0+}^{\alpha} \psi_n^{\epsilon} d\mathbf{x} \right] \le \frac{c_{\alpha}}{4} I_{0+}^{\alpha} \left[ \| \partial_{0+}^{\alpha} \psi_n^{\epsilon} \|_{L^2(\Omega_H)}^2 \right] + c_7 \| f_H \|_{L^{2/\alpha_1}(0,T,L^2(\Omega_H))}^2. \quad (4.36) $$

According to (4.35), (4.36), (4.19), (4.21), and the fact, from (4.9) and $\epsilon \in (0, 1)$, the quantities $\psi_n^{\epsilon}(0)$, $\phi_n^{\epsilon}(0)$ and $\varphi_{i,n}^{\epsilon}(0)$ are uniformly bounded in $H^1$, we can derive from (4.25) the following estimate (for all $t \in (0, T)$)

$$ \begin{aligned} & \frac{c_{\alpha}}{2} I_{0+}^{\alpha} [\| \partial_{0+}^{\alpha} \phi_n^{\epsilon} \|_{L^2(\Omega_H)}^2](t) + \epsilon (I_{0+}^{\alpha} [\| \partial_{0+}^{\alpha} \psi_n^{\epsilon} \|_{L^2(\Omega)}^2](t) + I_{0+}^{\alpha} [\| \partial_{0+}^{\alpha} \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2](t)) \\ & + \frac{v_i}{2} \| \nabla \varphi_{i,n}^{\epsilon}(t) \|_{L^2(\Omega_H)}^2 + \frac{v_g}{2} \| \nabla \psi_n^{\epsilon}(t) \|_{L^2(\Omega)}^2 \\ & \le c_8 (1 + \| \varphi_{i,n}^{\epsilon}(0) \|_{H^1(\Omega_H)}^2 + \| \psi_n^{\epsilon}(0) \|_{H^1(\Omega)}^2 + \| \phi_n^{\epsilon}(0) \|_{H^1(\Omega_H)}^4) \\ & + c_9 (\| f_H \|_{L^{2/\alpha_1}(0,T,L^2(\Omega_H))}^2) \\ & + c_{10} (\| u_n^{\epsilon} \|_{L^{\infty}(0,T,L^2(\Omega_H))}^2 + \| \phi_n^{\epsilon} \|_{L^{\infty}(0,T,L^2(\Omega_H))}^2 + I_{0+}^{\alpha} \| \phi_n^{\epsilon} \|_{L^4(\Omega_H)}^4) \\ & + c_{11} I_{0+}^{\alpha} [\| \phi_n^{\epsilon} \|_{H^1(\Omega_H)}^2] \| u_n^{\epsilon} \|_{L^{\infty}(0,T,L^3(\Omega_H))}^2 \le C. \end{aligned} \quad (4.37) $$

We can deduce that (from (4.27) and (4.37))

$$ \begin{aligned} & I_{0+}^{\alpha} [\| \partial_{0+}^{\alpha} \phi_n^{\epsilon} \|_{L^2(\Omega_H)}^2](t) + I_{0+}^{\alpha} [\| \sqrt{\epsilon} \partial_{0+}^{\alpha} \psi_n^{\epsilon} \|_{L^2(\Omega)}^2](t) + I_{0+}^{\alpha} [\| \sqrt{\epsilon} \partial_{0+}^{\alpha} \varphi_{i,n}^{\epsilon} \|_{L^2(\Omega_H)}^2](t) \\ & + \| \nabla \varphi_{i,n}^{\epsilon}(t) \|_{L^2(\Omega_H)}^2 + \| \nabla \psi_n^{\epsilon}(t) \|_{L^2(\Omega)}^2 \le C, \end{aligned} \quad (4.38) \\ \| \partial_{0+}^{\beta} u_n^{\epsilon} \|_{L^2(0,t,L^2(\Omega_H))}^2 \le C. $$

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_{T_k}^\gamma$ (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 argument and, it follows that there exist $(X^\epsilon; u^\epsilon) = (\phi^\epsilon, \phi_i^\epsilon, \psi^\epsilon; u^\epsilon)$ and $(\tilde{X}^\epsilon; \tilde{u}^\epsilon) = (\tilde{\phi}^\epsilon, \tilde{\phi}_i^\epsilon, \tilde{\psi}^\epsilon; \tilde{u}^\epsilon)$ such that there exists a subsequence of $(X_n^\epsilon; u_n^\epsilon) = (\phi_n^\epsilon, \phi_{i,n}^\epsilon, \psi_n^\epsilon; u_n^\epsilon)$ also denoted by $(X_n^\epsilon; u_n^\epsilon)$, such that

$$ 
\begin{align} 
(X_n^\epsilon; u_n^\epsilon) &\to (X^\epsilon; u^\epsilon) \quad \text{weakly in } L^\infty(0, T; \mathbb{V}_H \times \mathbb{V}_H \times \mathbb{V}) \times L^\infty(0, T; L^3(\Omega_H)), \nonumber \\
(X_n^\epsilon; u_n^\epsilon) &\to (X^\epsilon; u^\epsilon) \quad \text{strongly in } (L^2(\mathcal{Q}_H)^2 \times L^2(\mathcal{Q}) \times L^2(\mathcal{Q}_H), \nonumber \\
(\partial_{0+}^\alpha X_n^\epsilon; \partial_{0+}^\beta u_n^\epsilon) &\to (\partial_{0+}^\alpha \tilde{X}^\epsilon; \partial_{0+}^\beta \tilde{u}^\epsilon) \quad \text{weakly in } (L^2(\mathcal{Q}_H)^2 \times L^2(\mathcal{Q}) \times L^2(\mathcal{Q}_H)). \tag{4.39}
\end{align} 
$$

First we show that $(\partial_{0+}^\alpha X_n^\epsilon, \partial_{0+}^\beta u_n^\epsilon)$ exists in the weak sense and that $(\partial_{0+}^\alpha X^\epsilon; \partial_{0+}^\beta u^\epsilon) = (\tilde{X}^\epsilon; \tilde{u}^\epsilon)$. Indeed, we take $\omega \in C_0^\infty(0, T)$ and $\nu \in \mathbb{V}_H$ (then $\omega\nu \in D_4(0, T)$), by the weak convergence and Lebesgue's dominated convergence arguments we have

$$ 
\begin{align} 
\langle \tilde{\phi}^\epsilon, \omega\nu \rangle_{D'_4, D_4} &= \lim_{n \to \infty} \int_0^T \int_{\Omega_H} \omega(t) \partial_{0+}^\alpha \phi_n^\epsilon(\mathbf{x}, t) v(\mathbf{x}) d\mathbf{x} dt \nonumber \\
&= \lim_{n \to \infty} \int_{\Omega} \left( \int_0^T \omega(t) \partial_{0+}^\alpha \phi_n^\epsilon dt \right) v(\mathbf{x}) d\mathbf{x} \nonumber \\
&= - \lim_{n \to \infty} \int_0^T D_{T-}^\alpha \omega(t) \left( \int_{\Omega} \phi_n^\epsilon(\mathbf{x}, t) v(\mathbf{x}) d\mathbf{x} \right) dt - I_{0+}^{1-\alpha} \omega(0^+) \int_{\Omega} \phi_n^\epsilon(\mathbf{x}, 0) v(\mathbf{x}) d\mathbf{x} \tag{4.40} \\
&= - \int_0^T D_{T-}^\alpha \omega(t) \left( \int_{\Omega} \phi^\epsilon(\mathbf{x}, t) v(\mathbf{x}) d\mathbf{x} \right) dt - I_{0+}^{1-\alpha} \omega(0^+) \int_{\Omega} \phi^\epsilon(\mathbf{x}, 0) v(\mathbf{x}) d\mathbf{x} \nonumber \\
&= \int_0^T \int_{\Omega_H} \omega(t) \partial_{0+}^\alpha \phi^\epsilon(\mathbf{x}, t) v(\mathbf{x}) d\mathbf{x} dt. \nonumber
\end{align} 
$$

Consequently, $\tilde{\phi}^\epsilon = \partial_{0+}^\alpha \phi^\epsilon$ in the weak sense. In the same way we prove, in the weak sense, that $(\partial_{0+}^\alpha \phi_i^\epsilon, \partial_{0+}^\alpha \psi^\epsilon, \partial_{0+}^\beta u^\epsilon) = (\tilde{\phi}_i^\epsilon, \tilde{\psi}^\epsilon, \tilde{u}^\epsilon)$.

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

$$ 
\begin{align} 
\zeta_\alpha \int_0^T \kappa(t) \langle \partial_{0+}^\alpha \phi_n^\epsilon, \mathbf{h}_n \rangle_{\mathbb{V}'_H, \mathbb{V}_H} dt + \epsilon \int_0^T \kappa(t) \langle \partial_{0+}^\alpha \phi_{i,n}^\epsilon, \mathbf{h}_n \rangle_{\mathbb{V}'_H, \mathbb{V}_H} dt \nonumber \\
= - \int_0^T \kappa(t) \int_{\Omega_H} \mathcal{I}(\phi_n^\epsilon, u_n^\epsilon) \mathbf{h}_n d\mathbf{x} dt - \int_0^T \kappa(t) A_i(\phi_{i,n}^\epsilon, \mathbf{h}_n) dt + \int_0^T \kappa(t) \int_{\Omega_H} f_H \mathbf{h}_n d\mathbf{x} dt, \nonumber \\
\zeta_\alpha \int_0^T \kappa(t) \langle \partial_{0+}^\alpha \phi_n^\epsilon, \mathbf{k}_n \rangle_{\mathbb{V}'_H, \mathbb{V}_H} dt - \epsilon \int_0^T \kappa(t) \langle \partial_{0+}^\alpha \psi_n^\epsilon, \mathbf{k}_n \rangle_{\mathbb{V}'_H, \mathbb{V}_H} dt \nonumber \\
= - \int_0^T \kappa(t) \int_{\Omega_H} \mathcal{I}(\phi_n^\epsilon, u_n^\epsilon) \mathbf{k}_n d\mathbf{x} dt + \int_0^T \kappa(t) A_g(\psi_n^\epsilon, \mathbf{k}_n) dt + \int_0^T \kappa(t) \int_{\Omega_H} f_H \mathbf{k}_n d\mathbf{x} dt, \nonumber \\
\int_0^T \kappa(t) \langle \partial_{0+}^\beta u_n^\epsilon, \rho_n \rangle_{\mathbb{V}'_H, \mathbb{V}_H} dt = - \int_0^T \kappa(t) \int_{\Omega_H} \mathcal{G}(\phi_n^\epsilon, u_n^\epsilon) \rho_n d\mathbf{x} dt. \nonumber
\end{align} 
$$

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 (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 $(h, k, \rho) \in \mathbb{V}_H \times \mathbb{V} \times \mathbb{V}_H$)

$$ 
\begin{aligned}
& \langle c_\alpha \partial_{0^+}^\alpha \phi^\epsilon, h \rangle_{\mathbb{V}'_H, \mathbb{V}_H} + \epsilon \langle \partial_{0^+}^\alpha \varphi_i^\epsilon, h \rangle_{\mathbb{V}'_H, \mathbb{V}_H} + \int_{\Omega_H} \mathcal{I}(\phi^\epsilon, u^\epsilon)hdx + A_i(\varphi_i^\epsilon, h) = \int_{\Omega_H} f_H hdx, \\
& \langle c_\alpha \partial_{0^+}^\alpha \phi^\epsilon, k \rangle_{\mathbb{V}'_H, \mathbb{V}_H} - \epsilon \langle \partial_{0^+}^\alpha \psi^\epsilon, k \rangle_{\mathbb{V}'_H, \mathbb{V}_H} + \int_{\Omega_H} \mathcal{I}(\phi^\epsilon, u^\epsilon)kdx - A_g(\psi, v) = \int_{\Omega_H} f_H kdx, \\
& (\partial_{0^+}^\beta u^\epsilon, \rho)_{L^2(\Omega_H)} + \int_{\Omega_H} \mathcal{G}(\phi^\epsilon, u^\epsilon)\rho dx = 0.
\end{aligned}
\quad (4.41) $$

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.

<figure><img src="image_6.png" /></figure>

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 (H2)$_4$ hold, $f_H \in L^{2/\alpha_1}(0, T; L^2(\Omega))$ and $u_0 \in L^3(\Omega_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)$)

$$ 
\begin{aligned}
& \| \phi^\epsilon(t) \|_{L^2(\Omega_H)}^2 + \| \sqrt{\epsilon} \psi^\epsilon(t) \|_{L^2(\Omega)}^2 + \| \sqrt{\epsilon} \varphi_i^\epsilon(t) \|_{L^2(\Omega_H)}^2 \\
& \qquad + I_{0^+}^\alpha [ \| \phi^\epsilon \|_{L^4(\Omega_H)}^4 ](t) + I_{0^+}^\alpha [ \| \nabla \varphi_i^\epsilon \|_{L^2(\Omega_H)}^2 ](t) + I_{0^+}^\alpha [ \| \nabla \psi_n^\epsilon \|_{L^2(\Omega)}^2 ](t) \leq C, \\
& \| u^\epsilon \|_{L^\infty(0,t,L^3(\Omega_H))}^2 + \| \partial_{0^+}^\beta u^\epsilon \|_{L^2(0,t,L^2(\Omega_H))}^2 \leq C.
\end{aligned}
\quad (4.42) $$

(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)$)

$$ 
\begin{aligned}
& I_{0^+}^\alpha [ \| \partial_{0^+}^\alpha \phi^\epsilon \|_{L^2(\Omega_H)}^2 ](t) + I_{0^+}^\alpha [ \| \sqrt{\epsilon} \partial_{0^+}^\alpha \psi^\epsilon \|_{L^2(\Omega)}^2 ](t) + I_{0^+}^\alpha [ \| \sqrt{\epsilon} \partial_{0^+}^\alpha \varphi_i^\epsilon \|_{L^2(\Omega_H)}^2 ](t) \\
& \qquad + \| \nabla \varphi_i^\epsilon(t) \|_{L^2(\Omega_H)}^2 + \| \nabla \psi^\epsilon(t) \|_{L^2(\Omega)}^2 \leq C, \\
& \| u^\epsilon \|_{L^\infty(0,t,L^3(\Omega_H))}^2 + \| \partial_{0^+}^\beta u^\epsilon \|_{L^2(0,t,L^2(\Omega_H))}^2 \leq C.
\end{aligned}
\quad (4.43) $$

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

**Theorem 4.2.** Let assumptions (H1) and (H2)$_4$ be fulfilled and $\beta \ge \alpha$. Then for $(\phi_0, u_0)$ and $f_H$ given such that $(\phi_0, u_0) \in \mathbb{V}_H \times L^3(\Omega_H)$ and $f_H \in L^{\frac{2}{\alpha_1}}(0, T; L^2(\Omega_H))$ with $0 < \alpha_1 < \alpha$, there exists a weak solution $(\varphi_i, \phi, \psi, u)$ to problem (1.1) verifying

$$ 
\begin{aligned}
& \phi \in L^\infty(0, T; \mathbb{V}_H), \quad \partial_{0^+}^\alpha \phi \in L^2(0, T; L^2(\Omega_H)), \\
& \psi \in L^\infty(0, T; \mathbb{V}), \quad \psi \in L^\infty(0, T; \mathbb{V}_H), \\
& u \in L^\infty(0, T; L^3(\Omega_H)), \quad \partial_{0^+}^\beta u \in L^2(0, T; L^2(\Omega_H)).
\end{aligned}
\quad (4.44) $$
Moreover, if $\beta > 1/2$ then $u \in C([0, T]; L^2(\Omega_H))$ and $u(0^+) = u_0$ and if $\alpha > 1/2$ then $\phi \in C([0, T]; L^2(\Omega_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

$$ 
\begin{aligned}
(X^\epsilon; u^\epsilon) &\to (X; u) \text{ weakly in } L^\infty(0, T; \mathbb{V}_H \times \mathbb{V}_H \times \mathbb{V}) \times L^\infty(0, T; L^3(\Omega_H)), \\
(X^\epsilon; u^\epsilon) &\to (X; u) \text{ strongly in } (L^2(\mathbb{Q}_H))^2 \times L^2(\mathbb{Q}) \times L^2(\mathbb{Q}_H), \\
(\partial_{0+}^\alpha X^\epsilon; \partial_{0+}^\beta u^\epsilon) &\to (\partial_{0+}^\alpha \phi, 0, 0; \partial_{0+}^\beta u) \text{ weakly in } (L^2(\mathbb{Q}_H))^2 \times L^2(\mathbb{Q}) \times L^2(\mathbb{Q}_H).
\end{aligned}
\quad (4.45) $$

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. $\square$

**Theorem 4.3.** Let assumptions (H1), (H2)$_4$ and (H3) be fulfilled and $\beta \ge \alpha$. Then the solution $(\varphi_i, \phi, \psi, u)$ to system (1.1) is unique. Moreover, let $(\phi_0^{(i)}, u_0^{(i)}, f_H^{(i)})$ be given such that $(\phi_0^{(i)}, u_0^{(i)}) \in \mathbb{V}_H \times L^3(\Omega_H)$ and $f_H^{(i)} \in L^{\frac{2}{\alpha_1}}(0, T; L^2(\Omega_H))$, for $i = 1, 2$, with $0 < \alpha_1 < \alpha$. For $(\varphi_i^{(j)}, \phi_i^{(j)}, \psi_i^{(j)}, u_i^{(j)})$ a weak solutions to system (1.1), which corresponds to data $(\phi_0^{(j)}, u_0^{(j)}, f_H^{(j)})$ (for $j = 1, 2$), the following Lipschitz continuity relationship is satisfied.

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

$$ 
\begin{aligned}
& \| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2 + I_{0+}^\alpha \left[ \| \nabla \varphi_i \|_{L^2(\Omega_H)}^2 + \| \nabla \psi \|_{L^2(\Omega_H)}^2 \right] (t) \\
& \le C (\| \phi_0 \|_{L^2(\Omega_H)}^2 + \| u_0 \|_{L^2(\Omega_H)}^2 + \| f_H \|_{L^{2/\alpha_1}(0,T; L^2(\Omega_H))}^2).
\end{aligned}
\quad (4.46) $$

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

$$ | \mathcal{I}_2(\phi_1) - \mathcal{I}_2(\phi_2) | \le C_I | \phi | (1 + | \phi_1 | + | \phi_2 |), \quad (4.47) $$

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_i^{(j)}, \phi_i^{(j)}, \psi_i^{(j)}, u_i^{(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

$$ 
\begin{aligned}
& c_\alpha(\partial_{0+}^\alpha \phi, \varphi_i)_{L^2(\Omega_H)} + \int_{\Omega_H} (\mathcal{I}(\phi_1, u_1) - \mathcal{I}(\phi_2, u_2))\varphi_i d\mathbf{x} + A_i(\varphi_i, \varphi_i) = \int_{\Omega_H} f_H \varphi_i d\mathbf{x}, \\
& c_\alpha(\partial_{0+}^\alpha \phi, \psi)_{L^2(\Omega_H)} + \int_{\Omega_H} (\mathcal{I}(\phi_1, u_1) - \mathcal{I}(\phi_2, u_2))\psi d\mathbf{x} - A_g(\psi, \psi) = \int_{\Omega_H} f_H \psi d\mathbf{x}, \\
& (\partial_{0+}^\beta u, u)_{L^2(\Omega_H)} + \int_{\Omega_H} (\mathcal{G}(\phi_1, u_1) - \mathcal{G}(\phi_2, u_2))ud\mathbf{x} = 0, \\
& u(0) = u_0, \quad \phi(0) = \phi_0,
\end{aligned}
\quad (4.48) $$

where $(\phi_0, u_0, f_H) = (\phi_0^{(1)} - \phi_0^{(2)}, u_0^{(1)} - u_0^{(2)}, f_H^{(1)} - f_H^{(2)})$.
Then (since $\psi = \varphi_e$ and $\phi = \varphi_i - \varphi_e$ in $\Omega_H$),

$$
\begin{align}
& c_{\alpha}(\partial_{0^+}^{\alpha} \phi, \phi)_{L^2(\Omega_H)} + \int_{\Omega_H} (\mathcal{I}(\phi_1, u_1) - \mathcal{I}(\phi_2, u_2))\phi d\mathbf{x} + A_i(\varphi_i, \varphi_i) + \mathbf{A}_g(\psi, \psi) = \int_{\Omega_H} f_H \phi d\mathbf{x}, \nonumber \\
& (\partial_{0^+}^{\beta} u, u)_{L^2(\Omega_H)} + \int_{\Omega_H} (\mathcal{G}(\phi_1, u_1) - \mathcal{G}(\phi_2, u_2))ud\mathbf{x} = 0, \tag{4.49} \\
& u(0) = u_0, \quad \phi(0) = \phi_0. \nonumber
\end{align}
$$

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

$$
\begin{align}
& \frac{1}{2} \min(\mu c_{\alpha}, 1) \partial_{0^+}^{\alpha} (\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2) \nonumber \\
& \quad + \mu \min(v_i, v_g) (\| \nabla \varphi_i \|_{L^2(\Omega_H)}^2 + \| \nabla \psi \|_{L^2(\Omega_H)}^2) \tag{4.50} \\
& \le - \int_{\Omega_H} (F_{\mu}(\phi_1, u_1) - F_{\mu}(\phi_2, u_2)).(\phi, u) d\mathbf{x} + \int_{\Omega_H} f_H \phi d\mathbf{x} \nonumber \\
& \le c_0 (\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2 + \| f_H(t) \|_{L^2(\Omega_H)}^2) \nonumber
\end{align}
$$

and then

$$
\begin{align}
& \partial_{0^+}^{\alpha} (\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2) + (\| \nabla \varphi_i \|_{L^2(\Omega_H)}^2 + \| \nabla \psi \|_{L^2(\Omega_H)}^2) \nonumber \\
& \le c_1 (\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2) + c_2 \| f_H(t) \|_{L^2(\Omega_H)}^2. \tag{4.51}
\end{align}
$$

So (using Hölder inequality),

$$
\begin{align}
& \| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2 + I_{0^+}^{\alpha} \left[ \| \nabla \varphi_i \|_{L^2(\Omega_H)}^2 + \| \nabla \psi \|_{L^2(\Omega_H)}^2 \right] (t) \nonumber \\
& \le \| \phi_0 \|_{L^2(\Omega_H)}^2 + \| u_0 \|_{L^2(\Omega_H)}^2 + c_3 \| f_H \|_{L^{2/\alpha_1}(0,T;L^2(\Omega_H))} \, t^{\alpha-\alpha_1} \tag{4.52} \\
& \quad + c_1 I_{0^+}^{\alpha} \left[ \| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2 \right] (t). \nonumber
\end{align}
$$

By Lemma 3.5, we can deduce that first

$$
\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2 \le c_4 (\| \phi_0 \|_{L^2(\Omega_H)}^2 + \| u_0 \|_{L^2(\Omega_H)}^2 + \| f_H \|_{L^{2/\alpha_1}(0,T;L^2(\Omega_H))}^2)
$$

and then, from (4.51)

$$
\begin{align}
& \| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2 + I_{0^+}^{\alpha} \left[ \| \nabla \varphi_i \|_{L^2(\Omega_H)}^2 + \| \nabla \psi \|_{L^2(\Omega_H)}^2 \right] (t) \tag{4.53} \\
& \le C (\| \phi_0 \|_{L^2(\Omega_H)}^2 + \| u_0 \|_{L^2(\Omega_H)}^2 + \| f_H \|_{L^{2/\alpha_1}(0,T;L^2(\Omega_H))}^2). \nonumber
\end{align}
$$

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

$$
\begin{align}
& \partial_{0^+}^{\alpha} \| \phi(t) \|_{L^2(\Omega_H)}^2 + \frac{2}{c_\alpha} \min(v_i, v_g) (\| \nabla \varphi_i \|_{L^2(\Omega_H)}^2 + \| \nabla \psi \|_{L^2(\Omega_H)}^2) \nonumber \\
& \le -\frac{2}{c_\alpha} \frac{1}{\mu} \int_{\Omega_H} F_\mu(\phi_1, u_1) - F_\mu(\phi_2, u_2) \cdot (\phi, u) d\mathbf{x} \nonumber \\
& \quad + \frac{2}{c_\alpha} \frac{1}{\mu} \int_{\Omega_H} (\mathcal{G}(\phi_1, u_1) - \mathcal{G}(\phi_2, u_2)) u d\mathbf{x} + \frac{2}{c_\alpha} \| f_H(t) \|_{L^2(\Omega_H)} \| \phi(t) \|_{L^2(\Omega_H)} \tag{4.54} \\
& \le c_5 (\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2) + c_6 \| f_H(t) \|_{L^2(\Omega_H)}^2 \nonumber \\
& \quad + \frac{2}{c_\alpha} \frac{1}{\mu} \int_{\Omega_H} (\mathcal{G}(\phi_1, u_1) - \mathcal{G}(\phi_2, u_2)) \phi d\mathbf{x}, \nonumber \\
& \partial_{0^+}^{\beta} \| u(t) \|_{L^2(\Omega_H)}^2 \le -2 \int_{\Omega_H} (\mathcal{G}(\phi_1, u_1) - \mathcal{G}(\phi_2, u_2)) u d\mathbf{x}. \nonumber
\end{align}
$$
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_H))$, for $i = 1, 2$, we have

$$
\begin{aligned}
& \left| \int_{\Omega_H} (\mathcal{G}(\phi_1, u_1) - \mathcal{G}(\phi_2, u_2)) u d\mathbf{x} \right| \le c_7 (\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2) \\
& \qquad + c_8 \int_{\Omega_H} |\phi| \| u \| (|\phi_1| + |\phi_2|) d\mathbf{x} \\
& \le c_7 (\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2) \\
& \qquad + c_9 \| u(t) \|_{L^2(\Omega_H)} \| \phi(t) \|_{L^4(\Omega_H)}.
\end{aligned}
$$

So

$$
\begin{aligned}
& \frac{2}{c_\alpha} \frac{1}{\mu} \int_{\Omega_H} (\mathcal{G}(\phi_1, u_1) - \mathcal{G}(\phi_2, u_2)) u d\mathbf{x} \le c_{10} (\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2) \\
& \qquad + \eta_1 (\| \nabla \varphi_i \|_{L^2(\Omega_H)}^2 + \| \nabla \psi \|_{L^2(\Omega_H)}^2), \tag{4.55} \\
& -2 \int_{\Omega_H} (\mathcal{G}(\phi_1, u_1) - \mathcal{G}(\phi_2, u_2)) u d\mathbf{x} \le c_{11} (\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2) \\
& \qquad + \eta_2 (\| \nabla \varphi_i \|_{L^2(\Omega_H)}^2 + \| \nabla \psi \|_{L^2(\Omega_H)}^2),
\end{aligned}
$$

with the constants $\eta_1$ and $\eta_2$ chosen appropriately. According to estimates (4.55) and by choosing $\eta_1 = \frac{1}{c_\alpha} \min(v_i, v_g)$, estimate (4.54) becomes

$$
\begin{aligned}
& \partial_{0^+}^{\alpha} \| \phi(t) \|_{L^2(\Omega_H)}^2 + \zeta \Psi(t) \le c_{12} (\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2) + c_{13} \| f_H(t) \|_{L^2(\Omega_H)}^2, \\
& \partial_{0^+}^{\beta} \| u(t) \|_{L^2(\Omega_H)}^2 \le c_{14} (\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2) + \eta_2 \Psi(t),
\end{aligned}
\quad (4.56)
$$

with $\Psi(t) = \| \nabla \varphi_i \|_{L^2(\Omega_H)}^2 + \| \nabla \psi \|_{L^2(\Omega)}^2$ and $\zeta = \frac{1}{c_\alpha} \min(v_i, v_g)$.

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

$$
\| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2 \le c_{15} (\| \phi_0 \|_{L^2(\Omega_H)}^2 + \| u_0 \|_{L^2(\Omega_H)}^2 + \| f_H \|_{L^{2/\alpha_1}(0,T;L^2(\Omega_H))}^2)
$$

and then, from (4.56)

$$
\begin{aligned}
& \| \phi(t) \|_{L^2(\Omega_H)}^2 + \| u(t) \|_{L^2(\Omega_H)}^2 + I_{0^+}^{\alpha} \left[ \| \nabla \varphi_i \|_{L^2(\Omega_H)}^2 + \| \nabla \psi \|_{L^2(\Omega_H)}^2 \right] (t) \\
& \le C (\| \phi_0 \|_{L^2(\Omega_H)}^2 + \| u_0 \|_{L^2(\Omega_H)}^2 + \| f_H \|_{L^{2/\alpha_1}(0,T;L^2(\Omega_H))}^2).
\end{aligned}
\quad (4.57)
$$

Finally, if the difference between data is null i.e., $\phi_0 = 0$, $u_0 = 0$ and $f_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. $\square$

**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)

$$
\begin{equation}
\begin{aligned}
\mathcal{I}(\phi, u) &= -\frac{u}{\tau_{in}}\phi^2(1-\phi) - \frac{\phi}{\tau_{out}}, \\
\mathcal{G}(\phi, u) &= \begin{cases} \frac{u-1}{\tau_{open}} & \text{if } \phi < \phi_{gate}, \\ \frac{u}{\tau_{close}} & \text{if } \phi \ge \phi_{gate}, \end{cases}
\end{aligned}
\tag{4.58}
\end{equation}
$$

where $0 < \phi_{\text{gate}} < 1$, $\tau_{\text{in}}$, $\tau_{\text{out}}$, $\tau_{\text{open}}$ and $\tau_{\text{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_{\text{in}} \ll \tau_{\text{out}} \ll \min(\tau_{\text{open}}, \tau_{\text{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}$

$$
\mathcal{G}_{\zeta}(\phi, u) = \left( \frac{1}{\tau_{\text{close}}} - \left( \frac{1}{\tau_{\text{close}}} - \frac{1}{\tau_{\text{open}}} \right) h_{\zeta}(\phi) \right) (u - h_{\zeta}(\phi)), \quad (4.59)
$$

where the differentiable function $0 \le h_{\zeta} \le 1$ is given by

$$
h_{\zeta}(\phi) = \frac{1}{2} \left( 1 - \tanh \left( \frac{\phi - \phi_{\text{gate}}}{\zeta} \right) \right),
$$

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

$$
\begin{align}
\mathcal{I}_{2,\zeta}(\phi) &= -\left(\frac{1}{\tau_{\text{close}}} - \left(\frac{1}{\tau_{\text{close}}} - \frac{1}{\tau_{\text{open}}}\right)h_{\zeta}(\phi)\right)h_{\zeta}(\phi), \tag{4.60} \\
\hbar_{\zeta}(\phi) &= \frac{1}{\tau_{\text{close}}} - \left(\frac{1}{\tau_{\text{close}}} - \frac{1}{\tau_{\text{open}}}\right)h_{\zeta}(\phi). \notag
\end{align}
$$

According to the definition of tanh, we can deduce that $\lim_{\zeta \to 0} h_\zeta(\phi) = \begin{cases} 1 & \text{if } \phi < \phi_{\text{gate}}, \\ 0 & \text{if } \phi > \phi_{\text{gate}} \end{cases}$ and then
$\lim_{\zeta \to 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 \mathbb{R}) \times \mathbb{R}^2$ into $\mathbb{R}$ and locally Lipschitz continuous function on $(\phi, u)$ and, $\mathcal{I}_2$ and $\hbar$
sufficiently regulars.
$\square$
## 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 \ge \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, \dots$), with $t_k = k\Delta t$ for $k = 0, 1, \dots$ and given time step $\Delta t$)

$$ 
\begin{aligned}
\partial_{0^+}^{\gamma} [\psi](\mathbf{x}, t) &= \frac{1}{\Gamma(1-\gamma)} \int_{0}^{t} (t-\tau)^{-\gamma} \frac{\partial \psi}{\partial t}(\mathbf{x}, \tau) d\tau, \\
&= \frac{1}{\Gamma(1-\gamma)} \sum_{m=0}^{n-2} \int_{t_m}^{t_{m+1}} (t-\tau)^{-\gamma} \frac{\partial \psi}{\partial t}(\mathbf{x}, \tau) d\tau + \frac{1}{\Gamma(1-\gamma)} \int_{t_{n-1}}^{t} (t-\tau)^{-\gamma} \frac{\partial \psi}{\partial t}(\mathbf{x}, \tau) d\tau, \quad (5.1) \\
&\approx -I_{\text{memory}}^{\gamma}(\psi_{\text{history}})(\mathbf{x}, t) + \frac{\Delta t^{1-\gamma}}{\Gamma(2-\gamma)} \frac{\partial \psi}{\partial t}(\mathbf{x}, t),
\end{aligned}
$$

where

$$ I_{\text{memory}}^{\gamma}(\psi_{\text{history}})(\mathbf{x}, t) = -\frac{\Delta t^{-\gamma}}{\Gamma(2-\gamma)} \sum_{m=0}^{n-2} [(n-m)^{1-\gamma} - (n-m-1)^{1-\gamma}][\psi(\mathbf{x}, t_{m+1}) - \psi(\mathbf{x}, t_m)]. \quad (5.2) $$

It is clear that $I_{\text{memory}}^1(\psi_{\text{history}})$, for $\gamma = 1$, is a null function. The state function $\psi_{\text{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_{\text{memory}}^\gamma$ representing hypothetical memory effects.

The term $I_{\text{memory}}^\gamma$ 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_{\text{memory}}^\gamma$ 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 $\ddot{\gamma}$ and then generates changes in myocardial polarization.

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

<footnote>‡Thus the depolarizing effect of the stimulus becomes repolarizing and vice versa</footnote>
results corresponding to the middle point of heart domain $\Omega_H$ for $\phi$, $\phi_e$ and $u$, and in the middle point of torso domain $\Omega_B$ for $\phi_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$) : $|x - a|^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$.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Heart-torso domain.</figcaption></figure>

We use some biophysical parameters which are presented in Table 1 and we fix $L_H = 15cm$ and

<table><caption>**Table 1.** Cell membrane parameters.</caption><thead><tr><th>Description</th><th>name</th><th>value (unit)</th></tr></thead><tbody><tr><td>Cell surface to volume ratio</td><td>κ</td><td>200 (cm<sup>-1</sup>)</td></tr><tr><td>transmembrane capacitance</td><td>C<sub>α</sub></td><td>10<sup>-3</sup> (F/cm<sup>2</sup>)</td></tr><tr><td>Depolarization length</td><td>T<sub>in</sub></td><td>4.5 (ms)</td></tr><tr><td>Repolarization length</td><td>T<sub>out</sub></td><td>90 (ms)</td></tr><tr><td>Opening time constant</td><td>T<sub>open</sub></td><td>100 (ms)</td></tr><tr><td>Closing time constant</td><td>T<sub>close</sub></td><td>130 (ms)</td></tr><tr><td>Change-over voltage</td><td>φ<sub>gate</sub></td><td>-67 (mV)</td></tr><tr><td>Resting potential</td><td>φ<sub>min</sub></td><td>-80 (mV)</td></tr><tr><td>Maximum potential</td><td>φ<sub>max</sub></td><td>20 (mV)</td></tr><tr><td>Activation time</td><td>T<sub>act</sub></td><td>10 (ms)</td></tr></tbody></table>

$L_T = 45cm$. In those both applications, we impose the following initial conditions

$$ \phi(x, y, 0) = \phi_{min}, \quad U(x, y, 0) = \frac{1}{(\phi_{max} - \phi_{min})^2}. \qquad (5.3) $$

Finally, we consider $T = 1$, $K_i = K_e = 0.003 \text{ s.cm}^{-1}\text{I}_d$ (with $\text{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

$$ f_H(\mathbf{x}, t) = I_{i,app}\chi_H(\mathbf{x})\chi_{[a_i, T_{act}+a_i]}(t)\chi_{prop}(\mathbf{x}, t)\Phi(\mathbf{x}), \quad \text{on } (a_i, a_i + T_{act}) \qquad (5.4) $$
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

$$ \chi_H(\mathbf{x}) = \begin{cases} 1 & \text{if } \mathbf{x} \in \Omega, \\ 0 & \text{else} \end{cases} \qquad (5.5) $$

$$ \chi_{[a_i, T_{act}+a_i]}(t) = \begin{cases} 1 & \text{if } a_i < t < T_{act} + a_i, \\ 0 & \text{else} \end{cases} \qquad (5.6) $$

$$ \chi_{prop}(\mathbf{x}, t) = \begin{cases} 1 & \text{if } x + y < 2t/T_{act}, \\ 0 & \text{else.} \end{cases} \qquad (5.7) $$

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

$$ \Phi(\mathbf{x}) = 1 - \frac{1}{L_H} \left(x - \frac{L_H}{2}\right)^2 - \frac{1}{L_H} \left(y - \frac{L_H}{2}\right)^2. \qquad (5.8) $$

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:

$$ \mathcal{I}(\phi, u) = -0.0004 \frac{\phi - \phi_{min}}{\phi_{max} - \phi_{min}} \left( U(\phi_{max} - \phi_{min})^2 + \frac{(\phi - \phi_{gate})(\phi - \phi_{max})}{(\phi_{max} - \phi_{min})^2} \right), \qquad (5.9) $$

and

$$ \mathcal{G}(\phi, u) = 0.63 \frac{\phi - \phi_{min}}{\phi_{max} - \phi_{min}} - 0.013 U(\phi_{max} - \phi_{min})^2. \qquad (5.10) $$

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$, $\phi_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 $\phi_e$, the effect fractional-order dynamic is negligible except for the extremum of $\phi_e$ (as the value of $\alpha$ becomes smaller, the value of the minimum point becomes smaller too).
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Response to a single stimulus. Time evolution of the transmembrane potential $\phi$ at midpoint of heart with $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Response to a single stimulus. Time evolution of the extracellular potential $\varphi_e$ at midpoint of heart with $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Response to a single stimulus. Time evolution of the ionic variable $u$ at midpoint of heart $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Response to a multiple stimuli (two stimuli in sequence). Time evolution of the transmembrane potential $\phi$ at midpoint of heart $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Response to a multiple stimuli (two stimuli in sequence). Time evolution of the extracellular potential $\varphi_e$ at midpoint of heart with $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Response to a multiple stimuli (two stimuli in sequence). Time evolution of the ionic variable $u$ at midpoint of heart with $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>

## 5.2. Mitchell-Schaeffer ionic 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:

$$ \begin{aligned} \mathcal{I}(\phi, u) &= -\frac{u(\phi - \phi_{min})^2(\phi_{max} - \phi)}{\tau_{in}(\phi_{max} - \phi_{min})} - \frac{1}{\tau_{out}}\frac{\phi - \phi_{min}}{(\phi_{max} - \phi_{min})}, \\ \mathcal{G}(\phi, u) &= \begin{cases} \frac{u}{\tau_{open}} - \frac{1}{\tau_{open}(\phi_{max} - \phi_{min})^2} & \text{if } \phi < \phi_{gate}, \\ \frac{u}{\tau_{close}} & \text{if } \phi \ge \phi_{gate}. \end{cases} \end{aligned} \quad (5.11) $$

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

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Response to a single stimulus. Time evolution of the transmembrane potential $\phi$ at midpoint of heart with $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Response to a single stimulus. Time evolution of the extracellular potential $\varphi_e$ at midpoint of heart with $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Response to a single stimulus. Time evolution of the ionic variable $u$ at midpoint of heart $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Response to a single stimulus. Time evolution of the torso potential $\varphi_s$ at midpoint of torso with $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Response to a multiple stimuli (two stimuli in sequence). Time evolution of the transmembrane potential $\phi$ at midpoint of heart with $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 15.** Response to a multiple stimuli (two stimuli in sequence). Time evolution of the extracellular potential $\varphi_e$ at midpoint of heart with $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 16.** Response to a multiple stimuli (two stimuli in sequence). Time evolution of the ionic variable $u$ at midpoint of heart with $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 17.** Response to a multiple stimuli (two stimuli in sequence). Time evolution of the torso potential $\varphi_s$ at midpoint of torso with $\beta = 1$ and $\alpha = -1, -0.925, -0.9, -0.85$.</figcaption></figure>

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.

<figure><img src="image_22.png" /><figcaption>**Figure 18.** Response to a single stimulus. Time evolution of the transmembrane potential $\phi$ at midpoint of heart with $\alpha = 1$ and $\beta = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 19.** Response to a single stimulus. Time evolution of the extracellular potential $\varphi_e$ at midpoint of heart with $\alpha = 1$ and $\beta = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 20.** Response to a single stimulus. Time evolution of the ionic variable $u$ at midpoint of heart with $\alpha = 1$ and $\beta = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 21.** Response to a single stimulus. Time evolution of the torso potential $\varphi_s$ at midpoint of torso with $\alpha = 1$ and $\beta = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure 22.** Response to a multiple stimuli (two stimuli in sequence). Time evolution of the transmembrane potential $\phi$ at midpoint of heart with $\alpha = 1$ and $\beta = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 23.** Response to a multiple stimuli (two stimuli in sequence). Time evolution of the extracellular potential $\varphi_e$ at midpoint of heart with $\alpha = 1$ and $\beta = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Figure 24.** Response to a multiple stimuli (two stimuli in sequence). Time evolution of the ionic variable $u$ at midpoint of heart with $\alpha = 1$ and $\beta = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 25.** Response to a multiple stimuli (two stimuli in sequence). Time evolution of the torso potential $\varphi_s$ at midpoint of torso with $\alpha = 1$ and $\beta = -1, -0.925, -0.9, -0.85, -0.8$.</figcaption></figure>

## 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 [1, 10, 12, 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.).

**Conflict of interest**

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

**References**

1. R. A. Adams, *Sobolev spaces*, Academic Press, New-York, 1975.

2. O. P. Agrawal, Fractional variational calculus in terms of Riesz fractional derivatives, *J. Phys. A*, **40** (2007), 6287–6303.

3. R. Almeida, N. R. Bastos, M. T. Monteiro, Modeling some real phenomena by fractional differential equations, *Math. Methods Appl. Sci.*, **39** (2016), 4846–4855.

4. T. Arpadffy-Lovás, I. Baczko, B. Balati, M. Bitay, N. Jost, C. Lengyel, et. al., Electrical Restitution and its modifications by antiarrhythmic drugs in undiseased human ventricular muscle, *Front. Pharmacol.*, **11**, (2020), 479.

<footnote>§from some observations, cardiac MRI-Measurement or desired target</footnote>
5. T. Ashihara, J. Constantino, N. A. Trayanova, Tunnel propagation of postshock activations as a hypothesis for fibrillation induction and isoelectric window, *Circ. Res.*, **102** (2008), 737-745.

6. O. V. Aslanidi, A. P. Benson, M. R. Boyett, H. Zhang, Mechanisms of defibrillation by standing waves in the bidomain ventricular tissue with voltage applied in an external bath, *Physica D*, **238** (2009), 984-991.

7. G. W. Beeler, H. Reuter, Reconstruction of the action potential of ventricular myocardial fibres, *J. Physiol.*, **268** (1977), 177-210.

8. A. Belmiloudi, Time-varying delays in electrophysiological wave propagation along cardiac tissue and minimax control problems associated with uncertain bidomain type models, *AIMS Mathematics*, **4** (2019), 928-983.

9. A. Belmiloudi, S. Corre, Mathematical modeling and analysis of dynamic effects of multiple time-varying delays on electrophysiological wave propagation in the heart, *Nonlinear Anal-Real*, **47** (2019), 18-44.

10. A. Belmiloudi, Mathematical modeling and optimal control problems in brain tumor targeted drug delivery strategies, *Int. J. Biomath.*, **10** (2017), 1750056.

11. A. Belmiloudi, Dynamical behavior of nonlinear impulsive abstract partial differential equations on networks with multiple time-varying delays and mixed boundary conditions involving time-varying delays, *J. Dyn. Control Syst.*, **21** (2015), 95-146.

12. A. Belmiloudi, Robust control problem of uncertain bidomain models in cardiac electrophysiology, *Journal of Coupled Systems and Multiscale Dynamics*, **1** (2013), 332-350.

13. A. Belmiloudi, *Stabilization, optimal and robust control. Theory and applications in biological and physical sciences*, Springer-Verlag, London, 2008.

14. A. Belmiloudi, Robust control problems associated with time-varying delay nonlinear parabolic equations, *IMA J. Math. Control I.*, **20** (2003), 305-334.

15. A. Belmiloudi, Regularity results and optimal control problems for the perturbation of Boussinesq equations of the ocean, *Numer. Func. Anal. Opt.*, **21** (2000), 623-651.

16. M. Bendahmane, K. H. Karlsen, Analysis of a class of degenerate reaction-diffusion systems and the bidomain model of cardiac tissue, *Netw. Heterog. Media*, **1** (2006), 185-218.

17. G. A. Bocharov, F. A. Rihan, Numerical modelling in biosciences using delay differential equations, *J. Comput. Appl. Math.*, **125** (2000), 183-199.

18. M. Boulakia, M. A. Fernandez, J. F Gerbeau, N. Zemzemi, A coupled system of PDEs and ODEs arising in electrocardiogramms modeling, *Applied Math. Res. Exp.*, **28** (2008).

19. Y. Bourgault, Y. Coudiere, C. Pierre, Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology, *Nonlinear Anal-Real*, **10** (2009), 458-482.

20. N. Buric, K. Todorovic, N. Vasovic, Dynamics of noisy FitzHugh-Nagumo neurons with delayed coupling, *Chaos, Solitons and Fractals*, **40** (2009), 2405-2413.

21. J. O. Campos, R. S. Oliveira, R. W. dos Santos, B. M. Rocha, Lattice Boltzmann method for parallel simulations of cardiac electrophysiology using GPUs, *J. Comput. Appl. Math.*, **295** (2016), 70-82.
22. M. Caputo, Linear models of dissipation whose Q is almost frequency independent–II, *Geophys. J. Int.*, **13** (1967), 529–539.

23. E. M. Cherry, F. H. Fenton, T. Krogh-Madsen, S. Luther, U. Parlitz, Introduction to Focus Issue: Complex Cardiac Dynamics, *Chaos*, **27** (2017), 093701.

24. P. Colli Franzone, G. Savaré, *Degenerate evolution systems modeling the cardiac electric field at micro- and macroscopic level*. In: Evolution Equations, Semigroups and Functional Analysis (eds.A. Lorenzi, B. Ruf), Birkhauser, Basel, 49–78, 2002.

25. P. Colli Franzone, L. Guerri, S. Tentoni, Mathematical modeling of the excitation process in myocardial tissue: Influence of fiber rotation on wavefront propagation and potential field, *Math. Biosci.*, **101** (1990), 155–235.

26. T. Comlekoglu, S. H. Weinberg, Memory in a fractional-order cardiomyocyte model alters properties of alternans and spontaneous activity, *Chaos*, **27** (2017), 093904.

27. S. Corre, A. Belmiloudi, Coupled lattice Boltzmann simulation method for bidomain type models in cardiac electrophysiology with multiple time-delays, *Math. Model. Nat. Pheno.*, **14** (2019), 207.

28. S. Corre, A. Belmiloudi, Coupled lattice Boltzmann method for numerical simulations of fully coupled Heart and Torso bidomain system in electrocardiology, *Journal of Coupled System and Multiscale Dynamics*, **4** (2016), 207–229.

29. S. Corre, A. Belmiloudi, *Coupled Lattice Boltzmann Modeling of Bidomain Type Models in Cardiac Electrophysiology*, Mathematical and Computational Approaches in Advancing Modern Science and Engineering (eds. J. Bélair, et al.), Springer, 209–221, 2016.

30. H. Dal, S. Goktepe, M. Kaliske, E. Kuhl, A fully implicit finite element method for bidomain models of cardiac electromechanics, *Comput. Method. Appl. M.*, **253** (2013), 323–336.

31. S. Das, K. Maharatna, Fractional dynamical model for the generation of ECG like signals from filtered coupled Van-der, *Comput. Methods Programs Biomed.*, **122** (2013), 490–507.

32. K. Diethelm, *Fractional differential equations*, Springer, Berlin, 2010.

33. K. Diethelm, N. J. Ford, Analysis of fractional differential equations, *J. Math. Anal. Appl.*, **265** (2002), 229–248.

34. S. Dipierro, E. Valdinoci, A simple mathematical model inspired by the Purkinje cells: from delayed travelling waves to fractional diffusion, *Bull. Math. Biol.*, **80** (2018), 1849–1870,

35. M. Dupraz, S. Filippi, A. Gizzi, A. Quarteroni, R. Ruiz-Baier, Finite element and finite volume-element simulation of pseudo-ECGs and cardiac alternans, *Math. Method. Appl. Sci.*, **38** (2015), 1046–1058.

36. D. Baleanu, A. M. Lopes, (eds.) *Applications in engineering, life and social sciences*. In: Handbook of fractional calculus with applications, De Gruyter, 2019.

37. K. A. Ellenbogen, B. L. Wilkoff, G. N. Kay, C. P. Lau, A. Auricchio, *Clinical cardiac pacing, defibrillation and resynchronization therapy*, Elsevier - Health Sciences Division; 5th Revised edition, 2017.

38. F. Fenton, A. Karma, Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation, *Chaos*, **8** (1998), 20–47.
39. R. Fitzhugh, Impulses and physiological states in theoretical models of nerve membrane, *Biophys. J.*, **1** (1961), 445-465.

40. L. Glass, Synchronization and rhythmic processes in physiology, *Nature*, **410** (2001), 277-284.

41. J. M. Gomes, R. Weber dos Santos, E. M. Cherry, Alternans promotion in cardiac electrophysiology models by delay differential equations, *Chaos*, **27** (2017), 093915.

42. G. H. Hardy, J. E. Littlewood, Some properties of fractional integrals I, *Math. Z.*, **27** (1928).

43. D. A. Israel, D. J. Edell, R. G. Mark, Time delays in propagation of cardiac action potential, *Am. J. Physiol.*, **258** (1990), H1906-17.

44. D. Jeyaraj, M. Ashwath, D. S. Rosenbaum, Pathophysiology and clinical implications of cardiac memory, *Pacing Clin Electrophysiol*, **33** (2010), 346-352.

45. S. J. Kalbfleisch, J. Sousa, R. el-Atassi, H. Calkins, J. Langberg, F. Morady, Repolarization abnormalities after catheter ablation of accessory atrioventricular connections with radiofrequency current, *J. Am. Coll. Cardiol.*, **18** (1991), 1761-1766.

46. J. Keener, J. Sneyd, *Mathematical Physiology*, Springer-Verlag, New York, 2009.

47. A. A. Kilbas, H. M. Srivastava, J. J. Trujillo, *Theory and applications of fractional differential equations*, 2006.

48. A. Kubica, M. Yamamoto, Initial-boundary value problems for fractional diffusion equations with time-dependent coefficients, *Fract. Calc. Appl. Anal.*, **21** (2018), 276-311.

49. G. W. Leibniz, Letter from Hanover, Germany, to G.F.A. L'Hôpital, September 30; 1695, *Mathematische Schriften*, **2** (1849), 301-302.

50. L. Li, J. G. Liu, A generalized definition of Caputo derivatives and its application to fractional ODEs, *SIAM J. Math. Anal.*, **50** (2018), 2867-2900.

51. G. Lines, M. Buist, P. Grottum, A. J. Pullan, J. Sundnes, A. Tveito, Mathematical models and numerical methods for the forward problem in cardiac electrophysiology, *Computing and Visualization in Science*, **5** (2003), 215-239.

52. J. L. Lions, E. Magenes, *Problèmes aux limites non homogènes et applications*, Tome 1&2, Dunod, Paris, 1968.

53. C. H. Luo, Y. Rudy, A model of the ventricular cardiac action potential. Depolarization, repolarization, and their interaction, *Circ. Res.*, **68** (1991), 1501-1526.

54. C. H. Luo, Y. Rudy, A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes, *Circ. Res.*, **74** (1994), 1071-1096.

55. R. L. Magin, Fractional calculus models of complex dynamics in biological tissues, *Comput. Math. Appl.*, **59** (2010), 1586-1593.

56. C. C. Mitchell, D. G. Schaeffer, A two-current model for the dynamics of cardiac membrane, *Bull. Math. Biol.*, **65** (2003), 767-793.

57. J. Roger, A. McCulloch, A collocation-Galerkin finite element model of cardiac action potential propagation, *IEEE T. Biomed. Eng.*, **41** (1994), 743-757.
58. T. Nakagawa, T. Yagi, A. Ishida, Y. Mibiki, Y. Yamashina, H. Sato, et al., Differences between cardiac memory Twave changes after idiopathic left ventricular tachycardia and ischemic T wave inversion induced by acute coronary syndrome, *J. Electrocardiol.*, **49** (2016), 596-602.

59. N. Ozgen, Z. Lu, G. J. Boink, D. H. Lau, I. N. Shlapakova, Y. Bobkov, et al., Microtubules and angiotensin II receptors contribute to modulation of repolarization induced by ventricular pacing, *Heart Rhythm*, **9** (2012), 1865-72.

60. L. Padeletti, C. Fantappie, L. Perrotta, G. Ricciardi, P. Pieragnoli, M. Chiostri, et al., Cardiac memory in humans: vectocardiographic quantification in cardiac resynchronization therapy, *Clin. Res. Cardiol.*, **100** (2011), 100-51.

61. A. Panfilov, R. Aliev, A simple two-variable model of cardiac excitation, *Chaos Solitons and Fractals*, **7** (1996), 293-301.

62. I. Petras, *Fractional-order chaotic systems*, Fractional-order nonlinear systems, Springer, 2011.

63. A. N. Plotnikov, A. Shvilkin, W. Xiong, J. R. de Groot, L. Rosenshtraukh, S. Feinmark, et al., Interactions between antiarrhythmic drugs and cardiac memory, *Cardiovasc. Res.*, **50** (2001), 335-344.

64. A. J. Pullan, M. L. Buist, L. K. Cheng, *Mathematically modelling the electrical activity of the heart-from cell to body surface and back again*, World Scientific, Singapore, 2005.

65. A. Quarteroni, R. Sacco, F. Saleri, *Numerical Mathematics*, Vol. 37 of Texts in Applied Mathematics (2nd ed.), Springer-Verlag, Berlin, 2007.

66. M. B. Rosenbaum, H. H. Blanco, M. V. Elizari, J. O. Lazzari, J. M. Davidenko, Eletrotonic modulation of the T-wave and cardiac memory, *Am. J. Cardiol.*, **50** (1982), 213-222.

67. Y. Sakamoto, Y. Inden, H. Okamoto, K. Mamiya, T. Tomomatsu, A. Fujii, et al., T-wave changes of cardiac memory caused by frequent premature ventricular contractions originating from the right ventricular outflow tract, *J. Cardiovasc. Electrophysiol.*, **30** (2019), 1549-1556.

68. S. G. Samko, A. A. Kilbas, O. I. Marichev, *Fractional integrals and derivatives*, Gordon and Breach Science Publishers, Yverdon, 1993. Theory and applications; Edited and with a foreword by S. M. Nikolski; Translated from the 1987 Russian original; Revised by the authors.

69. R. F. Sandra, M. A. Savi, An analysis of heart rhythm dynamics using a three-coupled oscillator model, *Chaos, Solitons and Fractals*, **41** (2009), 2553-2565.

70. T. Seidman, H. Z. Zhou, Existence and uniqueness of optimal controls for a quasilinear parabolic equation, *SIAM J. Control Optim.*, **20** (1982), 747-762.

71. A. Shvilkin, H. D. Huang, M. E. Josephson, Cardiac memory: diagnostic tool in the making, *Circ-Arrhythmia Elect.*, **8** (2015), 475-482.

72. A. Shvilkin, K. K. Ho, M. R. Rosen, M. E. Josephson, T-vector direction differentiates post-pacing from ischemic T-wave inversion in precordial leads, *Circulation*, **111** (2005), 969-974.

73. E. A. Sosunov, E. P. Anyukhovsky, M. R. Rosen, Altered ventricular stretch contributes to initiation of cardiac memory, *Heart Rhythm*, **5** (2008), 106-113.

74. J. Sundnes, G. Lines, X. Cai, B. F. Nielsen, K. A. Mardal, A. Tveito, *Computing the electrical activity in the heart*, Springer, Berlin, 2006.
75. M. C. Suran, A. D. Margulescu, R. Bruja, D. Vinereanu, Surface ECG criteria can discriminate post-septal pacing cardiac memory from ischemic T wave inversions, *J. Electrocardiol.*, **58** (2020), 10-17.

76. W. Tan, C. Fu, C. Fu, W. Xie, H. Cheng, An anomalous subdiffusion model for calcium spark in cardiac myocytes, *Appl. Phys. Lett.*, **91** (2007), 183901.

77. R. Temam, *Navier-Stokes equations*, North-Holland, Amsterdam, 1984.

78. L. Tung, *A bi-domain model for describibg ischemic myocardial d-c potentials* [Thesis/Dissertation], Massachusetts Institute of Technology, 1978.

79. K. H. Ten Tusscher, A. V. Panfilov, Alternans and spiral breakup in a human ventricular tissue model, *Am. J. Physiol-Heart C.*, **291** (2006), H1088-H1100.

80. M. Veneroni, Reaction-diffusion systems for the macroscopic bidomain model of the cardiac electric field, *Nonlinear Anal-Real*, **10** (2009), 849-868.

81. E. J. Vigmond, R. W. dos Santos, A. J. Prassl, M. Deo, G. Plank, Solvers for the cardiac bidomain equations, *Prog. Biophys. Mol. Bio.*, **96** (2008), 3-18.

82. J. W. Waks, D. A. Steinhaus, A. Shvilkin, D. B. Kramer, Post-pacemaker T-wave Inversions: Cardiac Memory, *Am. J. Med.*, **129** (2016), 170-172.

83. B. J. West, M. Turalska, P. Grigolini, *Networks of echoes: imitation, innovation and invisible leaders*, Springer Science & Business Media, 2014.

84. S. Westerlund, L. Ekstam, Capacitor theory, *IEEE T. Dielect. El. In.*, **1** (1994), 826-839.

85. H. Ye, J. Gao, Y. Ding, A generalized Gronwall inequality and its application to a fractional differential equation, *J. Math. Anal. Appl.*, **328** (2007), 1075-1081.

86. Y. Zhou, *Basic theory of fractional differential equations*, World Scientific, Singapore, 2014.

<figure><img src="image_30.png" /></figure>

AIMS Press

©2021 the Author(s), licensee AIMS Press. This
is an open access article distributed under the
terms of the Creative Commons Attribution License
(http://creativecommons.org/licenses/by/4.0)