Methods



Fig. 2.1
Equivalent circuit diagram for the cell membrane. The cell membrane is modeled as a capacitor that is charged to the membrane voltage $$V_m$$ by the Nernst potentials ($$V_{{\text {Na}}^+}$$,…) of different ion species. Variable resistors indicate varying individual permeability of the membrane for each ion species. The current source corresponds to ion transporters in the membrane necessary for ion homeostasis




2.1.1 Single Cell Dynamics


Following Hodgkin’s and Huxleys’s scheme, the electrical behavior of the cell membrane can be modeled by an equivalent circuit as shown in Fig. 2.1: The (impermeable part of the) cell membrane is represented by a capacitor with capacity $$C_m$$. It can be charged by the accumulation of charge on its two plates (the inner and outer surfaces of the cell membrane) and the potential difference $$V_m=v_i-v_e$$ between the intracellular and the extracellular potential is called membrane potential, as explained in Sect. 1.​2.​1. This happens through passive and active ion channels that determine the permeability and thus the electrical conductivity of the membrane for specific ionic species ($$g_\text {x}$$ for ion species $$\text {x}$$). Denoting the current through the membrane of ion species $$\text {x}$$ by $$I_\text {x}$$ and the charge on the capacitor by $$Q$$, Kirchhoff’s current law implies that the sum of all currents (through the vertical elements in Fig. 2.1) is zero


$$\begin{aligned} 0&=\frac{\mathrm {d} Q}{\mathrm {d}t}+ I_\text {pump}+I_{\text {Na}^{+}}+I_{\text {K}^{+}}+I_{\text {Ca}^{2+}}+I_{\text {Cl}^{-}}+I_{\dots }\end{aligned}$$

(2.1)



$$\begin{aligned} \Rightarrow \quad \frac{\mathrm {d} Q}{\mathrm {d}t}&=C_m\frac{\mathrm {d} V_m}{\mathrm {d}t} = - I_\text {pump}-I_{\text {Na}^{+}}-I_{\text {K}^{+}}-I_{\text {Ca}^{2+}}-I_{\text {Cl}^{-}}-I_{\dots } \end{aligned}$$

(2.2)
resulting in the above equation for the membrane potential $$V_m$$. In the above equation, the convention was used that positive currents point out of the cell, i.e. from intracellular to extracellular, which implies that the charge $$Q$$ is the one on the intracellular “plate” of the capacitor. The right hand side is sometimes summarized in a term $$I_\text {ion}$$. Each term $$I_\text {x}$$, according to the circuit diagram in Fig. 2.1 is of the form


$$\begin{aligned} I_\text {x} = g_\text {x}(V_m-V_\text {x}), \end{aligned}$$

(2.3)
where $$V_\text {x}$$ is the reversal potential of the ion species $$\text {x}$$. The dots in Eq. (2.2) and Fig. 2.1 represent additional ion species and dynamic elements that could be included in the model. Due to the existence of multiple ion channel types for the same ionic species, the currents associated with each species in Fig. 2.1 have to be interpreted as net currents and $$g_\text {x}$$ as a net (average) conductivity. Thus, the model could be refined by treating each ion channel type as a separate circuit element. This is indeed done in physiologically detailed models such as the Wang-Sobie model [15], but it is not necessary for the purpose of this thesis.

In contrast, an important aspect that has to be taken into account and that is at the basis of the nonlinear phenomenon of action potential formation is the variability of the conductances $$g_\text {x}$$. These are functions of time as well as the quantity gating the channel (membrane potential, individual ion concentrations or mechanical stress, cf. Sect. 1.​2.​1 and Ref. [16]). As mentioned in the introduction of this chapter, the development of the Hodgkin-Huxley formalism preceded the discovery of the molecular functioning of ion channels. Still, from the measured current traces in voltage clamp experiments, it was already clear that the dynamics of the individual ion conductances can themselves only be modeled by dynamical systems for the so-called gating variables. These describe the transitions between different states of the ion channel or its subunits, the main parameters being the transition rates that are, in turn, determined by the quantity responsible for gating the channel (e.g., $$V_m$$ for voltage-gated channels). Depending on the absolute number of channels, a stochastic or deterministic model for the gating variables is required to accurately model the overall membrane conductance for a particular ionic species [17, 18]. As this thesis deals with macroscopic systems only, deterministic dynamics modeled by differential equations will be assumed in the following. In general, each current term in Eq. (2.2) is thus supplemented by a number of additional ordinary differential equations describing its (in-)activation dynamics. The full system of equations describing a single cell electrophysiological model in its most abstract form thus reads


$$\begin{aligned} C_m\frac{\mathrm {d} V_m}{\mathrm {d}t}&= \overbrace{-I_\text {pump}(V_m,\mathbf {h})-\sum _\text {x} I_\text {x}(V_m,\mathbf {h})}^{-I_\text {ion}(V_m,\mathbf {h})}\\ \nonumber \frac{\mathrm {d} \mathbf {h}}{\mathrm {d}t}&= \mathbf {H}(V_m,\mathbf {h}), \end{aligned}$$

(2.4)
where the vector $$\mathbf {h}$$ consists of the gating and all further variables necessary to model the individual components (this may also include ion concentrations, etc.). All current electrophysiological cell models have the form of Eq. (2.4) and the number of terms on the right hand side varies with model complexity, either absorbing many summands into one term or splitting some of them up into individual contributions (see Sect. 2.1.6 for details).


2.1.2 Bi-domain Description of Cardiac Tissue


The single-cell model introduced in Sect. 2.1.1 can be thought of as a description for an average piece of the cell membrane without any spatial extent. In this case, $$C_m$$, $$Q$$ and the current $$I_\text {ion}$$ in Eq. (2.4) have to be interpreted as capacitance, charge and current per unit membrane area. This can be used to construct a model of a $$d$$-dimensional domain of cardiac tissue by extending the so-called core conductor model of neuronal cable theory [19], an idea that emerged during efforts in the 1970s to model surface electrocardiograms [20–22]. The resulting mathematical framework is currently viewed as the standard and most accurate model of cardiac tissue as a continuum [23–25]. If an extended slab of cardiac tissue is considered, the intracellular potentials, $$v_i(\mathbf {x},t)$$ and $$v_e(\mathbf {x},t)$$, respectively, become functions of the position $$\mathbf {x}$$ in space, i.e. $$\mathbf {x}\in \mathbb R^d,t\in \mathbb {R}$$. The idea behind this is that both intracellular and extracellular spaces are defined in the whole tissue domain (thus occupying the same physical space) and have conductivity tensors $${\underline{\sigma }}_i$$ and $${\underline{\sigma }}_e$$, respectively, that are averaged quantities over multiple cells. The relative volumes of the physical intracellular and extracellular spaces and the geometrical arrangement of cells are implicitly contained in these averaged conductivities, as is the additional resistance introduced by gap junctions in the intracellular domain. The two domains represent spatially extended versions of the capacitor plates considered in the single cell model. Current densities within each domain caused by the corresponding potentials are:


$$\begin{aligned} \begin{aligned} \mathbf {j}_i&= -{\underline{\sigma }}_i\nabla v_i\\ \mathbf {j}_e&= -{\underline{\sigma }}_e\nabla v_e\end{aligned} \end{aligned}$$

(2.5)
The current densities enter continuity equations for the charge densities $$q_i$$ and $$q_e$$ in the intracellular and extracellular domain, respectively. It is assumed that no charge accumulates anywhere in the tissue, such that $$\partial q_i/\partial t=-\partial q_e/\partial t=:\partial q/\partial t$$ and


