Chronotropic Modulation of the Source-Sink Relationship of Sinoatrial-Atrial Impulse Conduction and Its Significance to Initiation of AF: A One-Dimensional Model Study

Initiation and maintenance of atrial fibrillation (AF) is often associated with pharmacologically or pathologically induced bradycardic states. Even drugs specifically developed in order to counteract cardiac arrhythmias often combine their action with bradycardia and, in turn, with development of AF, via still largely unknown mechanisms. This study aims to simulate action potential (AP) conduction between sinoatrial node (SAN) and atrial cells, either arranged in cell pairs or in a one-dimensional strand, where the relative amount of SAN membrane is made varying, in turn, with junctional resistance. The source-sink relationship between the two membrane types is studied in control conditions and under different simulated chronotropic interventions, in order to define a safety factor for pacemaker-to-atrial AP conduction (SASF) for each treatment. Whereas antiarrhythmic-like interventions which involve downregulation of calcium channels or of calcium handling decrease SASF, the simulation of Ivabradine administration does so to a lesser extent. Particularly interesting is the increase of SASF observed when downregulation G Kr, which simulates the administration of class III antiarrhythmic agents and is likely sustained by an increase in I CaL. Also, the increase in SASF is accompanied by a decreased conduction delay and a better entrainment of repolarization, which is significant to anti-AF strategies.


Introduction
Atrial fibrillation (AF) is the most common cardiac arrhythmia, characterized by high morbidity and mortality; the mechanisms underlying its initiation, known to be complex and multifactorial, are still largely unexplained [1]. The most generally recognized causes of initiation and maintenance of AF are conduction abnormalities along interatrial accessory pathways [2,3], abnormal interaction between sinoatrial node (SAN) cells and cardiac tissue nearby the insertion of pulmonary veins [4], increased fibrotic deposition within atrial tissue [5], electrical and calcium handling remodeling per se [6,7] or secondary to atrial tachycardia (AT) [8,9], and altered connexins ratio (Cx40/Cx43) leading to heterogeneity of conduction velocity [10,11]. This complex scenario basically underlies two main pro-AF mechanisms, that is, triggered activity and reentry [12]. Furthermore, the crucial role of bradycardia in initiating and maintaining AF is well documented, either when associated with pathological conditions, like the sick sinus syndrome [13,14], or brought about experimentally by cholinergic hyperactivation [15,16]. Importantly, SAN spontaneous activity plays an active role in controlling atrial arrhythmias by its ability to terminate or convert atrial flutter to AF during cholinergic withdrawal [17].
A factor that makes the treatment of AF particularly complex is the paradoxical proarrhythmic effect of some of the most commonly used antiarrhythmic drugs. A typical case is the administration of Adenosine which, meant to terminate AT, frequently triggers AF by increasing potassium conductance via K,ACh [18]. Analogous adverse effect is commonly found with Dobutamine and with other antiarrhythmic agents, whose mechanism of action involves shortening of ECG Effective Refractory Period (ERP) and/or of atrial Action Potential Duration (APD) by increasing repolarizing potassium currents [18][19][20][21].
An increasingly adopted and promising bradycardic agent is Ivabradine (Iva), which slows heart rate without significantly affecting inotropy [22], and thus it is widely used in the treatment of angina [23,24]. Unlike beta-blockers, Iva acts by directly closing channels, the main responsible for cardiac membrane pacemaker depolarization [25]. Despite that, its administration is associated with a 14% increased risk of AF, when compared with other bradycardic agents [26], with the involved mechanism resembling that of the sick sinus syndrome [27].
Finally, we note that not only class IV antiarrhythmic drugs ( CaL blockers) are occasionally associated with a higher risk of AF [28,29], but also CaL downregulation is a common finding in conditions when AT precipitates to AF [30].
Most of the substrates leading to AF include decrease in heart rate suggesting a common, even if not exhaustive, underlying mechanism. Despite this fact, evidences explaining why bradycardia should favor initiation of AF have not been provided so far. The aim of the present simulation study is to explore, at the most elementary tissue level, that is, interaction within cell pairs and a one-dimensional strand, the role of pacemaker source, atrial sink relationship in establishing more or less safe impulse conduction under chronotropic conditions typically involved in the treatment of AF. Safety factor measurements are based on changes in intercellular coupling and in the SAN/atrial membrane surface ratio. Mathematical modeling allows quantification of such a factor for SAN-atrial drive under ion channels modulation, which can better direct pharmacological and clinical antiarrhythmic strategies.

