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

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. 105) 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.


Introduction
Mathematical and computational models have become an increasingly popular tool for investigating biological and physiological systems. The quantitative capabilities of models can provide both unique insights into the mechanism of a problem and predictive power beyond experimental or clinical preparations. Once developed, a model can be used to test and generate future hypotheses in a way that may not be possible in experimental settings. The holy grail of computational biology is to develop comprehensive models that describe both mechanistic properties-for example, detailed molecular dynamics of biochemical interactions in a living organism-and subsequent emergent phenomena.
Development of comprehensive models is constrained by current computational power as well as lack of data. Instead of building comprehensive models, a more adaptable approach is to select a relevant spatial and temporal scale of the phenomenon and devise models that suit a particular research question. Thus, for the same biological or physiological systems, a wide spectrum of models may coexist, which aim to explain and predict the physiological process at different length or time scales. By analysing these models separately, researchers can gain a deeper understanding of the regimes where the different quantitative models are adequate. Bridging these models provides a way to propagate outputs derived from one model into inputs for another model.
In computational cardiac electrophysiology, there exist a range of models, which have been used to examine how subcellular electrical processes influence the diffusion of activation wavefronts across the heart [1]. Computing detailed biophysical models involves solving large systems of coupled ordinary or partial differential equations, which is computationally demanding. This limits both the number of simulations that can be run as well as their duration. It is therefore difficult to explore the sensitivities of a given model to input parameters and initial conditions. This, in turn, means that model outputs cannot easily be translated into inputs for models at other scales, such as those describing progression of patients through care pathways [2].
Atrial fibrillation (AF) is a cardiac arrhythmia that remains poorly understood despite progress in the development of detailed cardiac electrical models, experimental work and clinical studies. AF presents a prevalent heart rhythm disorder that significantly increases stroke risk in humans [3]. Improving identification, management and treatment of AF remains an important challenge [4]. AF consists of episodes of rapid and self-sustaining electrical excitation in the atrium of the heart, which punctuate periods of normal sinus rhythm when activation is driven by the heart's natural pacemaker. As the disease develops, episodes of AF become longer and more frequent until AF becomes permanent. Episode duration can vary between seconds and weeks, and constitutes the basic clinical marker to classify AF progression in patients. The ability to model and predict episode duration for a given patient would therefore be of significant clinical interest. In a previous publication [5], we described a biophysically motivated agentbased stochastic model to simulate progression of AF in a patient from first diagnosis, based on generating a time series of AF episodes with varying durations. The model parameters, which predicted episode start times and durations, were estimated from the literature where possible.
The duration of an AF episode is determined by mechanisms underlying its initiation and termination. AF is thought to be driven by rapid and self-sustaining electrical activity predominantly occurring in the left atrium [4,6]. An important mechanism is re-entry, in which a circulating activation wave continually propagates into recovering tissue. Several mechanisms have been associated with re-entry initiation, including atrial fibrosis [7], pulmonary vein triggers [6], and action potential and conduction velocity restitution [8]. Meanwhile, mechanisms of AF termination remain a poorly researched area, in part, owing to the computational cost of evaluating complex cardiac electrophysiology models over long periods. Initiating and maintaining an AF episode up to its termination in a simulation of a biophysically detailed model requires significant computational resources [9] if the episode lasts for more than a few seconds and/or a detailed atrial geometry is used.
Cellular automata (CA) models of the electrical activity on the surface of the heart are a simplified representation of cardiac electrophysiology, and the very first computer model of AF was a CA model [10]. They provide an intuitive way of describing how cardiac cells activate (depolarize) and deactivate (repolarize) by using simple update rules for the state of a single cell. These are usually based on the present states of the cell itself, and of its nearest neighbours. Since the original five-state Moe model [10], a series of studies have established that CA models of electrophysiology can represent behaviours seen with more detailed biophysical representations of excitable media [11][12][13][14]. CA models are simple to program and computationally cheap to run, allowing large numbers of simulations for little cost, and more detailed CA models have recently been devised [15,16].
The motivation of this study was to use a CA model as a computational platform to investigate stochastic initiation and termination of AF episodes. Our contribution can broadly be summarized as follows. First, in contrast to previous CA models [15,16], which used a simplified geometry such as a two-dimensional sheet with periodic boundary conditions, we generalized to a geometry representing the anatomical topology of the left atrium. While this is still stylized, we think this is a step towards reality. We show that the model is capable of inducing and terminating AF episodes stochastically. These phenomena are in line with the predictions of current state-of-art mechanistic models, and we are confident the CA model captures the essential dynamics of the real physiological system. Thus, we propose that CA models are a reasonable compromise between reality and computational efficiency when large numbers of long duration simulations are required. Second, we present a framework to analyse and infer the rate of stochastic initiation and termination of AF episodes. With the ability to run large numbers of simulations over long durations, we were able to accurately quantify these rates. This is necessary in order to be able to predict-in a statistical sense-the future progression of patients at a longer time scale. For example, these rates can be used to connect the CA model to the model we proposed to represent AF progression over years and decades [5]. We propose a framework of statistical analysis of patient trajectories, and apply it to a set of patient trajectories, generated from the CA model. We believe the ideas suggested may also be applicable to data from mechanistic models of other physiological systems, when computational resources are available to generate sufficiently many sample paths from such models.

