Synergy between Membrane Currents Prevents Severe Bradycardia in Mouse Sinoatrial Node Tissue

Bradycardia is initiated by the sinoatrial node (SAN), which is regulated by a coupled-clock system. Due to the clock coupling, reduction in the ‘funny’ current (If), which affects SAN automaticity, can be compensated, thus preventing severe bradycardia. We hypothesize that this fail-safe system is an inherent feature of SAN pacemaker cells and is driven by synergy between If and other ion channels. This work aimed to characterize the connection between membrane currents and their underlying mechanisms in SAN cells. SAN tissues were isolated from C57BL mice and Ca2+ signaling was measured in pacemaker cells within them. A computational model of SAN cells was used to understand the interactions between cell components. Beat interval (BI) was prolonged by 54 ± 18% (N = 16) and 30 ± 9% (N = 21) in response to If blockade, by ivabradine, or sodium current (INa) blockade, by tetrodotoxin, respectively. Combined drug application had a synergistic effect, manifested by a BI prolonged by 143 ± 25% (N = 18). A prolongation in the local Ca2+ release period, which reports on the level of crosstalk within the coupled-clock system, was measured and correlated with the prolongation in BI. The computational model predicted that INa increases in response to If blockade and that this connection is mediated by changes in T and L-type Ca2+ channels.


Introduction
The sinoatrial node (SAN) maintains the heart rate at rest in a range that allows for an instantaneous increase when a person performs work or generates a fight or flight response. While low heart rate at rest, also termed sinus bradycardia, is not a risk factor per se, when induced by specific drugs, it is associated with incident cardiovascular diseases and elevated mortality [1]. Severe symptomatic bradycardia can lead to atrial fibrillation and a decrease in oxygen supply [2]. The body may engage a "fail-safe" mechanism to prevent such conditions. SAN function is maintained by a coupled-clock system that consists of membrane ion channels, exchangers, and pumps (M clock), and an internal Ca 2+ clock, namely, the sarcoplasmic reticulum (SR). Both clocks communicate via Ca 2+ and Ca 2+ -activated adenylyl cyclase (AC)-cAMP-PKA signaling, with local Ca 2+ releases (LCRs) from the SR indicating the degree of coupling [3]. Due to the connectivity of the clocks, a change in a single membrane component can affect the others, and impair or upregulate clock function [4]. Thus, if one component in the coupled-clock function fails, another component can compensate for the reduced pacemaker function.
Here, we focus on I pacemaker ("funny") current (I f ), which is one of the elements that contribute to the generation of spontaneous diastolic depolarization and an important stabilizer of heart rhythm [5]. Reduction in I f was previously shown in heart failure conditions [6,7], aging [8], pulmonary hypertension [9], diabetes [10], and atrial fibrillation [11]. However, in such clinical scenarios, the SAN is still beating, and most patients show no signs of severe bradycardia. This may be ascribable to I f , which, when in synergy with other clock components, may upregulate their function and prevent severe bradycardia.
We hypothesize that SAN pacemaker cells bear an inherent fail-safe mechanism driven by synergy between ion channels and I f , and that this synergy is mediated by a coupledclock mechanism and specifically by Ca 2+ -dependent channels. To prove these hypotheses, we isolated SAN tissues from C57BL mice ( Figure 1A) and used confocal microscopy to measure Ca 2+ signaling in pacemaker cells within the tissue ( Figure 1B,C). We applied ivabradine (IVA), which reduces I f , alone or in combination with tetrodotoxin citrate (TTX), or tetramethrin (TMR), blockers of sodium current (I Na ) and T-type Ca 2+ current (I CaT ), respectively, which are two potential synergy components. A computational model of SAN cells ( Figure 1D) was used to understand the internal interactions between cell components. coupled-clock mechanism and specifically by Ca 2+ -dependent channels. To prove these hypotheses, we isolated SAN tissues from C57BL mice ( Figure 1A) and used confocal microscopy to measure Ca 2+ signaling in pacemaker cells within the tissue ( Figure 1B,C). We applied ivabradine (IVA), which reduces If, alone or in combination with tetrodotoxin citrate (TTX), or tetramethrin (TMR), blockers of sodium current (INa) and T-type Ca 2+ current (ICaT), respectively, which are two potential synergy components. A computational model of SAN cells ( Figure 1D) was used to understand the internal interactions between cell components.
Here, we experimentally show synergy between If and INa. Our model predicted that synergy exists between If and INa and that this synergy is mediated by the ICaT and L-type Ca 2+ current (ICaL). We also found that synergy between If and INa exists at the level of the LCR period ( Figure 1C), and promotes the positive feedback between ICaL, SR Ca 2+ , and Na + -Ca 2+ exchanger current (INCX) [12]. Taken together, our data show that a fail-safe mechanism is driven by a connection between ion channels and is mediated by membrane Ca 2+related mechanisms.   Here, we experimentally show synergy between I f and I Na . Our model predicted that synergy exists between I f and I Na and that this synergy is mediated by the I CaT and L-type Ca 2+ current (I CaL ). We also found that synergy between I f and I Na exists at the level of the LCR period ( Figure 1C), and promotes the positive feedback between I CaL , SR Ca 2+ , and Na + -Ca 2+ exchanger current (I NCX ) [12]. Taken together, our data show that a fail-safe mechanism is driven by a connection between ion channels and is mediated by membrane Ca 2+ -related mechanisms.

