Self-organization of conducting pathways explains electrical wave propagation in cardiac tissues with high fibrosis

Cardiac fibrosis occurs in many forms of heart disease and is considered to be one of the main arrhythmogenic factors. Regions with a high density of fibrosis are likely to cause blocks of wave propagation that give rise to dangerous cardiac arrhythmias. Therefore, studies of the wave propagation through these regions are very important, yet the precise mechanisms leading to arrhythmia formation in fibrotic cardiac tissue remain poorly understood. Particularly, it is not clear how wave propagation is organized at the cellular level, as experiments show that the regions with a high percentage of fibrosis (65-75%) are still conducting electrical signals, whereas geometric analysis of randomly distributed cells predicts connectivity loss at 40% at the most (percolation threshold). To address this question, we used a joint in vitro-in silico approach, which combined experiments in neonatal rat cardiac monolayers with morphological and electrophysiological computer simulations. We have shown that the main reason for sustainable wave propagation in highly fibrotic samples is the formation of a branching network of cardiomyocytes. We have successfully reproduced the morphology of conductive pathways in computer modelling, assuming that cardiomyocytes align their cytoskeletons to fuse into cardiac syncytium. The electrophysiological properties of the monolayers, such as conduction velocity, conduction blocks and wave fractionation, were reproduced as well. In a virtual cardiac tissue, we have also examined the wave propagation at the subcellular level, detected wavebreaks formation and its relation to the structure of fibrosis and, thus, analysed the processes leading to the onset of arrhythmias. Author summary Cardiac arrhythmias are one of the major causes of death in the industrialized world. The most dangerous ones are often caused by the blocks of propagation of electrical signals. One of the common factors that contribute to the likelihood of these blocks, is a condition called cardiac fibrosis. In fibrosis, excitable cardiac tissue is partially replaced with the inexcitable connective tissue. The precise mechanisms leading to arrhythmia formation in fibrotic cardiac tissue remain poorly understood. Therefore, it is important to study wave propagation in fibrosis from cellular to tissue level. In this paper, we study fibrosis of high density in experiments and computer simulations. We have observed a paradoxical ability of the tissue with extremely high fibrosis (up to 75% of fibroblasts) to conduct electrical signals and contract synchronously, whereas geometric analysis of randomly distributed cells predicted connectivity loss at 40% at the most. To explain this phenomenon, we have studied the patterns that cardiac cells form in the tissue and reproduced their self-organisation in a computer model. Our virtual model also took into account the polygonal shapes of the spreading cells and explained high arrhythmogenicity of fibrotic tissue.

The contraction of the heart is controlled by propagating waves of excitation. Abnormal 2 regimes of the wave propagation may cause cardiac arrhythmia, asynchronous 3 contractions of the heart and even lead to cardiac arrest and sudden cardiac death. 4 Cardiac arrhythmias often originate from blocks of propagation [1]. In that case, the 5 wave goes around the block, reenters the same inceptive region and, thus, forms a 6 persistent rotational activity called cardiac reentry, which is one of the main 7 mechanisms of lethal cardiac arrhythmias. 8 A normal heart has a complex structure, which is composed of the bundles of 9 elongated cardiac cells. Apart from premier excitable cardiac cells, there are also 10 inexcitable cells of connective tissue: cardiac fibroblasts. Their role is to maintain the 11 structural integrity of the heart [2] and repair injuries [3]. Fibroblasts outnumber 12 cardiomyocytes in a healthy human heart although occupying a much smaller total 13 volume [2,4]. However, many pathological conditions are associated with an excessive 14 growth of the fibrous tissue, called cardiac fibrosis, which is, therefore, considered to be 15 one of the major arrhythmogenic factors [5,6]. 16 The mechanisms of arrhythmia onset in fibrotic tissue remain poorly understood but 17 generally believed to be associated with the increased probability of waveblock 18 formation. It is a well-established fact that fibrosis slows down wave propagation [7] and 19 can completely block it if the fibroblasts density is high. The critical density of fibrosis 20 above which the conduction terminates is called percolation threshold. This concept 21 originates from the percolation theory, and, by definition, it specifies the point of 22 long-range connectivity loss/formation in random systems. Connectivity here refers to 23 electrical synchronisation in the tissue, or, in other words, the ability to transmit 24 electrical waves of excitation. Percolation threshold, i.e. the critical density of fibrosis 25 which breaks long-range connectivity, plays an important role in arrhythmogenicity. It 26 was shown that cardiac tissue is most susceptible to arrhythmias if the density of 27 fibrosis is only slightly (∼2-3% [8]) below the percolation threshold. There are two main 28 factors that may facilitate reentry formation. First, a large amount of fibroblasts 29 (acting as heterogeneities) increases the probability for waveblock formation [9]. Second, 30 a high fraction of fibrosis creates a 'maze' that effectively lengthens the travel distance 31 for the waves as they follow a longer zig-zag path [10] and, thus, provides sufficient 32 room for the emerging reentrant loops. As a result, high density of fibrosis both 33 facilitates the initiation and creates conditions for the existence of reentrant cycles, 34 resulting in a highly arrhythmogenic substrate. 35 Up to now wave propagation failure was studied only in generic mathematical 36 representations of cardiac fibrosis with each cell randomly chosen to be either myocyte 37 or fibroblast (see e.g. [9,11]). In this kind of 2D computer models, the propagation of 38 excitation failed at 40% of fibrosis [12], which is also within the range of values 39 predicted by classic mathematical models (e.g. 37-44% of the area uncovered by 40 conducting elongated ellipses with the shape similar to cardiomyocytes [13]). However, 41 experimental measurements [14] indicate that wave propagation and synchronous contraction in 2D cardiac monolayers is observed for up to 75% percentage of fibrosis.

