Spatially heterogeneous ionic remodelling promotes wavefront instability in simulated atrial tissue

Electrical remodelling is often seen as a key determinant in the development of atrial fibrillation. From a computational electrical modelling point of view, remodelling is usually incorporated by a global (within a specific atrial region) change of model parameters, e.g., ionic conductances and calcium handling properties. Given the spatially heterogeneous nature of cardiac tissue, one may expect that remodelling develops in a non-uniform manner. Computational models sometimes reflect this when fibrotic regions are considered; however, in the case of spatially heterogeneous electrical remodelling the literature is much sparser. Here we demonstrate that a spatially heterogeneous mixture of remodelled and non-remodelled regions can cause a destabilizing effect in simulated tissue paced at a sufficiently high frequency, indicating a potential mechanism for the initiation of fibrillatory behaviour. In particular, we use an established model of human atrial action potential with parameters based on sinus rhythm and chronic atrial fibrillation, and we distribute these two possible states according to realisations of the 2-dimensional Ising model. We show that there is a window of pacing frequencies where turbulent propagation behaviour (via the break up of incoming wavefronts) is sustained, and that a long turbulent transient can be sustained when adapting to fast pacing from a lower pacing frequency. These results might help us think about the mechanisms embodied in the aphorism ‘atrial fibrillation begets atrial fibrillation’, namely, that partially remodelled tissue could emerge in an intermediate stage of disease progression, promoting the conditions for the tissue to progress to becoming globally remodelled.

The beating of the heart is maintained by a complex system of interacting components at varying scales, which interact in such a way that a regular rhythm emerges. In the complex, non-stationary environment that the heart finds itself, the maintenance of a regular heart rhythm demands that the constituent components of the organ at smaller scales be able to adapt their behaviour in order to maintain the necessary collective behaviour. This can lead to pathological changes, as can be seen in the phenomenon of atrial remodelling, which is associated with atrial fibrillation (AF) (Nattel and Harada, 2014). This phenomenon encompasses a series of changes in ion channels, calcium handling and structural properties that occur under long durations of high-frequency pacing, which may originate from protective mechanisms, but ultimately result in a landscape more conducive to AF (by the shortening of the action potential (AP) for example) (Nattel, 1999). AF itself creates conditions of repeated high frequency pacing near re-entrant circuits, hence the phrase ‘AF begets AF’ (Wijffels et al., 1995), which reflects the fact that AF reproduces the conditions necessary for its own maintenance through progressive remodelling.

In computational studies, one can mimic remodelling by changing the model parameters. At the tissue level setting, the atria are typically partitioned into regions (e.g., left atrium, right atrium, pulmonary veins, etc.) and the parameters are typically uniform within these regions. The assumption of a uniform parameter regime within a given region is sometimes broken, e.g., when fibrotic regions are considered (Clayton, 2018) or if a spatial gradient in the parameters is used (Ten Tusscher and Panfilov, 2003). When electrical remodelling is considered in a given region, typically the parameters within that region are altered uniformly, hence there is an implicit assumption that the said region is either globally ‘normal’ or globally remodelled 1. Current studies therefore miss the possibility that remodelling could develop in a spatially heterogeneous fashion at a finer scale. Due to the complex structure of the tissue, one may think it plausible that such a heterogeneous emergence of remodelling would be possible. Indeed, there is some experimental evidence that ionic remodelling does develop non-uniformly in some cases, e.g., tissue from subjects with paroxysmal AF (Ramdat Misier et al., 1992) and atrial tachycardia (Fareh et al., 1998) displays intra-regional heterogeneity in refractory period durations. While certain aspects of structural remodelling (which necessarily introduces heterogeneity) have been considered in computational studies (Alonso and Bär, 2013; Vigmond et al., 2016; Clayton, 2018), there is little work in the context of spatially heterogeneous ionic remodelling, which is the topic of this paper.

Concretely, in this manuscript we will consider the effects of pacing a simulated patch of tissue that is comprised of regions in one of two states: the normal state and the remodelled state. The difference between the two states is reflected in the model parameters; in particular we use the established Grandi model of atrial AP (Grandi et al., 2011), in which the authors have specified sets of parameters for cells in sinus rhythm and cells in chronic AF, which we take as representative of the normal and remodelled states respectively. In addition to the usual ion-channel remodelling, we also assume a smaller conductivity in the remodelled regions, in accordance with alterations in gap junctions observed in remodelled myocytes (Nattel and Harada, 2014). To distribute the heterogeneity, we use the well known Ising model from statistical physics. Traditionally a model of ferromagnetism, in our case the Ising model is best understood as a random cluster model (Fortuin and Kasteleyn, 1972) (so is bond percolation, which has been used in other studies, e.g., Vigmond et al. (2016)), where the clusters correspond to remodelled regions randomly distributed in an ambient space of normal tissue. Broadly speaking, our results demonstrate three general points: At sufficiently high pacing frequencies, even sparsely distributed remodelled regions can destabilise wavefront propagation, causing spatiotemporally chaotic propagation to ensue.

Depending on the initial condition, a long transient period of chaotic propagation may precede an eventual settling down to a periodic rhythm. This can occur when pacing frequency suddenly increases for example.

The dynamics depends upon how the heterogeneity is distributed, which we model by varying the coupling parameter of the Ising model.

The proposed importance of tissue heterogeneity in relation to AF can be traced back to at least the classic work of Moe and Abildskov (1959) in their elucidation of the multiple wavelets hypothesis. But since then the literature on computational models incorporating (non-structural) tissue heterogeneity has been rather sparse. Xie et al. (2001) studied the effects of a model of ventricular action potential (AP) with heterogeneous properties in the I K 1 current by randomly sampling the conductance G K 1 at each point in the finite differences lattice. They studied the effects on spiral waves, finding that the AP heterogeneity can induce spiral wave breakup. Ten Tusscher and Panfilov (2003) also use heterogeneity in ion channel conductances to study the effect on spiral wave dynamics, but their heterogeneity is based upon the assumption of a spatial gradient, rather than distributed in a random fashion. Elliott et al. (2022) incorporate cellular heterogeneity by calibrating a large set of parameters to match certain physiological biomarkers to data, and then randomly sample these to set the parameters of the individual cells. In contrast, our work only uses two different sets of parameters (normal and remodelled), which act as representatives of two polar opposites states, rather than variation within the norm. Gussak et al. (2019) study parasympathetic nerve remodelling in the left atrium, and include a computational study that is similar in flavour to this work. They show that in a model including spatially uniform parasympathetic and sympathetic signalling effects, then the two have a tendency to counteract one another allowing for stable propagation. However, introducing spatial heterogeneity in the strength of the parasympathetic effects (which is suggested by their experimental data) causes exacerbation of the pro-arrhythmic effects caused by the sympathetic signalling.