$$\begin{aligned} \begin{aligned} \nabla \cdot \mathbf {j}_i +\frac{\partial q}{\partial t}&= -i_\text {ion}\\ \nabla \cdot \mathbf {j}_e -\frac{\partial q}{\partial t}&= i_\text {ion}, \end{aligned} \end{aligned}$$

(2.6)
where the lowercase quantities $$q$$ and $$i_\text {ion}$$ are charge and membrane current per unit volume, respectively. Otherwise, the same sign conventions as in Sect. 2.1.1 apply. The membrane currents appear as sources and sinks for the extracellular and intracellular domain, respectively, as crossing the cell membrane is the only way for charge to leave one of the two domains (in the absence of current injection via electrodes). As the corresponding quantities $$Q$$ and $$I_\text {ion}$$ are given per unit membrane area, these quantities have to be converted using the surface-to-volume ratio $$\beta $$ (units: area per volume) which specifies how much membrane surface there is on average for a given volume of tissue. Combining Eqs. (2.5) and (2.6), the Bi-domain equations for the tissue domain $${\mathcal {D}}$$ are obtained:


$$\begin{aligned}&\displaystyle \nabla \cdot {\underline{\sigma }}_i\nabla v_i\,=\,-\nabla \cdot \mathbf {j}_i\,=\beta \left( \frac{\partial Q}{\partial t}+I_\text {ion}\right) \,=\,\beta \left( C_m\frac{\partial (v_i-v_e)}{\partial t}+I_\text {ion}\right)&\text {in}\,{\mathcal {D}}\end{aligned}$$

(2.7a)



$$\begin{aligned}&\displaystyle \nabla \cdot {\underline{\sigma }}_e\nabla v_e\,=\,-\nabla \cdot \mathbf {j}_e\,=\,-\beta \left( \frac{\partial Q}{\partial t}+I_\text {ion}\right) \,=\,-\beta \left( C_m\frac{\partial (v_i-v_e)}{\partial t}+I_\text {ion}\right)&\text {in}\,{\mathcal {D}} \end{aligned}$$

(2.7b)
These equations have to be completed by boundary conditions for $$v_i$$ and $$v_e$$ at the boundary $${\partial {\mathcal {D}}}$$ of the tissue domain $${\mathcal {D}}$$. Assuming the space outside the tissue domain $${\mathcal {D}}$$ is a mono-domain with conductivity $${\underline{\sigma }}_o$$ and a corresponding outside potential $$v_o$$, these boundary conditions are


$$\begin{aligned}&\,\,\mathbf {n}\cdot {\underline{\sigma }}_i\nabla v_i= 0 \text {on}\, {\partial {\mathcal {D}}}\end{aligned}$$

(2.8a)



$$\begin{aligned}&\,\qquad \quad \,\,\,\,v_e= v_o \text {on}\, {\partial {\mathcal {D}}}\end{aligned}$$

(2.8b)



$$\begin{aligned}&\,\mathbf {n}\cdot {\underline{\sigma }}_e\nabla v_e= \mathbf {n}\cdot {\underline{\sigma }}_o\nabla v_o\quad \text {on}\, {\partial {\mathcal {D}}} \end{aligned}$$

(2.8c)
where $$\mathbf {n}$$ is a local unit vector perpendicular to the boundary $${\partial {\mathcal {D}}}$$. Equation (2.8a) is a no-flux boundary condition, expressing the fact that the intracellular space ends at the tissue boundary and thus no current in the intracellular domain can cross this boundary. Equations (2.8b) and (2.8c) ensure the continuity of the extracellular potential and currents at the interface between the extracellular space inside the tissue and the mono-domain outside the tissue (e.g., the bath).

For practical calculations, a different linear combination of the system (2.7) is used which corresponds to a change of variables from $$v_i$$ and $$v_e$$ to $$V_m=v_i-v_e$$ and $$v_e$$:


A328193_1_En_2_Equ130_HTML.gif
Reordering the terms in both equations yields


$$\begin{aligned}&\qquad \qquad \quad C_m\frac{\partial V_m}{\partial t}=\displaystyle \frac{1}{\beta }\nabla \cdot {\underline{\sigma }}_i\nabla (V_m+v_e)-I_\text {ion}(V_m,\mathbf {h})&\text {in}\,\, {\mathcal {D}}\end{aligned}$$

(2.9a)



$$\begin{aligned}&\qquad \qquad \qquad \quad \frac{\partial \mathbf {h}}{\partial t}\, = \mathbf {H}(V_m,\mathbf {h})&\text {in}\,\, {\mathcal {D}}\end{aligned}$$

(2.9b)



$$\begin{aligned}&\displaystyle \nabla \cdot ({\underline{\sigma }}_e+{\underline{\sigma }}_i)\nabla v_e= -\nabla \cdot {\underline{\sigma }}_i\nabla V_m&\text {in}\,\, {\mathcal {D}}. \end{aligned}$$

(2.9c)
Having been omitted so far, Eq. (2.9b) was reinserted from Eq. (2.4) for completeness. Equation (2.9c) does not contain $$t$$ and thus states that the extracellular potential can be determined from the membrane potential at any given instant. In Eq. (2.9a), the term $$\nabla \cdot {\underline{\sigma }}_i\nabla (V_m+v_e)$$ is more favorable than the seemingly simpler $$-\nabla \cdot {\underline{\sigma }}_e\nabla v_e$$, because in this way, the equation stays valid in the limit $${\underline{\sigma }}_e=\alpha {\underline{\sigma }}_e^*$$ for $$\alpha \rightarrow \infty $$. This will be useful in the next section. For the boundary conditions (2.8), the change of variables can be achieved by replacing $$V_m$$ in Eq. (2.8a):


$$\begin{aligned}&\mathbf {n}\cdot {\underline{\sigma }}_i\nabla (V_m+v_e) = 0&\text {on}\,\, {\partial {\mathcal {D}}}\quad \quad \quad \quad \end{aligned}$$

(2.10a)



$$\begin{aligned}&\qquad \qquad \qquad \qquad v_e= v_o&\text {on}\,\, {\partial {\mathcal {D}}}\end{aligned}$$

(2.10b)



$$\begin{aligned}&\qquad \qquad \mathbf {n}\cdot {\underline{\sigma }}_e\nabla v_e= \mathbf {n}\cdot {\underline{\sigma }}_o\nabla v_o&\text {on}\,\, {\partial {\mathcal {D}}} \end{aligned}$$

(2.10c)


2.1.3 Mono-domain Descriptions of Cardiac Tissue


For full numerical simulations using the bi-domain equations (2.9), one would have to time-step the membrane potential $$V_m$$ using Eq. (2.9a) from a given present state and calculate the corresponding extracellular potential by solving Eq. (2.9c) for $$v_e$$. The latter requires solving a large system of linear equations, an operation that is computationally very costly. Because of this, there exists an approximation that only involves a single partial differential equation (PDE) for $$V_m$$ that can be solved by time stepping, the so-called mono-domain equation. Although the mono-domain model is used extensively in the literature and is regarded a valid model in its own right without further justification, it is important to know its connection to the full and more realistic bi-domain formulation. In particular, this is true for understanding the effects of electric fields near tissue boundaries, which is one of the aims of this thesis. In the following two sub-sections, two different ways of deriving the mono-domain equation with different physical interpretation will be presented. When working with a mono-domain equation, one has to ensure that modifications to the equations are consistent with the choice of interpretation.


2.1.3.1 Neglecting the Extracellular Potential


