2.4 Numerical Simulations
2.4.4 Boundary Conditions
In section 2.4.1 on page 44 the underlying equations were derived, including boundary conditions for the boundaries between the cardiac tissue and an outer medium (e.g. a bath). In the absence of an external electrical field, a Neumann boundary condition (no-flux) was derived for the membrane potential Vm (Eq. (2.46)). Since secondary variables of the ionic cell model are not diffusive (they are only local), the boundary condition only needs to be specified for the membrane potential. In the following, we present how the no-flux boundary condition for the membrane potential is implemented for two cases. In the first part we discuss how it is defined for the outer boundary of the simulation grid, which is used when the simulated cardiac tissue occupies the whole simulation domain.
Furthermore, when we consider cardiac tissue of a non-trivial geometry embedded into a bath, we use the so called “phase field method”, in order to establish the no-flux boundary condition between the arbitrarily shaped cardiac tissue and the surrounding bath. While numerical simulations of arbitrarily shaped simulation domains can be solved in a more straightforward way with “finite element methods” (FEM), the implementation demands a more sophisticated approach in the case of finite differences. The phase field method is shown in the second part.
Outer Boundaries of the simulation grid
The calculation of the diffusive part, using e.g. the nine-point stencil (Eq. (2.66)) from the previous section, can not be performed in the conventional way at the grid nodes at the edges of the grid, since no neighboring points are available here which are required for the application of the stencil. Here, additional information need to be included, which specify the type of boundary conditions. For the desired no-flux boundary condition required for the membrane potentialVm, so called “ghost points” are used, which are (imaginary) added at the edges of the grid. Figure 2.23 depicts this technique for the one-dimensional case.
Figure 2.23: The use of “ghost points” in the context of the numerical realization of no-flux boundary conditions. In the (exemplary) one-dimensional case, a ghost point with a new grid node is added at the edge of the grid. The virtual membrane potential of this node is denoted byVm0. In the case of a no-flux boundary conditionVm0 takes the value ofVm2, in order to preserve a derivative of zero atVm1.
For the calculation of the diffusive part of the underlying equations for the membrane potential at the edge of the grid, denoted byVm1, a ghost point is added, which comprises an additional grid node with the virtual membrane potential at this point Vm0. In case of the no-flux boundary condition, the value ofVm0 is obtained by mirroringVmalong the edge, thus,Vm0 =Vm2. In this way, the no-flux boundary condition can be achieved by a derivative of zero at the grid pointVm1. Thus, at the edges of the simulation grid the general expression for the diffusive term in one dimension ∆Vmx ≈ Vmx−1−2hV2mx+Vmx+1
x reduces to
∆Vm1 ≈ Vm0 −2Vm1 +Vm2
h2x = 2Vm2−2Vm1
h2x . (2.68)
The subsequent generalization of this technique to two or three dimensional simulation domains is straightforward.
Inner Boundaries Using the Phase Field Method
When simulating arbitrarily shaped pieces of cardiac tissue (e.g. a whole heart in a sur-rounding bath), the use of ghost points is not feasible, due to the complex shape of the boundary. For this objective, the use of the “Phase Field method” is a very elegant way, to implement no-flux boundary conditions while keeping the finite difference scheme. The central object of this method is a phase field φ(r), which is defined at every grid node of the simulation domain, and has values between zero and one. In the simulations of this the-sis, the phase field distinguishes between cardiac tissue (φ= 1) and the surrounding bath (φ= 0), and interpolates smoothly the boundaries between the two domains, with a certain transition zone with the width ξ. X. Li et al. showed, that if this phase field is included into the differential equations in a specific way, the desired no-flux boundary condition is automatically (approximately) fulfilled [79]. In the simulations performed in this thesis, the phase field enters the monodomain equation in the following way:
φ∂Vm
∂t =∇ ·* Dφ∇V* m−φIIon
Cm . (2.69)
Flavio H. Fenton et al. demonstrated how this approach can be used for modeling
electri-2.4. Numerical Simulations
Figure 2.24: The Phase Field method. Subplot(a) shows an exemplary phase field along a cross section of a realistic heart geometry shown in (b). The section of the phase field shown in (a), corresponds to a cross section through the ventricular wall, marked by a red line in (b). Subplot (c) explains how the width of the transition zone between cardiac tissue and bath can be determined by the choice of the ξ parameter (ξ = 0.075 (diamonds), ξ = 0.05 (circles) and ξ = 0.075 (stars), respectively). Reprinted from [80], with the permission of AIP Publishing.
cal wave propagation in realistic heart geometries [80] and used the following differential equation for the creation of the phase field:
∂φ
∂t =ξ2∆φ−∂G(φ)
∂φ , (2.70)
G(φ) =−(2φ−1)2
4 +(2φ−1)4
8 . (2.71)
Figure 2.24 shows an example of how the phase field method is applied using a realistic heart geometry. The phase field is shown in subplot (a), where the smoothing procedure (solving Eq. (2.70)) has already been performed. The representative section shown in (a) is marked by a red line in (b), which shows a cross section of the heart geometry, with the left and right ventricle. That means, the value around zero atx= 0 cm is located outside the tissue.
When passing the epicardium, the phase field increases smoothly to one, before afterwards going back to zero, when leaving the ventricular wall. Subplot (c) shows a magnification of the transition zone. The phase field is plotted here using different values for ξ, which clarifies how the width of the transition zone can be modulated with the choice ofξ. In practice, the phase field φ is initialized with the value zero (surrounding bath) and one (cardiac tissue). Using this initial field, Eq. (2.70) is solved numerically, (with the diffusive term modeled as described in the previous section). In the study of this thesis, where the phase field method was used (section 3.3 on page 102), Eq. (2.69) was solved with the parameters dt= 0.1,h= 0.5 andξ = 0.5.