Results
2.1. I CaT and I Na Are Upregulated When I f Is Inhibited: Computational Evidence I f is one of the regulators of the spontaneous beating of the SAN. Inhibition of I f by pharmacological drugs or by genetic inhibition of HCN 4 does not stop the spontaneous beating. To explore whether a "fail-safe" mechanism is engaged to prevent severe bradycardia when I f is reduced, we used our previously published computational model of the single mouse SAN cell [13] ( Figure 1D). Inhibition of I f was simulated by reducing the I f maximal conductance coefficient (g If , see supplements) to 10% of its basal value ( Figure 1E). Figure 2 shows the effect of I f inhibition on main coupled-clock mechanisms. The model predicted that a reduction in g If prolonged the beat interval (BI), calculated as the time interval between two peaks of action potential (AP), and also associated it with increased amplitudes of I CaT (+114%), I Na (+48%), and Ca 2+ flux through the ryanodine receptors (RyR (+32%)), I NCX (+18%), I CaL (+10%), and other transmembrane currents and of Ca 2+ cycling parameters (<10%). The increase in I CaT and I Na in response to I f inhibition implies a connection between these currents, engaged in a potential mechanism to restrain the prolongation of BI.
previous Ca 2+ transient peak to the LCR onset. Moreover, 50% relaxation (blue line) is the time from the previous minimal Ca 2+ transient to the 50% relaxation point of the current Ca 2+ transient. (D) Scheme of SAN cell model parameters including membrane currents, ion concentrations, Ca 2+ cycling, and drug perturbations. (E) Simulated "funny" current (If) in the basal state (black) and with If inhibition, to simulate the ivabradine (IVA) treatment effect (yellow). Simulated sodium current (INa) in the basal state (black) and with INa inhibition, to simulate the tetrodotoxin (TTX) treatment effect (red).

ICaT and INa Are Upregulated When If Is Inhibited: Computational Evidence
If is one of the regulators of the spontaneous beating of the SAN. Inhibition of If by pharmacological drugs or by genetic inhibition of HCN4 does not stop the spontaneous beating. To explore whether a "fail-safe" mechanism is engaged to prevent severe bradycardia when If is reduced, we used our previously published computational model of the single mouse SAN cell [13] ( Figure 1D). Inhibition of If was simulated by reducing the If maximal conductance coefficient (gIf, see supplements) to 10% of its basal value ( Figure  1E). Figure 2 shows the effect of If inhibition on main coupled-clock mechanisms. The model predicted that a reduction in gIf prolonged the beat interval (BI), calculated as the time interval between two peaks of action potential (AP), and also associated it with increased amplitudes of ICaT (+114%), INa (+48%), and Ca 2+ flux through the ryanodine receptors (RyR (+32%)), INCX (+18%), ICaL (+10%), and other transmembrane currents and of Ca 2+ cycling parameters (<10%). The increase in ICaT and INa in response to If inhibition implies a connection between these currents, engaged in a potential mechanism to restrain the prolongation of BI.

