Unlike myocytes, fibroblasts are not electrically excitable, however they can maintain a resting membrane potential (RMP), and hence affect and be affected by the tissue electrophysiology. The electrophysiological function of fibroblasts is determined by the presence of a range of ionic channel currents in their cell membrane. It is known that both cultured and freshly isolated fibroblasts express voltage dependant potassium channels (I_Kv) (Chilton et al., 2005). A wide variety of currents have now been discovered in the cardiac fibroblast, predominantly potassium (I_Kur, I_Kir, I_to), but also non-specific ionic currents (I_ns) which can be mechanically and non-mechanically activated chlorine currents (I_Cl). The sources of data for these various currents are diverse, and are summarized in Table 1.

The general approach to modeling atrial fibroblast electrophysiology was based on a well-established ventricular fibroblast model by MacCannell et al. (2007). However, a set of ionic channel currents specific to atrial fibroblasts was considered, based on the CRN formulation of the currents for human atrial myocytes and experimental data from atrial fibroblasts. Data for the atrial fibroblast model was obtained from a study (Wu et al., 2014) on ionic currents in freshly isolated canine atrial fibroblasts. An optimal fit between the modeled ionic currents and experimental data was achieved by varying the current conductances (see Results below). Note that the dog has been widely used as an experimental model for AF (Schotten et al., 2011), and in the absence of ionic current data from human fibroblasts it is sensible to propose that the canine data provides the closest match.

The 4-aminopyridine (4-AP) sensitive current was modeled as I_{to, f}, which is selectively blocked by 4-AP. Fitting the current-voltage (I–V) relationship for this current to experimental data by Wu et al. (2014) involved a 90% reduction in the current conductance compared to its original formulation in the CRN model, but no changes to the CRN current kinetics was introduced.

The S9941 sensitive current was modeled as I_{Kur, f}, as S9941 is a blocker of ionic channels containing the subunit Kv1.5, which forms the channel conducting this current. Fitting the I-V relationship for this current to experimental data by Wu et al. (2014) involved a 40% reduction in the current conductance compared to its original formulation in the CRN model.

The end pulse current was measured by Wu et al. (2014) as the total voltage sensitive current. We considered the remainder of this current after subtracting both I_{to, f} and I_{Kur, f} as a sum of a further two currents: a linear non-specific ionic channel current, I_{ns, f}, and the inward-rectifier potassium channel current, I_{K1, f} (see Results below). Both I_{Kur, f} and I_{ns, f} are major contributors to the fibroblast RMP. The separation of these two currents ensured that I_{K1, f} reversed at the equilibrium Nernst potential for potassium ions, E_K, as is observed in electrophysiological recordings of this current in multiple cell types. The voltage and time dependences for these currents are based on the CRN formulation (index _f denotes their relevance to fibroblasts):

Here, V_f is the membrane voltage of fibroblasts, where o_{a, f,} u_{a, f} and o_{i, f,} u_{i, f} are activation and inactivation gates, respectively, and g_Kur and g_to are conductances of the respective channels.

The other two currents in the developed atrial fibroblast model are based on the MacCannell et al. (2007) formulation. Specifically, there are currents not measured by Wu et al. (2014) but assumed present in the atrial fibroblast based on experimental studies with other cardiac cells. These are the sodium-potassium exchange current (I_NaK) and the background sodium current (I_bNa).

Note that in the MacCannell et al. (2007) model and its derivatives, only I_K1 and the time and voltage dependant potassium current I_Kv are included. The latter current was shown to be sensitive to 4-AP and S9941 in their experimental work, (Shibukawa et al., 2005) and hence is likely to present a sum of I_{to, f} and I_{Kur, f}. However, unlike our novel model, MacCannell et al. (2007) modeled it as a single current.

M-F electrotonic coupling mechanism is based on the presence of M-F gap junctions (Rook et al., 1992; Gaudesius et al., 2003; Camelliti et al., 2004), in the same manner as inter-myocyte coupling. In order to simulate the M-F coupling, we used a similar formalism to that developed by MacCannell et al. (2007). Specifically, the CRN atrial myocyte (with the total membrane ionic current I_tot) was coupled with the novel atrial fibroblast (with a total membrane ionic current I_{tot, f}) via a linear M-F gap junctional conductance (G_gap):

