Deriving the Bidomain Model of Cardiac Electrophysiology From a Cell-Based Model; Properties and Comparisons

The bidomain model is considered to be the gold standard for numerical simulation of the electrophysiology of cardiac tissue. The model provides important insights into the conduction properties of the electrochemical wave traversing the cardiac muscle in every heartbeat. However, in normal resolution, the model represents the average over a large number of cardiomyocytes, and more accurate models based on representations of all individual cells have therefore been introduced in order to gain insight into the conduction properties close to the myocytes. The more accurate model considered here is referred to as the EMI model since both the extracellular space (E), the cell membrane (M) and the intracellular space (I) are explicitly represented in the model. Here, we show that the bidomain model can be derived from the cell-based EMI model and we thus reveal the close relation between the two models, and obtain an indication of the error introduced in the approximation. Also, we present numerical simulations comparing the results of the two models and thereby highlight both similarities and differences between the models. We observe that the deviations between the solutions of the models become larger for larger cell sizes. Furthermore, we observe that the bidomain model provides solutions that are very similar to the EMI model when conductive properties of the tissue are in the normal range, but large deviations are present when the resistance between cardiomyocytes is increased.

1. Introduction

Mathematical models are indispensable for understanding the complex processes underlying cardiac electrophysiology. A wide variety of models have been developed for the key processes going on across the membrane of cardiomyocytes (see, e.g., Rudy and Silva, 2006; Rudy, 2012; Qu et al., 2014; Amuzescu et al., 2021), where the latter paper presents a comprehensive overview of the evolution of these models. The models of the membrane dynamics have also been extended to yield descriptions of the electrophysiological properties of cardiac tissue, commonly represented by the bidomain model or the somewhat simpler monodomain model (see Tung, 1978; Neu and Krassowska, 1993; Sundnes et al., 2007; Clayton and Panfilov, 2008; Vigmond et al., 2008; Linge et al., 2009; Niederer et al., 2011a; Franzone et al., 2014). The use of mathematical models for understanding the properties of the cardiac action potential (AP) across the membrane of cardiomyocytes is very widespread, and so is the use of the bidomain/monodomain models for understanding the properties of the excitation wave traversing cardiac tissue during each heartbeat. However, the spatial bidomain/monodomain models have two inherent limitations. The main limitation is that the extracellular space, the membrane of the myocyte, and the intracellular space are all assumed to be present everywhere. This assumption is indeed courageous but has provided surprisingly accurate results and presently underpins the understanding of cardiac conduction. The second limitation is that convergence is obtained using a relatively coarse mesh (Δx ~ 0.25 mm, see Xie et al., 2004; Clayton and Panfilov, 2008; Niederer et al., 2011b) and thus a typical mesh block contains several hundred cardiomyocytes (see, e.g., Jæger et al., 2021a,b). Therefore, understanding of the conduction properties (see, e.g., Henriquez, 2014; Veeraraghavan et al., 2014) close to the myocytes cannot be achieved using these models (see, e.g., Jæger et al., 2021a).

These limitations of the homogenized (bidomain/monodomain) models are well known and several authors have developed alternatives where all individual cells are explicitly represented in the models (see, e.g., Spach et al., 2007; Jacquemet and Henriquez, 2009; Hubbard and Henriquez, 2014; Lin and Keener, 2014; Tveito et al., 2017a; Weinberg, 2017; Jæger et al., 2019, 2021a; Domínguez et al., 2021; Jæger and Tveito, 2021). Here, we will apply the EMI model where both the extracellular space (E), the cell membrane (M) and the intracellular space (I) are explicitly represented in the model (see, e.g., Tveito et al., 2017a,b; Jæger and Tveito, 2021), and compare properties with the homogenized bidomain model. First, we will show how the bidomain model can be derived from the more accurate EMI model. Earlier derivations of the bidomain equations (see, e.g., Neu and Krassowska, 1993; Franzone et al., 2014; Henriquez and Ying, 2021) relies on homogenization of cardiac tissue, whereas the derivation given here follows directly from the EMI model. As part of this derivation, we can identify the main sources of deviations between the models.

Next, we will compare the properties of the bidomain model and the EMI model using numerical simulations. We first show that the deviations between the results obtained by the bidomain model and the EMI model become small as the cell size is reduced. This property is consistent with the error introduced in the derivation of the bidomain model. Secondly, we demonstrate that, for conduction properties providing a normal excitation wave with a conduction velocity of about 50 cm/s, the solutions of the EMI model and the bidomain model are very similar. However, as the resistance between the myocytes (through the gap junctions) is increased, the deviation between the solutions increases considerably.

It should be noted that the representation of all individual cardiomyocytes implies a significant increase in the computation load since the mesh resolution needs to be reduced from about Δx ~ 0.25 mm for a finite difference method (FDM) of the bidomain model to about δx ~ 10 μm for a finite element method (FEM) code solving the EMI model (see Jæger et al., 2021a,b). The number of mesh blocks is Δx³/δx³ = 15, 600 times larger for the EMI model than for the bidomain model, and, therefore, the computational load increases significantly when every myocyte in the tissue is resolved.

The choice of using either an averaged model like the bidomain model or a cell-based model like the EMI model, depends on the application under consideration. The bidomain model is very useful for simulating large scale problems, whereas EMI is better suited when the dynamics close to individual myocytes, or even inside individual myocytes, are of importance.

2. Methods

In this section we will derive the bidomain model commonly used to model the electrical activity of the heart from a more detailed model where each cell is represented. This cell-based model is referred to as the EMI model and is derived from Maxwell's equations of electromagnetism in Agudelo-Toro (2012) and Jæger and Tveito (2021). We will start by introducing the equations of the EMI model before we describe the derivation of the bidomain model from these equations. Finally, we discuss how the bidomain model parameters can be defined using the parameter values and tissue geometry of the EMI model.

2.1. The EMI Model

Consider a domain consisting of a single cell, Ω_i, surrounded by an extracellular space, Ω_e, with a cell membrane, Γ, separating the two spaces Ω_i and Ω_e. For such a domain, the electrical activity may be modeled by the EMI model (see, e.g., Roberts et al., 2008; Stinstra et al., 2010; Tveito et al., 2017a; Jæger et al., 2019), given by the equations

$\begin{array}{ll}\nabla ·{\sigma }_i\nabla u_i=0,\text{ in}{\Omega }_i& (1)\end{array}$

$\begin{array}{ll}\nabla ·{\sigma }_e\nabla u_e=0,\text{ in}{\Omega }_e& (2)\end{array}$

$\begin{array}{ll}{\text{n}}_e·{\sigma }_e\nabla u_e=-{\text{n}}_i·{\sigma }_i\nabla u_i\equiv I_m,\text{ at }\Gamma ,& (3)\end{array}$

$\begin{array}{ll}{u}_i-u_e=v\text{ at }\Gamma ,& (4)\end{array}$

$\begin{array}{ll}{I}_m=C_m\frac{\partial v}{\partial t}+I_{\text{ion}}\text{ at }\Gamma ,& (5)\end{array}$

$\begin{array}{ll}{u}_e=0\text{ at }\partial {\Omega }_e^D,& (6)\end{array}$

$\begin{array}{ll}\frac{\partial u_e}{\partial {\text{n}}_e}=0\text{ at }\partial {\Omega }_e^N.& (7)\end{array}$

Here, u_i, u_e, and v are the intracellular, extracellular and membrane potentials (in mV) defined in Ω_i, Ω_e and at Γ, respectively, n_i and n_e are the outward pointing normal vectors of the intracellular and extracellular spaces, respectively, C_m is the specific membrane capacitance (in μF/cm²), I_ion is the ionic current density across the membrane (in μA/cm²), I_m the sum of the capacitive and ionic current densities (in μA/cm²), and σ_i and σ_e are the intracellular and extracellular conductivities, respectively (in mS/cm). The Equations (6) and (7) are Dirichlet and Neumann boundary conditions, respectively, for the outer boundary of the extracellular space.

2.1.1. Extension to Cells Connected by Gap Junctions

To model collections of connected cardiac cells, e.g., like illustrated in Figure 1A, the EMI model for a single cell may be extended to include a model for the currents through the gap junctions connecting neighboring cells (see, e.g., Tveito et al., 2017a; Jæger and Tveito, 2021; Jæger et al., 2021c). For example, for two connected cells 1 and 2, the EMI model can be extended to include equations of the form

