Caveolar Compartmentalization of Pacemaker Signaling is Required for Stable Rhythmicity of Sinus Nodal Cells and is Disrupted in Heart Failure

Background: Heart rhythm relies on complex interactions between electrogenic membrane proteins and intracellular Ca2+ signaling in sinoatrial node (SAN) myocytes; however, mechanisms underlying the functional organization of proteins involved in SAN pacemaking and its structural foundation remain elusive. Caveolae are nanoscale, plasma membrane pits that compartmentalize various ion channels and transporters, including those involved in SAN pacemaking, via binding with the caveolin-3 scaffolding protein, but the precise role of caveolae in cardiac pacemaker function is unknown. Our objective was to determine the role of caveolae in SAN pacemaking and dysfunction (SND). Methods: Biochemical co-purification, in vivo electrocardiogram monitoring, ex vivo optical mapping, in vitro confocal Ca2+ imaging, and immunofluorescent and electron microscopy analyses were performed in wild type, cardiac-specific caveolin-3 knockout, and 8-weeks post-myocardial infarction heart failure (HF) mice. SAN tissue samples from donor human hearts were used for biochemical studies. We utilized a novel 3-dimensional single SAN cell mathematical model to determine the functional outcomes of protein nanodomain-specific localization and redistribution in SAN pacemaking. Results: In both mouse and human SANs, caveolae compartmentalized HCN4, Cav1.2, Cav1.3, Cav3.1 and NCX1 proteins within discrete pacemaker signalosomes via direct association with caveolin-3. This compartmentalization positioned electrogenic sarcolemmal proteins near the subsarcolemmal sarcoplasmic reticulum (SR) membrane and ensured fast and robust activation of NCX1 by subsarcolemmal local SR Ca2+ release events (LCRs), which diffuse across ~15-nm subsarcolemmal cleft. Disruption of caveolae led to the development of SND via suppression of pacemaker automaticity through a 50% decrease of the L-type Ca2+ current, a negative shift of the HCN current (If) activation curve, and a 40% reduction of Na+/Ca2+-exchanger function, along with ~2.3-times widening of the sarcolemma-SR distance. These changes significantly decreased the SAN depolarizing force, both during diastolic depolarization and upstroke phase, leading to bradycardia, sinus pauses, recurrent development of SAN quiescence, and significant increase in heart rate lability. Computational modeling, supported by biochemical studies, identified NCX1 redistribution to extra-caveolar membrane as the primary mechanism of SAN pauses and quiescence due to the impaired ability of NCX1 to be effectively activated by LCRs and trigger action potentials. HF remodeling mirrored caveolae disruption leading to NCX1-LCR uncoupling and SND. Conclusions: SAN pacemaking is driven by complex protein interactions within a nanoscale caveolar pacemaker signalosome. Disruption of caveolae leads to SND, potentially demonstrating a new dimension of SAN remodeling and providing a newly recognized target for therapy.

Heart failure mouse model was done in the cardiovascular physiology core facility at the University of Wisconsin Madison.HF was induced by 8-weeks of myocardial infarction (MI) in C57BL/6J male mice by ligating the left anterior descending coronary artery as previously described (17).Specifically, mice were anesthetized and ventilated with 2% isoflurane and oxygen.Followed by a left thoracotomy, the left anterior descending coronary artery (LAD) was ligated with 8-0 nylon suture.Successful acute myocardial ischemia was confirmed when left ventricular wall became pale.Mice then were kept on heating pad until awake.
HF progression was verified by echocardiographic analysis at 8 weeks post-MI.Parameters including ejection fraction (EF) and left ventricular (LV) mass are shown in Data Supplement Figure S17.Specifically, systolic and diastolic function was assessed by echocardiographic examination in lightly anesthetized (1.5% isoflurane) mice using a 17.5-MHz transducer (RMV 707B; Visual Sonics, Toronto, ON, Canada), as previously described (18).Two-dimensional M-Mode images were recorded both in the short and long axes.Average values were obtained from the measurement of 3-5 consecutive cardiac cycles.