Model geometry
Electrical activation was modelled on an idealized spherical geometry, representing the left atrium of the human heart. We did not include the right atrium, because the main drivers of AF are believed to originate in the left atrium. In order to approximate the dimensions of the human left atrium, the volume of the sphere representing the left atrium was set to 40 ml [17], corresponding to a radius of 21.2 mm. We rescaled this geometry to obtain a unit sphere.
We implemented a Moe-type CA [10] in which the dynamics take place on discrete nodes on the surface of the sphere. To place the nodes on the spherical surface as uniformly as possible, we used an icosahedral dissection [18] to distribute 10 242 points evenly on the sphere, as visualized in figure 1. We also investigated an alternative way to distribute nodes using an rsif.royalsocietypublishing.org J. R. Soc. Interface 14: 20160968 Archimedian spiral [19]; this method can be generalized to nonspherical surfaces.
Polar coordinates (u, f ) were used to specify the locations of the nodes. We defined the anterior and posterior direction to be (u, f) ¼ ( p/2, 0) and ( p/2, p), respectively.
On the geometry, the anatomical objects-four pulmonary veins (PVs) and the mitral valve (MV)-were set to be electrically inactive. The MV, modelled as a circular area centred on the south pole (u, f) ¼ ( p, 0), was estimated to have circumference of 85 mm [20]. The four PVs were modelled as circular areas with a base radius of 5 mm [21] corresponding to 0.236 scaled units. The PVs were placed at (u, f ) ¼ (2p/5+p/10, +p/3). Nodes in these areas were permanently removed, and the remaining nodes comprised the substrate for the CA to take place.

Fibrosis
Fibrosis on the posterior atrial wall is thought to play an important role of inducing AF re-entry [1,7,22], and was modelled by removing nodes in the corresponding area. To model the spatial heterogeneity of fibrosis, we removed nodes according to a probability distribution set to be normally distributed, centred at (u,f ) ¼ (0.65p, 0), with a standard deviation equal to 0.4 sphere radii. The number of nodes removed (denoted FC) quantified the severity of fibrosis. Time-dependent fibrosis was not investigated in this study, as structural modelling of atrial tissue with fibrosis occurs at a time scale much slower than that of re-entrant activity [23].

Dynamics of the cellular automaton model
A multi-state Moe-type CA was used to represent electrical excitation in each node (or 'cell'). Each node on the sphere could be in one of a number of discrete states, labelled 0, 1, 2, . . . . In this type of discrete-time model, an action potential is represented by a time delay, during which an excited cell may trigger neighbouring cells within an interaction radius but cannot itself be reexcited. In our model, the cell was deemed 'at rest' at state 0 and 'activated' if its state was greater than 0. A cell at rest would become excited if the number of 'recently excited' neighbours in a local radius exceeded a threshold, upon which it changes from state 0 to state RP, the RP or action potential duration, at the next time step. A neighbour was considered to be 'recently excited' if it had been activated in the past four time steps. This number was chosen to achieve realistic spread of excitation (see below). Following excitation, the activated cells reduced their state by one each time step until the state reached 0, i.e. the 'rest' state. Each discrete time step in our simulation corresponds to approximately 2.5 ms in real time. In sinus rhythm, RP took values of about 120 time steps in the model (variations are described below), this representing a physiological RP of 300 ms.
To avoid grid discretization effects on the simulations owing to non-uniformities of the icosahedral mesh, the interaction radius between cells on the sphere was set to be greater than the length scale of the typical internode spacing (for complications, see Ventrella [24]). The speed at which an excitation wavefront could propagate (conduction velocity) was determined by two free parameters: the search radius and the threshold of number of active neighbours. We carefully calibrated both the active neighbour thresholds (¼ eight nodes) and local search radius (2.544 mm), corresponding to a region containing %36 nodes to achieve a baseline conduction velocity across the sphere of 0.5 m s 21 . Thus, the total time taken to travel across the unit sphere (defect-free) from north pole to south pole was %133 ms.

