Introduction

Computer models of electrical propagation in the heart represent a powerful tool to understand the mechanism of arrhythmias. Their major drawback is their heavy computational load, which imposes restrictions for realistic modeling. Models are designed as tradeoff between membrane complexity, tissue structure, tissue size and duration of simulation. The most effective computer models of the heart have up to 1'000'000 spatial nodes, thus heavily restraining the size of tissue that can be simulated in 2D or 3D. Stronger restrictions occur when the objective is the realistic simulation of arrhythmia, because of the important amount of cardiac cells  involved (the whole ventricular or atrial myocardium), and the several seconds of time required to reach a fully developed arrhythmic state. Models of conduction in the cardiac tissue based on membrane ionic currents have been developed over the last years, mostly for two dimensional tissue [1], but very few models take into account the geometry of the heart, which is a key factor for the development of arrhythmias. To overcome this problem, models incorporating anatomical 3D properties mainly use simplified membrane kinetics, such as the FitzHugh-Nagumo model [2]. While in the ventricles, reentry is known to be a 3D process involving transmural activity, major atrial activations involve only the thin walls of the atria and therefore the arrhythmic process could be simplified to a 2D problem in a 3D structure [2].

Based on this assumption, we describe in this paper a computer model of the human atria, using a simplified 3D geometry with a realistic atrial topology and size and a physiological cellular model (Beeler-Reuter or Luo-Rudy). The objective of this paper is to show how this anatomical computer model of human atria can contribute to the understanding of atrial arrhythmias and provide a tool for the development of therapeutic strategies.

Methods

Cardiac Tissue Model

The cardiac tissue consists of a grid of cardiac cells interconnected via gap junctions. Electrical propagation in this tissue can be modeled using a cable equation. If we neglect current flow in the extracellular domain, we obtain the monodomain formulation of action potential propagation, which is described by the following reaction-diffusion equation:

(1)

where S_v is the cell membrane surface to cell volume ration, D the conductivity tensor, V_m the transmembrane potential, C_m the membrane capacitance, I_ion the sum of the transmembrane ionic currents and I_stim a stimulus current. All simulations where conducted with C_m = 1 µF/cm² and S_v = 0.24 µm^-1. Intracellular resistivities can be individually programmed in the longitudinal and transversal directions via the conductivity tensor D (D is assumed to be diagonal), allowing introduction of heterogeneities in the tissue. In the presented results, the cardiac tissue is assumed to be globally homogeneous and isotropic. Heterogeneities have been introduced only to model major conductive obstacles like veins and valves or local slow conduction regions (between tricuspid valve and inferior vena cava for example). The sum of the transmembrane ionic currents I_ion in Equation (1) is computed either with the Beeler-Reuter [3] or the Luo-Rudy [4] models. These models are described by non-linear equations involving voltage dependent closing and opening gating variables for the different ion channels. In this paper we present simulations obtained only with the Beeler-Reuter model but similar results have been obtained with the Luo-Rudy formulation for I_ion.

Equation (1) has been discretized  by a finite differences method in both time and space, and solved in two steps: firstly an explicit computation of the membrane ionic currents I_ion, using Runge-Kutta integration and lookup tables to speed up gating variable computation, and secondly a semi-implicit computation of the remaining diffusive terms through a classical ADI scheme. The semi-implicit computation has been preferred to a simpler fully explicit method due to its stability even if a comparison between ADI and explicit methods show that it introduces a computational load increase of 30% approximately.

The choice of time and space discretizations of Equation (1) is the result of a tradeoff between computation speed and accuracy. Computer simulations of electrical propagation and especially reentry simulations impose strong requirements on the spatial discretization, which in turn has a major impact on the computational requirements. Typical signatures of incorrect spatial discretization are [5]: reduced propagation velocity (sometimes leading to artificial block of propagation), artificial anisotropy depending on the grid directions with slower propagation along the grid directions than in the diagonal directions (square waves patterns), and alteration of the reentry dynamics. The two main restrictions to spatial refinements are computation time and memory requirements (and numerical stability if explicit solving schemes were used). Based on these considerations, we have used space steps of 200 µm and fixed time steps of 25 µs. The resistivities (inverse of conductivity) have been varied between 80 Ohm·cm to 800 Ohm·cm, corresponding to propagation velocities ranging from approximately c = 90 cm/s to 30 cm/s. From a computational point of view, a reduction of propagation velocity by a resitivity increase is equivalent to a dilation of both atria while keeping a normal propagation velocity.