Assume a very well-conducting extracellular space,2 i.e. $$\alpha \gg 1$$ for $${\underline{\sigma }}_e=\alpha {\underline{\sigma }}_e^*$$. In the limit $$\alpha \rightarrow \infty $$, Eq. (2.9c) reduces to $$\nabla \cdot {\underline{\sigma }}_e^*\nabla v_e=0$$. This means that for large $$\alpha $$ the source term for the extracellular potential becomes negligible and one possible form of solution is $$v_e= -\mathbf {E}\cdot \mathbf {x}+c$$ with some fixed electric field vector $$\mathbf {E}$$ and an arbitrary constant $$c$$ (if boundary conditions  (2.10b) and (2.10c) permit).3 For this solution, it has to be assumed that $${\underline{\sigma }}_e^*$$ is spatially constant. If we assume the extracellular space is insulated ($$\mathbf {n}\cdot \sigma _e^*\nabla v_e=0$$) or grounded ($$v_e=0$$) on $${\partial {\mathcal {D}}}$$, this solution reduces even further to $$v_e= c$$ or $$v_e\equiv 0$$, respectively. In contrast, Eq. (2.9a) is still valid. The mono-domain equation in this scenario therefore is


$$\begin{aligned}&\frac{\partial V_m}{\partial t}=\nabla \cdot {{\underline{D}}}\nabla (V_m+v_e)-I_\text {ion}/C_m\quad \text {in} \,\,{\mathcal {D}}, \end{aligned}$$

(2.11)
where the diffusion tensor is given by $${{\underline{D}}}={\underline{\sigma }}_i/(\beta C_m)$$ and $$v_e$$ is merely a parameter for the problem, which is usually set to zero following the line of thoughts above. To make this problem well-defined, the boundary condition (2.10a) is sufficient:


$$\begin{aligned} \mathbf {n}\cdot {{\underline{D}}}\nabla (V_m+v_e) = 0\quad \text {on}\,\,{\partial {\mathcal {D}}} \end{aligned}$$

(2.12)


2.1.3.2 Equivalent Mono-domain Description


In Sect. 2.1.3.1, a mono-domain formulation of the bi-domain model was obtained by changing the physics of the problem. The extent to which the dynamics is altered can be significantly reduced in the case were the intracellular and extracellular conductivity tensors are scaled versions of the same tensor: $${\underline{\sigma }}_e=\alpha {\underline{\sigma }}_i$$. In this case, Eq. (2.7b) can be subtracted from Eq. (2.7a):


$$\begin{aligned} \nabla \cdot {\underline{\sigma }}_i\nabla V_m= \left( 1+\alpha ^{-1}\right) \beta \left( C_m\frac{\partial V_m}{\partial t}+I_\text {ion}\right) \end{aligned}$$
Rearranging for $${\partial V_m}/{\partial t}$$ yields the mono-domain equation


$$\begin{aligned}&\displaystyle \frac{\partial V_m}{\partial t}=\nabla \cdot {{\underline{D}}}\nabla V_m-I_\text {ion}(V_m,\mathbf {h})/C_m&\quad \text {in}\,\,{\mathcal {D}}\quad \quad \quad \end{aligned}$$

(2.13a)



$$\begin{aligned}&\displaystyle \frac{\partial \mathbf {h}}{\partial t}=\mathbf {H}(V_m,\mathbf {h})&\quad \text {in}\,\,{\mathcal {D}} \end{aligned}$$

(2.13b)
where now $${\underline{D}}={\underline{\sigma }}_i/[\beta C_m(1+\alpha ^{-1})]$$. In the limit $$\alpha \rightarrow \infty $$, the definition of $${\underline{D}}$$ is equal to that in Eq. (2.11). Similarly, the $$v_e$$ term in Eq. (2.11) becomes irrelevant if $${\underline{\sigma }}_e^*$$ is a scaled $${\underline{\sigma }}_i$$, because then $$\nabla \cdot {\underline{D}}\nabla v_e=0$$ follows from $$\nabla \cdot {\underline{\sigma }}_e^*\nabla v_e=0$$. Equations (2.11) and (2.13) are thus consistent in the sense that they are equivalent when the assumptions for both mono-domain descriptions are fulfilled.

Note that the mono-domain description (2.13) is completely equivalent to the full bi-domain model, but numerically cheap compared to Eq. (2.9). Consequently, for the case $${\underline{\sigma }}_e=\alpha {\underline{\sigma }}_i$$, numerical simulation of the bi-domain model is usually superfluous. There is, however, one caveat: from the set of boundary conditions (2.10), no boundary condition for the membrane potential alone can be derived. Also, unlike in Eq. (2.11), $$v_e$$ is not decoupled from the membrane potential dynamics and so Eq. (2.10a) cannot be enforced unless $$v_e$$ is calculated, which would defeat the original purpose of a mono-domain formulation. The usual (somewhat nonphysical) solution to this dilemma is to apply the ad-hoc no-flux boundary condition


$$\begin{aligned} \mathbf {n}\cdot {{\underline{D}}}\nabla V_m= 0\quad \text {on}\,\, {\partial {\mathcal {D}}}. \end{aligned}$$

(2.14)
This is sufficient for modeling the dynamics inside a tissue domain $${\mathcal {D}}$$, where the boundary of the domain is not meant to be a physical part of the system but rather a computational necessity. However, it is only a valid simplification, if one can rule out any influence of the boundary condition on the effects one wants to study. Sometimes, the no-flux boundary condition is physically justified, for example, if the extracellular space is known to be insulated electrically at the tissue boundary. However, in particular in this thesis, investigating effects at the boundaries of the tissue due to the conductive connection of the extracellular space to the surrounding environment is the original purpose of the modeling. In this case, an alternative solution is to assume that $$\alpha $$ and the outside conductivity $${\underline{\sigma }}_o$$ are large, which decouples the corresponding potentials from the dynamics of the membrane potential and makes them parameters as seen for Eq. (2.11). This way, an extracellular or outside potential (e.g., a uniform electric field) can be assumed rather than calculated, which makes it possible to use Eq. (2.10a) itself or a boundary condition obtained by subtracting Eq. (2.10c) from Eq. (2.10a):


$$\begin{aligned} \mathbf {n}\cdot {\underline{D}}\nabla (V_m+v_e)&= 0\qquad \qquad \qquad \,\text {on}\,\, {\partial {\mathcal {D}}}\\ \nonumber \text {or}\end{aligned}$$

(2.15)



$$\begin{aligned} \mathbf {n}\cdot \underbrace{\alpha {\underline{\sigma }}_i}_{={\underline{\sigma }}_e}\nabla V_m&= -\mathbf {n}\cdot {\underline{\sigma }}_o\nabla v_o\quad \text {on}\,\,{\partial {\mathcal {D}}} \end{aligned}$$

(2.16)
The second boundary condition (2.16) will be used in Sect. 3.​2 to study the effect of electric fields near boundaries of the tissue. As for the bi-domain model (2.9), the mono-domain equation (2.13) has to be supplemented by a concrete model $$\mathbf {H}$$ for the local variables $$\mathbf {h}$$ and the corresponding ionic currents $$I_\text {ion}(V_m,\mathbf {h})$$. The models used in this thesis will be introduced in Sect. 2.1.6.


2.1.4 Anisotropy


Although all simulations in this thesis are carried out assuming isotropic diffusion, for completeness of the mathematical framework, the construction of the corresponding tensors for the specific requirements of cardiac tissue shall be mentioned here briefly: The anisotropic conductivity of cardiac tissue due to the fiber orientation and sheet structure of the muscle (see Sect. 1.​3) can be implemented in the bi-domain and in the mono-domain model by choosing the diffusion tensor $${\underline{D}}$$ or equivalently the conductivity tensors. $${\underline{D}}$$ can be constructed from its three eigenvectors, each of which represents a direction in which a potential gradient leads to a current density in the same direction, with a conductivity equal to the corresponding eigenvalue. The fiber orientation determines one of these directions, the other two are uniquely defined by the normal vector of the laminar sheets (the remaining vector thus corresponds to the direction perpendicular to the fibers within the sheet). These three eigenvectors form an orthonormal basis of $$\mathbb R^3$$, implying that $${\underline{D}}$$ must be a real, symmetric matrix. Using these properties, the entire tensor $${\underline{D}}$$ is determined by two unit vectors in space defining the fiber direction $$\mathbf {f}=(f_x,f_y,f_z)$$ and the direction normal to the sheet plane $$\mathbf {n}=(n_x,n_y,n_z)$$ and by diffusion constants $$D_\parallel $$ along the fibers, $$D_\perp $$ perpendicular to the fibers (within the sheet) and $$D_{\perp \perp }$$ normal to the sheets (see Eq. 17 in [25]):


