# Computational modelling of mouse atrio ventricular node action potential and automaticity

Chiara Bartolucci<fnref n="1" /> ID, Pietro Mesirca<fnref n="2" />,<fnref n="3" />, Eugenio Ricci<fnref n="1" /> ID, Clara Sales-Bellés<fnref n="4" /> ID, Eleonora Torre<fnref n="2" />,<fnref n="3" /> ID, Julien Louradour<fnref n="2" />,<fnref n="3" /> ID, Matteo Elia Mangoni<fnref n="2" />,<fnref n="3" /> ID and Stefano Severi<fnref n="1" /> ID

### Handling Editors: Natalia Trayanova & Brian Delisle

The peer review history is available in the Supporting Information section of this article (https://doi.org/10.1113/JP285950#support-information-section).

<figure><img src="image_1.png" /></figure>

**Abstract** The atrioventricular node (AVN) is a crucial component of the cardiac conduction system. Despite its pivotal role in regulating the transmission of electrical signals between atria and ventricles, a comprehensive understanding of the cellular electrophysiological mechanisms governing AVN function has remained elusive. This paper presents a detailed computational model of mouse AVN cell action potential (AP). Our model builds upon previous work and introduces several key refinements, including accurate representation of membrane currents and

<footnote n="1">¹Computational Physiopathology Unit, Department of Electrical, Electronic and Information Engineering 'Guglielmo Marconi', University of Bologna, Cesena, Italy</footnote>

<footnote n="2">²Institut de Génomique Fonctionnelle, Université de Montpellier, CNRS, INSERM, Montpellier, France</footnote>

<footnote n="3">³LabEx Ion Channels Science and Therapeutics (ICST), Montpellier, France</footnote>

<footnote n="4">⁴BSICoS group, I3A Institute, University of Zaragoza, IIS Aragón, Zaragoza, Spain</footnote>
exchangers, calcium handling, cellular compartmentalization, dynamic update of intracellular ion concentrations, and calcium buffering. We recalibrated and validated the model against existing and unpublished experimental data. In control conditions, our model reproduces the AVN AP experimental features, (e.g. rate = 175 bpm, experimental range [121, 191] bpm). Notably, our study sheds light on the contribution of L-type calcium currents, through both $Ca_v$1.2 and $Ca_v$1.3 channels, in AVN cells. The model replicates several experimental observations, including the cessation of firing upon block of $Ca_v$1.3 or $I_{Na,r}$ current. $I_f$ block induces a reduction in beating rate of 11%. In summary, this work presents a comprehensive computational model of mouse AVN cell AP, offering a valuable tool for investigating pacemaking mechanisms and simulating the impact of ionic current blockades. By integrating calcium handling and refining formulation of ionic currents, our model advances understanding of this critical component of the cardiac conduction system, providing a platform for future developments in cardiac electrophysiology.

(Received 10 November 2023; accepted after revision 8 August 2024; first published online 13 September 2024)

**Corresponding author** S. Severi: Department of Electrical, Electronic and Information Engineering ‘Guglielmo Marconi,’ University of Bologna, Via dell’Università 50, 47521 Cesena (FC), Italy. Email: stefano.severi@unibo.it

**Abstract figure legend** The abstract figure illustrates the process of creating and validating a new numerical model for a single atrioventricular node (AVN) cell. The steps are as follows: (1) cells from the AVN were isolated from the mice hearts; (2) they underwent patch-clamp experiments to measure the L-type calcium current ($I_{Cav1.3}$); (3) experimental data and the new numerical model were fitted to generate the I-V (current-voltage) relationship; (4) the new numerical model was built incorporating various ionic currents and was able to simulate the electrical activity (membrane potential, $V_m$) and intracellular calcium dynamics ($Ca_i$) in the AVN cells; (5) the model was validated by comparing simulated action potentials and calcium transients under control conditions and after currents block and by generating a population of AVN models.

## Key points

* This paper introduces a comprehensive computational model of mouse atrioventricular node (AVN) cell action potentials (APs).

* Our model is based on the electrophysiological data from isolated mouse AVN cells and exhibits an action potential and calcium transient that closely match the experimental records.

* By simulating the effects of blocking specific ionic currents, the model effectively predicts the roles of L-type $Ca_v$1.2 and $Ca_v$1.3 channels, T-type calcium channels, sodium currents (TTX-sensitive and TTX-resistant), and the funny current ($I_f$) in AVN pacemaking.

* The study also emphasizes the significance of other ionic currents, including $I_{Kr}$, $I_{to}$, $I_{Kur}$, in regulating AP characteristics and cycle length in AVN cells.

* The model faithfully reproduces the rate dependence of action potentials under pacing, opening the possibility of use in impulse propagation models.

* The population-of-models approach showed the robustness of this new AP model in simulating a wide spectrum of cellular pacemaking in AVN.

Chiara Bartolucci is a Senior Researcher (RTD-A) in Biomedical Engineering at the Department of Electrical, Electronic and Information Engineering ‘Guglielmo Marconi’ - DEI - of the University of Bologna, Italy. She obtained a PhD in Biomedical Engineering in May 2015. The focus of her research is the development of computational cardiac models to understand the mechanisms of cardiac electrophysiology. She works on different cell types (e.g. ventricular, atrial, sinoatrial, and atrioventricular) and species (e.g. human, mouse, pluripotent stem cell).

<figure><img src="image_2.png" /></figure>
## Introduction

Since its discovery by Sunao Tawara and Karl Albert Aschoff in 1906 (Cavero & Holzgrefe, 2022), the atrioventricular node (AVN) (named, at that time, ‘Aschoff-Tawara node’) has been identified as a key component of the cardiac pacemaker-conduction system in mammalian hearts. The primary function of the AVN is to conduct the action potential from the atria to the ventricles, being normally the sole pathway for impulse propagation between them. Its slow conduction properties set the proper delay between atrial and ventricular activation and can also serve a protective function during some supraventricular tachyarrhythmias (Meijler & Janse, 1988). However, it also provides back-up pacemaker activity in case of failure of the sino-atrial node (Dobrzynski et al., 2003). The AVN is also an important pharmacological target for controlling ventricular rate upon atrial tachyarrhythmias. Despite its physiological importance, our current understanding of the cellular electrophysiological basis of the pacemaker activity of the AVN is still incomplete. In particular, besides well-known pacemaker mechanisms such as the ‘funny’ ($I_f$) current and ryanodine receptor 2 (RyR2) -dependent $Ca^{2+}$ release (Ridley et al., 2008), the L-Type $Ca_v$1.3, T-Type $Ca_v$3.1 calcium currents and the TTX-sensitive ($I_{Na,s}$, Nav1.1) and TTX-resistant ($I_{Na,r}$, Nav1.5) sodium currents have been shown to be important for cellular automaticity in the AVN (Marger et al., 2011a). Notably, the quantitative contributions of these additional ionic currents still need to be elucidated.

Computational models are well-established and powerful tools for the quantitative analysis of cardiac cell electrophysiology. However, although many AP models have been developed for most heart regions (Amuzescu et al., 2021; Trayanova et al., 2023), the AVN has not been widely studied. In 2009, the first biophysically detailed model of the AVN was published by Inada et al., based on experimental data from rabbit AVN (Inada et al., 2009). To identify and distinguish general and species-specific properties, however, it is desirable to have AVN cellular electrophysiology data from additional model species, accessible to gene targeting techniques. In 2011, Marger et al. (2011b) proposed a preliminary numerical model of mouse AVN cell automaticity, based on the mouse sinoatrial node (SAN) cell model (Mangoni et al., 2006; Zhang et al., 2000) to provide insights into the roles of L-type ($Ca_v$1.3), T-type, and $I_f$ currents. This model simulated the AP of isolated mouse AVN myocytes but did not integrate intracellular compartments preventing simulation of intracellular $Ca^{2+}$ and ionic dynamics along the pacemaker potential and the action potential upstroke and repolarization phases.

The purpose of this work was to obtain a more detailed, complete, and species-specific computational model of

the AP of mouse single AVN cells, by building from the Marger et al. model (Marger et al., 2011b) and fitting both old and previously unpublished experimental data. The introduction of calcium handling is paramount to generating a full set of AVN single-cell APs, due to the role that calcium has in cellular pacemaker activity. Implementation of cellular compartments and computation of intracellular $Ca^{2+}$ dynamics allowed us to create a second-generation model of mouse AVN pacemaking. The model reproduced experimental AVN AP and automaticity hallmarks. We used the model to simulate experimentally observed genetic deletion or blockade of ionic currents involved in pacemaker activity of AVN cells. In addition, we computed the rate-dependence of AVN APs by simulating the cell's stimulation at physiological pacing rates. Finally, we provide an in-depth characterization of model behaviour by applying the population-of-models approach.

## Methods

### Ethical approval

The study is in accordance with the *Guide for the Care and Use of Laboratory Animals* published by the US National Institute of Health (NIH Publication No. 85-23, revised 1996) and European directives (2010/63/EU). Experimental procedures for isolating AVN myocytes were approved by the Ethical Committee of the University of Montpellier and the French Ministry of Agriculture (protocol no: 30176) and comply with the policies of *The Journal of Physiology* regarding the use of animals in research Grundy (2015). Animals were housed in the Institute of Functional Genomics animal facility with free access to food and water and were exposed to 12 h light/dark reverse cycles (light, 08.00 to 20.00) in a thermostatically controlled room (21°C–22°C).

### Isolation of mouse AVN cells

AVN myocytes were isolated from 2- to 4-month-old wild-type, $Ca_v$3.1$^{-/-}$ (Kim et al., 2001), or $Ca_v$1.2$^{DHP-/-}$ mice, in which dihydropyridine (DHP) sensitivity in $Ca_v$1.2 $\alpha$1 subunits was eliminated by knock-in of a point mutation in the DHP binding site to the channel $\alpha$1 subunit (Sinnegger-Brauns et al., 2004).

Mice were deeply anesthetized by intraperitoneal injection of 0.3 ml of solution constituting ketamine (0.1 mg/g, Imalgène) and xylazine (0.01 mg/g, Rompun 2%, Bayer), followed by a second injection of pentobarbital (150 µL Euthasol Vet in 10 ml NaCl physiological solution). To avoid blood clots in the heart, we injected 0.5 ml of NaCl solution containing heparin (25 000 I.U.). Hearts were removed via thoracotomy when the tail sensitivity test became negative.
AVN pacemaker cells were isolated as previously described (Marger et al., 2011b). Briefly, the AV node region was excised in warmed (35°C) Tyrode solution containing (in mM): NaCl, 140; KCl, 5.4; CaCl₂, 1.8; MgCl₂, 1; Hepes, 5; and D-glucose, 5.5 (adjusted to pH 7.4 with NaOH). AV node tissue was then transferred into a 'low-Ca²⁺-low-Mg²⁺' solution containing (in mM): NaCl, 140; KCl, 5.4; MgCl₂, 0.5; KH₂PO₄, 1.2; taurine, 50; D-glucose, 5.5; bovine serum albumin (BSA), 1 mg/ml; HEPES, 5 (adjusted to pH 6.9 with NaOH). AV tissue was digested by Liberase TH (229 U/ml, Roche) and elastase (1.9 U/ml, Boehringer) for 22–24 min. Digestion was stopped by washing tissue strips in a modified 'Kraftbrühe' (KB) medium containing (in mM): L-glutamic acid potassium salt, 100; potassium aspartate, 10; KCl, 25; KH₂PO₄, 10; MgSO₄, 2; taurine, 20; creatine, 5; EGTA, 0.5; D-glucose, 20; BSA, 1 mg/ml; and HEPES, 5 (adjusted to pH 7.4 with KOH).

Single AV node cells were isolated by manual agitation using a flame-forged Pasteur's pipette, in KB solution at 35°C for 30–50 s. To recover the automaticity of the AVN cells, Ca²⁺ was gradually reintroduced in the cell's storage solution to a final concentration of 1.8 mM (Mangoni et al., 2001). AVN cells were then left to rest for 1 h before recording.

## Electrophysiological recordings

All electrophysiological recordings were carried out at 36°C, by employing an Axopatch 700B (Axon Instruments Inc.) patch-clamp amplifier connected to a Digidata 1550B interface (Molecular Devices).

The automaticity of AVN cells was recorded by the perforated patch-clamp technique with β-escin (50 µM) in basal extracellular Tyrode solution. Patch-clamp electrodes had a resistance of 4–5 MΩ when filled with an intracellular solution containing (mM): potassium aspartate, 130; NaCl, 10.0; Na-ATP, 2.0; creatine phosphate, 6.6; Mg-GTP, 0.1; CaCl₂, 0.04 (pCa 7.0); and HEPES-KOH, 10.0 (adjusted to pH 7.2 with KOH).

The whole-cell variation of the patch-clamp technique was used to record Ca_v 1.3-mediated ($I_{Cav1.3}$) and Ca_v 1.2-mediated ($I_{Cav1.2}$) $I_{Ca,L}$. The extracellular solution contained (in mM): tetraethylammonium-chloride (TEA-Cl), 135; CaCl₂, 2; MgCl₂, 1; 4-amino-pyridine, 10; D-glucose, 5.5; HEPES, 10; (adjusted to pH = 7.4 with TEAOH). Tetrodotoxin (TTX, 30 µM) was added to the extracellular solution to block $I_{Na}$ current. TTX was supplied by Latoxan (Portes Lès Valences, France). We resuspended lyophilized toxin into a H₂O/20% citrate buffer to obtain a mother solution of 6 mM TTX. The mother solution was then diluted in extracellular solution to about 1/250 (1/1000 citrate buffer, 30 µM TTX; $I_{Na,r}+I_{Na,s}$ block) to 1/5000 (1/10,000 citrate buffer,

100 nM TTX; $I_{Na,s}$ block). Electrodes had a resistance of 3–4 MΩ when filled with a solution containing (in mM): CsOH, 125; TEA-Cl, 20; CaCl₂, 1.2; Mg-ATP, 5; Li₂-GTP, 0.1; EGTA, 5 and 10 HEPES (pH adjusted to 7.2 with CsOH). Nifedipine (3 µM) was added to the extracellular solution with TTX to selectively inhibit $I_{Cav1.3}$ in isolated Ca_v 1.2<sup>DHP−/−</sup> atrioventricular cells. Nifedipine-sensitive $I_{Cav1.3}$ was obtained by subtracting the current recorded under nifedipine ($I_{Cav1.2}$) from the total $I_{CaL}$ recorded in the control TEA-Cl and TTX solution.

## Confocal line-scan imaging of Ca²⁺ in AVN cells

Spontaneous [Ca²⁺]ᵢ transients in isolated AVN cells were measured with confocal microscopy (Zeiss LSM 780). AVN isolated cells were seeded in a Fluoro-dish (FD3510-100, WPI), coated overnight with laminin (1 mg/ml, Sigma-Aldrich, St Louis, MO, USA), for 1 h before recordings. Cells were loaded for 20 min with a Tyrode solution containing the Ca²⁺ indicator Fluo4-AM (5 µM). After washing with Fluo4-AM two times, imaging of AVN cells was performed in a standard Tyrode bath solution pre-warmed at 37°C. Images were acquired by scanning AVN cells with an argon laser in line-scan mode. Fluorescence was excited at 488 nm and emissions were collected at >505 nm. A 63-oil immersion objective was used to record [Ca²⁺]ᵢ in AVN cells. Spontaneous [Ca²⁺]ᵢ transients were analysed using ImageJ software.

## Model development

For developing this new model, we adapted the parent mouse AVN single-cell model published by Mangoni in 2011 (MM2011) (Marger et al., 2011b). Since the latter did not include computation of intracellular calcium handling, we have taken equations governing intracellular ionic fluxes and calcium handling from the mouse SAN AP model published by Kharche et al. (2011).

Figure 1 shows a schematic representation of the different elements included in the present model of AVN automaticity. Four compartments are considered: the sarcoplasmic reticulum (SR), divided into junctional and network spaces (JSR and NSR, 0.12% and 1.16% of the whole cell volume, respectively), the subsarcolemmal calcium subspace, and the cytosol. The calcium subspace is considered as an independent intracellular compartment (1% of the total cell volume) because of the transient local calcium accumulation taking place upon Ca²⁺ entering through membrane Ca²⁺ channels and release from the SR.

The cell geometry was kept as a cylindrical structure where the length and radius of the cell were derived from the total membrane capacitance. This was accomplished by fixing both the membrane capacitance to 22 pF, as
experimentally measured in AVN cells, (Marger et al., 2011a; Mesirca et al., 2014) and the cell radius to 5 µm. The cell's length (65 µm) was calculated using the cylindrical approximation and the specific membrane capacitance set to 1 µF/cm².

Formulations of membrane currents were obtained by fitting the experimental I-V curves from the voltage-clamp activation protocol, and they are reported in the Results section. Currents that were not fitted to experimental data were derived from the parent model (MM2011) (Marger et al., 2011b) with the following adaptations: (1) the sustained inward current $I_{st}$ was reduced by 80%, (2) the maximal sodium-potassium pump current $i_{p,max}$ was increased 3-fold, (3) the calcium background current $I_{bCa}$ was reduced by 50%, (4) the diffusion time constant between the cytosol and sub-sarcolemma space $\tau_{diffCa}$ was made slower: 0.4 ms vs. 0.04 ms, (5) the sodium background conductance $g_{bNa}$ was adjusted to 0.0053 nS/pF, (6) the slow potassium conductance $g_{Ks}$ was modified to 1.88 nS/pF, (7) the formulation of $P_{up}$ was updated to Kharche's formulation (kept to the value set in Kharche's model).

