Numerical Algorithms in
Computational Biology:
Focus on Excitable Media and
Cardiac Electrophysiology
3d Spacetime Adaptive Mesh Refinement. I have for the first time implemented
a spacetime 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
locallyuniform patches of ddimensional
Cartesian meshes in a ddimensional
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 cellcentered, the data
points on different levels never align. When combined, the
multiple grid levels give an effectively nonuniform grid.
Each level l
has its own spatial resolution Δx_{l} and time
step Δt_{l}. 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 Δt_{l}
is an integer multiple r_{t} of the time step used for the next coarser
level. The ratio of Δt_{l} to Δx_{l}^{2}
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 (nonphysical) 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 LuoRudyI model using spacetime
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 spacetime
resolution, with white indicating coarse grids, yellow intermediate, and
green fine. Note that the fine grids closely follow the wave fronts. 
A 
B 
Phasefield Method for Irregularshaped
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 phasefield 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 timeconsuming task of
threedimensional mesh generation and, by allowing moving boundaries,
anticipation of the need to
consider boundary motion when electromechanic coupling effects are
incorporated in the future.
To apply this method to cardiac electrical
propagation or to other reactiondiffusion 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 zeroflux (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 3variable 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
quarterannulus. 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 xy 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 zeroflux 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 20ms intervals on four different domains using the phase field
method in Cartesian coordinates (red) and a reference solution (black).
(A) Propagation initiated offcenter in a quarterannulus 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 (AC) are
calculated in polar coordinates and plotted in the xy plane, while the reference solution for (D) is
calculated in Cartesian coordinates. 
A

B

C

D

Adaptive Mesh Refinement in Irregularshaped Domains. By combining the phasefield method with spacetime 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 threelevel 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 2d 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) Threelevel 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