$\begin{array}{ll}{\text{n}}_i^2·{\sigma }_i\nabla u_i^2=-{\text{n}}_i^1·{\sigma }_i\nabla u_i^1\equiv I_{1,2},\text{ at }{\Gamma }_g,& (8)\end{array}$

$\begin{array}{ll}{u}_i^1-u_i^2=w,\text{ at }{\Gamma }_g,& (9)\end{array}$

$\begin{array}{ll}{I}_{1,2}=C_g\frac{\partial w}{\partial t}+I_{\text{gap}},\text{ at }{\Gamma }_g,& (10)\end{array}$

where Γ_g is the interface between the two cells (i.e., the intercalated disc). Furthermore, $u_i^1$ and $u_i^2$ are the intracellular potentials of the two cells, w is the potential difference between the two cells, and ${\text{n}}_i^1$, and ${\text{n}}_i^2$ are the outward pointing normal vectors of the cells. In addition, C_g is the specific capacitance of the intercalated discs (in μF/cm²), I_gap is the current density through the gap junction proteins located at the intercalated discs (in μA/cm²), and I_1,2 is the sum of the capacitive current density over the intercalated discs and the current density through the gap junction proteins connecting the two cells. The current density through the gap junction proteins, I_gap, is commonly modeled using the simple passive model

$\begin{array}{ll}{I}_{\text{gap}}=\frac{1}{R_g}w=G_{g}w,& (11)\end{array}$

where R_g is the specific resistance of the gap junctions (in kΩcm²) and G_g is the corresponding specific conductance (in mS/cm²). A further explanation of the coupling between two adjacent cells is given in section 1.2.4 of Jæger and Tveito (2021).

FIGURE 1

Figure 1. (A) Illustration of an EMI model domain for four connected cells. The intracellular space (orange) is denoted by Ω_i and the extracellular space (red) is denoted by Ω_e. The cell membrane is defined as the interface between the intracellular and extracellular spaces and is denoted by Γ. Similarly, the intercalated discs (purple) are denoted by Γ_g and are defined as the interface between neighboring cells. Each cell is shaped as a cylinder with a diameter increasing slightly toward the center of the cell. (B) Illustration of the finite element mesh used to represent a single cell in the EMI model simulations.

2.2. Derivation of the Bidomain Model From the EMI Model

Instead of using the detailed model (Equations 1–11), modeling of the electrical activity of cardiac tissue is usually performed using the homogenized bidomain and monodomain models. In these models, the detailed geometry of the individual cells and intercalated discs do not have to be represented in the computational mesh because the intracellular space, the extracellular space and the cell membrane are all assumed to exist everywhere in the tissue. We will now describe a possible derivation of the homogenized bidomain model from the EMI model equations described above. Note, however, that more rigorous versions of this derivation, using mathematical two-scale homogenization, have also been presented (see, e.g., Neu and Krassowska, 1993; Franzone et al., 2014; Henriquez and Ying, 2021).

2.2.1. Starting Point of the Derivation

Assume that we have a relatively large collection of cells, and consider a small volume, Δ, in this cell collection, as illustrated in Figure 2A. We assume that this volume contains a number of cells with an associated surrounding extracellular space, and that the EMI model equations apply in the extracellular domain, in the intracellular domain, at the cell membrane and at the intercalated discs in this small block of tissue.

FIGURE 2

Figure 2. Illustration of a tissue block, Δ, containing a number of cells (orange) and a surrounding extracellular space (red). In Step 1 of the derivation we approximate the discontinuous intracellular space (A) consisting of individual cells connected by gap junctions by a continuous intracellular space (B).

Step 1: Approximating the Intracellular Conductivity

As a first step in the derivation, we wish to approximate the intracellular conductivity to take both the purely intracellular space and the gap junctions between neighboring cells into account. In other words, we wish to reformulate the full EMI model (Equations 1–11) to a system of the form

$\begin{array}{ll}\nabla ·{\stackrel{-}{\sigma }}_i\nabla u_i=0,\text{ in}{\Omega }_i& (12)\end{array}$

$\begin{array}{ll}\nabla ·{\sigma }_e\nabla u_e=0,\text{ in}{\Omega }_e& (13)\end{array}$