All intracellular ionic concentrations ($Na_i^+, Ca_i^{2+}, K_i^+$) have been described with dynamic equations and calcium buffering (troponin $Ca^{2+}$ site (TC) and troponin $Mg^{2+}$ site (TMC) in the cytosol, calmodulin located in the subspace ($CM_s$) and cytosol ($CM_i$) and calsequestrin (CQ) in

JSR) was included and taken from the model by Kharche (Kharche et al., 2011). Because of the existence of subspace and SR, four fluxes controlling calcium handling and balance have been added: the Sarco-Endoplasmic Reticulum Calcium ATPase pump (SERCA, $J_{up}$), the flux between NSR and JSR ($J_{tr}$), the release from JSR to the sub-sarcolemmal subspace ($J_{rel}$), and the $Ca^{2+}$ diffusion from the subspace to the whole cytosol ($J_{diffCa}$).

## Apparent input resistance estimation

We have estimated the apparent input resistance in basal conditions by stepping the voltage to ±1,±2 mV around the MDP and fitting with a linear (ohmic) relationship the total membrane current vs. voltage data points. The same procedure was repeated in the three conditions for which the cell becomes quiescent ($I_{Cav1.3}/I_{Na,r}$ and $I_{Kr}$ block).

## Population of models

The developed model was used to construct a population of *in silico* models based on the population of models methodology (Gemmell et al., 2014; Gong & Sobie, 2018; Muszkiewicz et al., 2016), as previously done by our group (Paci et al., 2018, 2020). We sampled a total of 18 parameters in the [80–120]% range, namely all the main ionic conductances governing the AP ($g_{Na,s}$, $g_{Na,r}$, $g_{Cav1.2}$,

<figure><img src="image_3.png" /><figcaption>**Figure 1. Schematic diagram of the mouse AVN single-cell model**</figcaption></figure>

$Ca^{2+}$ currents are depicted in red, $Na^+$ currents in purple, and $K^+$ currents in green. The model features three compartments: the cytoplasm (CYTO), the subspace (Sub-Sarcolemma), and the sarcoplasmic reticulum (SR), divided into network (NSR) and junctional (JSR). Buffers are present in the different compartments: calmodulin (CMi), troponin $Ca^{2+}$ site (TC) and troponin $Mg^{2+}$ site (TMC) in the cytosol; calmodulin in the subspace (CMs); calsequestrin (CQ) in the JSR. See Table A1 in the Appendix for currents legend. [Colour figure can be viewed at wileyonlinelibrary.com]
$g_{\text{Cav1.3}}$, $g_{\text{Cav3.1}}$, $g_{\text{K1}}$, $g_{\text{Kr}}$, $g_{\text{Ks}}$, $g_{\text{to}}$, $g_{\text{Kur}}$, $g_{\text{f}}$, $g_{\text{bNa}}$, $g_{\text{bK}}$, $g_{\text{bCa}}$, $g_{\text{st}}$, $g_{\text{sK}}$) and the maximal currents $I_{\text{NaK,max}}$ and $I_{\text{NaCa,max}}$. An initial population of 3000 AVN cells was generated and then calibrated based on available experimental data (MDP, $dV/dt_{\text{max}}$, APD$_{90}$, and rate in bpm from Choisy et al. (2015); Marger et al. (2011a); Mesirca et al. (2014, 2021)). Only if the cell biomarkers fell into two standard deviations around the average experimental values, was the model selected. Biomarkers were computed at steady state (300 s) as the average over the last 2 s. The population was calibrated in a non-paced configuration since all of the experimental data were acquired in non-stimulated cells showing spontaneous activity. Block of ionic currents was tested on the calibrated cells as was done for the basal model to evaluate the percentage of cells composing the population showing cessation of spontaneous activity.

## Hardware and software

Model differential equations were implemented in Matlab (Mathworks Inc., Natick, MA, USA) and solved with a variable order solver (ode15s), based on numerical differentiation formulas (Shampine & Reichelt, 1997), with a maximum time step of 1 ms. Automaticity was calculated for 500 s to ensure reaching the steady-state. Figures show the last seconds of the simulation. The AP biomarkers are: beating rate (beats per minute, bpm), maximum diastolic potential (MDP), AP amplitude (APA), and AP duration (APD$_{90}$ measured once the membrane voltage reached 90% of the resting voltage). AP biomarkers were calculated as averaged values from the last 2 s of the steady-state simulation at 500 s.

## Results

### Membrane currents

In this section, we specify sarcolemmal currents flowing through ion channels, ionic pumps, and exchangers (see Fig. 1), whose formulations have been tuned during model development. In addition, we present results of fitting experimental data derived from the literature and of original data we collected in the framework of model development.

The starting point for the Hodgkin-Huxley type formulations of the currents is the preliminary mouse AVN single-cell MM2011 model (Marger et al., 2011b). To reformulate ionic currents in the new model and obtain an improved description of AVN pacemaker activity, we have used already published experimental data in the form of current-to-voltage (I-V) curves from mouse AVN single cells from (Choisy et al., 2015; Marger et al., 2011a) and new data collected for this study. The full set of equations of the new model is reported in Appendix 1.

### L-type Ca²⁺ currents ($I_{\text{Cav1.2}}$ and $I_{\text{Cav1.3}}$)

The MM2011 model did not independently compute $I_{\text{Cav1.2}}$ and $I_{\text{Cav1.3}}$. We thus employed AVN cells from Ca$_v$1.2<sup>DHP-/-</sup> mice to experimentally separate the I-V curve of $I_{\text{Cav1.2}}$ from that of $I_{\text{Cav1.3}}$ (see Methods). Results are shown in panels Aa and Ab in Fig. 2. The light grey squares show the new experimental data, while filled black dots show the fit with the new model. The new model satisfactorily fitted experimental data. Panels B show the model activation ($dL_∞$, continuous lines) and inactivation ($fL_∞$, dashed lines) gates for $I_{\text{Cav1.2}}$ (Ba) and $I_{\text{Cav1.3}}$ (Bb), respectively. The fit with the MM2011 current formulation is also shown for $I_{\text{Cav1.3}}$: in the new model the maximal conductance is reduced (0.1365 vs. 0.1818 nS/pF), the current reversal potential is computed with the Nerst formula and the activation slope is set to 4.3 mV instead of 6.3 mV (Fig. 2Bb). For $I_{\text{Cav1.2}}$, the fit with the parent model is not shown because this current was not computed in the MM2011 model. For both $I_{\text{Cav1.2}}$ and $I_{\text{Cav1.3}}$ the voltage dependence of inactivation time constants was shifted of 40mV to hyperpolarized potentials respect the previous formulation.

### T-type Ca²⁺ current ($I_{\text{Cav3.1}}$)

To fit the T-type Ca²⁺ current encoded by Ca$_v$3.1 ($I_{\text{Cav3.1}}$), we have used the experimental data from (Marger et al., 2011b). Experimental and computed I-V curve are shown in Fig. 3A together with the activation ($dT_∞$) and inactivation ($fT_∞$) gates for $I_{\text{Cav3.1}}$ (Fig. 3B). Figure 3 also reports the fit with the original model (dark grey traces): in the new model the maximal conductance is reduced (0.2658 vs. 0.3091 nS/pF), the reversal potential is computed with the Nerst formula, the activation curve (Fig. 3B) is shifted to the right (half maximal activation -40 vs. -45 mV) and its slope is incremented (6 vs. 4.64 mV).

### Voltage gated sodium currents – TTX-resistant ($I_{\text{Na,r}}$) and TTX-sensitive components ($I_{\text{Na,s}}$)

Mouse SAN and AVN cells express both TTX-sensitive ($I_{\text{Na,s}}$) and TTX-resistant ($I_{\text{Na,r}}$) Na⁺ currents (Lei et al., 2005; Marger et al., 2011a). In Fig. 4 I-V curves presenting the fit of experimental data from (Marger et al., 2011a) (open light grey dots) with the new and original MM2011 model (black vs. dark grey filled markers) are reported, with TTX- sensitive component in panel Aa and the resistant component in panel Ab. Panels B show the model activation ($m_∞$) and inactivation ($h1_∞$ and $h2_∞$) gates for $I_{\text{Na,s}}$ (Ba) and $I_{\text{Na,r}}$ (Bb). Moreover, for both currents, the inactivation time constant formulations have been substituted with the one from Kharche (2011) with a small
change in τ<sub>h2</sub> in order to speed up the inactivation (see equation in the Appendix). Regarding I<sub>Na,s</sub>, the maximal conductance is reduced with respect to the MM2011 model, the activation curve (Fig. 4Ba) is shifted to the left (half-maximal activation –38 vs. –31 mV) and its slope increased (6 vs. 4.4 mV). For I<sub>Na,r</sub>, the maximal conductance is reduced (4.25e<sup>-5</sup> vs. 5e<sup>-4</sup> nS/pF), the activation curve (Fig. 4Bb) is slightly shifted to the right (half-maximal activation –41 vs. –43 mV), and the slope reduced (4 vs. 7 mV). Finally, the inactivation slope of h<sub>s1</sub> is increased (6 vs. 4 mV). Since experimental data for sodium currents were obtained at room temperature (22°C) we used a Q<sub>10</sub> factor equal to 2.79 for all time constants (multiplicative factor of 0.2378) and Q<sub>10</sub> = 1.5 for conductances from Kernick et al. (Kernick et al., 2019).

Hyperpolarization-activated funny current (I_f)

I_f was present in all spontaneously beating AVN cells recorded in Marger et al. (2011a). The corresponding