Sinus rhythm
The sinus node (SN) is located in the right atrium, so in our model, sinus rhythm was represented by the regular activation of a region of cells proximal to the right pulmonary veins (a circular area centred at (u, f) ¼ (5p/12, p/2) with radius 1.696 mm), which is typically the site of earliest activation in the left atrium following right atrial activation. The sinus period was set at 1 Hz for all simulations, unless otherwise specified.

Pulmonary vein triggers
Bursts of spontaneous activation near the PVs are thought to be triggers of re-entry [6]. To model PV bursts, a 2 mm annulus around each of the four PVs was set to be capable of autoexcitation. In each time step with probability p, one node in these annulus regions and its surrounding nodes (set as those points within 2.12 mm to the selected node) spontaneously fired to its maximal state. The location of this spontaneous firing was chosen uniformly on the annuli. The probability, p, quantifies how often these bursts occur; the corresponding burst rate in a continuous-time setting can be computed using p/(time step) ¼ continuous-time bursting rate BR, which is set as a model parameter. Note that triggers were stochastic and occurred on average BR times per second, rather than occur periodically every one per BRs.

Restitution
To model the effect of restitution where the RP (i.e. action potential duration) of a cell shows sensitivity to its previous rate of excitation, we implemented the following (non-dimensionalized) formula [25] where DI is the diastolic interval (the quiescent interval between the end of one activation and the following beat), and B, K are parameters controlling the steepness of the curve. K had units of discrete time (¼2.5 ms), and B was dimensionless. The floor function enforced that RP was an integer, which in combination with scale factor 121, allows a maximal RP of 120 time steps. RP was subject to a minimum of 64 time units, i.e. RP ¼ max(equation (2.1),64). This equates to 160 ms (considered the shortest physiologically relevant RP). We investigated the dependence of the transition rate into AF episodes on parameters B and K.

Implementation
The model was implemented with custom code written in Cþþ, and is publicly available on Github at https://github.com/ dblueeye/atrial-fibrillation-cellular-automata. Links to sample movies may also be found. The simulation ran at 16Â speedup, i.e. 16 s of simulated time could be evaluated in 1 s. The skeleton code of the simulation is detailed below to clarify implementation steps: (1) Initiation: set up location of the nodes on the sphere. Remove nodes on areas of the PV and MV. For each sample run, model fibrosis by removing a fixed number of nodes according to a spatial probability distribution. Briefly, assign a probability to each node, generated from a normal distribution centred at (u, f) ¼ (0.65p, 0) with a standard deviation 8.48. Then, arrange the probabilities into a list and compute the cumulative probability distribution F(i) with respect to the list. The inverse transform sampling was applied to the discrete distribution to select the node to be taken out. Repeat the procedure until FC numbers of nodes were taken out. Generate and store the nodes representing SN breakthrough. Generate and store a list of possible PV bursting locations and the nodes which would burst in a group. Generate a neighbourhood map between the nodes. (3) PV bursts: with probability p, there will be a PV burst. If this happens, choose one of the locations where PV bursts can take place. As described above, a group of nodes in that region is activated to their maximal state, and the new RP is computed and updated using equation (2.1). (4) Rest of the nodes: for the remaining nodes, check if any neighbours in the interaction range have been activated in the past four time steps (10 ms). If so, this node is activated to its maximal state, RP is again updated according to equation (2.1). Otherwise, the state of the node is reduced by 1. (5) Repeat from 1 until end of simulation.

Stochastic initiation of atrial fibrillation
Simulation results were visualized using an equal-area Mollweide projection [26], shown in figure 2. During sinus rhythm without PV bursts (figure 2a), activation of the left atrium began by SN breakthrough near the right PVs; wavefronts passed around the larger PVs smoothly with a conduction speed of 0.5 m s 21 . When wavefronts passed through areas of fibrosis, conduction slowing and conduction block were observed occasionally when the number of activated nearest neighbours remained subthreshold. When PV bursts were introduced, triggers initiated activation near single PVs at a constant rate; in some simulations, this led to transient re-entrant wavefronts forming, and in certain cases, these became permanent re-entrant wavefronts (figure 2b). Movies were uploaded to YouTube, and can be found on the Github project page, see the supporting information.
While the complete course of the stochastic process (for each node) could be stored, the resulting data file would be impractically large. Instead, we evolved the CA without exporting the dynamic states at all time steps. As our aim was to investigate statistical properties of the system initiating and terminating AF (defined as self-sustained activity differing from sinus rhythm), two AF classifiers were developed. We stored only the seeds of the pseudo-random number generator of those sample paths, which were classified as 'in AF' (details described below). If needed, the collected seeds could recreate the sample paths for subsequent analyses. This procedure permitted generation and storage of up to 10 4 -10 6 simulation runs, necessary to accurately compute the statistics of AF episodes, including the sampling of rare events such as termination.