$$\begin{aligned} D_{ij} = \delta _{ij}D_\perp + (D_\parallel -D_\perp )f_{i}\,f_{j} + (D_{\perp \perp }-D_\perp )n_in_j,\qquad i,j\in \{x,y,z\}\qquad \qquad \end{aligned}$$

(2.17)
A common simplification is to ignore the sheet structure and thus assume $${D_{\perp \perp }=D_\perp }$$, which eliminates one of the terms in Eq. (2.17). Note also that, in general, both the directions and diffusion constants may vary spatially.

The so-called anisotropy ratio is defined as the quotient $$D_\parallel /D_\perp $$. In the bi-domain equations, the two conductivity tensors $${\underline{\sigma }}_e$$ and $${\underline{\sigma }}_i$$ can be chosen independently. While the fiber orientation for both is the same, because it is given by the geometrical arrangement of the cells, the anisotropy ratio can be different. The case of equal anisotropy ratio is the one considered in the mono-domain approximation of Sect. 2.1.3.2.


2.1.5 The Phase-Field Method


As outlined in Sect. 2.1.3, boundary conditions constitute an important part of any problem described by a PDE on a spatial domain $${\mathcal {D}}$$. It was mentioned that a widely-used boundary condition for mono-domain cardiac models is no-flux. Another possibility which was derived for the case of a predetermined extracellular/outside potential is that of Eq. (2.16) in the same section. Both boundary conditions are of the Neumann type, i.e. they fix the normal fluxes $$\mathbf {n}\cdot {\underline{D}}\nabla u$$ on the boundary $${\partial {\mathcal {D}}}$$. As will be shown in Sect. 2.2.3, for a finite-difference numerical scheme, no-flux boundary conditions are straightforward to implement, if the boundaries of the domain $${\mathcal {D}}$$ are parallel to the coordinate axes, e.g., rectangles in cartesian coordinates or (fractional) annuli in polar coordinates. For arbitrarily shaped domains $${\mathcal {D}}$$, there is no such natural way of imposing zero flux on the domain boundary. In this case, one usually resorts to more complex and computationally intensive numerical schemes like finite-volume or finite-element methods (see Sect. 6.3 of Ref. [25] and references therein).

In Ref. [26], Fenton et al. introduce a phase-field method to impose no-flux boundary conditions on arbitrary geometries using finite differences on a simple (e.g., rectangular) computational domain $${\mathcal {D}}_{c}$$ which contains the actual physical domain $${\mathcal {D}}$$. In the following, this method will be explained in detail and extended to arbitrary Neumann boundary conditions of the type of Eq. (2.16). Furthermore, mathematical improvements to the convergence proof will enable a better assessment of the method’s numerical accuracy. The generalized method will be used extensively in Sect. 3.​2 to study the effect of electric fields at boundaries of curved tissue domains.


2.1.5.1 Description of the Method


The phase-field method is based on defining a phase field $$\phi $$, which is a smoothed version of the characteristic function of $${\mathcal {D}}$$, i.e. $$\phi \approx 1$$ inside the domain $${\mathcal {D}}\subset {\mathcal {D}}_{c}$$ and $$\phi \approx 0$$ otherwise, with a smooth transition in between. In this thesis, a different dynamical system than suggested in Ref. [26] shall be used for obtaining $$\phi $$ from the sharp, original characteristic function $$\phi _0$$:


$$\begin{aligned} \frac{\partial \phi }{\partial t}=(\phi _0-\phi )+\xi ^2\Delta \phi . \end{aligned}$$

(2.18)
In this equation, $$\phi $$ approaches a smooth phase field for $$t\rightarrow \infty $$, which stays close to $$\phi _0$$ and has a diffusive boundary with a characteristic length scale $$\xi $$ of the resulting transition at the interface between $$\phi =1$$ and $$\phi =0$$. It does not suffer from the stability issues in higher dimensions mentioned in [26] for the double-well potential method defined therein. Hence, any method that finds the equilibrium of Eq. (2.18) may be used instead of explicit time-stepping. In the following, this steady-state solution will be referred to as $$\phi $$. An example for $$\phi _0$$ and the resulting $$\phi $$ is depicted in Fig. 2.2a, b respectively. One property that will prove to be important for the line of thoughts below is that the normal slope of $$\ln \phi $$ decreases rapidly and monotonically within the interface region (in the outward direction) until reaching an approximately constant value on the order of $$-1/\xi $$, as shown in Fig. 2.2c. Hence, at large normal distances from the interface, $$\phi $$ decays approximately exponentially and $$\partial ^2\ln \phi /\partial n^2\approx 0$$. As illustrated in Fig. 2.2d, the normal direction is defined by the direction of the gradient $$\nabla \ln \phi $$ itself.

A328193_1_En_2_Fig2_HTML.gif


Fig. 2.2
Phase field behavior. a Initial binary phase field $$\phi _0$$, i.e. the characteristic function of $${\mathcal {D}}$$. b Steady-state solution of Eq. (2.18) with $$\phi _0$$ from (a) and $$\xi =2$$. c $$|\nabla \ln \phi |$$ which is equal to $$-\partial \ln \phi /\partial n$$ if one defines the normal vector $$\mathbf {n}$$ to be anti-parallel to $$\nabla \ln \phi $$. The logarithmic normal derivative increases quickly from 0 inside $${\mathcal {D}}$$ to an approximately constant value $$\approx $$ $$1/\xi $$ outside $${\mathcal {D}}$$. d $$\ln \phi $$ shown as a color plot and the corresponding gradient shown as arrows for the region indicated by a white rectangle in (c)

The phase field method then consists of substituting the original diffusion term (e.g., in Eq. (2.13)) by a modified term that involves the phase field $$\phi $$:


$$\begin{aligned} \nabla \cdot ({\underline{D}}\nabla u) \quad \rightarrow \quad \frac{1}{\phi }\nabla \cdot (\phi {\underline{D}}\nabla u) \end{aligned}$$

(2.19)
The resulting PDE after the substitution is assumed to be defined on the whole domain $${\mathcal {D}}_{c}$$, but only the part of the solution in $${\mathcal {D}}\subset \mathcal D_c$$ is taken to have physical meaning. Fenton et al. present a mathematical argument in Ref. [26] for a one-dimensional tissue domain which suggests that the substitution (2.19) for $$\mathcal D_c$$ leads to a solution in $$\mathcal D$$ that obeys a no-flux boundary condition


$$\begin{aligned} \mathbf {n}\cdot ({\underline{D}}\nabla u)=0\qquad \text {on }{\partial {\mathcal {D}}}. \end{aligned}$$

(2.20)
A similar proof for higher dimensions based on the same idea was provided in Ref. [27]. In the following, this idea and its shortcomings in terms of the necessary mathematical assumptions will be presented and an argument will be given how these deficits may be resolved.


2.1.5.2 Original Convergence Proof


Omitting any additional non-spatial degrees of freedom for the reaction term, the general PDE for a reaction-diffusion system with the replacement of Eq. (2.19) is