In the case of structural remodelling, e.g., fibrosis, some authors have considered a percolation type-process to represent the fibrotic regions (Cherry et al., 2007; Alonso and Bär, 2013; Vigmond et al., 2016), which relates to our approach in that bond percolation and the Ising model are both random cluster models. Fibrosis corresponds to a more advanced stage of remodelling though, and this type of approach with percolation has an effect more akin to the creation of fixed topological obstacles in the domain rather than cellular heterogeneity. Clayton (2018) also considers structural remodelling but uses a Gaussian random field to generate the remodelled regions, which are modelled as non-excitatory regions (but are still conductive via diffusion). The use of a Gaussian random field allows them to consider the effect of the scales of the heterogeneity, which is similar to this work with the Ising model.

So while spatial variability in ionic properties exist in some studies, it is typically limited to spatial gradients or the assumption of noise in specific ion channel conductances. In contrast, we consider the case in which two opposing states (remodelled and non-remodelled) coexist in randomly distributed regions.

By linking these microscopic cellular kinetics to mesoscopic structural arrangements and their combined effect on macro-scale tissue stability, this approach highlights the multiscale dynamics inherent to AP propagation in cardiac tissue. Characterizing how these structural and functional components interact to govern the transition to arrhythmia aligns with the broader principles of Network Physiology, emphasizing that tissue-level vulnerability emerges from complex interactions across distinct physiological sub-systems.

To setup the general problem considered in this manuscript, we introduce the following four assumptions relating to the structure of the model of atrial tissue: The tissue can be described by a 2-dimensional grid, where each region corresponds to a cell, or the average properties of a small volume of cells in close vicinity.

Each region can exist in one of two states, normal or remodelled, and the difference between these two states manifests itself in the ionic properties and the diffusivity at that particular site.

The tissue is in an equilibrium state (or a quasi-equilibrium state) so that the spatial distribution of remodelling remains constant over the time scales considered; i.e., the remodelling distribution does not evolve over the simulation.

Each site is correlated with neighbouring sites, so that remodelling at a site is more likely if remodelling has occurred at adjacent sites.

To resolve a concrete form for the spatial distribution of remodelling satisfying the above assumptions, we have chosen to use the well known Ising model from statistical physics.

The Ising model

The Ising model, which was originally formulated as a model of ferromagnetism, has found wide ranging applications across physics (Brush, 1967) and biology (Matsuda, 1981). Consider a 2-dimensional lattice L indexed by sites (i, j), where each site is to be considered as normal or remodelled. We assign to each site (i, j) a state, denoted r i j, which is given by

The Ising model describes a probability distribution over the whole configuration of sites r ≡ { r i j } (i, j) ∈ L. To describe the probability distribution, we must first assign each site an energy where K > 0 is a coupling parameter and h is an external field strength parameter. One can easily check that with the above definition, the energy at the site (i, j) is minimised when r i j has the same sign as h and r i ± 1, j = r i, j ± 1 = r i j. In fact, each neighbour of r i j in the opposite state incurs an energy of K units. We then define the energy of the total configuration of the lattice to be the sum of energies over individual sites, E (r) = ∑ i, j E i j, and the probability of the configuration r is a Boltzmann (or Gibbs) distribution (Landau et al., 1978) based on the energy, i.e., it has the form where Z is a normalising constant and β is the inverse-temperature, which for our purposes can be assumed to be 1 2. The characteristics of realisations of Equation 1 depend upon the coupling parameter K and the ‘external field’ parameter h, which in this scenario could be thought of as ‘remodelling propensity’ 3.

It can be shown that in two dimensions there is a so-called critical phase of the Ising model when K = K c ≡ log (1 + 2) / 2 and h = 0 (Gallavotti, 1999). When K < K c and h = 0, the coupling effect is small enough to yield significant spatial variation in typical realisations of the state, and in particular we will observe small clusters of remodelled/normal tissue (the typical size of which scale depends on how close K is to K c) interwoven with one another throughout the lattice (see the left of Figure 1). If one were to zoom out far enough a realisation of the Ising model when K < K c would appear to resemble uncorrelated noise 4. When K > K c and h = 0, the coupling is strong enough to cause a single state to engulf almost the entire lattice, with only some small clusters of the opposite state sparsely distributed across the lattice (see the right of Figure 1). If one were to zoom out in this case, the configuration would appear to be completely dominated by one particular state. When K = K c and h = 0, clusters of each state exist at all scales, so zooming out reveals a kind of self-similarity in the configuration, and never resembles uncorrelated noise nor the domination by a single state (see the middle of Figure 1).

In this work, we explore the effects of remodelling when it is spatially distributed according to the Ising model with h = 0 and K ∈ { 0.9 K c, K c, 1.1 K c }. These were chosen as representative structural regimes; see Figure 1. At criticality, i.e., when K = K c, we are approximating the effect of large, well defined regions of remodelling that co-exist with normal tissue. When K = 1.1 K c, we are approximating the effect of remodelling emerging within a sea of normal tissue, which has not yet grown into significant remodelled regions. When K = 0.9 K c, we are approximating the effect of remodelled regions that have emerged in a more spatially uniform fashion compared to K = K c. So for K = 1.1 K c, we are modelling something that one may expect early on in the remodelling timeline, while the other two cases can be thought of as two different ways that remodelling could appear midway through the development of remodelled tissue.