Caveolin-binding motif (CBM) analysis
Proteins' CBMs were identified based on the following motifs: ΦXXXXΦXXΦ or ΦXΦXXXXΦ where Φ is an aromatic amino acid and X is a nonaromatic amino acid (Data Supplement Figure S8).In brief, Python was used to convert protein amino acid sequences into binary and scanned for the aforementioned motifs.The python code is: # first string firstString = "FWHYARNDCQEGILKMPSTV" secondString = "11112222222222222222" string = "input amino acid sequence here" print("Original string:", string) translation = string.maketrans(firstString,secondString) # translate string print("Translated string:", string.translate(translation))

In vivo ECG Recording
Surface ECG recordings were done as previously described (19).Briefly, mice were lightly anesthetized with isoflurane vapor titrated to maintain the lightest anesthesia possible.On average, 1.5% volume isoflurane vapor was required to maintain adequate anesthesia.Loss of toe pinch reflex and respiration rate was used to monitor levels of anesthesia.Average respiration rate was not different between mouse groups.ECG was recorded for 5 minutes, and heart rate was measured as the average over a 30-second interval.
For isolated SAN preparations, after cannulation ventricles were dissected away, and the atria were flattened and then pinned to the bottom of a Sylgard-coated chamber and superfused with Tyrode solution at a constant rate of ~15 ml/min.The medial limb of the crista terminalis (CT) was cut in order to open the right atria appendage (RAA), and the pacing electrode was placed on the edge of the RAA.Two Ag/AgCl electrodes were immersed into the superfusion solution and placed on the sides of preparation to document pseudo-ECG.The inter-atrial septum (IAS) tissue was partially removed to reduce scattering of the optical signal from tissue that was not in focus.A rim of ventricular tissue was preserved for pinning the preparation in order to prevent damage to the atria.Both the left and right atria (LA and RA) as well as the atrio-ventricular junction (AVJ) were accessible in the preparation.

Imaging system
After isolation, the SAN preparations were then immobilized by infusion of Blebbistatin (10 μM, Tocris Bioscience, USA) in the perfusion media in order to suppress motion artifacts in optical recordings.Then, the preparations were co-stained with voltage-sensitive dye RH-237 (ThermoFisher Scientific, USA) (1.25 mg/ml in dimethyl sulfoxide) by direct application of the dye on the tissue.Excitation light (520/44 nm) was generated by a 150-W halogen lamp with a constant-current, low-noise, power supply (MHAB-150W, Moritex USA Inc., CA, USA).A flexible light guide directed the band-pass-filtered light onto the preparation and a shutter was used to ensure that the preparation was exposed to light only during image acquisition.The fluorescent light emitted from the preparation was long-pass (>715 nm) filtered using an edge pass filter (Thorlabs, NJ) before reaching the camera.Emitted light was directed towards a MiCAM Ultima-L CMOS camera (SciMedia, CA, USA) with high spatial (100 x 100 pixels, 60 ± 10 μm per pixel) and temporal (1,000 -3,000 frames/sec) resolution.The acquired fluorescent signal was digitized, amplified, and visualized using custom software (SciMedia, CA, USA).