Integrating Single Cell and Cell Pairs Action Potentials.
All simulations reported in the present study have been performed by means of two mathematical models of rabbit cardiac action potentials (APs): the Severi et al. sinoatrial model [31] and the Lindblad et al. atrial model [32]. Both models were downloaded and recompiled in their MATLAB version by means of COR facility at http://cor.physiol.ox.ac.uk/. The "ode15s" solver built into the R2010b version of MATLAB (The MathWorks, Inc., USA) was used to integrate the models equations. All simulations were run on a PC with Intel Core 2, 3 GHz CPU. SAN-atrial electrical cell coupling was simulated by solving the following system of differential equations: ,at − ,sa = ( ,sa ,sa + ion,sa ) , ,sa − ,at = ( ,at ,at + ion,at ) , where is the global electrical gap junctional resistance (MΩ) between the SAN and atrial elements, ,sa and ,at are the electrical capacitance of the original sinoatrial and atrial models (32 and 100 pF, resp.), ion,sa and ion,at are the total ionic current (nA), and ,sa and ,at are the membrane potential (mV) in both models. The left term of both equations represents the electrotonic current (nA) flowing, through the and with opposite sign, across SAN and atrial membrane, respectively. is a scaling factor that allows the control of the size (membrane surface) of the SAN cell model. = 2, for example, simulates electrical coupling of two parallel connected SAN cells firing synchronously, that is, a single SAN cell with twice the membrane surface and same channel densities, to the same atrial cell. was made varying throughout all simulations between 1 and 20, in order to specifically investigate an ( and ) region where the 3 pacing conditions (NP&ND, P&D, and P&ND) were all represented together with all the possible transitions between them. Sinoatrial AP waveforms, either in control or under any simulated intervention, were integrated starting with initial conditions taken after 120 s of spontaneous beating. Similarly, atrial AP waveforms were obtained by using initial conditions taken after 100 simulated electrically driven APs at a pacing cycle length (CL) intrinsic for the rabbit heart (355 ms).

Integrating AP Propagation in a One-Dimensional Cell
Strand. Action potential propagation along a strand made of a variable number ( SA ) of SAN cells and 10 atrial cells (Severi and Lindblad model, resp.), longitudinally connected with each other with an electrical resistance , was simulated by solving, for each th cell, the following differential equation: Source-sink properties were investigated by varying SA between 1 and 10 and from 5 up to 200 MΩ, step 5 MΩ. In some cases, in the last cell of the atrial side of the cable, a 2 ms suprathreshold current injection was simulated in order to elicit an AP propagated in atrial-SAN direction.  [31] for the description of both treatments on ion channels, membrane transporters, and calcium dynamics.

Modulation of the Membrane and Calcium
Clocks. As in the Severi et al. original paper, we simulated the application of 3 M Iva, which corresponds to the block of 66% of conductance [31,33]. Since a univocal formulation for Ryanodine application in order to silence calcium clock [34,35] is not provided for the Severi et al. model (whose CL barely changes after complete block of sarcoplasmic reticulum (SR) calcium release), we followed the indications by Maltsev and Lakatta [36] and simulated a Ryanodine-like (Rya * ) application by turning off SR calcium up-take and simultaneously downregulating SAN CaL by 34%, the aim being, as in all bradycardic treatments under study, to match the 28% CL prolongation found with Iva.

Effects of Classes III and IV Antiarrhythmic Agents.
The action of a class III antiarrhythmic agent like Dofetilide was simulated in single SAN cells and in cell pairs by a 74% Kr downregulation [37] and that of a typical class IV agent like Verapamil via 61.5% reduction of CaL conductance [38]. In cable simulations only 10% downregulation of both conductances was applied.