To generate realisations of the Ising model, we use the Swendsen-Wang algorithm (Swendsen and Wang, 1987). Due to finite size effects, some realisations at criticality will include a single cluster that engulfs the bulk of the lattice, and such realisations fail to capture the full richness of the scale-free effects. To prevent such realisations from being generated, when sampling configurations at criticality we condition on the events that the proportion of tissue that is remodelled falls between 30% and 70%, and that the largest cluster of remodelled tissue is at least 5% and at most 20% of the lattice. The former condition ensures that the realisation is not almost entirely remodelled, and the latter biases the realisations towards resembling configurations where remodelled tissue has emerged from normal tissue rather than the other way around. For K = 1.1 K c, we flip the configuration (swap the labels of the normal/remodelled states) if the dominant state of the realisation is the remodelled state, ensuring that the non-remodelled state is the dominant one. The configurations when K = 0.9 K c all resemble the left of Figure 1, so no additional conditioning is required.

Action potential model

We use an adaptation of the Grandi model of human atrial action potential (Grandi et al., 2011), as described by Voigt et al. (2013). To model the normal state we use the typical parameters given in Grandi et al. (2011), and to model the remodelled state we use the chronic AF parameters (cAF) given in the same work. Compared to the original Grandi model, the adaptation uses the Courtemanche model (Courtemanche et al., 1998) variant of I Na, which has faster recovery variables than the Grandi formulation, allowing it to respond to faster pacing without propagation failure occurring.

We do not fix the sodium concentrations 5 in the model but do fix the intracellular potassium concentration, which slowly drifts towards 0 if it is not fixed. The significance of the sodium concentration variables is that they are slowly evolving, approaching equilibrium over hundreds of beats (while the other variables typically adjust over only a few beats), and have a strong effect on the action potential duration (APD). The sodium concentration tends to increase upon faster pacing, which causes the APD to shorten and therefore allows for the faster pacing to be maintained. Abrupt increases in pacing frequency can therefore result in a long period of non-stationary behaviour as the sodium concentration adapts. See the Supplementary Material for a more detailed explanation of the dynamics of the sodium concentration in this model, and see Figure 5C in Voigt et al. (2013) for the APD rate dependence curves for both normal and remodelled parameters. At a given pacing period T, we numerically find the limit cycle for the single cell model for both sets of parameters, and set the state of the normal/remodelled regions accordingly. To find the limit cycle, we use Newton’s method on the Poincaré map of the single cell model for the two sets of parameters, which is explained in more detail in the Supplementary Material.

For normal tissue, we have assumed an effective diffusivity of D normal = 0.1 mm 2 /ms and an anisotropy ratio of 1. In remodelled tissue, we have assumed that some level of structural remodelling and/or gap junction downregulation has occurred and set D remodelled = 0.05 mm 2 /ms with an anisotropy ratio of 2, which give a similar reduction in diffusivity/increase in anisotropy ratio in Varela et al. (2016). Note that real tissue is not isotropic, however some studies have found that the effective anisotropy ratio is close to 1 in healthy atrial tissue (Houben et al., 2004). Moreover, patients with long term persistent AF have been measured to exhibit a greater degree of anisotropy than those with acute AF (van Schie et al., 2024), suggesting that a larger anisotropy ratio in the remodelled tissue serves as a reasonable phenomenological baseline. With entirely non-remodelled tissue, these parameters yield a conduction velocity of approximately 0.5 m/s (pacing at T = 1000 ms), and when heterogeneity is sampled from the critical Ising model (K = K c), the conduction velocity only slightly decreases to approximately 0.45 m/s. These conduction velocities seem reasonable in the current context given that they are within a physiological range (they are on the lower end of what one might observe in induced AF and within the normal range for persistent AF (van Schie et al., 2024)). Details of the numerical methods used for simulation are outlined in the Supplementary Material.

Simulation details and metrics

For all experiments we use a 400 × 400 grid with a spatial step of Δ x = 0.15 mm, and a time step of Δ t = 5 × 1 0 − 3 ms. We used the monodomain model to model propagation, which was implemented numerically using a finite volume method (see Supplementary Material for details). To stimulate wave propagation, a stimulus current was periodically applied at the bottom edge of the grid so that wavefronts travel upwards (in the absence of condition blocks at least). Unless otherwise stated, for each experiment we run the simulation for 2,000 beats. We use 10 instances of heterogeneity when considering K = K c for each set of pacing conditions (with the exception of the one reported in Section 3.3), as the critical Ising model displays significant variability at the macroscopic level, as shown in Figure 1. For K ≠ K c, we only use one realisation since the distribution of remodelled regions is essentially isotropic away from the critical point. We record the action potential V m and intracellular sodium concentration [ N a ] i at the start of each beat and halfway between two beats, i.e., if the pacing period is T ms, we collect the two quantitites every T / 2 ms. We also record the spatial mean of V m every ms, which works well as a one dimensional summary of the global behaviour of propagation.

At low pacing frequencies, we get a stable, periodic propagation regime, while at very high pacing frequencies, we get a period doubling effect 6, as the stimulus arrives before the tissue has fully repolarised. It can be difficult to assess periodic behaviour of the system visually due to it is high dimensionality, so we must use an appropriate metric to assess this. Note that for a system in discrete time, say y t, k − periodic behaviour implies that | y t + k − y t | = 0 for all t, and if we perturb the state slightly (and the k − periodic behaviour is stable), we would expect that | y t + k − y t | ∼ e − λ t for some λ > 0. We can construct a discrete time system from the data collected by taking V m (t, x) at the time of each stimulus, i.e., y n = V m (n T, x). We then use the metric given by: where d n (k) denotes the distance metric for the n th cycle assuming period k, and | D | denotes the volume of the domain. Note that for k = 1, this reduces to the spatial average of the difference between V m at successive beats, which will decay to 0 exponentially quickly if a stable periodic limit cycle is approached. For k = 2 there are two beats in a cycle, so the measure takes the distance between the first beats of two successive cycles and the distance between the second beats of the two successive cycles and averages them. For larger values of k, the interpretation is similar: there are k beats within a cycle, so we average over the distances obtained by choosing each of the k possible beats as representative of the cycle. For experiments with K = K c, the values of d n (k) that we report correspond to the average over the 10 instances of heterogeneity unless otherwise stated.

