Reentry and Ectopic Pacemakers Emerge in a Three-Dimensional Model for a Slab of Cardiac Tissue with Diffuse Microfibrosis near the Percolation Threshold

Abstract

Arrhythmias in cardiac tissue are generally associated with irregular electrical wave propagation in the heart. Cardiac tissue is formed by a discrete cell network, which is often heterogeneous. Recently, it was shown in simulations of two-dimensional (2D) discrete models of cardiac tissue that a wave crossing a fibrotic, heterogeneous region may produce reentry and transient or persistent ectopic activity provided the fraction of conducting connections is just above the percolation threshold. Here, we investigate the occurrence of these phenomena in three-dimensions by simulations of a discrete model representing a thin slab of cardiac tissue. This is motivated (i) by the necessity to study the relevance and properties of the percolation-related mechanism for the emergence of microreentries in three dimensions and (ii) by the fact that atrial tissue is quite thin in comparison with ventricular tissue. Here, we simplify the model by neglecting details of tissue anatomy, e. g. geometries of atria or ventricles and the anisotropy in the conductivity. Hence, our modeling study is confined to the investigation of the effect of the tissue thickness as well as to the comparison of the dynamics of electrical excitation in a 2D layer with the one in a 3D slab. Our results indicate a strong and non-trivial effect of the thickness even for thin tissue slabs on the probability of microreentries and ectopic beat generation. The strong correlation of the occurrence of microreentry with the percolation threshold reported earlier in 2D layers persists in 3D slabs. Finally, a qualitative agreement of 3D simulated electrograms in the fibrotic region with the experimentally observed complex fractional atrial electrograms (CFAE) as well as strong difference between simulated electrograms in 2D and 3D were found for the cases where reentry and ectopic activity were triggered by the micro-fibrotic region.

About the Authors:

Sergio Alonso

* E-mail: [email protected]

Affiliations Department of Physics, Universitat Politècnica de Catalunya, Barcelona, Spain, Physikalisch-Technische Bundesanstalt, Berlin and Braunschweig, Germany

ORCID http://orcid.org/0000-0002-3989-8757

Rodrigo Weber dos Santos

Affiliation: Graduate Program in Computational Modeling, Federal University of Juiz de Fora, Juiz de Fora, Brazil

Markus Bär

Affiliation: Physikalisch-Technische Bundesanstalt, Berlin and Braunschweig, Germany

Introduction

The synchronous contraction of the heart is a consequence of the highly organized propagation of the membrane action potential through the whole cardiac tissue. Several cardiac arrhythmias correspond to conditions where this synchronization fails and cardiac contraction is compromised. One of the most relevant arrhythmias is fibrillation which corresponds to a state where the contraction of many cardiac myocytes are strongly desynchronized in the atria, i.e. atria fibrillation, or, with more dangerous consequences, in the ventricles, i.e. ventricular fibrillation. Reentry is a circulating action potential wave, and believed to be a necessary but not a sufficient condition for fibrillation. Fibrillation is often attributed to a succession of multiple breakups and pairwise annihilations of spiral waves.

Spiral generation and spiral breakup [1–4] have been extensively studied by computational modeling, for reviews see [5, 6]. The generation of a spiral wave may result from the impact of an ectopic beat: an additional, premature stimulus (S2) which generates a second wave that interacts with the sinus wave (S1) giving rise to a spiral wave [5]. A second common mechanism for spiral wave generation is cardiac alternans, the alternation of short and large pulses in response to a high frequency pacing [6, 7].

Even though atrial fibrillation (AF) is the most common type of sustained cardiac arrhythmia [8], the mechanisms that trigger it are still poorly understood. One hypothesis supported by clinical studies [9] and animal experiments [10] is that AF can be initiated by one or more regions that generate ectopic beats. For example, six consecutive ectopic beats separated by 130 ms may induce atrial fibrillation [11]. Structural heterogeneities may also produce persistent atrial fibrillation [12].

In this work, we study how an heterogeneous region with the presence of both full-connected and non-conducting (fibrosis) areas can generate microreentry in 3D computational models of cardiac tissue. This type of fibrosis is known as microfibrosis due to the microscopic scale of the heterogeneities [13–15].

In 2007, Spach et al. [16] could reproduce their original experiments [17] using a 2D in-silico model that also accounted for fibrosis. A single cycle of microreentry was generated after using a standard S1–S2 protocol. They concluded that fibrosis increases the vulnerability to fast pacing. Similar results using multiple stimulus protocols were obtained by different groups [16, 18–22] and confirmed that fibrosis or other micro-structural heterogeneities may facilitate the initiation of reentries and help to sustain reentry or fibrillation. The use of the words facilitate and help appears naturally, since it is well known that S1–S2 or fast-pacing protocols can per se generate reentries. Therefore, there exists a complicated interaction among the different waves (S1–S2) and the heterogeneity of the tissue (fibrosis).

However, it was only recently shown [23, 24] that we can generate microreentries and ectopic pacemakers in 2D computer simulations of cardiac electrophysiology without the use of a second premature stimulus (S2) [16] or fast pacing [18]. If only a single stimulus is used, the usual interactions between subsequent waves are not present; and therefore, one can relate the genesis of the microrentries inside the fibrosis to the specific heterogeneity of the tissue, i.e. to its micro-structure rather than to functional features (like wave interactions).

The fibrosis forms a maze where waves fractionate and follow zig-zag pathways. Therefore, the topology of the maze is crucial for the emergence of a reentry inside the fibrotic region [23, 24]. This wave activity may continuously or transiently generate ectopic beats when it crosses the border of the microfibrosis region [18] into the healthy tissue areas. Such zig-zag dynamics was previously observed in experimental studies [25].

A straightforward approach to the modeling of microfibrosis is the use of simple 2D discrete and isotropic models of cardiac tissue, wherein a prescribed percentage of connections between neighboring myocytes is randomly removed [23, 26]. The probability of microreentry depends on the degree of fibrosis specified by the fraction ϕ of non-conducting links between neighboring cells in the tissue. By analyzing the probability distribution of tissues with reentry as a function of ϕ, it was found [23] that the maximum of this distribution was close to the percolation threshold ϕ_c [27, 28] of the conducting links in the associated discrete network. The value of ϕ_c provides a reference parameter that depends only on the topology, and, therefore, on the spatial dimension of the lattice used in the simulation. A similar relation was verified later using two-dimensional microscopic anisotropic models of the tissue together with more detailed physiological descriptions [24]. These results indicate that the degree of fibrosis as well as the topology of connections between myocytes are important factors when evaluating the pro-arrhythmic behavior of a particular region of the heart.

In recent years, related studies were carried out by Christensen et al. [29] for an anisotropic tissue model, by Kazbanov et al. [30] for tissues with spatially varying degrees of microfibrosis and by Vigmond et al. [19] for the modeling of complex fractionated atrial electrograms (CFAE) and reentry in a fibrosis model based on imaging data.

Propagation near fibrotic area has been associated with CFAE, in both animal [31] and clinical studies [32]. These findings have remarkably changed the management and prognostics of patients with atrial fibrillation via a new procedure of catheter ablation of atrial fibrillation that is guided by CFAE mapping [33]. Recent clinical studies of guided CFAE ablation revisited the benefits and risks of this procedure and clearly suggest that further investigation on this topic is required [34–36]. Different computer simulations have also already given valuable insights into how the microstructure of cardiac tissue, specifically microfibrosis, and CFAE are related [13, 15].