I_f experimental I–V curve is shown in Fig. 5A with open light grey dots, while the new and original model fits are shown in black and light grey filled symbols, respectively. Figure 5B reports the activation (y∞) gate for I_f. The I_f formulation has been updated with the one used by Kharche (Kharche et al., 2011), enabling a much better fit of experimental data. Nevertheless, some modifications have been made to the parameters. In particular: the maximal conductance has been reduced by 20%, the half-activation voltage has been shifted to –97 mV (from –106.8 mV), and its slope factor has been increased from 16.3 to 20. Furthermore, although Kharche's formulation presented different contributions from sodium and potassium currents to I_f, we used the same proportion for both as there is no evidence for different contributions of the two components in the AVN. The time constant (τ<sub>P</sub>) was also modified from Kharche's formulation to give a peak at around –50 mV and a slower activation in the diastolic phase.

<figure><img src="image_4.png" /><figcaption>Figure 2. L-type Ca²⁺ currents I-V curve fitting and activation/inactivation gates</figcaption></figure>

Voltage-dependent, L-type Ca²⁺ currents (I<sub>Cav1.2</sub>, panel Aa) and (I<sub>Cav1.3</sub>, panel Ab) I–V curve fits with old (dark grey) and new models (black) against new experimental data (light grey, *n* = 14, mean ± SD). The inset in panel Aa represents the voltage clamp protocol adopted for both currents. Extracellular ion concentrations (in mM): Ca²⁺, 2. Intracellular ion concentrations (in mM): Ca²⁺, 90e⁻³. Panels Ba and Bb report the model activation (dL<sub>∞</sub>) and inactivation (fL<sub>∞</sub>) gates for I<sub>Cav1.2</sub> and I<sub>Cav1.3</sub>, respectively, with the old (grey) and new formulations (black; n.c., normalized conductance). For I<sub>Cav1.2</sub>, the old fit is not shown because this current was not present in the original model.
## Rapid delayed rectifier K⁺ current (I<sub>Kr</sub>)

$I_{Kr}$ was found in all AVNCs tested in the previous work by Marger et al. (Marger et al., 2011a). We have fitted the new formulation of $I_{Kr}$ using experimental data by Marger et al. Comparison of I-V curves of the former and current models (dark grey vs. black filled markers), and experimental data is shown in Fig. 6A. Simulations of activation (Paf<sub>∞</sub>) and inactivation (Pi<sub>∞</sub>) gates are shown in Fig. 6B. In the new model, $I_{Kr}$ formulation maximal conductance has been reduced by 70% (0.068 vs. 0.227 nS/pF), only the fast activating component has been included, and the activation time constant has been modified (see equation in the Appendix), in order to speed up deactivation to obtain a less negative MDP closer to the experimental values of spontaneous action potentials.

## Ultra-rapid delayed rectifier K⁺ current (I<sub>Kur</sub>) and transient outward K⁺ current (I<sub>to</sub>)

The 4-AP-sensitive current, sum of $I_{to}$ and $I_{Kur}$, is shown in Fig. 7. Panel A reports the I-V curve and Panel B the activation ($r_{∞}$, continuous line) and inactivation ($q_{∞}$, dashed line) gates. The maximal conductance of both currents has been increased: for $I_{to}$ the new value is 0.04 nS/pF (vs. 0.022), for $I_{Kur}$ the new value is 0.015 nS/pF (vs. 0.003).

## Sodium/calcium exchange current (I<sub>NaCa</sub>)

To model the sodium/calcium exchanger, we fitted original data published by Choisy et al. (Choisy et al., 2015) with equations taken from the human pluripotent stem cell model from Paci et al. (Paci et al., 2012). Results of the new fit (filled black dots) are shown in Fig. 8

with the MM2011 formulation (filled dark grey dots) and experimental data (open light grey dots) for comparison.

## Mouse AVN model behaviour in basal conditions

Figure 9 shows the behaviour of the present model in terms of simulation of spontaneous AP, intracellular Ca²⁺ transients (Fig. 9B, C and G) and main ionic currents ($I_{Cav1.2}$, $I_{Cav1.3}$, $I_{Cav3.1}$, $I_{Na,r}$, $I_{Na,s}$, $I_{NaCa}$) (Fig. 9D, E) in comparison with different experimental examples (Fig. 9A, F). Simulated AP and Ca²⁺ dynamics hallmarks are compared to experimental data in the table embedded in Fig. 9.

Figure 10 shows the time course of individual currents that play a relevant role during diastolic depolarization (DD). In the early diastolic phase, $I_{Kr}$ slowly deactivates from 2.7 pA/pF to about 0.1 pA/pF within the first 40 ms. This slow deactivating behaviour allows the MDP to be set at the end of repolarization phase and diastolic inward currents then slowly take over the net diastolic current, allowing diastolic depolarization. In addition, the model predicts residual $I_{Kr}$ in the second half of the diastolic depolarization to the action potential upstroke phase, providing a background outward current component slowing basal diastolic depolarization. The model also predicts tiny $I_{Kur}$ background component. In the first part of diastolic depolarization $I_f$ and $I_{Cav3.1}$ supply the predominant fraction of inward current component to promote membrane depolarization. In the second part of diastolic depolarization, the TTX-resistant and TTX-sensitive sodium current and $I_{Cav1.3}$ activate and become more relevant. Notably, the model predicts $I_f$ to constantly increase in diastolic depolarization.

The model was also tested under paced conditions at different CLs: results are shown in Fig. 11.

<figure><img src="image_5.png" /><figcaption>Figure 3. T-type Ca²⁺ current I-V curve fitting and activation/inactivation gates</figcaption></figure>

A, T-type Ca²⁺ current I-V curve fitting with old (dark grey) and new model (black) against experimental data (light grey, n = 4) from Marger et al. (2011a). The inset represents the voltage clamp protocol adopted. Extracellular ion concentrations (in mM): Ca²⁺, 2. Intracellular ion concentrations (in mM): Ca²⁺, 200e⁻³. B, model activation (dT<sub>∞</sub>) and inactivation (fT<sub>∞</sub>) gates for I<sub>Cav3.1</sub> with the old (dark grey) and new formulation (black).
Population of models

Of the 3000 cells comprising the initial population, 370 passed the calibration. The remaining ones were discarded due to the absence of spontaneous oscillation or because

their AP markers were outside the experimental ranges (see Methods). Figure 12 shows the APs of the basal model and of the cells comprising the population, while Table 1 reports the comparison between their AP features. The majority of the population models shows a slightly lower

<figure><img src="image_6.png" /><figcaption>Figure 4. Na⁺ currents I-V curve fitting and activation/inactivation gates</figcaption></figure>

Na⁺ currents (I<sub>Na,s</sub> TTX-sensitive, panel Aa) and I<sub>Na,r</sub> TTX-resistant, panel Ab) I-V curves fits with old (dark grey) and new model (black) against experimental (light grey, *n* = 8) data from Marger et al. (2011a). The inset in panel Aa represents the voltage clamp protocol adopted for both currents. Extracellular ion concentrations (in mM): Na⁺, 45. Intracellular ion concentrations (in mM): Na⁺, 19. Ba and Bb report the model activation (m<sub>∞</sub>) and inactivation (h1<sub>∞</sub> and h2<sub>∞</sub>) gates for I<sub>Na,s</sub> and I<sub>Na,r</sub>, respectively, with the old (dark grey) and new formulation (black).

<figure><img src="image_7.png" /><figcaption>Figure 5. Funny current I-V curve fitting and activation gate</figcaption></figure>

A, *I*<sub>f</sub> I–V curve fitting with MM2011 (dark grey) and new model (black) against experimental data (light grey, *n* = 5, mean ± SEM) from Marger et al. (2011a). The inset represents the voltage clamp protocol adopted. Extracellular ion concentrations (in mM): Na<sup>+</sup>, 145; K<sup>+</sup>, 5.4. Intracellular ion concentrations (in mM): Na<sup>+</sup>, 8; K<sup>+</sup>, 140. *B*, model activation (y<sub>∞</sub>) gate for *I*<sub>f</sub> with the old (grey dashed) and new formulation (black).
rate than the basal model, but all the other features fall inside the population ranges.

## Model validation through the simulation of ion channel blocks

We have employed our new model to simulate the effects of blockade of ionic currents predicted to be involved in AVN automaticity. In Table 2 we summarized the simulation results of ionic current blockade in the basal model and in the model's population and compared them with experimental data. The results show agreement between simulations and experimental data and highlight the main currents determining the AP (Table 2).

$I_f$ is considered as one of the key pacemaker mechanisms in the AVN (Marger et al., 2011b; Mesirca

et al., 2021; Saeed et al., 2018; Temple et al., 2016). We thus simulated partial (32%) blockade of $I_f$ and compared it to experimental observations reported by Marger et al. (2011a) (Fig. 13C and Table 2). In the basal model, partial blockade of $I_f$ slows down the rate of pacemaker activity (155 vs. 175 bpm), to a similar percentage as experimentally reported (11% vs. 16%). Simulating total (100%) block of $I_f$ resulted in very slow spontaneous activity (89 bpm) but without arrest of pacemaking (Fig. 13C, dashed line and Table 2).

We then simulated 100% block of the two $I_{Na}$ components, $I_{Na,r}$ and $I_{Na,s}$. For $I_{Na,r}$, the model predicted arrest of pacemaker activity (Fig. 13A), consistent with the results obtained experimentally with 20 µM TTX by Marger et al. (Marger et al., 2011a). However, for $I_{Na,s}$, an 11% rate reduction is obtained (156 bpm) (Fig. 13B),

<figure><img src="image_8.png" /><figcaption>**Figure 6.** $I_{Kr}$ I-V curve fitting and activation/inactivation gates</figcaption></figure>

$A$, $I_{Kr}$ current fitting with MM2011 (dark grey) and new model (black) against experimental data (light grey, $n = 20$) from Marger et al. (2011a). The inset shows the voltage clamp protocol adopted. Extracellular ion concentrations (in mM): $K^+$, 5.4; Intracellular ion concentrations (in mM): $K^+$, 140. $B$, model activation ($Paf_\infty$) and inactivation ($Pif_\infty$) gates for $I_{Kr}$.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** 4-AP sensitive current I-V curve fitting and activation/incativation gates</figcaption></figure>

$A$, 4-AP-sensitive current (sum of $I_{to}$ and $I_{Kur}$) I-V curve fitting with MM2011 (dark grey) and new model (black) against experimental data (light grey, $n = 5$, mean ± SEM) from Marger et al. (2011a). The inset represents the voltage clamp protocol adopted. Extracellular ion concentrations (in mM): $K^+$, 5.4. Intracellular ion concentrations (in mM): $K^+$, 140. $B$, model activation ($r_\infty$) and inactivation ($q_\infty$) gates for $I_{4AP}$.
<table><caption>**Table 1.** Comparison of average and (min:max) interval AP features values of the calibrated population (n = 370) and the basal model.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Population</th><th rowspan="2">Present model</th></tr><tr><th>Mean</th><th>Min</th><th>Max</th></tr></thead><tbody><tr><td>Rate (bpm)</td><td>148</td><td>121</td><td>191</td><td>175</td></tr><tr><td>CL (ms)</td><td></td><td></td><td></td><td>342</td></tr><tr><td>MDP (mV)</td><td>-65</td><td>-67</td><td>-63</td><td>-65</td></tr><tr><td>APA (mV)</td><td>105</td><td>99</td><td>110</td><td>101</td></tr><tr><td>E<sub>th</sub> (mV)</td><td>-28</td><td>-31</td><td>-26</td><td>-27</td></tr><tr><td>dV/dt<sub>max</sub> (mV/ms)</td><td>8.6</td><td>7.1</td><td>10.3</td><td>7.3</td></tr><tr><td>APD<sub>90</sub> (ms)</td><td>90</td><td>84</td><td>97</td><td>89</td></tr><tr><td>CaTD<sub>90</sub> (ms)</td><td>217</td><td>162</td><td>273</td><td>177</td></tr></tbody></table>

milder than available experimental data, which show remarkable slowing of pacemaker activity by 53 ± 13.6% (Nikmaram et al., 2008).

We then simulated total blockade of $I_{Cav1.3}$. Under this condition, our model predicts the arrest of spontaneous AP discharge to resting membrane potential of -33.7 mV (Fig. 14A). This result is consistent with experimental evidence by Marger et al. (Marger et al., 2011b), which showed that $Ca_v1.3^{-/+}$ AVN cells had no spontaneous

activity and stabilized to a depolarized resting potential close to -35 mV. Simulating blockade of $I_{Cav1.2}$ (Fig. 14B) predicted an increase in spontaneous AP rate of +40% (245 vs. 175 bpm). Since no experimental data are available showing selective $I_{Cav1.2}$ block, we analysed the reasons for such an apparently counterintuitive numerical prediction. During the diastolic phase, most currents are not modified (Fig. A1), but we observe that upon $I_{Cav1.2}$ block $I_{Ca3.1}$ increases due to a slower inactivation (see the time constant in Fig. A2) which is caused by a reduced AP amplitude (lower OS). As a consequence, the $I_{Ca3.1}$ increase produces a steeper diastolic phase and thus a shorter CL (Fig. A3). Concerning the T-type calcium current, we also simulated $I_{Cav3.1}$ block (Fig. 14C). The model shows a 9% rate reduction (160 vs. 175 bpm), slightly lower than in our preliminary experimental data (WT: 191 ± 8, $n=6$; $Ca_v3.1^{-/+}$: 156 ± 16, $n=9$. P = 0.08. 18% reduction).

