A Numerical Method for Cardiac Mechanoelectric Simulations

Much effort has been devoted to developing numerical techniques for solving the equations that describe cardiac electrophysiology, namely the monodomain equations and bidomain equations. Only a limited selection of publications, however, address the development of numerical techniques for mechanoelectric simulations where cardiac electrophysiology is coupled with deformation of cardiac tissue. One problem commonly encountered in mechanoelectric simulations is instability of the coupled numerical scheme. In this study, we develop a stable numerical scheme for mechanoelectric simulations. A number of convergence tests are carried out using this stable technique for simulations where deformations are of the magnitude typically observed in a beating heart. These convergence tests demonstrate that accurate computation of tissue deformation requires a nodal spacing of around 1 mm in the mesh used to calculate tissue deformation. This is a much finer computational grid than has previously been acknowledged, and has implications for the computational efficiency of the resulting numerical scheme.     

Activation Dynamics in Anisotropic Cardiac Tissue via Decoupling

Bidomain theory for cardiac tissue assumes two interpenetrating anisotropic media--intracellular (i) and extracellular (e)--connected everywhere via a cell membrane; four local parameters sigma(i,e)(l,t) specify conductivities in the longitudinal (l) and transverse (t) directions with respect to cardiac muscle fibers. The full bidomain model for the propagation of electrical activation consists of coupled elliptic-parabolic partial differential equations for the transmembrane potential upsilon(m) and extracellular potential phi(e), together with quasistatic equations for the flow of current in the extracardiac regions. In this work we develop a preliminary assessment of the consequences of neglecting the effect of the passive extracardiac tissue and intracardiac blood masses on wave propagation in isolated whole heart models and describe a decoupling procedure, which requires no assumptions on the anisotropic conductivities and which yields a single reaction-diffusion equation for simulating the propagation of activation. This reduction to a decoupled model is justified in terms of the dimensionless parameter epsilon = (sigma(i)(l)sigma(e)(t) - sigma(i)(t)sigma(e)(l))/(sigma(i)(l) + sigma(e)(l))(sigma(i)(t) + sigma(e)(t)). Numerical simulations are generated which compare propagation in a sheet H of cardiac tissue using the full bidomain model, an isolated bidomain model, and the decoupled model. Preliminary results suggest that the decoupled model may be adequate for studying general properties of cardiac dynamics in isolated whole heart models.     

Two-dimensional Chebyshev pseudospectral modelling of cardiac propagation

Bidomain or monodomain modelling has been used widely to study various issues related to action potential propagation in cardiac tissue. In most of these previous studies, the finite difference method is used to solve the partial differential equations associated with the model. Though the finite difference approach has provided useful insight in many cases, adequate discretisation of cardiac tissue with realistic dimensions often requires a large number of nodes, making the numerical solution process difficult or impossible with available computer resources. Here, a Chebyshev pseudospectral method is presented that allows a significant reduction in the number of nodes required for a given solution accuracy. The new method is used to solve the governing nonlinear partial differential equation for the monodomain model representing a two-dimensional homogeneous sheet of cardiac tissue. The unknown transmembrane potential is expanded in terms of Chebyshev polynomial trial functions and the equation is enforced at the Gauss-Lobatto grid points. Spatial derivatives are obtained using the fast Fourier transform and the solution is advanced in time using an explicit technique. Numerical results indicate that the pseudospectral approach allows the number of nodes to be reduced by a factor of sixteen, while still maintaining the same error performance. This makes it possible to perform simulations with the same accuracy using about twelve times less CPU time and memory.     

Mathematical Modeling of Electrocardiograms: A Numerical Study

This paper deals with the numerical simulation of electrocardiograms (ECG). Our aim is to devise a mathematical model, based on partial differential equations, which is able to provide realistic 12-lead ECGs. The main ingredients of this model are classical: the bidomain equations coupled to a phenomenological ionic model in the heart, and a generalized Laplace equation in the torso. The obtention of realistic ECGs relies on other important features—including heart–torso transmission conditions, anisotropy, cell heterogeneity and His bundle modeling—that are discussed in detail. The numerical implementation is based on state-of-the-art numerical methods: domain decomposition techniques and second order semi-implicit time marching schemes, offering a good compromise between accuracy, stability and efficiency. The numerical ECGs obtained with this approach show correct amplitudes, shapes and polarities, in all the 12 standard leads. The relevance of every modeling choice is carefully discussed and the numerical ECG sensitivity to the model parameters investigated.     