43
In this paper, we study the phenomenon of wave propagation in cardiac tissue with a 44 high density of fibrosis using a joint in vitro-in silico approach. We performed 45 experiments in 25 monolayers with various percentages of fibroblasts and detected wave 46 propagation to determine the percolation threshold. We have found that, indeed, the 47 experimentally measured threshold (75% of the area covered by fibroblasts) is 48 substantially higher than what was predicted in conventional computer modelling 49 (40% [12]) or classic mathematical models. Further morphological examination revealed 50 that the key mechanism of conduction in highly fibrotic tissue is likely to be tissue 51 patterning. The cardiomyocytes were not located randomly but organised in a 52 branching network that wired the whole sample. 53 Next, in order to explain cardiac network formation, we applied a virtual cardiac 54 monolayer framework developed in [15], based on the Cellular Potts Models [16][17][18]. We 55 proposed a hypothesis that such self-organisation occurs due to cytoskeletons' alignment. 56 Based on this hypothesis, we were able to obtain branching patterns, as well as 57 reproduce the decrease in conduction velocity and wave percolation observed in 58 experiments. We have further studied in silico the process of formation of the 59 wavebreaks leading to reentry formation and analysed the tissue structures that caused 60 them.

61
This paper is organized as follows. First, we describe the experiments conducted 62 with the neonatal cell cultures. Second, we analyse the patterns and formulate the 63 hypothesis on the mechanism of their formation. Third, we implement the hypothesis 64 and reproduce pattern formation in silico in a Cellular Potts Model (CPM) of cardiac 65 tissue [15]. Next, we compare electrical signal propagation in computer modelling and in 66 experiments. Finally, we show a spontaneous formation of uni-directional blocks in the 67 model resulting in reentry formation and locate and analyse the structures causing it. We cultured neonatal rat ventricular cardiomyocytes mixed with cardiac fibroblasts in 72 variable proportion (20-88% fibroblasts) and studied electrical activity in these 73 monolayers using optical mapping.

74
In Figure 1a and in S1 Video the wave propagation in a sample with high fibrosis is 75 shown (66% fibroblasts). In spite of high percentage of fibrosis, this sample is still 76 conducting electrical waves, however the wave propagation pattern is complex. The 77 wave originates from the stimulation point marked by a yellow spike at the bottom of 78 the tissue and initially spreads in all directions. However, due to a large number of 79 fibroblasts, the wave is blocked at multiple sites. After a short delay (in areas outlined 80 with dashed ellipses), the wave propagates further into the left and the right parts of 81 the sample, and again spreads in various directions. This process repeats multiple times, 82 resulting in a complex, fractionated wave pattern containing narrow pathways and some 83 bulk excited regions. The conduction blocks in Figure 1a are shown in red, which are 84 the places where the wave propagation was blocked and the wave had to go around.