Further details with regards to numerical methods are given in the Supplementary Material.

General behaviour

First consider Figure 2, which displays example snapshots of AP for a particular instance of K = K c and cycle lengths T = 300 ms, T = 290 ms and T = 280 ms. This broadly displays the types of qualitative behaviour that we see throughout the experiments, namely, regular periodic propagation (for T = 300 ms, panel A), irregular but still relatively laminar propagation (for T = 290 ms, panel B), and irregular, turbulent behaviour (T = 280 ms, panel C). Also indicated are some regions of conduction block that emerge in the T = 290 ms experiment shown. These regions, which are not fixed in space, regenerate themselves by destabilizing the shape of the wavefront, leading to unpredictability in the AP at the local level (but globally the behaviour is quite coherent). In Figure 2B (T = 290 ms), we see that when the wavefront encounters the block a break occurs, but it is sufficiently small as to allow the wavefront to subsequently re-merge into a single coherent front. For the T = 280 ms experiment shown (Figure 2C), it is clear that the conduction blocks occupy much larger regions with a less regular shape, thus creating pathways for turbulent re-entrant circuits to emerge.

Effect of pacing period and coupling

All instances of experiments with T = 300 ms show the emergence of a period 1 limit cycle. This is made evident in the upper plot in Figure 3A which shows the time series of d n (1) values defined in Equation 2. Note that because periodic behaviour is obtained for all experiments with T = 300 ms, we have only run the simulation for 1,000 beats rather than 2,000. An example of the propagation can be seen in Figure 2A. We can see that after an initial transient period (corresponding to the tissue adapting to pacing) the lag 1 distances decay exponentially, which is what we would expect for a system approaching a stable attractor. Moreover, the lag 1 distance decays sharply from the outset, showing that there is a lack of chaotic transient behaviour.

The behaviour when T = 290 ms can be inferred from Figure 3B, which shows the lag 1 and lag 2 distances d n (1) and d n (2), and Figure 4, which shows the envelope of the spatial average of V m. In all cases, the lag 1 distances (Figure 3B) are on the order of 1–10 mV (except for a brief period near the beginning of the simulation when K = 1.1 K c). On the other hand, we can see that for all values of K, there is a period where the lag 2 distances d n (2) are small (on the order of 0.1 mV), indicating period 2 like behaviour in all experiments. Qualitatively, this corresponds to laminar wavefronts with a slightly different shape on alternating beats. When K = 0.9 K c, one can see from Figure 3B that this behaviour persists for the entire 2000 beats, however the characteristics of the dynamics are not the same as the exponential decay seen in d n (1) and d n (2) when T = 300 ms: a slowing in the rate of decay is observed along with small oscillations. The same type of oscillatory behaviour is observed when K = K c over the first 250 beats, before irregular behaviour sets in. These small scale irregular oscillations are manifested in the spatial average of V m as well, as shown in Figure 4. This suggests that the apparent period 2 behaviour for K = 0.9 K c may be a transient phenomenon, which would not be sustained if we were to run the simulation for a longer period of time. After the transient period 2 behaviour ceases in the K = K c example, we observe an intermediate state between laminarity and turbulence, characterised by small wave breaks but still relatively steady propagation as in Figure 2B, while the dynamics of propagation for K = 1.1 K c consists of larger wave breaks and irregular patterns of propagation. One can infer that there is a qualitative difference in propagation between K = 1.1 K c and K = K c by the fact that there is roughly an order of magnitude difference between the respective lag 1 and 2 distances in Figure 3B, indicating significantly greater spatial variability in the AP at each beat for K = 1.1 K c. This is also clearly manifested in Figure 4, which shows large, aperiodic fluctuations in the amplitude of the spatial mean of V m for K = 1.1 K c, but only small fluctuations for K = K c. We have included additional AP heatmaps for K = K c and K = 1.1 K c in the Supplementary Figures to visually demonstrate the difference between the propagation behaviours.

As a point of comparison, a completely homogeneous model with non-remodelled parameters transitions to period 2 behaviour where a uniform conduction block prevents propagation every other beat. Initially, all beats propagate due to an initially high sodium concentration, which slowly drops across the tissue causing a slow lengthening of the APD. This causes the tissue to exhibit behaviour in which propagation fails in a quasi-periodic manner, which drives further depletion of sodium (via failure to activate I N a) and eventual uniform conduction block; see Figure 5. From Figure 5A, it is evident that episodes of propagation failure occur intermittently for K = 1.1 K c (shown by the drops in d n (2)) from approximately cycle 600 onward, but they are far less regular in the homogeneous model before the transition to period 2 behaviour. One can also see a difference in the dynamics of the average intracellular sodium concentration in Figure 5B, with the homogeneous model experiencing a linear decline up to the transition (after which the decline accelerates), while for K = 1.1 K c, the sodium concentration is slower to decline and appears to plateau, suggesting that the irregular behaviour may persist indefinitely.

For T = 280 ms, we observe turbulent, aperiodic behaviour in all cases, with the exception of when K = 1.1 K c, where we observe a transient period of turbulent behaviour followed by a doubling of the period, corresponding to the stimulus failing to generate a propagating wavefront every other beat. This can be inferred from Figure 6A, where it is evident that the lag 2 distances decay exponentially for K = 1.1 K c (but remain large and unpredictable when K < 1.1 K c). We see similar behaviour when T = 270 ms (shown in Figure 6B), only with a shorter transient period of turbulent behaviour when K = 1.1 K c. When T = 260 ms, all experiments experienced period doubling within 1,000 beats (Figure 6C). Once period doubling sets in, the effective pacing period becomes 2 T, e.g., 520 ms for T = 260 ms, and as we have seen, stable propagation can be achieved at periods as low as T = 300 ms, so the alternation between successful and failed propagation ultimately induces stable behaviour. Figure 6D shows the lag 8 distance d n (8) for one particular experiment with K = K c and T = 270 ms (compared to the average over the 10 instances when K = K c). We show this because the experiment in question demonstrated an interesting effect, where an apparent limit cycle with period 8 was approached on two separate occasions, and not only this but it appeared to be relatively stable. However, it ultimately proved to be unstable as the sudden return to large values of d n (8) in Figure 6D shows. This demonstrates a type of intermittent behaviour in the model; a visualization of the period 8 behaviour is shown in the Supplementary Figures.