Results
Aim of the present study is to compare the effect of different simulated chronotropic interventions on the source-sink relationship between a SAN and an atrial rabbit AP model, either arranged in cell pairs or in a cells strand, in order to study differences in the electrical strength of the pacemaking source and discuss the corresponding significance to the initiation of AF.  correspond to a condition when SAN pacing occurs but does not provide enough electrotonic current to drive the atrial cell at its own frequency ("pace and not drive, " P&ND, red in panel (d)). This is shown in panels (c1) and (c2): in (c1) is set to be infinitely high, cells are practically uncoupled, and SAN cell displays its intrinsic spontaneous rhythm whereas atrial cell is quiescent at its resting potential ( = −75 mV). In (c2) is high but small enough to allow depolarizing electrotonic current to flow across it and induce subthreshold depolarizations into the atrial model cell. By further decreasing , we find a window of its values where, despite the electrical load exerted by the more polarized atrial membrane, SAN cell keeps firing and provides enough electrotonic current to synchronously drive the atrial cell ("pace and drive, " P&D, blue in figure). Finally, when is decreased below a certain value, the electrical load of the atrial cell prevails, preventing SAN firing and, with that, its own pacing ("not pace and not drive, " NP&ND, white in figure). It is found experimentally that these three conditions are not necessarily present for any cell pair but depend on the relative amount of membrane surface of the source and the sink cells [39]. This latter property is reproduced by the scaling parameter of the first equation of System (1) reported in Methods. Thus the source-sink behavior of the SAN-atrial cell pair can be summarized in a graphic like that reported in panel (d), previously described by Joyner and van Capelle [40], and obtained here by varying (step = 1) and (step = 50 MΩ). Given an excitable pacemaking source, electrically coupled to an excitable resting sink, there is, in general, an value separating NP&ND from P&ND or from P&D conditions. Also, only above a given value of , P&D condition can take place ( > 4 in panel (d)). The ensemble of ( and ) values which allow P&D condition fills the blue area, whose surface, as we will show, is a measure of the safety factor with which the SAN membrane can drive the atrial sink; we refer to the numerical value of such a factor as SASF. SASF surface scales dimensionally as MΩ. The conductance value of a single gap junctional channel (50 pS [41]) sets a nonlinear discretization for the -axes of Figure 1(d), which would allow us to convert the measured surface into a given number of gap junctional channels. Nevertheless, for the sake of the relative quantification of the strength of the source in different pacing conditions, our assumption of continuity for the -axes (and therefore for SASF) serves well our scope. The ability of an excitable current source to sustain a safe AP conduction is usually measured as safety factor (SF), and slightly different methods have been proposed in order to estimate it [42]. For each value of we could derive thedependency of SF, which we calculated according to the Leon and Roberge formulation: with ion the ionic current and the total membrane current flowing, respectively, across the atrial cell membrane and integrated over the time when they are negative ( and time windows, resp.). In doing so, we found that, for each value of the scaling factor , SF is greater than 1 only for values falling within the P&D area (see representative example in lower inset of Figure 1(d)). SASF surface gives, in other words, a compact representation of source-sink properties when both and membrane scaling ( ) are subject to changes. Figure 2 we summarize a number of different maneuvers (see Methods) we have performed on the SAN AP model in order to achieve increase (only one case) and decrease of its intrinsic beating rate. The simulated application of 1 M Iso leads to a 22% decrease in pacing CL by dramatically increasing the rate of diastolic depolarization (DDR, V/s) and leaving AP amplitude (APA, mV) unaltered (panels (b) and (h)). In contrast, the simulated application of 11 nM ACh leads to a 28% increase in CL (panel (c)), which is very close to that expected from the linear ACh-dose dependency for CL predicted by Rocchetti et al. in their study on real rabbit sinoatrial cells [43] and recently confirmed by our own study on guinea pig SAN cells [35]. We then simulated other interventions, like downregulating CaL by 61%, downregulating Kr by 74%, downregulating conductance by 66% via simulation of 3 M Iva application [31,33], and simultaneously downregulating SR rate of calcium uptake by 100% and CaL by 34% (panels (d)-(g)) in order to simulate Rya application (Rya * , see Methods). In all these instances we fine-tuned parameters in order to achieve exactly the same bradycardic effect, that is, the same 28% CL increase obtained with ACh. Corresponding relative changes in APA and DDR are shown in the histogram of panel (h). Absolute values of AP parameters can be found in Table 1.    Beating parameters were measured after 100 s of simulated spontaneous firing, with AP waveform in its steady state conditions. Cycle length (CL), maximum diastolic potential (MDP), upstroke value ( peak ), action potential amplitude (APA), action potential duration as measured at −40 mV (APD −40 mV ), diastolic depolarization rate (DDR), maximum value of the first time derivative of AP ( / max ), and threshold value as take-off potential (TOP) are reported.  (1) reported in Methods and varying, in turn, coupling resistance and scaling factor , as explained for Figure 1(d). Each simulation resulted, as shown above, in NP&ND, P&D, and P&ND configurations, all summarized in the color panels of Figure 3. In the case of downregulation of Kr and CaL , we also report in Figure 4, superimposed to the surface profiles shown in Figure 3 (black dotted lines), the color contour plots of the corresponding CL values. Our hypothesis that the P&D area of each graphic (SASF) is a measure of the strength of the SAN source is confirmed by the simulated experiment reported in Figure 5. In each beating SAN cell, taken in steady state conditions (after 120 s, see Methods), we simulated a hyperpolarizing constant current injection of increasing amplitude in order to find the amount of current needed to stop pacemaking within an arbitrarily chosen interval of time (5 s). Longer intervals were also tried and gave qualitatively identical results (not shown). The correlation between the constant current value needed to switch off pacemaking and the corresponding SASF value derived from Figure 3 clearly appears in the histogram of Figure 5(b) and is statistically quantified by the correlation analysis in panel (c) of the same figure.