Probability of initiating spontaneous re-entry and atrial fibrillation
An exploration of the model parameter space was undertaken to determine the primary mechanisms of re-entry initiation. Each simulation was started in sinus rhythm (by setting p ¼ 0), then PV bursts of varying time durations were initiated by setting p;BR/400, to simulate PV triggers on the domain. Following a period of time with PV bursts, p was reset to 0, and the simulation evolved for a further 10 s observation    sample paths leading to re-entrant wavefronts determined the probability of a given parameter set causing re-entry. We varied the following parameters: number of fibrotic cells (FC), PV bursting rate (BR) and restitution steepness (B, K). Baseline simulation parameters were: FC ¼ 300 (points), BR ¼ 20 Hz, B ¼ 1.0 and K ¼ 40 (discrete-time unit, ¼ 100 ms), and we point-mutated the parameters (FC,BR,B, K)-see table 1 for the full range. For each parameter set, PV burst duration was varied from 1 to 5 s, and 10 5 sample paths were generated to compute the probability of inducing re-entry. Results are summarized in figure 3b,c.
We found the probability of AF re-entry depended linearly on the duration of the PV bursting when this duration was less than or equal to 5 s. This suggests AF initiation may be modelled by a simple coarse-grained model in continuous time, in which initiation occurs with constant rate, written as follows: Sinus rhythm ! AF with a rate r 1 , ð3:1Þ where the transition rate r 1 is the slope of the linear response shown in the left panel of figure 3b,c. We performed a linear fit to the numerical data, and found the rate was monotonically dependent on fibrosis: r 1 is larger for higher FC. Estimated values for r 1 are reported in table 1.
We found a non-monotonic relation between re-entry probability and PV burst rate BR, seen in the middle panels of figure 3b,c. A PV burst was most likely to induce AF when BR was between 20 and 40 Hz. This could be due to 'crowding' effects in the case of BR, reducing excitable regions or increasing the likelihood of wavefront collision and termination of reentry. Similarly, the rate into re-entry had a non-monotonic response to the restitution parameters, increasing AF initiation rate then decreasing as the parameters were increased, steepening the restitution slope, as shown in the right panels of figure 3b,c. For increased restitution steepness (larger (B, K)), a potential mechanism could be that increased spatial heterogeneity of RPs is mediated by shorter refractory tails. These observations suggest that the CA model is able to capture complex interplay between the mechanisms inducing AF.

