BioMed Research International

Volume 2015 (2015), Article ID 137482, 14 pages

http://dx.doi.org/10.1155/2015/137482

## Adaptive Mesh Refinement and Adaptive Time Integration for Electrical Wave Propagation on the Purkinje System

^{1}Department of Mathematics, MOE-LSC and Institute of Natural Sciences, Shanghai Jiao Tong University, Minhang, Shanghai 200240, China^{2}Departments of Biomedical Engineering and Computer Science, Duke University, Durham, NC 27708-0281, USA

Received 19 November 2014; Accepted 4 February 2015

Academic Editor: Rodrigo W. dos Santos

Copyright © 2015 Wenjun Ying and Craig S. Henriquez. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

A both space and time adaptive algorithm is presented for simulating electrical wave propagation in the Purkinje system of the heart. The equations governing the distribution of electric potential over the system are solved in time with the method of lines. At each timestep, by an operator splitting technique, the space-dependent but linear diffusion part and the nonlinear but space-independent reactions part in the partial differential equations are integrated separately with implicit schemes, which have better stability and allow larger timesteps than explicit ones. The linear diffusion equation on each edge of the system is spatially discretized with the continuous piecewise linear finite element method. The adaptive algorithm can automatically recognize when and where the electrical wave starts to leave or enter the computational domain due to external current/voltage stimulation, self-excitation, or local change of membrane properties. Numerical examples demonstrating efficiency and accuracy of the adaptive algorithm are presented.

#### 1. Introduction

One of the long-recognized challenges in modeling cardiac dynamics [1–4] is developing efficient and accurate algorithms that can accommodate the widely varying scales in both space and time [5]. The electrical wave fronts typically occupy only a small fraction of the domain, are very sharp (in space), and change very rapidly (in time) while the electrical potential, in the region away from the wave fronts, is spatially broad and changes more slowly. With standard numerical methods on uniform grids, very small mesh parameters and very small timesteps must be used to correctly resolve the fine details of the sharp and rapidly changing wave fronts. These discretization parameters are often chosen heuristically and are fixed throughout the simulation, even if conditions change. Adaptive mesh refinement (AMR) methods have been proposed as a solution, in which coarse grids and large timesteps are used in the area where the electrical potential is changing slowly and fine grids and small timesteps are applied only in the region where the sharp electrical waves are located and the action potential changes very rapidly. Using this approach, the numbers of grid nodes and timesteps used with the adaptive algorithm are to some extent optimized. The original AMR algorithm was first proposed by Berger and Oliger for hyperbolic equations [6] and shock hydrodynamics [7]. The methods have been applied to cardiac simulations by Cherry et al. [8, 9] and Trangenstein and Kim [10].

The Berger-Oliger AMR algorithm is a hierarchical and recursive integration method for time-dependent partial differential equations. It starts time integration on a relatively coarse grid with a large timestep. The coarse grid is locally refined further if the computed solution at part of the domain is estimated to have large errors. Better solutions are obtained by continuing time integration on the fine grid with a smaller timestep until both coarse and fine grids reach the same time, called synchronization of levels. The fine grid may be locally refined further and is dynamically changing, which leads to both space and time adaptive algorithm.

The standard implementation of Berger-Oliger’s AMR algorithm uses block-structured grids and assumes that the underlying grids are logically rectangular and can be mapped onto a single index space. While the method has been shown to provide computation and accuracy advantages, it is challenging to apply to domains with complex geometry.

In this paper, we present an AMR algorithm that can be used for unstructured grids. The method is applied in both idealized and realistic tree-like domains similar to that found in the His-Purkinje system. The algorithm is based on that proposed by Trangenstein and Kim [10] where operator splitting technique is used to separate the space-independent reactions part from the linear diffusion part during time integration. Both the reactions and the diffusion parts are integrated with an implicit scheme. This allows larger and adaptive timesteps. As the AMR algorithm naturally provides a hierarchy of multilevel grids, the linear systems resulting from space discretization of the linear diffusion on the adaptively refined grids are solved by a standard geometric multigrid solver. The results show that uniform coarse discretization can lead to conduction failure or changes in dynamics in some parts of the branching network when compared to a uniform fine grid. The results also show that AMR scheme with adaptive time integration (AMR-ATI) can yield results as accurate as the uniform fine grid but with a speedup of 15 times.

The remainder of the work is organized as follows. Section 2 describes the partial differential equations, which model electrical wave propagation on the Purkinje system. Sections 3 and 4, respectively, present the time integration and space discretization for the reaction-diffusion equations. The adaptive mesh refinement and adaptive time integration procedures are outlined in Sections 5 and 6. In Section 7, some simulation results are presented with the AMR-ATI algorithm.

#### 2. Differential Equations

Suppose that we are given a fiber network of the Purkinje system with vertices and edges. See Figure 1 for an idealized branch of the Purkinje system. Let be the number of edges connected to vertices , for each . We call the* degree* of the vertex . A vertex with is called a* leaf-vertex*. Otherwise, it is called a* nonleaf-vertex*. Denote the th edge in the fiber network by . Denote the edges that have vertex as the common endpoint by . Denote by the endpoints of the adjacent edges, which overlap with vertex . The subscript is either 0 or 1 for .