Chronotropic Interventions on the SAN Firing. In
Whether we consider, for each treatment in the uncoupled condition and for the corresponding coupled configuration, the hyperpolarizing current needed to turn off pacemaking and the SASF value, respectively, we find that both parameters significantly correlate with DDR and with no other AP parameter (see panel (d)).

Transitory Changes in Electrical
Coupling. We hypothesize here that the source-sink SAN-atrial system is set, in physiological conditions, at a given value of ( , ), say (11,590 MΩ); that is, 11 SAN cells are spatially arranged around a single atrial cell and are cumulatively coupled to it with a total of 34 gap junctional channels (assumed to have, as pointed out above, a single channel conductance of 50 pS). We further assume that none of the antiarrhythmic interventions that we simulate changes this set point (none of them has known effects on geometrical nor on gap junctional coupling). On the other hand, we know from the literature [10,11] that AF is often accompanied by gap junctional remodeling. We show in Figure 6 what a transitory (500 ms) closure of only 2 gap junctional channels (yellow path in panels (a) and (c) and in the magnified details of panels (b) and (d)) is expected to cause on the safety factor of SAN pacing. Whereas in control conditions (NT) such slight junctional change leaves SANatrial intercellular conduction within the "safe" area (blue) of P&D behavior, the same modification transiently (horizontal double arrow in panel (d)) brings the CaL downregulated system into the P&ND region (red), where a single beat fails to be conducted from the SAN to the atrial membrane (bottom of panel (d)).

Relative Changes in Ion Currents for Each Treatment.
In order to identify the ion currents primarily involved in determining the strength of the pacemaking source when coupled to an atrial sink, we performed a series of simulations like that reported in Figure 7 Figure 8 shows how SAN-atrial conditions characterized by a larger "safety factor" correspond to a larger increase in when passing from uncoupled to coupled conditions (panel (a)). Similarly, they correspond to a larger decrease of Ks when coupled (panel (b)). Correlation coefficients were 0.77 and 0.93, respectively. Positive significant correlations were also found for CaL and NaK , though involving much smaller changes (data not shown).