Estimating time of atrial fibrillation initiation using a dynamic classifier
In the previous section, we investigated the hypothetical case where we controlled PV burst duration independently and subsequently observed for AF. In reality, PV bursts occur at random and cannot be simply turned off physiologically-AF may have initiated before the end of the burst period. Thus, the previous classifier is insufficient for estimating the true time of AF initiation. An alternative classifier to observe, record and track re-entry was thus proposed to estimate AF initiation time.
To model this, we again used the proportion of activated cells to be our 'signal', and defined an alternative AF classifier: tracking the proportion of activated cells out of total (non-fibrotic) cells, above a non-zero threshold for a period of time. We considered this analogous to clinical monitoring methods such as the electrocardiogram, which detect absence of regular peaks (e.g. P waves) for defined periods. Similar methods have been adopted by Manani [16]. In the following analysis, we set the non-zero threshold to be 0.5 and the time period to be 2 s. Using this definition, the classifier operates without perturbing the CA, and the onset time of AF reentry is a random variable: in different sample paths, the first time the classifier is triggered is random.
We refer to the first time the classifier indicates AF as t. This differs from simulation to simulation, and is random. We simulated 10 5 samples for selected sets of parameters to compute the cumulative distribution function s(t) ¼ Prob[t . t], the probability that the classifier is not activated before time t. This monotonically decreasing function quantified the statistics of the random transitions into the first re-entrant episode: the quicker the cumulative distribution function decays, the faster the system transitions to AF on average. The results are Table 1. Numerical results of best fits to the model results using the parameter sets (first column), for simulations up to 5 s using a two state linear transition model (second column), and up to 300 s using a two state nonlinear model (third and fourth columns). r 1 is the rate of AF initiation in both models, and r 2 is the rate at which initiation is inhibited for the nonlinear model. In specific parameter sets, r 1 , r 2 which suggests a very low rate of AF initiation.
parameter (FC,BR,B,K ) best fit r 1 in figure 3b,c best fit r 1 in figure 5a,b best fit r 2 in figure 5a,b presented in figure 4. The numerical results suggest that the cumulative distribution function is exponential, a signature that the waiting time distribution is also exponential, confirming the simple coarse-grained model with constant rate in equation (3.1). This analysis suggested a monotonic relation to all parameters, which differs from the first classifier, where burst rate BR showed non-monotonic dependence. We observed that the increased number of PV bursts at high bursting rates raised the proportion of activated cells, triggering the second classifier, but without leading to AF under the definition of the first classifier. To test this observation, we performed the following simulation: after the second classifier identified a re-entrant episode, we turned off the PV bursts and evolved the system for another 10 s. We excluded the sample paths which did not exhibit re-entry at the simulation endpoint, following the first classifier. For all parameter sets except the high BR ¼ 80 Hz case, more than 94% of re-entrant episodes identified by the second classifier led to self-perpetuating re-entry. In the BR ¼ 80 Hz case, only 16% led to re-entry.
This observation showed the second classifier, albeit realistic in practice, overestimated the transition rate into AF.   To test this assertion, we extended the analysis in figure 3b,c with a longer PV burst duration. For each parameter set, we simulated 5000 sample paths to compute the probability that the sample had transitioned into AF. The results are shown in figure 5. For some parameter sets, after a long period of PV bursting, the probability did not converge to 1.0 (e.g. for FC ¼ 500). In other words, the coarse-grained model equation (3.1) did not sufficiently capture AF dynamics when the PV burst duration was increased.
We thus generalized the coarse-grained model into a two-state model with a stochastic initiation and termination of AF under conditions of PV bursting, where SR represents sinus rhythm: SR ! AF with rate r 1 ð3:2aÞ and AF ! SR with rate r 2 : ð3:2bÞ Because we start in sinus rhythm, the initial probability of AF at t ¼ 0 is always 0 (and 1 for sinus rhythm). The temporal behaviour of the probability to be in AF can be calculated using standard methods [27], and we find Prob½in AF at time t ¼ r 1 r 1 þ r 2 ½1 À e Àðr1þr2Þt : ð3:3Þ A two-parameter fit was performed for each simulated parameter set, and the best fit is displayed in figure 5. Corresponding values of r 1 and r 2 are reported in table 1. The value r 2 quantifies the time scale at which stochastic AF initiation is inhibited by constant PV bursts. Comparing the relative values r 1 and r 2 , with high FC or low BR, inhibition of AF initiation dominated the process (r 1 , r 2 ) and the response of the termination rate to the parameters was also non-trivial.     Otherwise, the state of the system remains in SR/AF, respectively. Here, 1/k 1 and 1/k 2 quantify the average duration of the resting state (no PV bursts) and active state (with PV bursts), respectively. Short trains of bursting mean that k 2 ) k 1 . Selected parameter regimes were tested (data not shown) and preliminary results showed the coarse-grained model equations (3.4) and (3.5) faithfully projects the progression of the CA model for a range of parameter regimes. However, at much longer time scales ) Oð1=k 1 ,1=k 2 Þ, there were notable discrepancies. We investigated these differences in the following section, which suggests existence of higher-order states of AF dynamics.

Fourier analysis revealed higher order dynamics of atrial fibrillation
Our numerical simulations yielded many ( 10 4 ) sample paths which ended in AF re-entry. Fourier analysis was applied to sample paths from the baseline parameter set where re-entry was initiated (corresponding to the 'observation' phase of the time series shown in figure 3a, top panel). Results for 200 sample paths are shown in figure 6. For 5 s PV burst duration (figure 6a), we visualized sample paths along with the Fourier analysis, and observed that the dominant mode %5 Hz corresponds to the period of a single re-entrant wavefront. Subdominant half modes %2.5 Hz corresponded to the period of points that experienced 2 : 1 conduction block, e.g. points near one of the PVs which have previously been fast paced. There also existed higher harmonics, to which we did not seek to fit a physiological interpretation.
There was a notable variability in the Fourier spectrum for each sample path. This reflected the stochasticity of the system-including the quenched heterogeneity of RP, fibrosis and dynamical randomness from PV bursts-which propagated to the dynamics of re-entry modes. As the speed of the travelling wave is fixed at conduction speed 0.5 m s 21 , the dominant frequency is inversely proportional to the path length the wavefront travelled in one cycle. Both the duration of the re-entry and the length of cycle path exhibited %20% variability. We also examined the case when PV bursts lasting 150 s were applied (figure 6b), observing that the variability of the spectrum appeared smaller compared with the 5 s case. This suggests longer duration of PV bursts tend to drive the system into a stable dynamical mode that is hard to perturb. By comparing visualizations alongside the Fourier spectrum, we also identified that multiple 2 : 1 conduction blocks formed more frequently, and higher-order rotors were identified. Two representative snapshots are presented in figure 6a,b, top panel.
This analysis shows that even when the model state was classified as 'in AF', there can be multiple modes. The follow on question is whether the complexity of an AF episode affects its stability and its likelihood to terminate, either spontaneously or following intervention. For a single re-entrant wavefront, a short PV burst at the right time and location terminated AF (movie on YouTube, see the supporting information). This led to an investigation into spontaneous termination of AF in §3.7, comparing termination rates for different AF modes, to infer likely mechanisms of termination.