Here, we extend our studies with 2D simulations presented before in [23, 24] to 3D simulations of thin cardiac tissue represented by a discrete model composed of up to twenty, two-dimensional layers arranged in a cubic lattice with an embedded region of microfibrosis that will be shown: i) to reproduce microreentry inside the microfibrosis tissue, ii) to form an ectopic pacemaker and iii) to reproduce qualitatively the signature of CFAE.

In addition, the simulations permit to quantify the influence in the activity of the ectopic pacemaker of important physiological parameters such as electrical remodeling, the fraction of fibrosis (i. e. the fraction of broken links), the size of the micro-fibrotic region and the tissue thickness.

Many simulations of the atria are based on the monolayer hypothesis [11, 37], however, tissue thickness of the atria is known to vary between 0.5 and 3 mm [38]. In addition, some recent studies have suggested the possibility of transmural reentries in the atria [39], a phenomenon that could only be captured with a 3D model. Here, we investigate how the number of layers in the third dimension influences the occurrence of microreentries in a fibrotic region of the atria.

Electrophysiological properties such as APD and conduction velocity (CV) vary substantially within the heart. In addition, different pathologies, such as atrial fibrillation, ischemia or a substantial increase of myocyte to fibroblast coupling, may significantly decrease the values of both APD and CV [11, 40–43]. To assess how different electrophysiological features influences the probability to generate microreentries in a fibrotic region we test three levels of electrophysiological remodeling: normal remodeling (APD = 135ms, CV = 44 cm/s); intermediate remodeling (APD = 120ms, CV = 39 cm/s), and strong remodeling (APD = 76ms, CV = 34 cm/s).

Finally, the degree, type and pattern of fibrosis varies considerably for different kinds of diseases or during different stages of one particular pathology [44]. Therefore, we have investigated how the percentage of fibrosis, ϕ, and the size of a fibrotic region influence the probability to generate microreentries and create ectopic pacemakers in cardiac tissue.

Materials and Methods

Continuous model of cardiac tissue

The propagation of the action potential is described by the cable equation [45]. The equation relates the variation of the membrane potential (V) with the total ion current through the ion channels across the membrane (I), and the conduction of potential along the cell membrane. The cable equation reads:(1)where D = 1 cm²/s is the effective diffusion coefficient of the action potential [6]. The above equation is solved in its discrete form, as will be explained in the next sections. Moreover, homogeneous Neumann boundary conditions are used. We do not consider anisotropy or any detailed geometry and anatomy of the tissue to avoid extensive calculations. This amounts to a qualitative description of the tissue which nevertheless captures the main properties of the dynamics as will be demonstrated below.

The total current I is the sum of different types of ion currents I = ∑_i I_i. The number and form of these currents is determined by the specific ionic model employed [46]. Here, we use a simplified and generic cardiac model, which consists of only three ion currents [4, 47]:(2)where the currents do not exactly correspond to an explicit ion channel, but are controlled by parameters which can be fitted to particular experimental measurements or to more detailed models. The ionic currents read:(3)corresponding, to the fast inward (I_fi), slow outward (I_so) and slow inward (I_si) currents. These currents are controlled by the gating variables v and w. The evolution of such gate variables depends on the membrane potential:(4)where with p = θ(V − V_c) and q = θ(V − V_v), where θ(x) is the Heaviside step function. We have used the same parameter values as described before in [23]. One of the main advantages of using this simplified and generic cardiac cell model is the ability to test different scenarios in terms of electrophysiological remodeling via the modification of a single parameter of the model, τ_d.

We employ finite differences for the spatial discretization (Δx = 100 μm) and a second order Runge-Kutta method for the temporal integration of the equations (Δt = 0.0165 ms) for the numerical integration of Eqs (2–4).

Excitability and level of remodeling

Electrophysiological properties such as APD and CV vary substantially within the heart. In addition, under different pathologies, such as atria fibrillation, ischemia or a substantial increase of myocyte to fibroblast coupling, both APD and CV values can be significantly decreased. For instance, whereas APD values can be as short as 140 ms in healthy atria, it is reduced to 56 ms under chronic atrial fibrillation [11, 40, 41]. Likewise, CV values may be as slow as 44 cm/s within normal atrial, whereas under chronic atrial fibrillation CV can be even slower, with reported values at 37 cm/s [11, 40–43].

Numerical simulation of electrical remodeling has been used in many studies of atria fibrillation. Action potential duration and conduction velocities decreased under remodeling to: APD = 100 ms and CV = 66 cm/s in [48], APD = 144 ms and CV = 54 cm/s [38], and APD = 75 ms and CV = 61 cm/s [11]. AP remodeling can be due to completely different reasons, and the level of remodeling depends not only on the type of pathology but also on its stage.

In this work, we use a simplified approach to assess how different electrophysiological features influence the probability to generate microreentries in a fibrotic region. Despite the multitude of possibilities in terms of APD and CV values reported, we simulate only three levels of electrophysiological remodeling: normal remodeling (APD = 135ms, CV = 44 cm/s); intermediate remodeling (APD = 120ms, CV = 39 cm/s), and strong remodeling (APD = 76ms, CV = 34 cm/s).

These three cases of remodeling were simulated by using different values of the parameter τ_d of the Fenton-Karma model [47]: 0.25 ms, 0.30 ms and 0.35 ms. The action potential duration (APD) and the conduction velocity (CV) of our simplified cell model decrease with the parameter τ_d, and consequently, the wave length (WL) given by the product of these two quantities (WL = CV ⋅ APD) also decreases. Thus, waves are smaller and slower for larger values of τ_d, implying that the excitability decreases with this parameter. Therefore we make the qualitative relation between decreased excitability and electrical remodeling, i.e. higher values of τ_d corresponds to higher level of remodeling, as presented in Fig 1.

[Figure omitted. See PDF.]

Fig 1. Transmembrane voltage evolution for different remodeling levels.

Shape of the action potential duration for different values of τ_d employed in the three-dimensional simulations: τ_d = 0.25 ms (A), τ_d = 0.30 ms (B), and τ_d = 0.35 ms (C). Parameters used are: 3.33ms, = 19.6 ms, = 1000 ms, = 667 ms, = 11 ms, τ_o = 8.3 ms, τ_r = 50 ms, τ_si = 45 ms, k = 10, = 0.85, V_c = 0.13 and V_v = 0.055, for more details see [4]. These parameter values are kept constant for the rest of simulations.

https://doi.org/10.1371/journal.pone.0166972.g001

Discrete heterogeneous cardiac model

Randomly distributed inactive gap junctions, inert tissue cells among the cardiac muscle cells or the presence of unexcitable cells, e.g. fibroblast, in the tissue break the homogeneous propagation of action potential through the cardiac tissue. The resulting waves are expected to be corrugated and noisy. We model this disruption in the propagation by the introduction of random heterogeneities in the network of cells using the connectivity parameter η_ij:(5)where ξ = D/ℓ² with ℓ = 100 μm. If the connectivity parameter η_ij = 1 for all i and j we recover the homogeneous discrete network limit of Eq (2). The connection between two cells (i and j) is removed (η_ij = 0) with a certain probability ϕ. These abnormal connections are randomly scattered in the tissue. Our approach is a simple and qualitative model of reactive interstitial fibrosis [49].

Definition and calculation of the percolation threshold