Conduction Delay and Entrainment of Repolarization.
We have shown, up to this point, that the area of the blue surfaces in Figure 3 measures the safety factor SASF associated with spread of pacemaker activity from SAN to atrial cells. We have also observed that SASF is different among bradycardic interventions which lead to the same decrease in SAN beating rate. We also wanted to test how versus plots, like that reported in Figure 1, were derived for each chronotropic intervention described in Figure 2. Iso corresponded to a 22% CL decrease, whereas all bradycardic maneuvers increased CL by 28%. Color code as in Figure 1.  safety factor of pacemaking is going to affect AP conduction delay and entrainment of repolarization in our cell pair model. We did that for identical coupling conditions, that is, for SAN and atrial cells coupled with = 14 and = 550 MΩ, which correspond to a P&D state for all treatments.
Coupled conditions were simulated for 5 s and the last conducted APs reported in Figure 9(a). Panel (b) shows how AP conduction delay decreases as safety factor increases, being their relationship well fitted ( = 0.89) by a hyperbolic function. Same for relative APD −40 mV difference between SAN and atrial APs ( = 0.87) (panel (c)), where APD −40 mV measures the time lapse between the peak of the initial time derivative and the time when reaches the value of −40 mV. Very similar results (data not shown) were found when and were chosen in order to roughly match the center of the P&D area for each treatment. In other words, a more robust pacemaker source not only guarantees a much higher safety factor for the spread of excitation but also leads to a shorter AP conduction delay and a better entrainment of repolarization as well.

Simulations of a Linear SAN-Atrial Cell Strand.
In order to test whether SASF evaluation as a measure of source strength could be applied also in a cable-like AP propagation, we performed simulations on a SAN-atrial chain (see scheme, top panel of Figure 10) of cells, electrically connected with a variable junctional resistance (see Methods). For small values, the more polarized atrial side of the cable completely depressed spontaneous firing in the left SAN side (NP&ND condition in Figure 10(a1) and white area in panel (a)). As increased, and in analogy with what is observed for cell pairs reported in Figure 3, spontaneous APs were generated in the SAN side of the cable and conducted to the atrial side (P&D condition in Figure 10(a2) and blue area in panel (a)). A further increase prevented APs generated in the SAN side from being conducted to the atrial side (P&ND condition in Figure 10(a3) and red area in panel (a)). The same is shown in panels (b) and (c) of Figure 10 in the case of 10% downregulation of Kr and CaL , respectively. These results confirm, in a cable-like model, the effects reported above on cell pairs, where Kr downregulation increases and CaL downregulation decreases the strength of the source, that is, the safety factor for AP conduction (see also Figure 10(d)).
Given the relevance of unidirectional block (UB) in the initiation of AF, we tested for it the cable-like model described in Figure 10, with the additional aim of relating UB to source-sink properties as measured with SASF surfaces. UB developed when AP propagation was simulated in control conditions for SA = 10 and = 110 MΩ (yellow circle in the SASF representation); that is, conduction failed when the SAN side of the cable was let free of beating at its own intrinsic frequency (CL = 341 ms) but succeeded when the end of the atrial side of the cable was electrically paced at the same frequency (Figure 11(a)). This was indeed expected, since the yellow circle is located into the P&ND (red) region of the graph. As we have shown in Figure 10 and reported here in the inset of panel (b), Kr downregulation extends the P&D (blue) area, which now includes the yellow circle. Accordingly, conduction becomes permissive also in the SAN-atrial direction, making Kr downregulation effective, in this case, also in removing UB.