Here C_m is the membrane capacitance, subscripts m and f relate to a myocyte and fibroblast, and n to the number of fibroblasts coupled to one myocyte. G_gap was set as 0.5 nS (Jacquemet and Henriquez, 2008) and n was varied between 0 and 6. Unless specified otherwise, a myocyte was coupled to 2 fibroblasts.

The images with distributions of atrial fibrosis used in this study were obtained from 3 patients suffering persistent AF and recommended for routine, first time pre-ablation imaging. These data were randomized to ensure patient privacy. A Philips 1.5T MRI scanner was used, and the images were obtained at 1.4 × 1.4 × 4 mm³ resolution and saved in DICOM format. Two MRI modalities were used: LGE MRI (Figures 1A–C) and early gadolinium (Gd) angiographic imaging (gated MRA) (Figure 1D). The angiography images were utilized for an improved segmentation of LGE MRI data (see details below). These were taken in the same cardiac phase (atrial diastole) with respiratory gating.

A novel image processing pipeline was developed to reconstruct and register the LGE MRI data on the fibrosis distributions with the 3D atria. The pipeline utilizes fibrosis segmentation tools developed by the CARMA group at the University of Utah, such as plugins to 3D Slicer™. This enables the reconstruction of gradual fibrosis changes as a function of image intensity.

The extraction of fibrosis data from MR images requires thresholding of the intensity of the atrial wall. Obtaining the atrial wall in MRI data is non-trivial (Varela et al., 2015) and detailed segmentation can only be performed manually due to the relatively low resolution and signal to noise ratio (SNR). This process is extremely time consuming and may provide discrepancies between segmenters. However, the LGE MR images are acquired in the same session and at the same cardiac phase as a gated MRA acquisition, used as a navigation image for planning ablation procedures. The gated MRA delineates (Figure 1D) the blood pool which has a high intensity due to the contrast agent (Gd) and allows for accurate semi-automatic segmentation of the endocardial surface. This is achieved using a region growing algorithm within ITK-Snap, which utilizes seed points that iteratively dilate. The growth is limited by image intensity gradients. The segmented endocardial surface then is dilated to provide an epicardial surface. In some cases it required manual rectification in highly convoluted areas. The dilation is usually between 3 and 5 mm, which corresponds to the variable atrial wall thickness (Varela et al., 2015). The endocardial volume is subtracted from the dilated epicardial volume creating a 3D wall mask.