$$\begin{aligned} \frac{\partial u}{\partial t}&=f(u)+\frac{1}{\phi }\nabla \cdot (\phi {\underline{D}}\nabla u), \end{aligned}$$

(2.21)
where $$f$$ is assumed to be bounded for bounded $$u$$. The multi-dimensional equivalent of the mathematical argument given in [26] starts by integrating Eq. (2.21) over a small volume $$V$$ across the boundary $$\partial D$$, which yields


$$\begin{aligned} \int \limits _V \phi \frac{\partial u}{\partial t}\,\mathrm {d}V&= \int \limits _V \phi f(u)\,\mathrm {d}V + \int \limits _V\nabla \cdot (\phi {\underline{D}}\nabla u)\,\mathrm {d}V = \int \limits _V \phi f(u)\,\mathrm {d}V + \oint \limits _{\partial V}\phi {\underline{D}}\nabla u\cdot \,\mathrm {d}\mathbf {A}. \end{aligned}$$

(2.22)


A328193_1_En_2_Fig3_HTML.gif


Fig. 2.3
Volume integration at the phase-field interface. Sketch of the integration volume for Eq. (2.22) parallel to the boundary $$\partial \mathcal D$$ which is defined by the border between $$\phi _0=1$$ (green) and $$\phi _0=0$$ (gray). Arrows indicate normal vectors on the surfaces of the integration volume corresponding to the individual terms in Eq. (2.23)

The geometrical situation is illustrated in Fig. 2.3, with the volume $$V$$ shaded gray. A similar strategy for a different integration path is given in Ref. [28], but the following argumentation applies universally. Denoting the individual contributions to the surface integral separately leads to


$$\begin{aligned} \int \limits _V \phi \frac{\partial u}{\partial t}\,\mathrm {d}V&= \int \limits _V \phi f(u)\,\mathrm {d}V\nonumber \nonumber \\&\quad + \int \limits _{\partial V_n^1}\phi {\underline{D}}\nabla u\cdot \,\mathrm {d}\mathbf {n}_1 \nonumber \\&\quad + \int \limits _{\partial V_n^2}\phi {\underline{D}}\nabla u\cdot \,\mathrm {d}\mathbf {n}_2\nonumber \\&\quad + \int \limits _{\partial V_t^1}\phi {\underline{D}}\nabla u\cdot \,\mathrm {d}\mathbf {t}_1\nonumber \\&\quad + \int \limits _{\partial V_t^2}\phi {\underline{D}}\nabla u\cdot \,\mathrm {d}\mathbf {t}_2\nonumber \\&\quad +\dots \quad . \end{aligned}$$

(2.23)
The dots indicate that, in the case of three dimensions, there are two more terms for another direction parallel to the boundary $$\partial \mathcal D$$. However, for everything presented in the following, these terms can be treated analogously to the two explicitly mentioned in Eq. (2.23). For a given interface width $$\xi $$ of $$\phi $$, the normal width of the volume $$V$$ is chosen such that the two boundaries of $$V$$ parallel to the boundary $$\partial \mathcal D$$ are at iso-$$\phi $$ levels of $$\phi =1-\varepsilon $$ and $$\phi =\varepsilon $$ for an arbitrarily small $$\varepsilon >0$$” src=”/wp-content/uploads/2016/11/A328193_1_En_2_Chapter_IEq177.gif”></SPAN>.</DIV><br /> <DIV id=Par24 class=Para>By assuming smooth and bounded solutions <SPAN id=IEq178 class=InlineEquation><IMG alt=$$u$$ src= and choosing sufficiently small $$V$$ (corresponding to small enough $$\xi $$ for fixed $$\varepsilon $$), the individual integrands of Eq. (2.23) can be approximated by constants, except for the variation due to the factor $$\phi $$. By then letting $$\xi $$ tend to zero, i.e. reducing the width of $$V$$ normal to the boundary $$\partial \mathcal D$$ and making $$\phi $$ steeper and steeper, the integral on the left hand side and the integrals of the tangential components $$\partial V_t^1$$ and $$\partial V_t^2$$ on the right hand side of Eq. (2.23) vanish, because they are integrated over a volume and a surface, respectively, whose size tends to zero. The same is true for the reaction term, which is assumed to be bounded for bounded $$u$$. Using the fact that $$\phi \approx 1$$ on $$\partial V_n^2$$ and $$\phi \approx 0$$ on $$\partial V_n^1$$, only one of the two remaining integrals actually has a non-zero contribution, resulting in


$$\begin{aligned} 0=\int \limits _{\partial V_n^2}\phi {\underline{D}}\nabla u\cdot \,\mathrm {d}\mathbf {n}_2\approx |\partial V_n^2|{\underline{D}}\nabla u\cdot \mathbf {n}\end{aligned}$$

(2.24)



$$\begin{aligned} \Rightarrow \mathbf {n}\cdot {\underline{D}}\nabla u = 0, \end{aligned}$$

(2.25)
where $$|\partial V_n^2|$$ denotes the surface area of $$\partial V_n^2$$.

If the above reasoning was justified, Eq. (2.25) would indicate that no-flux boundary conditions are actually implemented by the phase field method. However, it is not known a priori, whether the integrands are bounded independent of the choice of $$\xi $$. As the integrands contain derivatives of $$u$$, it is not even sufficient to assume that $$u$$ is bounded. Indeed, it is quite natural to expect that by letting $$\xi $$ tend to zero and thus approaching the limit of a non-differentiable $$\phi $$, the gradients of the solution $$u$$ could blow up. Thus, if the crucial assumption of boundedness is not fulfilled, it can neither be concluded that integrals tend to zero because the integrated volume does, nor that for $$\partial V_n^1$$ the approximation $$\phi {\underline{D}}\nabla u=\varepsilon {\underline{ D}}\nabla u\approx 0$$ holds (for sufficiently small $$\varepsilon $$). In fact, it becomes obvious that the line of thoughts above is not sufficient, if one notes that it would also work without the factor $$\phi $$ in the volume integrals of Eq. (2.23), which corresponds to the additional prefactor $$1/\phi $$ in Eq. (2.19). As will be shown below, this prefactor is vital for the method to yield the desired result.

In a more recent publication [29], Li et al. showed that phase-field methods can be used to model a number of different boundary conditions, no-flux being a special case which can indeed implemented by Eq. (2.21). However, the mechanisms described below by which the phase-field method achieves its goal also explain some important details of the numerical implementation not provided by the very abstract proof in Ref. [29].


2.1.5.3 Mechanism of Boundary Condition Enforcement


Rewriting Eq. (2.21), the terms of the original PDE and the influence of the phase field can be separated:


$$\begin{aligned} \frac{\partial u}{\partial t}&=f(u)+\frac{1}{\phi }\nabla \cdot (\phi {\underline{D}}\nabla u)\nonumber \\&=f(u)+\nabla \cdot ({\underline{D}}\nabla u)+\frac{1}{\phi }(\nabla \phi )\cdot ({\underline{D}}\nabla u) \end{aligned}$$

(2.26)



$$\begin{aligned}&=f(u)+\nabla \cdot ({\underline{D}}\nabla u)+(\nabla (\ln \phi ))\cdot ({\underline{D}}\nabla u), \end{aligned}$$

(2.27)
where “$$\cdot $$” in the last two equations is the dot product between two vectors. In Eq. (2.27), the third term only becomes active, if there are components of the current density $$-{\underline{D}}\nabla u$$ perpendicular to the boundary (i.e. violating the intended no-flux boundary condition), as $$\nabla (\ln \phi )$$, by definition of the phase field, is just a scaled normal vector:


