Wenjun Ying

Hybrid Finite Element Method for Describing the Electrical Response of Biological Cells to Applied Fields

A novel hybrid finite element method (FEM) for modeling the response of passive and active biological membranes to external stimuli is presented. The method is based on the differential equations that describe the conservation of electric flux and membrane currents. By introducing the electric flux through the cell membrane as an additional variable, the algorithm decouples the linear partial differential equation part from the nonlinear ordinary differential equation part that defines the membrane dynamics of interest. This conveniently results in two subproblems: a linear interface problem and a nonlinear initial value problem. The linear interface problem is solved with a hybrid FEM. The initial value problem is integrated by a standard ordinary differential equation solver such as the Euler and Runge-Kutta methods. During time integration, these two subproblems are solved alternatively. The algorithm can be used to model the interaction of stimuli with multiple cells of almost arbitrary geometries and complex ion-channel gating at the plasma membrane. Numerical experiments are presented demonstrating the uses of the method for modeling field stimulation and action potential propagation

Efficient Fully Implicit Time Integration Methods for Modeling Cardiac Dynamics

Implicit methods are well known to have greater stability than explicit methods for stiff systems, but they often are not used in practice due to perceived computational complexity. This paper applies the backward Euler (BE) method and a second-order one-step two-stage composite backward differentiation formula (C-BDF2) for the monodomain equations arising from mathematically modeling the electrical activity of the heart. The C-BDF2 scheme is an L-stable implicit time integration method and easily implementable. It uses the simplest forward Euler and BE methods as fundamental building blocks. The nonlinear system resulting from application of the BE method for the monodomain equations is solved for the first time by a nonlinear elimination method, which eliminates local and nonsymmetric components by using a Jacobian-free Newton solver, called Newton--Krylov solver. Unlike other fully implicit methods proposed for the monodomain equations in the literature, the Jacobian of the global system after the nonlinear elimination has much smaller size, is symmetric and possibly positive definite, which can be solved efficiently by standard optimal solvers. Numerical results are presented demonstrating that the C-BDF2 scheme can yield accurate results with less CPU times than explicit methods for both a single patch and spatially extended domains.

A kernel-free boundary integral method for elliptic boundary value problems

This paper presents a class of kernel-free boundary integral (KFBI) methods for general elliptic boundary value problems (BVPs). The boundary integral equations reformulated from the BVPs are solved iteratively with the GMRES method. During the iteration, the boundary and volume integrals involving Greens functions are approximated by structured grid-based numerical solutions, which avoids the need to know the analytical expressions of Greens functions. The KFBI method assumes that the larger regular domain, which embeds the original complex domain, can be easily partitioned into a hierarchy of structured grids so that fast elliptic solvers such as the fast Fourier transform (FFT) based Poisson/Helmholtz solvers or those based on geometric multigrid iterations are applicable. The structured grid-based solutions are obtained with standard finite difference method (FDM) or finite element method (FEM), where the right hand side of the resulting linear system is appropriately modified at irregular grid nodes to recover the formal accuracy of the underlying numerical scheme. Numerical results demonstrating the efficiency and accuracy of the KFBI methods are presented. It is observed that the number of GM-RES iterations used by the method for solving isotropic and moderately anisotropic BVPs is independent of the sizes of the grids that are employed to approximate the boundary and volume integrals. With the standard second-order FEMs and FDMs, the KFBI method shows a second-order convergence rate in accuracy for all of the tested Dirichlet/Neumann BVPs when the anisotropy of the diffusion tensor is not too strong.

A kernel-free boundary integral method for implicitly defined surfaces

The kernel-free boundary integral (KFBI) method is a structured grid method for general elliptic partial differential equations. Unlike the standard boundary integral method, it avoids direct evaluation of volume and boundary integrals, which needs to know analytical expressions for the integral kernels. To evaluate a boundary or volume integral, the KFBI method first solves a corrected interface problem on a structured grid and then the numerical solution on the structured grid is interpolated to get approximate values of the integral at points on the boundary. Selection of control points of the boundary plays a key role in the KFBI method since both the correction for the interface equations and the interpolation with the structured grid based solution involve calculation of tangential derivatives of boundary data while stability and efficiency of the numerical differentiation critically depend on the distribution of control points. This work proposes a new point selection method, based on an overlapping surface decomposition of the boundary, which is implicitly defined by a level set function. The points selected are intersection points of the boundary with the grid lines of an underlying Cartesian grid. By the method, the interpolation stencils can be easily chosen to be locally uniform along a coordinate axis in two space dimensions and locally uniform on a coordinate plane in three space dimensions, which allows efficient numerical differentiation and boundary reconstruction/representation. An additional equilibrating process of boundary data further guarantees stable numerical differentiation. Numerical results demonstrating the method with examples in both two and three space dimensions are presented.

Thresholds for Transverse Stimulation: Fiber Bundles in a Uniform Field

Cable theory is used to model fibers (neural or muscular) subjected to an extracellular stimulus or activating function along the fiber (longitudinal stimulation). There are cases however, in which activation from fields across a fiber (transverse stimulation) is dominant and the activating function is insufficient to predict the relative stimulus thresholds for cells in a bundle. This work proposes a general method of quantifying transverse extracellular stimulation using ideal cases of long fibers oriented perpendicular to a uniform field (circular cells in a 2-D extracellular domain). Several methods are compared against a fully coupled model to compute electrical potentials around each cell of a bundle and predict the magnitude of applied plate potential (Phip) needed to activate a given cell (Phipact). The results show that with transverse stimulation, the effect of cell presence on the external field must be considered to accurately compute Phipact. They also show that approximating cells as holes can accurately predict firing order and Phipact of cells in bundles. Potential profiles from this hole model can also be applied to single cell models to account for time-dependent transmembrane voltage responses and more accurately predict Phipact. The approaches used herein apply to other examples of transverse cell stimulation where cable theory is inapplicable and coupled model simulation is too costly to compute.

A Simple Method for Computing Singular or Nearly Singular Integrals on Closed Surfaces

We present a simple, accurate method for computing singular or nearly singular integrals on a smooth, closed surface, such as layer potentials for harmonic functions evaluated at points on or near the surface. The integral is computed with a regularized kernel and corrections are added for regularization and discretization, which are found from analysis near the singular point. The surface integrals are computed from a new quadrature rule using surface points which project onto grid points in coordinate planes. The method does not require coordinate charts on the surface or special treatment of the singularity other than the corrections. The accuracy is about

Field Stimulation of Cells in Suspension: Use of a Hybrid Finite Element Method

O(h^3)

The Bidomain Model of Cardiac Tissue: From Microscale to Macroscale

, where

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

h

A semi-discrete tailored finite point method for a class of anisotropic diffusion problems

is the spacing in the background grid, uniformly with respect to the point of evaluation, on or near the surface. Improved accuracy is obtained for points on the surface. The treecode of Duan and Krasny for Ewald summation is used to perform sums. Numerical examples are presented with a variety of surfaces.