A Three-Dimensional Human Atrial Model with Fiber Orientation. Electrograms and Arrhythmic Activation Patterns Relationship

The most common sustained cardiac arrhythmias in humans are atrial tachyarrhythmias, mainly atrial fibrillation. Areas of complex fractionated atrial electrograms and high dominant frequency have been proposed as critical regions for maintaining atrial fibrillation; however, there is a paucity of data on the relationship between the characteristics of electrograms and the propagation pattern underlying them. In this study, a realistic 3D computer model of the human atria has been developed to investigate this relationship. The model includes a realistic geometry with fiber orientation, anisotropic conductivity and electrophysiological heterogeneity. We simulated different tachyarrhythmic episodes applying both transient and continuous ectopic activity. Electrograms and their dominant frequency and organization index values were calculated over the entire atrial surface. Our simulations show electrograms with simple potentials, with little or no cycle length variations, narrow frequency peaks and high organization index values during stable and regular activity as the observed in atrial flutter, atrial tachycardia (except in areas of conduction block) and in areas closer to ectopic activity during focal atrial fibrillation. By contrast, cycle length variations and polymorphic electrograms with single, double and fragmented potentials were observed in areas of irregular and unstable activity during atrial fibrillation episodes. Our results also show: 1) electrograms with potentials without negative deflection related to spiral or curved wavefronts that pass over the recording point and move away, 2) potentials with a much greater proportion of positive deflection than negative in areas of wave collisions, 3) double potentials related with wave fragmentations or blocking lines and 4) fragmented electrograms associated with pivot points. Our model is the first human atrial model with realistic fiber orientation used to investigate the relationship between different atrial arrhythmic propagation patterns and the electrograms observed at more than 43000 points on the atrial surface.


Introduction
The most common sustained cardiac arrhythmias in humans are related to the atria. Different atrial arrhythmias, mainly atrial fibrillation (AF), often provoke disabling symptoms and severe complications such as heart failure and stroke [1]. There are several experimental observations regarding the important role of the anatomical structure and electrophysiological heterogeneity on atrial electrical activity both in physiological and pathological conditions [2][3][4]. Additionally, it has been shown that atrial tachyarrhythmias produce a set of changes in atrial properties that lead to their perpetuation. These changes, denoted atrial remodeling, include alterations in the electrical cellular activity and in the anatomical structure. They have been described in several animal models [5] and in humans [6,7]. Changes in electrical activity cause a decrease in refractoriness subsequent to a significant action potential duration (APD) shortening [6][7][8], which may help the initiation and maintenance of multiple reentrant waves, as suggested by experimental studies [5,6].
It is well known that atrial tachyarrhythmias can be caused by various mechanisms, including rapid local ectopic activity, singlecircuit reentry and multiple-circuit reentry [9][10][11][12]. It is important to know the fundamental mechanisms underlying atrial arrhythmias since they have implications in antiarrhythmic therapy. It is thought that different mechanisms lead to changes in the characteristics of spatiotemporal organization of atrial tachyarrhythmias. For the study of the spatiotemporal organization of atrial tachyarrhythmias, different signal analysis techniques are currently being used to analyze the electrograms (EGM) recorded at different points of the atrial surface. These include the analysis of EGM morphology [13][14][15], dominant frequency (DF) [16][17][18][19][20][21] and regularity [19] or organization index (OI) [22,23].
Experimental and clinical studies have shown various degrees of spatiotemporal organization during sustained AF [10,18,24]. Recent evidence from high-density mapping and spectral analysis has shown that AF is associated with rotors and regular repetitive activation in a portion of the atrium. In addition, areas of complex fractionated atrial electrograms (CFAE) and high DF have been proposed as critical regions for maintaining AF [10,24,25]. Thus, identification of areas with these features may allow targeted ablation procedures to improve success rates. The OI, defined by Everett et al [22] for the quantification of the organization of atrial activity from spectral analysis of EGM is another parameter that is used to identify critical areas for maintaining atrial arrhythmias. However, there is a paucity of data on the relationship between CFAE, DF and OI. Even more, despite the important experimental and clinical evidence obtained in recent years, the relationship between the different characteristics of EGM and the propagation pattern that underlies them are still far from being completely understood.
In this study, a realistic three-dimensional (3D) computer model of the human atria is used to investigate the relationship between different atrial arrhythmic propagation patterns and the EGM observed at different points on the atrial surface. Firstly, an anatomically accurate 3D model of human atria was developed. Our model includes a realistic geometry with fiber orientation, anisotropic conductivity and electrophysiological heterogeneity for different atrial tissues. Then, different tachyarrhythmic episodes in the electrically remodeled model of human atria were simulated and the EGM were analyzed to characterize the arrhythmic patterns. This model is one of the most complete 3D human atrial models ever developed and it provides a useful tool for investigating the complex phenomena involved in atrial tachyarrhythmias by simulating experimental and clinical situations difficult to perform in vitro or in vivo.

