International Journal of Bifurcation and Chaos

Intracellular cardiac Ca²⁺ handling involves interactions of numerous distinct cellular and macromolecular structures. Such interactions coordinate the complicated behaviors of individual processes into spatial and temporal coherent patterns of Ca²⁺ sparks, oscillations and traveling waves. Understanding the dynamical behaviors and functional roles of intracellular Ca²⁺ handling requires not only detailed experimental information about subcellular structures and dynamics, but also useful tools to integrate and analyze this information. This requires interactions between experimental and simulation results within a virtual cell.

Modeling of protein interactions responsible for cardiac tension development can enhance the understanding of physiological and pathophysiological phenomena of the heart. Principal components of muscular tension development are the proteins actin, myosin, troponin and tropomyosin. The tension is produced by cross-bridge cycling of actin and myosin using adenosine triphosphate as energy source. The cross-bridge cycling is initiated by binding of intracellular calcium to troponin, resulting in configuration changes of tropomyosin.

In this work a hybrid model of protein interactions in cardiac tension development is derived on basis of recent measurements and descriptions on protein level. Dependencies on intracellular calcium concentration, sarcomere stretch and stretch velocity as well as cooperativity mechanisms are incorporated. The model quantifies the tension development by states associated to configurations of the involved proteins. The model enables in conjunction with electrophysiological models of cardiac myocytes the reconstruction of electro-mechanical phenomena. Numerical simulations with the hybrid model were performed, which illustrated the reconstruction of steady state and length switches experiments. The steady state experiments describe the force-cytosolic [Ca²⁺] relationship in intact rat cardiac trabeculae. The length switch experiments provide data on the redevelopment of force after sudden stretch in rabbit right ventricular papillary muscles. Results of the numerical simulations show quantitative agreement with experimental studies.

The hybrid model of cardiac tension development offers interfaces to further models of cardiac electro-mechanics. The hybrid model can be coupled with models of cellular electrophysiology and passive mechanics of myocardium allowing the inclusion of mechano-electrical feedback mechanisms. The hybrid model can be applied to elucidate cooperativity mechanisms, pathophysiological changes and metabolism of tension development.

Computer models have proved to be invaluable to experimental and theoretical research in cardiac electrophysiology, and several environments have been designed to aim at easing the development of such models. Recent advances in computing have made the development of graphical user interfaces much simpler, which has lead to the arrival of new software packages that are used primarily for teaching, not research.

Here, we introduce Cellular Open Resource (COR), a modeling environment that runs under Microsoft Windows^® and that can be used for both research and teaching. It is built around CellML^TM, which provides COR with "out of the box" access to an increasing database of single cell models. Though COR was initially developed to handle cardiac modeling, it may also be applied to other types of reaction (and diffusion) problems. The interface is designed with user friendliness in mind (e.g. an equation viewer can be used to graphically visualize an equation as it would appear in a publication). All the cell models used in a single cell or a multicellular problem are checked for correctness, before being dynamically compiled and converted into machine code, in order to maximize computing efficiency. The computation of a problem is both interactive and event-driven. One can, thus, modify the properties of a problem at any point in time and/or subject to a particular condition.

COR is freely available for academic use from .

The complete virtual engineering of the entire beating heart requires the support of sophisticated visualization tools that can effectively extract meaningful visual information from a combination of data sets and depict complex spatial, temporal and physical properties encoded in the data sets. Traditional surface graphics techniques have not been able to serve such requirements adequately, largely because of the necessity of transforming field-based data sets (simulated computationally and recorded clinically or experimentally) to surface representations, and because of the information loss during the process. Constructive volume geometry is a new computer graphics technique, and like the virtual tissue engineering of the beating heart, it employs field-based data types as its intrinsic primitives. In this paper, we describe the application of constructive volume geometry to the visualization of cardiac anatomy and electrophysiology. We demonstrate the capability of this technique in generating visualizations that can depict multiple structures in a combination of data sets and heterogenous interior structure through effective use of opacity and combinational operators. This is applied to tracking and mapping intramural bundles of muscle fibres, and to map voltage isosurfaces of reentrant activity, and the filaments around which reentrant waves propagate.