Experimental protocol
After isolation, motion suppression and staining, preparations were equilibrated for 20 min before imaging.Next, control maps of atrial activation during spontaneous rhythm were made as previously described (21).Location of the leading pacemaker was visualized and attributed to one of the four anatomical regions as previously described based on the location and the presence of HCN4 and connexin43 protein expression (19): (1) primary SAN (as characterized by HCN4positive and Cx43-negative expression pattern), (2) lower SAN (as characterized by both HCN4and Cx43-positive expression pattern and located along the CR between the SVC and IVC), (3) intra-atrial septum (IAS), and (4) latent pacemakers (as characterized by both HCN4-and Cx43positive expression pattern and located outside of SAN region and ICR).One mouse could have several sites of the leading pacemaker location because of the presence of multiple competing pacemakers during the optical recording.
Heart rhythm stability was calculated as shown in Figure 2 -average CL and CL lability.Average CL was calculated for 150 continuous beats per mouse at baseline condition after a 30-min of stabilization period.CL lability is calculated as the standard deviation of 150 continuous beats.
To estimate SAN recovery time (SANRT), the preparation was paced (S1S1 = 100 ms) for at least 30 secs through a pacing electrode located at the RA appendage (22).SANRT was measured as the time interval between the last pacing beat and the first spontaneous beat.Corrected SANRT (SANRTc) was calculated as the difference between the SANRT and the resting cycle length (CL) measured before the SANRT pacing protocol.

Optical mapping data processing
A customized Matlab-based computer program was used to analyze the optical signals (20,21,23).Signals were filtered using the low-pass Butterworth algorithm at 256 Hz.Maximum upstroke derivative (dV/dtmax) was calculated for each action potential using the normalized optical signal and its derivatives.Activation maps were constructed from activation times, which were determined from the dV/dtmax.

SAN Cell Isolation
Single SAN cells were isolated from mice using a modified method based on previously published protocols (19,(24)(25)(26).Cells exclusively from the region which corresponds to the primary pacemaker area characterized previously in the mouse heart by histology and immunolabeling of connexin43, connexin45, and HCN (19,27,28) as well as functionally by microelectrodes (29,30) and optical mapping (19)(20)(21)31).Briefly, after heart isolation and cannulation, the isolated SAN preparations were prepared as described above.Subsequently, the SAN region, bordered by the crista terminalis, atrial septum, and orifice of superior vena cava, was dissected from the heart.The SAN area was cut into strips which then was washed three times with 2 mins each time in the 'low Ca 2+ , Mg 2+ free' solution (in mmol/L: 140 NaCl, 5.4 KCL, 1.2 KH2PO4, 0.066 CaCl2, 50 taurine, 18.5 D-glucose, 5 HEPES and 1mg/ml bovine serum albumin (BSA), with pH adjusted to 6.9 with NaOH).The SAN tissue pieces were digested for 30 mins at 35°C in a cocktail of 1mg/ml Collagenase (Sigma) and 0.2mg/ml Elastase (Worthington) with gentle shaking every 5 min.The digested tissue was carefully rinsed with Kraft-Brühe solution (in mmol/L: 100 potassium glutamate, 10 potassium aspartate, 25 KCL, 10 KH2PO4, 2 MgSO4, 20 taurine, 5 creatine, 0.5 EGTA, 20 glucose, 5 HEPES, and 0.1% BSA, with pH adjusted to 7.4 with KOH) and gently triturated with a pair of fire polished, wide-bore Pasteur pipettes to release SAN cells.Isolated cells were then plated on Laminin (Sigma) coated coverslips for 30 mins, and readapted to normal extracellular Ca 2+ concentration of 1.8 mmol/L.SAN myocytes were identified by their small spindle shape and ability to beat spontaneously in the recording chamber when superfused with normal Tyrode's solution.We confirmed pacemaker phenotype in isolated cells by HCN4-positive and connexin43-negative immunofluorescence staining (in contract to HCN4-negative and connexin43-positive atrial cells) as shown in Figure 2A.Finally, HCN4-positive expression in studied SAN pacemaker cells is further confirmed by IF staining Figures 1C and 4D, Western blot (Figure 1E), RT-qPCR (Figure 6A), patch clamp recording of If current (Figure 6B) and sensitivity to a selective If current inhibitor Ivabradine (Figure 7A) as well as unique morphological organization on TEM photographs, including tack of transversal tubules, disorganized sarcomere organization and a specific pattern of mitochondria distribution (Figure 1A and Figure 4A, B vs. TEM photographs of atrial and ventricular myocytes shown in Supplementary Figure 1).