3D model of human atria
An anatomically realistic 3D model of human atria that includes fiber orientation was developed. Our model comprises the following main anatomical structures: left atrium (LA) and right atrium (RA), twenty pectinate muscles (PM) in the RA free wall, the fossa ovalis (FO) in the septum with its limbus, Bachmann's bundle (BB), the crista terminalis (CT), left and right appendages (LAPG and RAPG), left and right pulmonary veins (LPV and RPV), superior and inferior caval veins (SCV and ICV), the isthmus of RA, atrioventricular rings (AVR) and the coronary sinus (CS). The sinoatrial node (SAN) is situated near the ostium of the SCV. The CT layer originates next to the SAN and crosses the RA towards the ICV while the BB starts in a region next to the SCV and extends through the LA. The CS starts in the RA and extends along the posterior atrioventricular sulcus. The model also includes three different pathways for the inter-atrial conduction of electrical propagation: the BB, the limbus of the FO and discrete sites of the CS. Figure 1 shows the frontal ( Figure 1A) and dorsal ( Figure 1B) views of the model.
The original set of surfaces of the model was obtained from Harrild and Henriquez [26]. These surfaces were compared with data from the literature of normal human atrial structures and histological observations [2,27,28], and were modified. The dimensions of the different structures of the model are near the experimentally recorded average in human atria (see Table 1). The atrial wall is a monolayer surface, except the BB and the PM, which are solid structures.

Fibers orientation and anisotropy
One of the key characteristics of the atrial model is the incorporation of realistic fiber orientation based on histological observations [2,27,[29][30][31]. Our model was divided into 42 areas according to the orientation of the main muscle bundles (circular, longitudinal, transverse or oblique) in order to assign a realistic fiber direction to each region. In Figure 2, the main areas of the model and their fiber orientation can be seen. It shows circulating muscle bundles [2,27,31] around the CS and orifices of the pulmonary veins (PV), the SCV and ICV, the AVR, and both appendages (LAPG and RAPG). Fibers of BB, CT and PM have a longitudinal orientation. The free walls of both atria (RA and LA) are smooth and primarily have vertical fibers. The right free wall includes the CT, a fibrous bundle that runs vertically between both caval veins. The fibers between the superior and inferior PV are horizontally positioned, whereas the ostium of the PV contains a complex arrangement of vertical, horizontal and circular fibers [30]. This model is an improved version of the atrial models previously developed by our group [32][33][34].

Electrophysiological model
Nygren's model [35] of the human atrial action potential (AP) was implemented to reproduce the cellular electrical activity. In this study, the electrophysiological heterogeneity was introduced to reproduce APs in different parts of the atria [36]: CT, PM, APG, AVR and AWM (atrial working myocardium), which includes the remaining atrial structures. The maximum conductance values of I t , I Kr and I CaL were modified in Nygren's cellular model to obtain different AP models (see control values in Table 2).
The basic equation to calculate transmembrane voltage V m is: where C m is the specific membrane capacitance (50 pF), I ion is the total membrane ionic current, V m is the membrane potential and I stim is the stimulus current.
To reproduce atrial electrical remodeling, changes in the maximum conductance and kinetics of different ionic channels of human atrial cells observed in experimental studies of chronic AF [7] have been incorporated into the AP models. The following parameters were altered [37]: maximum conductance of I K1 was increased by 250% while the maximum conductance values of I CaL and I t were decreased by 74% and by 85%, respectively (see remodeled values in Table 2), the time constant of the fast inactivation of I CaL was increased by 62%, the activation curve of I t was shifted by +16 mV, and the inactivation curve of I Na was shifted by +1.6 mV. Figure 3 depicts the APs for the different atrial cellular models considered, under both physiological (control) and remodeling conditions ( Figure 3A and 3B, respectively). In these figures, we present the last AP obtained when a train of 10 stimuli at a basic cycle length of 1000 ms was applied. The corresponding APD 90 (APD to 90% of repolarization) values for the different atrial cells (both for control and remodeling conditions) are shown in Table 2. Under control conditions, APD 90 showed high values (ranged from 180 ms to 307 ms) in agreement with experimental data [7,36] and a great APD dispersion was observed (APD 90 max -APD 90 min = 127 ms). By contrast, under remodeling conditions, the APD decreased (ranged from 56 ms to 92 ms) and a smaller APD dispersion was observed (APD 90 max -APD 90 min = 36 ms). Figure 3C shows the APD 90 restitution curve for control and remodeled cells when the coupling interval (CI) between pulses is increased. Remodeling conditions not only induce a shortening in the APD 90 but also reduce the frequency dependent adaptation of APD 90 .

Action potential propagation and electrograms
In this study, the monodomain model of the electrical propagation of AP along the tissue model is described by the following reaction-diffusion equation [13,38,39]: where S v corresponds to the surface-to-volume ratio and D is the conductivity tensor. The monodomain equation was solved using a finite element method [40]. Our model includes regions with very high (CT), high (BB, limbus of the FO and PM), low (PV and isthmus of the RA) and very low conductivity (SAN). The AWM has mid-range conductivity (see areas with different colors in Figure 1).
A simplified two-dimensional model of 1506150 cells was used in order to select the values of tissue conductivity in the longitudinal and transversal direction of the fibers for the different regions. The conductivity values were fixed in order to obtain conduction velocities within the ranges reported in the literature [41][42][43][44][45], from 25 cm/s in very slow regions (SAN) to 143 cm/s in very fast regions (CT bundles) in normal (without remodeling) atria (see Table 3). Table 3 also shows that atrial remodeling provokes a reduction in the conduction velocities, between 9% and 24%. Anisotropic conduction was also assumed with a transversal to longitudinal ratio of conductivity of 1:9 for CT, which is in accordance with experimental studies [44] that have reported an anisotropy ratio in CT of 1:10. The isthmus of the RA and the SAN were set as isotropic, while an anisotropic ratio of 1:2 was used for BB, limbus of FO, PV and AWM [46]. Unipolar EGM at more than 43000 points of the atria surface (0.2 mm from the surface) were simulated. The extracellular potential (Fe) was computed using the large volume conductor approximation [38,47,48]: where ='V m is the spatial gradient of transmembrane potential V m , K is a constant that includes the ratio of intracellular and extracellular conductivities, r is the distance from the source point (x, y, z) to the measuring point (x', y', z') and dv is the differential volume. Surface pseudo-unipolar EGM at different points were visually inspected in order to seek their morphologies as either single, double or fractionated potentials. Double potentials were defined as EGM with two negative deflections, and CFAE were defined as EGM either exhibiting multiple (more than two) deflections or EGM with continuous (more than 50 ms) electrical activity without an isoelectric line [49]. EGM were digitally processed (every millisecond) with a 40-250 Hz band-pass filter, then rectified and low-pass filtered at 20 Hz [50]. Next, spectral analysis of the signals was performed with a fast Fourier transform obtaining a spectral resolution of 0.12 Hz. The DF defined as the frequency corresponding to the highest peak of the power spectrum was calculated. In addition, to measure the periodicity of the signal and the variability of the frequency in the spectrum, the OI was also calculated. The spectral power of the DF and its three harmonic peaks were calculated by computing the area under their peaks over a 0.75 Hz window. Finally, the OI was obtained as the ratio of this spectral power to the total power of the spectrum [23,50].