An excitation wave in nerve or cardiac tissue may fail to propagate if the temporal gradient of the transmembrane voltage at the front becomes too small to excite the tissue ahead of it. A simplified mathematical model is suggested, that reproduces this phenomenon and has exact traveling front solutions. The spectrum of possible propagation speeds is bounded from below. This causes a front to dissipate if it is not allowed to propagate quickly enough. A crucial role is played by the Na inactivation gates, even if their dynamics are by an order of magnitude slower than the dynamics of the voltage.

The sinoatrial node, the pacemaker of the heart, is a structurally and functionally complex and heterogeneous structure. Histology, immunohistochemistry, electrophysiology and mathematical modeling of the sinoatrial node are reviewed to reveal the structure-function relationships of the sinoatrial node. It is argued that contact between the sinoatrial node and surrounding atrial muscle is limited to permit driving of the atrial muscle by the sinoatrial node but not suppression of the sinoatrial node by the atrial muscle (which is more hyperpolarized). It is argued that a protective conduction block zone on one side of the sinoatrial node is the result of a lack of myocytes. Finally, it is argued that the cellular organization of the sinoatrial node is best described by the gradient model.

Multicellular models of homogeneous and isotropic human atria have been developed by incorporating cellular models of membrane electrical activity of single human atrial myocyte into a parabolic partial differential equation. These models are used to study the rate dependent conduction velocity of excitation wave, vulnerability of tissue to reentry and dynamical behaviors of reentry. Bidomain models were also developed to study the actions of a large and brief external electrical stimulus on wave propagation in human atria. These studies provide basic insights to understand the onset and termination of atrial arrhythmias in the human heart.

The dissipation of the excitation wavefronts is a specific mechanism of propagation failure if the sharp gradient of the transmembrane voltage at the wavefront smears out and spread of voltage becomes diffusive, as the main excitation current becomes inactivated. This is produced by the normal kinetics of the ionic currents underlying the action potential. Here we demonstrate that the dissipation of the excitation wavefront can cause arrhythmia as well as lead to its self-termination. We use Courtemanche et al. model of human atrial action potential to demonstrate how reentry creates dynamic electrophysiologic inhomogeneity of the tissue. Local dissipation of the excitation front causes wave breaks and instantaneous displacement of the tip of the reentry, and the same mechanism can lead to elimination of all wavelets, as the inhomogeneity creates conditions for simultaneous dissipation of their excitation fronts.

Our group has described previously the identification of arrhythmogenic pulmonary veins by rapid local electrical activations during atrial fibrillation. We have now investigated an algorithm for automated computer detection of this phenomenon from catheter electrodes in the upper pulmonary veins and assessed its performance in identifying arrhymogenic veins. Ten patients with persistent atrial fibrillation scheduled for pulmonary vein isolation at this hospital were studied. Electrogram recordings in the upper pulmonary veins were recorded and analyzed. Arrhythmogenic veins were identified by focal activity during sinus rhythm at electrophysiological studies. Recordings were visually assessed by a cardiologist for the presence of rapid repetitive electrical activations during atrial fibrillation. An index of rapid repetitive electrical activity (RREA index), the ratio of the number of activations with cycle lengths in the range 50 ms to 100 ms to the number of activations with cycle lengths in the range 100 ms to 200 ms, was devised to describe the extent of such activity automatically. The index was assessed as a predictor of arrhythmogenic veins. Electrograms from 19 upper pulmonary veins were recorded. Rapid activity was evident in 15 veins by visual manual assessment. The mean (range) automatic RREA index was 0.07 (0 to 0.16) for those identified as having no such activity manually, and 0.83 (0.22 to 1.68) for those identified with rapid activity (p<0.0001). With a threshold of RREA index in the range 0.17 to 0.21, the identification of veins with rapid firing was exactly the same as for manual assessment. Eleven upper pulmonary veins were identified as arrhythmogenic during electrophysiological study, and the identification of these veins by both manual and automatic assessment of rapid repetitive electrical activations gave a sensitivity of 100% (11/11) and specificity of 50% (4/8). A technique for automatic characterization of electrogram cycle length has been demonstrated and could be used online as a tool for identifying candidate sites for pulmonary vein isolation in patients despite persistent atrial fibrillation.