We have simulated 4-AP-sensitive current block by blocking $I_{Kur}$ at 50% corresponding to 50 µM 4-AP used in the experiment by Nikmaram et al. (Nikmaram et al., 2008) (Fig. 15B). In this condition, we observe a negligible effect on the rate (179 vs. 175 bpm), in contrast with the slight rate increase reported experimentally (106 vs. 125 bpm) (Nikmaram et al., 2008). The last tested block is a 100% reduction in $I_{Kr}$ (Fig. 15A), for which the cell stops beating. The experimental data show a CL prolongation of 28 ± 3.4% (Nikmaram et al., 2008), and a similar CL prolongation (370 ms vs. 286 ms, +29%) was obtained by simulating a current block by 70%.

We have also estimated the apparent input resistance of the model: in basal condition its value is 550 MΩ, in line with experimental values found in SAN cells (Verheijck et al., 1998). To our knowledge, there are no data on AVN cells. Under sodium current block, the resistance was quite similar (584 MΩ); when $I_{Cav1.3}$ is blocked, it is slightly lower (327 MΩ), which seems to confirm that upon this block it will become more difficult to restore the cell's excitability. Instead upon $I_{Kr}$ block the resistance increases (2.25 GΩ).

As for the basal model, 100% $I_{Cav1.3}$ block leads to cessation of automaticity in every cell of the population (Table 2). With 50% $I_f$, a more pronounced drop in rate (33%) was obtained on average in the population, while upon 100% block the rate further dropped to 54 ± 23 bpm. Still, each cell composing the population showed spontaneous pacing. Complete block of $I_{Na,r}$ – while stopping the basal model – sees half of the population maintaining periodic oscillations. However, simultaneous 100% block in $I_{Na,r}$ and $I_{Na,s}$ lead to loss of pacemaking in almost every model (98%) of the population, in line with single cell experimental data (Marger et al., 2011a). Similarly, complete block of $I_{Kr}$ stops the beating in the basal model but only in half of the models of the population, which is more consistent with

<figure><img src="image_10.png" /><figcaption>**Figure 8.** $Na^+/Ca^{2+}$ exchange current I-V curve fitting</figcaption></figure>

$Na^+/Ca^{2+}$ exchange current I-V curve fitting with MM2011 (dark grey) and new model (black) against experimental data (light grey, $n=6$, mean ± SEM) from Choisy et al. (2015). The inset represents the voltage ramp adopted to compute the current. Extracellular ion concentrations (in mM): $Na^+$, 140; $Ca^{2+}$, 2. Intracellular ion concentrations (in mM): $Na^+$, 10; $Ca^{2+}$, $5e^{-5}$.
the experimental data reporting slowing of the beating rate (Nikmaram et al., 2008).

Discussion

In the pursuit of a comprehensive understanding of AVN cellular electrophysiology, the present study aimed to develop a novel computational model of mouse single cell AP and automaticity, starting from the mouse SAN model by Mangoni et al. (2006) and from a former formulation of an AVN model, the MM2011 model (Marger et al., 2011b). We introduced several new structural features into the model in comparison with the parent MM2011, namely: subcellular compartimentalization (Fig. 1), computation of calcium handling, dynamic computation of ionic intra-cellular concentrations and of calcium buffering. We have paid particular attention to the critical reviewing of the formulations of membrane currents and exchangers in light of all the available data from native AVN mouse cells, both from the literature and by including new *ad hoc*

experimental recordings of ionic currents. A stringent fit of I-V curves was obtained for all the considered currents (Figs 2–8). The model was then validated by comparing it with several experimental results in terms of block of individual ionic currents and by investigating their impact on the AP.

In control conditions, the model reproduces the main experimental AP features. The simulated AP automaticity is characterized by a CL of 286 ms, close to experimentally observed values (Choisy et al., 2015; Marger et al., 2011a; Mesirca et al., 2014, 2021). Critical parameters like MDP, APA and APD₉₀ lay within the range of the experimental data variability (Fig. 9 and Table within). Regarding the calcium transient, the model can reproduce its experimental time course (Fig. 9, panels F and G) as well as quantitative information such as CaTD at 30%, 50%, and 90% of repolarization (see Table in Fig. 9). Although a few specific AP features do not exactly match the experimental ones (upstroke steepness, dV/dt<sub>max</sub> and AP triangular shape, APD₅₀), the model notably allowed reproduction of almost all the experimentally reported

<figure><img src="image_11.png" /><figcaption>Figure 9. Example of experimental AP and cytosolic Ca²⁺ transient and the corresponding time courses simulated together with the ionic currents</figcaption></figure>

A–F, experimental AP and cytosolic Ca²⁺ transients. B–E and G, simulated action potentials of single-cell AVN model and Ca²⁺ transient in the cytosol. C, Ca²⁺ transients in the cellular submembrane space. D and E, calcium and sodium currents (I<sub>Cav1.2</sub>, I<sub>Cav1.3</sub>, I<sub>Cav3.1</sub>, I<sub>Na,r</sub>, I<sub>Na,s</sub>, I<sub>NaCa</sub>). Comparative tables between experimental and simulated values for AP and Ca²⁺ transient biomarkers are reported on the bottom right. AP biomarkers are reported in panel B: OS: overshoot; TOP: take-off potential; MDP: maximum diastolic potential; APD: action potential duration (at 90% of repolarization). Ca²⁺ transient biomarkers are reported in panel G: dCaᵢ/dt<sub>max</sub>: calcium transient maximum derivative; CaTD: calcium transient duration (at 30%, 50%, and 90% of repolarization). [Colour figure can be viewed at wileyonlinelibrary.com]
<table><caption>Table 2. Summary of the current block effects on the AP rate (in control condition 210 bpm) both in basal model and in the population</caption><thead><tr><th>Current block</th><th>Basal model rate (bpm) and relative change (%)</th><th>Non-beating cells in the population</th><th>Mean ± std rate (bpm) and average relative change (%) in the population</th><th>Experimental rate change (%)</th></tr></thead><tbody><tr><td>100% I<sub>Na,s</sub></td><td>156 (-11%)</td><td>0</td><td>135 ± 11 (-9%)</td><td>-33% (Nikmaram et al., 2008)</td></tr><tr><td>100% I<sub>Na,r</sub></td><td>No beat</td><td>187 (51%)</td><td>167 ± 25 (+13%)</td><td>—</td></tr><tr><td>100% I<sub>Na,r</sub> + I<sub>Na,s</sub></td><td>No beat</td><td>363 (98%)</td><td>125 ± 69 (-2%)</td><td>No beat (Marger et al., 2011a)</td></tr><tr><td>100% I<sub>Cav1.3</sub></td><td>No beat</td><td>370 (100%)</td><td>—</td><td>No beat (Marger et al., 2011a)</td></tr><tr><td>100% I<sub>Cav1.2</sub></td><td>245 (+40%)</td><td>0</td><td>199 ± 20 (+34%)</td><td>—</td></tr><tr><td>100% I<sub>Cav3.1</sub></td><td>160 (-9%)</td><td>0</td><td>139 ± 11 (-7%)</td><td>-18% (new exp data)</td></tr><tr><td>50% I<sub>kur</sub></td><td>179 (+2%)</td><td>0</td><td>152 ± 13 (+2%)</td><td>-15% (Nikmaram et al., 2008)</td></tr><tr><td>100% I<sub>Kr</sub></td><td>No beat</td><td>194 (52%)</td><td>133 ± 19 (-8%)</td><td>-23% (Nikmaram et al., 2008)</td></tr><tr><td>32% I<sub>f</sub></td><td>155 (-11%)</td><td>0</td><td>99 ± 15 (-33%)</td><td>-16% (Marger et al., 2011a)</td></tr><tr><td>100% I<sub>f</sub></td><td>89 (-49%)</td><td>0</td><td>54 ± 23 (-64%)</td><td>—</td></tr></tbody></table>

effects of current blockades, as evident in Figs 13–15. In addition, even though there are no experimental data available for the mouse AVN stimulated condition, we tested our model by pacing it at physiological rates and the rate dependence (shown in Fig. 11) expresses its well known behaviour, similar to what has been reported by Workman et al. for rabbit AVN cells (Workman et al., 2000). The results also show CL values for which APD alternans and 1:2 capture are obtained, demonstrating the model’s potential for investigating such phenomena.

Key among our refinements of membrane current formulations were the L-type Ca²⁺ currents. Two different L-type Ca²⁺ channels, namely Ca_v1.2 and Ca_v1.3, are differentially expressed in the atria and cardiac conduction system (Marionneau et al., 2005) and play a central role in the automaticity of SAN and AVN cells. The role of Ca_v1.2 in generating the cardiac AP upstroke and that of Ca_v1.3 in generating the SAN diastolic depolarization are well established (Mesirca et al., 2015), their function in generating AVN automaticity still needs to be fully elucidated (Marger et al., 2011b) and it is not fully quantified. Our new experimental data confirm the expression of Ca_v1.2 in AVN cells, despite I<sub>Cav1.2</sub> not being included in the parent model MM2011, as well as providing separate computations of I<sub>Cav1.2</sub> and I<sub>Cav1.3</sub> in generating total I<sub>Ca,L</sub> in the AVN (Fig. 2), with Ca_v1.3-mediated I<sub>Ca,L</sub> activating at slightly more negative voltages. Our new experimental data using Ca_v1.2<sup>DHP-/-</sup> AVN cells are in line with the previous observation by Marger et al. (2011b) that Ca_v1.3 channels constitute

the predominant L-type channel isoform in these cells. Nevertheless, a significant contribution of I<sub>Cav1.2</sub> was found, with a peak current density around 2.5 pA/pF, which is about half of I<sub>Cav1.3</sub> density (5 pA/pF) (see Fig. 2). This observation in AVN is also consistent with our previous estimation that in mouse SAN cells I<sub>Cav1.3</sub> accounts for about 70% of the total I<sub>Ca,L</sub> (Mangoni et al., 2003), which indicates that Ca_v1.3 expression is up-regulated in the primary (SAN) and secondary (AVN) pacemaker centres. Our successful fitting of new experimental data is shown in Fig. 2. In this study, we employed Ca_v1.2<sup>DHP-/-</sup> mice to quantitatively separate I<sub>Cav1.2</sub> from I<sub>Cav1.3</sub>. These mice present the advantage of limiting potential interference of cellular remodelling that might have been present in Ca_v1.3<sup>-/-</sup> AVN myocytes in the work by Marger et al. (Marger et al., 2011b). In Ca_v1.3<sup>-/-</sup> SAN myocytes, we did not observe remodelling of calcium-handling components. Indeed, SAN myocytes express similar densities of I<sub>f</sub> (Mangoni et al., 2003) and similar NCX/SERCA activity (Torrente et al., 2016) in comparison to wild-type counterparts. In addition, we did not observe overt remodelling in Ca_v1.2<sup>DHP-/-</sup> SAN myocytes in relation to diastolic calcium release and spontaneous firing (Louradour et al., 2022). In conclusion, even though we cannot exclude a differential degree of remodelling in SAN versus AVN myocytes, we think it unlikely that remodelling has affected our experimental data.

The importance of Ca_v1.3 in AVN physiology channels has been assessed experimentally by Marger et al. (Marger
et al., 2011a, 2011b) where perfusion of isradipine blocked AP discharge and $Ca_v1.3^{-/-}$ AVNCs did not display spontaneous beating: the model confirms this result by showing the cessation of firing upon 100% block of $I_{Cav1.3}$ (Fig. 14A). Our new model also predicts that this current is also relevant in the last part of the DD

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Close up view of computed diastolic membrane voltage (A) and main ionic currents (B) activated during diastolic depolarization</figcaption></figure>

The inset at the top left reports a zoom on the first 50 ms showing the amplitude of $I_{Kr}$, while the one on the top right reports the sum of the outward currents ($I_{out}$), the sum of the inward currents ($I_{in}$) and the total net current ($I_{tot}$) during diastole. [Colour figure can be viewed at wileyonlinelibrary.com]

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Rate dependence of the baseline model</figcaption></figure>

Light and dark grey areas mark CL values for which APD alternans and 1:2 capture are obtained, respectively.

phase, as highlighted in Fig. 10. With 100% block of $Ca_v1.2$ (Fig. 14B), the AP rate increases: since there is no experimental data to validate this result, the model was used to investigate the mechanisms behind such behaviour. $I_{Cav1.2}$ block reduces the AP amplitude (lower OS), which in turn causes $I_{Cav3.1}$ to deactivate more slowly (Figs A2, A3). As a result, $I_{Cav3.1}$ increases and produces a steeper diastolic phase and a shorter CL. Given the unexpected behaviour of a reduced CL following $I_{Cav1.2}$ block, experimental validation of this prediction is required to verify its validity.

Likewise, reformulation of T-type $Ca^{2+}$ current ($I_{Cav3.1}$) was also performed to achieve a more accurate representation of its behaviour and it is showcased in Fig. 3C. $Ca_v3.1$ channels are one of the relevant contributors to the diastolic phase (Fig. 10) and its block, simulated with the model, indicates a slowing of the rate in qualitative accordance with new data and in the same direction as previous experimental evidence (Marger et al., 2011a).

Sodium currents, both TTX-sensitive and TTX-resistant, underwent meticulous recalibration as detailed in Fig. 4. $I_{Na}$ current is considered a predominant mechanism underlying the AP upstroke phase (Marger et al., 2011a). AVNCs expressed both $I_{Na,r}$ and $I_{Na,s}$ isoforms, which were consequently included in the model: their role in the AP firing is predicted in Fig. 13. As expected, the model showed no AP discharge when $I_{Na,r}$ was blocked (Fig. 13A). The experimental evidence on $I_{Na,s}$ block did not show the complete arrest of pacemaker activity, but a considerable prolongation of the CL (Nikmaram et al., 2008), which is qualitatively reproduced by the model (11% longer CL). Given the relevant functional role of sodium current on AV node

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Action potential traces for the population</figcaption></figure>