Several experiments have been performed to show the realistic behavior of this two-dimensional tissue [6].

Human Atria Model

Based on the assumption that atria are constituted of thin walls [2], we have built our anatomic model of human right and left atrium by folding the described two-dimensional cardiac tissue model into a 3D structure of one layer of cells. The human atria are modeled by a simplified structure composed of two connected ellipsoids with a diameter of 4.5 cm and an overall length of 8 cm. Non conductive obstacles are placed to simulate the veins and the valves, as shown in Figure 1. The total area of the structure is approximately 100 cm², with the orifices accounting for 20% of this area. The surface of the mitral (MV) and tricuspid (TV) valves is 4 cm², 2 cm² for the superior vena cava (SVC), 2.5 cm² for the inferior vena cava (IVC) and 1 cm² for each pulmonary vein (PV). These values are comparable with those found in human atria.

In spite of its crude approximation, our atrial geometry presents all major anatomical obstacles important for the simulation of atrial arrhythmias. Furthermore, it has a realistic size with a reasonable spatial discretization. Further work will enable us to develop a more accurate anatomical structure and geometry of the atria including pectinate muscles and a septal structure.

(a)	(b)

Figure 1. Geometry of the proposed human atria model with the holes represented in black, (a) anterior view with Superior Vena Cava (SVC) and Sino-Atrial Node (SA Node) visible in the right atrium, (b) posterior view with Inferior Vena Cava (IVC) and Tricuspid Valve (TV) visible in the right atrium, Pulmonary Veins (PV) and Mitral Valve (MV) visible in the left atrium.

Electrical activation is initiated from selected regions by injecting an intracellular current (I_stim in Equation (1)) to simulate periodic depolarization of the sino-atrial (SA) node, or premature stimuli. Furthermore, different action potential modifications have been tested by modulation of the ionic currents, simulating membrane defects or electrical disturbances.

The model has been implemented in C++ with double precision and runs on multiple platforms, including Windows PC's, cluster of Linux PC's, and a SGI Cray Origin2000 with 38 64-bits processors. The computation of 1 ms of real time typically requires 1 min of CPU with a single Pentium II processor.

Results

Sinus Rhythm Propagation

Normal propagation was initiated from the SA node region. From there the activation wavefront propagates over the right to the left atrium. In the model of healthy tissue, the propagation of a normal sinus beat results in a total activation time for the atria of about 90 ms. The propagation speed is 90 cm/s and the action potential duration is 289 ms at a cycle length of 1000 ms. This propagation velocity is in the range of values reported for physiological tissues [7]. A simulation of 5 sinus rhythm beats at a rate of 60 bpm is presented real time in Video 1.

Video 1: Sinus rhythm propagation at 60 bpm in a healthy atrial tissue (propagation velocity c = 90 cm/s). The video is displayed real time and the total depolarization time for both atria is 90 ms. (The VLC Media Player is recommended for viewing the video.)

Atrial Arrhythmias

Several experiments of arrhythmia initiation have been performed, using a programmed stimulation protocol as in electrophysiological studies. S₁ was the basic drive impulse initiated from the SA node at a cycle length of 1000 ms. Ectopic beats S₂ and S₃ (area 1.5 mm²,  intensity equal to the double of the fully repolarized threshold) have been introduced with various coupling intervals and at several locations.  In the healthy tissue, defined with a normal anatomy and size, all initiation attempts have failed since no sustained arrhythmia could be seen. Sustained reentry could not exist in the healthy tissue model due to wavelength constraints: an action potential with a duration of 289 ms and a propagation velocity of 90 cm/s leads to a wavelength of 27 cm, which is much larger than the atrial dimensions. Therefore, at basic conditions, this model seems electrically stable. Arrhythmias could be initiated when the wavelength was decreased to a value lower than the atria dimensions. First, the propagation velocity was reduced from 90 cm/s to 30 cm/s, either globally in the whole tissue, or locally (to mimic slow conduction in the isthmus between IVC and TV for example) by increasing the resistivity from 80 Ohm·cm to 800 Ohm·cm. As already mentioned, this decreased propagation velocity can also be interpreted as a normal propagation velocity in dilated atria. Furthermore, a decrease of action potential duration by an inhibition of the I_s current in the Beeler-Reuter model facilitates reentry because of an action potential shortening.