An Efficient Technique for the Numerical Solution of the Bidomain Equations

Computing the numerical solution of the bidomain equations is widely accepted to be a significant computational challenge. In this study we extend a previously published semi-implicit numerical scheme with good stability properties that has been used to solve the bidomain equations (Whiteley, J.P. IEEE Trans. Biomed. Eng. 53:2139–2147, 2006). A new, efficient numerical scheme is developed which utilizes the observation that the only component of the ionic current that must be calculated on a fine spatial mesh and updated frequently is the fast sodium current. Other components of the ionic current may be calculated on a coarser mesh and updated less frequently, and then interpolated onto the finer mesh. Use of this technique to calculate the transmembrane potential and extracellular potential induces very little error in the solution. For the simulations presented in this study an increase in computational efficiency of over two orders of magnitude over standard numerical techniques is obtained.     

Cardiac propagation simulation.

We have completed a range of membrane-based simulations of action potential propagation in two- and three-dimensional models of ventricular myocardium. The two-dimensional simulations included a bidomain representation of the myocardium which explicitly characterized the component volume conductors in the intracellular, interstitial, and extracellular spaces. With these simulations, we studied the contribution of the extracellular volume conductor to transmural myocardial propagation during depolarization. We also used two-dimensional bidomain simulations to study the effect of the interstitial volume conductor in the setting of planar myocardial depolarization with nominal and extreme tissue conductivities. Our three-dimensional simulations included a monodomain representation of the myocardium which characterized the three component volume conductors as a single lumped conductor. With these simulations, we examined the effects of the intramural rotation of the fiber axes on the timing and pattern of activation. To achieve practical solution times, we extended numerical techniques from previous reports and developed a range of new techniques applicable to this class of problems. Simulations of the depolarization wavefront used the nonlinear Ebihara and Johnson membrane equations for the fast sodium current as the membrane model. Simulations of the full action potential cycle combined the Ebihara and Johnson fast sodium current with the Beeler and Reuter membrane equations. Our results demonstrated that the individual volume conductors and the rotation of fiber axes have unique and identifiable consequences on the electrical activation in models of ventricular myocardium.     

Development of a model for point source electrical fibre bundle stimulation

A model is presented for determining the excitation (transmembrane) potentials on nerve and muscle fibres in a cylindrical bundle from an external point source electrode at the surface and within the preparation. The fibre bundle is considered to be immersed in an infinite isotropic conductive medium and is idealised as an infinitely extending cylinder. This cylinder is initially represented as an isotropic monodomain. A subsequent degree of complexity introduces anisotropy in the monodomain, and finally the bundle is represented as an anisotropic bidomain comprised of the interstitial radial and longitudinal conductivities, the intracellular longitudinal conductivity and the fibre membrane between the two domains. In this latter model, electrical coupling from extracellular to intracellular space is included by means of the bidomain formulation. Computational aspects are discussed, and preliminary results for prescribed conditions are presented.     

The use of spectral methods in bidomain studies.

A Fourier transform method is developed for solving the bidomain coupled differential equations governing the intracellular and extracellular potentials on a finite sheet of cardiac cells undergoing stimulation. The spectral formulation converts the system of differential equations into a "diagonal" system of algebraic equations. Solving the algebraic equations directly and taking the inverse transform of the potentials proved numerically less expensive than solving the coupled differential equations by means of traditional numerical techniques, such as finite differences; the comparison between the computer execution times showed that the Fourier transform method was about 40 times faster than the finite difference method. By application of the Fourier transform method, transmembrane potential distributions in the two-dimensional myocardial slice were calculated. For a tissue characterized by a ratio of the intra- to extracellular conductivities that is different in all principal directions, the transmembrane potential distribution exhibits a rather complicated geometrical pattern. The influence of the different anisotropy ratios, the finite tissue size, and the stimuli configuration on the pattern of membrane polarization is investigated.     