The increase of the fraction of non-conducting links leads to isolated regions confined by removed links. Below a particular fraction of ϕ_c the network exhibits a percolation of conducting links [50]. Wave propagation is not possible for fractions above this percolation threshold. For ϕ just below the value of ϕ_c, the effective medium approach [51] that enables modelers to replace heterogeneous discrete models by homogeneous continuum models breaks down and a single excitation wave may break into pieces and produce reentries and persistent irregular activity [23].

The percolation threshold depends on several factors including the spatial dimension of the system. The value of the fraction ϕ that produces conduction block in a particular two-dimensional lattice may permit propagation in a different lattice with higher connectivity or in the corresponding 3D medium. Therefore, it is not clear that previous results in 2D can be extrapolated to 3D. Here, we want to evaluate this transition from 2D to 3D systems by using a discrete model of up to 20 layer thickness defined on a cubic lattice. The layers of this model correspond to the square lattices employed in early studies in 2D [23].

We systematically increase the number of layers (N) in a three-dimensional model starting from a monolayer, i.e. the two-dimensional system, assuming a layer thickness of 100 μm. Altogether, we vary from N = 1 (monolayer, 2D) to N = 20 layers (3D tissue with a thickness of 2 mm). For these parameters the percolation threshold can be easily calculated from one hundred randomly generated networks for different values of ϕ [23]. The change of the percolation threshold is plotted in Fig 2.

[Figure omitted. See PDF.]

Fig 2. Percolation threshold for different slab thickness.

Dependence of the percolation threshold in a slab of tissue for different thickness. For each thickness N we generate 100 realizations for each value of ϕ and check the probability of finding a pathway from one side of the system to the other. After a fitting process, we evaluate the fraction ϕ_c which correspond to a probability of 50%. The size of the system is 4 cm × 4 cm × N layers, where each layer corresponds to 0.01 cm.

https://doi.org/10.1371/journal.pone.0166972.g002

Simulation protocol

We consider that a small cylindrical region in the center of our simulation domain represents the fibrotic region and it is strongly heterogeneous. Inside this fibrotic region, a fraction ϕ of connections are non-conducting (η_ij = 0). Outside this fibrotic region all the connections between neighboring myocytes are normal and homogeneous (η_ij = 1). In addition, we model electrical remodeling varying the excitability of the fibrotic region, we set τ_d to , where can be 0.25 ms, 0.30 ms or 0.35 ms.

We evaluate the interaction between the wave propagating in an homogeneous region of the tissue with the fibrotic region with a fraction ϕ of non-conducting links. A single excitation wave of action potential is initially induced in a corner of the three-dimensional slab. The action potential propagates through the medium till it interacts with the fibrotic region. Depending on the values of ϕ and ϕ_c (percolation threshold for the slab used in simulation, see Fig 2) three behaviors are observed: For low fraction (ϕ ≪ ϕ_c) the wave propagates inside the fibrotic region and deforms, and although propagation inside is slower than outside, wave propagation is stable. For large fraction (ϕ ≫ ϕ_c) the wave does not enter into the fibrotic region and propagates only around this region. For intermediate fraction (ϕ ∼ ϕ_c) the wave propagates in the fibrotic region, however the propagation inside is highly irregular. We evaluate the different dynamics by simulating the activity during 1500 ms. The resulting numerical simulations are discussed in the next section.

Fibrotic border effects

We incorporate a fibrotic border (FB) that surrounds the region with microfibrosis. This border region models a smooth transition between non-fibrotic and fibrotic tissue in comparison with a stiff transition employed in our previous models [23, 24].

First, we compare the two different borders: 1) a localized fibrosis without FB, as presented in Fig 3(A). Here, a circular region of microfibrosis is immediately surrounded by homogeneous tissue; 2) fibrosis with a fibrotic border where the circular region of microfibrosis is immediately surrounded by a second zone where the probability of microfibrosis decreases with the radius to zero. Outside of these two regions there is no microfibrosis, see Fig 3(C).

[Figure omitted. See PDF.]

Fig 3. Effects of the fibrotic border in the fibrotic area.

Sketch of the fibrotic region without (a) and with (c) fibrotic border and probability of reentry without (b) and with (d) fibrotic border for N = 3 in a system 4 × 4 × 0.03 cm³ and R₀ = 1.40 cm for a total time of simulation of 1500 ms. Two types of dynamics are observed: reentry outside of the fibrotic region, Reentry (blue) or activity restricted to the interior of the fibrotic region, Activity (green).

https://doi.org/10.1371/journal.pone.0166972.g003

We run one hundred of simulations for different values of ϕ, see Fig 3, and calculate the probability of the different types of reentries observed in each set of simulations. We have also classified the observed reentries as activity or reentry [52], see sketches in Fig 3. The label activity describes the simulations that have resulted in sustained reentry but are confined to the region with microfibrosis. The label reentry characterizes the simulations where we observed ectopic beats leaving the microfibrosis and inducing sustained reentry all over the tissue, see a sketch in Fig 3. Simulations were monitored for 1500 ms and were classified as: (i) no reentry, (ii) non-sustained reentries, (iii) confined, sustained activity or (iv) global, sustained labeled reentry.

The total number of reentries in both cases are similar, compare Fig 3(B) and 3(D). However, the type of activity, confined activity vs. global reentry, varies systematically. The inclusion of a border region facilitated the transformation of confined activity to sustained, global reentry affecting also the homogeneous surrounding tissue. Therefore, we include the gradual fibrotic border in the next simulations.

Calculation of electrocardiograms

One way to indirectly assess the degree of fibrosis in cardiac tissue is the analysis of the waveform of electrograms measured by catheters tips located near cardiac tissue during a clinical electrophysiology study. Propagation near fibrosis is associated with fractionated electrograms [31, 32]. Two-dimensional computer simulations have already given valuable insights about the relation between microfibrosis and complex fractionated electrograms [13, 15, 19].

Here, we continue these studies by also simulating electrograms (EG) measured by electrodes near the 3D tissue surface (0.1 mm from the surface). The extracellular potential is computed using the large volume conductor approximation [15]. To this end, we integrate all the currents in the tissue with respect to the position of the electrode:(6)where is the distance from the electrode to the particular point of the tissue and A is a constant including the ratio of extracellular and the intracellular resistivities. For a list of different methods of calculations of EG see [53].

Results

The interaction of a wave with the fibrotic region was studied in two-dimensional models [23]. A remaining activity was observed inside the fibrotic area. Under certain conditions, this activity confined to the microfibrosis region may leave and re-excite the rest of the tissue, giving rise to the generation of ectopic beats. Here, we extend those results to 3D.

Reentry in a three-dimensional slab of tissue

In Fig 4 we show a typical example of a simulation leading to reentry due to localized fibrosis in 3D. The corresponding movie is provided as supplementary material (S1 Movie). First, the wave of action potential propagates around the circular region, see Fig 4(A) and 4(B). However, the action potential slowly enters into the fibrotic region, where it breaks down into pieces, see Fig 4(C)–4(E). One of these pieces arrives to the border where it induces the formation of a new reentrant wave, see Fig 4(F)–4(G). The action potential inside the fibrotic region induces new waves into the homogeneous tissue. Therefore, via the mechanism of microreentry inside a region with microfibrosis a localized ectopic pacemaker is generated.

[Figure omitted. See PDF.]

Fig 4. Reentry in a 3D slab of cardiac tissue.

