High-order discontinuous Galerkin methods for the monodomain and bidomain models

1.   Introduction

The heart activity is defined through the cardiac cycle which, at a first analysis, is characterized by two alternating phases: the diastole, the period of relaxation, and the systole, the period of contraction. The continuous sequence of contractions, formally called rhythmicity, is caused by the propagation of bio-electrical signals through the cells resulting in an active contraction of the cardiac muscle. Mathematical models are widely employed to understand and predicting the complex processes underlying this phenomenon. The most popular ones are the monodomain and the bidomain models, which mathematically correspond to systems of (non linear) reaction-diffusion equations coupled to a system of ordinary differential equations for the ionic currents ^{[6, 15]}.

The cardiac electrical propagation is represented by a wave characterized by a fast and steep front. Therefore, in order to accurately capture the proper dynamics, Finite Element methods usually require very refined computational grids. This clearly entails high computational costs to solve the electrophysiology problem in real scenarios. On the other hand, high-order spectral methods have become popular for their ability at capturing sharp parts such as shock waves. In this context, discontinuous Galerkin (DG) methods may provide an effective alternative for the solution of the electro-physiology problems, guaranteeing high-order accuracy and more flexibility, see ^{[10, 17]} and for the problem at hand ^[12].

The objective of this work is to present and test in practice a high-order DG method to discretize the monodomain and bidomain models. In particular, Section 2 is aimed at introducing the electrophysiology problem and the associated mathematical models. The DG formulation is illustrated in Section 3 while Section 4 provides the details of the high-order approximation. In order to obtain the final numerical algorithm, Section 5 incorporates the semi-implicit time discretization. Finally, Section 6 presents some numerical tests. Specifically, we perform a preliminary analysis in model problems to highlight the feasibility of this approach showing some convergence tests refining the mesh size and increasing the polynomial order. Then, the results of a realistic (although simplified) simulation of the propagation of the electric potential is shown to highlight the effectiveness of the proposed method. Final remarks are then discussed.

2.   The physical problem and the mathematical formulation

In what follows we present a brief introduction to the monodomain and bidomain equations. For a detailed derivation, we address the interested reader to ^[15].

The heart's active mechanical contraction is triggered by the cardiac cell's electrical activation. Cardio-myocytes are activated and deactivated at each heartbeat following a characteristic electrical cycle. The cell is initially at rest ($-90 \, mV$). At the beginning of the activation, its potential increases rapidly ($\approx2\, ms$) and reaches the value of $+20\, mV$. Later, a plateau near $0\, mV$ is observed, followed by a slow repolarization to the initial potential.

From a microscopical point of view, each single cell is involved in a passage of chemical ions through specific channels, e.g., calcium $Ca^{++}$, sodium $Na^+$ and potassium $K^+$. From a macroscopical point of view, one can describe the dynamics as a continuous electrical diffusion over the entire cardiac tissue driven by the directions of the muscle fibers ^[14].

Applying general electromagnetism laws, the bidomain model has been formulated for the macroscopic dynamics (see ^{[6, 15]} for more details). To complete the formulation, a mathematical model for the ionic current is required. In such context, it is noteworthy to mention the original formulation by Hodgkin and Huxley ^[11]. FitzHugh ^[9] and Nagumo ^[13] proposed later a simplification of the latter and complex models successively followed such as the ones proposed in ^{[2, 20, 21]}. In this work, the FitzHugh-Nagumo reduced ionic model (FHN) is considered. This simple model represents only one ionic channel o describe the ionic currents.

Given an open and bounded domain $\Omega \in \mathbb{R}^d$, $d = 2, 3$, and a final time $T > 0$, the unknowns of the bidomain model are: the trans-membrane potential $V_m = \phi_i - \phi_e$, where $\phi_i$ and $\phi_e$ are the intracellular and extracellular potentials, respectively, and the gating variable $w$ representing the percentage of opening of the ionic channel. The model parameters include: the positive constants $\chi_m, C_m$, representing the surface area-to-volume ratio and the membrane capacitance, the permeability tensors $\Sigma_i, \Sigma_e$ in the internal and external cellular field, the external applied currents $I_i^{ext}, I_e^{ext}$, and some known constants to tune the ionic model ($\kappa, a, \epsilon, \Gamma$). In particular, $\Sigma_i$ and $\Sigma_e$ account for the anisotropy given by the cardiac fibers. Furthermore, initial and Neumann boundary conditions are imposed through some known functions $\phi_{i, 0}, \phi_{e, 0}, w_{0}, b_i, b_e$. The former conditions assign the initial state of the system while the latter conditions prescribe the behaviour on the boundary of the domain in terms of inward or outward currents.