There Is No Synergy between If and ICaT
The model suggested that there is synergy between If and ICaT in single pacemaker cells. As shown before [14], pacemaker cells within the SAN demonstrate significant

There Is No Synergy between I f and I CaT
The model suggested that there is synergy between I f and I CaT in single pacemaker cells. As shown before [14], pacemaker cells within the SAN demonstrate significant heterogeneity in the densities of I f , I CaL , potassium currents (I K ), and other currents, which might be dictated by their localization in different pacemaker cells clusters. Thus, we chose to use pacemaker tissues that contain different clusters to assess the possibility of an association between the channels. IVA (3 µM) and TMR (10 µM) were applied to intact mouse SAN tissue to block I f and I CaT , respectively, and their individual and combined effects on the spontaneous BI were measured. Previous works demonstrated that 3 µM IVA blocks the HCN 4 channel [15] without affecting the T-type, L-type, delayed outward potassium current densities [16], or SR Ca 2+ content, while 10 µM IVA was shown to affect the T or L-type channels [16]. TMR at 0.1 µM was shown to block I CaT in single rabbit SAN cells, while 50 µM TMR abolished both I CaT and I CaL [17]. Yet, as we detected no change in the BI or LCR period of mouse SAN cells residing in the SAN tissue after administration of 0.1 µM TMR (N = 6, Figure S1), the current experiments were performed with 10 µM TMR. Figure 4A shows representative time courses of Ca 2+ signaling before (control) and after the administration of IVA and following the administration of IVA+TMR. Figure 4B shows representative time courses of Ca 2+ signaling before (control) and after the administration of TMR and following the administration of TMR+IVA. Note that the same cell was traced before and after each treatment (paired measurements). IVA increased the BI by 49 ± 20%, compared to its control, while TMR increased the BI by 28 ± 5% compared to its control. Administration of both IVA and TMR increased the BI by 64 ± 13% ( Figure 4C). Beats per minute (BPM, calculated as 60,000/BI per cell) decreased by 24 ± 7% with IVA, 20 ± 3% with TMR, and 34 ± 4% with both IVA and TMR, compared to the control. Namely, there was no additive effect upon application of both blockers together.
BI variability (BIV), estimated by the average standard deviation of the BIs in each cell, was then calculated to determine whether each blocker and blocker combination affect BI periodicity as well. Compared to the control, BIV was increased by 105 ± 47% with IVA, by 69 ± 29% with TMR, and by 462 ± 134% with both IVA and TMR ( Figure 4D). Namely, there was an additive effect of the two blockers on BIV.
To determine the effect of each blocker and blocker combination on Ca 2+ transient and LCR properties, global and local Ca 2+ signaling parameters were measured before and after each drug treatment. Ca 2+ transient amplitude (estimated by the fluorescence ratio [F/F0]) and 50% Ca 2+ transient relaxation time were not affected ( Figure 4E,F) by IVA, TMR, or by IVA+TMR. The LCR period, defined as the time from the previous Ca 2+ transient peak to the LCR onset (as illustrated in Figure 1C), was prolonged only upon administration of IVA (88 ± 42%) and IVA+TMR (79 ± 21%), but not upon treatment with TMR alone ( Figure 4G).