Investigation into stochastic atrial fibrillation termination suggests stable and unstable re-entry modes
Observation of simulations which generated figure 5 indicated stochastic termination of AF could be a direct result of PV bursts. To test this hypothesis, we randomly collected 500 sample paths ending in AF in previous experiments and performed 50 extended simulations on each. Recall that AF was induced by a set of PV bursts over some duration, say T 1 , in previous experiments. After AF was initiated, we waited a time window T 2 without PV bursting, and applied another set of PV bursts (1 s duration), and observed if re-entry was terminated after the second set of PV bursts had been applied. A schematic diagram of this is shown in figure 7a. Figure 7b shows the termination probability significantly depends on T 1 . For T 1 ¼ 1 s, it was very likely to terminate This analysis suggests two modes of AF: some re-entrant circuits can be terminated easily by PV bursts, and others cannot. In figure 7c, we show probability to termination, ordered by each sample path in the y-direction and each waiting window duration T 2 in the x-direction. A clear alignment in the x-direction of either blue or white stripes showed that if a sample path can be terminated, the probability of termination does not critically depend on T 2 ; if the sample path cannot be terminated, most likely, it cannot be terminated for any T 2 . The longer T 1 (the duration of the first set of PV bursts to induce re-entry), the smaller the proportion of unstable AF (episodes which can be terminated). Thus, the overall probability to terminate AF is orders of magnitude smaller than AF induced by shorter T 1 .
We therefore hypothesized that activation and termination can be modelled using the multiple-state model: and SR ! stable AF: ð3:6bÞ Results suggest that the transition rates are not constant and critically depend on the duration of PV bursts. To quantify transition rates, a classifier identifying the signal state must be developed; we aim to develop this in the future. Our presented framework can be applied to measure transition rates once a reliable classifier is implemented. We remark that the multiple-state system has a 'memory' for marginal observables (in AF or not) in line with our previously proposed hidden state binary model [5], which can be used to project the progression of AF over long time scales.

Discussion
In this study, we investigated stochastic onset and termination of AF episodes by using a CA model on a two-dimensional sphere, with an equivalent topology of the human left atrium. We demonstrated the capability of the model to generate large sets of sample paths to infer the statistical properties of initiation and termination of re-entry and AF (up to Oð10 6 Þ sample paths and for duration approximately Oð10 4 Þ s). Three potential arrhythmogenic mechanisms were investigated, fibrosis density (FC), pulmonary vein bursting rate (BR) and RP restitution steepness (B,K). By probing this parameter space, we investigated the probability of AF onset and termination resulting from PV bursts. We found a linear dependence between burst duration and probability of re-entry initiation for all parameters for short PV burst durations. Increased FC led to a linear increase in probability of initiating re-entry, but there was a non-monotonic relationship between probability of initiating re-entry and both BR and restitution steepness. One possible explanation is that while increased BR and steeper restitution both act to promote re-entry, these mechanisms may also act to promote the termination of re-entry by wavefront collision. When PV burst duration was increased, probability of re-entry at simulation endpoint did not increase linearly, such that at high FC and high BR, likelihood of re-entry remained constant. The behaviour of the model depends on parameters in a non-trivial way, which suggests the existence of complex mechanisms that inhibit or suppress AF initiation, and may even terminate re-entrant circuits before they have fully formed. By fitting a two state nonlinear model to our simulation outputs, we estimated initiation and inhibition rates r 1 and r 2 for given parameter sets. Finally, we analysed a subset of the sample paths in AF, and found existence of stable and unstable AF modes. A second set of PV bursts could spontaneously terminate a proportion of induced AF episodes, with termination probability reducing, subject to duration of the first PV bursts.
This study offers an alternative novel methodology and framework for investigating mechanisms of spontaneous AF, which differ from conventional modelling and experimental studies in its capability for rapid statistical sampling of longer duration episodes. Our findings and conclusions are set out and discussed in the subsequent paragraphs.