Intracellular confocal Ca 2+ recordings
SAN cells were loaded with the Ca 2+ indicator dye Fluo4-AM (Invitrogen) for 15 mins and then were washed with Tyrode's solution (in mmol/L: 140 NaCl, 5.4 KCL, 1.2 KH2PO4, 1 MgCl2, 5.55 D-glucose, 5 HEPES) for 15 mins at room temperature (22 ± 2°C).SAN cells with characteristic spindle and/or spider-like morphologies showing sustained spontaneous activity were selected for confocal Ca 2+ imaging using Leica SP5 confocal microscopy system following previously published methods (18,32).Confocal line scans were applied along the longest straight membrane of the cells to capture the local Ca 2+ release events with 2,000 Hz scanning speed under 63×/1.40NA oil immersion objective, pinhole of 1 airy unit.

Patch clamp studies
ICa,L and ICa,T: For calcium current recordings, electrophysiological recordings were carried out in the whole cell configuration of the patch-clamp technique at room temperature using the Axopatch 200B amplifier (Axon Instruments, Foster City, CA) with pCLAMP 10.7 software.Recording pipettes were pulled from thin-walled borosilicate glass capillaries (World Precision Instruments, Inc., Sarasota,FL) with pipette resistance of 3-5 MW.The recordings were filtered at 2 kHz and digitized at 20 kHz.ICa currents were recorded in a voltage clamp mode of the patchclamp technique.All solutions and buffers are indicated in mM/l.The external bath solution contained 120 Tetraethylammonium-Chloride, 10 CsCl, 10 Glucose, 10 HEPES, 1.5 MgCl2, 1 CaCl2, and pH 7.4 (CsOH).Intracellular pipette solution contained 100 Cs-methansulfonate, 30 CsCl, 10 HEPES, 5 EGTA, 5 Mg-ATP, 2 MgCl2, and pH 7.2 (CsOH).Whole cell currents were recorded from a holding potential -80 mV with 300-ms test pulses from -40 to 50 mV, in 10 mV increments for ICa.For ICa,L a holding potential was -40mV.We calculated ICa,T by subtracting ICa,L from ICa.Data were analyzed using Microcal Origin software (Origin Lab Corporation Northampton, MA).

If:
If current was recorded from spontaneously beating single cells superfused with a modified Tyrode solution containing (mM): NaCl, 115; KCl, 25; CaCl2, 1.8; MgCl2, 1; MnCl2, 2; BaCl2, 2; Hepes, 5; and glucose, 10 (pH 7.4, adjusted with NaOH).Isoproterenol was not added.Recordings were acquired using an Axopatch 200A amplifier (Molecular Devices), attached to a Digidata 1440A data acquisition system (Molecular Devices), and acquired using Clampex 10.0 at 20 KHz with a low-pass filter of 5 KHz.Electrodes were fabricated with patch glass (Warner Instruments) using a Sutter Instrument micropipette puller P-97 and were filled with an internal solution containing (mM): Potassium aspartate, 130; NaCl, 10; CaCl2, 2; Hepes, 10; EGTA-KOH, 5; MgCl2, 2; NaATP, 2, and NaGTP, 0.1 (pH 7.2, adjusted with KOH).The liquid junction potential resulting from the external and internal solutions used was 14 mV and was corrected after the recordings.Pipettes had resistances of 3-7 MW, before series resistance compensation.Steadystate current amplitudes were calculated at the end of a series of 3 s steps from -44 mV to -154 mV in 5 or 10 mV increments from a holding potential of -14 mV.Voltage steps were followed by a 1 s step to -114 mV to elicit tail currents and assess the voltage dependence of the activation.Steps were given at 10 s intervals.Some cells received 100 µM ZD 7288 to block HCN currents which eliminated most of the currents.No leak or blocker subtraction were used.The mean amplitudes of the tail currents at the end of the -114 mV step were plotted against the activation potential.Curves were fitted to a Boltzmann function, I = A2 + [(A1-A2)/(1+e (V-V1/2)/z )], where A2 is the maximum tail current amplitude, A1 is the current offset, V1/2 is the midpoint of activation, and z is the slope, using Origin 9.0 (OriginLab Corporation).
Fully activated If current-voltage (I/V) relationships were obtained according to a previously published protocol (33).