The anatomy of the atrioventricular conduction system was first described nearly a hundred years ago. Since then, it has been an occasional subject of controversy mainly through a lack of adherence to the original definitions based on histology. The gross landmarks for locating the atrial component of the conduction system are found in the right atrium. The components and structure of the system in human are comparable to that found in commonly used laboratory animals. The conduction system is composed of specialized myocytes. Its atrial component, the atrioventricular node, is in contact with atrial myocardium. Having penetrated the atrioventricular insulating plane, the major ventricular bundles are encased in fibrous sheaths that separate the specialized myocytes from the ordinary myocardium. Only at the distal ramifications of the bundle branches do the fibrous sheaths disappear, allowing continuity with ventricular myocardium. Being the only muscular pathway connecting atrial with ventricular myocardium, knowledge of its structure can help in developing potential therapies for some forms of cardiac arrhythmias.

The atrioventricular node (AVN) is a small but critically important component of the cardiac electrical conduction system and is located at the junction between right atrium and right ventricle of the heart. It plays important roles in normal and abnormal impulse propagation. Mathematical models of the conduction properties of the AVN have been made, but detailed in silico reconstruction of AVN electrophysiology lags behind that of other cardiac regions. One important facet of detailed reconstruction of AVN electrical activity is the development of comprehensive, ionic conductance-based models of single cell electrophysiology. With a view to facilitating the construction of such models, this article reviews progress made regarding single AVN cell electrophysiological data during the last decade, predominantly focusing on that derived from morphologically normal AVN cells. Properties and potential roles of a range of currents are discussed: including L-type and T-type calcium currents (I_Ca,L and I_Ca,T respectively), background current, hyperpolarization-activated current (I_f), delayed and transient outward potassium currents (I_Kr and I_to, respectively), sodium–calcium exchanger current (I_NaCa), the sustained inward current (I_st) and acetylcholine-activated potassium current (I_KACh).

We study wave propagation in a recently developed model, which reproduces geometry and fiber orientation in the right and left ventricles of the human heart. The cardiac tissue is represented using the previously developed γ-ionic model for human ventricular tissue using a spatial resolution of 0.5 mm. We simulate three-dimensional reentrant behavior resulting from a single vortex located in the free wall of the right, left ventricles and in the interventricular septum. We found that single reentrant scroll waves can generate V-shaped collision areas and in some cases, epicardial breakthrough patterns. The simulated ECGs of single spiral waves show similarities with monomorphic and polymorphic ventricular tachycardia, depending on the location of the reentrant sources. We model complex activation patterns resembling ventricular fibrillation by simulating the effects of an ATP-sensitive potassium channel opener and find that VF is, in that case, organized by a small number of vortices.

In this work, the authors use computer modeling to theoretically investigate the mechanisms involved in figure-of-eight reentry during acute regional myocardial ischemia, a pattern of excitation which may lead to ventricular fibrillation and sudden cardiac death. For this purpose, a modified version of the Luo–Rudy dynamic model for the action potential and ionic currents has been used, together with a two-dimensional model of the regionally ischemic ventricle. The virtual tissue comprises several realistically dimensioned and located transitional border zones for hyperkalemia, hypoxia and acidosis, simulating the substrate heterogeneity created by acute ischemia. Different types of patterns of excitation following the delivery of a premature stimulus were obtained, including figure-of-eight reentry. Action potentials and selected ionic currents which explain the reentry process are analyzed. The effect of the degree of ATP-sensitive current activation in the vulnerability to reentry is also studied. The results are in accordance with experimental observations, and demonstrate the ability of second-generation mathematical models to analyze and explain the mechanisms involved in ischemic reentry.