Utilising tools developed by the CARMA group and available as a 3D Slicer plugin (http://capulet.med.utah.edu/namic/cmrslicer/), provides two different methods for fibrosis segmentation. The first, named “automatic scar segmentation,” uses a K-means clustering algorithm. This produces a binary fibrosis model in the masked region, identifying the densest fibrosis areas. The second method (“threshold model”) collects the MR intensity data from the masked image applying no threshold or segmentation. Both methods project the data onto a surface mesh of the endocardium. This 3D image-based atrial mesh is registered to a respective 3D mesh of the Visible Human Female atria using the IRTK-toolkit (https://biomedia.doc.ic.ac.uk/software/irtk/) using a non-linear deformation. The registered 3D mesh is projected onto the Visible Human dataset using a “nearest neighbor” algorithm. Fibrosis datasets obtained using the “threshold” tool are segmented by LGE MRI threshold, with the threshold values defined on the range of the image intensity in the masked image.

The range by LGE MRI intensity values corresponding to fibrosis were calculated individually for each patient, due to the intensity variation between scans, using the maximum intensity (MI) as the upper limit (Figure 2B). We hypothesized that the histogram intensity (Figure 2B) was a result of an overlap between two intensity distributions corresponding to healthy and fibrotic tissues (McGann et al., 2008). Thus, the right tail of the histogram was associated with fibrosis, as described in previous studies (Figure 2A), and the left tail of the histogram was associated with healthy tissue.

The range of intensity values was divided into 7, with the lowest 2/7th being allocated to healthy tissue and each fibrosis level, 1–5 allocated as a following 1/7th of the range (moving toward higher values). Due to the densest fibrosis (level 5) covering less tissue than has previously been associated with fibrotic scars (Karim et al., 2013), the highest 1/7th of the range was expanded by 20% (moving toward lower values), with all other parts of the range evenly reduced. This provided significant regions of dense fibrosis in agreement with previous reconstructions. The resultant intensity distribution between levels was as follows: healthy tissue, 0 to 0.27 MI; level 1, 0.27–0.4 MI; level 2, 0.4–0.55 MI; level 3, 0.55–0.69 MI; level 4, 0.69–0.82 MI; level 5, 0.82–1.0 MI. The segmented fibrosis regions were considered to be transmurally uniform in the 3D atrial model.

To investigate the role of segmentation variation, we also created two alternative segmentations for each LGE MRI dataset, one with a higher and another with a lower threshold. The higher threshold case in described above, and the lower threshold case was created by shifting the threshold of dense fibrosis (dashed line in Figure 2B) leftward by 10 units of LGE MRI intensity.

For the whole atrial simulations, the monodomain equation was solved on the 3D human atrial geometry with zero-flux boundary conditions, using the forward Euler method with spatial and temporal steps of 0.3 mm and 0.005 ms, respectively. The equation is as follows:

Here ∇ is the gradient operator and D is the diffusion coefficient, which can be a tensor when the tissue is anisotropic. The model uses the 3D fiber vector array to calculate the diffusion tensor to each voxel dependant on the orientation of the fiber (Clayton et al., 2011; Colman et al., 2014). This calculation also requires the longitudinal and transverse diffusion coefficients, D_L and D_T, where D_L was constant at 0.3 mm²/ms and D_T was varied to change the anisotropy ratio, R = D_L /D_T, as described below.

The model settings were varied to simulate different conditions linked to fibrosis and AF. Four types of AF-associated remodeling were considered: (1) Ionic remodeling (Colman et al., 2013, 2014), but no structural remodeling (the standard values of D_L = 0.3 mm²/ms and R = 10:1), which is associated with early-stage AF. (2) Diffuse fibrosis, where instead of ionic remodeling the M-F coupling was added and inter-myocyte coupling was decreased (D_L = 0.12 mm²/ms, R = 10:1) throughout the atria, both changes associated with short-term AF (Rook et al., 1992; Li et al., 1999). (3) Interstitial fibrosis, where in addition to the presence of diffuse fibrosis atrial tissue anisotropy was increased (D_L = 0.12 mm²/ms, R = 16:1), which is associated with high level of collagenous depositions in long-term AF (Kawara et al., 2001; Krul et al., 2015). (4) Patient-specific patchy fibrosis (D_L ≤ 0.3 mm²/ms, R = 10:1) introduced as described below.

The segmented distinctive fibrosis levels, indexed 0–5 depending on the LGE MRI intensity, were applied to simulate the dense-to-diffuse fibrosis gradients in the 3D atrial model. Index 0 represents no fibrosis, and regions with indices 1–5 are modeled as increasingly severe fibrosis. Specifically, indices 0–5 corresponded to progressively increasing M-F coupling (0–5 fibroblasts per myocyte) and decreasing diffusive coupling between myocytes (100–20% of the standard D_L = 0.3 mm²/ms). Thus, 3 patient-specific fibrosis distributions were generated.

The tissue effective refractory period (ERP) was calculated in a quasi-1D tissue by applying a stimulating current at one end to initiate AP propagation to the other end. This was modeled as a long 3D slab of 50 × 7 × 7 cells with a space step of 0.3 mm and a diffusion coefficient of 0.3 mm² ms⁻¹. The stimulus of −2500 mA was applied in the first 10 layers of cells. This resulted in a 1D plane wave propagating along the slab. Stimuli were applied at a given BCL, followed by a last stimulus applied at a variable S2 interval. The latter was increased until an AP was observed at the non-stimulated end of the cable. The minimum S2 interval for which AP was still able to propagate through the cable was taken as ERP. These simulations also enabled the calculation of the conduction velocity (CV) by measuring conduction times along the slab.

Re-entry was initiated in the 3D model through fast pacing in the RAA or LAA, allowing for the wave to propagate transversely to the crista terminalis (CT) and generate a conduction block.

The fitting of ionic channel currents in the atrial fibroblast models to the respective patch-clamp data (Wu et al., 2014) can be seen in Figures 3A–D. I-V curves for I_Kur, I_to and I_K1 (Figures 3A–C) in the model all are in good agreement with the experimental data. After fitting the ionic currents to experimental data, the RMP in the resulting model for a single atrial fibroblast was −42.5 mV, which was in excellent agreement with the respective experimental value of −42.8 mV recorded from atrial fibroblasts (Wu et al., 2014). This value is more positive than the RMP of −48 mV in MacCannell et al. (2007) ventricular fibroblast model, although both are within the experimental data range (Table 1).

Coupling to fibroblasts significantly affected electrophysiological properties of the myocyte. Figures 4A,B, illustrates APs in the control CRN model (i) coupled to the developed atrial fibroblast model (ii), and the CRN model for a myocyte coupled to the MacCannell et al. (2007) ventricular fibroblast model (iii). In these simulations, a myocyte was coupled to 2 fibroblasts. The M-F model (ii) produced prolonged repolarisation compared to the uncoupled control model (i), while the M-F ventricular model (iii) produces a more rapid repolarization than in control. Hence, coupling to atrial (ii) and ventricular (iii) fibroblasts produces opposite effects on AP in the CRN atrial myocyte model (i). Figures 4C,D show the respective membrane potential in atrial and ventricular fibroblasts for the M-F coupling models (ii) and (iii).

The fibroblast potentials were purely electrotonic due to the coupling with the myocyte and had substantially lower amplitude than that of the myocyte. The atrial fibroblast had a slower repolarization and slightly lower amplitude than that of the coupled ventricular fibroblast. More importantly, it had a higher RMP, and hence a larger effect of the myocyte RMP.

Figure 5 compares electrophysiological characteristics of the cases (i)-(iii). The RMP (Figure 5A) was increased (became more positive) compared to control in both M-F coupling models (ii) and (iii). At a BCL of 500 ms, in the M-F model (ii) with atrial fibroblast the RMP was −76.5 mV (Figure 5A) and in the M-F model (iii) with ventricular fibroblast the RMP was −77.5 mV, compared to −79.8 mV in the uncoupled CRN model (i). RMP is lowest in case (ii) throughout all BCLs, which can be particularly important at fast rates associated with AF.

The dynamic effects of coupling a CRN myocyte with the fibroblast models significantly changed the refractoriness of the myocyte. There was little difference in APD between the M-F coupling model with atrial fibroblasts (ii) and the control (i) except for at a high pacing rate, e.g., at a BCL of 300 ms (Figure 5B). In comparison, the M-F coupling model with ventricular fibroblasts (iii) reduced the APD throughout the range of BCLs, with a maximum reduction of about 20 ms. Effects of M-F coupling on APD became more pronounced in both models (ii) and (iii) as the number of fibroblasts coupled to a myocyte was increased from 2 to 4 and 6.

The excitation threshold for the stimulus current was also investigated (Figure 5C) for the models. In the M-F coupling model (ii) with atrial fibroblasts, the threshold was −15.9 pA/pF at the BCL of 500 ms. It was significantly lower when compared to both the uncoupled myocyte model (i) (−19.2 mA/mS) and the M-F coupling model (iii) with ventricular fibroblasts (−17.5 mA/mS). Changes of the excitation threshold in models (ii) and (iii) compared to control CRN (i) correlated with the respective changes in the RMP: higher (more positive, but lower absolute value) RMP corresponded to the lower threshold (see Figures 5A,C).

Effects of M-F coupling on ERP were larger than those on APD in both M-F coupling models (ii) and (iii). ERP in the M-F coupling model with atrial fibroblasts (ii) was larger than that in control (i) by about 5–10 ms at BCLs above 300 ms (Figure 5D), but dropped off rapidly at BCL was decreased below 275 ms. ERP in the M-F coupling model with ventricular fibroblasts (iii) was substantially larger than that in control (i) by about 50–100 ms across the range of BCLs.

Despite decreasing the excitation threshold, both M-F coupling models (ii) and (iii) reduced the CV compared to the control (i), with the ventricular M-F model producing the greater reduction. This can be explained by a substantial reduction of the sodium current, I_Na, in case (ii) (Figure 5F) and by a large increase of the ERP in case (iii) (Figure 5D). The decrease of I_Na in the myocyte, particularly pronounced due to coupling with atrial fibroblasts at high rates (Figure 5F), can be explained by more positive RMP in atrial fibroblasts (Figure 5A). The latter results in a stronger electrotonic M-F current, I_Gap, (Figure 4) and partial inactivation of I_Na.

Note that RMP in the novel atrial fibroblast model (−42 mV) was more representative of atrial physiology (Table 1) than that the MacCannell et al. (2007) ventricular fibroblast model (−47 mV), as it was developed based on more comprehensive patch-clamp data from atrial fibroblasts. Moreover, M-F coupling with the ventricular fibroblast model produced non-physiologically high ERP (about 350 ms) and low CV (0.4 m/s) in atrial tissue at high rates (BCL of 300 ms). The high ERP in this case may be a result of slow repolarisation processes, characterized by long time constants, in the ventricular models–whereas our novel atrial fibroblast model uses formalisms, including kinetics inherited from the CRN model, specific to faster repolarisation in the atria. Hence, using the ventricular fibroblast model instead of the atrial fibroblast model in 3D human atrial models could lead to different behavior in AF, as shown in the next section.

In 3D atrial simulations, fast pacing in the anterior RA resulted in a conduction block leading to wave breaks and re-entry in all simulated conditions (see Methods). The precise location and timing of the wave-block varied depending on the condition considered and pacing rate used, but was always clearly observed between the CT and the pacing site (seen in Figure 6).

In the diffuse fibrosis conditions, M-F coupling with atrial and ventricular fibroblasts resulted in greatly varying duration of re-entry: about 3 and 12 s, respectively. A much longer duration of re-entry in the latter case can be explained to the wavelength reduction associated with extremely low CV values seen with this model (Figure 5D) at high rates. The short wavelength can be seen in Figure 6 as a larger excitable gap (between a re-entering wave and the preceding one) in the case of coupling with ventricular fibroblasts (Figures 6A,C), as compared to that with atrial fibroblasts (Figures 6B,D). This also points the sensitivity of 3D atrial behavior to the choice of a fibroblast model, especially at high rates typical of AF.

In all diffuse fibrosis conditions, re-entry was unstable and self-terminated after several seconds (Figure 7B). In the ionic remodeling condition, known to lead to a large decrease of APD and the wavelength (Aslanidi et al., 2011; Colman et al., 2013), a single stable rotor persisted through the entire simulation (Figure 7A). Another pattern was observed in the high anisotropy condition linked with interstitial fibrosis. Re-entrant waves initiated due to the wave break at the CT and a relatively short wavelength were afterwards broken down into multiple secondary wavelets (Figure 7C). The breakdown mechanism was due to the presence of a higher number of sites with slow transverse conduction (similar to the CT). The general 3D atrial activation pattern of multiple wavelets propagating quasi-chaotically and independently was typical of AF.

However, stabilization of sustained activity in specific atrial locations was not observed in any of the uniform conditions considered. In the next section, such a stabilization of re-entrant wave is linked with patient-specific areas of heterogeneous patchy fibrosis in the 3D atria.

The segmentations of LGE MR images yielded different fibrosis distributions in 3 AF patients. The reconstructed fibrosis distributions varied in the extent and severity, as seen in Figure 8. All datasets showed areas of more severe, dense fibrosis (corresponding to the highest LGE MRI intensity) in the posterior LA (Figure 8), which is a common ablation target. The relative quantized distributions of fibrosis in the thresholded datasets showed a higher distribution in the lower fibrosis levels. Significantly lower amounts of fibrosis were observed in the higher levels (4–5) despite the histogram bias for level 5 fibrosis. In all the datasets good connectivity of fibrotic regions was observed, with the fibrosis levels changing gradually (Figure 8).

In 3D atrial model simulations, variable duration pinning to the fibrotic region was typically seen for all patient-specific cases (Figure 9). Permanent rotor stabilization in the BZ was seen in two cases, where rotor was pinned until the end of simulation. One simulation produced rotor movement directly around an inner area of dense fibrosis (Figure 9A) and another within an adjacent border zone of less severe patchy fibrosis (Figure 9B). In both cases, the rotor core was localized within regions of broadly defined patchy fibrosis (levels 3–4). Spatio-temporal dynamics of the rotors is illustrated in Figures 9C,D and Supplementary Video 1.

In patient-specific 3D atrial model simulations with variation of LGE MRI intensity threshold, different behaviors were observed (Figure 10). In simulations with the higher thresholds (lower fibrosis level) no stable rotor pinning was observed, with only short lived (less than 1 s) meandering re-entry in the LA fibrotic regions in two patients (Figures 10A,B). Simulations of the same two patient-specific atrial models with the lower threshold (higher fibrosis level) showed severe wave breakdown, with multiple re-entrant wavelets in the LA and to a lesser extent the RA. Interactions of the wavelets in fibrotic regions prevented the generation of a single rotor (Figures 10B,D). The complex electrical activity in both patient-specific 3D atria was characteristic of AF and sustained for the entire simulation (10 s).

The 3D atrial simulations performed for this section provide further evidence to the importance of using realistic fibrosis patterns to study the dynamics of AF, and elucidate the role of fibrosis in arrhythmogenesis due to the unique wave interactions for different fibrotic distributions.

A novel model for atrial fibroblast electrophysiology was developed based on experimental data by Wu et al. (2014) and applied to study the effects of M-F coupling with the atrial myocyte. Atrial fibroblasts in the model had higher (more positive) RMP compared to the previous models of ventricular fibroblasts by MacCannell et al. (2007) and hence a stronger effect on the excitability of atrial myocytes. This is in agreement with previous studies showing that the RMP of fibroblasts can have significant effects on AP properties of the myocyte (Jacquemet and Henriquez, 2008; Maleckar et al., 2009).

During the early phases of AP, fibroblasts in the M-F model act as a current sink effectively reducing the upstroke of the myocyte, with the AP amplitude (APA) reduced by 5 mV with the atrial fibroblast model (Figure 4). Similar behavior has been shown in other studies (Jacquemet and Henriquez, 2008; Maleckar et al., 2009), but with larger M-F coupling conductances. The APA reduction can slow down the propagation in atrial tissue due to the reduced electrotonic current transferred to neighboring myocytes. During the later phases of the AP, the fibroblasts in the M-F model act as a current source, prolonging the repolarisation phase (Figure 4) and increasing ERP (Figure 5E) in the myocyte.

The latter result is similar to the findings by Jacquemet and Henriquez (2008) and Maleckar et al. (2009) who both used a modified MacCannell et al. (2007) model for ventricular fibroblasts. However, the greatly prolonged ERP in the atrial myocyte coupled to ventricular fibroblasts (Figure 5E) is non-physiological, and also results in non-physiologically large reductions in atrial CV (Figure 5D). Coupling to the novel atrial fibroblast model produced more physiological values of atrial ERP and CV (Figure 5). A relatively small CV decrease in this case also had a different mechanism. The excitation threshold in the myocyte coupled to atrial fibroblasts was decreased (Figure 5C) due to increased RMP in the myocyte (which was caused by inflowing I_Gap from the atrial fibroblasts with relatively higher RMP). Counter-intuitively, the decreased threshold did not result in increased atrial CV–in fact, CV was significantly decreased (Figure 5D). This was due to the inflow of I_Gap into the myocyte leading to a partial inactivation of I_Na (Figure 5F).

Thus, the developed model for atrial fibroblast is more relevant to atrial electrophysiology (Table 1) than the existing ventricular fibroblast models, and also provides more relevant insights into the arrhythmogenic effects of M-F coupling in the atria. The latter include CV reductions that can be responsible for regional conduction blocks and the generation of re-entrant waves. The fibroblast model provides a valuable tool for further studying AF in 3D atrial models. The model could be further validated if novel data from human atrial fibroblasts become available.

The 3D atria modeling enables exploring the interplay between effects of fibrosis on atrial structure and function and the resulting AF scenarios. Three uniform effects of AF-induced remodeling in the entire 3D atrial model were considered: ionic remodeling, diffuse fibrosis (M-F coupling + decreased M-M coupling) and interstitial fibrosis (diffuse fibrosis + increased atrial tissue anisotropy). Simulations showed a range of re-entrant behaviors (Figure 7).

Similar to previous studies (Colman et al., 2013), the 3D model was first validated against atrial activation data in sinus rhythm, demonstrating a good agreement with clinical measurements by Lemery et al. (2007) Further simulations of ionic remodeling effects also agreed with previous results (Aslanidi et al., 2011; Colman et al., 2013), showing a stable rotor in the 3D atrial model in the case when M-M coupling was normal (Figure 7A). This was due to the remodeling greatly decreasing atrial APD, and therefore the wavelength.

In case of diffuse fibrosis, simulations produced unstable re-entry that self-terminated after several seconds. This was due to the fact that APD and ERP in this case was relatively high (Figures 5A,E, 7B), which resulted in a long wavelength comparable to the size of the atria. Re-entry was sustained longer in the model with coupling to ventricular fibroblasts compared to atrial fibroblasts (12 vs. 3 s), due to low CV (Figure 5D), and hence the wavelength (Figure 6), in the former case. These results also point to the sensitivity of 3D atrial behavior to the choice of a fibroblast model, which should be relevant to atrial electrophysiology.

In case of interstitial fibrosis, AF-induced increase in tissue anisotropy (Kawara et al., 2001; Koura et al., 2002; Verheule et al., 2013; Angel et al., 2015; Krul et al., 2015) resulted in multiple conduction blocks and wave breakdown throughout the 3D atria (Figure 7C). This was the only condition where the 3D atria model produced a sustained AF-like patterns with multiple re-entrant wavelets. Overall, changes of the 3D wave dynamics due to the increased anisotropy were more significant than the effects of M-F coupling or wavelength reduction due to ionic remodeling (Figure 7). This suggests that the structural effects of fibrosis, which develop in later stages in AF progression, may play the most important role in AF sustenance. However, stabilization of sustained activity in specific atrial locations (Haissaguerre et al., 2016), was not observed in any of the uniform cases, and can only be seen with heterogeneous distributions of fibrosis.

The 3D atrial simulations demonstrated unique wave interactions with different fibrotic distributions, and hence the importance of using patient-specific LGE MRI data to study AF.

The reconstructed LGE MRI patient datasets showed increased fibrosis distribution on the LA posterior wall (Figure 8). This region may be statistically prone to fibrotic infiltration and is often a site where complex fractionated electrograms are identified, which are common targets for ablation. The lower extent of fibrosis on the LA roof and floor may be due to the lower Z axis resolution of the MRI scan. Regardless of limitations of the imaging method, LGE MRI intensity has been strongly linked to fibrosis, and the reconstructed fibrosis distributions agree with existing knowledge from various modalities.

In 3D atrial simulations with patient-specific fibrosis distribution, a clear stabilization of rotors in regions of patchy fibrosis was observed in 2 out of 3 cases. One case produced a larger circus movement around a dense fibrosis region (Figure 9C) and another produced a smaller extent movement along the border of a dense region (Figure 9D). In both cases, rotors moved in slow-conducting BZ of patchy fibrosis surrounding inner regions of dense fibrosis. Note that fiber orientation can also be important in determining the rotor dynamics (Colman et al., 2014). However, in our simulations slow conduction in the fibrotic BZ appeared to be the main determinant of the rotor location–without patchy fibrosis, rotors meandered freely as seen in Figures 6, 7.

Further simulations also showed that patient-specific fibrosis distributions can induce wave breakdown (Figure 10) which is symptomatic of AF. This was only observed in simulation of high levels of fibrosis (obtained with a low LGE MRI intensity threshold) and resulted in permanent AF for the duration of the simulation (10 s). Smaller re-entrant wavelets were observed in these simulations moving chaotically throughout the fibrotic region. Contrarily, simulations of lower levels of fibrosis in the same patient (obtained by using a higher LGE MRI intensity threshold) produced very short duration re-entry that was neither stable nor resulted in wave breakdown. Note that the overall amount of fibrosis was not very different between the high and low thresholds. Rather, fibrosis was redistributed toward a higher proportion of severe fibrosis (levels 3–5) with an overall longer slow-conducting BZ. Hence, the rotor dynamics in AF may be determined primarily by the presence of such distributions.

Our results are in agreement with previous computational studies of atrial fibrosis: McDowell et al. (2015) and later Zahid et al. (2016) have shown that rotors in 3D atrial models can be localized near areas of fibrosis introduced based on LGE MRI data. Similar observations were made in cases of discrete fibrotic zones (Zahid et al., 2016), statistical distributions of fibrotic properties from dense fibrosis to healthy tissue (McDowell et al., 2015), and the respective continuous distributions in the current study. Hence, the mechanism of rotor stabilization proposed in our study can be independent on a specific way of modeling fibrosis. The accuracy of modeling fibrosis can be important for finding precise location of the rotors, and our approach of reconstructing continuous areas of patchy fibrosis may provide an advantage over models based on binary (Zahid et al., 2016) or statistical (McDowell et al., 2015) reconstructions.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.