Isolated cardiomyocytes immunofluorescence labeling
Immunolabeling was performed as previously described (19,34) on mouse SAN cells using the following primary antibodies (Data Supplement Table I).Briefly, isolated cardiomyocytes were fixed with 4% buffered paraformaldehyde or ice-cold methanol for 10 min.Cells fixed with paraformaldehyde were permeabilized with Triton X-100 (0.1%) for 10 min.After washing with PBS containing 0.05 % Tween-80 (PBS-T) (three 5-min washes), cells were incubated with 1 ml blocking solution (2% BSA and 2% goat serum in PBS-T) for 1.5 h at room temperature to block nonspecific binding.Subsequently, cells were incubated overnight with respective primary antibodies in blocking solution at 4°C.Excess primary antibody was washed off with the use of blocking solution (three 5-min washes).The cells were then incubated overnight with Alexaconjugated secondary antibodies (Molecular Probes, Eugene, OR; 2 mg/ml) diluted 1:800 in blocking solution.Highly cross-absorbed goat anti-mouse Alexa Fluor 488, goat anti-rabbit Alexa Fluor 568 and goat anti rat Alexa Fluor 647 were used.The cells were then washed with PBS-T (three 5-min washes), and mounted on a coverslip.To determine nonspecific binding, control experiments with secondary antibody alone were also performed.Imaging was performed with a Leica SP5 Confocal microscope under 63×/1.40NA oil immersion objective, pinhole of 1 airy unit.A sequential imaging pattern were applied during double staining fluorescent imaging.
To further confirm the close location of Cav1.2, Cav-3 and HCN4 in one complex triple staining experiments were performed.Labeled cells were incubated with highly cross-absorbed goat antimouse Alexa Fluor 488, goat anti-rabbit Alexa Fluor 568 and goat anti-rat Alexa Fluor 647 were used.Imaging was performed similarly as for double stained cells.

