Mechanisms of stochastic onset and termination of atrial fibrillation studied with a cellular automaton model

Abstract

Mathematical models of cardiac electrical excitation are increasingly complex, with multiscale models seeking to represent and bridge physiological behaviours across temporal and spatial scales. The increasing complexity of these models makes it computationally expensive to both evaluate long term (more than 60 s) behaviour and determine sensitivity of model outputs to inputs. This is particularly relevant in models of atrial fibrillation (AF), where individual episodes last from seconds to days, and interepisode waiting times can be minutes to months. Potential mechanisms of transition between sinus rhythm and AF have been identified but are not well understood, and it is difficult to simulate AF for long periods of time using state-of-the-art models. In this study, we implemented a Moe-type cellular automaton on a novel, topologically equivalent surface geometry of the left atrium. We used the model to simulate stochastic initiation and spontaneous termination of AF, arising from bursts of spontaneous activation near pulmonary veins. The simplified representation of atrial electrical activity reduced computational cost, and so permitted us to investigate AF mechanisms in a probabilistic setting. We computed large numbers (approx. 10⁵) of sample paths of the model, to infer stochastic initiation and termination rates of AF episodes using different model parameters. By generating statistical distributions of model outputs, we demonstrated how to propagate uncertainties of inputs within our microscopic level model up to a macroscopic level. Lastly, we investigated spontaneous termination in the model and found a complex dependence on its past AF trajectory, the mechanism of which merits future investigation.

Keywords: atrial fibrillation; cellular automata; model; re-entry; termination.

Figure 1.

Visualization of the spherical geometry representing the left atrium, with nodes distributed regularly over the surface. Anatomical features (black) were rendered electrically inactive. LS/LIPV, left superior/left inferior pulmonary vein. RS/RIPV, right superior/right inferior pulmonary vein. MV, mitral valve. Fibrotic cells (red) were distributed randomly over a disc centred on the posterior atrial wall.

Figure 2.

Snapshots of the simulation using a Mollweide projection. The centre of the projection is at the posterior atrial wall. Fibrotic nodes (red), which cannot be initiated, are set to a constant state. (a) In sinus rhythm, SN breakthrough starts proximal to the right pulmonary veins (PVs), and cells are immediately excited from 0 to (maximal refractory period, RP) state 120, decreasing its state by 1 each time step until it reaches 0. Cells nearby are excited to 120 if the number of neighbouring cells which are excited exceeds 8, and this starts a wavefront of activation over the sphere, with slower activation through fibrotic areas and around PVs. (b) In re-entry, existing wavefronts self-perpetuate across the domain, and SN breakthrough does not initiate wavefronts of excitation. RP restitution has meant that cells excite to a lower state compared with sinus rhythm, and this also leads to a shorter wavetail.

Figure 3.

(a) Protocol for investigating AF initiation. Starting in sinus rhythm (SR), PV bursts of up to 5 s were initiated, after which a 10 s observation window with sinus pacing was simulated to probe existence of AF. (Top/bottom sample path: with/without AF. In the top sample path, the sinus breakthrough region cannot be excited by sinus pacing, because re-entrant waves keep re-exciting the region from state 0.) (b) Probability of re-entry (the proportion of simulations finishing in re-entry over 10⁵ sample paths) as a function of the PV burst duration, for the model parameters fibrosis density (FC), PV bursting rate (BR), and restitution steepness (B,K). The slope of the curves quantifies the continuous-time rate to induce AF re-entry. Discrete markers: simulation results; continuous lines: best linear fits. (c) Data replotted as a function of model parameters, with each line representing PV burst duration, which highlights the non-monotonic dependence on model parameters for BR (middle) and restitution steepness (B, K). We have joined the scatter plots in this panel for optimal visualization of the non-monotonicity in this panel.

Figure 4.

(a) The second classifier identifies the system as ‘in AF’ if the proportion of active cells is greater than 0.5 for more than 2 s. (b) The onset time of AF is a random time τ, and the cumulative distribution of τ is plotted for selected parameter sets. The data were fitted by an exponential function . Note that when we varied BR, the cumulative distribution was monotonically decreasing for any given time. This indicates that the second classifier identifies the transition rate to enter AF is a monotonic increasing function of BR, in contrast to figure 3 where AF is mostly induced at the intermediate regime of BR.

Figure 5.

A parallel analysis of figure 3 to analyse AF initiation rate for longer duration of PV bursts up to 300 s. (a) Re-entry probability as a function of PV burst duration, for different model parameters. Discrete markers: simulation results; continuous curves: best fits using equation (3.3). (b) Re-entry probability as a function of model parameters, for different PV burst durations. Re-entry probability is nonlinear for all parameters, and exhibits non-monotonic dependence on BR and restitution steepness. We have joined the scatter plots in this panel for optimal visualization of the non-monotonicity in this panel.

Figure 6.

Fourier analysis of the classifier signal in AF re-entry. Both panels used baseline parameters FC = 300, BR = 20 Hz, B = 1 and K = 40. In (a), the Fourier spectra of 200 sample paths ending with re-entry after 5 s PV bursts were overlaid. In (b), 200 sample paths with re-entry after 150 s PV bursts were plotted. We plotted 20 sample paths on the temporal domain in the insets. Above the main plots, we present typical snapshots of the visualization (movies provided in the supporting information), showing a more ‘homogeneous’ travelling wave in (a), and a more ‘fragmented’ wavefront in (b).

Figure 7.

(a) Schematic diagram of the protocol to probe spontaneous termination. We turned on PV bursts for a duration T₁ (unit: s), to obtain sample paths which initiated AF. For these, a second PV burst of duration 1 was applied, following a waiting window of duration T₂ (s). Termination probability was measured in the 10 s window after the second PV burst. (b) The termination probability across all sample paths depends on T₁; in contrast, the value of T₂ does not change the termination probability significantly. (c) We performed 50 trials for T₂ = [0.5, 1, 2, 4, 6, 8, 10], for each sample path in AF after the first PV burst. Probability of AF termination for each sample path (single line), as a function of the sample path and T₂, is shown in the heat map. Strong correlation in the horizontal direction suggests two subpopulations: stable AF which cannot be terminated, and unstable AF which can be terminated, regardless of the T₂ value. The results suggest the subpopulation of stable AF increases as T₁ increases.