I f and I Na Blockers Have a Synergistic Effect on BI
To determine whether synergy exists between I f and I Na , 3 µM IVA and/or 5 µM TTX were applied to pacemaker cells within intact mouse SAN tissue and the spontaneous BI was measured. Figure 3A shows representative time courses of Ca 2+ signaling before (control) and after the administration of IVA and following the additional administration of TTX. Figure 3B shows representative time courses of Ca 2+ signaling before (control) and after the administration of TTX and following the additional administration of IVA. Compared to untreated cells, IVA prolonged the BI by 54 ± 18% ( Figure 3C), while TTX prolonged the BI by 30 ± 9%. In contrast to TMR, a synergistic effect was observed upon administration of both IVA and TTX, which prolonged the BI by 143 ± 25%. BPM decreased by 25 ± 6% with IVA, 18 ± 4% with TTX, and 50 ± 6% with both IVA and TTX, compared to the control. No difference was found between the change in BI in cells treated with TTX following IVA (IVA+TTX) versus cells first treated with IVA and then with TTX (TTX+IVA). Taken together, the combined treatment enhanced the mono-drug blockade, suggesting that I Na can play a role in the fail-safe mechanism when I f is reduced.    BIV increased by 106 ± 40% with IVA, by 176 ± 83% with TTX, and by 100 ± 36% with the IVA and TTX combination ( Figure 3D). Namely, the combined drug treatment did not further enhance the BIV increase obtained with each blocker separately.
The Ca 2+ transient amplitudes did not change after the administration of IVA or IVA+TTX. In contrast, TTX alone decreased the Ca 2+ transient amplitude by 13 ± 3% compared to the control ( Figure 3E). In addition, TTX treatment led to a 5 ± 2% increase in the 50% Ca 2+ transient relaxation time (Figure 3F), while IVA and IVA+TTX had no effect.
The LCR period was prolonged by 81 ± 31% on application of IVA, by 49 ± 18% with TTX, and by 164 ± 46% with both IVA and TTX ( Figure 3G). Thus, the nonlinear change in BI in response to each blocker compared to their combination was also reflected in a nonlinear change in the LCR period.

The Molecular Mechanisms That Mediate the Synergy between I f and I Na : Computational Evidence
The nonlinear effect of I f and I Na blockers on BI and Ca 2+ dynamics is likely due to feedback in the pacemaker cell mediated by still unknown mechanisms. Because specific blockers do not exist for each ion channel type, and because the clock mechanisms are coupled, it is experimentally impossible to test the underlying mechanisms. Therefore, we used our model to predict the internal mechanisms mediating the crosstalk between I f and I Na that was found experimentally.
The TTX effect was simulated by reducing the channel maximal conductance coefficients (see supplements and Figure 1E); g Na1.1 , the TTX-resistant maximal channel conductance coefficient, was reduced to 80% of its basal value and g Na1. 5 , the TTX-sensitive maximal channel conductance coefficient, was reduced to 50% of its basal value. Lower values led to model instability. Figure 5A shows a simulation of the main coupled-clock mechanisms under basal conditions, and on application of IVA and/or TTX. The model predicted a BI increase of 32% with I f inhibition, 31% with I Na inhibition, and 70% with I f and I Na inhibition together ( Figure 5A, internal figure), supporting the nonlinear commutative effect of IVA and TTX shown experimentally. In parallel, I Na inhibition significantly decreased the amplitude of I CaT (−74.5%), Ca 2+ flux through the RyR (−48%), and I NCX (−37%), together with a mild decrease in the amplitudes of I f (−17%), I CaL (−17%) and intracellular Ca 2+ (−16%) ( Figure 5A). Taken together, a reduction in I f leads to a major increase in I Na , while a reduction in I Na does not increase I f .
To determine whether an increase in I Na is key to prevention of severe bradycardia, we simulated a 'current clamp' by fixing I Na to its basal state, while I f was blocked in parallel. In this way, the increase of I Na is blocked, and its effect on the BI can be evaluated. The model predicted that I f inhibition with I Na clamped causes a longer BI than I f inhibition with varying I Na , which indicates that the increase in I Na in response to I f inhibition restrains the increase in BI at the cellular level ( Figure 5B).
We then explored whether the synergy between I f and I Na also exists if I f is increased. Figure S2 shows that I Na decreased in response to a ten-fold increase in g If . Thus, the negative feedback between I f and I Na serves as a fail-safe mechanism for both bradycardia and tachycardia.
To uncover the internal mechanisms that mediate between the reduction in I f and the increase in I Na , we tested the individual effects of various cellular parameters by 'clamping' each parameter to its basal state, while blocking I f in parallel. Only fixation of both I CaT and I CaL restrained the indirect increase in I Na amplitude caused by I f inhibition, by bringing the I Na amplitude closer to its basal value and decreasing the BI (Figure 6). Fixation of I CaT ( Figure S3) or I CaL ( Figure S4) led to a smaller decrease in I Na amplitude, while fixation of I NCX ( Figure S5) or RyR flux ( Figure S6) did not reduce the indirect increase in I Na caused by I f inhibition. Taken together, I CaT and I CaL primarily mediate the increase in I Na in response to a decrease in I f .