Problem 1 (Bidomain model coupled with FitzHugh-Nagumo model). For each $t\in(0, T]$, find $\phi_i$, $\phi_e$ and $w$ such that:

Note that Problem 1 needs the following compatibility condition, which is necessary for the existence of the solution:

A simpler formulation can be derived assuming the intracellular and extracellular permeability tensors to be proportional ^{[6, 15]}. This assumption is valid when the propagation is regular and no chaotic processes such as fibrillation are taking place. In such a way, the unknowns reduce to $V_m$ and $w$. Therefore, the monodomain problem reads as follows:

Problem 2 (Monodomain problem coupled with FitzHugh-Nagumo model). For each $t\in(0, T]$, find $V_m$ and $w$ such that:

where

$\xi$ being the proportional factor: $\Sigma_e = \xi\Sigma_i$, and for suitable initial and boundary conditions which derive from the bidomain ones.

3.   Semi-discrete discontinuous Galerkin approximation

Starting from the strong form given by Problems 1 and 2, the next step is the pursuit of a DG semi-discrete formulation. Terms and symbols are defined similarly to the usual convention ^{[3, 4]}.

3.1.   Discontinuous Galerkin formulation

Let us introduce a shape-regular triangulation $\mathcal{T}_h$ of $\Omega$, where $\mathcal{F} _h = \mathcal{F} _h^I \cup \mathcal{F} _h^B$ is the set of the faces of the partition which includes the internal and boundary faces, respectively. Let the DG space be defined as

where $\mathbb{P}^p(\mathcal{K})$ is the space of polynomials of total degree less than or equal to $p\ge1$ over $\mathcal{K}\in \mathcal{T}_h$ and $\cdot|_\mathcal{K}$ is the restriction operator to the element $\mathcal{K}$. Moreover, we define $N_h = \dim(\Theta^{h, p}) < \infty$ as the dimension of the space $\Theta^{h, p}$ and the following bilinear forms.

Here $\Sigma_*$ is either $\Sigma_i$ or $\Sigma_e$, $\theta \in \{-1, 0, 1\}$, $\nabla_h$ denotes the element-wise gradient and we set

where $\overline{\gamma}$ is a stability parameter defined edge-wise, with $h > 0$ the mesh size (supposed to be quasi uniform) and $\alpha > 0$ a fixed parameter, and $\mathbf{n} \in \mathbb{R}^d$ is the outward normal unitary vector to the corresponding face $F$. Moreover, the jump and average operators $[\![\cdot]\!]$ and $\{\!\{\cdot \}\!\}$ are defined on $F\in \mathcal{F}_h^I$ in the standard way ^[4], i.e.,

where the subscript $\{1, 2\}$ indicates the evaluation of the variable $u$, the vector variable $\mathbf{v}$ or the normal vector $\mathbf{n}$ with respect to the two adjacent elements $\mathcal{K}_1$ and $\mathcal{K}_2$ such that

Then, we can write the semi-discretized formulation for the bidomain problem.

Problem 3 (Bidomain model - semidiscrete DG formulation). For any $t\in(0, T]$, find $\phi_i^h(t), \phi_e^h(t), w^h(t) \in \Theta^{h, p}$ such that:

for each $v_h\in \Theta^{h, p}$, where

$I_{ion}(\cdot, \cdot)$ is defined in Problem 1, $\phi_{i, 0}^h, \phi_{e, 0}^h, w^h_0 \in \Theta^{h, p}$ are suitable projections onto $\Theta^{h, p}$ of $\phi_{i, 0}, \phi_{e, 0}, w_0$.

Following the same notation, we also present the weak formulation for the monodomain problem.

Problem 4 (Monodomain model - DG semidiscrete formulation). For any $t\in(0, T]$, find $V_m^h(t) \in \Theta^{h, p}$ and $w^h(t) \in \Theta^{h, p}$ such that:

for every $v_h \in \Theta^{h, p}$, where $I_{ion}$ is defined in Problem 2 and $V_{m, 0}^h, w^h_0 \in \Theta^{h, p}$ are suitable projections onto $\Theta^{h, p}$ of $V_{m, 0}, w_0$.

In Problems 3 and 4, three different methods are employed according to the choice of the coefficient $\theta$:

● $\theta = 1$: Symmetric Interior Penalty method (SIP) ^[7];

● $\theta = 0$: Incomplete Interior Penalty method (IIP) ^[18];

● $\theta = -1$: Non Symmetric Interior Penalty method (NIP) ^[16].

Notice that $\alpha$ in the definition (3.1) needs to be large enough to guarantee coercivity of the SIP and IIP formulations.

3.2.   Algebraic formulation

Let $\{\varphi_j\}_{j = 1}^{N_h}$ be a basis of $\Theta^{h, p}$ and let $\boldsymbol{V}^h_m(t), \boldsymbol{\phi}^h_i(t), \boldsymbol{\phi}^h_e(t), \boldsymbol{w}^h(t) \in \mathbb{R}^{N_h}$ be the vectors containing the expansion coefficients of $V_m^h(t), \phi_i^h(t), \phi_e^h(t), w^h(t)$ with respect to such a basis for each instant of time $t\in (0, T]$.

For $i, j = 1, \cdots, N_h$ and $z = \{i, e, \emptyset\}$, we define the following matrices and right hand side vector:

The semi-discrete DG algebraic formulation of the bidomain problem leads to the following ODE system:

Problem 5 (Algebraic formulation of the bidomain model). Find

such that for any $t \in (0, T]$:

together with initial conditions.

The monodomain DG ODE system can be written similarly.

Problem 6 (Algebraic formulation of the monodomain model). Find $\boldsymbol{V}_m^h(t), \boldsymbol{w}^h(t) \in \mathbb{R}^{N_h}$ such that for any $t \in (0, T]$:

together with initial conditions.

4.   A note on the basis functions

The Dubiner spectral functions ^[8] are chosen as basis $\{\varphi_j\}_{j = 1}^{N_h}$ for the space $\Theta^{h, p}$. The advantage of using an $L^2$ orthogonal basis on simplices rather than a nodal basis relies on the feasibility of applying high-order discretizations without incurring in high costs constructing the discretized space and/or ill-conditioned matrices. Moreover, the mass matrix is in this case diagonal.

On the other hand, the projection of known functions on the spectral space is certainly more involved since it requires $L^2$ scalar products. This issue can be overcome by implementing fast evaluations of Jacobi polynomials and efficient quadrature formulae to be also used in the weak formulations of Problems 3 and 4.

In $d = 2$, the Dubiner basis is defined on the reference simplex as follows ^[3]. This is possible by collapsing the two-dimensional reference square

into the reference triangle

applying the following transformation (Figure 1):

Hence, the Dubiner basis functions are constructed as the transformations of suitable basis functions initially defined on $\hat{Q}$. More precisely, on the reference square $\hat{Q}$ the basis functions are defined as modified tensor products of Jacobi polynomials.

Definition 4.1 (Jacobi polynomials ^[19]). The $n$-th Jacobi polynomial of indices $\alpha, \beta \ge -1$ is defined as:

Definition 4.2 (2D Dubiner basis functions ^[3]). If $d = 2$, the Dubiner basis function indexed by $(i, j) \in \mathbb{N}^2, i+j\le p$ is defined as:

where

Remark 4.1. Despite the great potentiality of an orthogonal spectral basis, it is often necessary to pass to nodal representations. Let $\{\varphi_j\}_{j = 1}^{N_h}$ be the set of Dubiner basis functions and let

A nodal evaluation of such solution can be performed through the linear combination

while the converse operation can be performed through a $L^2$ scalar product thanks to the orthonormality property:

5.   Fully discrete formulation