85
The main propagation paths are shown with white and black arrows. In this sample, the 86 electrical wave propagation was still possible and long-range connectivity was still 87 present, however the amount of fibroblasts was close to the percolation threshold. 88 We have found that the percolation threshold for the neonatal rat cardiac 89 monolayers was 75 ± 2% of fibroblasts. We have also measured conduction velocities in 90 the samples below the percolation threshold ( Figure 1b). The measurements show that 91  the velocity decreased when approaching the percolation threshold. In the samples with 92 low fibrosis, the velocity was approximately 10-14 cm/s, and it decreased twofold in the 93 samples with 70% fibrosis. The number of conduction blocks was higher in samples with 94 high fibrosis. As a result, the mean conduction velocity decreased with the increase in a 95 percentage of fibroblasts.

96
After optical mapping of the wave propagation, we have fixated the samples and 97 studied their morphology using immunohistochemical labeling. We have found that the 98 cardiomyocytes in the samples have formed connected networks capable of electrical 99 wave propagation. In Figure 2 cardiomyocytes are shown in pink, and the cluster that 100 they have formed is outlined with a white contour. The cardiomyocytes were organised 101 in a branching structure that wired the whole sample. We have followed the pathway 102 using a confocal microscope and it was possible to find long-range connectivity in the 103 tissue. Therefore, we observed that the cardiomyocytes were organised in conduction 104 pathways and assumed that there must be a mechanism responsible for their 105 self-organisation.

107
Pattern formation in cell populations during development was extensively studied using 108 Cellular Potts Models [16,17,19]. However, after trying several approaches used before 109 (see the discussion section of our paper), we have not found an existing model that 110 could have been applicable to the cardiac tissue and at the same time could reproduce 111 the experimentally observed branching structures. Finally, after careful analysis of our 112 experimental preparation, we found that an essential feature of our structure was the  The precise biological mechanism of such alignment is unknown. However, in our 119 view, it can be derived from a well-known property of actin cytoskeleton reinforcement 120 in response to the external force [20]. Actin filaments, as well as the adaptors that link 121 them, remodel under applied tension. In case of contact of two cells, the actin filaments 122 are connected through adherens junctions, which transmit the tensile forces between the 123 filaments of the neighbouring cells. It was shown that higher tension stabilises the Cardiac syncytium in a 3-days culture of the neonatal rat cardiomyocytes. The cells have formed inercalated discs (ID, bright-green, highlighted with red arrows), aligned their cytoskeletons (the aligned strains on the both sides from ID are shown with white arrows) and formed a branching network. Nuclei are shown in blue (DAPI, labels DNA), and actin strands are shown in green (phalloidin, labels F-actin).
whole complex [20]. The tension is maximal, if the actin filaments are aligned with each 125 other, which gives a preference for intercellular alignment of the cytoskeletons. 126 Therefore, we have incorporated this mechanism into our Cellular Potts Model 127 (CPM) of cardiac cells [15]. This model was already adjusted to reproduce characteristic 128 shapes of the cardiac cells in virtual tissue model, and here we extended it with a new 129 energy term, that corresponded to the alignment of actin bundles in the neighbouring 130 cells.

131
In our CPM model, cell-substrate adhesion sites, to which actin bundles are 132 anchored, were represented as separate entities: specially labeled subcells of the lattice. 133 If two adhesion sites of two cells came into contact, we established the connection 134 between them. The connection means, that the new bond energy was applied to them. 135 This energy depended on the angle between the linked actin bundles (see Figure 4). The 136 minimum of the energy corresponded to smooth coupling between the bundles, or zero 137 angle. In this case, the bond was the most stable, but couplings with the non-zero angle 138 between the bundles had a tendency to break apart. With this new energy term that favours cytoskeletons alignment, the cardiomyocytes 140 in simulations created branching patterns. In Figure 5a a resulting simulated structure 141 of the sample with 70% fibroblasts is shown. We see that in this sample only 30% of 142 cardiomyocytes were able to build a network, and even with such a high density of 143 fibroblasts, this network was fully interconnected. Our further studies showed that such 144 interconnection was established for every sample with 30% of cardiomyocytes (n = 10) 145 of 1 cm × 1 cm size, which is close to the size of the Petri dish used in our   Figure 1a). One can see, that the number of percolation blocks per unit 165 area is similar to that of the experimental activation pattern, and the trajectories of the 166 waves share the same features. Note, that the scale of the simulated activation map is 167 slightly smaller than those of the experimental one.

