Numerical Algorithms in
Computational Biology:
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. A B 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. A B Back to Elizabeth Cherry's home page
A
B
C
D