$$\begin{aligned} \nabla (\ln \phi )=\frac{\partial \ln \phi }{\partial n}\mathbf {n} \end{aligned}$$

(2.28)
With this information, Eq. (2.27) reads:


$$\begin{aligned} \frac{\partial u}{\partial t}&=f(u)+\nabla \cdot ({\underline{D}}\nabla u)+\frac{\partial \ln \phi }{\partial n}\mathbf {n}\cdot ({\underline{D}}\nabla u) \end{aligned}$$

(2.29)
Note that the third term resembles an advection term in fluid dynamics, though the velocity field is not divergence free. Shifting $${\underline{D}}$$ to the other side of the scalar product, the velocity field for the quantity $$u$$ is given by


$$\begin{aligned} -\frac{\partial \ln \phi }{\partial n}{\underline{D}}\mathbf {n}. \end{aligned}$$

(2.30)
As $${\underline{D}}$$ is positive definite, $$\mathbf {n}\cdot {\underline{D}}\mathbf {n}\ge 0$$. Taking $$\mathbf {n}$$ to point out of the domain $$\mathcal D$$, according to Fig. 2.2c $$\partial \ln \phi /\partial n<0$$ holds, which means that the velocity (2.30) points into the same half-space as $$\mathbf {n}$$, although it might not be normal to $$\partial \mathcal D$$. Therefore, $$u$$ is “advected” outwards with a velocity proportional to $$|\,\mathrm {d}\ln \phi /\,\mathrm {d}n|$$.

To show that the assumed boundedness of the integrands in Eq. (2.23) is actually present, the so-called maximum principle will be used: If the local reaction term in Eq. (2.21) is ignored for the moment, a pure diffusion equation is obtained. The maximum principle states that the solution $$u(\mathbf {x},t)$$ for $$0\le t<T$$ and $$\mathbf {x}\in \mathcal D_c$$ is attained either on the spatial boundary $$\partial \mathcal D_c$$ or at some point $$\mathbf {x}\in \mathcal D_c$$ at $$t=0$$ [30]. Another way to say this is: The solution will never exceed the maximum of the initial field $$u$$ at $$t=0$$, provided that the values on the boundary $$\partial \mathcal D_c$$ stay below the maximum as well. The same applies to the minimum. For the purpose of the following reasoning, it is assumed that the boundary $$\partial \mathcal D_c$$ of the auxiliary domain is subject to a passive boundary condition which never attains a maximum or minimum, for example a Dirichlet condition $$u=\langle u(t=0)\rangle _{\mathcal D_c}$$ on $$\partial \mathcal D_c$$. Below, it will be argued why the boundary condition on the auxiliary domain is not important for the solution $$u$$ in the physical domain $$\mathcal D$$.

Having the maximum principle at hand, an upper bound for the left hand side of Eq. (2.23) can be obtained by differentiating Eq. (2.21) (still without the local reaction term) with respect to $$t$$ and noting that the time derivative $$\frac{\partial u}{\partial t}$$ itself therefore fulfills the diffusion equation


$$\begin{aligned} \frac{\partial \frac{\partial u}{\partial t}}{\partial t}&=\frac{1}{\phi }\nabla \cdot \left( \phi {\underline{D}}\nabla \frac{\partial u}{\partial t}\right) . \end{aligned}$$

(2.31)
The initial condition $$\left. \frac{\partial u}{\partial t}\right| _{t=0}$$ is connected to the initial condition $$u_0(\mathbf {x})=u(\mathbf {x},t=0)$$ of the original diffusion equation via Eq. (2.21) or (2.29). Given that $$u_0$$ fulfills the no-flux boundary condition in the interface region, the $$\phi $$-dependent terms in Eq. (2.29) drop out and thus $$\frac{\partial u}{\partial t}$$ is bounded independent of $$\phi $$ by the maximum principle for all $$t$$:


$$\begin{aligned} \left| \left. \frac{\partial u}{\partial t}\right| _{(\mathbf {x}, t)}\right| \le C \qquad \forall \mathbf {x}, t\quad \text {and}\quad \forall \xi . \end{aligned}$$

(2.32)
In the context of excitable media, with the source term $$f=I_\text {ion}$$ introducing additional $$\frac{\partial u}{\partial t}$$ contributions, usually, the fastest time scale of $$u$$ is given by the rise velocity upon activation. The maximum principle therefore ensures that $$\partial u/\partial t$$ is not increased further by the $$\phi $$-dependent diffusion dynamics and $$C$$ can be interpreted as the upstroke velocity of the action potential.

The boundedness of the integrand involving the reaction term $$f(u)$$ in Eq. (2.23) is seen immediately, because $$u$$ itself is subject to a diffusion equation and thus obeying the maximum principle, if $$f(u)$$ is omitted in Eq. (2.21). Again, this means that the $$\phi $$-dependent diffusion dynamics cannot increase or decrease $$u$$ further than determined by the reaction dynamics. Therefore $$u$$ is bounded and so is $$f(u)$$.

For the remaining terms on the right-hand side of Eq. (2.23), a simpler special case is considered, namely that of a diffusion tensor $${\underline{D}}$$ which is nearly spatially constant in the interface region and of an approximately straight boundary, such that spatial derivatives of $${\underline{D}}$$ and $$\mathbf {n}$$ can be neglected.4 Under these assumptions, differentiating Eq. (2.29) in the normal direction yields:


$$\begin{aligned} \frac{\partial \frac{\partial u}{\partial n}}{\partial t}&=\nabla \cdot \left( {\underline{D}}\nabla \frac{\partial u}{\partial n}\right) +\frac{\partial \ln \phi }{\partial n}\mathbf {n}\cdot \left( {\underline{D}}\nabla \frac{\partial u}{\partial n}\right) +\frac{\partial ^{2} \ln \phi }{\partial n^{2}}\mathbf {n}\cdot ({\underline{D}}\nabla u)\end{aligned}$$

(2.33)



$$\begin{aligned}&=\frac{1}{\phi }\nabla \cdot \left( \phi {\underline{D}}\nabla \frac{\partial u}{\partial n}\right) +\frac{\partial ^{2} \ln \phi }{\partial n^{2}}\mathbf {n}\cdot ({\underline{D}}\nabla u) \end{aligned}$$

(2.34)
If only the first term on the right hand side was present, $$\frac{\partial u}{\partial n}$$ would be bounded by the maximum principle. Due to the shape of the phase field (compare Fig. 2.2 and the beginning of this section), $$\frac{\partial ^{2} \ln \phi }{\partial n^{2}}$$ is negative everywhere in the interface region and approaching zero from below in the outward direction. Therefore, the second term on the right hand side is only present in the interface region and opposite in sign to the normal component of the flux $${\underline{D}} \nabla u$$. As $${\underline{D}}$$ is positive definite, the induced change of $$\partial u/\partial n=\mathbf {n}\cdot \nabla u$$ leads to a reduction of the normal component of $${\underline{D}}\nabla u$$, because in general, for a vector $$\mathbf {v}$$ with a normal component $$v_n$$


$$\begin{aligned} \frac{\partial \mathbf {n}\cdot {\underline{D}}\mathbf {v}}{\partial v_n}=\mathbf {n}\cdot {\underline{D}}\mathbf {n} \ge 0, \end{aligned}$$