A new, anatomically accurate, mathematical model of the right and left porcine ventricular myocardium is described based on measurements of the geometry and fibrous-sheet structure. Passive and active properties of the myocardium are calculated using an orthotropic constitutive law based on the fibrous-sheet structure and a biophysical cellular based model of cardiac contraction. Using Galerkin finite element techniques, the equations of finite deformation are solved to determine deformation and regional wall stress through the heart cycle. The mechanics model is coupled via myocardial wall stress, to a one-dimensional coronary blood flow model embedded in the myocardium. Bidomain electrical activation of the myocardium is also modeled, with ionic current based electrophysiological equations and reaction–diffusion equations based on orthotropic conductivity tensors referred to the fibrous-sheet material axes. Metabolic models are used to couple energy supply to contraction and excitation in the heart, and at the body surface, a framework for quantifying the effect of ischemic heart disease is developed.

Scroll waves of electrical excitation are self-organized reentrant waveforms circling a filament and are a suspected cause of lethal cardiac arrhythmias. The lack of information about the filaments inside the cardiac wall raises the question of what is the dynamics of scroll waves and what equilibrium shape, if any, can we expect for their filaments. Here, we analyze the filament configuration and describe their equilibration in a curved slab of twisted fibers, a reasonably close model of the ventricular wall. We then generalize our findings by the use of the recently developed minimal principle to show that the filaments are predictable; namely, they are aligned with minimal resistance paths (geodesics) determined by the fiber configuration.

Despite its importance, the mechanisms underlying ventricular fibrillation (VF) remain poorly understood. Experimental studies are mainly limited to the study of activity on the heart surface, yet the reentrant mechanisms believed to sustain VF are three-dimensional. In 3D, a reentrant wave rotates around a linear filament, and the properties of these filaments can contribute to the stability of reentry. In this paper we describe how filaments can be detected during reentry in cardiac virtual tissues, and how their length and twist can be extracted from numerical results. As an example, we present some results showing the behavior of reentry in 3D virtual tissues with both idealized slab geometry, and anatomically detailed right ventricular (RV) free wall geometry derived from the Auckland canine ventricles. In the slab geometry reentry was stable even with 120° rotational anisotropy, but with a transmural gradient of action potential duration (APD) to mimic the normal distribution of APD across the ventricular wall, reentry became unstable. In contrast, we found that reentry in the anisotropic RV free wall with uniform APD was unstable, possibly because the rotation of fibers is not uniform as in the slab geometry. With transmural differences in APD the instability was augmented, indicating that transmural differences can act synergistically with other mechanisms to increase instability of reentry. This study shows that numerical models can generate information about the behavior of filaments during reentry in 3D, and highlights the need for a theoretical foundation to explain this observed behavior.

Mathematical models of cardiac anatomy and physics provide information, which help to understand structure and behavior of the heart. Miscellaneous cardiac phenomena can only be adequately described by combination of models representing different aspects or levels of detail. Coupling of these models necessitates the definition of appropriate interfaces. Adequateness and efficiency of interfaces is crucial for efficient application of the combined models.

In this work an integrated model is presented consisting of several models interconnected by interfaces. The integrated model allows the reconstruction of macroscopic electro-mechanical processes in the heart. The model comprises a three-dimensional are of left ventricular anatomy represented as truncated ellipsoid. The integrated model includes electrophysiological, tension development and elastomechanical models of myocardium at levels of single cell, proteins, and tissue patches, respectively.

The model is exemplified by simulations of extracorporated left ventricle of small mammals. These simulations yield temporal distributions of electrophysiological parameters as well as descriptions of electrical propagation and mechanical deformation. The simulations show characteristic macroscopic ventricular function resulting from the interplay between cellular electrophysiology, electrical excitation propagation, tension development, and mechanical deformation.

Effects of cardiac mechanical heterogeneity on the electrical function of the heart are difficult to assess experimentally, yet they pose a serious (patho-)physiological challenge. Here, we present an in silico study of the effects of mechanical heterogeneity on action potential duration (APD) in mechanically interacting muscle regions and consequent effects on the dispersion of repolarization, a well-established determinant of cardiac arrhythmogenesis.