Steady-state action potentials of the basal model (black) and of the calibrated cell population (grey, *n* = 370).
myocytes we reconstructed the reactivation time courses
of both these currents along with the functional refractory
period of the *in silico* myocyte. We have tested the P1/P2
protocol and the result is shown in Fig. A4. It can be
seen that after 2 ms, both I<sub>Na</sub> components have recovered
from inactivation by more than 90%, and after 3 ms
by 100%. Due to this fast recovery, the impact on
the refractory period is minimal as shown in Fig. A5
where S1S2 protocol was performed. In this condition,

a full AP is obtained for diastolic intervals (DI) longer
than 7 ms.

The role of I<sub>f</sub> was similarly refined, incorporating
insights from Kharche's formulation to enhance
agreement with experimental data, as exemplified in
Fig. 5. HCN channels are important for basal excitability
not only in SANCs (Ricci et al., 2023) but also in AVNCs
of the mouse, as experimentally reported (Choisy et al.,
2015; Marger et al., 2011b). The contribution of the

<figure><img src="image_15.png" /><figcaption>Figure 13. Simulated AP in control condition (grey line) and under blockade of different currents as indicated (black line)</figcaption></figure>

A, 100% I<sub>Na,r</sub>. B, 100% I<sub>Na,s</sub> block. C, 32% (continuous lines) and 100% (dashed line) I<sub>f</sub> blocks.
funny current to the automaticity of AVNCs was assessed by Marger (Marger et al., 2011b), where the specific blocker ZD-7228 reduced the DD slope and slowed the pacemaker activity by about 16%. We have investigated the effect of the experimental block amount (32%) with the model, obtaining a slight rate reduction of 11%. We also tested the effect of a 100% $I_f$ block in the model (Fig. 13C). The results show dramatic CL prolongation (~49%), supporting a central role in the pacemaking mechanisms of AVNCs, as could be expected by looking

at its amplitude relative to other diastolic curves (Fig. 10). This result is also confirmed by Mesirca and co-workers (Mesirca et al., 2014) using conditional cardiac-specific transgenic mouse model expressing non-conducting hHCN4, which acted as dominant negative subunits to suppress $I_f$ in order to study the importance of this current in the generation and conduction of the cardiac impulse.

Furthermore, the model's representation of $I_{Kr}$, $I_{to}$, and $I_{Kur}$ was improved, as shown in Figs 6 and 7. Additionally, the sodium/calcium exchanger underwent rigorous fitting

<figure><img src="image_16.png" /><figcaption>**Figure 14.** AP in control condition (grey line), and in the presence of current blocks (black line)
A, 100% $I_{Cav1.3}$. B, 100% $I_{Cav1.2}$ block. C, 100% $I_{Cav3.1}$ block.</figcaption></figure>
efforts based on Choisy's data (Choisy et al., 2015) by
changing its formulation, as depicted in Fig. 8.

I<sub>Kr</sub> was found in all AVNCs tested by Marger et al. (Marger et al., 2011b), suggesting that this current is important also in this cell type. The application of 1 µM E-4031 in an AVN preparation increased the CL by 28 ± 3.4% from 452 ± 20 ms to 580 ± 35.8 ms (Nikmaram et al., 2008). We have simulated a 100% current block (Fig. 15A) and obtained cessation of beating, supporting the hypothesis of a central role for I<sub>Kr</sub> in AVNCs repolarization. The difference between simulation and experiment can be ascribed to both the longer basal CL measured by Nikmaram et al. in comparison to the model (452 vs. 286 ms) and to the repolarizing contributions of atrial cells expressing I<sub>K1</sub> which are present in the tissue preparation. In fact, by testing a block percentage lower than 80% we were also able to prolong the CL, without stopping the beating.

The 4-AP sensitive currents can be important in
regulating the APD in AVN cells, and their block was
experimentally tested by Nikmaram (Nikmaram et al.,
2008), who has shown a mild prolongation of the CL

upon 4-AP administration. Such prolongation might be
due to APD prolongation; however, no specific data
were available. Our simulations did not show any effects
on APD and consequently on CL (Fig. 15B). Finally,
we investigated the modulatory effect of SR calcium
release and uptake on beating rate and AP properties. In
particular, blocks of 50/100% and an increase of 50% were
tested for both Ca²⁺ release and uptake fluxes (J<sub>rel</sub> and J<sub>up</sub>):
there were only negligible effects on the AP whereas the
CaT90 was prolonged with both J<sub>rel</sub> and J<sub>up</sub> block (204 ms
and 250 ms). To our knowledge, there are no experimental
data on AVN mouse cells to confirm these results, which
remain a prediction of the model.

As far as we know, the only other model that describes
the AVN cell behaviour has been published by Inada et al.
(Inada et al., 2009) in rabbit species. The focus of Inada's
work was the propagation of electrical signals through
the AVN and for this reason, the authors developed a
one-dimensional multicellular model including the SAN,
AVN, and atrial cells. Their model can reproduce different
experimental findings, e.g. the conduction velocity of the
AVN, the effects of different drugs on AVN function, such

<figure><img src="image_17.png" /><figcaption>Figure 15. AP in control condition (grey line) and in the presence of current blocks (black line)
A, 100% I<sub>Kr</sub> block. B, 50% I<sub>Kur</sub>.</figcaption></figure>
as the slowing of AVN conduction by beta-blockers, and different types of arrhythmias, like atrial fibrillation and atrioventricular nodal reentry tachycardia. As a single-cell description of an AVN cell, the present model introduces new and different important elements. Indeed, Inada's model did not include Na⁺ current, since Munk et al. (Munk et al., 1996) and Ren et al. (Fu-Xian et al., 2006) showed experimentally that, in the rabbit, I<sub>Na</sub> is absent in putative AVN cells. Regarding potassium currents, it was assumed that I<sub>to</sub> is absent in the AVN region and I<sub>Kur</sub> current was not included. They considered only I<sub>Kr</sub> to be present in the AVN, while we included I<sub>Ks</sub> and I<sub>K1</sub> too. For the L-type calcium current, Inada's model did not distinguish between the two isoforms (Ca<sub>v</sub>1.3/Ca<sub>v</sub>1.2), while we account for both descriptions. They took into consideration a generic background current, I<sub>b</sub>, while we have distinguished the three ions (I<sub>bNa</sub>, I<sub>bCa</sub>, I<sub>bK</sub>). For simplicity, intracellular Na⁺ and K⁺ concentrations were assumed to be constant.

## Limitations and future developments

Like any mathematical model, our AVN cell model relies on simplifications and assumptions to represent the complex biological reality. These assumptions could potentially overlook intricate interactions and nuanced behaviours present in real cells. For example, detailed simulation of the effects of activation of β-adrenoceptors, which is not included in the present model, would be an important step forward. This model development would require previous descriptions of the effects of adrenergic activation on AVN ionic currents. The model primarily draws upon data from mouse AVN cells. However, interspecies differences could influence ion channel kinetics and cellular behaviour. Extrapolating findings from mouse cells to other species or human cells might require further validation and adjustments. The availability of high-quality experimental data for parameter fitting is critical. The lack of comprehensive data for all ionic currents, especially in specific conditions or disease states, could limit the accuracy of our model's predictions (e.g. the relative contribution of I<sub>to</sub>/I<sub>Kur</sub> should be obtained from the peak vs. end-pulse current, but there are no such data available yet). The model currently focuses on individual ionic currents, overlooking potential coregulation in expression of different currents in a cell-specific way or their modulation by intracellular signalling pathways and connectivity with atrial cells. These interactions could underlie emergent properties of AVN pacemaking and impulse conduction at the cellular and tissue-wide levels, properties that we did not attempt to study in the present work.

## Conclusion

In conclusion, the present model represents a significant step forward in understanding the complexity of the AVN physiology and pacemaker mechanism in the cardiac conduction system. By meticulously refining the formulations of sarcolemmal currents flowing through ionic channels, pumps, and exchangers, we have achieved a more accurate representation of AVN cell behaviour.

This model not only bridges historical insights, as exemplified by Mangoni et al.'s work, but also integrates contemporary experimental data from Marger et al. and Choisy et al. (Choisy et al., 2015; Marger et al., 2011a), bringing our understanding closer to reality. In particular, our computations highlight important roles for L-type Ca<sub>v</sub>1.3 and hyperpolarization-activated 'funny' channels underlying I<sub>f</sub> in the generation of pacemaking in the AVN. It also accounts for the experimentally observed role of TTX-resistant and TTX-sensitive Na⁺ currents I<sub>Na,r</sub> and I<sub>Na,s</sub> in AVN automaticity.

The model's accurate representation of ion channel behaviour and its influence on action potentials and pacemaker activity opens the possibility of obtaining predictive insights into cardiac dysfunctions and potential therapeutic interventions. Moving forward, efforts to refine the model with additional data, validate it across different species, and accounting for more intricate cellular dynamics will be essential for harnessing its full potential in advancing our understanding of AVN physiology and cardiac electrophysiology as a whole.

## Appendix

List of model equations

Cell dimension

$$C_m = 22 \ [pF]$$

$$R_{cell} = 0.00005 \ [dm]$$

$$Cm_{surf} = 1 \left[ \frac{\mu F}{cm^2} \right]$$

$$L_{cell} = \frac{1}{2\pi R_{cell}} \left( \frac{Cm}{Cm_{surf} \cdot 10^8} - 2\pi R_{cell}^2 \right) \ [dm]$$

$$V_{cell} = 3.1416 \cdot R_{cell}^2 \cdot L_{cell}$$

$$V_{sub} = 0.01 \cdot V_{cell}$$

$$V_i = 0.46 \cdot V_{cell} - V_{sub}$$

$$V_{JSR} = 0.0012 \cdot V_{cell}$$

$$V_{NSR} = 0.0116 \cdot V_{cell}$$

Ionic concentrations (mM)

$$Na_o = 140, K_o = 5.4, Ca_o = 2$$
Nernst potential (mV)

$$E_{Na} = \frac{RT}{F} \cdot \log\left(\frac{Na_o}{Na_i}\right)$$

$$E_{Ca} = 0.5 \frac{RT}{F} \cdot \log\left(\frac{Ca_o}{Ca_{sub}}\right)$$

$$E_{K} = \frac{RT}{F} \cdot \log\left(\frac{K_{o}}{K_{i}}\right)$$

$$E_{Ks} = \frac{RT}{F} \cdot \log\left(\frac{K_{o} + 0.1 \cdot Na_{o}}{K_{i} + 0.1 \cdot Na_{i}}\right)$$
$$I_{Nas}$$

$$I_{Nas} = g_{Nas} * m_{s}^{3} * h_{s} * Na_{o} * \frac{F^{2}}{RT}$$

$$* \frac{e^{(V-E_{Na})*\frac{F}{RT}} - 1}{e^{V*\frac{F}{RT}} - 1} * V \left[ \frac{pA}{pF} \right]$$

$$g_{Nas} = 1.5 * 1.9318e^{-05} \left[ \frac{nS}{pF} \right]$$

$$m_{s,\infty} = \frac{1}{\left(1 + e^{-\frac{(V+38)}{6}}\right)}$$

$$\tau_{ms} = 0.2378$$

$$* \left( \frac{0.6247}{0.832 \cdot e^{-0.335 \cdot (V+56.7)} + 0.627 \cdot e^{0.082 \cdot (V+65.01)}} + 0.04 \right)$$

$$\frac{dm_s}{dt} = \frac{m_{s,\infty} - m_s}{\tau_{ms}}$$

$$h_{s1,\infty} = \frac{1}{\left(1 + e^{-\frac{(V+56)}{3}}\right)}$$

$$\tau_{hs1} = 0.2378 * 0.113$$

$$* \left( \frac{1}{13475.066 \cdot e^{\frac{(V-59.398)}{15.645}}} + \frac{1}{1.113 + 0.044 \cdot e^{\frac{-(V-86.768)}{8.059}}} \right)^{-1}$$

$$\frac{dh_{s1}}{dt} = \frac{hs_{1\infty} - h_{s1}}{\tau_{hs1}}$$

$$h_{s2,\infty} = h_{s1,\infty}$$

$$\tau_{hs2} = 0.2378 * 0.125$$

$$* \left( \frac{1}{140557.232 \cdot e^{\frac{(V-59.455)}{12}}} + \frac{1}{2.471 + 0.767 \cdot e^{\frac{-(V-68.931)}{18.237}}} \right)^{-1}$$

$$\frac{dh_{s2}}{dt} = \frac{hs_{2\infty} - h_{s2}}{\tau_{hs2}}$$

$$F_{Nas} = \frac{0.0952 * e^{-0.063*(V+34.4)}}{1 + 1.66 + e^{-0.225*(V+63.7)}} + 0.0869$$

$$h_s = (1 - F_{Nas}) * h_{s1} + F_{Nas} * h_{s2}$$

$$I_{Nar}$$

$$I_{Nar} = g_{Nar} * m^3 * h * Na_o * \frac{F^2}{RT}$$

$$* \frac{e^{(V-E_{Na})*\frac{F}{RT}} - 1}{e^{V*\frac{F}{RT}} - 1} * V \left[ \frac{pA}{pF} \right]$$

$$g_{Nar} = 1.5 * 4.25e^{-05} \left[ \frac{nS}{pF} \right]$$

$$m_{\infty} = \left( \frac{1}{1 + e^{-\frac{V+41}{4}}} \right)$$

$$\tau_m = \frac{0.6247 * 0.2378}{0.832 \cdot e^{-0.335 \cdot (V+56.7)} + 0.627 \cdot e^{0.082 \cdot (V+65.01)}} + 0.04$$

$$\frac{dm}{dt} = \frac{m_{\infty} - m}{\tau_m}$$