Multicompartment diffusion and reaction modeling system applied to transmitter-mediated photoreception.

An integrated system of programs has been developed with broad applicability to the numerical solution of models with parabolic partial differential equations coupled to ordinary differential equations, as arise, for example, in diffusive transport bulk- or surface-limited by reaction rate processes. The programs have been designed to run optimally in various minicomputer environments and to be as portable as possible. The difference scheme for the parabolic equations is new, and competes favorably with several commonly used implicit schemes. Convergence is proved, and conditions on stability are given. To solve the two sets of coupled diffusion + reaction equations, a new numerical method is developed which digitally filters the second space differences of the diffusion difference equations, making the simple, explicit difference scheme convergent for arbitrary time step size. The method competes favorably with the Cranl-Nicolson scheme in speed and accuracy.The modeling system is applied to certain photoresponsive cells of the Aplysia californica (R₂ giant neuron and the ventral photoresponsive neuron) which hyperpolarize when illuminated, due to an increase of membrane potassium permeability. It has been hypothesized that light releases an internal transmitter from cytoplasmic granules. Three model compartments parallel cellular morphology: the first represents the granule component with bulk-limited diffusion; the second corresponds to the cytoplasm and involves simple diffusion; and the third is situated near the plasma membrane where the transmitter concentration is directly related to membrane conductance. The bimolecular binding succesfully predicts the dynamic non-linearity.     

Physiology Driven Adaptivity for the Numerical Solution of the Bidomain Equations

Previous work [Whiteley, J. P. IEEE Trans. Biomed. Eng. 53:2139–2147, 2006] derived a stable, semi-implicit numerical scheme for solving the bidomain equations. This scheme allows the timestep used when solving the bidomain equations numerically to be chosen by accuracy considerations rather than stability considerations. In this study we modify this scheme to allow an adaptive numerical solution in both time and space. The spatial mesh size is determined by the gradient of the transmembrane and extracellular potentials while the timestep is determined by the values of: (i) the fast sodium current; and (ii) the calcium release from junctional sarcoplasmic reticulum to myoplasm current. For two-dimensional simulations presented here, combining the numerical algorithm in the paper cited above with the adaptive algorithm presented here leads to an increase in computational efficiency by a factor of around 250 over previous work, together with significantly less computational memory being required. The speedup for three-dimensional simulations is likely to be more impressive.     

Mathematical modelling of the spatio-temporal response of cytotoxic T-lymphocytes to a solid tumour.

In this paper a mathematical model describing the growth of a solid tumour in the presence of an immune system response is presented. In particular, attention is focused upon the attack of tumour cells by so-called tumour-infiltrating cytotoxic lymphocytes (TICLs), in a small, multicellular tumour, without necrosis and at some stage prior to (tumour-induced) angiogenesis. At this stage the immune cells and the tumour cells are considered to be in a state of dynamic equilibrium--cancer dormancy--a phenomenon which has been observed in primary tumours, micrometastases and residual disease after ablation of the primary tumour. Nonetheless, the precise biochemical and cellular mechanisms by which TICLs control cancer dormancy are still poorly understood from a biological and immunological point of view. Therefore we focus on the analysis of the spatio-temporal dynamics of tumour cells, immune cells and chemokines in an immunogenic tumour. The lymphocytes are assumed to migrate into the growing solid tumour and interact with the tumour cells in such a way that lymphocyte-tumour cell complexes are formed. These complexes result in either the death of the tumour cells (the normal situation) or the inactivation (sometimes even the death) of the lymphocytes. The migration of the TICLs is determined by a combination of random motility and chemotaxis in response to the presence of chemokines. The resulting system of four nonlinear partial differential equations (TICLs, tumour cells, complexes and chemokines) is analysed and numerical simulations are presented. We consider two different tumour geometries--multi-layered cell growth and multi-cellular spheroid growth. The numerical simulations demonstrate the existence of cell distributions that are quasi-stationary in time and heterogeneous in space. A linear stability analysis of the underlying (spatially homogeneous) ordinary differential equation (ODE) kinetics coupled with a numerical investigation of the ODE system reveals the existence of a stable limit cycle. This is verified further when a subsequent bifurcation analysis is undertaken using a numerical continuation package. These results then explain the complex heterogeneous spatio-temporal dynamics observed in the partial differential equation (PDE) system. Our approach may lead to a deeper understanding of the phenomenon of cancer dormancy and may be helpful in the future development of more effective anti-cancer vaccines.     