Using a novel mathematical description of ventricular electromechanical activity (virtual muscle), we first assessed how differences in intrinsic contractile properties affect the electrical behavior of cardiac muscle representations. In spite of identical electrophysiological model descriptions in virtual muscle samples, faster muscle models show shorter APD than their slower counterparts. This is a consequence of Ca²⁺-mediated feedback from mechanical to electrical activity in the individual muscle models.

This mechano-electric feedback (MEF) is, of course, significantly more complex in native cardiac tissue, as the heterogeneous muscle elements interact both mechanically and electrically. Cardiac mechanical heterogeneity, in its most reduced form, can be represented by a duplex consisting of two mechanically interacting muscle segments. Our in silico model of heterogeneous myocardium therefore consists of two individual virtual muscles that are mechanically interconnected in-series to form a virtual heterogeneous duplex. During isometric contraction of the duplex (i.e. at constant external length), internal mechanical interactions affect Ca²⁺ handling and APD of muscle elements, resulting in an increased dispersion of repolarization beyond the intrinsic APD differences.

Duplex electromechanical activity is strongly affected by the activation sequence of its elements. Late activation of the faster (subepicardial type) duplex element, postponed by time-lags that correspond to normal transmural activation delays, optimizes duplex contractility and smoothes out intrinsic APD differences, thereby reducing dispersion in repolarization. This smoothing effect is not observed upon delayed activation of the slower (subendocardial type) duplex element. In both settings, changes in repolarization timing follow a nonlinear dependence of APD on activation delay. Furthermore, asynchronous activation of identical elements in a homogeneous duplex causes an impairment of contractile function and increases dispersion of repolarization. This suggests that the normal electrical activation sequence in the heart requires matching mechanical and electrical heterogeneity for optimal cardiac performance.

On the subcellular level, our results suggest that mechanical modulation of Ca²⁺ handling is a key mechanism of MEF in heterogeneous myocardium, which contributes to the matching of local mechanical and/or electrical activity to global hemodynamic demand.

Multichannel magnetocardiography measures the magnetic field distribution of the human heart noninvasively from many sites over the body surface. Multichannel magnetocardiogram (MCG) analysis enables regional temporal differences in the distribution of cardiac magnetic field strength during depolarization and repolarization to be identified, allowing estimation of the global and local inhomogeneity of the cardiac activation process.

The aim of this study was to compare the spatial distribution of cardiac magnetic field strength during ventricular depolarization and repolarization in both normal subjects and patients with cardiac abnormalities, obtaining amplitude measurements by magnetocardiography.

MCGs were recorded at 49 sites over the heart from three normal subjects and two patients with inverted T-wave conditions. The magnetic field intensity during depolarization and repolarization was measured automatically for each channel and displayed spatially as contour maps. A Pearson correlation was used to determine the spatial relationship between the variables.

For normal subjects, magnetic field strength maps during depolarization (R-wave) showed two asymmetric regions of magnetic field strength with a high positive value in the lower half of the chest and a high negative value above this. The regions of high R-wave amplitude corresponded spatially to concentrated asymmetric regions of high magnetic field strength during repolarization (T-wave). Pearson-r correlation coefficients of 0.7 (p<0.01), 0.8 (p<0.01) and 0.9 (p<0.01) were obtained from this analysis for the three normal subjects. A negative correlation coefficient of -0.7 (p<0.01) was obtained for one of the subjects with inverted T-wave abnormalities, suggesting similar but inverted magnetic field and current distributions to normal subjects. Even with the high correlation values in these four subjects, the MCG was able to identify differences in the distribution of magnetic field strength, with a shift in the T-wave relative to the R-wave.

The measurement of cardiac magnetic field distribution during depolarization and repolarization of normal subjects and patients with clinical abnormalities should enable the improvement of theoretical models for the explanation of the cardiac depolarization and repolarization processes.

This paper presents the development of a bidomain model of electrical defibrillation and post-shock arrhythmogenesis of the rabbit heart. The model incorporates a realistic representation of cardiac geometry, fiber structure, and membrane ionic processes. Fibrillation is induced in the rabbit heart with burst pacing. Defibrillation shocks of various strengths are delivered to the model in an attempt to terminate the arrhythmia. Two sets of defibrillation simulations are conducted based on the waveform of the defibrillation shock: monophasic and biphasic. The simulations presented here demonstrate the feasibility of modeling realistically the effect of the defibrillation shock on the heart.