Discussion
The present modeling study aims to investigate, at the cellular level, the source-sink SAN-atrial relationship in conditions which are critical for the development of AF. We first define a graphical representation which quantifies the safety factor for impulse conduction from SAN to atrial cell membrane into a single parameter (SASF). We then use SASF in order to compare identical bradycardic effects on SAN firing, which follow pacemaker autonomic regulation or antiarrhythmiclike treatments targeting AF.
Representations of the source-sink properties of SANatrial electrical coupling in graphics like those of Figures 1 and  3 have been previously described by Joyner and van Capelle [40], who showed that some degree of electrical uncoupling is   an essential design for proper SAN-atrial conduction. Several other simulation studies have been performed on this same issue [44][45][46], though not addressing its potential relevance for the initiation of AF. Basically it has been demonstrated that some degree of SAN to atrial interdigitation improves reliability of conduction and that the critical junctional resistance allowing AP entrainment in SAN-atrial cell pairs can easily fall within the GΩs range. Though the geometry of SAN is far to be fully elucidated and the possible arrangement of mutual intercellular coupling still needs to be clarified, nevertheless it is assumed that SAN ability to electrically drive the larger atrial volume is based on the presence of different cell types within its interior [47,48], on a complex gradient of gap junctional distribution [49,50], and on some favorable intercellular geometry [47]. This latter is achieved by interdigitation of different cell types [46,50], which will likely result into a many-to-one SAN-atrial cells connection. Graphics (Figures 1(d) and 3) reporting versus (number of SAN cells coupled to a single atrial cell) compactly summarize these properties, by defining, for each , therange that allows SAN membrane to electrically drive the atrial cell (P&D blue area). The same type of representations can be applied also to cable-like AP propagation (Figures 10  and 11).
Autonomic modulation, which is known to shift the leading pacemaker site within the SAN [51], is involved into the initiation and maintenance of AF [52]. Particularly, cholinergic stimulation is indicated as one of the main factors for the initiation of spontaneous AF, even though the exact knowledge of underlying mechanism is unknown [52]. Our simulations of cholinergic (equivalent to 11 nM ACh) and adrenergic (equivalent to 1 M Iso) modulation of the Severi  5(b)), which is likely underlying the increased propensity to AF under cholinergic hyperstimulation [53] and related to the dramatic differences in coupling induced-Δ in the two instances (Figure 7(b)). Indeed, the correlation between SASF and Δ (Figure 8(a)) suggests that coupling-induced changes in play a major role in determining SASF. Recent studies have attributed the molecular mechanism of SAN rhythm to the interplay between two clocks, one involving the hyperpolarization activated cation current (membrane clock) [54] and the second attributable to activation of the electrogenic NaCa exchanger by spontaneous SR releases of calcium (calcium clock) [36]. The two clocks have been shown to possess different intrinsic dynamic properties, which can contribute to the physiopathology of the heart rhythm [35]. Thus, it is of interest to investigate the source-sink changes induced into the SAN when the two clocks are separately downregulated in order to obtain the same bradycardic effect. Whereas both maneuvers lead to a decrease in SASF, we note that the simulated application of Iva reduces SASF to a lesser extent when compared to Rya * (Figures 3 and 5(b)). If, from one hand, the downregulation of the membrane clock appears to be safer than that of the calcium clock, on the other hand the Iva-induced reduction in SASF might be the mechanism underlying the slight increase in AF development found when Iva is used, as it is frequently, instead of beta-blockers in the treatment of ischaemic diseases [26]. Sinus node dysfunction is often associated with initiation and maintenance of AF [13], and an electrotonically weaker SAN, even in the absence of structural, electrical, or junctional remodeling, is more prone to transient SAN blocks, leading to a substrate favoring arrhythmias. This is shown, for instance, in the case of CaL downregulation (Figure 6), where even a very rapid and transient 5% increase in , ineffective in physiological conditions (panel (b)), brings the system out of the P&D area, leading to a 1-beat failure in SAN-atrial conduction (panel (d)). A comprehensive study of the SASF associated with pharmacological treatments of AF is beyond the scope of the present study. We limit our analysis to two ion channels modulatory effects involved in the action of classes III and IV antiarrhythmic drugs, like Dofetilide and Verapamil, respectively, known to decrease Kr and CaL conductance [38,55]. Furthermore, since our aim was to compare bradycardic effects of the same extent (28% reduction of CL), we did not try to match, in our simulations, the therapeutic doses clinically administered for each one of these agents.
CaL downregulation, which results from class IV antiarrhythmics administration, is also a common finding in conditions when atrial tachycardia develops into AF [30]. It is interesting, at this regard to observe that a decrease in CaL brings about a dramatic reduction of SASF both in cell pairs (Figures 3 and 5(b)) and in cable propagation (Figures 10(c) and 10(d)), which might explain the transition to AF due to (1) a decreased safety factor for AP conduction, (2) an increased conduction delay, and (3) a diminished control in the entrainment of repolarization (Figures 9(b) and 9(c)). This is in agreement with the experimental findings of Hondeghem and Hoffmann in isolated rabbit hearts [29], and with documented Verapamil-induced SAN asystole [28] as well. An opposite effect on SASF is brought about in SAN-atrial cell pairs by the downregulation (−74%) of Kr in order to simulate the action of pure anti-Kr class III antiarrhythmic agents, like Dofetilide and Nifekalant [56,57]. In this case, the decrease of Kr leads to (1) an increase in SASF, likely mediated by the coupling-induced greater increase in and in CaL (Figures 7(b) and 8(a)), (2) a decrease in conduction delay, and (3) a better entrainment of repolarization, all synergic with a reduced proneness to develop AF. Furthermore, the larger coupling-induced reduction in Ks (Figures 7(b) and 8(b)) is going to reduce the likelihood of AF initiation by further increasing APD [58].
The role of CaL in granting safety of conduction can be summarized as follows. The depolarizing (negative) ion exceeding the depolarizing (negative) electrotonic (see (3)), which makes AP spread safe (SF > 1) or, in other words, makes the source sink system fall into the P&D area, is almost entirely due to the contribution of CaL , which is far the larger depolarizing current in SAN membrane. When this current is reduced, SASF gets reduced as well ( Figure 3) and, with it, do the strength of the SAN source and the beating frequency (bottom panel of Figure 4). When, on the opposite,  Kr downregulation is applied in the same cable conditions (yellow dot in the inset. Note that now it falls into the blue P&D area), AP conduction succeeds in both ways. CaL increases, like during Kr reduction (see Figure 7(b)), the safety of conduction increases as SASF does ( Figure 3); beating frequency decreases (middle panel of Figure 4) more than it does under CaL downregulation.
Finally, in the cable simulation of Figure 11, SASF calculation predicts UB for a given set of passive parameters ( SA and , yellow dots in the insets) in control conditions, and its removal after Kr downregulation following SASF increases. Unidirectional block, together with decrease in conduction velocity, is a well-recognized cause of reentry, which, in turn, as mentioned above, is one of the two main mechanisms underlying initiation of AF [12]. The fact that, as we show, decreasing Kr prevents UB induced by electrical uncoupling, might point to an unexplored property of class III antiarrhythmic agents, which are usually mainly described for their AP prolonging action.
To summarize, it is accepted that the successful SAN drive of the atria is based on a combination of active and passive electrical properties, and on a given geometrical arrangement between the two cell types within the SAN node and at its border [47]. The three factors determine a set point which, in physiological conditions, belongs to the P&D area of the graphic representations of Figure 3. The surface of such area can be used to quantify (SASF) the strength of the SAN as an electrotonic source. We find that (1) paraand orthosympathetic modulations exert opposite effects on SASF, (2) class III antiarrhythmic agents increase it and class IV decreases it, and (3) the same bradycardic effect obtained by separately downregulating the membrane and calcium clock corresponds to a smaller and larger SASF reduction, respectively.