Even in a tissue with decreased wavelength, most programmed stimulation protocols have led to non sustained arrhythmias. Only critically timed and located premature beats could initiate atrial flutter or atrial fibrillation. Atrial flutter has been induced by applying one single ectopic beat S₂ between IVC and TV carefully timed after a sinus rhythm propagation. This leads to a single periodic reentrant wavefront anatomically anchored around IVC and TV. These results are similar to those reported in [2]. Atypical flutter can also be induced with premature beats in the PV region by using the same technique. The average rate of flutter in our model is 250 bpm and the periodic structure is clearly visible.
The initiation of atrial fibrillation in the model is more complex than flutter initiation. Several (two or more) carefully located and timed ectopic beats following a normal sinus rhythm depolarization are required to reach a fibrillatory state. Furthermore fibrillation could only be initiated in a tissue with global slow conduction. Reduced action potential duration facilitates arrhythmia perpetuation but is not required. Up to 30% inhibition of I_s current has been introduced in various simulations.  Although most attempts resulted in unstained wavefront perpetuation, by carefully timing the two ectopic beats, the model was able to initiate atrial fibrillation and sustain it for more than 40 seconds (limited only by the simulation time). Figure 2 shows an example of atrial fibrillation initiated in a tissue with 10% inhibition of the I_s current, a propagation velocity of c = 30 cm/s, and a stimulation protocol  involving two ectopic beats S₂ and S₃ located in the high right atrium with an S₁-S₂ interval of 350 ms and an S₂-S₃ interval of 225 ms. After an initiation period with multiple wave breaks and wave front interactions due to slow recovery fronts [8], a fibrillatory state appears in the model.
 

(a) 1650 ms

(b) 1700 ms

(c) 1750 ms

(d) 1800 ms

Figure 2: Snapshots of atrial fibrillation (from 1650 ms to 1800 ms after basic sinus rhythm initiation). AF was initiated by applying two ectopic beats S₂ and S₃ with a coupling interval of 350 ms and 225 ms respectively. The propagation velocity is c = 30 cm/s and the action potential duration is 260 ms.

In this example fibrillation last for about seven seconds and then is converted spontaneously to sinus rhythm. Compared to the periodic pattern of atrial flutter, atrial fibrillation is observed as multiple reentering wavelets traveling randomly and interacting with each others. Fibrillation involves both anatomically anchored and free spiral tips, thus being a complex mixture of anatomic and functional reentries, clearly visible in the PV region of Figure 2. The average rate of atrial fibrillation  in our model is 420 bpm.
Similarly to what we can observe in clinical experiments, most of the simulated atrial fibrillations convert spontaneously to atrial flutter or sinus rhythm. An example of initiation of atrial fibrillation followed by a spontaneous conversion to atypical flutter is presented in Video 2, which shows 10 seconds of real time propagation. Initiation involves a 3 ectopic beats protocol following a sinus rhythm depolarization, very similar to the one used for Figure 2, in a tissue also having similar properties (c = 30 cm/s, I_s inhibition of 10%). After 2 seconds, a fibrillatory state is reached and last for about 3 seconds. Finally, a spontaneous reorganization occurs and stabilizes the remaining wavefront anchored to IVC and TV, thus leading to a stable and periodic atypical flutter. The difference between the flutter and fibrillation patterns is clearly visible in Video 2.

It is interesting to observe the number of wavelets during the initiation, perpetuation and conversion of the atrial fibrillation described in Figure 2 (see Figure 3(a)). The stimulation protocol induces one main wavelet having one of its tip anchored around SVC, while its other free tip starts meandering. After a 1500 ms initiation phase, multiple wavefronts are produced from the free spiral tip and spread in the whole tissue, finally leading after 2500 ms to atrial fibrillation with four to seven wavefronts. Episodes with few wavefronts sometimes happen during fibrillation: two of such episodes can be seen in Figure 3(a) around 5000 ms, where only one single wavefront remains for few milliseconds, but because of inhomogeneities in tissue repolarization due to the previous fibrillatory state, multiple wave breaks rapidly lead back to fibrillation instead of sinus rhythm. The simulations also confirm the multiple wavelets hypothesis by Moe [9], which states that fibrillation is maintained by the presence of at least four independent wavelets.

 