We analyze the asymptotic structure of the Hodgkin–Huxley system of equations, in terms of the concepts of slow manifold and fast foliation, based on Tikhonov's theorem on asymptotics of solutions of slow–fast systems of differential equations. We test Zeeman's conjecture that the jump onset–slow return structure of the action potential in realistic equations of biological excitability may be due to a cusp singularity of the slow manifold with respect to the fast foliation. We find that although the cusp singularity can appear in such equations, the characteristic features in question cannot be reproduced within the Tikhonov scheme and require development of different asymptotic approaches.

Wolff–Parkinson–White syndrome is a disease where an arrhythmia is caused by the ventricles being electrically excited by an additional accessory pathway that links the atria to the ventricles. The spread of the activation wave from this pathway to the ventricles is modeled using a simplified model of Hodgkin–Huxley sodium channel kinetics, in a two ion-channel model. The model is investigated both analytically (using an asymptotic analysis) and numerically, and both methods are shown to give the same result. It is found that for a given width of the accessory pathway, there is a critical sodium channel density needed for the activation wave to spread from the pathway to the tissue. This result provides an explanation for the success of class-I anti-arrhythmic drugs in treating Wolff–Parkinson–White syndrome.

EADs were produced by lowered g_K,s in Luo–Rudy cell models, and could be triggered by the addition of a small noise current. Endogenous EADs, produced in a fully deterministic, or a noise driven, one-dimensional virtual tissue could initiate ectopic (unidirectional or bidirectional) propagation, and retrograde, unidirectional propagation that is the analogue of reentry in a two-dimensional tissue. Such potentially reentrant activity occurred with a low probability, for a very narrow parameter range.

Heterogeneity in the electrical activity across the ventricular wall might result from transmural differences in myocardial membrane ionic current densities. Computational models of endo- and epi-cardial action potentials of guinea-pig myocytes were developed, assuming transmural differences in the ionic current densities and kinetics of i_Kr and i_Ks. The cell models were able to reproduce the characteristics of action potentials of endo and epi cells. One- and two-dimensional models of ventricular tissue were also developed to study the effects of transmural heterogeneity on the propagation of excitation waves across the ventricle wall. It was shown that intercellular coupling reduces the transmural heterogeneity across the ventricle wall.

The vulnerable window in a heterogeneous virtual LRl cardiac tissue, with a linear gradient in G_K, is wider when following propagation down the gradient, towards tissue with longer action potential duration, than when following propagation up the gradient. Spiral wave solutions in a uniform linear gradient in G_K drift, with a velocity component along the gradient of the order of mm/s, towards tissue with a longer APD.

A simple two-variable model is used to replace the formulation of calcium dynamics in the Luo–Rudy ventricular cell model. Virtual ventricular cell and tissue are developed and validated to reproduce restitution properties and calcium-dependent voltage patterns present in the original model. Basic interactions between the membrane potential and Ca²⁺ dynamics in the virtual cell and a strand of the virtual tissue are studied. Intracellular calcium waves can be linked to both action potentials (APs) and delayed afterdepolarizations (DADs). An intracellular calcium wave propagating from the cell interior can induce an AP upon reaching the cell membrane. The voltage and the intracellular Ca²⁺ patterns within the same cell can be highly desynchronized. In a one-dimensional strand of the virtual tissue calcium motion is driven by the AP propagation. However, calcium release can be induced upon certain conditions (e.g. Na⁺ overload of the cells), which results in DADs propagating in the wake of AP. Such propagating DADs can reach the excitation threshold, generating a pair of extrasystolic APs. Collision of a propagating AP with a site of elevated intracellular Ca²⁺ concentration does not affect the propagation under the normal conditions. Under Na⁺ overload local elevation of the intracellular Ca²⁺ leads to generation of an extrasystolic AP, which destroys the original propagating AP.