Notice that for the K = 1.1 K c experiments in Figure 6A (T = 280 ms) and Figure 6B (T = 270 ms), before approaching the period 2 limit cycle, the system has several periods where the lag 2 distances drop significantly, but subsequently rise to large values rather quickly. These correspond to short bouts of periodic behaviour that quickly lose stability, and suggest a slow time scale in the system that controls the stability of the emergent limit cycle. We visualise this in Figure 7 by showing the evolution of the slowly evolving average N a i (partitioned into remodelled and non-remodelled areas). Figure 7A shows that the values of d n (2) for T = 280 ms appear to be correlated with the declining average N a i, and when the period doubling sets in, the rate of decay of N a i undergoes a sudden increase (until it eventually plateaus towards an equilibrium value). A similar analysis of the case shown in Figure 6D (T = 270 ms, K = K c and near period 8 behaviour) is shown in Figure 7B. In this case, N a i in the remodelled areas remains roughly steady over the first 900 beats or so, whilst it drops in the non-remodelled areas. But when the near periodic behaviour emerges, a sufficiently stable propagation regime is sustained as to cause average N a i to rise globally, which presumably lays the foundations for this cycle to become unstable once again by shortening the APD (and hence the length of the wave back) allowing propagation to break through in regions where it was previously blocked in the period 8 cycle.

Effect of initial condition

Here we consider the effect of an abrupt change in pacing period on the dynamics. In the experiments previously discussed, we have assumed that the tissue has already approximately adjusted to the pacing conditions by exploiting the limit cycle behaviour of the single cell model. In reality however, we might expect increases in pacing frequency to occur suddenly due to some external trigger, in which case the tissue would have to go through a transient period where it adapts to the new pacing frequency. We approximate this behaviour by initializing the state based on the single cell model paced every 600 m, and commence pacing the tissue with T = 300 ms.

Figure 8 displays the essence of the behaviour for the three couplings considered. Figures 8A,B show the envelope of the spatial average of V m over the entire simulations, the former showing K = 0.9 K c and K = K c, and the latter showing K = 1.1 K c. We have separated the latter instance as we have run this for a total of 5,000 beats as opposed to 2,000 beats. The reason for the extension of run time will become evident below. We also provide the plots of the distances d n (1) and d n (2) for each of these experiments in the Supplementary Figures. For K = 0.9 K c and K = K c, we can see that turbulent behaviour commences from the outset (evidenced by the large fluctuations in amplitude displayed in Figure 8A), but gradually subsides and is replaced by the periodic rhythm observed in Figure 2A. This shows that transient episodes of chaotic behaviour are possible even when a stable limit cycle exists, and can come about due to a sudden bout of fast pacing. The duration of the turbulent period is longer for K = K c than it is for K = 0.9 K c, and moreover the adaptation to periodic behaviour is more gradual.

For K = 1.1 K c the behaviour is rather different. In the first 300 s of the simulation there is propagation failure every other beat (evidenced by the constant envelope with large amplitude seen in Figure 8B). After the 300 s mark however, there is a sudden transition to chaotic behaviour, which then persists for 1,200 s (4,000 beats). This suggests the opposite of what we described above for K < 1.1 K c: initially the system is able to maintain a periodic rhythm at the original pacing period (T = 600 ms) and does not allow for propagation every second beat. But after repeatedly pacing the tissue for a sufficiently long duration, portions of the wavefront begin to break through the blocked region and turbulent propagation ensues, and evidently does not stop.

Figure 8C shows the dynamics of the average sodium concentration in the remodelled and non-remodelled areas. At a pacing period of T = 600 ms, the equilibrium sodium concentration is significantly lower than when T = 300 ms, particularly in the non-remodelled regions, which induces longer wave backs, explaining the initial turbulence for K = 0.9 K c and K = K c, and the initial propagation failure for K = 1.1 K c. We can also see that the rate of increase of sodium in the non-remodelled regions is significantly greater for K = 0.9 K c than it is for K = K c (though this difference is less pronounced for the remodelled areas), and appears to be approaching a plateau by beat number 2000, while it is still steadily increasing for K = K c. This would explain why the K = K c example takes much longer to settle down to the limit cycle compared to K = 0.9 K c. Presumably, the more equidistributed remodelling when K = 0.9 K c has the effect of shortening the APD in non-remodelled areas via diffusive effects, allowing depolarization to occur more uniformly over time (and hence increasing the average sodium concentration), whereas for K = K c, this would tend to be localised around large remodelled clusters. For K = 1.1 K c, we see a very slow increase in the sodium concentration over the first 1,000 beats in the non-remodelled regions, and a decrease in the remodelled regions 7, which experience a sudden increase in growth rate when turbulence sets in after 1,000 beats. This is not surprising: the slow increase in sodium concentration in the non-remodelled regions over the first 1,000 beats reflects an increase in sodium concentration in the immediate proximity of the edge where pacing is applied (but the rest of the domain experiences little movement in sodium concentration due to propagation failure), shortening the APD in this region and the surrounding regions due to diffusive effects. Evidently, the sodium concentration across this region increases up to a threshold upon which the APD in this region becomes just short enough to allow wavefronts to break through certain regions that were previously blocked by refractory tissue. The reason we extended the K = 1.1 K c experiment was in anticipation of an eventual transition to periodicity like the other two examples. However, the sodium concentration appears to have approached an equilibrium by beat number 5,000, suggesting that such a transition may not occur at all. Moreover, in the aforementioned region along the paced edge, the sodium concentration actually begins to slowly decrease after approximately 3,500 beats, which can be seen in the Supplementary Figures. This suggests that periodic behaviour may not be achieved at all in this example, and that the turbulent behaviour may be sustained indefinitely.