Numerical simulation showing the reentry from the fibrotic region into the whole slab of tissue. Nine snapshots of the tissue slab every t = 66 ms, therefore for times t = 66 ms, t = 132 ms, t = 198 ms, t = 264 ms t = 330 ms, t = 396 ms, t = 462 ms, t = 528 ms, and t = 594 ms. Parameter values: ϕ = 0.72, , system size is 7 × 7 × 0.1 cm³ (N = 10 layers), radius of the fibrotic region R_o = 1.4 cm, and total time of simulation t = 1500 ms.

https://doi.org/10.1371/journal.pone.0166972.g004

In order to understand the process of reentry in more detail we plot in Fig 5 the maps of activation times. The time of the first activation of each cell is plotted in Fig 5(A) and 5(B). While the wave takes around 200 ms to arrive to the opposite corner of the system, it takes more than 400 ms to reach the center of the fibrotic area. Even more interesting is the time of the second activation of each cell, which only occurs in the case of a reentry. We observe in Fig 5(C) and 5(D) that the second activation occurs first near the border of the fibrotic region at 300 ms and it produces a propagating wave even before the center of the fibrotic region is excited for the first time at 450 ms (Fig 5(A) and 5(B)).

[Figure omitted. See PDF.]

Fig 5. First and second activation time maps in a 3D slab of cardiac tissue.

Map of the first activation of the cells corresponding to simulation shown in Fig 4 (A) and the enlargement of the map in the fibrotic area (B). Map of the second activation of the cells corresponding to simulation shown in Fig 4 (C) and the enlargement of the map in the fibrotic area (D). Black arrow shows the origin of the second activation times which produces the first reentry, and white arrow shows premature activation which dies without producing reentry.

https://doi.org/10.1371/journal.pone.0166972.g005

The probability of reentry depends on the thickness

Two-dimensional studies with the same model [23] predicted the appearance of reentry for a intermediate value of τ_d inside the fibrotic region (). They are summarized in the first row in Fig 6. For larger values of waves propagate too slow to enter into the fibrotic region (results not shown). While reentry probability arrives to 40% in two-dimensional systems with strong remodeling, only very small probability is found for and no reentry is found for .

[Figure omitted. See PDF.]

Fig 6. Dependence of the reentry probability on thickness and excitability.

Columns correspond to three different values of the parameter inside the fibrotic region corresponding to three different excitabilities. Outside this region τ_d = 0.25ms. Rows correspond to five different thickness of the three-dimensional slab of tissue (N). For each value of and N the dependence of the probability of reentry as function of ϕ is plotted. For each value of ϕ two types of reentries are plotted: full Reentry (blue) and confined Activity (green). Dashed lines correspond to the percolation threshold value. Parameter values: system size varies from 4 × 4 × 0.01 cm³ (N = 1 layer) to 4 × 4 × 0.2 cm³ (N = 20 layers), radius of the fibrotic region R_o = 1.4 cm, and total time of simulation t = 1500 ms.

https://doi.org/10.1371/journal.pone.0166972.g006

If we add two layers to the single layer employed for the two-dimensional simulations (N = 3), two important changes are observed, see second row in Fig 6. First, for strong remodeling, the whole distribution of reentries shifts, to higher values of ϕ, in agreement with the dependence of the percolation threshold on the thickness, see Fig 2. Second, we observe an increase of the reentry probability for intermediate and large values of .

Next, we consider a three-dimensional slab with five layers (N = 5), see third row in Fig 6. The results for five layers confirms the shift to larger values of ϕ due to the change of the percolation threshold, and the further increase of the probability of reentry for intermediate values of . In particular, for intermediate remodeling and ϕ = 0.7 the probability of reentry is as high as 90%. We highlight also that 5 layers correspond to 0.5 mm of thickness and that thickness is already observed in some regions of the normal atria.

The probability of reentry for , is much smaller when 10 layers are employed (N = 10). On the other hand, a sizable non-zero probability for reentry for , exists only in a narrow window of fractions of ϕ. And, more important, activity is observed already under normal remodeling with 6% of probability around the percolation threshold value. Finally, for thick slabs of tissue (N = 20), no reentry is found for strong remodeling in comparison with the high probability in two dimensional system. For intermediate remodeling the results are similar to those obtained with half of the layers. For normal remodeling there is a 8% of probability of activity again around the percolation threshold value. The increase of the thickness therefore shifts the optimal value of for reentry to occur to values that corresponds to intermediate remodeling, whereas in the pure 2D case, the probability was largest for strong remodeling (). Altogether, the strongest occurrence of reentry does not occur for normal or strong but for intermediate remodeling. This is presumably because for normal conditions the excitation waves are more robust against perturbation by the heterogeneous distribution of available conducting links. In parallel, the tendency for rotor formation weakens as remodeling gets more pronounced.

The probability of reentry depends on the size of the fibrotic region

To study the influence of the size of the heterogeneous region, we consider a set of parameter values with high probability of reentries: ms and N = 10, see Fig 6. In Fig 7 we show the results of these simulations, where we change the radius of the fibrotic region from 0.4 cm to 1.4 cm. The higher the size of the region is, the more likely the reentries are.

[Figure omitted. See PDF.]

Fig 7. Dependence of the reentry probability on the size of the fibrotic region.

For weak excitability () and thickness N = 10 layers we change the radius of the fibrotic region. For the radius here considered we observe reentry only for ϕ = 0.71, 0.72 and 0.73. The probabilities for each value of ϕ are plotted as function of the radius R₀. Parameter values:system size is 4 × 4 × 0.1 cm³ (N = 10 layers), and total time of simulation t = 1500 ms.

https://doi.org/10.1371/journal.pone.0166972.g007

With these results we can also estimate the minimum size for a fibrotic region in a 1 mm thick three-dimensional system to generate reentries: this is around 1 cm of diameter, which is in the same order than the fibrotic spatial patterns derived from patients [54]. This value is smaller than for a pure two-dimensional system where it is around 1.5 cm of diameter [52].

Simulation of fractionated electrograms in 3D

We consider again the 3D numerical simulation shown in Fig 4, where we choose 7 observation sites (electrodes), see Fig 8. The first panels of this figure correspond to four snapshots of the action potential propagation (view from the top of the surface of the 3D simulation), where the location of the electrodes are also shown.

[Figure omitted. See PDF.]

Fig 8. Electrograms at different locations in the tissue during reentry.

The signals of different electrodes are plotted in arbitrary units for a numerical simulation showing reentry. Four snapshot of the top view of the tissue slab for times t = 16.5 ms, t = 500 ms, t = 1000 ms, and t = 1500 ms. The location of the electrodes is shown in the first snapshot by white circles. The electrodes are located at 0.01 cm from the slab of tissue. Signals calculated in each electrode, following Eq (6), are plotted. Black and red curves correspond to electrodes inside and outside of the fibrotic region, respectively. Note that vertical axis scale is modified for panels C-E due to the weak signal inside the fibrotic area. Parameter values: ϕ = 0.72, , system size is 7 × 7 × 0.01 cm³ (N = 10 layers), R_o = 1.4 cm, and total time of simulation t = 1500 ms.

https://doi.org/10.1371/journal.pone.0166972.g008

