A Simple Model for Identifying Critical Structures in Atrial Fibrillation

Atrial fibrillation (AF) is the most common abnormal heart rhythm and the single biggest cause of stroke. Ablation, destroying regions of the atria, is applied largely empirically and can be curative but with a disappointing clinical success rate. We design a simple model of activation wavefront propagation on a structure mimicking the branching network architecture of heart muscle and show how AF emerges spontaneously as age-related parameters change. We identify regions responsible for the initiation and maintenance of AF, the ablation of which terminates AF. The simplicity of the model allows us to calculate analytically the risk of arrhythmia. This analytical result allows us to locate the transition in parameter space and highlights that the transition from regular to fibrillatory behaviour is a finite-size effect present in systems of any size. These clinically testable predictions might inform ablation therapies and arrhythmic risk assessment.

AF. But why rotors should form and their relationship to the underlying structure and its arrhythmogenic characteristics are unknown. Heterogeneity in structural and electrophysiological tissue properties (e.g., fibrosis [6], orifices [7], action potential duration [8], restitution [9], tissue thickness [10]) have all been implicated in the perpetuation of AF.
Ablating atrial tissue extensively and empirically via catheter electrodes may cure AF in some cases [7,11], but an inability to identify the specific regions critical to the persistence of AF has resulted in the failure to improve on disappointing clinical outcomes. Thus, enhancing mechanistic understanding, identification of critical causative regions, and prediction of arrhythmic risk for preventative strategies requires a novel approach.
Heart muscle cells are discrete and so it is natural to use cellular automata to model activation wave fronts as was done in some of the earliest mathematical models of AF [12,13]. These studies investigated the effect of heterogeneity in functional parameters (e.g., refractory period) and macroscopic obstacles on the dynamics of wave fronts. In agreement with physiological experiments, the studies showed that an increase in heterogeneity in the refractory period resulted in greater wave front complexity. Bub, Shrier, and Glass devised cellular automata models of cardiac conduction and experimentally compared the behavior of the model to cardiac cell monolayers [14,15]. Bub et al. have also investigated the relationship between structural heterogeneity and the stability of activation wave fronts in isotropically coupled tissue for comparison with cardiac monolayers [16].
In contrast, we have investigated the effect of anisotropic structural heterogeneity in the underlying lattice of cells to mimic the lateral uncoupling of myocardial strands by fibrosis in cardiac muscle [3]. The level of structural heterogeneity is controlled by ν, a model parameter that determines the fraction of cells with lateral cell-to-cell couplings. Qualitatively, the model reproduces many of the known characteristics of real AF. For example, we find that with low levels of structural heterogeneity wave propagation remains planar. As the level of structural heterogeneity is increased beyond a threshold value, we observe spontaneous localized disruption of propagation and progressive degeneration to fibrillation, either self-terminating or sustained. Hence, the simple phenomenological model reproduces the observation that age related changes, that is, an increase in structural anisotropy [3,17], might spontaneously induce and sustain AF. We stress that AF behavior is initiated spontaneously in this model as a result of the pertinent underlying structure of the tissue as opposed to a combination of rapidly paced tissue and functional heterogeneity or macroscopic obstacles.
The simplicity of the model provides mechanistic insight into the initiation and maintenance of AF. We can characterize and identify the composition of local regions critical to the initiation and maintenance of AF and show that ablation of such critical regions may terminate AF. Furthermore, we predict the existence of a critical level of heterogeneity beyond which AF occurs. Because of the simplicity of the model, we can analytically calculate the risk of AF as a function of model parameters. This analytic expression reveals that the threshold value of structural anisotropy decreases with an increasing refractory period. Bringing the threshold ν ★ below ν by increasing the refractory period prevents microreentry circuits from persisting or forming because the wave front in the circuit collides with its unexcitable tail.
Our model of AF is complementary to biophysical models where subcellular processes, such as ion channel activity, are taken into account [18,19]. Biophysical models of cardiac tissue commonly consist of a system of partial differential equations (PDEs), known as the bidomain model, coupled nonlinearly to a system of ordinary differential equations modeling the cell membrane dynamics. With simplification such biophysical models are amenable to analytic calculations [20]. E.g., Keener [21,22] and Biktashev [23] have derived analytic calculations for spiral wave dynamics on a reduced form of these equations. PDE based models assume a continuous medium while our main investigation is on the effects of anisotropic structural heterogeneity for which the discreteness of the tissue may be important [24][25][26].
A healthy atrial area is L′ × L′ ≈ 20 cm 2 [27,28]. Atrial muscle tissue consists of nearly cylindrical myocytes of length Δx′ ≈ 100 μm and diameter Δy′ ≈ 20 μm packed in an irregular brick-wall-like pattern [29,30]. Cells are coupled primarily from end to end along the longitudinal axis of the cell with little side-to-side coupling, thus forming a branching cablelike structure [29,31]. The structure and strength of these connections can change; e.g., fibrosis can uncouple cells, primarily in a side-to-side fashion [29]. To represent this branching network of interacting cells, each cell is coupled to its neighbors longitudinally and to its neighbors transversally with probability ν. Structurally, the substrate can be pictured as a set of aligned longitudinal cables of single cell thickness which are coupled to neighboring cables with frequency ν. The structure is a 2D sheet. A 2D sheet is a good first approximation because the atrial muscle wall is relatively thin, 2.5 mm [31]. However, one might extend the model to 3D by coupling adjacent 2D sheets of tissue. To mimic the basic topological features of an atrium we apply periodic boundary conditions vertically and open boundary conditions horizontally to the 2D model. This provides a cylindrical topology which is a convenient substrate to study the evolution of wave fronts as our focus is on investigating the effects of increasing anisotropic structural heterogeneity, rather than the detailed effects of atrial anatomy. Once the substrate is generated, it remains fixed throughout the numerical simulation. Thus coupling related changes in the structure are incorporated into the substrate by the single variable ν.
The coordinated contraction of the heart is a result of a regular propagating wave front of activation originating at the heart′s natural pacemaker. The ionic currents which determine the activation of a cell are not explicitly considered in the model. In its simplest representation, a cell may be in one of three states: resting (repolarized), excited (depolarizing), or refractory, see Fig. 1. An excited cell induces neighboring resting cells to become excited and the wave front is a coherent propagation of this excitation throughout the tissue. The time taken for a cell to depolarize, Δt′ ≈ 0.6 ms, is much shorter than the refractory period, τ′ ≈ 150 ms. Cells along the left boundary are pacemakers and self-excite at a fixed period T′ = 300-1000 ms to mimic normal cardiac rhythm.
Cellular electrical dysfunction of any cause introduces a degree of noise into the excitation process. A fraction δ of randomly chosen cells in the substrate are assigned to be dysfunctional. Such cells have a finite probability ∊ of not exciting in response to an excited neighbor and hence can block the transmission of excitation.
Translating real tissue values into the model yields L = L′/Δx′ = 1000. We coarse grain the model by taking Δx′ → bΔx′, where b is the number of cells within a unit of space in the model. As the conduction velocity θ′ x = Δx′/Δt′ ≈0.2 ms −1 is fixed we also get Δt′ → bΔt′. We take b = 5 which gives L = 200 and τ=τ′/(bΔt′) = 50. The evolution of the system varies depending on the fractions of transverse connections ν and dysfunctional cells δ which are the only remaining parameters. Coarse graining the transversal connections gives ν → 1 − (1 − ν) b as this is the probability that a group of b cells have at least one vertical connection between them. The parameters of the model are system size L × L, refractory period τ, the fraction of dysfunctional cells δ, the probability of dysfunction ∊ = 0.05, the period of pacing T = T′/(bΔt′) = 220 using T′ = 660 ms, and the fraction of transverse connections ν.
When ν = 1, mimicking neonatal tissue in humans and other mammals, we observe regular wave fronts for all δ, consistent with the rarity of AF observed in children. We now investigate the effect of structural heterogeneity (ν < 1). For simplicity, we consider a fixed fraction of dysfunctional cells (δ = 0.05) leaving the fraction of transverse connections ν as the only control parameter. Hence, the parameter space ν may be divided into a rotor forming range and a stable range in which no rotors form, see Fig. 3. The frequency of formation and the temporal duration of rotors increases sharply from zero as the fraction of transverse connections ν decreases below a threshold value ν ★ ≈ 0.14. Only when ν is sufficiently small will the relevant path length of the reentrant circuits exceed the refractory period, which is a necessary condition to sustain a rotor, see Fig. 2(e). Decreasing ν results in rotors forming and disassociating more rapidly, resulting in a fibrillatory state, see Fig. 3. This mimics the natural onset of AF and its transition to paroxysmal and persistent atrial fibrillation. Moreover, this is consistent with the increase of AF incidence with advancing age. The critical regions forming a simple reentrant circuit as displayed in Fig. 2(e) are primarily responsible for the fibrillation observed in the model.
In our model, we can analytically calculate the risk of developing atrial fibrillation. We identify the risk of developing atrial fibrillation with the probability that the L × L system has at least one region with a reentrant circuit as displayed in red in Fig. 2(e). The probability that a cell has at least one transverse coupling is p ν = 1 − (1 − ν) 2 . Given an arbitrary dysfunctional cell i, let ℓ i denote the distance from that cell to its first neighbor to the right that has at least one transverse coupling. Hence, the probability that ℓ i < τ/2 (assuming τ is even) (1) The average number of dysfunctional cells is δL 2 and P risk is the complement of all dysfunctional cells having l i < τ/2: (2) The first expression is in dimensionless units, while the second expression is in units with dimensions where Δt′ = Δx′/θ′ x and θ′ x is the conduction velocity in the longitudinal direction.
The analytically calculated risk's fit to the data, see solid red line in Fig. 3(a), confirms that the critical regions shown in Fig. 2(e) drive fibrillation in the model. By rendering the region containing this structure unexcitable, mimicking ablation, we find that the rotor terminates and planar wave fronts reappear, representing a return to regular cardiac rhythm, see Figs.