168
The percolation threshold in virtual cardiac monolayers was equal to 71.5 ± 1.5% of 169 fibroblasts, meaning that 100% of the samples (n = 10) with 70% fibrosis were 170 interconnected, whereas samples with more than 73% fibrosis (n = 10) were never 171 conducting. For 72% fibrosis 20% of the samples were functional and for 71% fibrosis, 172 those were 80%. The conduction velocity dependence on the density of fibrosis is shown 173 in Figure 6b. For each simulated sample we have indicated the mean value of the 174 velocity, and the standard deviation of the velocity distribution was shown with error 175 bars. One can see, that the dependency of the velocity on the fibroblas density is similar 176 to one measured in experiments (see Figure 1b). Closer to the percolation threshold the 177 fluctuations in conduction velocities amplified due to the stochastic nature of the 178 percolation block. Moreover, variations in conduction velocities within one sample were 179 very high, because of a large number of conduction blocks.

180
The percolation threshold in our simulations was slightly lower (≈ 72%) than in 181 experiment (75%). However, both values are much higher than the predictions of the 182 models with random cells distribution (40%).  185 We have shown that in virtual tissues the spontaneous onset of reentry can occur. It 186 was observed in samples with a high level of fibrosis. Figure 7 shows such process for a 187 sample with 70% fibrosis. We see that after the first stimulus applied to the left border 188 of the sample (shown in yellow), the wave propagation was blocked at the lower part of 189 the sample, but continued propagation through the upper part. After reaching the right 190 effect is called source-sink mismatch [1]. When the source (a smaller cluster) is 205 insufficient compared to the sink (a larger inactive cell cluster), the wave propagation is 206 blocked. This effect was observed in chemical systems [22] and later in cardiac 207 monolayers [23].

208
In a sample in Figure 7 the "diode" may initiate a reentry if the sample is stimulated 209 from the left or from the top with a high frequency (4 Hz or more). 210 We have performed studies in 6 large samples with 70% of fibroblasts. All of these 211 samples had 2-5 areas of the uni-directional block and many bi-directional blocks, but 212 only one of these samples was arrhythmogenic. Thus, in addition to "diodes", some 213 extra geometrical conditions are required. The precise conditions are to be studied in 214 the future, however, they are related to the presence of the long conducting circuits in 215 the tissue. In fact, in Figure 7 we see that apart from the diode, there is also a loop, 216 which is shown with a red dashed curve in the bottom-right image. The diode is a part 217 of this loop, however, the loop is large enough to account for recovery of the bottom 218 part of the tissue after one rotation. The samples with a higher density of fibrosis were 219 even less likely to have long circuits, thus we didn't observe sustained reentry there. We 220 concluded, that reentry formation requires not only uni-directional blocks but also long 221 circuits, and the intermediate densities of fibrosis are the most arrhythmogenic ones.

222
This result was previously shown for random cell distributions, that reentry is most 223 likely to occur 5-10% below the percolation threshold [8]. The principle holds true in 224 our model, but quantitatively the densities of fibrosis are different.

225
Our model shows, that the areas of the uni-directional block can be naturally formed 226 during tissue growth. Every time one spreading cell comes in contact with the other cell 227 cluster, there is a chance that this connection will be asymmetric. It may explain the  We observed paradoxical electrical wave propagation in samples with up to 73-75% of 232 fibroblasts instead of mathematically predicted 40% for randomly distributed cells. We 233 have shown both in vitro and in silico, that electrical wave propagation was possible due 234 to the formation of the conduction pathways that rewired the whole monolayer. We have 235 proved the existence of this branching network with immunohistochemical images. We 236 have measured the conduction velocity, which decreased with the increase of the portion 237 of fibroblasts in the monolayer in a similar way in experimental and computational 238 studies. Assuming that cardiomyocytes align their cytoskeletons to fuse into cardiac 239 syncytium, the morphology of conductive pathways was successfully reproduced in 240 computer modeling as well as the electrophysiological properties of the monolayers.

241
To explain the formation of the pathways, we have considered several hypotheses.