In this section a first order semi-implicit time discretization is presented in order to obtain a fully discretized system of equations from Problems 3 and 4. Thus, we split the interval $(0, T]$ into $N$ uniform sub-intervals $(t^n, t^{n+1}]$ of length $\Delta t$ such that $t^n = n \Delta t \quad \forall n = 0, \cdots, N-1$. Then, we define the fully discretized solutions $V_m^{h, n}, w^{h, n}, \phi_i^{h, n}, \phi_e^{h, n} \in \Theta^{h, p}$ that are approximations of $V_m^h(t^n)$, $w^h(t^n)$, $\phi_i^h(t^n)$, $\phi_e^h(t^n)$, respectively, whose expansion coefficients vectors are denoted with $\boldsymbol{V}_m^{h, n}$, $\boldsymbol w^{h, n}$, $\boldsymbol{\phi}_i^{h, n}$, $\boldsymbol{\phi}_e^{h, n} \in \mathbb{R}^{N_h}$, respectively.

In the semi-implicit scheme, only the non-linear term is treated in explicit in order to take advantage of an implicit discretization while still having linearity in the time-step advancing scheme.

Problem 7 (Fully discretized formulation of the Bidomain model). Given $\phi_i^{h, 0}, \phi_e^{h, 0}, w^{h, 0} \in \Theta^{h, p}$, find

for $n = 0, \ldots, N-1$ such that:

Problem 8 (Fully discretized formulation of the Monodomain model). Given $V_m^{h, 0}, w^{h, 0} \in \Theta^{h, p}$, find $\boldsymbol{V}_m^{h, n+1}, \boldsymbol{w}^{h, n+1} \in \mathbb{R}^{N_h}$ for $n = 0, \ldots, N-1$ such that:

Remark 5.1. Problems 1 and 2 are known to be well-posed in the sense of existence, uniqueness and regularity of the solutions under suitable assumptions ^{[5, 6]}. More precisely, $V_m(t)$ is unique in $H^1(\Omega)$ $\forall t \in (0, T]$ while $\phi_i(t)$ and $\phi_e(t)$ are unique in $H^1(\Omega)/\mathbb{R}$, i.e., up to an additive constant. Since $V_m$ is defined as the difference between the two potentials, this constant is necessarily the same for both $\phi_i$ and $\phi_e$. Therefore, the numerical solution of Problem 7 requires to fix the constant, e.g., by imposing the value of $\phi_i$ or $\phi_e$ in a particular point of the domain or by imposing its mean value. In the former case, the condition

is applied for arbitrary $\mathbf{x}\in\Omega$, $c\in\mathbb{R}$, $\forall n = 1, \dots, N$. We opt instead for the latter strategy which is more consistent with the variational formulation of the problem in $H^1(\Omega)$ and is represented by the condition

or

where $d_j: = \int_\Omega \varphi_j$, for an arbitrary $c\in\mathbb{R}$, $\forall n = 1, \dots, N$.

6.   Numerical results

This section is devoted to some numerical tests aimed at assessing the quality of the proposed numerical strategies. Some convergence tests are discussed in Section 6.1 while a more realistic simulation is shown in Section 6.2. The code has been implemented in a small C++ library, called $\texttt{DUBeat}$^*, which is largely based on $\texttt{lifex}$ ^[1].

https://matteocalafa.com/DUBeat/

6.1.   Verification test cases

We consider a simple square domain $\Omega = (0, 1)^2 \subset \mathbb{R}^2$. We opt for the symmetric interior penalty method (SIP, $\theta = 1$) and the following parameters choices:

Moreover,

where $\mathcal{I}_2 \in M_{2, 2}(\mathbb{R})$ is the $2$-dimensional identity matrix, i.e., we assume that the coordinate system is already aligned with the principal fibers directions and there is no anisotropy.

Boundary conditions and applied currents are assigned assuming that the exact solutions are

with $\phi_e$ being shown in Figure 2.

In the following tests, the errors are computed with respect to the standard $L^2(\Omega)$ norms, the $H^1(\Omega)$ norm

and the DG norm

We tested the proposed formulations in two dimensions on a sequence of uniformly refined grids of granularity $h\approx2^{-\sigma}$, $\sigma = 0, 1, \dots, 5$, see Figure 3.