Co-localization analysis
Colocalization plots were used to determine the amount of different proteins that colocalized with cav-3 in the SAN cell membrane.In a 2-channel confocal image, each voxel has two intensity values (ranging from 0 to 1 after normalization in a 12-bit image), one for each red and green fluorescent signals.Voxels with normalized intensities <0.3 were considered background fluorescence and were excluded from colocalization analysis.A colocalization plot, generated with a custom-made Matlab-based computer program, displays these intensity values as a function of each other channel.By definition, two proteins are highly colocalized in a particular volume when fluorescence intensities corresponding to these two proteins are high in the voxel corresponding to this volume (1,18).Therefore, if two proteins are colocalized in many voxels, the colocalization plot will contain a significant diagonal distribution.Voxels with the highest degree of colocalization will be displayed in the upper right quadrant.In contrast, if the two proteins are not colocalized, the colocalization plot shows voxel values near each axis, with no diagonal elements present.The analysis was complemented with determination of percentage of colocalized voxels using BlobProb plugin in Fiji (http://www.fiji.sc)(35).Thresholds for each channel were determined as described above.

Proximity ligation labeling
Proximity ligation assay (PLA) was performed on methanol-fixed cardiomyocytes with the Duolink (inSitu) kit (Sigma-Aldrich).Cells were incubated overnight with primary antibodies against NCX and Cav-3.Secondary antibodies were PLA anti-Rabbit-plus and PLA anti-mouseminus.PLA ligation and amplification steps were performed with far-red Duolink PLA kit according to manufacturer's instructions (Sigma-Aldrich, United States).Cell imaging was performed with 60x oil immersion Leica SP5 confocal microscope using 63x/140 NA oil immersion objective.Cells were illuminated by 633 helium-neon laser and detection was performed within 640-700 nm range.PLA analysis was performed by automatically thresholding cell images in FIJI (http://www.fiji.sc)and applying build-in "Analyze particles" plugin in FIJI to count the number of particles per area.

Co-immunoprecipitation
Mouse and human SAN tissue (pre-validated by immunostaining of HCN4 positive and Cx43 negative area as shown in Data Supplement Figure S2 (36,37)) lysates as well as HEK293 cell line were used, and immunoprecipitations were carried out using anti-Cav-3 or anti-NCX, or anti-HCN4 antibodies as previously described (38); negative control of b-actin and phospholamban.Briefly, tissue samples from mice and human and HEK293 cells were flash frozen in liquid nitrogen, homogenized and lysed using buffer containing 150 mM NaCl, 25 mM Tris⋅HCl, 10 mM NaEGTA, 20 mM NaEDTA, and supplemented with 0.5% CHAPS and 1x protease inhibitor cocktail (B14001, bimake).The lysate (500 µg) was sonicated, rotated at 4ºC for 2 hours and centrifuged at 16,000 × g for 15 min at 4ºC.The soluble supernatant was incubated with aforementioned antibodies (20 µg) at 4ºC for 2 hours and then incubated with 25 µL (Protein A Mag Sepharose/Protein G Mag Sepharose, 28944006 /28944008, GE Life Sciences; 1:1 Protein A/G) at 4ºC for 2 hours, followed by incubation with anti-Cav-3, anti-NCX and anti-HCN4 (antibody details see Data Supplement Table I) in a total of 300-400 µL (depending on protein concentration) of lysate.After flowthrough collection, samples were eluted with 2.5% Acetic Acid.Laemmli Buffer (161-0747, Bio-Rad) and DTT were added and all samples were incubated at 65ºC for 10 min.Immune complexes were analyzed by SDS-PAGE (4-20% gradient gels, Bio-Rad) and Western blot by probing with appropriate antibodies.

Transmission Electron Microscopy
SAN tissue was isolated as previously described in SAN cell isolation section and fixed in the following fixative: 2.5% glutaraldehyde, 2.0% paraformaldehyde, and 0.2 mol/l cacodylate buffer for 24 -48 hours.The samples were rinsed in the same buffer, post-fixed in 1% osmium tetroxide, dehydrated in a graded ethanol series, rinsed in propylene oxide, and embedded in Epon 812 substitute.After resin polymerization, the samples were then sliced into 70-nm sections with a Leica EM UC6 ultramicrotome and placed on 200 mesh transmission electron microscopy grids.The samples were post-stained in 8% uranyl acetate in 50% EtOH and Reynold's lead citrate, viewed on a Philips CM120 transmission electron microscope, and documented with a SIS MegaView III digital camera.
Electron microscopy images were analyzed by using the NIH ImageJ software.A threshold size for individual caveolae was set between 50 and 100 nm.The number of caveolae was counted as per unit length (μm) of myocyte sarcolemmal membranes from a series of random electron microscopy micrographs.Images with recognized sarcolemmal reticulum (SR) from WT, HF and Cav-3KO SAN tissue were selected and the distance between SR membrane to the closest sarcolemmal membrane was calculated.

Statistical Analysis
Student t test was used in 2-group comparisons.Multiple groups of normally distributed data of similar variance were compared by 1-or 2-way ANOVA.For multiple comparisons, the Bonferroni corrected P value is shown.All statistical analyses were performed using GraphPad Prism 5 (GraphPad Software).P<0.05 was considered statistically significant.Values were presented as mean ± SEM.

Mathematical modeling and simulation
Subcellular Ca 2+ signaling model.We used an established 3D model of rabbit ventricular Ca 2+ signaling (39,40) as the basis for developing an analogous model of the murine SAN cell that describes subcellular stochastic properties of individual SR Ca 2+ release units (CRUs).Figure 5A depicts the model schematic: in each peripheral CRU, i.e., a CRU coupled to external membrane, there are 5 Ca 2+ compartments (cytosolic, submembrane, cleft space, network SR, and junctional SR), whereas central CRUs are not coupled to external membrane and have less RyRs vs. the peripheral CRUs (Data Supplement Table II).We assigned the first two CRUs close to the surface membrane as peripheral CRUs and the rest as central CRUs.CRUs are coupled by Ca 2+ diffusions in submembrane, cytosol and the network SR.We modified this model to describe the structural data, global cellular dimensions, and electrical capacitance of typical mouse sinoatrial node cells based on (41,42).The new 3D SAN cell model has a dimension of 60 μm × 7 μm × 7 μm with a capacitance of 25 pF and comprises 2176 (34 × 8 × 8) CRUs.Detailed model updates and parameters are in the Appendix and Data Supplement Table II.Briefly, we updated Ca 2+ current models in the original 3D rabbit Ca 2+ signaling model to account for the mouse SAN cellspecific properties (Data Supplement Table II).Specifically, we replaced the model of L-type Ca 2+ current (ICaL) in the original 3D Ca 2+ signaling model with a Markov-equivalent implementation of the ICaL formulation from (41) using the approach employed in (43).The Na + /Ca 2+ exchanger current (INCX) and T-type Ca 2+ current (ICaT) of the new model were simulated using the formula in (41).Furthermore, we updated the RyR release model to reflect enhanced SAN RyR Ca 2+ sensitivity compared to the original model for rabbit ventricle.Our model does not explicitly account for functional caveolar or extra-caveolar Ca 2+ compartments, which would be intermediate (in terms of size and Ca 2+ concentration levels) between the cleft and subsarcolemmal spaces.Localizations of L-type Ca 2+ channels, T-type Ca 2+ channels and NCX were based on colocalization data collected in Figure 1.
Coupled Ca 2+ and Vm model.To build a 3D mouse SAN myocyte model incorporating detailed local Ca 2+ signaling and electrophysiology, we coupled the new 3D subcellular Ca 2+ signaling model with a mouse SAN AP model (0D AP model).The 0D AP model was updated from (41) as previously described in (42) (details in the Appendix).For each time step, all Ca 2+ -dependent transmembrane currents (ICaL, ICaT, INCX, Sarcolemmal membrane background Ca 2+ current and Ca 2+ pump -ICaBk and ICaP) and intracellular Ca 2+ handling were calculated using the 3D Ca 2+ signaling model, and the remaining transmembrane currents (Na + and K + currents) were computed using the 0D AP model.The membrane voltage (Vm) was then computed from the total transmembrane current using the electrical circuit model as in (41,42).The intracellular Na + concentration was fixed to 10 mM as in (39).Of note, spatial models of rabbit SAN myocyte incorporating stochastic descriptions of Ca 2+ handling have been previously constructed, reveal key contributions of subcellular Ca 2+ release units and their coupling with sarcolemmal ionic proteins to maintaining normal and robust pacemaking activity (44)(45)(46)(47).
Simulation protocols.Our new 3D SAN cell model was implemented in C++.The ordinary differential equations are solved using the forward Euler method except that ICaL and RyR flux were solved stochastically as in (39,40,43).The time step was 0.01 ms as in (39).A period of 30 sec was simulated and the data from the last 20 sec was analyzed.

0D mouse SAN AP model
We used an updated version of the Kharche et al. model (41) model of mouse sinoatrial node myocyte (SAM), as described in (42) based on experimental data obtained in the Proenza lab in 2-3 months old male C58BL/6J mice at physiological temperature.In the updated model, the parameters for ICaL, If, Ito and Isus were initially manually tuned to match our voltage-clamp experiments; then, an optimization based on reverse multivariable regression (nonlinear iterative partial least squares method, as in ( 48)) was performed to identify model parameters (Data Supplement Table II) that best matched various recorded AP waveform parameters (including cycle length, APD90 and APD50, AP amplitude, diastolic depolarization rate and duration, upstroke velocity, maximum diastolic potential).
Modeling scheme of L-type Ca 2+ current.Left panel: a Hodgkin-Huxley scheme of ICaL modified from (41) with the addition of dl3 state; Right panel: the Markov-equivalent scheme for the Hodgkin-Huxley ICaL model.CDI -Ca 2+ -dependent inactivation.VDI -voltage-dependent inactivation.The conceptualization of conversion between the two schemes is detailed in (43).

L-type Ca 2+ current:
The L-type Ca 2+ channels are assumed to be located in the cleft area in the peripheral CRUs and closely coupled with RyRs.We modified the Hodgkin-Huxley type model of ICaL in (41) to a Markov-equivalent scheme using the approach described in (43).To reflect the difference in maximum open probability between the Hodgkin-Huxley type model (close to 100%) and experimental measurements (~5%), we introduced an additional state (dl3) in the Markovequivalent scheme above as detailed in (43).The model is then solved stochastically as in (43).The ICaL current is described by Ca 2+ buffering: Similar to the previous studies (39,40), the buffering of Ca 2+ in each compartment was modeled using an instantaneous buffering function, with an exception that troponin-Ca 2+ was modelled dynamically.In our 3D model of SAN, we removed buffering of SR, Myosin (Ca 2+ ), and Myosin (Mg) in the cleft and submembrane compartments as these buffers were absent in these compartments in the previous models (41,49).

type Ca 2+ current (ICaT): The
,,* ⋅  where PO is a binary variable describing whether a LTCC channel is open or closed, NL,O is the total number of open LTCC, and ica is the unitary current and computed as  = 4 ()  ⋅  !.⋅0.001 ⋅ exp(2 ⋅ ) −  ICaT formula from the Kharche et al. model (41) was coupled with the 3D Ca 2+ signaling model.We evenly distributed the Ca 2+ flux carried by ICaT to the peripheral CRUs, with 50% of the proteins in the cleft and 50% in the submembrane compartment for each CRU (Data * [] *

Supplement Table II). Na + /Ca 2+ exchanger (INCX):
(41)NCX model of the 3D Ca 2+ signaling model was replaced with that from(41).In the new baseline 3D Ca 2+ model of SAN, NCX was equally distributed among the peripheral CRUs, with 50% of the proteins in the cleft and 50% in the submembrane compartment for each CRU (Data Supplement

Table II). SR Ca 2+ release and leak:
(39)49)late the Ca 2+ release from SR via RyR, we developed a new model of RyR Ca 2+ release (schematic at right) based on the RyR models from(39,49).The new model adds [Ca]SR-dependence to closed-to-open transitions rates and higher sensitivity to cleft Ca 2+ compared with the Sato-Bers model(39).Also, in the new model, the SR Ca 2+ leak was redirected to the cleft space instead of the cytosolic space in(40).The new RyR model is described by ()/0 =  /0 −  /0 −  /0 5& =  *()/0  6  .(/:; and  (/:;2 are kept the same as in(39)

Table II). Data Supplement Table I. List of used primary antibodies.
(49)s Ca 2+ concentration in the cleft or submembrane space.The compartmentalization of ICaBk was kept the same as in (49) (Data Supplement TableII).Cx is the Ca 2+ concentration in the cleft or submembrane space.The compartmentalization of ICaP was kept the same as in(49).Parameters for both ICaBK and ICaP was tuned so that the two whole-cell currents from the 3D Ca 2+ signaling model were comparable to those from the updated 0D electrophysiological model (Data Supplement 2+ flux (ICaBk): Similar to Sato-Bers model (39), ICaBk was modeled as  ()<= =  ()<= ( −  () ) B where