Reentry near the percolation threshold in a heterogeneous discrete model for cardiac tissue

Reentry Near the Percolation Threshold in a Heterogeneous Discrete Model for Cardiac Tissue Sergio Alonso and Markus Ba ¨r Physikalisch-Technische Bundesanstalt, Abbestrasse 2-12, 10587 Berlin, Germany (Received 23 March 2012; revised manuscript received 10 August 2012; published 9 April 2013) Arrhythmias in cardiac tissue are related to irregular electrical wave propagation in the heart. Cardiac tissue is formed by a discrete cell network, which is often heterogeneous. A localized region with a fraction of nonconducting links surrounded by homogeneous conducting tissue can become a source of reentry and ectopic beats. Extensive simulations in a discrete model of cardiac tissue show that a wave crossing a heterogeneous region of cardiac tissue can disintegrate into irregular patterns, provided the fraction of nonconducting links is close to the percolation threshold of the cell network. The dependence of the reentry probability on this fraction, the system size, and the degree of excitability can be inferred from the size distribution of nonconducting clusters near the percolation threshold. DOI: 10.1103/PhysRevLett.110.158101 PACS numbers: 87.19.Hh, 82.39.Rt, 87.10.e, 89.75.Kd Introduction.—Large excitation followed by a refractory state is characteristic of excitable media. Such dynamics controls many biological systems, e.g., electrical spikes in neurons, calcium waves in cells, or electric propagation through cardiac tissue [1]. In particular, electrical waves synchronize the contraction of the heart [2]. Abnormal propagation in cardiac tissue may give rise to fibrillation, i.e., chaotic electrical activity in the tissue, which desynchronizes the whole heart. Cardiac tissue is highly heterogeneous. Heterogeneity sizes range from the cell size (100m) to several millimeters (large blood vessels). Large heterogeneities have been recently found to be crucial for defibrillation [3]. However, cardiac arrhythmias are often successfully modeled by homogeneous models [4–6]. Such models are appropriate if heterogeneities appear on scales much smaller than the wave front width [7]. Cardiac tissue is formed by excitable cells called myocytes. Myocytes are inhomogeneously distributed and are connected by gap junctions. The extracellular matrix, nonconducting gap junctions, blood vessels, or fibroblasts may disturb the gap-junction coupling among the myocytes [2]. The presence of nonconducting connections in cardiac tissue can strongly affect wave propagation [8–11] and increases the risk of arrhythmias [12]. Heterogeneities in cardiac tissue provide a possible anatomical mechanism for reentry [12], whereas the breakup of spiral waves in homogeneous systems [6,13] is a prominent candidate for a functional mechanism for the chaotic dynamics observed in fibrillation [12]. For example, a rapid pacing of the homogeneous tissue may produce cardiac alternans and can eventually induce reentry [14]. However, the mechanisms underlying reentry and fibrillation are still not fully understood. As a control parameter, we consider the fraction of nonconducting links in the cellular network. ¼0 corresponds to completely homogeneous tissue and ¼1 to a set of disconnected myocytes. While for small fraction, a wave propagates through the system; for large fraction, a wave stops. In between, there is an intermediate value of the fraction pwhere both behaviors interchange. This fraction defines the percolation threshold [15]. For fractions just below p, larger clusters of blocking connections appear [16]. The window between this point and percolation turns out to be a dangerous parameter region which can be modeled only in a discrete, heterogeneous model of cardiac tissue. We find that a wave propagating in homogeneous tissue and interacting with a heterogenous region, which mimics a damaged region on the cardiac muscle, can generate reexcitations. The reexcitations may appear in the form of ectopic beats or as reentries. To obtain such reexcitations, the combination of two factors inside the damaged area are found to be essential: a fraction of nonconducting links close to percolation and a change of the local properties towards weaker excitability. We show that the probability of reentry is strongly correlated with the probability of formation of clusters above a typical size. Such clusters appear close to percolation. Model.—We perform simulations using many different realizations of random heterogeneous tissue with a given fraction of nonconducting links. In order to obtain statistically representative results, we keep our model as simple as possible and employ the Fenton-Karma model for wave propagation in cardiac tissue [17]. Here, the parameters of the model are chosen [17,18] in order to fit the propagation properties of the modified Beeler-Reuter model. The equation for the propagation of the action potential reads @tVi¼ðIfi þIso þIsiÞþX n j¼1 ijðVjViÞ(1) for the cell i, where is proportional to the conductivity and nis the number of first neighbors. Two adjacent cells may not be connected because of the loss of the gap-junction connection or by the formation of excess PRL 110, 158101 (2013) PHYSICAL REVIEW LETTERS week ending 12 APRIL 2013 0031-9007=13=110(15)=158101(5) 158101-1 Ó2013 American Physical Society connective tissue (fibrosis). Two neighboring cells iand j are connected if ij ¼1and disconnected if ij ¼0[19]. The potential Viis conveniently scaled to dimensionless variables (0<V i<1)[17]. Three ion currents are considered here Ifi ¼vpðViVcÞð1ViÞ=d; Iso ¼Við1pÞ=oþp=r; Isi ¼wð1þtanhðkðViVsi cÞÞÞ=ð2siÞ; (2) mimicking fast inward sodium (Ifi), slow outward potassium (Iso), and slow inward calcium (Isi) currents. The variables vand ware gating variables [20]. For normal physiological conditions [21], the model yields an action potential duration (APD) of 150 milliseconds (ms) for a single pulse. Fixing ¼D=x2, where Dis proportional to the conductivity [17], the model correctly reproduces the characteristic scale of wave propagation, i.e., width of wave front 500 mand an action potential width of 6 cm. Interaction with a damaged region.—In Fig. 1(a) we show a typical process to generate reentry in a homogeneous excitable medium [22]. Figure 1(b) shows an alternative mechanism of wave breakup based on the interaction with an inhomogeneous region. We block a fraction of randomly picked connections in the heterogeneous region. The wave may break inside the region and send new waves into the normal and homogeneous tissue, see Fig. 1(b). For the complete subsequent spatiotemporal evolution along the diagonal of the simulation box see Fig. 1(c). We restrict our study to consider dynamics caused by a single action potential wave, in comparison with previous studies [9–11] where the interaction between consecutive waves produced the reentry. The wave dynamics depends on the random realization of the connectivity matrix ij. We carried out simulations for 100 realizations of the random heterogeneous network for a given value of . Histograms approximating the reentry probability PðÞare shown in Fig. 2. We classify as reentry any simulation with wave activity after 1.5 s whereas the first excitation wave typically leaves the system after 0.4 s. For small , the wave crosses the heterogeneous region with a small delay. For large , the wave does not enter into the heterogeneous region. Reentries are observed for intermediate values of , and the probability of their occurrence increases with the size of the damaged area, see Fig. 2(a). The reentry probability also increases with d, see Fig. 2(b). The increase of dis equivalent to a decrease of the maximal conductance for the sodium current in more detailed models [23], i.e., a decrease of excitability (see Fig. S1 [24]). This dependency of reentry probabilities on system size and parameter dmay be explained by the hypothesis that the occurrence of reentry requires a cluster larger than a certain critical size. With increasing size of the damaged region, the probability that such a cluster is contained in the system becomes larger. Moreover, a decrease of excitability, i.e., an increase of d, reduces the width of the action potential wave. Therefore, it is plausible that the critical size for a cluster to cause reentry also decreases. Propagation in completely heterogeneous systems.— Whileachangeofdmay correspond to an electrical remodeling of the properties of the tissue, the heterogeneities mimic a structural remodeling of the tissue. To study the effects of the nonconducting links more systematically, we focus on the dynamics in the damaged region and consider a uniformly heterogeneous system. An increase of  decreases the propagation velocity (see Fig. S2 [24]), while it does not substantially modify the APD. For a value of  above the percolation threshold p, the velocity of the waves is predicted to decrease to zero [7]. Only for < pthe wave propagates. However, for certain fractions < p, wave propagation through the system is irregular and local conduction block occurs. Wave activity then often gives rise to reentry and sustained irregular activity, see Fig. 3. FIG. 1 (color online). (a) Spiral wave generation in homogeneous media (77cm2) by an external perturbation. (b) Reentry by interaction of a wave under normal conditions (d¼ 0:25 ms) with a circular region with high fraction of nonconducting links indicated by a dashed circumference (radius R0¼0:8cm,¼0:49 and d¼0:35ms), total area 77cm 2. (c) Spatiotemporal plot along the diagonal line shown in first panel of (b). Full vertical lines indicate times corresponding to snapshots in panel (b); dashed lines indicate the region with nonconducting links. PRL 110, 158101 (2013) PHYSICAL REVIEW LETTERS week ending 12 APRIL 2013 158101-2 The appearance of reentries is related with the existence of large clusters of heterogeneities, compare Figs. 3(a) and 3(b) where a large nonconducting cluster is the cause of reentry. The formation of large nonconducting clusters requires a proximity to the percolation threshold. This qualitative argument provides a link between the occurrence of reentry and the percolation threshold, and related properties like the dependence of the maximum cluster on . To check this hypothesis, we have calculated PðÞfor bond and site percolation [19] in square and hexagonal topologies, see Fig. 4. In all four cases, we obtain qualitatively similar results. Reentry is observed for values of the fraction of nonconducting links or nonconducting sites  that are 5%–10% smaller than the respective percolation threshold p. Relation between size of nonconducting clusters and reentry.—We return to the argument that reentry is caused by clusters with a diameter larger than a critical length ‘. We generated different random realizations of the grid to calculate numerically the cluster size distribution. Using 103realizations for each value of , we calculate numerically the probability L ‘of finding a cluster larger than a critical length ‘in a system of size Land the probability of percolation (L L) corresponding to the case ‘¼L. Then, the function L xðÞ¼ð1þtanh½ð xÞ=xÞ=2is fitted to the results for x¼‘and x¼L (see Fig. S3 [24]). We hypothesize that the probability of reentry PðÞis strongly correlated with the probability 0.44 0.46 0.48 0.5 0.52 φ 0 0.1 0.2 0.3 0.4 0.5 0.6 P ( φ ) 0.44 0.46 0.48 0.5 0.52 φ 0 0.1 0.2 0.3 0.4 0.5 0.6 P ( φ ) P ( φ ) P ( φ ) P ( φ ) P ( φ ) 0.44 0.46 0.48 0.5 0.52 φ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.44 0.46 0.48 0.5 0.52 φ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.44 0.46 0.48 0.5 0.52 φ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.44 0.46 0.48 0.5 0.52 φ 0 0.1 0.2 0.3 0.4 0.5 0.6 R0=1.8cm R0=1.4cm R0=0.8cm τd=0.35ms τd=0.325ms τd=0.3ms (a) (b) FIG. 2 (color online). Dependence of the reentry probability on for different radius keeping d¼0:35 ms inside the region (a) and for different values of dinside the region keeping Ro¼1:8cm (b). Dashed lines indicate p. Solid lines are Gaussian fits. FIG. 3 (color online). Panel (a) shows clusters (dark) into the conducting area (bright). White square marks a large cluster responsible for the reentry. Panels (b)–(d) show the wave breakup in the heterogeneous network. Straight arrow in (b) marks the place of reentry and circular arrow in (c) follows the reentry direction. Parameter values: area 44cm2,¼0:49, and d¼0:3ms. 0.46 0.48 0.5 0.52 0.54 φ 0 0.05 0.1 0.15 0.2 0.25 P( φ ) P( φ )P( φ ) P( φ ) 0.6 0.62 0.64 0.66 0.68 0.7 φ 0 0.05 0.1 0.15 0.2 0.25 0.34 0.36 0.38 0.4 0.42 0.44 φ 0 0.05 0.1 0.15 0.2 0.25 0.44 0.46 0.48 0.5 0.52 φ 0 0.05 0.1 0.15 0.2 0.25 (b) (a) (c) (d) FIG. 4 (color online). (color online). Dependence of the reentry probability on the fraction of nonconducting links, keeping d¼0:3ms, for square (a) and hexagonal (b) networks, and nonconducting sites for square (c) and hexagonal (d) networks. Dashed lines indicate percolation threshold: p¼0:5(a), p¼0:41 (b), p¼0:65 (c), and p¼0:5(d). Solid lines are Gaussian fits. PRL 110, 158101 (2013) PHYSICAL REVIEW LETTERS week ending 12 APRIL 2013 158101-3 QðÞto find a cluster of size larger than ‘in the absence of percolation, which would completely block wave propagation and prevent reentry. We assume that QðÞ is well approximated by the probability of having a large enough cluster multiplied with the probability that the system is below percolation: QðÞ¼L ‘ðÞð1L LðÞÞ (see Fig. S3 [24]). In Fig. 5, we consider the case of bond percolation in a square network. The dependence of QðÞon for different values of ‘and Lis shown. Since we have stipulated that the minimal cluster size ‘to generate reentry monotonically decreases with the excitability ‘ðdÞ, Fig. 5shows that reentry becomes more likely for larger system size L [Fig. 5(a)] and for weaker excitability [Fig. 5(b)]. The results are consistent with the outcome of the numerical simulations in Fig. 2. Furthermore, Fig. 5(b) exhibits a shift of the maximum of QðÞto the left of the percolation threshold for decreasing ‘, which correlates with the equivalent shift of PðÞobtained for increasing d observed in Fig. 2(b). Discussion.—When electrical waves propagate through cardiac tissue, their interaction with heterogeneous damaged areas of the tissue may lead to reentry. The reentry probability depends on the local physiological properties as well as on the size of the damaged area. In the numerical simulations, a decrease of the excitability, i.e., increase of d, increases the probability of reentry (see Fig. S4 [24]). For normal propagation conditions, the probability of reentry strongly decreases and wave propagation is robust, as expected for healthy hearts. Cardiac tissue is heterogeneous and different types of heterogeneities coexist. Detailed descriptions are available for structural and electrical remodeling of the cardiac tissue [12]. Moreover, similar heterogeneities like the ones used here were previously employed, e.g., to model blocked gap junctions [9], fibrosis [10], intercalation of fibroblast [25], or current sinks [26]. Although a completely heterogeneous cardiac tissue where waves break continuously seems unlikely, the presence of a highly heterogeneous region surrounded by normal tissue is quite feasible. We propose that such an anatomical situation can induce a reentry in cardiac tissue. It may also explain the appearance of ectopic beats where a new excitation pulse appears spontaneously in the tissue. In the atria, such regions may cause complex fractionated atrial electrograms, where a localized zone of the atria shows electrograms with high frequencies and unusually short spatial correlations while the surrounding areas have much lower frequencies [27]. Reentries appear for nonconducting cell-to-cell link fractions around 5%–10% below the percolation threshold. This threshold itself depends on the topology of the network. We consider four different cases with different topologies and find consistent results. This suggests that the knowledge of percolation properties of real cardiac tissue is critical for the prediction of heterogeneity-related (anatomical) cardiac arrhythmias. The hypothesis that reentry requires a certain critical size of a cluster of nonconducting links could be validated by a study of clustersize statistics near the percolation threshold. This approach explains the changes in reentry probability related to changes in system size and excitability. A dynamical control of the reentry generated by a damaged region is complicated because the induction of new waves may promote even more reentries. Our results indicate that the elimination of the reentries can be achieved by an increase of the fraction of nonconducting links to values above the percolation threshold. This corresponds to the most common therapy for atrial fibrillation; the ablation of small pathological regions of cardiac tissue. Further studies on more realistic cardiac models and geometries should give more quantitative results. However, we expect that the qualitative picture that emerged from our study is generic. The results should also carry over to other excitable systems. Chemical waves that interact with reactive droplets in oil [28] or catalyst-loaded beads [29,30] can break generating spiral waves [31]. Such experiments may be used to confirm our predictions. FIG. 5 (color online). (color online). Dependence of QðÞon for different sizes LLkeeping ‘¼3mm (a), and for different ‘keeping LL¼44cm 2(b). QðÞis calculated from L ‘and L L, which are fitted for each value of to the results obtained from 103random grid realizations (x¼0:01 cm). PRL 110, 158101 (2013) PHYSICAL REVIEW LETTERS week ending 12 APRIL 2013 158101-4 In conclusion, the vulnerability of the tissue to the appearance of heterogeneity induced reentries is determined by the proximity to the percolation threshold of the heterogeneous cell network combined with a decrease of the excitability and a critical size of the heterogeneous region. We are grateful to H. Koch and R. Weber dos Santos for valuable discussions. We acknowledge financial support from the German Science Foundation DFG within the framework of SFB 910 ‘‘Control of self-organizing nonlinear systems.’’ [1] J. P. Keener and J. Sneyd, Mathematical Physiology (Springer, New York, 1998). [2] D. P. Zipes and J. Jalife, Cardiac Electrophysiology: From Cell to Bedside (Saunders, Philadelphia, 2004). [3] S. Luther et al.,Nature (London) 475, 235 (2011). [4] A. Panfilov and P. Hogeweg, Phys. Lett. A 176, 295 (1993). [5] M. Ba ¨r and M. Eiswirth, Phys. Rev. E 48, 1635(R) (1993). [6] A. Karma, Phys. Rev. Lett. 71, 1103 (1993). [7] S. Alonso, R. Kapral, and M. Ba ¨r, Phys. Rev. Lett. 102, 238302 (2009). [8] J. G. Stinstra, R. MacLeod, and C. S. Henriquez, Ann. Biomed. Eng. 38, 1390 (2010). [9] E. M. Cherry, J. R. Ehrlich, S. Nattel, and F. H. Fenton, Heart Rhythm 4, 1553 (2007). [10] K. H. W. J. Ten Tusscher and A. V. Panfilov, Europace 9, vi38 (2007). [11] Z. Qu, H. S. Karagueuzian, A. Garfinkel, and J. N. Weiss, Am. J. Physiol. 286, H1310 (2004). [12] J. N. Weiss, Z. Qu, P.-S. Chen, S.-F. Lin, H. S. Karagueuzian, H. Hayashi, A. Garfinkel, and A. Karma, Circulation 112, 1232 (2005). [13] M. Gerhardt, H. Schuster, and J. J. Tyson, Science 247, 1563 (1990). [14] Z. Qu, A. Garfinkel, P.-S. Chen, and J. N. Weiss, Circulation 102, 1664 (2000). [15] S. Torquato, Random Heterogeneous Materials (Springer, New York, 2002). [16] D. Stauffer, Phys. Rep. 54, 1 (1979). [17] F. H. Fenton and A. Karma, Chaos 8, 20 (1998). [18] F. H. Fenton, E. M. Cherry, H. M. Hastings, and S. J. Evans, Chaos 12, 852 (2002). [19] For bond percolation, corresponds to the fraction of nonconducting links and the values of ij are uncorrelated. For site percolation, corresponds to the fraction of disconnected cells, and if i¼0then ij ¼08j. [20] The evolution of the gate variables depends on the action potential @tv¼ð1pÞð1vÞ= vðVÞpv=þ vand @tw¼ð1pÞð1wÞ= wpw=þ w, where  vðVÞ¼ ð1qÞ v1þq v2with p¼1if VVc, and p¼0if V<V c, and q¼1if VVv, and q¼0if V<V v. [21] Parameter values: d¼0:25 ms, v1¼19:6ms, v2¼ 1000 ms,þ v¼3:33ms, w¼11 ms,þ w¼667 ms, o¼8:3ms,r¼50ms,si ¼45ms,k¼10,Vv¼0:85, Vc¼0:13, and D¼0:001 cm2=ms. [22] We employ x¼100 mfor spatial discretization and second order Runge-Kutta method for temporal integration (t¼0:023 ms) of Eqs. (1) and (2). [23] S. Alonso and A. V. Panfilov, Chaos 17, 015102 (2007). [24] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.110.158101 for further detailed numerical results. [25] Y. Xie, A. Garfinkel, P. Camelliti, P. Kohl, J. Weiss, and Z. Qu, Heart Rhythm 6, 1641 (2009). [26] B. E. Steinberg, L. Glass, A. Shrier, and G. Bub, Phil. Trans. R. Soc. A 364, 1299 (2006). [27] K. Nademanee, J. McKenzie, E. Kosar, M. Schwab, B. Sunsaneewitayakul, T. Vasavakul, C. Khunnawat, and T. Ngarmukos, J. Am. Coll. Cardiol. 43, 2044 (2004). [28] M. Toiya, H. O. Gonza ´lez-Ochoa, V. K. Vanag, S. Fraden, and I. R. Epstein, J. Phys. Chem. Lett. 1, 1241 (2010). [29] J. Maselko, J. S. Reckley, and K. Showalter, J. Phys. Chem. 93, 2114 (1989). [30] A. F. Taylor, M. R. Tinsley, F. Wang, Z. Huang, and K. Showalter, Science 323, 614 (2009). [31] G. Wang, Q. Wang, P. He, S. Pullela, M. Marquez, and Z. Cheng, Phys. Rev. E 82, 045201(R) (2010). PRL 110, 158101 (2013) PHYSICAL REVIEW LETTERS week ending 12 APRIL 2013 158101-5