In this work we investigated how the spatial organisation of electrical remodelling influences wave propagation in atrial tissue. While it is well established that heterogeneity can promote conduction block and arrhythmogenesis, remodelling is often represented as a property of entire regions rather than as a spatially emergent process occurring at mesoscopic scales. By considering remodelled tissue as clusters embedded within a sea of non-remodelled tissue, we were able to examine how the degree of spatial organisation influences the resulting dynamics. A key observation is that the spatial organisation of remodelling can have an effect on the propagation dynamics, and in particular it often manifests through the length of transient episodes of turbulent behaviour (as is made evident in, e.g., Figures 6, 8). Our analysis also revealed the importance of ultra-slow sodium-mediated adaptation (the dynamics of which are also dependent on the distribution of remodelling) in maintaining the different dynamical regimes observed. These results suggest that the spatial structure of remodelling may be as important as its overall extent when determining tissue-level dynamics.

The effects of diffusivity

We have seen that if carefully initialised, all experiments exhibited regular periodic propagation with pacing period T = 300 ms, but developed more complex behaviour at all the faster pacing frequencies. To explain why these different behaviours occurred, we first note some general dynamical features that we would expect to see based on the structure of the model. One thing that must be considered when comparing the results is that we halved the diffusivity in the remodelled areas, so the effective diffusivity for K = 1.1 K c is larger than for the other two cases, which has implications for the behaviour observed in the corresponding experiments. If we were to consider the model in the case of a (homogeneous) one-dimensional cable, then simple scaling relations imply that varying the diffusivity D will cause the conduction velocity and wavelength to scale like D. In the case of an sufficiently long cable, if periodic pacing results in stable behaviour, then the induced wave travels fast enough such that the refractory wave back drops to resting potential by the time the next stimulus arrives, allowing another wave to propagate. Increasing D increases the wavelength and hence the mass of refractory tissue, but the increased wave speed will compensate for this and the tissue will still drop to resting potential at the next stimulus. However, when D is large enough such that the wavelength becomes comparable with the domain size, then boundary effects become noticeable, and in particular, once the wavefront approaches the boundary, the forward movement of the wave stops and hence the rate at which the wave back drops to resting potential decreases, resulting in an increased likelihood of propagation failure. Hence for K = 1.1 K c (where only approximately 5% of the domain is remodelled), it is not surprising that the system experiences regular propagation failure for T ≤ 280 ms, whereas for the other two values of K (where 30%–50% of the domain is remodelled), it only occurs at T = 260 ms.

In addition to the effects of the altered diffusivity, the altered parameters reflecting ionic remodelling exert an effect on tissue level refractoriness via the cellular level effects on APD. In the single cell model, remodelling causes the APD to become relatively insensitive to changing pacing frequencies, and in particular it is much shorter than the APD of the non-remodelled model at low pacing-frequencies (Voigt et al., 2013). This then gives us an indication as to why in Section 3.3, K = 1.1 K c undergoes a long period of propagation failure, while the other two values of K do not: Since we initialise according to T = 600 ms, the non-remodelled tissue exhibits a significantly longer APD than it would if we initialised according to T = 300 ms, however for the remodelled tissue the increase in APD is much less pronounced, which allows propagation to occur in remodelled regions without a long period of adaptation. As T approaches 300 ms, the two APDs (measured as the APD90) become quite close (Voigt et al., 2013). However with the remodelled parameters, the AP shape is much more ‘triangular’ than in the non-remodelled case, i.e., in the remodelled cells the AP drops to resting potential in a linearly fashion, while the non-remodelled cells exhibit a slower plateau towards resting potential from the point it reaches its APD90 (see Supplementary Material to see how the shapes differ). If T = 290 ms, for example, with remodelled parameters the AP is about 5 mV lower than for the non-remodelled parameters at the time of stimulation. This is significant, because when diffusion is involved, the effect will be that on average, the AP in non-remodelled regions is lowered slightly (and conversely for remodelled tissue), accelerating repolarization. The distribution of remodelled tissue is significant in this sense too: the diffusive effects of a small cluster of remodelled tissue will be dominated by the influence of a large surrounding region of non-remodelled tissue, however a large cluster of remodelled tissue allows the characteristics of the remodelled AP to persist in this region, which can then exert an effect on non-remodelled tissue near the boundary of the cluster. It is in this sense that we may be able to explain the differences in behaviour between K = K c and K = 0.9 K c: the homogeneous distribution of remodelling for K = 0.9 K c creates a stronger synchronization effect between regions (although small scale spatial fluctuations still exist), while when K = K c, we get larger, more well defined regions with differing AP characteristics.

It is clear from the above discussion that some of the particularities of the behaviour observed may be peculiar to the Grandi model and/or the value of the diffusivity used. It would be conceptually straightforward to repeat the experiments performed here with a different ionic model such as the Courtemanche model (Courtemanche et al., 1998), or to test the results of varying the diffusivity/domain size. The same can be said of the anisotropy, which we have not investigated the effect of in any detail here. We leave a systematic exploration of the effects of these factors as future work.

The effects of increasing pacing frequency

When T = 290 ms we saw that experiments with K = K c and K = 1.1 K c differed qualitatively in how they deviated from regular periodic behaviour, with the latter displaying turbulent propagation and the former displaying chaotic behaviour, but not what we would think of as turbulence (e.g., large wave breaks, re-entrant circuits etc.,). For K = K c, we saw what mediates this behaviour in Figure 2B, namely, the development of small, irregular conduction blocks corresponding to regions of incomplete repolarization. When K = 1.1 K c, the emergence of aperiodic behaviour happens to occur in a similar fashion, however the faster wave speed (see Section 4.1) causes newly formed fronts to encounter the blocks faster, causing them to grow larger and rather quickly provoke more significant wave breaks than what is seen in Figure 2B. What is interesting is that the remodelled regions appear to have a dual effect on the dynamics: on the one hand sparsely distributed remodelled regions are sufficient to destabilise the dynamics enough to promote an environment of constantly reproducing conduction blocks, driving turbulent propagation. On the other hand, when remodelled regions are sufficiently large, the opposite effect is produced, i.e., the remodelled regions create a stabilizing effect by shortening the average wavelength and conduction velocity, which helps to control the size of the conduction blocks by reducing the mass of refractory tissue (see the discussion in Section 4.1).