Conclusions
Among the several conditions leading to AF, the role of safety in SAN-atrial AP conduction in terms of sourcesink relationship has not been previously explored. Our simulation study shows that such relationship should be taken into account, together with usually recognized active, passive, and geometrical factors, in order to reconstruct the mechanism underlying the initiation and maintenance of AF. Particularly, a parameter like SASF compactly describes the strength of SAN source and thus might draw light on further SAN modeling and on the mechanism of action of antiarrhythmic agents.

Limitations of the Study
Despite the critical role that the morphological structure of SAN-atrial arrangement, including the complex interaction between several recognized cell types [47,59,60], is going to play into the atrial transition to fibrillation, such information could not possibly be considered in our simulations, given the still incomplete and sparse literature on the subject, and the much simpler geometry of our simulation setting. On the other hand, our simplified approach is meant to provide a more general background on SAN source-atrial sink properties, which can be improved as the underlying structure will be elucidated. A limitation of the Severi et al.
SAN AP model is the lack of sensitivity of beating rate to reduction of SR calcium release. Thus, as mentioned above, in order to simulate Rya application, we followed a procedure already in use by other authors [36]. Finally, unfortunately the measure of SASF cannot possibly be pursued in any thinkable experimental preparation, which confines its use to the theoretical in silico explanation of experimentally observed phenomena.