Numerical and computational methods
A hexahedral mesh was built from the three-dimensional anatomical model using Femap from Siemens PLM software. This mesh includes 52906 elements and 100554 nodes with a  spatial resolution ranging from 300 to 700 mm, the average element size being approximately 530 mm. Similar spatial resolutions have been used in previous computational studies of atrial propagation [26,51]. Equation (2) was numerically solved using the software EMOS [52,53]. EMOS is a parallel code (mpibased) that implements the finite element method (FEM) [40] and the Operator Splitting [54] for solving the monodomain model (Eq. 2). The time step was fixed to 0.02 ms. Simulation of 10 seconds of atrial activity took 14 hours on a computing node with two 6-core Intel Xeon X5650 clocked at 2.66 GHz and 48GB DDR3 RAM.

Simulation protocols
SAN and ectopic stimulation were simulated by applying rectangular current pulses of 6 ms in duration and 30 mA in amplitude on an area of approximately 10 mm 2 . Atrial tachyarrhythmias were generated by the S1-S2 protocol, which consisted of a train of stimuli with a basic cycle length of 1000 ms applied during 10 seconds in the SAN area to simulate the sinus rhythm (S1). After the last beat of the SAN stimulation, an ectopic focus (S2) was delivered in six different locations in the atria: CT, SCV, isthmus of RA (near the CS), and ostium of the LPV and of the RPV. Two different ectopic foci S2 were used in our study: 1) a train of six stimuli at cycle length (CL) of 130 ms; and 2) a continuous focus also at CL of 130 ms.

Atrial action potential propagation
An atrial propagation pattern in sinus rhythm was simulated in the model by applying a periodic stimulus (10 beats) at a basic cycle length of 1000 ms in the anatomic location of the SAN. Figure 4 shows different snapshots of the propagation of the last beat applied for both physiological ( Figure 4A) and remodeling ( Figure 4B) conditions. A delayed activation front in the remodeled atria can be observed. In Figure 4, it is possible to note how the stimulus applied in the SAN region caused the initiation of a propagating wavefront that quickly spread to the ICV favored by the high conductivity of the CT arch. The anisotropy of the CT caused an almost triangular wavefront (see snapshots at 36 ms in Figures 4A and 4B). The depolarizing wave reached the BB after 22 ms and propagated to the anterior septal portion of the LA through interatrial BB, inducing the first activation of the LA at 46 ms (in normal atria) and 54 ms (in remodeled atria) after the SAN activation. Thereafter, at 53 ms and 62 ms from the SAN activation for normal and remodeled atria, respectively, the interatrial connection at the limbus of the FO contributed to left atrial septal activation. The activation wave also propagated to LA through the third interatrial connection, the RA-CS-LA pathway, at 87 ms and 98 ms, respectively, and the first LA inferior wall activation was observed at 117 ms and 132 ms for normal and remodeled atria, respectively.
In Figure 4 (see snapshots at 70 and 101 ms) the delay in the propagation induced by remodeling conditions can be observed. The entire atrial depolarization occurred at 120 ms in normal conditions and at 138 ms in remodeled atria (see snapshot at 138 ms in Figure 4B), finishing in the distal left APG. Moreover, the whole atrium repolarized faster in remodeling conditions than  under physiological conditions, at 234 ms and 432 ms after SAN activation, respectively. After the 10 th stimulus was applied in the SAN zone, two different ectopic foci (S2) were applied at six different points of the atria (see methods): 1) a burst of six beats, and 2) a continuous ectopic focus, both with a CL of 130 ms. When we applied the burst of six beats, we observed either atrial flutter (AFL) or reentrant tachycardia in five of the six locations tested (isthmus of RA near the CS, SCV, CT, ostium of the LPV and center of the LA posterior wall), and AF when the ectopic focus was applied in the ostium of the RPV. Different AF patterns were also observed when a continuous ectopic focus was applied in the six locations. Figure 5A depicts an example of AFL obtained by applying the burst of six ectopic beats in the isthmus of RA (near the CS), 231 ms after the last sinus stimulus. When the burst was applied, a conduction block was observed in the CS area which induced a macroreentry rotating counterclockwise around the tricuspid annulus (see the arrow in snapshot at 4100 ms in Figure 5A), that was maintained during the whole simulation (<8 seconds). The activation pattern of this macroreentry resembles an AFL (see video S1 in the Supporting Information). The direction of the reentrant circuit can also be observed when comparing the APs of the sites 1, 2 and 3 in Figure 5A (see dotted arrows). The macrorentrant wavefront depolarized the rest of the atria, including the posterior wall of the LA (see site 4 in Figure 5A) with a 1:1 activation pattern. The CL of this arrhythmic pattern was almost constant (<200 ms) in the whole atria, a characteristic of AFL. The EGM in different points of the atria (see sites 2 and 4 in Figure 5A) only displayed single potentials, showing a stable and regular atrial activation, which is also characteristic of AFL. This regularity is reflected as a single narrow DF peak of 5.0 Hz and high OI values, close to unity, in sites 2 and 4. Similar DF and OI values were observed in the entire atrial tissue. Figure 5B shows an example of reentrant tachycardia (see snapshots at 3760 ms and 3840 ms) obtained by applying a burst of six ectopic beats approximately in the middle of the CT, 160 ms after the last sinus stimulus. In this case, the ectopic focus initiated a reentry anchored to the RA, that propagated between the center and the top areas of the CT (see video S2 in the Supporting Information). This propagation pattern was favored by the greater anisotropy, fiber orientation and heterogeneity in the CT region. In fact, the CT area acts as a blocking line (see the dashed line in snapshot at 3840 ms in Figure 5B), maintaining the reentry activity and keeping the rest of the tissue in a reentrant tachycardia pattern. The blocking line is a source of electrotonic APs where functional conduction blocks are observed. The APs during tachycardia at five sites are also shown in Figure 5B. It is important to note that in the CT region, where the reentry turns, two impulses arrive to site 2 from sites 1 and 3 in each cycle. Site 3 shows electrotonic APs when the reentry passes through the region of site 1 and 2, but it is blocked before arriving to site 3. Likewise, site 1 shows electrotonic APs when the reentry turns and passes through site 3 and 2, but it is blocked before arriving to site 1. This reentrant activity spreads to the rest of the atria (see APs at sites 4 and 5 in Figure 5B).