Cellular automaton models of cardiac electrophysiology
Although CA models have been superseded in popularity by more biophysically detailed models [28], CA models show similar tissue scale behaviours [11][12][13][14], and thus are still valuable both in stand alone theoretical studies and combined with clinical investigations [15,16,29,30]. Our approach complements studies by Manani et al. [16], who used a CA formulation to investigate the effect of time-dependent fibrosis on arrhythmia susceptibility. Our model represents variability and uncertainty through its stochastic formulation and the large number of sample paths, and thus permits a systematic investigation within the model framework, while accepting model limitations. A major limitation of the CA model compared with continuum models is its inability to directly model the mechanism of conduction slowing and CV restitution, although most other potential AF mechanisms [31] may be handled with the CA formulation.

Model parameter space
In this study, we fixed the size and shape of the left atrium, and the size and location of anatomical objects. We did not include the left atrial appendage, and assumed that location of SN breakthrough into the left atrium was fixed. Heterogeneity was investigated by randomly varying initial RPs, rather than by region specific heterogeneity in parts of the left atrium. We recognize that these are all parameters which may vary between individuals, and may significantly impact probability of AF initiation and termination. We chose to focus on biophysical mechanisms rather than on anatomical variability, and recognize that not only is the potential parameter space vast, but also that additional investigations into the effects of these parameters are important.

Atrial fibrillation onset
Recent studies have investigated mechanisms related to electrical and structural remodelling, highlighting the importance of inter-patient variability. McDowell et al. [22,32] found that combinations of fibrosis subtypes were proarrhythmic and that patient-specific distribution of fibrosis had a major impact on AF initiation, and anchored wavefronts to specific atrial regions, with other electrophysiological changes not significantly altering this behaviour. Krummen et al. [33] reported that steepening AP restitution slope in patients initiated re-entry, with the associated computational study identifying specific ionic pathways responsible for restitution steepening. Regional electrical heterogeneity of the atria was investigated by Colman et al. [34], who found region-dependent action potential duration heterogeneity in the atrium increased susceptibility to AF onset and maintenance of reentrant circuits.
Our study has investigated these three mechanisms plus PV firing rate, albeit with a discrete rather than with continuous model, and different assumptions and formulations (we did not model fibrosis subtypes or include region-specific RPs for our cells). Our study results differ from the conclusions of these continuum studies, especially regarding the steepening of restitution slope, where we found a non-monotonic relationship between AF onset and restitution steepness not predicted by Krummen et al. There is no general consensus on whether a steep restitution slope is pro-or anti-arrhythmic [35], and our results showed there is a 'window' of steepness which maximizes probability of AF onset. This was also true for the other parameters, where excessive fibrosis and PV burst rate inhibited increased onset of AF. We comment that a PV burst rate up to 80 s 21 , while representing the number of triggers across all four PVs rather than a single focal source, may appear unphysiological, but it is also possible that many focal PV bursts go undetected.

Atrial fibrillation termination
Clinical studies predominantly investigate how targeted ablations terminate AF, and these have been explored theoretically in a number of studies [36]. However, few studies explore spontaneous termination owing to the difficulty of capturing such rare events. A few clinical studies have been documented: Ndrepepa [37] referred to generators of fibrillatory activity in the left atrium, and reported that AF termination was polymorphic in its mechanism. Alcaraz [38,39] analysed the atrial activity of patients during AF and immediately prior to termination, and found the existence of more organized atrial activity (measured by sample entropy) 1 min prior to termination, and that the late activity had a significantly lower dominant frequency mean value. Some studies of dominant frequency and harmonics have suggested Fourier analyses as useful predictors of termination [40].
Our study was inconclusive regarding termination. We found that PV bursts are a potential mechanism for terminating as well as initiating AF, and also act to inhibit initiation rate for longer durations of PV bursts. Fourier analysis of the sample paths revealed both stable and unstable modes of AF, but no clear trend was observed. We found, however, that the longer the period of PV bursting, the smaller the probability that induced AF will be terminated by future PV bursts. This suggests dynamical memory effects exist within the model caused by extended PV burst pacing, which influences the stability and robustness of the induced re-entry wavefronts. This agrees with the 'AF begets AF concept' [41], and recent studies of Uldry et al. [42,43], who reported an increase in AF complexity with duration, and that spontaneous termination mechanisms differed depending on dynamics of AF and its underlying complexity.
In other recent studies, Krogh-Madsen et al. [44] also found that remodelling maintained AF by shortening atrial wavelength (electrical by shortening action potential duration, structural by slowing conduction), which correlated with increased AF episode duration, with dynamics of re-entry differing between types of remodelling. Liberos et al. [9] suggested that cell-cell ionic differences as a mechanism of AF termination, by decelerating re-entrant activity and increase in rotor tip meandering. Our model did not include electrical remodelling similar to these studies, but our model is well placed to analyse atrial wavelength and track the rotor tips in future studies, to see if similar mechanisms exist within our formulation. The general consensus is that AF complexity increases over time together with AF episode durations, with size of atria and atrial obstacles thought to play a critical role in termination. Petrutiu et al. [45] found that non-terminating episodes exhibited larger dominant frequencies compared with spontaneously terminating episodes, and more abrupt changes in dominant frequency were observed prior to spontaneous termination. An open mechanistic question remains over whether spontaneous termination is preceded by a progressive fusion of wavelets or a simultaneous block of all wavelets in the tissue. We believe our work is well placed to evaluate these questions through the capabilities to run longer duration simulations.
Our study may additionally complement existing ablationbased termination studies by identifying similar mechanisms or proposing novel therapeutic studies. Rappel et al. [46] demonstrated that ablation caused termination in a heterogeneous domain by creating an excitable gap, dislodging a stable anchored wavefront or by closing critical isthmus channels. Uldry et al. [47] reported 10-20% success rate when using atrial septal pacing at alternating frequencies to pass the atria.