(2.35)
which means that the change of the normal component of $${\underline{D}}\mathbf {v}$$ is positive, if the normal component of $$\mathbf {v}$$ is increased. The second term in Eq. (2.34) can therefore be viewed as a feedback term that aims at reducing $$\mathbf {n}\cdot {\underline{D}}\nabla u$$ to zero by changing $$\mathbf {n}\cdot \nabla u$$ appropriately. This can in principle mean that $$|\mathbf {n}\cdot \nabla u|$$ is increased (if $${\underline{D}}$$ is anisotropic), but zero flux will always be reached at finite values, as indicated above, depending on the parallel components of $$\nabla u$$. As the phase field transition length $$\xi $$ tends to zero, the “feedback strength” is increased. Applying the same strategy for the derivative of $$u$$ parallel to the boundary yields equations for $$\mathbf {p}\cdot \nabla u=\partial u/\partial p$$ identical to Eq. (2.34), only without the last term. Therefore, this component is bounded by the maximum principle.

In summary, the gradient $$\nabla u$$ will thus be bounded independent of the steepness of the phase field at the boundary of $$\mathcal D$$ (at least in the approximations used here), which renders the line of thoughts after Eq. (2.23) correct. In addition to verifying the validity of the phase field approach, the above reasoning contains some useful information about the mechanisms, by which the no-flux boundary condition is achieved: Most importantly, one ambiguity of the ansatz (2.19) was removed. It is now clear why simply multiplying $${\underline{D}}$$ by the phase-field $$\phi $$ does not lead to the desired results, although Eq. (2.23) apparently suggested this. If this was done, Eq. (2.29) would contain $$\partial \phi /\partial n$$ instead of $$\partial \ln \phi /\partial n$$ (in addition to a factor $$\phi $$ before the diffusion term, which would, however, not have any impact). Replacing $$\ln \phi $$ by $$\phi $$ in Eq. (2.34) would destroy the gradient reduction effect of the second term. Therefore, it would not be possible to conclude that the gradient is bounded.

Furthermore, when the constant $$C$$ in Eq. (2.32) is interpreted as the upstroke velocity, as suggested, it can be seen from Eq. (2.23) that the accuracy of the method depends on this maximum upstroke velocity. This was also noted in Ref. [26], without, however, taking into account the possibility that this might be caused by the phase field method itself. This has now been ruled out by the line of thoughts above. An additional source of inaccuracy stems from the steepness of gradients arriving at the interface region from inside the domain $$\mathcal D$$. Equation (2.34) indicates that the boundary condition of the original problem formulated for $$\mathcal D$$ has been replaced by a dynamical process that aims at reducing $$\mathbf {n}\cdot ({\underline{D}}\nabla u)$$ to zero. As this does not happen instantaneously, larger initial gradients mean larger remainders of $$\mathbf {n}\cdot ({\underline{D}}\nabla u)$$ after the reduction process. In excitable media reaction-diffusion systems, the upstroke velocity and the steepness of gradients are, of course, closely related.

Finally, Eq. (2.29) together with the information illustrated in Fig. 2.2c indicates that the velocity, with which the quantity $$u$$ is transported outwards, increases steadily in the interface region and becomes $$\approx $$ $$ 1/\xi $$ after some distance. The smaller $$\xi $$, the less information can therefore travel back from $$\mathcal D_c\backslash \mathcal D$$ to $$\mathcal D$$. The only process available for this is the diffusion term in Eq. (2.29), which does not depend on $$\phi $$ and is therefore arbitrarily slow. This means that the auxiliary (arbitrary) boundary condition applied for numerical reasons at the boundary $$\partial \mathcal D_c$$ of the auxiliary domain becomes less and less important. This fact is useful for the numerical implementation of the phase field method (see Sect. 2.2.3).


2.1.5.4 Extension to Other Boundary Conditions


Although the phase-field method was designed by Fenton et al. to implement no-flux boundary conditions on arbitrary geometries, with the information above, it can be extended to a more general type of Neumann boundary conditions of the type


$$\begin{aligned} \mathbf {n}\cdot {\underline{D}}\nabla u = \mathbf {n}\cdot \mathbf {E}, \end{aligned}$$

(2.36)
where $$\mathbf {E}$$ is a fixed vector. This is a special case of Eq. (2.16), when $${\underline{\sigma }}_o\nabla v_o$$ is a constant electric field vector. It is seen immediately from the line of thoughts after Eq. (2.22) that the corresponding phase field replacement in this case is


$$\begin{aligned} \nabla \cdot ({\underline{D}}\nabla u) \quad \rightarrow \quad \frac{1}{\phi }\nabla \cdot (\phi ({\underline{D}}\nabla u - \mathbf {E})). \end{aligned}$$

(2.37)
With this substition, Eq. (2.21) is not altered for regions where $$\phi \approx 1$$, since $${-\nabla \cdot \mathbf {E}=0}$$. In the integrals of Eq. (2.23), however, the additional term appears, leading to $$\mathbf {n}\cdot ({\underline{D}}\nabla u-\mathbf {E})=0$$ in Eq. (2.25), which is the desired boundary condition. The same result can be obtained from Ref. [29] by starting from Eqs. (2.28) and (2.37) and substituting the correct boundary condition.


2.1.6 Models


So far, in this chapter, the cross-membrane current $$I_\text {ion}(V_m,\mathbf {h})$$ in, e.g., Eqs. (2.4), (2.9) or (2.13) has been treated as an unknown quantity whose dynamics is determined by the membrane potential $$V_m$$ and an additional set of equations for the local (i.e. spatially uncoupled) variables $$\mathbf {h}$$. As described in Sect. 2.1.1, one approach to construct this term is to empirically determine the contributions of all different ion species and membrane compounds that lead to transmembrane currents. Indeed, many models based on this bottom-up approach exist, although they are by definition based on incomplete data. A comprehensive list of physiological models, ordered by tissue type (atrial/ventricular) and species can be found in Ref. [31]. One advantage of this approach is that model parameters can directly be related to physiological properties, which enables investigations of genetic defects and facilitates the adaption of the model for disease conditions. However, a major drawback of these models is their computational cost, due to the large number of local variables used. One extreme example is the Wang-Sobie model for neonatal rat ventricular myocytes [15] with 35 dynamic variables and 74 parameters. Although the large number of parameters might be beneficial for the reasons mentioned above, it is at the same time hindering generic approaches such as parameter scans or fitting procedures. Also, for most physiological models, sensitivity of the dynamics to individual parameters is not documented, which further complicates the analysis.

On the other end of the spectrum, generic models are created using a top-down approach to the model construction problem. Rather than incorporating every possible detail of the system for which information exists, they try to mimic the basic properties of cardiac tissue (or even more generic: of excitable media) using the fewest number of variables and parameters possible to achieve a particular degree of realism. A very generic model for excitable media is the Barkley model [32], which consists of only two variables and three parameters which determine the basic properties of the medium such as the excitability, excitation threshold and duration. The Barkley model will be used in this work for numerical simulations, when values of physiological parameters are not important or unknown and when results are to be obtained for the large class of excitable media. In all other simulations, the somewhat more detailed three-variable Fenton-Karma model [33] will be used, which has 14 parameters that enable it to reproduce important aspects of the dynamics of more realistic cardiac models.

A328193_1_En_2_Fig4_HTML.gif


Fig. 2.4
Barkley model dynamics. Four attempts, marked by “$$\times $$”, to induce an excitation in the Barkley model (without diffusion) by instantaneously setting the model variable $$u$$. In the first case, the stimulation is subthreshold. The second and third stimulation produce almost identical excitations, whereas the fourth stimulation occurs at a time when excitability is not fully recovered yet, as $$v$$ has not returned to zero. Because of this refractoriness, the stimulation fails to produce an excitation. Model parameters: B1 (see Table A.1, p. 184)


2.1.6.1 The Barkley and the Bär-Eiswirth Model


This two-variable, generic model was developed by Barkley in the beginning of the 1990s in an effort to understand spiral wave dynamics in generic excitable media. Although already used in Ref. [32], it was newly introduced in its most efficient numerical form in Ref. [34]. Its equations are