During the simulation the EG measured by the different electrodes outside the fibrotic area are basically the same, see Fig 8(a), 8(b), 8(f) and 8(g), and equivalent to a EG measured in an homogeneous system. However, the EG measured by the three electrodes above the fibrotic region are different, see Fig 8. Inside this region the wave propagates slowly and it arrives to the middle with a large delay and pieces of waves arrive simultaneously from the border. The signals of the electrodes inside the fibrotic area are weaker, and look irregular and fractionated, see Fig 8(c)–8(e). Note that the amplitude of the EG inside the fibrotic region is weaker and the scale is one order of magnitude smaller. Actually the amplitude of the EG decreases linearly with the fraction ϕ due to the systematic loss of active tissue for the integral defined in Eq (6).

Two numerically calculated EG in 2D and 3D are compared in Fig 9. They correspond to conditions where microfibrosis induces reentries. The EG measured in the homogeneous part of the tissue is basically the same in 3D and in 2D, apart from an increase in the amplitude of the EG due to the larger number of cells in 3D, see Fig 9(A) and 9(B). However the EG in the fibrotic region looks different, see Fig 9(C) and 9(D). In the 2D case the EG is already fractionated but the waveform still resembles somewhat that one of the homogeneous EG. The 3D EG measured inside the microfibrosis is much more random and very different from the homogeneous EG, probably due to an averaging effect of uncorrelated signals from different layers of the simulation. This result shows that 3D effects on the measured EGs are not negligible and should be taken into account in future studies when comparing EG in more detailed and realistic physiological models with experimentally measured CFAEs.

[Figure omitted. See PDF.]

Fig 9. Comparison of Electrograms in 2D and 3D.

Electrograms calculated on two points corresponding to locations marked with (A) and (D) in Fig 8. Two-dimensional electrograms calculated under strong remodeling conditions outside (A) and inside (C) the fibrotic area (ϕ = 0.48) in comparison with three-dimensional electrograms calculated under intermediate remodeling conditions outside (B) and inside (D) the fibrotic area (ϕ = 0.72).

https://doi.org/10.1371/journal.pone.0166972.g009

Discussion

Percolation threshold and occurrence of microreentry in microfibrosis

Our previous simulations in 2D discrete and isotropic models of cardiac tissue with simple lattice topologies (square, hexagonal) exhibited microreentry inside a fibrotic region [23]. Equivalent results were obtained with two-dimensional models with sub-cellular resolution that took into account the anisotropy of the tissue and the gap junction distribution together with a much more complete physiological description [24]. In both studies, the probability of microreentry depended on the fraction of fibrosis, ϕ, and it was shown to be most prevalent just above the percolation threshold ϕ_c where conduction is lost. Since the percolation threshold depends only on the lattice topology, details of the electrophysiological should not matter much for the identification of the ϕ-values for which reentry is most likely.

In this paper, 3D simulations of homogeneous tissues with an embedded region of microfibrosis confirmed again that the probability to generate reentry from a microfibrotic region is higher just above the percolation threshold for conducting links in the discrete network. Reentry is likely to appear for percentages of non-conducting links between 50% and 70%. Note that myocyte connectivity decreases to a value of 65% in high-resolution 3D reconstruction of an infarct border zone [21]. We show a straightforward relation between the size of the fibrotic region and the probability of reentry: reentries are more likely to occur inside larger-sizes of fibrotic regions. The same relation was observed before in 2D simulations [23]. One possible explanation for this observation is that the increase of the size of the microfibrosis region leads to an increase of the length of a possible reentry circuit that is formed via a zig-zag propagation in the maze. Moreover, the 3D results reported here show that there are non-trivial and non-intuitive relations between the probability of reentries, the excitability of the tissue (APD and CV) and the thickness of the tissue.

The probability of reentry depends on the remodeling of the tissue and we observe a non-monotonic behavior: probability of reentry first increases when we move from parameter values near propagation boundary (strong remodeling) to weakly excitable tissues (intermediate remodeling); then, the probability of reentry decreases when we further move from weakly excitable to normal remodeling. Non-trivial relations are also observed varying the tissue thickness. For strong remodeling, probability of reentry first increases as we increase the thickness of the tissue from a monolayer to 0.3 mm, but decreases as we further enlarge the thickness from 0.5 mm to 2 mm.

Source-sink mismatches produces unidirectional propagation block [55] and have an important role in the generation of the zigzag reentry [18]. The probability of source-sink mismatch increases with size and thickness. Therefore, by further increasing the thickness of the system, the source-sink mismatches are too strong and instead of generating unidirectional blocks, generate full propagation block. The fine balance between having enough driving force for wave propagation, I_Na, and enough source-sink mismatches to generate reentry zigzag circuits explains the trends observed in the results presented in Fig 6. A similar balance or tradeoff between tissue structure and I_Na amplitude was also used earlier to explain slow conduction velocities observed in experiments [56–58].

The role of fibrosis: triggering vs. perpetuating arrhythmic activity

There are several attempts to model the interaction of waves within regions of heterogeneous cardiac tissue, where cells are partially disconnected [59–61]. Their main results are a decrease of the macroscopic conduction velocity [51, 62] and a decrease of excitability which may preclude spiral breakup [26]. Related results are the observations that propagating waves perform a zigzag motion associated with the presence of non-excitable tissue [63], scars [64], propagation inside vein tissue [20] or microfibrosis [13].

However, it is worth noting that the above mentioned studies consider the role of fibrosis in the maintenance or perpetuation of arrhythmia, e. g. atrial fibrillation. To achieve this goal, the arrhythmia is always induced by premature stimulation (S1–S2 protocol) or fast pacing in the simulations. In this work, our focus is on the role of fibrosis as a trigger for arrhythmia. Our simulations present how a region with microfibrosis can sustain microreentry and become an ectopic pacemaker. For instance, Fig 8 shows that after a single stimulus the 3D tissue is rhythmically paced by the region of microfibrosis, with a cycle length around 200 ms. Future work may compare the outcome of different simulation protocols for a given heterogeneous setup that models regions of fibrosis inside a healthy tissue. We expect that reentry occurs earlier in the S1–S2 and fast pacing protocols than in our single-stimulus approach. Nevertheless, the parameter regions where reentry is found to be most likely to occur in our study here should also be the ones where propagation of multiple excitation waves is most vulnerable to the emergence of reentries. This is the first 3D simulations that support the hypothesis that ectopic pacemakers can be explained by the mechanism of microreentry inside a region of fibrosis. Recently, the role of fibrosis and other tissue heterogeneities in the generation of ectopic foci was also studied using different mechanisms of arrhythmia, such as trigger activity [65] and abnormal automaticity [66].

We systematically change the thickness of our 3D model. Fig 6 shows that the results are qualitatively comparable when using the monolayer or a 0.3 mm thick tissue. However, both the qualitative and quantitative results obtained with 0.5 mm to 2 mm-thick tissues are substantially different from the results obtained with a monolayer model indicating a crossover from genuine 2D to 3D behavior. This finding is in accordance with recent experimental and simulation-based evidence of the impact of endo-epicardial electrical dissociation during atria fibrillation [67]. During atrial fibrillation, the electrical activity of endocadial and epicardial surfaces is found to be progressively uncoupled, most likely due to interstitial fibrosis. In addition, the uncoupling is correlated with increasing stability and complexity of the atrial fibrillation substrate. Therefore, one may conclude that atrial fibrillation progressively becomes a three-dimensional process.

Fractionated and asymmetric electrograms