Discussion
The present study investigated the synergy between ion channels in pacemaker cells within mouse SAN tissue. The experimental measurements showed that the change in the spontaneous BI when both If and INa were blocked was higher than the sum of changes induced by each blocker individually (synergistic effect). Together with the model predic tion that INa increases in response to If blockade, these results support the first hypothesis that a fail-safe mechanism is an inherent feature of SAN pacemaker cells and is driven by feedback among ion channels. The experimental effect on the LCR period when both I and INa were blocked compared to their individual effects and the model prediction tha the feedback between If and INa is eliminated when ICaT and ICaL are clamped support the second hypothesis that this effect is mediated by coupled-clock mechanisms, specifically ion channels that are affected by intercellular Ca 2+ dynamics.
Our combined experimental measurements and numerical model simulations sug gested that synergy exists between If and INa. The experiments showed that when both I and INa were blocked, the effect on the spontaneous BI was higher than the summed effec of the two blockers. The model showed a similar phenomenon when If and INa were inhib ited by reducing their maximal conductance coefficients. It predicted that when If is re duced, INa increases. An increase in INa eliminates further deceleration of diastolic depo

Discussion
The present study investigated the synergy between ion channels in pacemaker cells within mouse SAN tissue. The experimental measurements showed that the change in the spontaneous BI when both I f and I Na were blocked was higher than the sum of changes induced by each blocker individually (synergistic effect). Together with the model prediction that I Na increases in response to I f blockade, these results support the first hypothesis that a fail-safe mechanism is an inherent feature of SAN pacemaker cells and is driven by feedback among ion channels. The experimental effect on the LCR period when both I f and I Na were blocked compared to their individual effects and the model prediction that the feedback between I f and I Na is eliminated when I CaT and I CaL are clamped support the second hypothesis that this effect is mediated by coupled-clock mechanisms, specifically ion channels that are affected by intercellular Ca 2+ dynamics.
Our combined experimental measurements and numerical model simulations suggested that synergy exists between I f and I Na . The experiments showed that when both I f and I Na were blocked, the effect on the spontaneous BI was higher than the summed effect of the two blockers. The model showed a similar phenomenon when I f and I Na were inhibited by reducing their maximal conductance coefficients. It predicted that when I f is reduced, I Na increases. An increase in I Na eliminates further deceleration of diastolic depolarization and prevents further bradycardia. However, when I Na was reduced, the model predicted only a small change in I f . Thus, the negative feedback between I f and I Na is unidirectional. Furthermore, the model predicted that an increase in I f will reduce I Na . Thus, the predicted synergy between I f and I Na also protects against tachycardia. Note that at the tested IVA concentration, there is a negligible direct effect on I Na [18]. Although I f plays an important role in pacemaking in the SAN, spontaneous beating is still maintained upon its inhibition with even higher (and non-specific) IVA concentrations [15,16] or elimination by genetic manipulation of HCN 4 [11]. The synergy between I f and I Na may also act here as a "fail-safe" mechanism.
Our second main finding was that Ca 2+ -activated channels mediate the synergy between I f and I Na . Our experiments showed that the effect of simultaneous I f and I Na blockage on the LCR period was higher than the summed effect of each blocker. The LCR period is linked to changes in the BI in response to numerous perturbations that affect both the M and Ca 2+ clocks (i.e., a decrease in cAMP/PKA levels or β adrenergic receptor stimulation) [19] and is, thus, considered a readout of the degree of clock coupling [4]. An increase in the LCR period prolongs the diastolic depolarization through delayed activation of I CaL and I NCX (positive feedback between I CaL and I NCX [12]) and through direct activation of BK channels [20], which consequently lead to a longer BI. Thus, the LCR period itself reports on the synergy. The model predicted that when I CaT and I CaL are clamped, the I Na increase in response to a decrease in I f is moderated and restrains further bradycardia. Clamping I CaT and I CaL to their basal values blocks the positive feedback between I CaL -SR-I NCX . Our suggested theory of the synergy between I f and I Na is based on the regulation of pacemaker activity by I CaL and I CaT and their regulation of I CaL -SR-I NCX feedback. The increase in I CaL and I CaT in response to decreased I f maintains intracellular Ca 2+ levels. This, in turn, maintains Ca 2+ release from the SR (RyR flux), preserves the LCR period, and via maintenance of I NCX , accelerates the diastolic depolarization, which prevents further bradycardia. Note that clamping either I NCX or RyR flux did not prevent the feedback between I f and I Na . Our experiments showed that removal of cytosolic Ca 2+ , measured as the Ca 2+ transient relaxation time (T 90 ), did not significantly change in response to IVA or IVA+TTX, suggesting that the SERCA pump function does not directly control the feedback between I f and I Na either. Clamping I Na ( Figure 5) reduced I CaT , I CaL , I NCX , and RyR flux. Reduction in I Na compared to its higher unclamped value prolongs BI, and thus, affects the diastolic depolarization voltage and activation of I CaT and I CaL. These reductions eliminate the positive feedback between I CaL -SR-I NCX , and thus, lead to reduced I NCX and RyR flux. Considering both the experimental results and the numerical model simulation, it can be concluded that the coupled-clock system per se is the fail-safe mechanism that prevents bradycardia in response to a decrease in I f . Note that even under normal conditions, I f is not the only regulator of SAN automaticity. There are several other currents, including I Ca,T , I Ca,L , and I NCX , that contribute to spontaneous diastolic depolarization [21]. Thus, other coupled-clock mechanisms can act as fail-safe mechanisms even when in some, cell I f is close to zero [14].
The model predicted that in parallel to the increase in I Na , there was an increase in I CaT , I NCX , and RyR flux in response to I f inhibition. I f inhibition leads to prolonged BI through increases in I CaT and I CaL , which lead to increased outflux currents and subsequently to increased I NCX . Prolongation of BI allows the SR to slowly refill, resulting in an increase in Ca 2+ available to activate I CaT and I NCX . Note, however, that it also allows for less Ca 2+ release per time interval. Ion channels are also coupled through changes in voltage and other internal signals, and positive feedback has also been shown to exist between I NCX and I CaL [12]. Note, however, that this feedback was bidirectional and was measured if one of the currents increased or decreased. A positive feedback mechanism in cardiomyocytes was also described between the SK channel and I CaL [22] and together with RyR [23].
The model predicted that I CaT increases in response to a reduction in I f . However, experimentally, when I f was blocked together with I CaT , the effects on the spontaneous BI and LCR period were similar to the summed effect of each blocker. Thus, no synergy exists between I f and I CaT . A decrease in I f prolongs the BI and deaccelerates the diastolic depolarization. A longer diastolic depolarization phase allows for the activation of more T-type Ca 2+ channels, increasing their amplitude. It has been suggested that I CaT stabilizes the rate of depolarization when the maximal diastolic potential is more positive [17]. Such conditions were achieved here when IVA was applied. Note that TMR by itself prolonged BI to a similar degree as TTX. However, compared to TTX, it does not have synergistic effect with IVA. Thus, the prolongation of BI by one drug does not necessarily lead to a non-linear effect on BI when combined with another drug.
The BI of SAN is not constant and has some variability, which, itself, is considered an important index that correlates with various heart diseases [24]. Although the absolute change in the average spontaneous BI when both I f and I Na were blocked was different from the sum of changes induced by each blocker individually, the variability was similar between IVA, TTX, and IVA+TTX treatments. Because there is no linear relationship between BI and BIV [24][25][26], it is possible that cell perturbations differentially affect BI and its variability. In contrast, blockade of both I f and I CaT resulted in a change in BI variability significantly greater than the sum of the changes caused by each of the respective blockers individually. It was previously shown that T-type channels contribute to setting the SAN firing rate [27], and that high variability exists upon target inactivation of T-type Ca v 3.1 channels [28]. Thus, T-type and funny channels may be important stabilizers of heart rhythm and may have a synergistic impact on variability increase when both are disabled. Taken together, it is likely that both I f and I CaT play a significant role in pacemaker synchronization.
Voltage clamping is the conventional method used to study ion channels [29]. This method is useful for measurement of the direct effect of drugs on specific channel populations. However, because ion channels are coupled through changes in voltage and other internal signals, this method is not suitable for measuring the dynamic feedback between the cell parameters.
SAN bradycardia is often associated with a shift in the leading pacemaker, from the superior SAN inferiorly towards the subsidiary atrial pacemakers, including the atrioventricular node [30]. In this work, a cell in the central SAN area was imaged (see Figure 1A), and the same cell was traced before and after drug perturbation, thereby circumventing the impact of a potential shift in the leading pacemaker on our conclusions.
Two components of I Na have been reported in the mouse SAN: TTX-resistant and TTX-sensitive Na + channels. It was recently shown [31] that TTX-sensitive channels may contribute to SAN pacing, while TTX-resistant I Na is likely to be responsible for AP propagation from the SAN to the atrium. Note that we did not use cells in the periphery and only cells that responded to TTX were used. Moreover, at transmembrane potentials reported for mouse SAN cells [32], TTX-resistant Na + channels are expected to be inactive and non-contributing to pacemaker activity.
This work focused on inhibition of HCN 4 by IVA. However, 20% of HCN channels in the mouse are HCN 2 [33]. However, specific inhibitors of HCN 2 have not been tested on SAN cells yet and, thus, its potential interaction with other channels cannot be tested here.