$$\begin{aligned} \frac{\partial u}{\partial t}&=\frac{1}{\epsilon }u(1-u)\left( u-\frac{v+b}{a}\right) +\nabla ^2 u\end{aligned}$$

(2.38a)



$$\begin{aligned} \frac{\partial v}{\partial t}&=g(u,v), \end{aligned}$$

(2.38b)
with the parameters $$\epsilon $$, $$a$$ and $$b$$. For the original Barkley model, $$g(u,v)=u-v$$. The fast variable $$u$$ can be identified as a normalized $$V_m$$ in the mono-domain formalism of Sect. 2.1.3, while the local part on the right hand side of the $$u$$-equation represents the cross-membrane current term $$-I_\text {ion}/C_m$$. There is only one local variable, such that $$\mathbf {h}$$ is actually a scalar and $$h=v$$. The diffusion tensor is normalized to $${\underline{D}}=\text {Id}$$. An example of the local dynamics of the model, i.e. without the diffusion term, following a series of stimulations is shown in Fig. 2.4. It can be understood from the $$(u,v)$$ phase space structure, which is typical for simple excitable systems5 and is depicted in Fig. 2.5. It has a stable fixed point at $$(0,0)$$, whose basin of attraction is the whole square $$[0,1[\times [0,1]$$, but trajectories starting a considerable distance6 right of the nullcline $$u=(v+b)/a$$ make a large detour close to $$u=1$$, $$v=1$$ and $$v=0$$ before returning to the fixed point. In contrast, trajectories starting close to $$(u,v)=(0,0)$$ return directly. The first case corresponds to action potential following a superthreshold depolarization, whereas the second is a subthreshold depolarization not strong enough to evoke an action potential. The meaning of the three parameters in terms of the characteristics of excitable media can readily be read off the phase space diagram: $$b/a$$ is the excitation threshold during rest, whereas $$\epsilon $$ determines the time scale of $$u$$ compared to $$v$$. Its main impact is on the upstroke velocity of an action potential, which also determines the conduction velocity in the spatially extended case in conjunction with diffusion. $$a$$ on its own mainly influences the action potential duration, as it determines the slope of the above-mentioned nullcline and thus the value of $$v$$ at which $$u$$ starts to decrease.

A328193_1_En_2_Fig5_HTML.gif


Fig. 2.5
Barkley model: local phase space. Shown are two stimulation attempts from the fully recovered state: one is subthreshold (cyan), causing the system to return immediately to the stable fixed point at (0,0), whereas the other one is far enough in the $$\dot{u}>0$$” src=”/wp-content/uploads/2016/11/A328193_1_En_2_Chapter_IEq347.gif”></SPAN> region to produce an excitation (<SPAN class=EmphasisTypeItalic>blue</SPAN>). <SPAN class=EmphasisTypeItalic>Arrows</SPAN> on the trajectories mark equidistant time points, indicating the fast dynamics of <SPAN id=IEq348 class=InlineEquation><IMG alt=$$u$$ src= compared to $$v$$. The tilted $$u$$ nullcline has the equation $$v=au-b$$ and is responsible for the refractoriness shortly after an excitation, when still $$v>0$$” src=”/wp-content/uploads/2016/11/A328193_1_En_2_Chapter_IEq352.gif”></SPAN>, but <SPAN id=IEq353 class=InlineEquation><IMG alt=: in this case, the distance to the right necessary to cross the $$u$$ nullcline is considerably increased. When $$v>a-b$$” src=”/wp-content/uploads/2016/11/A328193_1_En_2_Chapter_IEq355.gif”></SPAN>, the system is in absolute refractoriness, because then for every <SPAN id=IEq356 class=InlineEquation><IMG alt= it is $$\dot{u}<0$$. Model parameters: B1 (see Table A.1, p. 184)

Despite its simplicity, the Barkley model (now, again, including the Laplacian term) supports a large number of different wave patterns, including pulses, plane waves and both stationary and meandering spiral waves. However, there are no (homogeneous) parameters for which the Barkley displays spiral wave breakup or other types of chaos. It was shown by Bär and Eiswirth that turbulence can be observed by delaying the increase of the slow variable $$v$$. This can be done by redefining $$g$$ in Eq. (2.38b) as follows [37]:


$$\begin{aligned} g(u,v)={\left\{ \begin{array}{ll} -v &{} u<\frac{1}{3}\\ 1-6.75u(u-1)^2 -v&{} \frac{1}{3} \le u \le 1\\ 1 -v &{} u \ge 1 \end{array}\right. } \end{aligned}$$

(2.39)
This modification will be referred to as the Bär-Eiswirth model in the rest of this thesis.


2.1.6.2 The Fenton-Karma Model


The equations for this generic model were developed by Fenton and Karma in 1998 with the aim of capturing the most important features specific to cardiac models (as opposed to general excitable media). Consequently, although the model has only three variables, it is considerably more tightly connected to the physiological mechanisms of excitability in the cardiac muscle. The model considers three cross-membrane currents, named fast inward (fi), slow outward (so) and slow inward (si) current. These roughly correspond to the sodium, potassium and calcium currents, respectively (cf. Sect. 1.​2.​1). However, the complexity of different ion channel types and other contributions to membrane potential dynamics has been condensed into these net currents and their activation and inactivation. Therefore, the parameters (most of which are timescales and thresholds) have a physiological meaning in that they reflect typical or net properties of these summarized currents, but nevertheless they cannot be interpreted as descriptions of the microscopic components responsible for membrane potential dynamics.

The equations of the Fenton-Karma model are


$$\begin{aligned} \frac{\partial u}{\partial t}&=\nabla \cdot {\underline{D}}\nabla u-J_\text {fi}(u,v)-J_\text {so}(u)-J_\text {si}(u,w)\end{aligned}$$

(2.40a)



$$\begin{aligned} \frac{\partial v}{\partial t}&= \Phi (u_c-u)(1-v)\left( \frac{\Phi (u-u_v)}{\tau _{v1}^-} + \frac{\Phi (u_v-u)}{\tau _{v2}^-}\right) - \Phi (u-u_c)\frac{v}{\tau _v^+}\qquad \end{aligned}$$

(2.40b)



$$\begin{aligned} \frac{\partial w}{\partial t}&= \Phi (u_c-u)\frac{1-w}{\tau _w^-}-\Phi (u-u_c)\frac{w}{\tau _w^+} \end{aligned}$$

(2.40c)
where the three cross-membrane currents are defined as


$$\begin{aligned} J_\text {fi}(u,v)&= -\frac{v}{\tau _d}\Phi (u-u_c)(1-u)(u-u_c)\end{aligned}$$

(2.41a)



$$\begin{aligned} J_\text {so}(u)&= \frac{u}{\tau _0}\Phi (u_c-u)+\frac{1}{\tau _r}\Phi (u-u_c)\end{aligned}$$

(2.41b)



$$\begin{aligned} J_\text {si}(u,w)&= -\frac{w}{2\tau _{si}}(1+\tanh [k(u-u_c^{si})]). \end{aligned}$$

(2.41c)
The sum of these currents corresponds to the term $$-I_\text {ion}/C_m$$ of Sect. 2.1.3. Again, $$u$$ is a normalized membrane potential $$V_m$$ and on the order 0 to 1, such that one unit of $$u$$ roughly corresponds to a 100 mV membrane potential change. The vector of local variables for the Fenton-Karma model therefore is $$\mathbf {h}=(v,w)$$.

In its original published form, the function $$\Phi (x)$$ in Eqs. (2.40) and  (2.41) is simply the step function
Tags:
Nov 21, 2016 | Posted by in CARDIOLOGY | Comments Off on Methods

Full access? Get Clinical Tree
Get Clinical Tree app for offline access