$$h_{1,\infty} = \frac{1}{\left(1 + e^{\frac{(V+65.11)}{6}}\right)}$$

$$\tau_{h1} = \tau_{hs1}$$

$$\frac{dh_1}{dt} = \frac{h_{1\infty} - h_1}{\tau_{h1}}$$

$$h_{2,\infty} = h_{1,\infty}$$

$$\tau_{h2} = \tau_{hs2}$$

$$\frac{dh_2}{dt} = \frac{h_{2\infty} - h_2}{\tau_{h2}}$$

$$F_{Na} = \frac{0.0952 * e^{-0.063*(V+34.4)}}{1 + 1.66 + e^{-0.225*(V+63.7)}} + 0.0869$$

$$h = (1 - F_{Na}) * h_1 + F_{Na} * h_2$$

$$\frac{dV}{dt} = -(I_{Nar} + I_{Nas} + I_{CaL} + I_{CaD} + I_{CaT} + I_{K1} + I_{Kr} + I_{to} + I_{sus} + I_{f} + I_{Ks} + I_{bNa} + I_{bK} + I_{bCa} + I_{NaK} + I_{NaCa} + I_{st} + I_{SK} + I_{KACh}) + I_{stim}$$

$$I_{CaD} = g_{CaD} * \left( f_D * d_D + \frac{0.006}{1 + e^{-\frac{V+14.1}{6}}} \right) * (V - E_{Ca}) \left[ \frac{pA}{pF} \right]$$

$$I_{CaD} = g_{CaD} * \left( f_D * d_D + \frac{0.006}{1 + e^{-\frac{V+14.1}{6}}} \right) * (V - E_{Ca}) \left[ \frac{pA}{pF} \right]$$

$$g_{CaD} = 0.1365 \left[ \frac{nS}{pF} \right]$$
$$d_{D,\infty} = \frac{1}{1 + e^{-\frac{V+22}{4.3}}}$$
$$\alpha_{d_D} = -\frac{0.01419 * (V + 35)}{e^{-\frac{V+35}{2.5}} - 1} - \frac{0.04245 * V}{e^{-0.208*V} - 1}$$

$$\beta_{d_D} = \frac{0.00571 * (V - 5)}{e^{0.4*(V-5)} - 1}$$

$$\tau_{d_D} = \frac{1}{\alpha_{d_D} + \beta_{d_D}}$$

$$\frac{dd_D}{dt} = \frac{d_{D,\infty} - d_D}{\tau_{d_D}}$$

$$f_{D,\infty} = \frac{1}{1 + e^{\frac{V+48}{5.4}}}$$

$$\alpha_{f_D} = \frac{0.00312 * (V + 68)}{e^{\frac{V+68}{4}} - 1}$$

$$\beta_{f_D} = \frac{0.025}{1 + e^{-\frac{V+68}{4}}}$$

$$\tau_{f_D} = \frac{1}{\alpha_{f_D} + \beta_{f_D}}$$

$$\frac{df_D}{dt} = \frac{f_{D,\infty} - f_D}{\tau_{f_D}}$$

$$I_{CaL} (Cav\ 1.2)$$

$$I_{CaL} = g_{CaL} * \left( f_L * d_L + \frac{0.006}{1 + e^{-\frac{V+14.1}{6}}} \right) * (V - E_{Ca}) \left[ \frac{pA}{pF} \right]$$

$$g_{CaL} = 0.07 \left[ \frac{nS}{pF} \right]$$

$$d_{L,\infty} = \frac{1}{1 + e^{-\frac{V+10}{4}}}$$

$$\alpha_{d_L} = -\frac{0.01419 * (V + 35)}{e^{-\frac{V+35}{2.5}} - 1} - \frac{0.04245 * V}{e^{-0.208*V} - 1}$$

$$\beta_{d_L} = \frac{0.00571 * (V - 5)}{e^{0.4*(V-5)} - 1}$$

$$\tau_{d_L} = \frac{1}{\alpha_{d_L} + \beta_{d_L}}$$

$$\frac{dd_L}{dt} = \frac{d_{L,\infty} - d_L}{\tau_{d_L}}$$

$$f_{L,\infty} = \frac{1}{1 + e^{\frac{V+36}{4.6}}}$$

$$\alpha_{f_L} = \frac{0.00312 * (V + 68)}{e^{\frac{V+68}{4}} - 1}$$

$$\beta_{f_L} = \frac{0.025}{1 + e^{-\frac{V+68}{4}}}$$

$$\tau_{f_L} = \frac{1}{\alpha_{f_L} + \beta_{f_L}}$$

$$\frac{d f_L}{dt} = \frac{f_{L,\infty} - f_L}{\tau_{f_L}}$$

$$I_{CaT} (Cav\ 3.1)$$

$$I_{CaT} = g_{CaT} * d_T * f_T * (V - E_{Ca}) \left[ \frac{pA}{pF} \right]$$

$$g_{CaT} = 0.2658 \left[ \frac{nS}{pF} \right]$$

$$d_{Tinf} = \frac{1}{1 + e^{-\frac{V+40}{6}}}$$

$$\alpha_{d_T} = 1.068 * e^{\frac{V+26.3}{30}}$$

$$\beta_{d_T} = 1.068 * e^{-\frac{V+26.3}{30}}$$

$$\tau_{d_T} = \frac{1}{\alpha_{d_T} + \beta_{d_T}}$$

$$\frac{dd_T}{dt} = \frac{d_{T,\infty} - d_T}{\tau_{d_T}}$$

$$f_{Tinf} = \frac{1}{1 + e^{\frac{V+71}{3.4}}}$$

$$\alpha_{f_T} = 0.0153 * e^{-\frac{V+71.7}{83.3}}$$

$$\beta_{f_T} = 0.015 * e^{\frac{V+71.1}{15.38}}$$

$$\tau_{f_T} = \frac{1}{\alpha_{f_T} + \beta_{f_T}}$$

$$\frac{df_T}{dt} = \frac{f_{T,\infty} - f_T}{\tau_{f_T}}$$

$$I_{K1}$$

$$I_{K1} = g_{K1} * \frac{V - E_K}{1 + e^{0.07*(V+80)}} \left[ \frac{pA}{pF} \right]$$

$$g_{K1} = 0.0693 \left[ \frac{nS}{pF} \right]$$

$$I_{Kr}$$

$$I_{Kr} = g_{Kr} * P_{af} * P_i * (V - E_K) \left[ \frac{pA}{pF} \right]$$

$$g_{Kr} = 0.0682 \left[ \frac{nS}{pF} \right]$$

$$P_{af,\infty} = \frac{1}{1 + e^{-\frac{(V+23)}{6.5}}}$$

$$\tau_{af} = \frac{0.3256}{0.0372 * e^{\frac{V-15}{20}} + 0.00096 * e^{-\frac{V+40}{5}}} + 5,$$
$$ \frac{dP_{af}}{dt} = \frac{P_{af,\infty} - P_{af}}{\tau_{af}} $$

$$ P_{i\infty} = \frac{1}{1 + e^{\frac{V+15}{6.5}}} $$

$$ \tau_i = 2 $$

$$ \frac{dP_i}{dt} = \frac{P_{i,\infty} - P_i}{\tau_i} $$

$$ I_{to} \& I_{Kur} $$

$$ I_{to} = g_{to} * q * r * (V - E_K) \left[ \frac{pA}{pF} \right] $$

$$ g_{to} = 0.0402 \left[ \frac{nS}{pF} \right] $$

$$ I_{Kur} = g_{Kur} * r * (V - E_K) \left[ \frac{pA}{pF} \right] $$

$$ g_{Kur} = 0.0151 \left[ \frac{nS}{pF} \right] $$

$$ q_{\infty} = \frac{1}{1 + e^{\frac{V+59.37}{13.1}}} $$

$$ \tau_q = 10.1 + \frac{65.17}{0.57 * e^{-0.08*(V+49)} + 0.000024 * e^{0.1*(V+50.93)}} $$

$$ \frac{dq}{dt} = \frac{q_{\infty} - q}{\tau_q} $$

$$ r_{\infty} = \frac{1}{1 + e^{\frac{-(V-10.93)}{19.7}}} $$

$$ \tau_r = 2.98 + \frac{15.59}{1.037 * e^{0.09*(V+30.61)} + 0.369 * e^{-0.12*(V+23.84)}} $$

$$ \frac{dr}{dt} = \frac{r_{\infty} - r}{\tau_r} $$

$$ I_{Ks} $$

$$ I_{Ks} = g_{Ks} * xs^2 * (V - E_{Ks}) \left[ \frac{pA}{pF} \right] $$

$$ g_{Ks} = 1.88 \left[ \frac{nS}{pF} \right] $$

$$ \alpha_{xs} = \frac{0.014}{1 + e^{-\frac{V-40}{9}}} $$

$$ \beta_{xs} = 0.001 * e^{-\left(\frac{V}{45}\right)} $$

$$ x * s_{\infty} = \frac{\alpha_{xs}}{\alpha_{xs} + \beta_{xs}} $$

$$ \tau_{xs} = \frac{1}{\alpha_{xs} + \beta_{xs}} $$

$$ \frac{dxs}{dt} = \frac{xs_{\infty} - xs}{\tau_{xs}} $$

$$ I_f \text{ (from Kharche with } \tau_P \text{ mod)} $$

$$ I_f = i_{fNa} + i_{fK} $$

$$ i_{fNa} = 0.3833 * gh * P * (V - E_K) \left[ \frac{pA}{pF} \right] $$

$$ i_{fK} = 0.3833 * gh * P * (V - E_{Na}) \left[ \frac{pA}{pF} \right] $$

$$ g_h = 0.2073 \left[ \frac{nS}{pF} \right] $$

$$ P_{\infty} = \frac{1}{1 + e^{\frac{V+97}{20}}} $$

$$ \tau_p = \frac{1.505}{e^{-0.0119(V+604.3)} + e^{\frac{V-41}{10}}} $$

$$ \frac{dP}{dt} = \frac{P_{\infty} - P}{\tau_p} $$

## Background currents

$$ I_{bNa} $$

$$ I_{bNa} = g_{bNa} * (V - E_{Na}) \left[ \frac{pA}{pF} \right] $$

$$ g_{bNa} = 0.0053 \left[ \frac{nS}{pF} \right] $$

$$ I_{bK} $$

$$ I_{bK} = g_{bK} * (V - E_K) \left[ \frac{pA}{pF} \right] $$

$$ g_{bK} = 0.0011 \left[ \frac{nS}{pF} \right] $$

$$ I_{bCa} $$

$$ I_{bCa} = g_{bCa} * (V - E_{Ca}) \left[ \frac{pA}{pF} \right] $$

$$ g_{bCa} = 3e^{-4} \left[ \frac{nS}{pF} \right] $$

$$ I_{NaK} $$

$$ I_{NaK} = \frac{1}{C_m} * I_{NaK_{max}} * \left( \frac{Na_i}{K_{mNa} + Na_i} \right)^3 * \left( \frac{K_o}{K_{mK} + K_o} \right)^2 * \frac{1.6}{1.5 + e^{-\frac{V+60}{40}}} \left[ \frac{pA}{pF} \right] $$

$$ I_{NaK_{max}} = 143.4 \left[ pA \right] $$
$$K_{mNa} = 5.64 \ [mM]$$

$$K_{mK} = 0.621 \ [mM]$$

$$I_{NaCa}$$

$$I_{SK} = g_{SK} * m_{cak}^2 * (V - E_{K}) \left[ \frac{pA}{pF} \right]$$

$$I_{NaCa} = K_{NaCa} * \frac{(e^{\gamma_{NaCa}*V*\frac{E}{RT}} * Nai^3 * Cao - e^{((\gamma_{NaCa}-1)*V*\frac{E}{RT})} * Nao^3 * Casub * \alpha)}{(K_{mNai}^3 + Nao^3) * (K_{mCa} + Cao) * (1 + K_{sat} * e^{((\gamma_{NaCa}-1)*V*\frac{E}{RT})})} \left[ \frac{pA}{pF} \right]$$

$$K_{NaCa} = 7000$$

$$\gamma_{NaCa} = 0.35$$

$$\alpha = 2.5$$

$$K_{mNai} = 87.5 \ (mM)$$

$$K_{mCa} = 1.38 \ (mM)$$

$$K_{sat} = 0.1$$

$$I_{st}$$

$$I_{st} = g_{st} * d_s * f_s * (V - E_{st}) \left[ \frac{pA}{pF} \right], g_{st} = 0.1545 \left[ \frac{nS}{pF} \right]$$

$$\alpha_{d_s} = \frac{1}{0.15 * e^{-\frac{V}{11}} + 0.2 * e^{-\frac{V}{700}}}$$

$$\beta_{d_s} = \frac{1}{16 * e^{\frac{V}{8}} + 15 * e^{\frac{V}{50}}}$$

$$d_{s\infty} = \frac{\alpha_{d_s}}{\alpha_{d_s} + \beta_{d_s}}$$

$$\tau_{d_s} = \frac{1}{\alpha_{d_s} + \beta_{d_s}}$$

$$\frac{dd_s}{dt} = \frac{d_{s,\infty} - d_s}{\tau_{d_s}}$$

$$\alpha_{f_s} = \frac{1}{3100 * e^{-\frac{V}{13}} + 700 * e^{-\frac{V}{70}}}$$

$$\beta_{f_s} = \frac{1}{95 * e^{-\frac{V}{10}} + 50 * e^{-\frac{V}{700}}} + \frac{0.000229}{1 + e^{-\frac{V}{5}}}$$

$$f_{s\infty} = \frac{\alpha_{f_s}}{\alpha_{f_s} + \beta_{f_s}}$$