4(a)-4(f).
The model displays a sharp transition from regular to fibrillatory behavior as a function of the control parameter ν. The existence of the transition is technically a finite-size effect because as the system size increases, the planar wave front range (C) of the phase diagram decreases and will eventually disappear for L → ∞. However, the threshold value ν ★ , defined by the point of steepest slope of P risk , increases extremely slowly with system size as ν* ≈ 1 − (δL 2 ) −(1/τ) ≈ 1 − (δL′ 2 /b 2 Δx′Δy′) −bΔx′/θx′τ′ for τ ⪢ 1 and δL 2 τ ⪢ 1. We note that ν ★ decreases with increasing refractory period. For L = 200, we predict that the transition would occur at ν ★ ≈ 0.14 using δ = 0.05 and τ = 50. The calculation of P risk can be extended to multiple cell layers (see Supplemental Material [32]) and the abrupt transition remains.
Antiarrhythmic drugs which prolong the refractory period increase the wavelength of reentry θ′ × τ′, thereby reducing the number of possible independent electrical wave fronts which can propagate in the atrium. This reduces the complexity of the fibrillation. However, in addition to this, our model also predicts that prolonging the refractory period reduces the structure-related risk of inducing fibrillation by rendering once arrhythmic regions, ν < ν ★ (τ), nonarrhythmic, by lowering ν ★ . This implies that functional changes can render structurally arrhythmic regions nonarrhythmic without actually modifying the structure.
The incidence of AF increases with age and is associated with the accumulation of fibrosis.
Our objective was to gain insight into the mechanisms responsible for the initiation and maintenance of AF in this context using a phenomenological model. This deliberately simple model incorporates pertinent basic properties of tissue architecture and electrical dysfunction, excluding electrophysiological properties normally thought to be important in fibrillation such as the dispersion of refractoriness and restitution. The most interesting aspects of this study are the following. (1) The mere simplicity of the model. (2) The integration of simplified dynamics on a relevant discrete medium mimicking the anisotropic branching architecture of heart muscle. (3) A single model that reproduces multiple wellknown observed features of AF. (4) The spontaneous emergence of AF (without external stimuli) when the transversal coupling is reduced, mimicking fibrosis, beyond a threshold value. (5) Identification of local critical regions initiating and maintaining AF. (6) Analytical tractable risk of AF. (7) The existence of an analytically tractable threshold value of heterogeneity below which AF emerges, providing insight into how to control the risk. (8) A targeted ablation procedure is possible for those patients at the early onset of AF driven by fibrosis due to the small number of critical regions.
It would be of interest to investigate a nonuniform distribution of vertical coupling (e.g., localized fibrosis) or refractory period [33]. The present model identifies specific targets for ablation based on the underlying structure of the myocardium. Patterned cultures of cardiac myocyte monolayers or cardiac tissue slices could test the primary predictions of the model as stated above. If the biological experiments agree with the model predictions then such targets within a biological system might be characterized electrophysiologically or architecturally with the aim of translating such insight into clinical applications and thereby, it is hoped, improving the success rate of treating AF with ablation.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.   When the wave front in the top cable reaches the rightmost vertical connection, it induces activity in the bottom cable that will also propagate retrograde. If the relevant path length (red line) of the circuit is greater than the refractory period, the cell to the left of the dysfunctional cell will be in a resting state when the retrograde wave front arrives and a reentrant circuit of activation forms. The reentrant circuit terminates when the dysfunctional cell blocks the retrograde wave front. This is the simplest example of a critical region which will induce fibrillationlike behavior in the model.  The ablation terminates the reentrant circuit and the system returns to a planar wave front phase. See Fig. 3 in the Supplemental Material [32] for the full image and animation.