$\begin{array}{ll}{\text{n}}_e·{\sigma }_e\nabla u_e=-{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_i\equiv I_m,\text{ at }\Gamma ,& (14)\end{array}$

$\begin{array}{ll}{u}_i-u_e=v\text{ at }\Gamma ,& (15)\end{array}$

$\begin{array}{ll}{I}_m=C_m\frac{\partial v}{\partial t}+I_{\text{ion}}\text{ at }\Gamma ,& (16)\end{array}$

$\begin{array}{ll}{u}_e=0\text{ at }\partial {\Omega }_e^D,& (17)\end{array}$

$\begin{array}{ll}\frac{\partial u_e}{\partial {\text{n}}_e}=0\text{ at }\partial {\Omega }_e^N,& (18)\end{array}$

where ${\stackrel{-}{\sigma }}_i$ is the average conductivity of the intracellular space including gap junctions. We wish to express ${\stackrel{-}{\sigma }}_i$ such that the total intracellular resistance of a tissue block is close to the total intracellular resistance of the tissue block when the full EMI model (Equations 1–11) applies. Such an expression for ${\stackrel{-}{\sigma }}_i$ is derived below in section 2.3.3.

In the remaining part of the derivation we will treat the intracellular domain as a continuous domain (see Figure 2B), and assume that the simplified EMI system (Equations 12–18) applies.

Step 2: Applying the Divergence Theorem

In the next step of the derivation, we consider the purely intracellular part of the tissue block Δ and apply the divergence theorem for ${\stackrel{-}{\sigma }}_i\nabla u_i$ to obtain,

$\begin{array}{ll}{\displaystyle {\int }_{\partial {\Omega }_i^{\Delta }}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS={\displaystyle {\int }_{{\Omega }_i^{\Delta }}}\nabla ·{\stackrel{-}{\sigma }}_i\nabla u_{i}dV,& (19)\end{array}$

where ${\Omega }_i^{\Delta }$ is the intracellular space contained in the tissue block and $\partial {\Omega }_i^{\Delta }$ is the boundary of the intracellular space contained in the tissue block. This boundary can be separated into the boundary between the intracellular space and the extracellular space contained in the tissue block, i.e. the cell membrane, and the intracellular part of the outer boundary of the tissue block in each spatial direction. Rewriting the surface integral, we obtain

$\begin{array}{l}{\displaystyle {\int }_{A_i^{x,+}}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS+{\displaystyle {\int }_{A_i^{x,-}}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS\\ +{\displaystyle {\int }_{A_i^{y,+}}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS+{\displaystyle {\int }_{A_i^{y,-}}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS\\ +{\displaystyle {\int }_{A_i^{z,+}}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS+{\displaystyle {\int }_{A_i^{z,-}}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS\\ +{\displaystyle {\int }_{{\Gamma }_{\Delta }}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS={\displaystyle {\int }_{{\Omega }_i^{\Delta }}}\nabla ·{\stackrel{-}{\sigma }}_i\nabla u_{i}dV,& (20)\end{array}$

where $A_i^{x,+}$ is the intracellular part of the boundary of the tissue block in the positive x-direction, $A_i^{x,-}$ is the intracellular part of the boundary of the tissue block in the negative x-direction, and the surfaces $A_i^{y,+}$, $A_i^{y,-}$, $A_i^{z,+}$, and $A_i^{z,-}$ are defined similarly for the intracellular part of the boundaries of the tissue block in the y- and z-directions. Furthermore, Γ_Δ is the membrane contained in the tissue block.

By applying the divergence theorem and similar definitions for the extracellular space, we likewise obtain

$\begin{array}{l}{\displaystyle {\int }_{A_e^{x,+}}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS+{\displaystyle {\int }_{A_e^{x,-}}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS\\ +{\displaystyle {\int }_{A_e^{y,+}}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS+{\displaystyle {\int }_{A_e^{y,-}}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS\\ +{\displaystyle {\int }_{A_e^{z,+}}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS+{\displaystyle {\int }_{A_e^{z,-}}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS\\ +{\displaystyle {\int }_{{\Gamma }_{\Delta }}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS={\displaystyle {\int }_{{\Omega }_e^{\Delta }}}\nabla ·{\sigma }_e\nabla u_{e}dV.& (21)\end{array}$

Step 3: Applying the EMI Model Equations (12–(14)

By inserting Equations (12) and (13) into Equations (20) and (21), we find that the right hand sides of Equations (20) and (21) are zero. Moreover, by inserting Equation (14), we get

$\begin{array}{l}{\displaystyle {\int }_{A_i^{x,+}}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS+{\displaystyle {\int }_{A_i^{x,-}}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS\\ +{\displaystyle {\int }_{A_i^{y,+}}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS+{\displaystyle {\int }_{A_i^{y,-}}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS\\ +{\displaystyle {\int }_{A_i^{z,+}}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS+{\displaystyle {\int }_{A_i^{z,-}}}{\text{n}}_i·{\stackrel{-}{\sigma }}_i\nabla u_{i}dS-{\displaystyle {\int }_{{\Gamma }_{\Delta }}}{I}_{m}dS=0& (22)\end{array}$

for the intracellular part, and

$\begin{array}{l}{\displaystyle {\int }_{A_e^{x,+}}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS+{\displaystyle {\int }_{A_e^{x,-}}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS\\ +{\displaystyle {\int }_{A_e^{y,+}}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS+{\displaystyle {\int }_{A_e^{y,-}}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS\\ +{\displaystyle {\int }_{A_e^{z,+}}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS+{\displaystyle {\int }_{A_e^{z,-}}}{\text{n}}_e·{\sigma }_e\nabla u_{e}dS+{\displaystyle {\int }_{{\Gamma }_{\Delta }}}{I}_{m}dS=0& (23)\end{array}$

for the extracellular part.

Step 4: Extending the Variables and Parameters to Be Defined Everywhere

In order to avoid having to represent the detailed geometry of the cell tissue, we now define some new variables U_i, U_e, and V that each are defined in the entire domain Ω = Ω_i∪Ω_e, and thus also in the entire tissue block, Δ. We want these variables to fulfill the integral conditions specified in Equations (22) and (23). In addition, we assume that the definitions of the membrane potential and I_m specified in Equations (15) and (16) apply in the entire domain. In other words, in an arbitrary tissue block, Δ, of Ω, we seek solutions U_i, U_e, and V such that

$\begin{array}{l}{\displaystyle {\int }_{A_i^{x,+}}}\text{n}·{\stackrel{-}{\sigma }}_i\nabla U_{i}dS+{\displaystyle {\int }_{A_i^{x,-}}}\text{n}·{\stackrel{-}{\sigma }}_i\nabla U_{i}dS\\ +{\displaystyle {\int }_{A_i^{y,+}}}\text{n}·{\stackrel{-}{\sigma }}_i\nabla U_{i}dS+{\displaystyle {\int }_{A_i^{y,-}}}\text{n}·{\stackrel{-}{\sigma }}_i\nabla U_{i}dS\\ +{\displaystyle {\int }_{A_i^{z,+}}}\text{n}·{\stackrel{-}{\sigma }}_i\nabla U_{i}dS+{\displaystyle {\int }_{A_i^{z,-}}}\text{n}·{\stackrel{-}{\sigma }}_i\nabla U_{i}dS-{\displaystyle {\int }_{{\Gamma }_{\Delta }}}{I}_{m}dS=0,& (24)\end{array}$

$\begin{array}{l}{\displaystyle {\int }_{A_e^{x,+}}}\text{n}·{\sigma }_e\nabla U_{e}dS+{\displaystyle {\int }_{A_e^{x,-}}}\text{n}·{\sigma }_e\nabla U_{e}dS\\ +{\displaystyle {\int }_{A_e^{y,+}}}\text{n}·{\sigma }_e\nabla U_{e}dS+{\displaystyle {\int }_{A_e^{y,-}}}\text{n}·{\sigma }_e\nabla U_{e}dS\\ +{\displaystyle {\int }_{A_e^{z,+}}}\text{n}·{\sigma }_e\nabla U_{e}dS+{\displaystyle {\int }_{A_e^{z,-}}}\text{n}·{\sigma }_e\nabla U_{e}dS+{\displaystyle {\int }_{{\Gamma }_{\Delta }}}{I}_{m}dS=0,& (25)\end{array}$

$\begin{array}{ll}V=U_i-U_e,& (26)\end{array}$

$\begin{array}{ll}{I}_m=C_m{V}_t+I_{\text{ion}}.& (27)\end{array}$

Here, n is the outward pointing normal vector of the tissue block, and ${\stackrel{-}{\sigma }}_i$, σ_e, C_m, I_m and I_ion have been extended to be defined in the entire domain.

Step 5: Approximate the Surface Integrals

Since I_m is now defined in the entire tissue block, and not just on the membrane, the surface integral over the membrane can be approximated by

$\begin{array}{ll}{\displaystyle {\int }_{{\Gamma }_{\Delta }}}{I}_{m}dS\approx {\displaystyle {\int }_{\Delta }}\chi I_{m}dV,& (28)\end{array}$

where Δ is the entire tissue block and χ is the membrane surface to volume ratio, i.e., the surface area of the membrane contained in Δ divided by the volume of Δ.

In addition, the integrals in Equations (24) and (25) over the outer boundary of the tissue block is separated into the intracellular and extracellular parts of the tissue block, and in this step of the derivation, we wish to approximate these integrals to be defined over the entire tissue block boundaries. In order to do this, we apply the approximation

$\begin{array}{ll}{\displaystyle {\int }_{A_i^{x,+}}}\text{n}·{\stackrel{-}{\sigma }}_i\nabla U_{i}dS\approx {Ā}_i^x{\displaystyle {\int }_{A^{x,+}}}\text{n}·{\stackrel{-}{\sigma }}_i\nabla U_{i}dS,& (29)\end{array}$

and similar approximations for the remaining surfaces. Here, ${Ā}_i^x$ is the average fraction of the cross-sectional area of the tissue block perpendicular to the x-direction that is occupied by the intracellular space and A^x,+ is the entire boundary of the tissue block in the positive x-direction. Inserting this type of approximation in all the integrals over the outer boundaries of the tissue block, Equations (24) and (25) can be approximated as

$\begin{array}{l}{\displaystyle {\int }_{A^{x,+}}}\text{n}·{Ā}_i^x{\stackrel{-}{\sigma }}_i\nabla U_{i}dS+{\displaystyle {\int }_{A^{x,-}}}\text{n}·{Ā}_i^x{\stackrel{-}{\sigma }}_i\nabla U_{i}dS\\ +{\displaystyle {\int }_{A^{y,+}}}\text{n}·{Ā}_i^y{\stackrel{-}{\sigma }}_i\nabla U_{i}dS+{\displaystyle {\int }_{A^{y,-}}}\text{n}·{Ā}_i^y{\stackrel{-}{\sigma }}_i\nabla U_{i}dS\\ +{\displaystyle {\int }_{A^{z,+}}}\text{n}·{Ā}_i^z{\stackrel{-}{\sigma }}_i\nabla U_{i}dS+{\displaystyle {\int }_{A^{z,-}}}\text{n}·{Ā}_i^z{\stackrel{-}{\sigma }}_i\nabla U_{i}dS-{\displaystyle {\int }_{\Delta }}\chi I_{m}dV=0,& (30)\end{array}$

$\begin{array}{l}{\displaystyle {\int }_{A^{x,+}}}\text{n}·{Ā}_e^x{\sigma }_e\nabla U_{e}dS+{\displaystyle {\int }_{A^{x,-}}}\text{n}·{Ā}_e^x{\sigma }_e\nabla U_{e}dS\\ +{\displaystyle {\int }_{A^{y,+}}}\text{n}·{Ā}_e^y{\sigma }_e\nabla U_{e}dS+{\displaystyle {\int }_{A^{y,-}}}\text{n}·{Ā}_e^y{\sigma }_e\nabla U_{e}dS\\ +{\displaystyle {\int }_{A^{z,+}}}\text{n}·{Ā}_e^z{\sigma }_e\nabla U_{e}dS+{\displaystyle {\int }_{A^{z,-}}}\text{n}·{Ā}_e^z{\sigma }_e\nabla U_{e}dS+{\displaystyle {\int }_{\Delta }}\chi I_{m}dV=0.& (31)\end{array}$

Furthermore, we may define a set of scaled bidomain conductivities,

$\begin{array}{ll}{\tilde{\sigma }}_i^x={Ā}_i^x{\stackrel{-}{\sigma }}_i^x,{\tilde{\sigma }}_i^y={Ā}_i^y{\stackrel{-}{\sigma }}_i^y,{\tilde{\sigma }}_i^z={Ā}_i^z{\stackrel{-}{\sigma }}_i^z,& (32)\end{array}$

$\begin{array}{ll}{\tilde{\sigma }}_e^x={Ā}_e^x{\sigma }_e^x,{\tilde{\sigma }}_e^y={Ā}_e^y{\sigma }_e^y,{\tilde{\sigma }}_e^z={Ā}_e^z{\sigma }_e^z,& (33)\end{array}$

where ${\stackrel{-}{\sigma }}_i^x$, ${\stackrel{-}{\sigma }}_i^y$ and ${\stackrel{-}{\sigma }}_i^z$ and ${\sigma }_e^x$, ${\sigma }_e^y$ and ${\sigma }_e^z$ refer to the possible directional dependence of ${\stackrel{-}{\sigma }}_i$ and σ_e. We also define the associated bidomain conductivity tensors,

$\begin{array}{ll}{M}_i=\left(\begin{array}{ccc}{\tilde{\sigma }}_i^x& 0& 0\\ 0& {\tilde{\sigma }}_i^y& 0\\ 0& 0& {\tilde{\sigma }}_i^z\end{array}\right),M_e=\left(\begin{array}{ccc}{\tilde{\sigma }}_e^x& 0& 0\\ 0& {\tilde{\sigma }}_e^y& 0\\ 0& 0& {\tilde{\sigma }}_e^z\end{array}\right).& (34)\end{array}$

By introducing these tensors, Equations (30) and (31) can be rewritten as

$\begin{array}{ll}{\displaystyle {\int }_{\partial \Delta }}\text{n}·(M_i\nabla U_i)dS-{\displaystyle {\int }_{\Delta }}\chi I_{m}dV=0,& (35)\end{array}$

$\begin{array}{ll}{\displaystyle {\int }_{\partial \Delta }}\text{n}·(M_e\nabla U_e)dS+{\displaystyle {\int }_{\Delta }}\chi I_{m}dV=0,& (36)\end{array}$

where ∂Δ represents the entire outer surface of the tissue block.

Step 6: Reapply the Divergence Theorem for the New Variables

We may now reapply the divergence theorem for the newly defined variables U_i and U_e defined in the entire tissue block. This yields

$\begin{array}{ll}{\displaystyle {\int }_{\Delta }}\nabla ·(M_i\nabla U_i)dV-{\displaystyle {\int }_{\Delta }}\chi I_{m}dV=0,& (37)\end{array}$

$\begin{array}{ll}{\displaystyle {\int }_{\Delta }}\nabla ·(M_e\nabla U_e)dV+{\displaystyle {\int }_{\Delta }}\chi I_{m}dV=0.& (38)\end{array}$

We also note that the volume Δ was chosen arbitrarily. Therefore, the more general relation

$\begin{array}{ll}\nabla ·(M_i\nabla U_i)-\chi I_m=0,& (39)\end{array}$

$\begin{array}{ll}\nabla ·(M_e\nabla U_e)+\chi I_m=0,& (40)\end{array}$

holds.

Step 7: Rearranging the Terms and Inserting Equations (26) and (27)

By rearranging Equation (39) and adding Equations (39) and (40), these equations may be rewritten as

$\begin{array}{ll}\nabla ·(M_i\nabla U_i)=\chi I_m,& (41)\end{array}$

$\begin{array}{ll}\nabla ·(M_i\nabla U_i)+\nabla ·(M_e\nabla U_e)=0.& (42)\end{array}$

Finally, by inserting Equations (26) and (27), we obtain the bidomain model equations

$\begin{array}{ll}\nabla ·(M_i\nabla V)+\nabla ·(M_i\nabla U_e)=\chi \left(C_m\frac{\partial V}{\partial t}+I_{\text{ion}}\right),& (43)\end{array}$

$\begin{array}{ll}\nabla ·(M_i\nabla V)+\nabla ·\left((M_i+M_e)\nabla U_e\right)=0,& (44)\end{array}$

where we recall that M_i and M_e are intracellular and extracellular conductivity tensors (in mS/cm) defined in Equations (32)–(34), χ is the membrane surface to volume ratio (in cm⁻¹), C_m is the specific membrane capacitance (in μF/cm²), I_ion is the current density through ion channels, pumps and exchangers on the cell membrane (in μA/cm²) and V and U_e (in mV) are the bidomain model membrane and extracellular potentials, respectively, defined in the entire domain. Furthermore, the intracellular potential (in mV) may be computed by

$\begin{array}{ll}{U}_i=V+U_e.& (45)\end{array}$

In addition, the boundary conditions

$\begin{array}{ll}{U}_e=0\text{ at }\partial {\Omega }^D,& (46)\end{array}$

$\begin{array}{ll}\frac{\partial U_e}{\partial \text{n}}=0\text{ at }\partial {\Omega }^N,& (47)\end{array}$

are assumed to hold at the boundary of the domain where Ω^D coincides with the EMI model boundary ${\Omega }_e^D$, and Ω^N coincides with the EMI model boundary ${\Omega }_e^N.$

2.3. Expressions for the Bidomain Model Parameters

The bidomain model as derived above introduces a set of new parameters, namely the conductivity tensors, M_i and M_e, and the surface to volume ratio, χ. Considering their definitions, values for these parameters may be derived from the geometry and parameters of the EMI model. In this subsection, we suggest an approach for making these definitions by considering an EMI model mesh of a volume Ω, containing an intracellular volume, Ω_i, and extracellular volume, Ω_e, a surface for the cell membranes, Γ, and a collection of surfaces for the intercalated discs, Γ_g. For simplicity, we assume that the value of all the EMI model parameters and the tissue geometry do not vary in different parts of the domain, so that the bidomain model parameters can be treated as constants throughout the domain. In addition, we assume that the total domain Ω = Ω_i∪Ω_e is shaped as a rectangular cuboid with lengths L_x, L_y and L_y in the x-, y- and z-directions, respectively. An alternative approach for setting up the bidomain model conductivities from the EMI model parameters and a simplified tissue geometry is presented in Henriquez and Ying (2021).

2.3.1. Surface to Volume Ratio, χ

In order to compute the surface to volume ratio from an EMI model mesh, we may simply compute

$\begin{array}{ll}{A}_{\Gamma ,\Omega }={\displaystyle {\int }_{\Gamma }}1dS,& (48)\end{array}$

$\begin{array}{ll}{V}_{\Omega }={\displaystyle {\int }_{\Omega }}1dS,& (49)\end{array}$

where A_Γ,Ω represents the total membrane area in the domain and V_Ω represents the volume of the domain. Assuming an even distribution of cells throughout the domain, the surface to volume ratio can then be defined as

$\begin{array}{ll}\chi =\frac{A_{\Gamma ,\Omega }}{V_{\Omega }}.& (50)\end{array}$

2.3.2. Average Cross-Sectional Area Fractions

We first consider the average intracellular fraction of the cross-sectional area perpendicular to the x-axis, ${Ā}_i^x$. Let A^x(x) be the cross-sectional area of Ω perpendicular to the x-axis, and let $A_i^x(x)$ be the fraction belonging to Ω_i. Then

$\begin{array}{ll}{V}_{{\Omega }_i}={\displaystyle {\int }_0^{L_x}}{A}_i^x(x)A^x(x)dx={Ā}_i^x{\displaystyle {\int }_0^{L_x}}{A}^x(x)dx={Ā}_i^x{V}_{\Omega }.& (51)\end{array}$

Hence,

$\begin{array}{ll}{Ā}_i^x=\frac{V_{{\Omega }_i}}{V_{\Omega }}.& (52)\end{array}$

Similar arguments yield

$\begin{array}{ll}{Ā}_i^x={Ā}_i^y={Ā}_i^z=\frac{V_{{\Omega }_i}}{V_{\Omega }},& (53)\end{array}$

$\begin{array}{ll}{Ā}_e^x={Ā}_e^y={Ā}_e^z=\frac{V_{{\Omega }_e}}{V_{\Omega }}=\left(1-\frac{V_{{\Omega }_i}}{V_{\Omega }}\right).& (54)\end{array}$

Here, V_Ωi can be computed from the EMI model mesh as

$\begin{array}{ll}{V}_{{\Omega }_i}={\displaystyle {\int }_{{\Omega }_i}}1dS.& (55)\end{array}$

2.3.3. Average Intracellular Conductivity

As described above, we wish to define an average conductivity ${\stackrel{-}{\sigma }}_i$ such that the simplified EMI model (Equations 12–18) is a good approximation of the full EMI model (Equations 1–11). In particular, we wish to find a ${\stackrel{-}{\sigma }}_i$ such that the total intracellular resistance of the simplified model is close to the total intracellular resistance of the full model. To simplify this argument, we assume that there is no capacitive current across the intercalated discs, i.e. that the current between two cells is given by I_gap (see Equation 11).

We start by considering the total resistance in the x-direction of the domain. In the full EMI model (Equations 1–11) with the capacitive current set to zero, this is given by the sum of the resistance over the purely intracellular space ($R_c^x$) and the resistance over the gap junctions ($R_j^x$) (Shaw and Rudy, 1997):

$\begin{array}{ll}{R}_i^x=R_c^x+R_j^x.& (56)\end{array}$

The total resistance in the purely intracellular space is given by (Plonsey and Barr, 2007)

$\begin{array}{ll}{R}_c^x=\frac{L_x}{{\sigma }_i{Ā}_i^x{L}_y{L}_z},& (57)\end{array}$

where ${Ā}_i^x{L}_y{L}_z$ is the average intracellular cross-sectional area of the domain perpendicular to the x-direction. Assuming that the cells are organized as a regular grid in the x-, y- and z-directions, the total resistance through gap junctions in the x-direction is given by

$\begin{array}{ll}{R}_j^x=\frac{(N_x-1)R_g}{N_y{N}_z{A}_j^x},& (58)\end{array}$

where R_g is the specific gap junction resistance (in kΩcm²), as it appears in the full EMI model, $A_j^x$ is the area of a single intercalated disc perpendicular to the x-direction and N_x, N_y, and N_z are the number of cells in the x-, y-, and z-directions, respectively. Thus, N_x − 1 is the number of intercalated disc collections along the length of the domain in the x-direction, N_yN_z is the number of intercalated disc for each such collection, and $N_y{N}_z{A}_j^x$ is the total cross-sectional area of each of the intercalated disc collections.

In the simplified model (Equations 12–18), the total resistance is given by (Plonsey and Barr, 2007)

$\begin{array}{ll}{R}_i^x=\frac{L_x}{{\stackrel{-}{\sigma }}_i^x{Ā}_i^x{L}_y{L}_z}.& (59)\end{array}$

Therefore, in order for the total resistance to be the same in the two formulations, we wish ${\stackrel{-}{\sigma }}_i^x$ to satisfy

$\begin{array}{ll}\frac{L_x}{{\stackrel{-}{\sigma }}_i^x{Ā}_i^x{L}_y{L}_z}=\frac{L_x}{{\sigma }_i{Ā}_i^x{L}_y{L}_z}+\frac{(N_x-1)R_g}{N_y{N}_z{A}_j^x},& (60)\end{array}$

which yields

$\begin{array}{ll}{\stackrel{-}{{\sigma }_i}}^x=\frac{{\sigma }_i}{1+\frac{{\sigma }_i{R}_g(N_x-1){Ā}_i^x{L}_y{L}_z}{L_x{N}_y{N}_z{A}_j^x}}.& (61)\end{array}$

From the EMI model mesh, we may compute

$\begin{array}{ll}{A}_{j,\Omega }^x={\displaystyle {\int }_{{\Gamma }_g^x}}1dS,& (62)\end{array}$

as the total area of all intercalated discs perpendicular to the x-direction, ${\Gamma }_g^x$. Since $A_j^x$ is defined as the area of a single intercalated disc perpendicular to the x-direction, we note that

$\begin{array}{ll}{A}_j^x=\frac{A_{j,\Omega }^x}{(N_x-1)N_y{N}_z},& (63)\end{array}$

where (N_x − 1)N_yN_z is the total number of intercalated discs in the x-direction. Inserting Equations (63) into Equation (61) yields

$\begin{array}{ll}{\stackrel{-}{\sigma }}_i^x=\frac{{\sigma }_i}{1+\frac{{(N_x-1)}^2{Ā}_i^x{L}_y{L}_z{\sigma }_i{R}_g}{L_x{A}_{j,\Omega }^x}}.& (64)\end{array}$

We also note that from Equation (52), we have that

$\begin{array}{ll}{Ā}_i^x=\frac{V_{{\Omega }_i}}{V_{\Omega }}=\frac{V_{{\Omega }_i}}{L_x{L}_y{L}_z}\Rightarrow {Ā}_i^x{L}_y{L}_z=\frac{V_{{\Omega }_i}}{L_x},& (65)\end{array}$

and inserting this into Equation (64), we obtain

$\begin{array}{ll}{\stackrel{-}{\sigma }}_i^x=\frac{{\sigma }_i}{1+\frac{{\sigma }_i{R}_g{V}_{{\Omega }_i}}{{\delta }_x^2{A}_{j,\Omega }^x}}.& (66)\end{array}$

where δ_x = L_x/(N_x − 1). Similar arguments for the y- and z-directions result in

$\begin{array}{ll}{\stackrel{-}{\sigma }}_i^y=\frac{{\sigma }_i}{1+\frac{{\sigma }_i{R}_g{V}_{{\Omega }_i}}{{\delta }_y^2{A}_{j,\Omega }^y}},& (67)\end{array}$

$\begin{array}{ll}{\stackrel{-}{\sigma }}_i^z=\frac{{\sigma }_i}{1+\frac{{\sigma }_i{R}_g{V}_{{\Omega }_i}}{{\delta }_z^2{A}_{j,\Omega }^z}}.& (68)\end{array}$

Since

$\begin{array}{ll}\frac{{\stackrel{-}{\sigma }}_i^x}{{\stackrel{-}{\sigma }}_i^y}=\frac{1+\frac{{\sigma }_i{R}_g{V}_{{\Omega }_i}}{{\delta }_y^2{A}_{j,\Omega }^y}}{1+\frac{{\sigma }_i{R}_g{V}_{{\Omega }_i}}{{\delta }_x^2{A}_{j,\Omega }^x}},& (69)\end{array}$

the anisotropy is governed by the difference between ${\delta }_x^2{A}_{j,\Omega }^x$ and ${\delta }_y^2{A}_{j,\Omega }^y,$ and similar for the other combination of axes.

2.3.4. Intracellular Conductivity Tensor

Inserting Equations (52) and (66)–(68) into Equation (32), we get

$\begin{array}{ll}{\tilde{\sigma }}_i^x={Ā}_i^x{\stackrel{-}{\sigma }}_i^x=\frac{V_{{\Omega }_i}}{V_{\Omega }}\frac{{\sigma }_i}{1+\frac{{\sigma }_i{(N_x-1)}^2{R}_g{V}_{{\Omega }_i}}{L_x^2{A}_{j,\Omega }^x}},& (70)\end{array}$

$\begin{array}{ll}{\tilde{\sigma }}_i^y={Ā}_i^y{\stackrel{-}{\sigma }}_i^y=\frac{V_{{\Omega }_i}}{V_{\Omega }}\frac{{\sigma }_i}{1+\frac{{\sigma }_i{(N_y-1)}^2{R}_g{V}_{{\Omega }_i}}{L_y^2{A}_{j,\Omega }^y}},& (71)\end{array}$

$\begin{array}{ll}{\tilde{\sigma }}_i^z={Ā}_i^z{\stackrel{-}{\sigma }}_i^z=\frac{V_{{\Omega }_i}}{V_{\Omega }}\frac{{\sigma }_i}{1+\frac{{\sigma }_i{(N_z-1)}^2{R}_g{V}_{{\Omega }_i}}{L_z^2{A}_{j,\Omega }^z}}.& (72)\end{array}$

2.3.5. Extracellular Conductivity Tensor

The extracellular conductivity tensor can be found directly from the cross-sectional area fractions and we get

$\begin{array}{ll}{\tilde{\sigma }}_e^x={\tilde{\sigma }}_e^y={\tilde{\sigma }}_e^z=\frac{V_{{\Omega }_e}}{V_{\Omega }}{\sigma }_e=\left(1-\frac{V_{{\Omega }_i}}{V_{\Omega }}\right){\sigma }_e.& (73)\end{array}$

3. Results

In order to compare the EMI model with the homogenized bidomain model, we set up a few example applications and perform numerical simulations of the two models. Note here that all EMI model simulations are performed in three dimensions (3D), whereas the bidomain model simulations are performed in two dimensions (2D) or one dimension (1D).

3.1. Simulation Set-Up

In our numerical simulations of the EMI model, we consider collections of cells shaped as cylinders with a slightly varying diameter. In all simulations, except for the ones where the cell length is varied and is explicitly specified, each cell is 120 μm long (in the x-direction) and has a radius varying from 6 μm at the cell ends to 7 μm at the center of the cell (see Figure 1). We let the distance from the boundary of the extracellular space to the cell collection be 2 μm in all spatial directions. The parameter values used in the simulations are specified in Table 1. The parameters used in the bidomain model are computed from the EMI model parameters and mesh as described in section 2.3.

TABLE 1

Table 1. Default parameter values used in the simulations.

All EMI model simulations are performed in 3D. However, for our example test cases with a 1D strand of cells and a 2D grid of cells, we use 1D and 2D versions, respectively, of the bidomain model. In the simulations of a 1D strand of cells, we apply homogenous Neumann boundary conditions on the outer boundary of the extracellular domain in the y- and z-directions and homogenous Dirichlet boundary conditions in the x-direction. In the simulations of a 2D grid of cells, we apply homogenous Neumann boundary conditions on the outer boundary of the extracellular domain in the z-direction and homogenous Dirichlet boundary conditions in the x- and y-directions.

3.2. Numerical Methods

The EMI model simulations are performed using the operator splitting procedure described in Jæger et al. (2021c,d), the numerical methods applied to (Jæger et al., 2021d) and the MFEM C++ finite element method library (Anderson et al., 2020; MFEM, 2021). For details on the numerical methods applied to solve the EMI model, we refer to Jæger et al. (2021a,c,d). The bidomain model simulations are performed in Matlab using a first-order temporal operator splitting procedure as described in Sundnes et al. (2006), where the ordinary differential part of the equations is solved using forward Euler and the partial differential part of the equations is solved using an implicit finite difference scheme. Unless otherwise specified, we use a time step of Δt = 0.001 ms in the simulations of both models. In the bidomain model simulations, we use a spatial discretization of Δx = Δy = 10 μm, roughly matching the typical edge length in the applied EMI model finite element mesh.

3.3. 1D Strand of Cells With a Passive Membrane Model

We first consider an example with a 1D strand of cells connected in the longitudinal direction (x-direction). The total length of the cell strand is 2 mm and we consider a number of different choices for the length of a single cell (and the associated total number of cells). In addition, we vary the value of the specific gap junction resistance, R_g. The membrane dynamics, I_ion, is modeled by a simple passive membrane model

$\begin{array}{ll}{I}_{\text{ion}}=\frac{1}{R_m}(v-v_0),& (74)\end{array}$

where R_m = 5 kΩcm² is the specific membrane resistance and v₀ = −80 mV is the resting membrane potential. We stimulate the first (leftmost) 400 μm of membrane in the x-direction by a constant stimulus current of size −10 μA/cm². Figure 3 shows the intracellular potential at time t = 20 ms along a line in the x-direction in the center of the domain for the EMI model and the associated solution of the bidomain model for a few combinations of cell length and R_g values. We observe that for small cells (lower panels), the solution of the bidomain model is in very good agreement with the results of the EMI model. However, if the cell size is increased, and the gap junction resistance is increased, there is a significant difference between the results of the EMI model and the bidomain model.

FIGURE 3

Figure 3. Intracellular potential, u_i, at time t = 20 ms in EMI model and bidomain model simulations using a passive membrane model (Equation 74) and the default parameter values specified in Table 1, except that the value of R_g is increased by the factor indicated by the column titles. In addition, the cell length (L_cell) is varied as described for each row of plots.

3.4. 1D Strand of Cells With an Active Membrane Model

Next, we consider an example with a 1D strand of 20 cells of length 120 μm with an active membrane model, modeled by the human left atrial base model from Jæger et al. (2021b). We initiate a traveling wave by stimulating the first 360 μm of cell membrane in the x-direction (corresponding to three cardiomyocytes) by a 1 ms long constant stimulus current of size −40 μA/cm². We measure the conduction velocity as the distance between a point a in the center of the domain in the x-direction and a point b located at 4/5 of the total domain length, divided by the difference in time between when the membrane potential in these two points reach a value above −20 mV. Using the default parameter values specified in Table 1, we get a conduction velocity of 50.8 cm/s in the bidomain model simulation. This is close to the value found in the EMI model simulation, which is 53.3 cm/s.

In Figure 4, we further investigate the relationship between the conduction velocity found in the bidomain and EMI model simulations when R_g is increased, representing reduced cell coupling. We consider two different discretization resolutions, the default resolution of Δt = 0.001 ms and Δx ~ 10 μm and a refined resolution of Δt = 0.0005 ms and Δx ~ 5 μm. We observe that for values of R_g relatively close to the default value, the conduction velocities found in simulations of the bidomain and EMI models are very similar. However, when R_g is considerably increased, the difference between the two model formulations appears to be more significant and the conduction velocity is considerably higher in the bidomain model simulations than in the corresponding EMI model simulations. Furthermore, we observe that for the EMI model, conduction is blocked when R_g is increased by a factor larger than about 2,000, whereas for the bidomain model, R_g can be increased by a factor of about 20,000 before conduction is blocked for the default resolution and conduction is not blocked for the considered values of R_g for the refined resolution. In addition, we note that the simulations of refined resolution appears to give very similar conduction velocities as for the default resolution in the EMI model simulations. For the bidomain model simulations, the two resolutions give very similar results for the first range of R_g values, but as the R_g value is severely increased, we can observe a difference between the two resolutions.

FIGURE 4

Figure 4. Conduction velocity as R_g is increased in EMI model and bidomain model simulations of a strand of 20 connected cells with an active membrane model (Jæger et al., 2021b). The values on the x-axis represent the factor with which the default R_g value in Table 1 is multiplied. The remaining parameter values are as specified in Table 1. Note that the plot is separated into three panels in order to improve the visibility of the data. Note also that we consider two different discretization resolutions in the simulations of each model. In the default case, Δt = 0.001 ms and Δx in the bidomain model and the typical edge length in the EMI model is 10 μm, whereas in the refined case (dashed lines), both the spatial and temporal discretization steps are reduced to half of the default values.

3.5. 2D Grid of Cells With an Active Membrane Model

Next, we consider a case of a grid of 25×25 connected cells with the same active membrane model as for the 1D strand simulations. We stimulate the membrane of an area corresponding to the 5×5 cells in the lower left corner by the same stimulation current as in the 1D case. Figure 5 shows the membrane potential, v, and the extracellular potential, u_e, from the EMI model and bidomain model simulations using the default parameter values specified in Table 1 at time t = 5 ms. The solution of the two models appears to be very similar. However, in Figure 6, we have performed a similar simulation where R_g is increased by a factor of 200. We consider the solution at t = 20 ms and observe that the traveling excitation wave has clearly traveled faster and reached further in the bidomain model simulation than in the EMI model simulation, consistent with the results of Figure 4.

FIGURE 5

Figure 5. Membrane potential, v, and extracellular potential, u_e, at time t = 5 ms in EMI model and bidomain model simulations using the default parameter values specified in Table 1 and an active membrane model (Jæger et al., 2021b). This simulation required a CPU time of 132 min for the EMI model and 2 min for the bidomain model.

FIGURE 6

Figure 6. Membrane potential, v, and extracellular potential, u_e, at time t = 20 ms in EMI model and bidomain model simulations using and active membrane model (Jæger et al., 2021b) and the default parameter values specified in Table 1, except that the value of R_g is increased by a factor of 200.

4. Discussion

The bidomain model continues to provide essential insights into cardiac conduction and how the electrochemical dynamics of the heart is affected by blocking ion channels (see, e.g., Zemzemi et al., 2013; Sharifi, 2017), increasing gap junction resistance (see, e.g., Roth, 1988; Bruce et al., 2014), introducing ischemia (see, e.g., Stinstra et al., 2004, 2005; Heidenreich et al., 2012) or performing defibrillation (see, e.g., Skouibine et al., 2000; Trayanova et al., 2006, 2011; Quarteroni et al., 2017). However, as almost any model, its utility is limited by the inherent resolution of the model. It is useful for understanding cardiac conduction at the tissue level, but it cannot be applied for analyses of conduction close to individual cardiomyocytes. Therefore, detailed models representing individual myocytes have been developed.

Here, we show that the bidomain model can be derived directly from the cell-based EMI model. Classically, the bidomain model is derived using elegant homogenization techniques (see Neu and Krassowska, 1993; Henriquez and Ying, 2021). In the derivation presented here, the deviation between the properties of the bidomain model and the EMI model is seen directly as part of the derivation. In short, the advantage of the present derivation is that it is more straightforward to follow and that it gives indications of where the deviations in the results between the two models stem from.

4.1. Source of Difference Between EMI and Bidomain Solutions

There are essentially three steps in the derivation of the bidomain model where approximations are introduced and thus, most likely, are responsible for the difference in the solutions of the two models. First; the resistance of the intracellular space and the gap junctions are combined into one common and averaged resistance. Second, the average of a function over a volume is approximated by the average of the same function over the surface of the volume. Third, the average of a function on a surface is approximated by the average of the same function on an extended surface. It is beyond the scope of this paper to perform a detailed analysis of these deviations, but based on these observations, it comes as no surprice that the error becomes smaller when the cell size is reduced.

4.2. Differences and Similarities

The EMI model and the bidomain model provide remarkably similar results when the parameters of importance for the conduction velocity are in the normal range. It is safe to claim that the bidomain model represents normal cardiac conduction very well if the scale of interest contains many cells. Certainly, the bidomain model cannot be used to study conduction in the vicinity of individual cells, and it also runs into difficulties for large cells combined with high values of resistance across the gap junctions. It is observed that the bidomain model consistently overestimates the conduction velocity. For normal parameters, the difference is small, but for strongly increased resistance across the gap junctions, the bidomain model significantly overestimates the conduction velocity.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author Contributions

KH and AT: concept and writing. KH: derivation of models and simulations. AT: definition of test problems. Both authors contributed to the article and approved the submitted version.

Funding

This project was supported by the Norwegian Research Council through the EMIx project; 324239.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Agudelo-Toro, A. (2012). Numerical simulations on the biophysical foundations of the neuronal extracellular space (Ph.D. thesis). Niedersächsische Staats-und Universitätsbibliothek Göttingen.

Google Scholar

Amuzescu, B., Airini, R., Epureanu, F. B., Mann, S. A., Knott, T., and Radu, B. M. (2021). Evolution of mathematical models of cardiomyocyte electrophysiology. Math. Biosci. 334:108567. doi: 10.1016/j.mbs.2021.108567

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, R., Andrej, J., Barker, A., Bramwell, J., Camier, J.-S., Dobrev, J. C. V., et al. (2020). MFEM: a modular finite element library. Comput. Math. Appl. 81, 42–74. doi: 10.1016/j.camwa.2020.06.009

CrossRef Full Text | Google Scholar

Bruce, D., Pathmanathan, P., and Whiteley, J. P. (2014). Modelling the effect of gap junctions on tissue-level cardiac electrophysiology. Bull. Math. Biol. 76, 431–454. doi: 10.1007/s11538-013-9927-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Clayton, R., and Panfilov, A. (2008). A guide to modelling cardiac electrical activity in anatomically detailed ventricles. Prog. Biophys. Mol. Biol. 96, 19–43. doi: 10.1016/j.pbiomolbio.2007.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Domínguez, S., Reimer, J., Green, K. R., Zolfaghari, R., and Spiteri, R. J. (2021). “A simulation-based method to study the LQT1 syndrome remotely using the EMI model,” in Emerging Technologies in Biomedical Engineering and Sustainable TeleMedicine (Cham: Springer), 179–189.

Google Scholar

Franzone, P. C., Pavarino, L. F., and Scacchi, S. (2014). Mathematical Cardiac Electrophysiology, Vol. 13. Cham: Springer.

Google Scholar

Heidenreich, E., Ferrero, J., and Rodríguez, J. (2012). “Modeling the human heart under acute ischemia,” in Patient-Specific Computational Modeling. Lecture Notes in Computational Vision and Biomechanics, Vol. 5 (Dordrecht: Springer), 81–103.

Google Scholar

Henriquez, C. S. (2014). A brief history of tissue models for cardiac electrophysiology. IEEE Trans. Biomed. Eng. 61, 1457–1465. doi: 10.1109/TBME.2014.2310515

PubMed Abstract | CrossRef Full Text | Google Scholar

Henriquez, C. S., and Ying, W. (2021). “The bidomain model of cardiac tissue: from microscale to macroscale,” in Cardiac Bioelectric Therapy (Boston, MA: Springer), 211–223.

Google Scholar

Hubbard, M. L., and Henriquez, C. S. (2014). A microstructural model of reentry arising from focal breakthrough at sites of source-load mismatch in a central region of slow conduction. Am. J. Physiol. Heart Circ. Physiol. 306, H1341–H1352. doi: 10.1152/ajpheart.00385.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Jæger, K. H., Edwards, A. G., Giles, W. R., and Tveito, A. (2021a). From millimeters to micrometers; re-introducing myocytes in models of cardiac electrophysiology. Front. Physiol. 12:763584. doi: 10.3389/fphys.2021.763584

PubMed Abstract | CrossRef Full Text | Google Scholar

Jæger, K. H., Edwards, A. G., Giles, W. R., and Tveito, A. (2021b). Mutations change excitability and the probability of re-entry in a computational model of cardiac myocytes in the sleeve of the pulmonary vein. bioRxiv. doi: 10.1101/2021.09.24.461636

CrossRef Full Text | Google Scholar

Jæger, K. H., Edwards, A. G., McCulloch, A., and Tveito, A. (2019). Properties of cardiac conduction in a cell-based computational model. PLoS Comput. Biol. 15:e1007042. doi: 10.1371/journal.pcbi.1007042

PubMed Abstract | CrossRef Full Text | Google Scholar

Jæger, K. H., Hustad, K. G., Cai, X., and Tveito, A. (2021c). Efficient numerical solution of the EMI model representing the extracellular space (E), cell membrane (M) and intracellular space (I) of a collection of cardiac cells. Front. Phys. 8:579461. doi: 10.3389/fphy.2020.579461

CrossRef Full Text | Google Scholar

Jæger, K. H., Hustad, K. G., Cai, X., and Tveito, A. (2021d). “Operator splitting and finite difference schemes for solving the EMI model,” in Modeling Excitable Tissue (Cham: Springer), 44–55.

Google Scholar

Jæger, K. H., and Tveito, A. (2021). “Derivation of a cell-based mathematical model of excitable cells,” in Modeling Excitable Tissue (Cham: Springer), 1–13.

Google Scholar

Jacquemet, V., and Henriquez, C. S. (2009). Genesis of complex fractionated atrial electrograms in zones of slow conduction: a computer model of microfibrosis. Heart Rhythm 6, 803–810. doi: 10.1016/j.hrthm.2009.02.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, J., and Keener, J. P. (2014). Microdomain effects on transverse cardiac propagation. Biophys. J. 106, 925–931. doi: 10.1016/j.bpj.2013.11.1117

PubMed Abstract | CrossRef Full Text | Google Scholar

Linge, S., Sundnes, J., Hanslien, M., Lines, G. T., and Tveito, A. (2009). Numerical solution of the bidomain equations. Philos. Trans. R. Soc. Lond. A 367, 1931–1950. doi: 10.1098/rsta.2008.0306

PubMed Abstract | CrossRef Full Text | Google Scholar

MFEM (2021). MFEM: Modular Finite Element Methods [Software]. Available online at: www.mfem.org.

Neu, J., and Krassowska, W. (1993). Homogenization of syncytial tissues. Crit. Rev. Biomed. Eng. 21, 137–199.

Google Scholar

Niederer, S., Mitchell, L., Smith, N., and Plank, G. (2011a). Simulating human cardiac electrophysiology on clinical time-scales. Front. Physiol. 2:14. doi: 10.3389/fphys.2011.00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer, S. A., Kerfoot, E., Benson, A. P., Bernabeu, M. O., Bernus, O., Bradley, C., et al. (2011b). Verification of cardiac tissue electrophysiology simulators using an n-version benchmark. Philos. Trans. R. Soc. A 369, 4331–4351. doi: 10.1098/rsta.2011.0139

PubMed Abstract | CrossRef Full Text | Google Scholar

Plonsey, R., and Barr, R. C. (2007). Bioelectricity: A Quantitative Approach. Boston, MA: Springer Science & Business Media.

Google Scholar

Qu, Z., Hu, G., Garfinkel, A., and Weiss, J. N. (2014). Nonlinear and stochastic dynamics in the heart. Phys. Rep. 543, 61–162. doi: 10.1016/j.physrep.2014.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Quarteroni, A., Manzoni, A., and Vergara, C. (2017). The cardiovascular system: mathematical modelling, numerical algorithms and clinical applications. Acta Numerica 26, 365–590. doi: 10.1017/S0962492917000046

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, S. F., Stinstra, J. G., and Henriquez, C. S. (2008). Effect of nonuniform interstitial space properties on impulse propagation: a discrete multidomain model. Biophys. J. 95, 3724–3737. doi: 10.1529/biophysj.108.137349

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth, B. J. (1988). The electrical potential produced by a strand of cardiac muscle: a bidomain analysis. Ann. Biomed. Eng. 16, 609–637. doi: 10.1007/BF02368018

PubMed Abstract | CrossRef Full Text | Google Scholar

Rudy, Y. (2012). From genes and molecules to organs and organisms: heart. Comprehensive Biophys. 268–327. doi: 10.1016/B978-0-12-374920-8.00924-3

PubMed Abstract | CrossRef Full Text

Rudy, Y., and Silva, J. R. (2006). Computational biology in the study of cardiac ion channels and cell electrophysiology. Q. Rev. Biophys. 39, 57–116. doi: 10.1017/S0033583506004227

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharifi, M. (2017). Computational approaches to understand the adverse drug effect on potassium, sodium and calcium channels for predicting tdp cardiac arrhythmias. J. Mol. Graphics Modell. 76, 152–160. doi: 10.1016/j.jmgm.2017.06.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaw, R. M., and Rudy, Y. (1997). Ionic mechanisms of propagation in cardiac tissue: roles of the sodium and L-type calcium currents during reduced excitability and decreased gap junction coupling. Circ. Res. 81, 727–741. doi: 10.1161/01.RES.81.5.727

PubMed Abstract | CrossRef Full Text | Google Scholar

Skouibine, K., Trayanova, N., and Moore, P. (2000). A numerically efficient model for simulation of defibrillation in an active bidomain sheet of myocardium. Math. Biosci. 166, 85–100. doi: 10.1016/S0025-5564(00)00019-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Spach, M. S., Heidlage, J. F., Dolber, P. C., and Barr, R. C. (2007). Mechanism of origin of conduction disturbances in aging human atrial bundles: experimental and model study. Heart Rhythm 4, 175–185. doi: 10.1016/j.hrthm.2006.10.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Stinstra, J., MacLeod, R., and Henriquez, C. (2010). Incorporating histology into a 3D microscopic computer model of myocardium to study propagation at a cellular level. Ann. Biomed. Eng. 38, 1399–1414. doi: 10.1007/s10439-009-9883-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Stinstra, J. G., Hopenfeld, B., and Macleod, R. S. (2004). Using models of the passive cardiac conductivity and full heart anisotropic bidomain to study the epicardial potentials in ischemia. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2, 3555–3558. doi: 10.1109/IEMBS.2004.1403999

PubMed Abstract | CrossRef Full Text | Google Scholar

Stinstra, J. G., Shome, S., Hopenfeld, B., and MacLeod, R. S. (2005). Modelling passive cardiac conductivity during ischaemia. Med. Biol. Eng. Comput. 43, 776–782. doi: 10.1007/BF02430957

PubMed Abstract | CrossRef Full Text | Google Scholar

Sundnes, J., Lines, G. T., Cai, X., Nielsen, B. F., Mardal, K.-A., and Tveito, A. (2007). Computing the Electrical Activity in the Heart, Vol. 1. Berlin: Springer Science & Business Media.

Google Scholar

Sundnes, J., Nielsen, B. F., Mardal, K. A., Cai, X., Lines, G. T., and Tveito, A. (2006). On the computational complexity of the bidomain and the monodomain models of electrophysiology. Ann. Biomed. Eng. 34, 1088–1097. doi: 10.1007/s10439-006-9082-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Trayanova, N., Constantino, J., Ashihara, T., and Plank, G. (2011). Modeling defibrillation of the heart: approaches and insights. IEEE Rev. Biomed. Eng. 4, 89–102. doi: 10.1109/RBME.2011.2173761

PubMed Abstract | CrossRef Full Text | Google Scholar

Trayanova, N., Plank, G., and Rodríguez, B. (2006). What have we learned from mathematical models of defibrillation and postshock arrhythmogenesis? application of bidomain simulations. Heart Rhythm 3, 1232–1235. doi: 10.1016/j.hrthm.2006.04.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Tung, L. (1978). A bi-domain model for describing ischemic myocardial dc potentials (Ph.D. thesis). Massachusetts Institute of Technology.

Google Scholar

Tveito, A., Jæger, K. H., Kuchta, M., Mardal, K.-A., and Rognes, M. E. (2017a). A cell-based framework for numerical modeling of electrical conduction in cardiac tissue. Front. Phys. 5:48. doi: 10.3389/fphy.2017.00048

CrossRef Full Text | Google Scholar

Tveito, A., Jæger, K. H., Lines, G. T., Paszkowski, Ł., Sundnes, J., Edwards, A. G., et al. (2017b). An evaluation of the accuracy of classical models for computing the membrane potential and extracellular potential for neurons. Front. Comput. Neurosci. 11:27. doi: 10.3389/fncom.2017.00027

PubMed Abstract | CrossRef Full Text | Google Scholar

Veeraraghavan, R., Gourdie, R. G., and Poelzing, S. (2014). Mechanisms of cardiac conduction: a history of revisions. Am. J. Physiol. Heart Circ. Physiol. 306, H619–H627. doi: 10.1152/ajpheart.00760.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Vigmond, E., Dos Santos, R. W., Prassl, A., Deo, M., and Plank, G. (2008). Solvers for the cardiac bidomain equations. Prog. Biophys. Mol. Biol. 96, 3–18. doi: 10.1016/j.pbiomolbio.2007.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinberg, S. (2017). Ephaptic coupling rescues conduction failure in weakly coupled cardiac tissue with voltage-gated gap junctions. Chaos 27, 093908. doi: 10.1063/1.4999602

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, F., Qu, Z., Yang, J., Baher, A., Weiss, J. N., Garfinkel, A., et al. (2004). A simulation study of the effects of cardiac anatomy in ventricular fibrillation. J. Clin. Invest. 113, 686–693. doi: 10.1172/JCI17341

PubMed Abstract | CrossRef Full Text | Google Scholar

Zemzemi, N., Bernabeu, M. O., Saiz, J., Cooper, J., Pathmanathan, P., Mirams, G. R., et al. (2013). Computational assessment of drug-induced effects on the electrocardiogram: from ion channel to body surface potentials. Br. J. Pharmacol. 168, 718–733. doi: 10.1111/j.1476-5381.2012.02200.x

PubMed Abstract | CrossRef Full Text | Google Scholar