Three-dimensional pseudospectral modelling of cardiac propagation in an inhomogeneous anisotropic tissue

Various investigators have used the monodomain model to study cardiac propagation behaviour. In many cases, the governing non-linear parabolic equation is solved using the finite-difference method. An adequate discretisation of cardiac tissue with realistic dimensions, however, often leads to a large model size that is computationally demanding. Recently, it has been demonstrated, for a two-dimensional homogeneous monodomain, that the Chebyshev pseudospectral method can offer higher computational efficiency than the finite-difference technique. Here, an extension of the pseudospectral approach to a three-dimensional inhomogeneous case with fibre rotation is presented. The unknown transmembrane potential is expanded in terms of Chebyshev polynomial trial functions, and the monodomain equation is enforced at the Gauss-Lobatto node points. The forward Euler technique is used to advance the solution in time. Numerical results are presented that demonstrate that the Chebyshev pseudospectral method offered an even larger improvement in computational performance over the finite-difference method in the three-dimensional case. Specifically, the pseudospectral method allowed the number of nodes to be reduced by ≈85 times, while the same solution accuracy was maintained. Depending on the model size, simulations were performed with ≈18–41 times less memory and ≈99–169 times less CPU time.     

Stiffness Analysis of Cardiac Electrophysiological Models

The electrophysiology in a cardiac cell can be modeled as a system of ordinary differential equations (ODEs). The efficient solution of these systems is important because they must be solved many times as sub-problems of tissue- or organ-level simulations of cardiac electrophysiology. The wide variety of existing cardiac cell models encompasses many different properties, including the complexity of the model and the degree of stiffness. Accordingly, no single numerical method can be expected to be the most efficient for every model. In this article, we study the stiffness properties of a range of cardiac cell models and discuss the implications for their numerical solution. This analysis allows us to select or design numerical methods that are highly effective for a given model and hence outperform commonly used methods.     

Effect of a perfusing bath on the rate of rise of an action potential propagating through a slab of cardiac tissue

Experiments show that the rate of rise of the action potential depends on the direction of propagation in cardiac tissue. Two interpretations of these experiments have been presented: (i) the data are evidence of discrete propagation in cardiac tissue, and (ii) the data are an effect of the perfusing bath. In this paper we present a mathematical model that supports the second interpretation. We use the bidomain model to simulate action potential propagation through a slab of cardiac tissue perfused by a bath. We assume an intracellular potential distribution and solve the bidomain equations analytically for the transmembrane and extracellular potentials. The key assumption in our model is that the intracellular potential is independent of depth within the tissue. This assumption ensures that all three boundary conditions at the surface of a bidomain are satisfied simultaneously. One advantage of this model over previous numerical calculations is that we obtain an analytical solution for the transmembrane potential. The model predicts that the bath reduces the rate of rise of the transmembrane action potential at the tissue surface, and that this reduction depends on the direction of propagation. The model is consistent with the hypothesis that the perfusing bath causes the observed dependence of the action-potential rate of rise on the direction of propagation, and that this dependence has nothing to do with discrete properties of cardiac tissue.     

Multiscale modeling of gastrointestinal electrophysiology and experimental validation