Figure 4 shows the computed errors obtained with the monodomain and bidomain problems in two dimensions. Except for the pre-asymptotic regime when $h\approx 1$, the errors for $d = 2$ shows the expected theoretical convergence orders that are shown in the bottom triangles. If $p\in \mathbb{N}$ is the polynomial degree, the order is indeed $O(h^p)$ for the $H^1$ and $DG$ errors and $O(h^{p+1})$ for the $L^\infty$ and $L^2$ errors. We also observe that the curves are more flat for $h\ll1, p = 2$ because the error due to the grid refinement is very small and therefore comparable to the error caused by the time advancing scheme.

Correct convergence orders are also verified for $p$-refinements in Figure 5 where the grid size $h = 2^{-3}$ is fixed and the polynomial degree varies from 1 to 5.

6.2.   Towards a realistic simulation

In this section a first pseudo-realistic simulation is performed in order to assess the coherence between the monodomain model, the numerical discretization and the physical phenomenon. We aim at analysing the behaviour of a localized applied current on a bounded portion of heart tissue. The domain $\Omega \subset \mathbb{R}^2$ is defined as a reference square as in Section 6.1 and boundary conditions are homogeneous so that the only source of potential is represented by $I^{ext}$ while the domain $\Omega$ is isolated from the outer regions. The applied current is defined as

where $\mathcal{I}_{[a, b]}(x)$ is the indicator function, i.e., $\mathcal{I}_{[a, b]}(x) = 1$ if $x \in [a, b]$ and 0 otherwise. This definition of $I^{ext}$ represents a temporary and localized electric shock in the center of the domain. Conductivity and ionic model parameters are chosen as in Section 6.1 except for $C_m = 10^{-2}$ and $\epsilon = 40$. This latter choice allows to obtain electric activation and repolarization in reasonable times. Then, the time step is $\Delta t = 10^{-3}$, second order polynomials ($p = 2$) are employed as well as a uniformly refined grid with $h = 2^{-6}$.

Finally, results are shown in Figure 6 where the activation and subsequent propagation of the electric potential are well visible. The maximum value $V_m = 1$ is achieved on the crest while the internal part undergoes repolarization where the potential falls below the rest value and then slowly returns to zero.

7.   Conclusions

The numerical discretizations of the monodomain and bidomain models have been successfully carried out through a high-order DG formulation and the use of the Dubiner spectral basis. Furthermore, the space discretization has been complemented with a semi-implicit time integration. Finally, the numerical tests in Section 6 have from one side verified the convergence properties of the numerical solution and, on the other hand, compared the numerical simulation with the physiological propagation of the electric potential. In particular, the space discretization error decreases very fast as the grid is refined, meeting the expected theoretical orders of convergence, and soon becoming comparable with the time-step error. Also the physiological test, despite the anisotropy due to the mesh orientation, exhibits good reconstruction of both the depolarization and repolarization phases so that the electric wave propagation is correctly computed.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

PCA acknowledges the consortium iNEST (Interconnected North-East Innovation Ecosystem), Piano Nazionale di Ripresa e Resilienza (PNRR) – Missione 4 Componente 2, Investimento 1.5 – D.D. 1058 23/06/2022, ECS00000043, supported by the European Union's NextGenerationEU program and the INdAM - GNCS Project "Algoritmi efficienti per la gestione e adattazione di mesh poligonali", codice CUP_E53C22001930001. PFA has been partially funded by the research grant PRIN2020 n. 20204LN5N5 funded by MUR. PFA has been partially supported by ICSC–Centro Nazionale di Ricerca in High Performance Computing, BigData, and Quantum Computing funded by European Union—NextGeneration EU. PCA, PFA, CV are members of the INdAM group GNCS "Gruppo Nazionale per il Calcolo Scientifico" (National Group for Scientific Computing). CV has been partially supported by the Italian Ministry of University and Research (MIUR) within the PRIN (Research projects of relevant national interest) MIUR PRIN22-PNRR n. P20223KSS2 Machine learning for fluid-structure interaction in cardiovascular problems: efficient solutions, model reduction, inverse problems, and by the Italian Ministry of Health within the PNC PROGETTO HUB LIFE SCIENCE - DIAGNOSTICA AVANZATA (HLS-DA) "INNOVA", PNC-E3-2022-23683266-CUP: D43C22004930001, within the "Piano Nazionale Complementare Ecosistema Innovativo della Salute" - Codice univoco investimento: PNC-E3-2022-23683266

Conflict of interest

The authors declare no conflicts of interest.