Limitations
Aside from its inhibitory effect on voltage-dependent sodium channels, TTX has also been shown to inhibit Ca 2+ currents [34]. However, this inhibition requires TTX doses higher than those used in this work. If it did inhibit Ca 2+ channels, then the magnitude of the synergy between I f and I Na would be expected to be higher since the I Na effect is mediated by increased activity of Ca 2+ channels.
High doses of TMR have been found to reduce I CaL together with I CaT [17]. While high TMR doses were not applied here, some nonspecific drug effects may have occurred.
Differences in the values of global parameters in the model compared to the experiments may have arisen due to use of SAN tissue in the experiments versus single SAN cells in the model. They also may be ascribable to the technical limitations of the model to generate extremely short or long BIs, and to the lack of certain cell mechanisms, such as metabolic pathways, bioenergetics, and BIV, in the models. Yet, the experimental and computational trends were similar.
In this work, BI was calculated from Ca 2+ measurements of SAN tissue that was treated with fluo-4. It is known that when compared to patch or external electrodes, this technique yields prolonged BI in SAN tissue and cells. Yet, our BIs were in the range of previous publications of mouse SAN tissue measurements [6,[35][36][37]. Others have published shorter BIs [38,39], but based on their illustration, the tissue was stretched, which could shorten the BI [40]. Moreover, their tissue had multiple rhythms while we only used tissues with synchronically beating cells. Note that since blebbistatin or any other drug that eliminates contraction was not used, we were able to clearly see the spontaneously beating area and to determine if there was more than one rate of the pacemaker (multiple sites).
Note that I Na inhibition may also suppress SAN-to-atrium propagation. However, this suppression is not affected by I f . The experiments showed that when both I f and I Na were blocked, the effect on the spontaneous BI was higher than the summed effect of each blocker. Thus, synergy must exist between the currents that are involved in SAN-to-atrium propagation.
Because Na 1.5 is TTX-resistant [41], a higher concentration of TTX is needed to block these channels. However, 10 µM TTX completely stopped the SAN electrical activity in our experiments.
Pacemaker shifts can occur with changes in BI [42]. However, we found no movement of the cell before and after drug perturbation and the same cell was traced before and after treatment. Moreover, we did not use blebbistatin or any other drug that eliminate contraction, which enabled us to clearly see the pacemaker contraction. Note that TMR by itself prolongs BI to a similar degree as TTX. When TMR was applied together with IVA, the effect on the spontaneous BI and LCR period was similar to the summed effect of each blocker (IVA or TMR). Thus, the prolongation of BI by one drug does not necessarily lead to a pacemaker shift that leads to a non-linear effect on BI when combined with another drug.

SAN Confocal Ca 2+ Imaging
To measure Ca 2+ signals, intact SAN preparations were loaded with fluo-4-AM (Ther-moFisher, 30 µmol/L) over 1 h, at 37 • C, on a shaker set to 60 RPM. The tissues were washed twice with Tyrode solution before imaging. Ca 2+ fluorescence was imaged using an LSM880 confocal laser scanning microscope (Zeiss) with a 40×/1.2 N.A. water immersion lens ( Figure 1B). Baseline recordings were performed after 30 min rest at 37 • C. Tissues were excited with a 488 nm argon laser and emission was collected with a low-pass 505 nm filter. Images were acquired in line scan mode (1.22 ms per scan; pixel size, 0.01 µm) along the pacemaker cells. The same cell was imaged before (baseline) and after drug perturbation(s). Each recording lasted at least 3 s.

Ca 2+ Analysis
Ca 2+ signaling was analyzed using a modified version of the software "Sparkalyzer" [43]. The fluorescence signal (F) was normalized by the minimal value between beats (F 0 ). Ca 2+ transients were semi-automatically detected and Ca 2+ sparks were manually marked. BI was calculated as the average time between Ca 2+ transient peaks, and the BIV was calculated as its standard deviation. The Ca 2+ transient amplitude, Ca 2+ transient 50% relaxation time, and LCR period were automatically calculated by the software as described before ( Figure 1C) [43].

Drugs
I f was blocked with 3 µmol/L IVA (Toronto Research Chemicals, North York, ON, Canada). I Na was blocked with 5 µmol/L TTX (Alomone Labs, Jerusalem, Israel). I CaT was blocked with 10 µmol/L TMR (Sigma-Aldrich, Saint Louis, MO, USA). All drugs were initially dissolved in dimethyl sulfoxide (DMSO). Images were recorded at least 5 min after TTX or TMR application, and at least 15 min after IVA application.

Statistics
Experimental results are presented as mean ± SEM. Statistical comparisons between baseline and post-treatment were performed with one-way ANOVA and paired t-tests. Differences were considered statistically significant at p < 0.05. In each experiment, N refers to the number of SAN preparations.

Computational Model
To investigate the internal dynamics of SAN cells treated with different blockers, we used our previously published computational model of the mouse SAN cell [13], which itself based on previous publications [5,[44][45][46][47][48][49][50][51][52]. The model contains 43 state variables and differential equations describing the dynamics of cell membrane potential, ion currents, Ca 2+ cycling, and post-translational modifications during an AP ( Figure 1D). The model initial conditions (Table S1) and constants (Table S2) were based on the experimental results. The IVA effect was simulated by reducing the I f maximal conductance coefficient, g If , to 10% of its value ( Figure 1E). The TTX effect was simulated by reducing g Na11 , the TTX-resistant channel maximal conductance coefficient, to 80% of its value and g Na15 , the TTX-sensitive channel maximal conductance coefficient, to 50% of its value ( Figure 1F). Current clamps were simulated by bringing the specific current to its basal state by an alternation of its maximum conductance coefficient. The software was run in MATLAB (The MathWorks, Inc., Natick, MA, USA). Numerical integration was performed using the MATLAB ode15s stiff solver, and the model simulations were run for 200 s to ensure that a steady state was reached. Computation was performed on an Intel(R) Core(TM) i7-4790 CPU @ 3.60 GHz machine with 8 GB of RAM.