Normal gastrointestinal (GI) motility results from the coordinated interplay of multiple cooperating mechanisms, both intrinsic and extrinsic to the GI tract. A fundamental component of this activity is an omnipresent electrical activity termed slow waves, which is generated and propagated by the interstitial cells of Cajal (ICCs). The role of ICC loss and network degradation in GI motility disorders is a significant area of ongoing research. This review examines recent progress in the multiscale modeling framework for effectively integrating a vast range of experimental data in GI electrophysiology, and outlines the prospect of how modeling can provide new insights into GI function in health and disease. The review begins with an overview of the GI tract and its electrophysiology, and then focuses on recent work on modeling GI electrical activity, spanning from cell to body biophysical scales. Mathematical cell models of the ICCs and smooth muscle cell are presented. The continuum framework of monodomain and bidomain models for tissue and organ models are then considered, and the forward techniques used to model the resultant body surface potential and magnetic field are discussed. The review then outlines recent progress in experimental support and validation of modeling, and concludes with a discussion on potential future research directions in this field.     

Segmentation of 3D cell membrane images by PDE methods and its applications

We present a set of techniques that enable us to segment objects from 3D cell membrane images. Particularly, we propose methods for detection of approximate cell nuclei centers, extraction of the inner cell boundaries, the surface of the organism and the intercellular borders—the so called intercellular skeleton. All methods are based on numerical solution of partial differential equations. The center detection problem is represented by a level set equation for advective motion in normal direction with curvature term. In case of the inner cell boundaries and the global surface, we use the generalized subjective surface model. The intercellular borders are segmented by the advective level set equation where the velocity field is given by the gradient of the signed distance function to the segmented inner cell boundaries. The distance function is computed by solving the time relaxed eikonal equation. We describe the mathematical models, explain their numerical approximation and finally we present various possible practical applications on the images of zebrafish embryogenesis—computation of important quantitative characteristics, evaluation of the cell shape, detection of cell divisions and others.     

Nanoparticle Deposition onto Biofilms

We develop a mathematical model of nanoparticles depositing onto and penetrating into a biofilm grown in a parallel-plate flow cell. We carry out deposition experiments in a flow cell to support the modeling. The modeling and the experiments are motivated by the potential use of polymer nanoparticles as part of a treatment strategy for killing biofilms infecting the deep passages in the lungs. In the experiments and model, a fluid carrying polymer nanoparticles is injected into a parallel-plate flow cell in which a biofilm has grown over the bottom plate. The model consists of a system of transport equations describing the deposition and diffusion of nanoparticles. Standard asymptotic techniques that exploit the aspect ratio of the flow cell are applied to reduce the model to two coupled partial differential equations. We perform numerical simulations using the reduced model. We compare the experimental observations with the simulation results to estimate the nanoparticle sticking coefficient and the diffusion coefficient of the nanoparticles in the biofilm. The distributions of nanoparticles through the thickness of the biofilm are consistent with diffusive transport, and uniform distributions through the thickness are achieved in about four hours. Nanoparticle deposition does not appear to be strongly influenced by the flow rate in the cell for the low flow rates considered.     

Extracellular Measurement of Anisotropic Bidomain Myocardial Conductivities. I. Theoretical Analysis

The passive electrical properties of cardiac tissue, such as the intracellular and interstitial conductivities along the longitudinal and transverse axes, have not been often measured because intracellular electrodes are usually needed for these measurements. In this paper, we present a theoretical analysis of two myocardial models developed to estimate these properties by analyzing potentials recorded with a pair of extracellular electrodes while injecting alternating current between another pair of electrodes. First, the cardiac tissue is represented by a standard bidomain model which includes a membrane capacitance; second, this model is modified by adding an intracellular capacitance representing the intercalated disks. Numerical solutions are computed with a fast Fourier transform algorithm without constraining the anisotropy ratios of the interstitial and intracellular domains. We systematically investigate the effects of changes in the bidomain parameters on the voltage-to-current ratio curves. We also demonstrate how the bidomain parameters can be theoretically estimated by fitting, with a modified Shor's r algorithm, the simulated potentials along the longitudinal and transverse axes for different frequencies between 10 and 10000 Hz. An important finding is that the interelectrode distance must be similar to the myocardial space constant so as to obtain frequency dependent measurements. © 2001 Biomedical Engineering Society. PAC01: 8719Nn, 8719Hh, 8716Uv, 0230Uu, 8716Ac