Focus on Excitable Media and Cardiac Electrophysiology
3d Space-time Adaptive Mesh Refinement. I have for the first time implemented a space-time adaptive mesh refinement algorithm for excitable media, and for models of cardiac electrical activity in particular. My adaptive mesh refinement algorithm is an extension of an AMRA developed by Marsha Berger and collaborators to integrate hyperbolic systems of conservation laws, such as the Euler equations of fluid dynamics. The algorithm is founded on the use of Cartesian grids and approximates a given continuous field on a set of nested locally-uniform patches of d-dimensional Cartesian meshes in a d-dimensional Cartesian box. Every grid patch is defined as a separate data structure, independent of other patches, and has associated with it a level of resolution l. Grid points align from one level to the next, but because the grid spacing varies and the fields are cell-centered, the data points on different levels never align. When combined, the multiple grid levels give an effectively non-uniform grid.
Each level l has its own spatial resolution Δxl and time step Δtl. Although multiple grid patches may exist at a given level of spatial resolution, the same time step is used for all of them to facilitate synchronizing data at different resolutions. The time step Δtl is an integer multiple rt of the time step used for the next coarser level. The ratio of Δtl to Δxl2 for all levels is fixed, which allows the same explicit difference scheme to be stable on all grid patches.
The AMRA assumes that some explicit finite difference scheme (specified by the user) has been chosen to represent both space and time derivatives. Each grid patch is defined separately and maintains its own solution vector, so that grid patches can be integrated independently of other patches, except for the determination of boundary data. Grids of different spatial resolutions are integrated from coarse to fine levels to ensure that internal boundary data for fine grids always can be interpolated from data already computed on coarser grids. Communication among grids is limited to two procedures: updating values at coarse cells with averaged values from the finer cells they overlay and determination of internal (non-physical) boundary values for fine grid patches, which either are provided directly from neighboring grid patches at the same level (if available) or are interpolated from neighboring grid patches at the next coarsest level.
The power of the AMRA arises from its ability to refine or to coarsen the representations of fields automatically and efficiently by varying the number of grid points locally. Grid patches with higher resolution in space and time are created when an estimate of the local truncation error on a coarser mesh exceeds a specified tolerance and are deleted when no longer required.
Figure 1. Representative spiral wave breakup dynamics of the Luo-Rudy-I model using space-time adaptive mesh refinement with three levels of resolution. (A) Spiral wave propagation and breakup. Excited tissue is shown in blue and red and quiescent tissue in yellow. (B) Corresponding grid structure using three levels of space-time resolution, with white indicating coarse grids, yellow intermediate, and green fine. Note that the fine grids closely follow the wave fronts.
Phase-field Method for Irregular-shaped Domains. Complex geometrical structures such as realistic models of the heart's atria or ventricles frequently are handled using finite element or finite volume methods. The use of finite elements to describe the irregular shapes of these models increases the complexity of the numerical scheme and makes parallelization more difficult. Together with my collaborators, I have adapted a technique called the phase-field method used in crystal growth physics to handle the irregular boundary conditions related to cardiac electrophysiology applications. This method provides several benefits, including avoidance of the complex and time-consuming task of three-dimensional mesh generation and, by allowing moving boundaries, anticipation of the need to consider boundary motion when electro-mechanic coupling effects are incorporated in the future.
To apply this method to cardiac electrical propagation or to other reaction-diffusion equations, the diffusion tensor D is multiplied by a phase matrix, defined initially to be the identity where tissue is present and 0 where tissue is absent, and then interpolated smoothly between these values on tissue boundaries. This procedure automatically gives the required zero-flux (Neumann) boundary condition at the cost of slightly blurring boundaries. We have verified that the effective blurring necessary to reproduce the correct solutions is comparable to the error obtained experimentally when measuring the tissue edges. Figure 2 demonstrates the accuracy with which the method handles irregular boundaries in four different cases by comparing the phase field results with those obtained using exact representations of the boundaries. In all cases, the 3-variable cell model of Fenton and Karmas was used with parameters set to reproduce the modified BR model in a mesh of 200x200 elements. Figure 2A and 2B show isochrones of a propagating wave front initiated in the center (A) and on one edge (B) of a quarter-annulus. The solution obtained using the phase field method is compared with the solution obtained using exact representations of the boundaries in polar coordinates. Black contours represent the reference solution calculated in polar coordinates (plotted in the x-y plane) and red contours represent the solution using the phase field in Cartesian coordinates using dt=0.5ms, dx=dy=dr=0.025cm, and r*dΘ between 0.011 and 0.039cm, depending on the radius r, to ensure adequate spatial resolution. Figure 2C is the same as in 2B except that the structure is more complicated since a hole is included in the center. Figure 2D shows propagation on a square with a rectangular hole in a domain with an anisotropy ratio of 5:1 with fibers oriented at an angle of –53 degrees. Both simulations were obtained using Cartesian coordinates, one using direct zero-flux boundary conditions on the edges and the other using the phase field method. In all cases, excellent agreement is obtained between the two methods, which justifies our use of this method for handling complex boundary conditions.
Figure 2. Isochrones of a propagating wave front at 20-ms intervals on four different domains using the phase field method in Cartesian coordinates (red) and a reference solution (black). (A) Propagation initiated off-center in a quarter-annulus domain. (B) Propagation initiated at the right corner. (C) Same as (B) but with a hole in the domain. (D) Propagation initiated on one side of a square domain with a rectangular hole and conduction anisotropy at an angle. Reference solutions for (A-C) are calculated in polar coordinates and plotted in the x-y plane, while the reference solution for (D) is calculated in Cartesian coordinates.
Adaptive Mesh Refinement in Irregular-shaped Domains. By combining the phase-field method with space-time adaptive mesh refinement, it is possible to extend the adaptive algorithm described above to handle domains with irregular shapes. To do so, it is necessary to compute the phase field at the finest spatial resolution to be used and then to use coarsened representations at the other spatial resolutions used, so that all levels of mesh have access to the phase field at their resolution. After this step, the phase field is added to the Laplacian term as described above. Figure 3 shows an example of the voltage and corresponding three-level grid structure in a 2D slice of the canine ventricular anatomy 10 ms after application of a small stimulus. Note that the finer grid structure closely follows the expanding wave front.
Figure 3. Adaptive mesh refinement in a 2-d slice of canine ventricles. (A) Voltage plot 10 ms after a stimulus applied in the lower left portion of the tissue. Higher voltages are shown in blue. (B) Three-level grid structure corresponding to the voltage plot in (A). Coarse resolution regions are shown in white, intermediate in yellow, and fine in green.