$$\tau_{f_s} = \frac{1}{\alpha_{f_s} + \beta_{f_s}}$$

$$\frac{df_s}{dt} = \frac{f_{s,\infty} - f_s}{\tau_{f_s}}$$

$$I_{SK}$$

$$g_{SK} = 0.004 \left[ \frac{nS}{pF} \right]$$

$$m_{cak_{\infty}} = \frac{Ca_r}{1 + Ca_r}$$

$$\tau_{cak} = \frac{0.001}{\text{beta } a_k * (1 + Car)}, \text{ bet } a_k = 100$$

$$Ca_r = \left( \frac{Ca_{sub}}{Ca_c} \right)^{n_{SK}}, \quad C a_c = 0.7e^{-5} \ (mM), \ n_{SK} = 1.7$$

$$I_{KACh}$$

$$I_{KACh} = g_{KACh} (K_i - (K_o * e^{-\frac{VF}{RT}}) \left[ \frac{pA}{pF} \right])$$

$$g_{KACh} = 0.0 \cdot K_o^{0.41} \left[ \frac{nS}{pF} \right]$$

Intracellular calcium dynamics

$$j_{up} = P_{up} \frac{(Ca_i/K_{mf})^{n_{up}} - (Ca_{NSR}/K_{mr})^{n_{up}}}{1 + (Ca_i/K_{mf})^{n_{up}} - (Ca_{NSR}/K_{mr})^{n_{up}}} \left[ \frac{mM}{ms} \right],$$

$$P_{up} = 0.04 \left[ \frac{mM}{ms} \right], \ K_{mf} = 0.000246 \ [mM],$$

$$K_{mr} = 3.29 \ [mM], \ n_{up} = 2$$

$$j_{tr} = \frac{Ca_{NSR} - Ca_{JSR}}{\tau_{tr}} \left[ \frac{mM}{ms} \right], \ \tau_{tr} = 60 \left[ \frac{1}{ms} \right]$$

$$j_{Ca_{dif}} = \frac{Ca_{sub} - Ca_i}{\tau_{Ca_{dif}}} \left[ \frac{mM}{ms} \right], \ \tau_{Ca_{dif}} = 0.4 \left[ \frac{1}{ms} \right]$$

$$j_{rel} = ks * O1 * (Ca_{JSR} - Ca_{sub}) \left[ \frac{mM}{ms} \right], \ ks = 250 \left[ \frac{1}{ms} \right]$$

$$k_{CaSR} = Max_{SR} - \left( \frac{(Max_{SR} - Min_{SR})}{1 + \left( \frac{EC_{50SR}}{Ca_{JSR}} \right)^{HSR}} \right)$$

$$Max_{SR} = 15$$

$$Min_{SR} = 1$$

$$EC_{50SR} = 0.45 \ [mM]$$
$$HSR = 2.5$$

$$ko_{SRCa} = \frac{ko_{Ca}}{k_{CaSR}} \left[ \frac{1}{mM^2 \cdot ms} \right]$$

$$ko_{Ca} = 1.5 \left[ \frac{1}{mM^2 \cdot ms} \right]$$

$$ko_m = 0.06 \left[ \frac{1}{ms} \right]$$

$$ki_{Ca} = 0.05 \left[ \frac{1}{mM \cdot ms} \right]$$

$$ki_m = 0.005 \left[ \frac{1}{mM \cdot ms} \right]$$

$$ki_{SRCa} = ki_{Ca} * k_{CaSR} \left[ \frac{1}{mM \cdot ms} \right]$$

$$\frac{dR1}{dt} = ki_m \cdot RI - ki_{SRCa} \cdot Ca_{sub} \cdot R1 - (ko_{SRCa} \cdot Ca_{sub}^2 \cdot R1 - ko_m \cdot O1)$$

$$\frac{dO1}{dt} = ko_{SRCa} Ca_{sub}^2 \cdot R1 - ko_m \cdot O1 - (ki_{SRCa} \cdot Ca_{sub} \cdot O1 - ki_m \cdot I1)$$

$$\frac{dI1}{dt} = ki_{SRCa} \cdot Ca_{sub} \cdot O1$$

$$-ki_m \cdot I1 - (ko_m \cdot I1 - ko_{SRCa} \cdot Ca_{sub}^2 \cdot RI)$$

Intracellular ionic concentrations

$$\frac{d f_{TC}}{dt} = k_{fTC} [Ca^{2+}]_i (1 - f_{TC}) - k_{bTC} f_{TC}$$

$$\frac{d f_{TMC}}{dt} = k_{fTMC} [Ca^{2+}]_i (1 - f_{TMC} - f_{TMM}) - k_{bTMC} f_{TMC}$$

$$\frac{d f_{TMM}}{dt} = k_{fTMM} [Mg^{2+}]_i (1 - f_{TMM} - f_{TMM}) - k_{bTMM} f_{TMM}$$

$$\frac{d f_{CMi}}{dt} = k_{fCM} [Ca^{2+}]_i (1 - f_{CMi}) - k_{bCM} f_{CMi}$$

$$\frac{d f_{CMs}}{dt} = k_{fCM} [Ca^{2+}]_{sub} (1 - f_{CMs}) - k_{bCM} f_{CMs}$$

$$\frac{d f_{CQ}}{dt} = k_{fCQ} [Ca^{2+}]_{NSR} (1 - f_{CQ}) - k_{bCQ} f_{CQ}$$

$$k_{fTC} = 88.8 \left[ \frac{1}{mM \cdot ms} \right], k_{bTC} = 0.446 \left[ \frac{1}{ms} \right]$$

$$k_{fTMC} = 227.7 \left[ \frac{1}{mM \cdot ms} \right], k_{bTMC} = 0.00751 \left[ \frac{1}{ms} \right]$$

$$k_{fTMM} = 2.277 \left[ \frac{1}{mM \cdot ms} \right], k_{bTMM} = 0.751 \left[ \frac{1}{ms} \right]$$

$$k_{fCM} = 227.7 \left[ \frac{1}{mM \cdot ms} \right], k_{bCM} = 0.542 \left[ \frac{1}{ms} \right]$$

$$k_{fCQ} = 0.534 \left[ \frac{1}{mM \cdot ms} \right], k_{bCQ} = 0.445 \left[ \frac{1}{ms} \right]$$

$$[Mg^{2+}]_i = 2.5 (mM)$$

$$\frac{d[Na_i^+]}{dt} = -(I_{fNa} + I_{st} + I_{Nas} + I_{Nar} + I_{bNa} + 3 * I_{NaK} + 3 * I_{NaCa}) * \frac{C_m}{F * (V_i + V_{sub})}$$

$$\frac{d[K_i^+]}{dt} = -(I_{Kr} + I_{Ks} + I_{to} + I_{sus} + I_{SK} + I_{K1} + I_{fK} + I_{KACh} - 2I_{NaK} + I_{bK} - I_{stim}) * \frac{C_m}{F * (V_i + V_{sub})}$$

$$\frac{d[Ca_i^{2+}]}{dt} = \frac{(j_{Ca_{dif}} * V_{sub} - j_{up} * V_{NSR})}{V_i} - (d f_{CMi} \cdot CM_{tot} + d f_{TC} \cdot TC_{tot} + f_{TMC} \cdot TMC_{tot})$$

$$\frac{d[Ca_{sub}^{2+}]}{dt} = -\frac{((I_{CaL} + I_{CaT} + I_{bCa} + I_{CaD} - 2I_{NaCa}) * \frac{C_m}{2*F} + j_{rel} * V_{JSR})}{V_{sub}} - (j_{Ca_{dif}} + CM_{tot} \delta a_{fCMs})$$

$$\frac{d[Ca_{NSR}^{2+}]}{dt} = j_{up} - j_{tr} * \frac{V_{JSR}}{V_{NSR}}$$

$$\frac{d[Ca_{JSR}^{2+}]}{dt} = j_{tr} - (j_{rel} + CQ_{tot} * \delta a_{fCQ})$$

$$CM_{tot} = 0.045 [mM]$$

$$TC_{tot} = 0.031 [mM]$$

$$TM C_{tot} = 0.062 [mM]$$

$$CQ_{tot} = 10 [mM]$$
<figure><img src="image_18.png" /><figcaption>Figure A1. Comparison of computed diastolic membrane voltage and main ionic currents</figcaption></figure>

Comparison of computed diastolic membrane voltage and main ionic currents in control (Aa and Ab, respectively) and in 100% I<sub>Cav1.2</sub> block (Ba and Bb) conditions. [Colour figure can be viewed at wileyonlinelibrary.com]

<table><caption>Table A1. List of currents abbreviations</caption><tr><td>I<sub>Na,s</sub></td><td>Neuronal tetrodotoxin (TTX)-sensitive sodium current</td></tr><tr><td>I<sub>Na,r</sub></td><td>Cardiac TTX-resistant voltage-dependent sodium current</td></tr><tr><td>I<sub>st</sub></td><td>Sustained inward current</td></tr><tr><td>I<sub>bNa</sub></td><td>Sodium background current</td></tr><tr><td>I<sub>f</sub></td><td>Funny current</td></tr><tr><td>I<sub>NaCa</sub></td><td>Sodium/calcium exchange current</td></tr><tr><td>I<sub>NaK</sub></td><td>Sodium/potassium ATPase current</td></tr><tr><td>I<sub>Kr</sub></td><td>Rapid delayed rectifier potassium current</td></tr><tr><td>I<sub>Ks</sub></td><td>Slow delayed rectifier potassium current</td></tr><tr><td>I<sub>K1</sub></td><td>Inward rectifier potassium current</td></tr><tr><td>I<sub>to</sub></td><td>Transient outward potassium current</td></tr><tr><td>I<sub>Kur</sub></td><td>Ultra-rapid delayed rectifier potassium current</td></tr><tr><td>I<sub>SK</sub></td><td>Calcium-dependent-potassium current</td></tr><tr><td>I<sub>bK</sub></td><td>Potassium background current</td></tr><tr><td>I<sub>Cav1.3</sub></td><td>L-type 1.3 calcium current</td></tr><tr><td>I<sub>Cav1.2</sub></td><td>L-type 1.2 calcium current</td></tr><tr><td>I<sub>Cav3.1</sub></td><td>T-type calcium current</td></tr><tr><td>I<sub>bCa</sub></td><td>Calcium background current</td></tr><tr><td>J<sub>rel</sub></td><td>SR calcium release flux</td></tr><tr><td>J<sub>tr</sub></td><td>Calcium translocation from NSR to JSR flux</td></tr><tr><td>J<sub>up</sub></td><td>Calcium uptake via SERCA pump flux</td></tr><tr><td>J<sub>diffCa</sub><sup>2+</sup></td><td>Calcium diffusion flux</td></tr></table>

<figure><img src="image_19.png" /><figcaption>Figure A2. Voltage dependence of the model's I<sub>Cav3.1</sub> inactivation time constant</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure A3. Simulations results in control and 100% I<sub>Cav1.2</sub> block</figcaption></figure>

Simulations results in control (black) and 100% I<sub>Cav1.2</sub> block (grey dashed line) conditions: membrane voltage (top), T-type Ca<sup>2+</sup> current (second panel), T-type Ca<sup>2+</sup> current inactivation gate (third panel) and T-type Ca<sup>2+</sup> current inactivation gate time constant (bottom).
<figure><img src="image_21.png" /><figcaption>**Figure A4.** P1/P2 protocol results on sodium currents</figcaption></figure>

After 2 ms, both $I_{Na}$ components have recovered from inactivation by more than 90%, and after 3 ms by 100%. The inset shows the protocol adopted. [Colour figure can be viewed at wileyonlinelibrary.com]

<figure><img src="image_22.png" /><figcaption>**Figure A5.** S1S2 protocol results</figcaption></figure>

The reference action potential was obtained in a stimulated condition at 450 bpm. A current stimulus of 400 nA for 1 ms was employed to excite the model. In this condition, a full AP is obtained for diastolic intervals (DI) longer than 7 ms. [Colour figure can be viewed at wileyonlinelibrary.com]

## References

Amuzescu, B., Airini, R., Epureanu, F. B., Mann, S. A., Knott, T., & Radu, B. M. (2021). Evolution of mathematical models of cardiomyocyte electrophysiology. *Mathematical Biosciences, **334***, 108567.

Cavero, I., & Holzgrefe, H. (2022). Remembering the canonical discoverers of the core components of the mammalian cardiac conduction system: Keith and Flack, Aschoff and Tawara, His, and Purkinje. *Advances in Physiology Education, **46*** (4), 549–579.

Choisy, S. C., Cheng, H., Orchard, C. H., James, A. F., & Hancox, J. C. (2015). Electrophysiological properties of myocytes isolated from the mouse atrioventricular node: L-type ICa, IKr, $I_f$, and Na-Ca exchange. *Physiological Reports, **3*** (11), 1–18.

Dobrzynski, H., Nikolski, V. P., Sambelashvili, A. T., Greener, I. D., Yamamoto, M., Boyett, M. R., & Efimov, I. R. (2003). Site of origin and molecular substrate of atrioventricular junctional rhythm in the rabbit heart. *Circulation Research, **93*** (11), 1102–1110.

Fu-Xian, R., Niu, X.-L., Yan, O. U., Zhen-Hua, H., Ling, F.-D., Zhou, S.-S., & Li, Y.-J. (2006). Morphological and electro-physiological properties of single myocardial cells from Koch triangle of rabbit heart. *Chinese Medical Journal (England), **119**(24), 2075–2084.

Gemmell, P., Burrage, K., Rodriguez, B., & Quinn, T. A. (2014). Population of computational rabbit-specific ventricular action potential models for investigating sources of variability in cellular repolarisation. *PLoS ONE, **9**(2), e90112.