Video 2: Initiation of atrial fibrillation and spontaneous conversion to atypical flutter. The video is displayed real time and the total duration of the simulation is 10 seconds. The propagation velocity is c = 30 cm/s and the action potential duration is 260 ms.
(The VLC Media Player is recommended for viewing the video.)

In addition to the number of wavelets, the amount of excited tissue, defined here as a transmembrane potential V_m > -70 mV, has also been computed and is shown in Figure 3(b). The initiation phase covering the first 2500 ms of simulation is characterized by important variations in the percentage of excited tissue . This reflects the global tissue behavior driven by a unique reentrant wavefront. The periodicity of this wavefront anchored around SVC is responsible for the periodic variation of the amount of excited tissue. Then a significant drop of these variations can be observed during the first fibrillatory episode characterized by a high number of wavelets (from 2500 ms to 3500 ms). This can be explained by the fact that the periodic behavior is replaced by multiple simultaneous depolarization and repolarization processes each related to an independent wavelet, globally leading to significantly reduced variations in the amount of excited tissue (stabilized between 50% to 60%).  As the number of wavefront decreases, an augmentation of the variations is observed. It can be noted that between 5500 ms and 6500 ms, there is a gradual reorganization of the fibrillation process that will finally lead to sinus rhythm.
 
 

(a)

(b)

Figure 3: Initiation of atrial fibrillation and spontaneous conversion to sinus rhythm: (a) development of the number of independent wavelets, (b) development of the percentage of excited tissue (V_m > -70 mV) . The stimulation protocol is described in Figure 2.

Conclusions

We have built a virtual atrium using validated membrane channel physiology and conduction features as described in physiological experiments. As the size and conduction velocity are within realistic values, the compromise made so far in this model is that the atrial walls are considered to be a single cell layer only. With few exceptions, this is acceptable for the study of atrial arrhythmias, while it would not be the case on the ventricular level. Another limitation of our model is that the whole atrial tissue (except the anatomical holes) has uniform properties (action potential duration and propagation velocity). Future work will consist in the development of a more accurate structure and geometry and the introduction of anisotropy and heterogeneity based on data from research in basic electrophysiology and clinical experiments on atrial arrhythmias.
Simulations of up to 40 seconds of real time have been performed, while in the literature, results in realistic models with a high number of spatial nodes are limited to about one second [1]. Functional and anatomic reentries were initiated only for appropriate basic conditions, such as slowed conduction or reduced action potential duration. During our virtual electrophysiological experiments, we have been able to initiate atrial flutter-like reentry showing a single macro reentrant circuit with a periodic pattern, as well as atrial fibrillation. The number of independent wavelets is comparable to those observed in humans [7,10].

These simulations allow us to observe in great detail how atrial arrhythmias originate from the prematurely stimulated tissue and how they are perpetuated or spontaneously converted to atrial flutter or sinus rhythm. One important observation is that in our model atrial arrhythmias are a combination of functional and anatomic reentries and therefore the geometry plays an important role. Indeed, the effect of anatomy of the heart on the mechanisms of arrhythmia is clearly visible in our experiments. In the mechanism of fibrillation, we have observed that veins and valves have two contradictory effects: on one hand they act as anchors for the waves [11], having therefore a stabilization effect and determining the reentrant pathways as well as the rate of arrhythmia. On the other hand they also induce multiple wave breaks mainly in the PV region, thus destabilizing propagation as shown in the snapshot sequence of Figure 2.  Pulmonary veins were known to be arrhythmogenic [12], but from our experiments they also seem to play a major role in fibrillation perpetuation because of their wave breaking and anchoring capacities. Another critical region for the initiation of atrial arrhythmias is the isthmus between IVC and TV.The transition from fibrillation to a periodic and more stable reentrant wave is facilitated due to the anchoring phenomenon. It seems that major anatomical obstacles like IVC, SVC and both valves play a more important role in this anchoring process than smaller holes like PV.

In conclusion, our model of human atria, even if it is based on simplified assumptions, can reproduce atrial arrhythmias showing a similar behavior than reported in clinical experiments. The advantage of the computer model is that it can show details difficult to study in nature and the experiments are reproducible. It will be used as a tool to study the potential and theoretical impact of the following therapeutic strategies: anti-tachy pacing, defibrillation, ablation and even pharmacological interventions.

References