Under normal conditions, surface EG show first a positive deflection, when the AP wave is arriving, followed by a negative deflection, indicating the wave is moving away. The amplitudes of these reflections are similar or symmetric. Such typical waveforms are observed, for instance, in our results presented in Fig 8. However, previous computational studies have already characterized how changes in electrogram waveform are associated with complex and fractionated propagation in the presence of fibrosis, arrhythmia, and complex 3D geometry and anisotropy [11, 13, 15, 18]. For instance, in Fig 8(f) and 8(g), the waveforms are asymmetric, with a positive deflection bigger than the following negative one. This indicates the collision of the two propagating waves, as described before in [68]. The two waves have to propagate around the fibrotic region before they collide. A negative monophasic electrogram, as in the electrode E of Fig 8(e), is usually associated with an AP focus [11, 18], since this would reflect that most of the waves propagate away for the electrode position.

Fractionated electrograms are carefully studied and related to fractionated propagation in microfibrosis in [13, 15, 18]. In our simulations we observe fractionated and complex electrograms near the fibrotic region, as presented by the waveforms in Fig 8(c)–8(e).

Limitations, future works and conclusions

The ionic model used in the numerical simulations is simple and rather qualitative. In order to have a direct control of the parameters related with ion channels and their possible pathological or pharmacokinetic modifications, more elaborate and detailed ionic models should be employed in future studies. However previous results in two dimensions with the same model [23] were later found to be very much in line with the outcome of studies with a more realistic and physiologically detailed model [24]. Therefore we are confident that the three-dimensional results presented here should server as a good starting point and provide guidance for studies in the near future involving more realistic models incorporating detailed electrophysiology and anatomy as well as experimentally obtained characterization of fibrosis.

The geometry employed here, a three-dimensional thin slab of tissue, is relatively simple and we have already ongoing works that consider anisotropy and more accurate geometric descriptions. The number of neighboring cells in the cubic lattice employed here is six while a more realistic organization of the cells is probably closer to a three-dimensional close-packed hexagonal organization in which each cell has twelve nearest neighbors. Since the percolation threshold in such a lattice is substantially larger than the one in cubic lattice, reentries will also occur at a much larger value of ϕ in the hexagonal lattice than in the cubic lattice.

The development of more elaborate descriptions of heterogeneous cardiac tissue including specific information of this type is an actual research direction of several groups with a lot of implications in cardiac modeling for the correct description of fibrotic regions in cardiac tissue. For example, reactive interstitial fibrosis has been modeled by introducing non-conductive materials that disconnect a fraction ϕ of neighboring myocytes, however, fibroblast are also frequently present in such regions which can be electrically coupled to myocytes and affect their properties [43, 69–72].

In summary, we find that a confined region with microfibrosis may generate microreentries. Such active regions may act as ectopic pacemakers for the surrounding tissue, exhibit the electrogram signature of CFAEs and are intimately linked to the fraction of fibrosis, ϕ, and topology of connections between myocytes. The essential parameter characterizing the topology of the lattice of discrete nodes in the simulation (or the network of cells in experiments) is the associated percolation threshold ϕ_c. The phenomena reported here: microreentry inside fibrotic regions, macroreentry and ectopic pacemakers triggered from fibrotic areas and the associated complex fraction electrograms all occur if the fraction of broken links representing the degree of fibrosis is such that the network of the remaining conducting links is just above the percolation threshold.

Supporting Information

[Figure omitted. See PDF.]

S1 Movie. Reentry in a 3D slab of cardiac tissue.

Numerical simulation showing the reentry from the fibrotic region into the whole slab of tissue.

https://doi.org/10.1371/journal.pone.0166972.s001

(AVI)

Author Contributions

1. Conceptualization: SA RS MB.

2. Data curation: SA.

3. Formal analysis: SA.

4. Funding acquisition: SA RS MB.

5. Investigation: SA RS MB.

6. Methodology: SA RS MB.

7. Writing – original draft: SA RS MB.

8. Writing – review & editing: SA RS MB.

Citation: Alonso S, dos Santos RW, Bär M (2016) Reentry and Ectopic Pacemakers Emerge in a Three-Dimensional Model for a Slab of Cardiac Tissue with Diffuse Microfibrosis near the Percolation Threshold. PLoS ONE 11(11): e0166972. https://doi.org/10.1371/journal.pone.0166972

References

1. Panfilov A, Hogeweg P. Spiral breakup in a modified FitzHugh-Nagumo model. Physics Letters A. 1993;176(5):295–299.

2. Karma A. Spiral breakup in model equations of action potential propagation in cardiac tissue. Physical review letters. 1993;71(7):1103. pmid:10055449

3. Bär M, Eiswirth M. Turbulence due to spiral breakup in a continuous excitable medium. Physical Review E. 1993;48(3):R1635. pmid:9960866

4. Fenton FH, Cherry EM, Hastings HM, Evans SJ. Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity. Chaos: An interdisciplinary journal of Nonlinear Science. 2002;12(3):852–892. pmid:12779613

5. Qu Z, Garfinkel A, Weiss JN. Vulnerable window for conduction block in a one-dimensional cable of cardiac cells, 1: single extrasystoles. Biophysical journal. 2006;91(3):793–804. pmid:16679367

6. Alonso S, Bär M, Echebarria B. Nonlinear physics of electrical wave propagation in the heart: a review. Reports on Progress in Physics. 2016;79(9):096601. pmid:27517161

7. Pastore JM, Girouard SD, Laurita KR, Akar FG, Rosenbaum DS. Mechanism linking T-wave alternans to the genesis of cardiac fibrillation. Circulation. 1999;99(10):1385–1394. pmid:10077525

8. Gaztañaga L, Marchlinski FE, Betensky BP. Mechanisms of cardiac arrhythmias. Revista Española de Cardiología (English Edition). 2012;65(2):174–185.

9. Narayan SM, Krummen DE, Rappel WJ. Clinical mapping approach to diagnose electrical rotors and focal impulse sources for human atrial fibrillation. Journal of cardiovascular electrophysiology. 2012;23(5):447–454. pmid:22537106

10. Ryu K, Shroff SC, Sahadevan J, Martovitz NL, Khrestian CM, Stambler BS. Mapping of Atrial Activation During Sustained Atrial Fibrillation in Dogs with Rapid Ventricular Pacing Induced Heart Failure: Evidence for a Role of Driver Regions. Journal of Cardiovascular Electrophysiology. 2005;16(12):1348–1358. pmid:16403068

11. Tobón C, Ruiz-Villa CA, Heidenreich E, Romero L, Hornero F, Saiz J. A three-dimensional human atrial model with fiber orientation. Electrograms and arrhythmic activation patterns relationship. PLOS One. 2013; p. e50883. pmid:23408928

12. Haissaguerre M, Shah AJ, Cochet H, Hocini M, Dubois R, Efimov I, et al. Intermittent drivers anchoring to structural heterogeneities as a major pathophysiologic mechanism of human persistent atrial fibrillation. The Journal of physiology. 2016;. pmid:26890861

13. Jacquemet V, Henriquez CS. Genesis of complex fractionated atrial electrograms in zones of slow conduction: A computer model of microfibrosis. Heart Rhythm. 2009;6(6):803–810. http://dx.doi.org/10.1016/j.hrthm.2009.02.026. pmid:19467508

14. Kim AM, Olgin JE. Microfibrosis and complex fractionated atrial electrograms. Heart Rhythm. 2009;6(6):811–812. http://dx.doi.org/10.1016/j.hrthm.2009.03.006. pmid:19467509