Gong, J. Q. X., & Sobie, E. A. (2018). Population-based mechanistic modeling allows for quantitative predictions of drug responses across cell types. *NPJ Systems Biology and Applications, **4*** (1), 1–11.

Grundy, D. (2015). Principles and standards for reporting animal experiments in The Journal of Physiology and Experimental Physiology. *Journal of Physics, **593**(12), 2547–2549.

Inada, S., Hancox, J. C., Zhang, H., & Boyett, M. R. (2009). One-dimensional mathematical model of the atrio-ventricular node including atrio-nodal, nodal, and nodal-His cells. *Biophysical Journal, **97**(8), 2117–2127.

Kernik, D. C., Morotti, S., Wu, H., Garg, P., Duff, H. J., Kurokawa, J., Jalife, J., Wu, J. C., Grandi, E., & Clancy, C. E. (2019). A computational model of induced pluripotent stem-cell derived cardiomyocytes incorporating experimental variability from multiple data sources. *The Journal of Physiology, **597**(17), 4533–4564.

Kharche, S., Yu, J., Lei, M., & Zhang, H. (2011). A mathematical model of action potentials of mouse sino-atrial node cells with molecular bases. *American Journal of Physiology-Heart and Circulatory Physiology, **301**(3), H945–H963.

Kim, D., Song, I., Keum, S., Lee, T., Jeong, M.-J., Kim, S.-S., Mcenery, M. W., & Shin, H.-S. (2001). Lack of the burst firing of thalamocortical relay neurons and resistance to absence seizures in mice lacking α1G T-type $Ca^{2+}$ channels. *Neuron, **31**(1), 35–45.

Lei, M., Goddard, C., Liu, J., Léoni, A.-L., Royer, A., Fung, S. S.-M., Xiao, G., Ma, A., Zhang, H., Charpentier, F., Vandenberg, J. I., Colledge, W. H., Grace, A. A., & Huang, C. L.-H. (2005). Sinus node dysfunction following targeted disruption of the murine cardiac sodium channel gene *Scn5a*. *The Journal of Physiology, **567**(2), 387–400.

Louradour, J., Bortolotti, O., Torre, E., Bidaud, I., Lamb, N., Fernandez, A., Le Guennec, J.-Y., Mangoni, M. E., & Mesirca, P. (2022). L-Type Cav1.3 calcium channels are required for beta-adrenergic triggered automaticity in dormant mouse sinoatrial pacemaker cells. *Cells, **11*** (7), 1114.
Mangoni, M. E., Couette, B., Bourinet, E., Platzer, J., Reimer, D., Striessnig, J., & Nargeot, J. (2003). Functional role of L-type Ca<sub>v</sub> 1.3 Ca<sup>2+</sup> channels in cardiac pacemaker activity. *Proceedings of the National Academy of Sciences*, **100**(9), 5543–5548.

Mangoni, M. E., & Nargeot, J. (2001). Properties of the hyperpolarization-activated current (I(f)) in isolated mouse sino-atrial cells. *Cardiovascular Research*, **52**(1), 51–64.

Mangoni, M. E., Traboulsie, A., Leoni, A.-L., Couette, B., Marger, L., Le Quang, K., Kupfer, E., Cohen-Solal, A., Vilar, J., Shin, H.-S., Escande, D., Charpentier, F., Nargeot, J., & Lory, P. (2006). Bradycardia and slowing of the atrioventricular conduction in mice lacking CaV3.1/α1G T-type calcium channels. *Circulation Research*, **98**(11), 1422–1430.

Marger, L., Mesirca, P., Alig, J., Torrente, A., Dübel, S., Engeland, B., Kanani, S., Fontanaud, P., Striessnig, J., Shin, H.-S., Isbrandt, D., Ehmke, H., Nargeot, J., & Mangoni, M. E. (2011a). Pacemaker activity and ionic currents in mouse atrioventricular node cells. *Channels*, **5**(3), 241–250.

Marger, L., Mesirca, P., Alig, J., Torrente, A., Dübel, S., Engeland, B., Kanani, S., Fontanaud, P., Striessnig, J., Shin, H.-S., Isbrandt, D., Ehmke, H., Nargeot, J., & Mangoni, M. E. (2011b). Functional roles of Cav1.3, Cav3.1 and HCN channels in automaticity of mouse atrioventricular cells: Insights into the atrioventricular pacemaker mechanism. *Channels*, **5**(3), 251–261.

Marionneau, C., Couette, B., Liu, J., Li, H., Mangoni, M. E., Nargeot, J., Lei, M., Escande, D., & Demolombe, S. (2005). Specific pattern of ionic channel gene expression associated with pacemaker activity in the mouse heart. *The Journal of Physiology*, **562**(1), 223–234.

Meijler, F. L., & Janse, M. J. (1988). Morphology and electro-physiology of the mammalian atrioventricular node. *Physiological Reviews*, **68**(2), 608–647.

Mesirca, P., Alig, J., Torrente, A. G., Müller, J. C., Marger, L., Rollin, A., Marquilly, C., Vincent, A., Dubel, S., Bidaud, I., Fernandez, A., Seniuk, A., Engeland, B., Singh, J., Miquerol, L., Ehmke, H., Eschenhagen, T., Nargeot, J., Wickman, K., ... Mangoni, M. E. (2014). Cardiac arrhythmia induced by genetic silencing of “funny” (f) channels is rescued by GIRK4 inactivation. *Nature Communications*, **5**(1), 4664.

Mesirca, P., Nakao, S., Nissen, S. D., Forte, G., Anderson, C., Trussell, T., Li, J., Cox, C., Zi, M., Logantha, S., Yaar, S., Cartensen, H., Bidaud, I., Stuart, L., Soattin, L., Morris, G. M., da Costa Martins, P. A., Cartwright, E. J., Oceandy, D., ... D'souza, A. (2021). Intrinsic electrical remodeling underlies atrioventricular block in athletes. *Circulation Research*, **129**(1), E1–E20.

Mesirca, P., Torrente, A. G., & Mangoni, M. E. (2015). Functional role of voltage gated Ca²⁺ channels in heart automaticity. *Frontiers in Physiology*, **6**, 19.

Munk, A. A., Adjemian, R. A., Zhao, J., Ogbaghebrael, A., & Shrier, A. (1996). Electrophysiological properties of morphologically distinct cells isolated from the rabbit atrioventricular node. *The Journal of Physiology*, **493**(3), 801–818.

Muszkiewicz, A., Britton, O. J., Gemmell, P., Passini, E., Sánchez, C., Zhou, X., Carusi, A., Quinn, T. A., Burrage, K., Bueno-Orovio, A., & Rodriguez, B. (2016). Variability in

cardiac electrophysiology: Using experimentally-calibrated populations of models to move beyond the single virtual physiological human paradigm. *Progress in Biophysics and Molecular Biology*, **120**(1–3), 115–127.

Nikmaram, M. R., Liu, J., Abdelrahman, M., Dobrzynski, H., Boyett, M. R., & Lei, M. (2008). Characterization of the effects of Ryanodine, TTX, E-4031 and 4-AP on the sino-atrial and atrioventricular nodes. *Progress in Biophysics and Molecular Biology*, **96**(1–3), 452–464.

Paci, M., Casini, S., Bellin, M., Hyttinen, J., & Severi, S. (2018). Large-scale simulation of the phenotypical variability induced by loss-of-function long QT mutations in human induced pluripotent stem cell cardiomyocytes. *International Journal of Molecular Sciences*, **19**(11), 3583.

Paci, M., Passini, E., Klimas, A., Severi, S., Hyttinen, J., Rodriguez, B., & Entcheva, E. (2020). All-optical electrophysiology refines populations of in silico human iPSC-CMs for drug evaluation. *Biophysical Journal*, **118**(10), 2596–2611.

Paci, M., Sartiani, L., del Lungo, M., Jaconi, M., Mugelli, A., Cerbai, E., & Severi, S. (2012). Mathematical modelling of the action potential of human embryonic stem cell derived cardiomyocytes. *Biomedical Engineering Online [Electronic Resource]*, **11**(1), 61.

Ricci, E., Bartolucci, C., & Severi, S. (2023). The virtual sino-atrial node: What did computational models tell us about cardiac pacemaking? *Progress in Biophysics and Molecular Biology*, **177**, 55–79.

Ridley, J. M., Cheng, H., Harrison, O. J., Jones, S. K., Smith, G. L., Hancox, J. C., & Orchard, C. H. (2008). Spontaneous frequency of rabbit atrioventricular node myocytes depends on SR function. *Cell Calcium*, **44**(6), 580–591.

Saeed, Y., Temple, I. P., Borbas, Z., Atkinson, A., Yanni, J., Maczewski, M., Mackiewicz, U., Aly, M., Logantha, S. J. R. J., Garratt, C. J., & Dobrzynski, H. (2018). Structural and functional remodeling of the atrioventricular node with aging in rats: The role of hyperpolarization-activated cyclic nucleotide-gated and ryanodine 2 channels. *Heart Rhythm*, **15**(5), 752–760.

Shampine, L. F., & Reichelt, M. W. (1997). The MATLAB ODE suite. *SIAM Journal on Scientific Computing*, **18**(1), 1–22.

Sinnegger-Brauns, M. J., Hetzenauer, A., Huber, I. G., Renström, E., Wietzorrek, G., Berjukov, S., Cavalli, M., Walter, D., Koschak, A., Waldschütz, R., Hering, S., Bova, S., Rorsman, P., Pongs, O., Singewald, N., & Striessnig, J. (2004). Isoform-specific regulation of mood behavior and pancreatic β cell and cardiovascular function by L-type Ca²⁺ channels. *Journal of Clinical Investigation*, **113**(10), 1430.

Temple, I. P., Logantha, S. J. R. J., Absi, M., Zhang, Y., Pervolaraki, E., Yanni, J., Atkinson, A., Petkova, M., Quigley, G. M., Castro, S., Drinkhill, M., Schneider, H., Monfredi, O., Cartwright, E., Zi, M., Yamanushi, T. T., Mahadevan, V. S., Gurney, A. M., White, E., ... Dobrzynski, H. (2016). Atrioventricular node dysfunction and ion channel transcriptome in pulmonary hypertension. *Circulation: Arrhythmia and Electrophysiology*, **9**(12), e003432.
Torrente, A. G., Mesirca, P., Neco, P., Rizzetto, R., Dubel, S., Barrere, C., Sinegger-Brauns, M., Striessnig, J., Richard, S., Nargeot, J., Gomez, A. M., & Mangoni, M. E. (2016). L-type Ca v 1.3 channels regulate ryanodine receptor-dependent Ca²⁺ release during sino-atrial node pacemaker activity. *Cardiovascular Research, **109**(3), 451-461.

Trayanova, N. A., Lyon, A., Shade, J., & Heijman, J. (2023). Computational modeling of cardiac electrophysiology and arrhythmogenesis. *Physiological Reviews, **104**(3), 1265-1333.

Verheijck, E. E., Wilders, R., Joyner, R. W., Golod, D. A., Kumar, R., Jongsma, H. J., Bouman, L. N., & Ginneken, A. C. G. V. (1998). Pacemaker synchronization of electrically coupled rabbit sinoatrial node cells. *Journal of General Physiology, **111*(1), 95.

Workman, A. J., Kane, K. A., & Rankin, A. C. (2000). Rate-dependency of action potential duration and refractoriness in isolated myocytes from the rabbit AV node and atrium. *Journal of Molecular and Cellular Cardiology, **32**(8), 1525-1537.

Zhang, H., Holden, A. V., Kodama, I., Honjo, H., Lei, M., Varghese, T., & Boyett, M. R. (2000). Mathematical models of action potentials in the periphery and center of the rabbit sinoatrial node. *American Journal of Physiology-Heart and Circulatory Physiology, **279**(1), H397-H421.

## Additional information

### Data availability statement

The original contributions presented in the study are included in the article; further inquiries can be directed to the corresponding author.

### Competing interests

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

### Author contributions

C.B., P.M., M.E.M., and S.S. were responsible for the study conception and design. C.B., C.S.B., and E.R. were responsible

for simulations. C.B., E.R., and S.S. were responsible for data analysis. C.B., P.M., E.R., M.E.M., and S.S. were responsible for data interpretation. C.B., P.M., E.R., M.E.M., and S.S. were responsible for writing the manuscript. P.M., E.T., J.L., M.E.M. were responsible for carrying the experiments. C.B., P.M., E.R., C.S.B., E.T., J.L., M.E.M., and S.S. were responsible for the final approval of the manuscript submitted for publication. All authors agree to be accountable for all aspects of the work. All persons designated as authors qualify for authorship, and all those who qualify for authorship are listed.

### Funding

C.B. and S.S. have received funding from the European Union - NextGenerationEU through the Italian Ministry of University and Research under PNRR - M4C2-I1.3 Project PR_00000019 'HEAL ITALIA' to S.S. CUP J33C22002920006. The views and opinions expressed are those of the authors only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. F.R. and L.D. are the members of the INdAM group GNCS (National Group for the Scientific Computing). Work in the M.E.M. team was supported by Fondation Leducq TNE FANTASY 19CV03.

### Acknowledgements

Open access publishing facilitated by Universita degli Studi di Bologna, as part of the Wiley — CRUI-CARE agreement.

### Keywords

atrio-ventricular node, action potential, cardiac electrophysiology, cardiomyocyte, computer simulation, ion channels, mouse, pacemaker activity

### Supporting information

Additional supporting information can be found online in the Supporting Information section at the end of the HTML view of the article. Supporting information files available:

### Peer Review History