Clinical relevance
We have presented a methodology to translate emergent simulated behaviour in a greatly simplified computational model of atrial electrophysiology into transition rates into and out of AF. Determining rate transitions between sinus rhythm and AF is an important step forward because these rates can be used in models of AF progression [5], can be evaluated against clinical data, and in the future could be used to predict AF onset and progression.

Future work
Our work in this article focuses on the framework of the stochastic analysis. We acknowledge that CA models are a simplified representation of reality. However, this approach permits large numbers of simulations to obtain probability distributions and probe particular mechanisms. We propose several areas in which the model could be investigated further.

Geometry
We adopted a simplified quasi-spherical geometry to model the left atrium. Because the dynamical rules of the CA only involve the neighbourhood relations between nodes, it is straightforward to construct a CA model on any two-dimensional surface embedded in three-dimensional space. The difficulty of evenly distributing the nodes may be overcome by using the Archimedean spiral [48]. It may also be possible to extend this to three dimensions.

Directional fibrosis
In this work, we modelled fibrosis by setting nodes to be electrically active, whereas fibrosis may act to promote faster propagation in certain directions within cardiac tissue [16]. This could be achieved in our CA model by assigning weights to neighbouring nodes.

Representing interventions
As the computational cost of a CA model is lower than that of biophysically detailed models, it is an ideal platform to develop and evaluate effects of intervention strategies such as ablation or external pacing. However, as a coarse-grained approach, the CA model is unlikely to capture detailed biochemical or biophysical effects within these, or within other interventions such as pharmacological modification of cell and tissue electrophysiology.

Restitution and remodelling
Restitution changes may not be instantaneous. One way to model restitution with memory would be to replace equation (2.1) by where a measures the strength of the 'memory'. When a ¼ 1, there exists no restitution, and when a ¼ 0, it converges to our proposed model (2.1).
Only initial state structural remodelling was investigated in this study. Additional structural and electrical remodelling may be implemented in the CA framework, both as an initial condition and as a transient process (e.g. with ageing). For example, removing cells (and adding them back) from the domain of excitable cells could model acute scar formation or recovery from ischemia.

Pattern recognition of the re-entrant wavefronts
Our analysis revealed stable and unstable modes of AF. Visualization of selected sample paths suggested some characteristic differences between these modes: just prior to termination, unstable AF terminated via conduction block through fibrosis regions or pulmonary veins. This often included spontaneous PV bursts at the channel isthmus in a short excitable window. In comparison, stable (did not spontaneously terminate) wavefronts appeared to have more complex pathways of propagation.
While the Fourier spectrum suggested potential differences, the analysis was inconclusive as there was a large variability over sample paths in a given parameter set. As our classifiers contain only the temporal information, we could additionally use spatial information of the re-entrant wave front to (e.g. rotor tip tracking, local electrogram) to inform our analysis.

Supporting information
Please visit Github project /dblueeye/atrial-fibrillationcellular-automata for a working implementation and for movie URLs on YouTube. Competing interests. We declare we have no competing interests. Funding. We acknowledge support from the UK Engineering and Physical Sciences Research Council (www.epsrc.ac.uk) grant no. EP/ K037145/1.