15. Campos FO, Wiener T, Prassl AJ, Weber dos Santos R, Sanchez-Quintana D, Ahammer H, et al. Electroanatomical Characterization of Atrial Microfibrosis in a Histologically Detailed Computer Model. Biomedical Engineering, IEEE Transactions on. 2013;60(8):2339–2349. pmid:23559023

16. Spach MS, Heidlage JF, Dolber PC, Barr RC. Mechanism of origin of conduction disturbances in aging human atrial bundles: Experimental and model study. Heart Rhythm. 2007;4(2):175–185. http://dx.doi.org/10.1016/j.hrthm.2006.10.023. pmid:17275753

17. Spach MS, Dolber PC, Heidlage JF. Influence of the passive anisotropic properties on directional differences in propagation following modification of the sodium conductance in human atrial muscle. A model of reentry based on anisotropic discontinuous propagation. Circulation Research. 1988;62(4):811–32. pmid:2450697

18. Hubbard ML, Henriquez CS. A microstructural model of reentry arising from focal breakthrough at sites of source-load mismatch in a central region of slow conduction. American Journal of Physiology-Heart and Circulatory Physiology. 2014;306(9):H1341–H1352. pmid:24610922

19. Vigmond E, Pashaei A, Amraoui S, Cochet H, Hassaguerre M. Percolation as a mechanism to explain atrial fractionated electrograms and reentry in a fibrosis model based on imaging data. Heart Rhythm. 2016;. pmid:26976038

20. Cherry EM, Ehrlich JR, Nattel S, Fenton FH. Pulmonary vein reentry—properties and size matter: insights from a computational analysis. Heart Rhythm. 2007;4(12):1553–1562. pmid:18068635

21. Rutherford SL, Trew ML, Sands GB, LeGrice IJ, Smaill BH. High-Resolution 3-Dimensional Reconstruction of the Infarct Border Zone: Impact of Structural Remodeling on Electrical Activation. Circulation Research. 2012;111(3):301–311. pmid:22715470

22. Smaill BH, Zhao J, Trew ML. Three-Dimensional Impulse Propagation in Myocardium: Arrhythmogenic Mechanisms at the Tissue Level. Circulation Research. 2013;112(5):834–848. pmid:23449546

23. Alonso S, Bär M. Reentry Near the Percolation Threshold in a Heterogeneous Discrete Model for Cardiac Tissue. Physical Review Letters. 2013;110(15):158101. pmid:25167313

24. Gouvêa de Barros B, Weber dos Santos R, Lobosco M, Alonso S. Simulation of Ectopic Pacemakers in the Heart: Multiple Ectopic Beats Generated by Reentry inside Fibrotic Regions. BioMed Research International. 2015;2015:713058. pmid:26583127

25. de Bakker JM, van Capelle FJ, Janse MJ, Tasseron S, Vermeulen JT, de Jonge N, et al. Slow conduction in the infarcted human heart.’Zigzag’ course of activation. Circulation. 1993;88(3):915–26. pmid:8353918

26. Panfilov A. Spiral breakup in an array of coupled cells: the role of the intercellular conductance. Physical Review Letters. 2002;88(11):118101. pmid:11909429

27. Smith J, Ritzenberg A, Cohen R. Percolation theory and cardiac conduction. Computers in Cardiology. 1984;11.

28. Kaplan DT, Smith JM, Saxberg BE, Cohen RJ. Nonlinear dynamics in cardiac conduction. Mathematical biosciences. 1988;90(1):19–48. pmid:11539069

29. Christensen K, Manani KA, Peters NS. Simple model for identifying critical regions in atrial fibrillation. Physical review letters. 2015;114(2):028104. pmid:25635565

30. Kazbanov IV, Ten Tusscher KH, Panfilov AV. Effects of Heterogeneous Diffuse Fibrosis on Arrhythmia Dynamics and Mechanism. Scientific reports. 2016;6. pmid:26861111

31. Kalifa J, Tanaka K, Zaitsev AV, Warren M, Vaidyanathan R, Auerbach D, et al. Mechanisms of Wave Fractionation at Boundaries of High-Frequency Excitation in the Posterior Left Atrium of the Isolated Sheep Heart During Atrial Fibrillation. Circulation. 2006;113(5):626–633. pmid:16461834

32. Nademanee K, McKenzie J, Kosar E, Schwab M, Sunsaneewitayakul B, Vasavakul T, et al. A new approach for catheter ablation of atrial fibrillation: mapping of the electrophysiologic substrate. Journal of the American College of Cardiology. 2004;43(11):2044–2053. http://dx.doi.org/10.1016/j.jacc.2003.12.054. pmid:15172410

33. Nademanee K, Lockwood E, Oketani N, Gidney B. Catheter ablation of atrial fibrillation guided by complex fractionated atrial electrogram mapping of atrial fibrillation substrate. Journal of Cardiology. 2010;55(1):1–12. http://dx.doi.org/10.1016/j.jjcc.2009.11.002. pmid:20122543

34. Kumagai K, Sakamoto T, Nakamura K, Nishiuchi S, Hayano M, Hayashi T, et al. Combined Dominant Frequency and Complex Fractionated Atrial Electrogram Ablation After Circumferential Pulmonary Vein Isolation of Atrial Fibrillation. Journal of Cardiovascular Electrophysiology. 2013;24(9):975–983. pmid:23691953

35. Chen J, Lin Y, Chen L, Yu J, Du Z, Li S, et al. A decade of complex fractionated electrograms catheter-based ablation for atrial fibrillation: Literature analysis, meta-analysis and systematic review. IJC Heart and Vessels. 2014;4(0):63–72. http://dx.doi.org/10.1016/j.ijchv.2014.06.013.

36. Magnani S, Muser D, Chik W, Santangeli P. Adjunct ablation strategies for persistent atrial fibrillation—beyond pulmonary vein isolation. Journal of thoracic disease. 2015;7(2):178. pmid:25713734

37. Podziemski P, Zebrowski JJ. A simple model of the right atrium of the human heart with the sinoatrial and atrioventricular nodes included. Journal of clinical monitoring and computing. 2013;27(4):481–498. pmid:23430363

38. Dössel O, Krueger MW, Weber FM, Wilhelms M, Seemann G. Computational modeling of the human atrial anatomy and electrophysiology. Medical & biological engineering & computing. 2012;50(8):773–799. pmid:22718317

39. Hansen BJ, Zhao J, Csepe TA, Moore BT, Li N, Jayne LA, et al. Atrial fibrillation driven by micro-anatomic intramural re-entry revealed by simultaneous sub-epicardial and sub-endocardial optical mapping in explanted human hearts. European Heart Journal. 2015;36(35):2390–2401. pmid:26059724

40. Colman MA, Aslanidi OV, Kharche S, Boyett MR, Garratt C, Hancox JC, et al. Pro-arrhythmogenic effects of atrial fibrillation-induced electrical remodelling: insights from the three-dimensional virtual human atria. The Journal of physiology. 2013;591(17):4249–4272. pmid:23732649

41. Dobrev D, Ravens U. Remodeling of cardiomyocyte ion channels in human atrial fibrillation. Basic research in cardiology. 2003;98(3):137–148. pmid:12883831

42. Campos FO, Prassl AJ, Seemann G, Weber dos Santos R, Plank G, Hofer E. Influence of ischemic core muscle fibers on surface depolarization potentials in superfused cardiac tissue preparations: a simulation study. Medical & Biological Engineering & Computing. 2012;50(5):461–472. pmid:22411321