A possible dynamical systems interpretation for the turbulent behaviour and chaotic non-turbulent behaviour comes from the notion of an interior crisis (Grebogi et al., 1983). This occurs when an unstable limit cycle collides with a chaotic attractor, causing a sudden widening of the attractor. We suggest this because on the one hand, if we use a homogeneous, non-remodelled model with T = 290 ms, there is a stable limit cycle in which propagation fails every other beat after a transient period of quasi-periodic propagation failure 8, and we expect that a slight perturbation in the parameters 9 only slightly perturbs the limit cycle, however it need not be stable after the perturbation. On the other hand, when K = K c we see unpredictable but globally coherent behaviour, which could plausibly correspond to motion along a low dimensional chaotic attractor. One could then imagine that interpolating between the two regimes would lead the unstable limit cycle and chaotic attractor to collide, resulting in dynamics that exhibit large fluctuations between beats. The result would then be intermittent episodes of motion that follows the unstable limit cycle, or in the terms of the system at hand, intermittent episodes of propagation failure. Indeed, this type of effect can be observed in some of the T = 290 ms experiments with K = K c (though not all of them, it appears to be dependent on the Ising model realisation), where the dominant dynamics are like those in Figure 2B, but with occasional episodes of re-entrant like behaviour that causes the subsequent beat to be blocked. We show an example of this in the Supplementary Figures.

For T = 280 ms we see that K = K c (and K = 0.9 K c) both undergo a transition to turbulent behaviour, so the aforementioned stabilizing effect of remodelling only manifests itself over a small window of pacing periods. The mechanism for this is a straightforward extension of what we described above, namely, that wave fronts encounter conduction blocks faster, causing a reverberation effect to take place and eventually spark turbulent behaviour. The only difference in this case is that a shorter interval between stimuli provokes the behaviour, rather than an increased conduction velocity. When K = 1.1 K c, we see a long transient period of turbulent behaviour, which proves to be unstable and eventually gives way to period doubling, where the stimulus fails to provoke a wavefront every other beat. When T = 270 ms, we largely see the same behaviour as when T = 280 ms, the main difference of note being that when K = 1.1 K c, the transient period of turbulent behaviour is shorter. And finally when T = 260 ms, all experiments undergo a transition to propagation failure, showing that in all cases there is a threshold pacing that cannot be maintained. One could interpret this as a kind of ‘protective’ effect, since the threshold (in terms of T) for repeated propagation failure is higher when the proportion of remodelled tissue is lower, which one would logically expect since remodelled tissue is thought to stabilise high frequency re-entrant circuits (Nattel et al., 2008). Nevertheless, the transient periods for sparse remodelling (K = 1.1 K c) can last for long durations (For T = 280 ms it lasts for approximately 1,400 beats, for example), so one could imagine that a growing fraction of remodelled tissue results in the lowering of this pacing threshold, and the prolongation of the average length of the transient until the turbulent behaviour becomes stable. In the experimental work of Wijffels et al. (1995), there is a similar transition in the duration of AF episodes relative to the amount of time that AF was induced beforehand using a pacemaker (which one would expect to be correlated with the degree of remodelling present), where the duration of triggered AF episodes rise sharply from a few seconds up to weeks as the duration connected to the AF pacemaker rises.

We note that the aforementioned notion of an interior crisis could plausibly apply to explain the behaviour of the transition to turbulent behaviour for K < 1.1 K c as well, see, for example, the intermittent drops in d n (2) for T = 270 ms, K = 0.9 K c in Figure 6B. However, we want to note that the exceptional case that we described when period 8 behaviour emerges intermittently (T = 270 ms, K = K c) might be more resemblant of Pommeau-Manneville intermittency (Pomeau and Manneville, 1980), where trajectories slowly diverge from an unstable limit cycle, but ultimately return to it is vicinity after a chaotic passage 10. As for the transient episodes of turbulent behaviour when K = 1.1 K c, inspection of the lengths of the transients suggests that a so-called boundary crisis (Grebogi et al., 1983) could be a plausible explanation. When a boundary crisis occurs, a chaotic attractor is destroyed upon collision with the basin of attraction of a stable limit cycle, however there are residual effect on the dynamics in the sense that trajectories can become ‘tangled up’ for a lengthy period near the chaotic set. Moreover, if the boundary crisis occurs at some critical parameter T c, then the average length of the transient chaotic episode, say τ, usually follows a power law: τ ∼ | T − T c | − γ for some exponent γ (Grebogi et al., 1987). This would explain why the duration of the transient episode is significantly shorter for T = 270 ms than for T = 280 ms (and even shorter for T = 260 ms), although we do not have enough data to estimate the value of T c nor to test whether there is truly a power law relationship. However, this would be a conceptually straightforward addition to this study that we leave as future work. We also note that similar scaling relationships can be derived for times between bursts in crisis induced intermittency (and Pommeau-Manneville intermittency), but again we do not have sufficient data to assess the validity of these relations here and leave it as potential future work.

Adaptation to fast pacing