242
First, we suggested that differential adhesion together with cell elongation may be 243 enough to explain this patterning. A similar system with autonomously elongating cells 244 was successfully used to explain vasculogenesis [24]. However, we did not impose 245 obligatory elongation on the cardiac cells, because they are not necessarily elongated 246 according to our experimental observations. Cardiomyocytes obtain their typical 247 brick-like shape with the guidance of the extracellular matrix over the course of 248 development. However, there is no evidence for any internal autonomous mechanism for 249 elongation. In experiments on the glass, cells were not only bipolar but also tripolar and 250 multipolar. In our model, this feature was reproduced with explicit introduction of the 251 actin bundles. Similar approach was previously used to describe the shapes of dendric 252 cells [25]. Cooperation between aligned actin bundles pulling in one direction shaped 253 the cell more efficiently, which resulted in clustering of these bundles and multipolar cell 254 shapes. Therefore, since the elongation was not imposed, there was no mechanism 255 forcing cardiomyocytes out of clusters to search for new connections.

256
Next, we considered various mechanisms that were previously used to explain similar 257 branching patterns in angiogenesis. The main sources of instability in those models were 258 chemotaxis and contact-inhibition [26]. The percolation in the networks formed due to 259 sharp gradients of chemoattractants was also studied previously in a continuous 260 model [27]. However, there was no evidence for directed migration of cardiac cells and 261 for type-specific contact-inhibition, similar to those that select a tip cell in a growing 262 blood vessel. Therefore, we discarded this hypothesis either.

263
This mechanism is similar to diffusion limited aggregation: a process in which 264 random-walking particles form fractal aggregates. In our model, a cell does not only 265 stick to the growing pattern but also aligns its orientation with it. The mechanism for it 266 actually exists as a part of syncytium formation, when cardiomyocytes fuse and align 267 their cytoskeletons. Once a cell aligns cytoskeleton with its neighbours, this structure 268 maturates and fixes the cell in place. In some sense, this process causes 269 contact-inhibition only between cardiomyocytes and not between cardiomyocytes and 270 fibroblasts, which results in the branching structure like one shown in angiogenesis [26]. 271 Formation of conduction pathways and complex texture of the tissue may be 272 important not only in terms of arrhythmogenicity. For example, it was shown that 273 texture of cardiac tissue at subcellular level can substantially affect the propagation of 274 external current during defibrillation [28,29]. It would also be interesting to quantify 275 the excitation patterns in terms of the number of re-entrant sources and wavefront 276 complexity [30]. Therefore, it will be interesting to perform a similar study for the 277 textures generated with our model.

278
There are several limitations of the methods used in this study. First of all, we have 279 conducted experiments with cell cultures, which are different from the real cardiac tissue. It would be interesting, yet more complicated, to study patterning of the real 3D 281 cardiac tissue. However, the mechanisms of patterning that we discovered here could be 282 universal and might be applicable to the real cardiac development as well. Second, 283 fibrosis is a more complex condition than just excessive growth of fibroblasts, which also 284 involves collagen deposition. The collagen electrically insulates the cardiac fibers from 285 one another and is also considered to increase arrhythmogenicity. In this study, we did 286 not take the extracellular matrix into account, but it would be also interesting to 287 measure its effect on the percolation threshold in the future studies. Finally, percolation 288 depends on the size of the sample and scales with this size. For random systems, the 289 scaling laws are known. In our case, we did not consider scaling and used similar size of 290 the samples in silico as in the experiment. However, it would be interesting to study 291 scaling of the percolation threshold and see how the size can affect the probability of 292 wavebreak formation. 293 We conclude, that the cardiomyocytes in fibrotic areas can form a connected network 294 and allow electrical signal propagation in monolayers containing up to 75% of 295 connective tissue. for CM-specific labelling, Alexa Fluor 488 phalloidin (Molecular Probes, USA, A12379) 319 for F-actin non-specific staining and DAPI for labelling cell DNA. Pictures were taken 320 with an inverted fluorescence microscope (Axio Imager with ApoTome optical sectioning 321 module, Zeiss). Immunofluorescent staining of the CMs was performed with the use of 322 secondary and primary antibodies according to a previously described protocol 2 . Signals 323 from each fluorescent label were recorded in a corresponding wavelength range.