43. Brown TR, Krogh-Madsen T, Christini DJ. Computational Approaches to Understanding the Role of Fibroblast-Myocyte Interactions in Cardiac Arrhythmogenesis. BioMed Research International. 2015;2015. pmid:26601107

44. Angel N, Li L, Macleod RS, Marrouche N, Ranjan R, Dosdall DJ. Diverse Fibrosis Architecture and Premature Stimulation Facilitate Initiation of Reentrant Activity Following Chronic Atrial Fibrillation. Journal of Cardiovascular Electrophysiology. 2015;26(12):1352–1360. pmid:26249367

45. Keener JP, Sneyd J. Mathematical physiology. vol. 1. Springer; 1998.

46. Clayton R, Bernus O, Cherry E, Dierckx H, Fenton F, Mirabella L, et al. Models of cardiac tissue electrophysiology: progress, challenges and open questions. Progress in Biophysics and Molecular Biology. 2011;104(1):22–48. pmid:20553746

47. Fenton F, Karma A. Fiber-rotation-induced vortex turbulence in thick myocardium. Physical Review Letters. 1998;81(2):481.

48. Krueger MW, Rhode KS, O‘Neill MD, Rinaldi CA, Gill J, Razavi R, et al. Patient-specific modeling of atrial fibrosis increases the accuracy of sinus rhythm simulations and may explain maintenance of atrial fibrillation. Journal of electrocardiology. 2014;47(3):324–328. pmid:24529989

49. Oliveira RS, de Barros BG, Gomes JM, Lobosco M, Alonso S, Bär M, et al. Reactive Interstitial and Reparative Fibrosis as Substrates for Cardiac Ectopic Pacemakers and Reentries. In: Bioinformatics and Biomedical Engineering. Springer; 2016. p. 346–357. https://doi.org/10.1007/978-3-319-31744-1_31

50. Torquato S. Random heterogeneous materials: microstructure and macroscopic properties. vol. 16. Springer Science & Business Media; 2002. https://doi.org/10.1115/1.1483342

51. Alonso S, Kapral R, Bär M. Effective medium theory for reaction rates and diffusion coefficients of heterogeneous systems. Physical Review Letters. 2009;102(23):238302. pmid:19658977

52. Alonso S, Bär M. Reentry produced by small-scale heterogeneities in a discrete model of cardiac tissue. In: Journal of Physics: Conference Series. vol. 727. IOP Publishing; 2016. p. 012002.

53. Malmivuo J, Plonsey R. Bioelectromagnetism: principles and applications of bioelectric and biomagnetic fields. Oxford university press; 1995.

54. Zahid S, Cochet H, Boyle PM, Schwarz EL, Whyte KN, Vigmond EJ, et al. Patient-derived models link re-entrant driver localization in atrial fibrillation to fibrosis spatial pattern. Cardiovascular research. 2016;110(3):443–454. pmid:27056895

55. Prudat Y, Kucera JP. Nonlinear behaviour of conduction and block in cardiac tissue with heterogeneous expression of connexin 43. Journal of molecular and cellular cardiology. 2014;76:46–54. pmid:25128085

56. Rohr S, Kucera JP, Kléber AG. Slow Conduction in Cardiac Tissue, I: Effects of a Reduction of Excitability Versus a Reduction of Electrical Coupling on Microconduction. Circulation Research. 1998;83(8):781–794. pmid:9776725

57. Fast VG, Kléber AG. Cardiac tissue geometry as a determinant of unidirectional conduction block: assessment of microscopic excitation spread by optical mapping in patterned cell cultures and in a computer model. Cardiovascular Research. 1995;29(5):697–707. pmid:7606760

58. Auerbach DS, Grzeda KR, Furspan PB, Sato PY, Mironov S, Jalife J. Structural heterogeneity promotes triggered activity, reflection and arrhythmogenesis in cardiomyocyte monolayers. The Journal of Physiology. 2011;589(9):2363–2381. pmid:21486795

59. Bub G, Shrier A, Glass L. Spiral wave generation in heterogeneous excitable media. Physical review letters. 2002;88(5):058101. pmid:11863783

60. Tanaka K, Zlochiver S, Vikstrom KL, Yamazaki M, Moreno J, Klos M, et al. Spatial distribution of fibrosis governs fibrillation wave dynamics in the posterior left atrium during heart failure. Circulation research. 2007;101(8):839–847. pmid:17704207

61. Mendonca Costa C, Campos FO, Prassl AJ, Weber dos Santos R, Sánchez-Quintana D, Ahammer H, et al. An efficient finite element approach for modeling fibrotic clefts in the heart. Biomedical Engineering, IEEE Transactions on. 2014;61(3):900–910. pmid:24557691

62. Ten Tusscher KH, Panfilov AV. Influence of diffuse fibrosis on wave propagation in human ventricular tissue. Europace. 2007;9(suppl 6):vi38–vi45. pmid:17959692

63. Casaleggio A, Hines ML, Migliore M. Computational model of erratic arrhythmias in a cardiac cell network: the role of gap junctions. PLOS One. 2014; p. e100288. pmid:24941068

64. de Sa DDC, Thompson N, Stinnett-Donnelly J, Znojkiewicz P, Habel N, Müller JG, et al. Electrogram fractionation the relationship between spatiotemporal variation of tissue excitation and electrode spatial resolution. Circulation: Arrhythmia and Electrophysiology. 2011;4(6):909–916. pmid:21984446

65. Zimik S, Nayak AR, Pandit R. A Computational Study of the Factors Influencing the PVC-Triggering Ability of a Cluster of Early Afterdepolarization-Capable Myocytes. PloS one. 2015;10(12):e0144979. pmid:26675670

66. Pumir A, Arutunyan A, Krinsky V, Sarvazyan N. Genesis of ectopic waves: role of coupling, automaticity, and heterogeneity. Biophysical journal. 2005;89(4):2332–2349. pmid:16055545

67. Eckstein J, Maesen B, Linz D, Zeemering S, van Hunnik A, Verheule S, et al. Time course and mechanisms of endo-epicardial electrical dissociation during atrial fibrillation in the goat. Cardiovascular research. 2010; p. cvq336. pmid:20978006

68. Jacquemet V, Virag N, Ihara Z, Dang L, Blanc O, Zozor S, et al. Study of unipolar electrogram morphology in a computer model of atrial fibrillation. Journal of cardiovascular electrophysiology. 2003;14(s10):S172–S179. pmid:14760921

69. McDowell KS, Arevalo HJ, Maleckar MM, Trayanova NA. Susceptibility to Arrhythmia in the Infarcted Heart Depends on Myofibroblast Density. Biophysical Journal. 2011;101(6):1307–1315. http://dx.doi.org/10.1016/j.bpj.2011.08.009. pmid:21943411

70. Miragoli M, Salvarani N, Rohr S. Myofibroblasts Induce Ectopic Activity in Cardiac Tissue. Circulation Research. 2007;101(8):755–758. pmid:17872460

71. Jacquemet V, Henriquez CS. Loading effect of fibroblast-myocyte coupling on resting potential, impulse propagation, and repolarization: insights from a microstructure model. American Journal of Physiology—Heart and Circulatory Physiology. May 2008;294(5):H2040–H2052. pmid:18310514

72. Xie Y, Garfinkel A, Camelliti P, Kohl P, Weiss JN, Qu Z. Effects of fibroblast-myocyte coupling on cardiac conduction and vulnerability to reentry: a computational study. Heart Rhythm. 2009;6(11):1641–1649. pmid:19879544