Atrial flutter and reentrant tachycardia
The EGM observed in areas of conduction-block (see EGM at site 2) are characterized by double potentials, one short negative component associated with the electrotonic APs elicited when the reentry is passing at a close distance, and another larger and steeper negative component generated when the reentry passing through the electrode area. The EGM at site 5 displays single potentials with only negative deflections of short-time when the curved front passes on this site (see EGM at site 5 in Figure 5B). The EGM observed in the rest of the atria (out of the blocking line area) showed single potentials.
Only slight variations in the CL of APs were observed between different points, which suggest a periodic and regular reentrant activity. The entire atrial tissue presents similar DF (around  5.5 Hz) and high OI values (see sites 2 and 5 in Figure 5B, with values of 0.87 and 0.99, respectively). Reentrant tachycardias were also observed when the ectopic focus was applied in SCV, LPV and the LA posterior wall (data not shown).
Atrial fibrillation Figure 6 shows two AF episodes triggered by an ectopic focus placed in the ostium of the RPV. In the first case, a train of six ectopic beats (CL = 130 ms) provoked an AF episode that was maintained by different reentrant mechanisms (see snapshots in Figure 6A). Indeed, rotors, figure of eight reentries, macrorreentries, fragmentations and wavefront collisions were observed in  both atria (see video S3 in the Supporting Information). This irregular activity caused the appearance of electrotonic APs and irregular APs in different areas and at different times (see APs at sites 2 and 3 in Figure 6A).
The EGM primarily show irregular and polymorphic potentials, but double potentials and CFAE were also observed in minor areas. The first potential of the EGM calculated at site 1 (see + in EGM at site 1 in Figure 6A) only shows a long-lasting negative deflection when a spiral wavefront passes over this point and away from it (see the black narrow in snapshot at 2160 ms in Figure 6A). At this site, potentials with positive deflection greater than negative ones (see N in EGM at site 1 in Figure 6A) were observed when wave collisions occurred on the recording point (see snapshot at 2355 ms in Figure 6A).
The EGM calculated at site 2 shows potentials with a small additional positive deflection generated when the wavefront is fragmented (see & in EGM at site 2 in Figure 6A). Then one of the two fronts is shifted to the top of the CT returning to depolarize the site 2 (see the arrow in snapshot at 3220 ms in Figure 6A), which gives rise to potentials with double positive deflections.
CFAE at site 3 were observed when the tip of the rotor turned on this point (see the arrow in snapshot at 8675 ms in Figure 6A). The wavefront surrounds the pivot point, without depolarizing it completely, which results in multiple low amplitude deflections in the EGM (see * in EGM at site 3 in Figure 6A).
Variations in the CL of the APs were observed between different points of the atria, which exceeded 20 ms in many cases. This fact is reflected in the appearance of multiple frequency peaks and low OI values in the spectral analysis. Heterogeneity in DF values and low OI values were observed in the entire atrial tissue (see spectral analysis of EGM at sites 1, 2 and 3 in Figure 6A). Figure 6B shows an example of a focal AF episode triggered by a continuous ectopic focus applied in the ostium of the RPV. This ectopic activity also induced multiple reentrant waves that collided and fragmented with each other (see snapshots in Figure 6B and video S4 in the Supporting Information). Rotors with irregular trajectories that collided and fragmented creating new activation fronts were mainly observed in the RA. Although LA also presented fragmentation and collisions in the superior and anterior walls, the posterior wall near RPV was activated at a higher frequency by the wavefront generated by the ectopic focus.
In site 1 of Figure 6B, a reentry turns around a functional block for several milliseconds. Electrotonic APs reflected as deviations from the baseline at the EGM, CFAE, (see x in EGM at site 1 in Figure 6B) were observed when the reentry turned very close to site 1 (see snapshots at 4545 ms and 4650 ms in Figure 6B). In other instants, the spiral wavefront moved over the site 1 generating potentials with a single negative deflection in the EGM (snapshots not shown).
The APs and the EGM calculated in site 2, near the ectopic focus, show a rapid and regular atrial activation with similar CL values. The EGM has potentials with only negative deflections when the curved wavefront passes through this point (see snapshot at 7625 ms and % in EGM at site 2 in Figure 6B). The stable and regular activity in this site is reflected in the spectral analysis as a single narrow DF peak of 7.7 Hz (value of the ectopic focus frequency) and high OI value close to unity.
Site 3 in Figure 6B shows APs and EGM for a point on the blocking line, which induced reentry anchored to the CT in certain instants (see the dashed line in snapshot at 7625 ms in Figure 6B). The reentry propagates between the center and the top areas of the CT, similarly to that observed during reentrant tachycardia in Figure 5B. Site 3 is depolarized by the wavefront and it is again activated, before repolarizing, by the front turning on line blocking. This behavior is reflected in EGM with double potentials (see # in EGM at site 3 in Figure 6B).
The only zone of the atria that shows a high frequency of activation (similar to the focus) and high OI values is the region of the posterior wall of LA near to the ectopic focus. The remaining tissue is unable to depolarize at this high frequency and it shows irregular CL, which is reflected in the spectral analysis with multiple frequency peaks and lower DF and OI values (DF of 5.6 Hz and OI of 0.36, and DF of 5.2 Hz and OI of 0.48, at sites 1 and 2, respectively, in Figure 6B). Every AF rhythms simulated when the continuous focus was applied in the other five locations (not shown) displayed stable and regular activity near the area of the ectopic focus and high DF values similar to the focus frequency. In the rest of the tissue, multiple reentrant waves generating irregular activity were observed.

Validation of 3D model of human atria with complete fiber orientation
A 3D model of the human atria with realistic anatomy and fiber orientation was developed. Our model also includes real anisotropy characteristics and electrophysiological heterogeneities in the main atrial structures. To our knowledge, this model is one of the most complete 3D human atrial models ever developed.
Our model is based on the surfaces provided by Harrild and Henriquez [26] that were improved using detailed anatomical and histological observations [2,[27][28][29][30][31][55][56][57][58]. We modified the geometry and dimensions of the PV, APG, PM and AVR rings, and the CS was created, in accordance with human data [27,28,[55][56][57][58]. Specifically, the diameters antero-posterior and medial-lateral of our model are 4.5 cm and 4.3 cm for the RA, respectively, and 4.3 cm and 3.7 cm for the LA, respectively. These values fall within the ranges reported by Cohen et al [28] in humans (2.8-5. Our model also includes anisotropy and fiber orientation. It has been observed that atrial anatomy, anisotropy and fiber orientation are critically important in determining the spread and direction of the activating wave front [2,3,4,29,31]. There are a number of well-defined bundles in the atrial muscle where the AP propagates faster [29,[41][42][43]. These bundles include the CT, BB, the PM and the limbus of the FO. All of them are considered in our model. The atrial model was divided into 42 regions and a detailed fiber orientation was included in each region based on histological observations [2,[27][28][29][30][31]59,60]. This model goes beyond previous models, as it not only considers a detailed atrial anatomy, regional heterogeneity and anisotropy, but also accurate fiber orientation for the whole atria. Jacquement et al [13] developed an atrial model of simple geometry with a coarse fiber structure, which was manually introduced. Seeman et al [61] and recently Aslanidi et al [62] included only the fiber orientation of the main bundles (BB, PM and CT) in their human atrial model. Krogh-Madsen et al implemented the Harrild and Herriquez human atrial model [26] without structural modifications, anisotropy and fiber orientation to simulate arrhythmic episodes [63]. Ho et al [31] argued that the heterogeneous myoarchitecture of RA and LA and the inter-atrial septum must be considered in computer models that seek to investigate mechanisms of atrial arrhythmias, as our model does. An image-based anatomical model of the sheep atrial fiber orientation has been recently presented [64]. The model reproduces the whole atria with highly detailed myofiber architecture while our model simplifies the complex myocardial structure by dividing the whole atria into 42 regions. Fiber orientation in both models has some similarities, for example fibers at the base of SVC and around PV have circumferential orientation or highly organized tracts in the CT and PM. However, it is not clear how similar are fiber architecture in the whole atria between sheets and human.
Our model also includes a realistic RA-CS-LA interatrial connection; where the CS connects electrically the RA and the LA through some fiber tracks. Different anatomic studies [65] have shown that the CS musculature is continuous with LA myocardium in the proximal portion of the CS and with RA myocardium at the CS orifice. Therefore, the CS musculature may form an electrical connection between the RA and the LA.
The electrophysiological heterogeneity also plays an important role during atrial conduction. In order to reproduce the AP morphologies in different parts of the atria, the human atrial model developed by Nygren was modified based on the electrophysiological heterogeneity observed in canine atria [36] as previously suggested [61]. The effect of electrical remodeling was introduced by altering different currents (I K1 , I CaL , I t and I Na ) in accordance with experimental studies of chronic AF [7]. In this study, under normal conditions, AP showed morphological differences in the different electrophysiological zones, consistent with those reported experimentally [36]. It is noteworthy that, under electrical remodeling conditions, the APs showed similar morphologies for different atrial zones with only slight APD variations between them.
Atrial action potential propagation. Our model reproduces the normal atrial activation originated in the SAN, both under physiological and remodeling conditions. Under physiological conditions, the activation initiated in the SAN region reached the BB after 22 ms, which follows experimental studies that have reported the first activation of BB at approximately 19 ms [59]. The propagating wavefront quickly spreads to the ICV, becoming almost triangular, favored by the high conductivity of the longitudinal fiber orientation of the CT. These characteristics of the wavefronts have been reported for normal sinus rhythm in humans [43].
In our simulation study, the RA free wall was completely depolarized after 60 ms. This result is in agreement with experimental observations of the complete RA wall activation at 57 ms [43]. Interestingly, the LA was also activated through the limbus of FO at 53 ms, and through the RA-CS-LA connection at 117 ms, which is consistent with recent studies that have suggested that the activation of the posterior epicardial LA takes place at 54610 ms [66] and that a third interatrial electrical connection exists in the region of the CS [65,66].
The complete LA activation was observed in the distal LAPG at 120 ms, which is also in agreement with experimental results. Canavan et al [67] showed that the last activation in the atrial tissue occurs just before 120 ms. Lemery et al [66] reported the latest LA activation at 116618 ms. In addition, it has also been documented that atrial activation ends at the LAPG [68].
Under remodeling conditions, the latest atrial activation also occurs in our model in the distal LAPG, but after 138 ms, showing a conduction delay. The delay is due to an overall reduction of the conduction velocity that changes from 9% to 24% depending on the zone [7,69]. Experimental studies in dogs [70] observed a mean reduction of 25% in conduction velocity values measured in remodeling conditions. Even more, under remodeling conditions a faster repolarization of the whole atria was observed due to a reduced APD.
Simulated tachyarrhythmic episodes resulted from the application of ectopic foci with both transient and continuous activity. It has been demonstrated that extrastimuli, either both transitory or continuous, can act as triggers and, in some cases, they may be responsible for the initiation and maintenance of AF episodes [9,71]. In our study, the foci were applied in six different locations in the atria: CT, SCV, isthmus of RA (near the CS), center of the LA posterior wall and ostium of the LPV and RPV. Although tachyarrhythmic episodes are initiated by focal triggers most commonly localized in areas near the PV [9,71], ectopic foci have also been recorded in other regions. In an extensive cases report [72] the ectopic foci were located in the posterior wall of the LA including VP (38.3%), in the SCV (37%), in the CT (3.7%) in the ligament of Marshall (8.2%), in the CS (1.4%) and in the interatrial septum (1.4%).
In our simulations, one AFL, four reentrant tachycardias and one AF episode were obtained by applying transient ectopic foci, but only AF episodes were observed when a continuous ectopic focus was applied in the same locations.

Flutter and reentrant tachycardia
The application of the train of six ectopic beats in the isthmus of RA (near the CS) triggered an episode of AFL. It was maintained by a macroreentry located in the RA, turning counterclockwise around the tricuspid annulus. In patients, it has been observed that a high percentage (85%) of AFL is due to reentrant excitation travelling around the tricuspid valve ring [73]. Notably, our model can reproduce this kind of AFL episode, where the circular fiber orientation around the tricuspid annulus (as in the real atrial musculature) plays an important role in the AFL reentrant circuit. The macroreentrant circuit activated the rest of the atria with a 1:1 pattern of activation and functional block was not found in the CT. This supports similar observations reached in a previous analysis of human AFL [74], where the atrial activation in the right lateral wall was uniform, without blocking in the region of the CT. In our simulation of AFL, no variation was observed in the CL (<200 ms) in the different regions of the atria. This result is also consistent with the mean values of CL (238 ms and 245 ms) obtained in experimental and clinical studies [75,76].
The EGM recorded in different points of the atria presented high organization, with uniform and regular single potentials, which reflected in the spectral analysis as a single narrow DF peak. This is in agreement with EGM recorded in vivo [77] in both LA and RA which also showed uniform and regular potentials. In our model, the spectral analysis showed DF of 5 Hz during AFL, which falls within the range reported in humans (4.2-5.8 Hz) [78]. Additionally, EGM also showed high OI values according to the high regularity in the atrial activation [23].
Application of the ectopic focus in the CT, SCV, ostium of LPV and the center of the LA posterior wall triggered reentrant tachycardia episodes. When the ectopic focus was applied in the CT, the reentrant tachycardia was maintained by a reentry anchored to the CT. The reentry moved between the central and the superior part of the CT, showing a block line in this region. Experimental studies have shown a great number of the reentrant tachycardias in absence of structural heart diseases occurring along the CT [60,79], due to the anisotropy in this region that favors the role of CT as a natural barrier to the atrial activation. Even more, it has also been demonstrated in experimental and clinical studies that the CT is an anatomical substrate underlying atrial arrhythmias [74,80].
The reentrant tachycardia produced EGM with mostly single, uniform and regular potentials, linked to stable and regular atrial activation. This organized activity is reflected in the spectral analysis with a single narrow DF peak (consistent with the CL of atrial tachycardia) and a high OI value in most areas of the atria. A 1:1 uniform activation pattern across both atria and a high regularity in the atrial activation were observed. Experimental and clinical studies [77,79] have shown stable and regular atrial activation with mostly single and regular EGM during reentrant tachycardias. Everett et al [23] have linked high OI values with high regularity in atrial activation. It is noteworthy that double potentials were observed in the EGM along the CT when it acts as a blocking line. Interestingly, our results also confirm that double potentials (or at least part of them) could be related to electrotonic potentials provoked by sequential activation on both sides of the blocking line. In an experimental study on atrial EGM during AF, Konings et al [15] registered double potentials along the lines of conduction block, similar to the results of our simulations.

Atrial fibrillation
When a train of six ectopic beats in the ostium of the RVP was applied, an episode of AF was initiated ( Figure 6A). Multiple reentrant waves, rotors, fragmented fronts and collisions were observed, which are characteristic of fibrillatory activity [10,15,77]. CL variations and EGM with single potentials double potentials and CFAE were observed. Experimental and clinical studies of AF [16,17,77] have recorded intracardiac polymorphic EGM, and irregular and disorganized activation patterns in areas with fibrillatory conduction.
In our study, EGM showed only a long-lasting negative deflection when a spiral wavefront passed over the recording point, and potentials with positive deflection greater than the negative ones when wave collisions occurred. These results are consistent with previous studies [13,81] which have shown that the EGM morphology is related to the shape of the wavefront and to its curvature, suggesting that single potentials may have different morphologies. Our results confirm the results obtained by Jacquemet et al [13] in simulation studies using a model with simplified fiber structure. He observed potentials with a positive deflection much lower than the negative deflection during either curved or spiral fronts. He also obtained potentials with a much larger positive deflection than negative deflection during a collision between two wavefronts. However, Konings et al [15] reported short-double potentials along either side of the line of collision, unlike observations in our study. The differences could be due to the collisions between several wavefronts (.2) observed by Konings et al, whereas in our study the collisions occur mainly between two wavefronts.
In our study, CFAE were observed when the tip of the rotors turned on the recording point (pivot point). During the last decades experimental and clinical studies have demonstrated that the maintenance of the AF in many cases depends on these rotors [17,82,83]. Recent observational studies have demonstrated that substrates serving as ''AF perpetuators'' can be identified by searching for areas that have CFAE, but the underlying etiology of CFAE and their relationship with rotors has not yet been elucidated [84]. Konings et al [15] recorded CFAE during AF in humans, suggesting that CFAE can indicate pivot points, slow conduction and complex reentrant patterns. Zlochiver et al [14] demonstrated that rotor meandering might also underlie, at least in part, the CFAE that they observed close to the driver. Our results reveal a direct relationship between the tip of the rotor (pivot point) and the CFAE.
The spectral analyses display multiple frequency peaks around the DF peaks and lower OI values than those observed during either reentrant tachycardia or AFL. These results are in agreement with experimental studies [15][16][17]77], reporting that spectral analysis of the EGM recorded during AF episodes shows multiple frequency peaks in areas with irregular and unstable activation. Indeed, conduction blocks and wave collisions increase the irregularity and variability in the frequency [15,16,82,85]. Everett et al [23] associated low OI values with the presence of EGM with double and fragmented potentials in canine models of AF.
In our study, the focal AF triggered by an ectopic focus with continuous activity between RPV showed substantial differences from the previous AF. Nearby zones and the location of the ectopic focus showed APs and EGM with rapid and regular activation with single potentials and similar CL values between them. This electrical activity gives rise to a single frequency peak with the highest DF value, similar to the frequency of the ectopic focus and with the highest OI values, in accordance with experimental observations. Lin et al [16] have reported that EGM recorded near the ectopic sources show rapid, regular and stable activity. The frequency spectrum obtained from sources maintaining AF usually shows a single frequency peak with a narrow morphological base [85]. Moreover, Takahashi et al [50], in a study with 25 patients, associated high OI values with sources that maintain the arrhythmia. Due to the high frequency of the ectopic focus, the rest of the atria cannot follow a 1:1 activation, therefore, the conduction block results in a reduction of DF in remote areas of the focal source, as suggested by different authors [16,19,20]. Additionally, several studies in which AF is maintained by a focal source in the LA have reported higher DF values in the LA than in the RA [16,18,25].
Our results show EGM with potentials that display a pronounced negative deflection when the curved front passes on the recorded points. Jacquemet et al [13] also obtained potentials with a negative deflection much larger than the positive deflection in the presence of curved fronts.
In areas of the atria away from the ectopic focus, multiple reentrant waves that collided and fragmented were observed, which produced irregular and polymorphic EGM with single, double and CFAE; a hallmark of fibrillatory conduction [10,15,77]. Indeed, double potentials are related with wave fragmentations and conduction block and CFAE are related with pivot points, as previously mentioned. Finally, CL variations are reflected in the spectrum as multiple frequency peaks and low OI values, indicating an irregular conduction pattern with high frequency variability.

Limitations of the study
Some of the limitations of this study include: The Nygren model, despite being a fairly complex and detailed model, has some limitations. It lacks a detailed dynamic intracellular Ca 2+ cycling, which can play a relevant role in some cardiac arrhythmias. Several electrophysiologically detailed human atrial cell models have been developed. Although all of them resemble the AP of human atria, they have different morphologies (see [86] for review). We used the modifications on I K1 , I CaL and I t reported by Zhang et al [37] to obtain electrically remodeled atrial cells models. In that work, the authors modified two atrial cellular models (Nygren et al [35] and Courtemanche et al [87]) and observed that the remodeling conditions (called in this work AF-1) provoked similar APD 90 reductions (62% versus 68%). In another study, Sánchez et al [88] compared the results of two atrial cell models (Courtemanche et al [87] and Maleckar et al [89]) using slightly different remodeling conditions. They observed that the APD 90 reduction was also similar in both models. Even more, although the restitution curve in control presented some differences between the cellular models, under remodeling conditions both restitution curves flattened and were qualitatively very similar. They concluded that their study of AF dynamics was independent of the model used. The aforementioned results suggest that the main findings of our work are also independent of the cellular model used.
The remodeling conditions used in our study are based on the experimental work of Bosh et al [7]. However, another experimental study [6] has reported shorter reductions of APD during AF (around 50%) and several computational studies [37,88,90] have introduced the effect of electrical remodeling mainly using quantitatively different modifications in the same channels. These different remodeling conditions may produce different patterns of atrial arrhythmias and requires further investigation. Another limitation is that the heterogeneity of our model is based on canine atria data [36]. New experimental data about the effect of atrial electrical remodeling on the anatomical structures of human atria would allow more realistic simulations of AF dynamic under chronic AF. Although the used remodeling conditions can reproduce the action potential phenotype observed in patients with permanent AF, electrical remodeling is not the only process accompanying permanent AF. Indeed, many of these patients have significant structural remodeling with fibrosis, which contributes to short APD [91,92] and increases the complexity of the arrhythmia [93]. Our 3D anatomical model of human atria does not take into account the real thickness of the atrial walls. Several studies have shown very complex fiber structures in areas of human atria [2,27,30,31]. Our results were obtained using a specific virtual atria model. Although our model includes a great number of anatomical and morphological details, it corresponds to a particular set of parameters (electrophysiology, anatomy, fiber direction, anisotropy and heterogeneity, among others). For example, our model includes four pulmonary veins, which is the most common anatomical structure in humans though other patterns have been observed [94]. In addition, although there are also inter-subject differences in fiber orientation, we have tried to model the most common fiber orientation observed experimentally for the different parts of the atrial model. The influence of aging on the anatomy of the atria is also relevant [95]. Therefore, more comprehensive results could be obtained when considering inter-subject variability. Indeed, Romero et al [96] analyzed the impact of the ionic current variability in line with biological intersubject differences, on the APD 90 in a model of human ventricular cells. Recently, a semi-automatic method to incorporate atrial anisotropy, heterogeneities and fiber orientation into patient specific models has been developed [97]. Although the fiber orientation was also based on data from literature and from different human atria preparations, the specific anatomy of patients is reproduced in 3D models using data from CT and MRI. This method could be used to further extend this study for personalized anatomical models.
Finally, the mesh used in the present study is relatively coarse. It comprises of 100554 nodes with a spatial resolution ranging from 300 to 700 mm (530 mm in average). Although this discretization is similar to the used in previous simulation studies of AF [26,51] recent works have used higher resolutions. Aslanidi et al [62] used a spatial resolution of 33063306300 mm 3 and showed that decreasing the space step form 330 mm to 250 mm resulted in very small (,3%) changes of the conduction velocity of AP. Although a higher spatial resolution could slightly affect to AF dynamics observed in our study, we would expect similar patterns of activation. Moreover, these changes would not affect the relationship between atrial activation patterns and simulated EGMs.

Conclusions
We have developed a realistic 3D model of human atria with a very detailed fiber orientation, real anisotropy and electrophysiological properties. Additionally, it includes all main interatrial connections. To our knowledge, this model is one of the most complete human atrial models developed.
In our study using an electrically remodeled atria model, a transient focus induced atrial flutter, reentrant tachycardias and AF; whereas a continuous ectopic focus provoked AF in all of the tested sites.
Our model can reproduce the stability and regularity that cause arrhythmias as AFL, reentrant tachycardia and areas with focal activity. This electrical activity is reflected in EGM with single potentials and spectra with a narrow frequency peak and high OI values. Our model can also reproduce characteristic patterns of AF, which were maintained by multiple reentrant waves, rotors, fragmentation and wave collision, which are reflected in polymorphic EGM and spectra with multiple frequency peaks and low OI values. Interestingly, EGM potentials with only a negative deflection are related to spiral and curved wavefronts that pass and move away. In addition, potentials with a much greater positive deflection are related with wave collisions. Moreover, double potentials are related with either wave fragmentations or blocking lines, while CFAE are related to pivot points.
Finally, this is the first work that uses a 3D human atrial to investigate the relationship between different atrial arrhythmic propagation patterns and the EGM observed at more than 43000 points on the atrial surface.

Supporting Information
Video S1 Atrial Flutter (supplements Fig. 5A in the manuscript). Simulation of an AFL episode, triggered by six ectopic beats applied at the isthmus of the RA (near the CS), and maintained by a macroreentry rotating counterclockwise around the tricuspid annulus. (AVI) Video S2 Reentrant tachycardia (supplements Fig. 5B in the manuscript). Simulation of reentrant tachycardia episode, triggered by six ectopic beats applied at the middle of the CT, and maintained by a reentry anchored to the RA, moving between the center and the top areas of the CT. (AVI) Video S3 Atrial fibrillation induced by a transitory focus (supplements Fig. 6A in the manuscript). Simulation of AF episode, triggered by six ectopic beats applied in the ostium of the RPV, and maintained by multiple reentrant waves. Rotors, figure of eight reentries, macrorreentries, fragmentations and wavefront collisions can be observed in both atria.

(AVI)
Video S4 Atrial fibrillation induced by a continuous focus (supplements Fig. 6B in the manuscript). Simulation of focal AF episode, triggered by a continuous ectopic focus applied in the ostium of the RPV. Multiple reentrant waves are observed in the RA and in the superior and anterior walls of LA. Posterior wall near to RPV show stable activity at high frequency due to the wavefront generated by the ectopic focus. (AVI)