324
Channels were then preudo-coloured and merged together. In this study, we have used a mathematical model that was previously developed for 346 cardiac monolayers formation [15]. It is based on the Cellular Potts Model (CPM) 347 formalism, which describes cells as the domains of the regular lattice and assigns energy 348 to the system of cells to describe their growth and motility. In our previous work [15], we 349 have selected the main features of the cardiac cells (such as area, number of protrusions, 350 etc.) and parametrised our model to reproduce the cell shapes. In this study, we 351 extended the model with a new energy term responsible for syncytium formation.

352
The evolution of the cardiac monolayer is described by the Hamiltonian: where H adhesive + H elastic is the basic CPM model and H protr is the term describing the 354 protrusion dynamics of the cardiac cells, which produces a characteristic polygonal  The spreading of virtual cells is shown in S5 Video. One of the extra features of our 365 approach, is that elongation of the cells is not imposed. Cardiac cells in a real heart are 366 indeed elongated, as they are guided by the extracellular matrix, whereas in a petri dish 367 cells can spread in any direction with no preference. However, this does not result in a 368 circular shape due to the limited number of attachment sites where spreading occures. 369 Therefore, instead of explicit declaration of elongation, we suggested that virtual cells 370 have only a small ( 10) number of well-developed mature attachment sites. Moreover, 371 these sites in a model tend to cluster stontaneously (see S5 Video), as the actin strands 372 pull the cell more efficiently acting together, rather than separately. Therefore, these 373 clustered attachment sites turn cell shapes into bipolar, tripolar or multipolar states.

374
The bipolar state is the most stable one for isolated cells, however tripolar cells are also 375 relatively common. This correctly represents the population of cells observed in 376 experiments.

377
H protr is an energy term, that decreases with the distance from the centre of mass of 378 the cell, and which is applied only to these labeled attachment sites. As a result, the 379 attachment sites spread out and reproduce a characteristic polygonal shape of the 380 cardiac cells. Here, in this model, we omit the details of the spreading process, which 381 involves polymerization and depolymerisation of the actin, attachment/detachment, etc., 382 but we mimic the overall dynamics of the protrusions. Also, we assume, that for every 383 attachment site a corresponding actin bundle exists, which stretches from the 384 attachment site towards the proximity of the nuclei.

385
If two attachment sites in the neighbouring cells come into contact, the bond between these attachment sites can be formed. In the algorithm, this bond appears if one attachment site attempts to move over the other. In such case, instead of copying of the subcell the connection establishes. A new energy term H junctions applies to the subcells, which are involved in a newly established junction. This term determines the stability of the cell-to-cell junction and depends on the angle between the actin bundles, associated with the attachment sites involved (see Figure 4). It is determined as follows: where α i,j (also shown in Figure 4) is the angle between two cytoskeleton bundles of 386 two neighbouring cells which ends at points i and j which are labeled as the parts of the 387 junction. The higher is the angle, the less stable is the junction. Therefore, the 388 junctions with continuous actin bundles on both sides persist, but the kinked bundles 389 tend to lose the connection. As a result, the actin bundles of the neighbouring cells tend 390 to align and stay in the aligned state.
The parameters of the model used in simulations (see Table 1) in this paper were 392 adjusted to compensate the additional energy term H junctions . The most of them are 393 close to the parameters used in our original paper [15]. The value of the new energy 394 term was set to E bond = 5.0, which provided enough stability for the junctions to 395 maintain the branching structure but at the same time not too much stability to allow 396 cells to search for possible new connections. Addition of H junctions effectively increased 397 the adhesion between cardiomyocytes. Therefore, the differential adhesion was toned 398 down in a model to allow cardiomyocytes to migrate randomly before they maturate 399 and stick to the pattern. In this study we have changed type-specific adhesion  The virtual samples generated with the Cellular Potts Model were then used as a map 411 for electrophysiological studies. The methodology was already described in detail in our 412 previous paper [15]. The parameters of the electrophysiological model were the same as 413 in [15], except for the coupling coefficients that were slightly adjusted (D in = 0.4 cm 2 /s, 414 D L = 0.4 × 10 −2 cm 2 /s, D T = 0 cm 2 /s) to reproduce the maximal conduction velocity 415 in the control samples with low fibrosis. The same coefficients were then used for all