In Section 3.3 we explored how the model would behave in response to a sudden increase in pacing period by setting the initial condition to approximate the state based on a pacing period of T = 600 ms, and then proceeded to pace at double the frequency, i.e., with T = 300 ms. When K = 1.1 K c, the model was unable to accommodate the faster pacing for approximately 1,000 beats (or approximately 300 s), i.e., propagation continued with a frequency of 1/600 ms, before giving way to a persistent epoch of turbulent behaviour. It does not appear to approach the limit cycle that we know exists from Figure 3A, while the experiments for the other values of K eventually approach their respective limit cycles, i.e., they adapt to the pacing frequency. The K = 1.1 K c experiment produces intermittency: Notice that in Figure 8B, towards the end of the simulation (from approximately 1,000 s onward) there are periods where the upper envelope gravitates around −30 mV and the lower envelope slightly above −60 mV (shown by the dashed black lines). These correspond to periods of relative order in which propagation succeeds at every beat without any major wave breaks or re-entry circuits forming, signalling motion that approaches the limit cycle, but the system always returns to chaotic behaviour. This either occurs by the upper and lower envelope moving closer together (e.g., right at the end of the simulation), corresponding to re-entrant activity 11, or by a large drop in the lower envelope, corresponding to propagation failure and hence a domain wide descent towards the resting potential.

When K < 1.1 K c there is a transient episode of turbulent behaviour, which eventually gives way to a periodic rhythm. Compared to K = 1.1 K c, where remodelling is sparse, here the remodelled regions play a key role in sustaining propagation as they have a lower resting potential. However, the low sodium concentration present in the initial state causes a lengthened APD in the non-remodelled regions, and therefore although a wavefront generally forms upon stimulation of the domain, it faces irregular conduction blocks in the non-remodelled regions. But after a sufficient duration of pacing, the sodium concentration rises to a level as to allow for coherent propagation, and the system approaches the limit cycle that we observed in the earlier experiments. The remodelled and non-remodelled tissue therefore act in opposition to one-another in this scenario: the non-remodelled tissue has a tendency to adapt slowly and block propagation instigated at high pacing frequencies, whereas the remodelled tissue displays a short and inflexible APD, permitting propagation under a sudden increase in pacing frequency. The result of this opposition is a transient period of spatially disorganised propagation while the non-remodelled tissue adapts.

One theme throughout this work is that the experiments require simulating hundreds or thousands of beats to determine whether the behaviour is stable or not. In one sense, that is a strength in that it could connect to phenomena such as intermittency. In another it is a weakness in that it is highly computationally demanding. This would make our study difficult to adapt to experiments involving anatomically realistic atrial domains. We would argue that the primary reason that these long transients arise is the fact that the APD adapts on an ultra-slow time scale (mediated by sodium concentration). Hence, if one wished to perform similar experiments to the ones in this work on more complex domains, the most appropriate approach would be to use a significantly simpler ionic model (or a low dimensional phenomenological model like the Fenton Karma model (Fenton and Karma, 1998)) with the ultra-slow dynamics abstracted to a discrete per beat map. It is plausible that techniques that exploit separation of time scales (such as averaging and homogenisation (Pavliotis and Stuart, 2008)) could be used to derive the sodium dynamics (or equivalently a parameter that phenomenologically represents sodium) in terms of a discrete map, which could be coupled with the propagation model and updated once per beat. Moreover, due to the strong separation of time scales, it is plausible that the sodium dynamics could be accelerated without qualitatively changing the dynamics, which would allow the computational load to be reduced by an order of magnitude. The details of such an approach would require rigorous investigation that we leave as future work.

The role of heterogeneity

In this work we have explicitly assumed that atrial tissue is heterogeneous in that certain regions exhibit remodelled characteristics and others do not. It is notable however that real tissue can be assumed to have heterogeneities at multiple scales without necessarily assuming that it arises from a ‘pathological’ state (i.e., remodelling in this case) (Bueno-Orovio et al., 2014; Captur et al., 2017). The monodomain equations with homogeneous conduction properties can be derived using homogenisation (Colli Franzone, 2014), which essentially averages over small scales to yield a macroscopic effective conductivity. However, as we have seen with the K = 1.1 K c experiments, even a small amount of (non-uniform) contamination at microscopic scales can yield effects at the macroscopic scale. As we have seen though, the effects at the macroscopic scale are only revealed when pacing frequency is sufficiently high, and can be exacerbated by slow adaptation to an increasing pacing frequency. That is to say, the homogenised model may work perfectly well in representing the behaviour under typical conditions, but become a poor representation under specific conditions such as high frequency pacing.

One possible direction for future work would be to alter the model structure used here by relaxing the assumption of a strict ‘normal’ and ‘remodelled’ binary. It is plausible that regions may exist in different stages of remodelling, e.g., paroxysmal AF is usually associated with remodelling of calcium handling, but not with the more extensive ion-channel remodelling observed in chronic AF (Voigt et al., 2014). Another possible direction for future work would be to explore the use of a more realistic data generation process for the remodelling configurations. Here we have used the Ising model, which is an equilibrium model, but in reality an active model (e.g., a cellular Potts modelling framework like in Kudryashova et al. (2019)) would be more appropriate, since in reality the tissue itself is not a static medium. Moreover, the Ising model may not be sufficiently flexible to capture the higher-order statistical structure observed in real biological data, e.g., in Kudryashova et al. (2019), suggesting that a more complex data generation process may be required if fitting to experimental data were pursued.

The key notion conceptually in this study is not necessarily the specificity of remodelling, but rather the opposition between cells with different AP characteristics, which theoretically could arise from other factors, e.g., distribution of parasympathetic fibres (Gussak et al., 2019). One could compare this to the technique used for modelling fibrotic tissue, in which it is often the case that tissue remains conductive but not excitable (e.g., Clayton, 2018), i.e., where all ionic currents are set to 0. This would be the extreme case, where the opposition is between excitable and non-excitable regions, in which the latter can essentially operate as topological obstructions, creating fixed conduction blocks. We should note that fibrotic tissue is not always modelled as non-excitable, as recent studies have found that myofibroblasts can exhibit a fast sodium current (Koivumäki et al., 2014). Introducing these into models can provide a continuum between a regime where fibrotic regions support propagation and one where they behave more like topological obstacles (Sánchez et al., 2019). Our case is most similar to the latter scenario, where there are no fixed conduction blocks, but conduction blocks can spontaneously arise due to lack of adequate synchronization in repolarization. In this sense, what we have produced here is more similar in spirit to the classic model developed by Moe and Abildskov (1959), in which non-uniform repolarization properties is explicitly the cause of fibrillatory behaviour.