The model consists of a square matrix of 200 × 200 cells representing the discrete, simplified Kirchhoff network model framework (Jæger and Tveito, 2023). The sinoatrial node is located at the centre of the tissue and is modelled as an ellipse having semi-major and semi-minor axes of 75 and 12 cells, respectively (Figure 1). The SAN is surrounded by an insulating border ( R gap = ∞ MΩ ) that prevents interactions with the atrial tissues aside from the five SEPs; these have dimensions of 20 × 11 cells. Considering the Fabbri–Wilders–Severi (FWS) model cell length of 67 µm (Fabbri et al., 2017a), the above dimensions correspond to 1.34 mm, which is compatible with other experimental reports (Csepe et al., 2016) and previous computational investigations (Amsaleg et al., 2022). In addition to the FWS model, we adopted the model by Koivumäki et al. (2011) for the atrial cells (dubbed as the ‘K’ model hereafter) and the model by Morgan et al. (2016) for the fibroblasts. In addition, the behaviour of the Mazhar–Bartolucci–Severi (MBS) atrial model (Mazhar et al., 2023) recently developed by our group was tested. In this work, the latter was simplified by removing the equations describing the mechanoelectric feedback and calmodulin-kinase II activity, reducing it to 29 ordinary differential equations (ODEs) for the description of the atrial action potentials (APs).

A linear gradient in the fibroblasts density was implemented from the left side of the SAN to the atrial tissue. In particular, the central SAN had 45.7% fibroblasts (Csepe et al., 2015, 2017; Kalyanasundaram et al., 2021; Li et al., 2023) substituting the SAN cells. The gradient allowed the SEPs to have a lower fibroblast density of 27.7 ± 16.3% vs 45.7 ± 10.9% for the central SAN, as recently reported (Li et al., 2023). The probability that a fibroblast would substitute a SAN cell depended on the gradient and on a random number between 0% and 100% extracted from a uniform distribution. If the random number exceeded the probability of not having a fibroblast, as given by the gradient (i.e. N rand > (100% − 45.7%) in the central SAN), then the SAN cell was replaced with a fibroblast. The linear gradient dropped to zero 10 cells to the right of the SEP exit, and the fibroblasts were allowed to replace the atrial cells in these areas. However, the remaining atrial tissue was considered to be completely free of fibroblasts. With regard to cellular heterogeneity, the maximal conductances of the main ionic currents ( P CaL , P CaT , g K r , K NaCa , i NaK,max , g N a , g K s , g f , g t o , and g Kur ) along with the maximal uptake rate of the SERCA pump ( P u p , b a s a l ) of the FWS model were randomised by sampling a log-normal distribution of width 0.2. No heterogeneity was considered in the atrium.

To improve the generalisability of these results, five different virtual SAN tissues (#1–5) were implemented for each setup (H, UF, and HF). Each of these distributions of heterogeneity, fibroblasts, or both were generated using five different random seeds for sampling the log-normal (heterogeneity) and uniform (fibroblast) distributions.

Coupling resistance values of 1 G Ω were each chosen for the SAN cells and fibroblasts in accordance with experimental measures and previous computational works (Verheule et al., 2001; Anumonwo et al., 1992; Rook et al., 1992; Oren and Clancy, 2010; Karpaev et al., 2018). In the atrium, R gap was set to 1 M Ω to allow conduction velocities of 70 cm s⁻¹ (Fedorov et al., 2010; Ravelli et al., 2011; Kojodjojo et al., 2006). However, we doubled g N a (maximal conductance of the fast sodium current) in the Koivumäki model for the whole atrial tissue to achieve conduction velocities of 100 cm s⁻¹, which are typical of crista terminalis tissue (Ferrer et al., 2015; Amsaleg et al., 2022). In addition, a sigmoidal gradient of the coupling conductivity was implemented, similar to that by Amsaleg et al. (2022): S = 1 R R A + 1 R SAN − 1 R R A 1 1 + e − α x + β , where R R A and R SAN are the respective gap-junctional resistances of the atrial and SAN tissues, α = 1 , and β = 37 to centre the gradients in the SEPs (Figure 2). Importantly, it was necessary to spread the gradient outside the SEPs in a semicircular shape (Amsaleg et al., 2022) to achieve driving (Figure 2). None of the parameter configurations ( R gap in the SAN and atrium, α and β of the sigmoidal gradient) allowed the atrium to be excited if the change in conductivity occurred exclusively inside the SEPs. Additionally, the no-flux boundary conditions were applied.

Simulations were performed starting from single cell steady-state conditions (500 s) and lasted for 50 s. The last 5 s of the output was saved and analysed. To simulate atrial tachycardia and to investigate the responses of the SAN to overdrive suppression, a high-frequency pacing protocol was applied (Li et al., 2020). A 15 × 15 cluster of atrial cells in the top-right corner of the tissue was stimulated at 2 Hz with a current of 5 nA for 1 ms from 32 to 41.5 s before stopping the pacing. In this case, the last 10 s of the simulations output was exported to compute the sinus node recovery time ( S N R T ) and C L .

To assess the effects of autonomic control over the atrial driving capability and to test the SAN responses under challenging conditions, the effects of 25 nM of acetylcholine (ACh) and 1 µM isoproterenol (ISO) continuous infusion were assessed under both baseline and stimulated conditions. Moreover, 50% and 25% blocks in I f and I CaL of the SAN cells were tested to investigate the robustness of atrial propagation with respect to loss of key diastolic and AP ionic currents, respectively. Finally, additional SEP widths of 7, 15, and 19 cells were investigated.

For clarity, we define ‘U’ as the uniform setup (no heterogeneity or fibroblasts), ‘H’ as the setup with SAN heterogeneity, ‘UF’ as the setup with fibroblasts and without heterogeneity, and ‘HF’ as the setup including both SAN heterogeneity and fibroblasts.

The C L was computed as the time difference between two consecutive maximum points of the first derivative of the membrane voltage for the atrium or between two consecutive overshoots for the SAN, as these two methods proved to be more robust in the two cases. For the SAN, the last C L value was averaged between the cells to evaluate SAN synchronization at the end of the simulation. For the atrium, the C L values were first averaged in time and then between the cells. The safety factor ( S F ) for conduction was computed as the ratio of the net charge received by the neighbouring cells to the charge required to obtain a full upstroke in single-cell simulations, as in Boyle and Vigmond (2010): S F = C m Δ V m − Q ion − Q stim Q thr = Q gap Q thr ,

where Q gap is computed as Q gap = ∫ t a I gap d t .

Here, I gap is the net current exchanged with the neighbours (a negative current is considered to be depolarising), and t a is the time interval over which I gap is negative, as computed from the time instant at which I gap is 1% of its negative peak to when I gap changes sign. Because the atrial cells could be stimulated by the SAN cells very slowly, the linear dependence of Q thr with respect to t a was saturated at 4 ms. Under this condition, a non-negligible charge used for depolarisation would indeed fall below the 1% of I gap threshold, thus resulting in Q net being smaller than Q thr despite the cell reaching the threshold for AP firing. The S F values were computed only for the cells at the leading SEP interface, before being averaged and compared across different conditions. This was because Li et al. (2020) showed that the SAN/atrial interface provides the lowest S F . If the SAN excited the atrium across all SEPs (as in the absence of cellular heterogeneity and fibroblasts), the values were averaged between all the SEPs.

In the externally paced model tissues, the S N R T was computed as the time difference between the first atrial AP occurring after the last stimulus and the time instant of the last stimulus (Li et al., 2017, 2020). Finally, the simulated electrograms and activation time maps were generated (see the Supplementary Material for details).

Statistical differences in the C L , S F , or S N R T obtained under the different setups (H, UF, and HF, n = 5) were assessed via paired t-tests or one-sample t-tests in case of comparison with the U setup (n = 1). All data analyses were performed in MATLAB R2019b using the built-in functions, and p-values of 0.05 or lower were deemed significant.

The simulations were based on high-performance computing approaches. We developed an efficient CUDA/C program to update the state variables of the cells; the program implemented the forward Euler method to solve a set of ODEs by computing the current values of the state variables in parallel at each time step. Euler’s method is an instance of the bulk-synchronous parallel programming pattern (Valiant, 1990) comprising a sequence of parallel steps, where each parallel step updates the state variables. The bulk-synchronous pattern naturally maps to the architecture of a modern graphics processing unit (GPU) with a large number of computing cores that can operate in parallel. GPUs are particularly effective for implementing numerical algorithms over regular domains, such as vectors of matrices; indeed, we obtained a ∼ 150 × speed-up with the parallel operations compared to serial implementation.

The CUDA/C programs were written in a lightweight extension of the standard C++ programming language that included a few additional keywords and syntactic elements. A typical CUDA/C program includes ‘normal’ functions and data that are handled by the CPU as well as special functions (called kernels) that are executed by the GPU. Two different compilers were used to build the CPU and GPU executables, which were then linked in a single program.

The resulting CUDA program was implemented on a Linux server equipped with an NVidia Titan V GPU with 5,120 cores and 12 GB of device RAM, along with an AMD Ryzen Threadripper 2950x CPU operating at 3.5 GHz. A time step of 5 µs was used, and the data were undersampled at 200 µs; smaller time steps of 1 µs did not provide different results.

The behaviour of the model under the control conditions was tested first to explore the roles of heterogeneity and fibroblasts in the absence of extrinsic modifications. As shown in Figure 3, the presence of fibroblasts (UF setup) only minimally prolonged the average atrial C L compared to the uniform setup (817 ± 1 ms vs 814 ms, p < 0.01). Notably, the random fibroblast distribution showed negligible atrial σ C L (standard deviation of the C L ), while different instances of cellular heterogeneity produced widely variable C L values (817 ± 181 ms). The combination of heterogeneity and presence of fibroblasts showed a marked reduction of the C L (683 ± 32 ms) compared to those with uniform tissue and presence of fibroblasts alone (both p < 0.001); furthermore, it showed a lower C L variability than with heterogeneity alone. The large variability observed in the presence of heterogeneity (Supplementary Figure S2A) is attributed to the fact that exit blocks occur in the SEPs, with consequent alternation in the leading SEP and delay in atrial activation. The presence of fibroblasts completely abolished this phenomenon, allowing the atrial tissue to be paced at a constant steady-state rate ( σ C L = 0 , Supplementary Figure S2A). This can be ascribed to two main reasons. First, the fibroblasts help the myocytes to achieve better CL synchronisation (lower σ C L in the SAN) by blunting the intrinsic rate differences among the cells caused by cellular heterogeneity (Supplementary Figure S2A). Second, the fibroblasts act on the pacing and driving abilities of the SAN by increasing the S F with respect to the uniform tissue condition (1.62 ± 0.02 vs 1.45, p < 0.001). Figure 3 shows that heterogeneity alone increases the S F (1.58 ± 0.06 vs 1.45, p < 0.01) and that its combination with fibroblasts produces similar effects (1.60 ± 0.05, p < 0.01 vs U).

Next, we investigated the ability of the SAN to recover from external high-frequency pacing by delivering stimuli at 2 Hz to allow comparisons with prior experiments (Li et al., 2017, 2020). Excitation of the atrial tissue at this rate resulted in overdrive suppression of the SAN, with all the paced stimuli entering the SAN from all five SEPs. Once the pacing was stopped, the tissue showed C L s and σ C L similar to those under the control condition (compare the left panels of Figures 3 and 4). As in the control condition, the HF setup showed tachycardia ( C L = 588 ± 51 ms) compared to the U (815 ms), H (820 ± 218 ms), and UF (813 ± 1 ms) setups. The right panel of Figure 4 presents the S N R T values after pacing, with the heterogeneity condition showing a slightly lower value than the uniform setup (1,047 ± 18 ms vs 1,102, p < 0.01). However, in this condition, the first recovered beat originates from the SEPs, and the macro re-entry circuits with complex patterns (e.g. 2:1 exit block in the SEPs) are then stabilised in most of the tissues (Supplementary Movie S1). These observations, together with the fact that the C L s have not yet reached steady-state values after pacing, explain the high C L dispersion (Supplementary Figure S2B). Exit blocks were not observed in the UF (1,124 ± 9 ms, p < 0.01 vs U) and HF (1,007 ± 163 ms, p > 0.05 vs U) setups, but all five tissues showed re-entry activities under both cases. Thus, the fibroblasts appear to not provide faster recovery under the control conditions. Indeed, the lower average S N R T value and higher standard deviation of the HF setup are explained uniquely by the fact that tissue HF #1 had a shorter S N R T (718 ms). This occurs in all cases because the SNRT is given by the time taken by the last (nth) paced stimulus (at 41.5 s) to enter the SEPs, travel through the SAN, and exit through the SEPs (Supplementary Movie S1). However, this is not observed in HF tissue #1 (Supplementary Movie S2), where the (n-1)th stimulus excites the atrium (at 42.2 s). This is because the nth stimulus (at 41.5 s) collides with the (n-2)th stimulus (fired at 40.5 s) at the exit of SEP4. Therefore, the (n-1)th stimulus (fired at 41.00 s that has reached the midpoint between SEP3 and SEP4 in the meantime) does not collide with an incoming depolarisation in SEP4 and is free to excite the atria several hundreds of milliseconds before the nth stimulus can arrive and set the re-entry.

To study the effects of parasympathetic stimulation in the model, we simulated the administration of 25 nM of ACh to the SAN. The results in Figure 5 show profound bradycardia under the uniform (1,448 ms vs 814 ms without ACh, +78% equal to the responses of single cells), heterogeneous (1,489 ± 287 ms vs 817 ± 181 ms), and fibroblast (1,277 ± 197 ms vs 817 ± 1 ms) conditions. The combination of heterogeneity and fibroblasts provides C L s closer to the physiological resting values (991 ± 76 ms, equal to an average heart rate of 61 bpm; p < 0.05 vs H and p < 0.001 vs U). However, this is still larger than the control condition (683 ± 32 ms). Moreover, both heterogeneity and fibroblasts provide strength to atrial driving in this case ( S F = 1.60 ± 0.05, 1.76 ± 0.08, and 1.54 ± 0.09 for the H, UF, and HF cases vs 1.44 for the uniform case, with p < 0.01, p < 0.001, and p > 0.05, respectively). Here, the atrial σ C L in the H setup is still large owing to the occurrence of exit blocks and leading SEP alternation (Supplementary Figure S2C).

We next tested a challenging condition in which the SAN underwent simultaneous 25 nM of ACh infusion and external 2 Hz pacing, as in case of the ex vivo human hearts in Li et al. (2017, 2020). In this setting, the CLs showed values of 2,102 ms, 1,333 ± 152 ms, 1,169 ± 71 ms, and 952 ± 176 ms for the U, H, UF, and HF cases (p < 0.001 for all setups vs U; p < 0.05 for HF vs H). The SAN was again completely suppressed by the depolarisations travelling through all the SEPs. For the S N R T , Figure 6 interestingly indicates that the H, UF, and HF cases showed strikingly fast recoveries of pacemaking, contrary to that in the control case: 1,923 ± 325 ms for the heterogeneous tissue, 1,440 ± 269 ms with the fibroblasts, and 1,206 ± 141 ms in presence of both heterogeneity and fibroblasts vs 3,020 ms for the uniform model (p < 0.01, p < 0.001, and p < 0.001, respectively). The decrease in S N R T by the addition of fibroblasts to the tissue is larger than that given by the presence the presence of cellular heterogeneity alone. Compared to the results without ACh administration, the uniform setup showed S N R T prolongation of + 174 % ; with this setup, the SAN depolarisations originated in the middle of the SEPs as in the case without ACh, but no re-entry was established even though 2:1 exit blocks occurred for the first two beats (Supplementary Figure S4 and Supplementary Movie S3, σ C L = 950 ms). In presence of heterogeneity, re-entry circuits were formed in all tissues, but the administration of ACh increased the number of exit blocks after pacing cessation. Given the longer C L , was larger, there were clusters of spontaneous cells that fired APs before the arrival of the re-entry stimulus, later merging with it at the entrance of the leading SEP. This occurred in the presence of fibroblasts as well, but only tissue #3 of the HF setup showed exit blocks from the SEPs. Supplementary Figure S2D indeed shows very large atrial σ C L in the H setup, with a substantially lower value for the UF setup and almost no variability in the HF case. This again suggests that fibroblasts help the SEPs in capturing the atrium (Supplementary Figure S4).

After the experiments with ACh infusion, we evaluated the model responses to sympathetic stimulation. In this case, given the on/off behaviour of the FWS model with respect to adrenergic activation (calibrated to 1 µM release), it was not possible to use other concentrations, as reported in previous experiments (1 nM in Li et al. (2020)). We consequently simulated the infusion of 1 µM ISO, which shortened the C L in all cases with respect to the control: U = −21% (639 vs 814 ms), H = −30% (531 ± 13 ms vs 764 ± 183 ms), UF = −19% (663 ± 2 ms vs 817 ± 1 ms), and HF = −18% (548 ± 16 ms vs 665 ± 55 ms). The H setup did not show any exit block, with consequently reduced C L and C L variability (Figure 7 and Supplementary Figure S2E); this is attributed to the remarkable increase ( > 20%) in the S F in both the H and HF cases compared to the uniform one (1.88 ± 0.08 and 1.86 ± 0.08 vs 1.53, respectively; p < 0.001 in both cases). Although the UF setup showed a milder increase to 1.63 ± 0.07 (p > 0.05 vs U), it did not show any exit blocks.

We next added external 2 Hz pacing along with 1 µM ISO infusion. The models showed average C L values similar to those of the non-paced conditions (Figure 8). Despite the S N R T in the U setup being shorter than that without ISO (841 ms vs 1,102 ms, −24%), the H and HF cases showed much faster recoveries ( S N R T = 775 ± 18 ms and 790 ± 31 ms, p < 0.01 and p < 0.05 vs U, respectively) compared to their control values. The UF setup showed the largest S N R T (862 ± 4 ms, p < 0.001 vs U, H, and HF), which is still less than that without ISO.

To evaluate if the model was robust to loss of depolarising current during diastole and if the heterogeneity and fibroblasts contributed to its strength, the funny current I f was blocked by 50% in the SAN (Figure 9). As a consequence, a 12% C L prolongation was observed in the uniform tissue compared to the control tissue (912 ms vs 814 ms), just as in the single-cell simulations with the FWS model (Fabbri et al., 2017a). Similar results were obtained in the HF case (average + 10.7 % , 745 ± 27 ms vs 665 ± 55 ms without blocking); however, a counterintuitive tendency toward C L shortening was noted in the H setup ( − 9 % , 741 ± 25 ms vs 817 ± 181 ms). This can be explained by the presence of exit blocks in the control H setup, which induce the longer average C L in the atrium despite the fact that the CL is longer with the I f block within the SAN (741 ± 25 ms vs 662 ± 26 ms), as expected.

Both setups show lower C L (741 ± 25 ms and 745 ± 27 ms vs 912 ms, with p < 0.001, p < 0.001, and p < 0.001 for H, UF, and HF, respectively) and higher S F (1.62 ± 0.04 and 1.63 ± 0.02 vs 1.39, with all p < 0.001) values compared to the uniform SAN case. With the UF setup, only three tissues achieved atrial driving; these showed C L s similar to that of the U case (907 ± 2 ms, +11% prolongation compared to control), along with similar S F s (1.58 ± 0.13, p > 0.05).

The L-type calcium current determines the slow SAN AP. Therefore, we simulated a 25% reduction in I CaL to assess tissue robustness to a loss of the main current during the upstroke phase. This block leads to loss of driving in the U and in all H models. Intriguingly, three out of the five UF models and four out of the five HF models showed pace-and-drive behaviours. The UF tissue #1 as well as HF tissues #3 (Supplementary Movie S4) and #4 showed 2:1, 3:1, and 2:1 exit blocks, respectively, with consequent extreme atrial bradycardia ( C L = 2,209 ms or 27 bpm with only two APs in 5 s, C L = 2,583 ms or 23 bpm, and C L = 1,797 ms or 33 bpm). The HF tissue #5 also showed complex electrical activity, with re-entry excitation showing one exit block in the last 5 s of the simulation and a resulting C L of 1,142 ± 463 ms. Only HF tissue #2 showed a regular physiological rate ( C L = 806 ± 0 ms, 74 bpm), but its activation started from the centre of #SEP2 (Supplementary Figure S5). The UF tissues #2 and #4 show stable driving but at a slightly bradycardic rate ( C L = 1,073 ± 0 ms and 1,076 ± 0 ms, both equal to 56 bpm). In all these cases, the SF appeared to be slightly reduced compared to that of the HF model under other conditions (Supplementary Figure S3).

Supplementary Figure S3 summarises the S F results for the four setups under the different conditions. In the UF condition (panel C), there was frequency dependence of the S F , which increased with longer C L s (e.g. with ACh or I f block). In addition, SAN frequency synchronisation was worse in the HF case than the H and UF cases with both I f and I CaL current blocks (Supplementary Figure S2, panels G and H). In all other cases, the combination of heterogeneity and fibroblasts produced lower or equal σ C L values within the SAN.

To gain further understanding of the driving mechanisms, two additional model parameters were investigated, namely the width of the exit pathways and cellular organisations at their interfaces. Table 1 reports the effects of the SEP width on the driving capability of the SAN, showing that there is a lower limit (width = 11 cells) for propagation. The results are similar for both the K and MBS models, with no modulations of the C L values with respect to the SEP width. However, the conduction velocity along the leading SEP increases proportionally with the width and a corresponding decrease in the conduction time. The MBS model required a mosaic configuration in addition to the design described previously to be able to pace the atrium; this included a gradient of atrial cells interspersed in the SEPs, without which it was not possible to obtain propagation in the uniform SAN condition. When considering heterogeneity (without mosaic), only one of the five models achieved very poor driving, with a 3:1 exit block and consequent pronounced bradycardia ( C L = 1,932 ms). However, even with the mosaic (uniform) configuration, the atrium was effectively activated in only three of the five models ( C L = 814 ± 0 ms).

A final set of simulations including heterogeneity, fibroblasts (HF tissue #3 with K model), and simultaneous administration of 25 nM of ACh as well as 1 µM of ISO (hereafter called the ‘full’ model), was implemented to model a condition as close to the SAN physiology as possible. To reproduce the physiological activation sequence, R gap was lowered to 10 and 100 M Ω for the SAN cells and fibroblasts, respectively. In addition, the mosaic model was implemented according to the MBS model simulations, while the FWS model was modified to account for the simultaneous adrenergic and cholinergic effects in an additive manner (Supplementary Material). This setting provides a C L of 1,192 ms, with the activation starting from the left central part of the SAN and exciting the atrium through SEP3 after roughly 80 ms (Figure 10 and Supplementary Movie S5). The propagation wave then re-enters the SAN from SEP1 and SEP5, where it collides with the depolarisations arriving from the head and tail of the SAN. The S F in this condition was 1.56.

When the R gap of the fibroblasts was retained at its original value of 1 G Ω , a more physiological C L of 961 ms and higher S F of 1.75 were achieved; however, atrial excitation occurred after ∼ 120 ms through SEP#5, and the slow CV inside the SAN allowed entry of the depolarisations from the other SEPs (Supplementary Movie S6). Interestingly, pacemaking activity was completely suppressed in the absence of the fibroblasts.

The important and novel roles of cellular heterogeneity and fibroblasts displayed by the proposed model include increasing the ability of the human SAN to drive the atrium (Figures 3, 5, 7, and 9). In all conditions, the H configuration showed substantially higher S F values ( > + 10 %); this is likely due to the presence of stronger cells (i.e. expressing, for example, more I CaL ) that deliver more current at the SEP interfaces. The UF and HF setups showed similar behaviours under most conditions since the fibroblasts acted as current sources and allowed the SAN cells to provide more charge to the atrial tissue, as described below. The information in Supplementary Table S2 supports this observation by reporting how the heterogeneity and fibroblasts determine more depolarised maximum diastolic potentials ( M D P s) as well as higher upstroke velocities, especially in the leading SEP, compared to when they are absent. The contributions of cellular heterogeneity to the ionic properties have long been deemed central to robust pacemaking (Boyett et al., 2000). With regard to fibroblasts, experimental and computational studies on the ventricular tissues (Fahrenbach et al., 2007; Xie et al., 2009) have allowed hypothesising their ability to influence the source–sink relationship for effective stimulus propagation in the SAN (Csepe et al., 2015). Indeed, a possible first validation of the present results can be obtained from co-cultures of fibroblasts and ventricular myocytes (Miragoli et al., 2007): in these experiments, the percentage of myocyte preparations showing spontaneous activity was above zero only when a minimum amount of fibroblasts (20%) was present, and this percentage was maximum ( > 80%) with 50% fibroblast density. Considering that fibroblast–myocyte coupling has been demonstrated in vivo recently in mouse infarcted ventricular tissue (Wang et al., 2023), our simulations point toward the role of the depolarised resting potential of the fibroblasts in robust pacemaking.

The administration of 1 µM of ISO showed synergistic effects with both the H and HF configurations, producing greater increases in the S F than that for the uniform case (Supplementary Figure S3). In fact, if we exclude the S F increase in the UF setup in the presence of ACh (likely due to the myocytes undergoing a longer and hyperpolarised diastolic phase, that allows fibroblasts to supply them more charge, only ISO provided higher S F values than the other conditions (control, ACh, and current blocks) for the same configurations (Supplementary Figure S3).

The full model showed an S F value (1.56) closer to the average obtained with ACh administration (1.54 ± 0.09) than that with ISO infusion (1.86 ± 0.08). When considering a weaker coupling of the fibroblasts (original R gap value of 1 G Ω ), this is sufficient to substantially increase the S F toward the ISO value (1.75); this is consistent with the results obtained from coupled rabbit SAN and atrial cell models, which showed a substantial dependence of the SF on R gap (Cacciani and Zaniboni, 2015). An additional beneficial consequence of the increased R gap in the full model is the avoidance of bradycardia (62 vs 50 bpm). Nevertheless, both effects are obtained at the cost of losing phase synchronisation inside the SAN (lower apparent conduction velocity), a fact that may facilitate re-entry onset (Supplementary Movie S6). Most interestingly, when the fibroblasts were removed from the full model, the SAN tissue lost automaticity. This supports the importance of fibroblasts in mitigating atrial hyperpolarisation, thereby protecting depolarisation of the SAN myocytes and allowing atrial driving.

The simulations with I f and I CaL blocks showed that heterogeneity protected the SAN from loss of key ionic currents when coupled to the atrium (Figure 9), confirming its beneficial actions commonly found in isolated pacemaker tissues (Campana et al., 2022; Maltsev et al., 2022). For fibroblasts, stronger beneficial effects were observed, allowing the SAN to pace the atrium despite a 25% I CaL block with the UF and HF setups, while the U and H tissues remained quiescent. Similar to the increase in S F , this is attributed to the depolarisations of the MDPs and especially the increase in the d V / d t max of the SAN cells in the leading SEP owing to the synergism of heterogeneity and fibroblasts (Supplementary Table S3).

In this work, the indirect S N R T , obtained as the time interval between the first atrial AP and last paced beat (Li et al., 2017, 2020), was used to quantify the ability of the SAN to recover its driving capability after overdrive suppression. In the heterogeneous tissue under the control conditions (no blocks or drug administration), the first spontaneous beat propagates from the middle of the SEP slightly before than it does in the uniform tissue. After this beat, re-entry is established in all five models, but they do not stabilise on specific SEPs as in the case with fibroblasts. The latter configuration does not have exit blocks, in line with the observations of Zhao et al. (2023), who observed lower exit block occurrences upon sodium current block or ACh administration compared to when the fibroblasts were not considered in remodelling owing to heart failure. This suggests that even though fibroblasts provide strength to atrial driving, they also lead to macro re-entry stabilisation in the present model.

After pacing at 2 Hz, all setups showed SNRT values close to experimental observations in two explanted human hearts: U = 1,102 ms, H = 1,047 ± 18 ms, UF = 1,124 ± 9 ms, and HF = 1,063 ± 163 ms vs the values of 1,020 and 1,459 ms reported by Li et al. (2020). When corrected for the basal C L , the S N R T of the uniform setup (288 ms) was in agreement with the extant computational reports (167 ms with pacing at 950 ms; Li et al. (2020)). However, after external pacing, no exit blocks were observed in the experiments of Li et al. (2020) and in the computational model reported by Zhao et al. (2023). In both these cases, the exit blocks were obtained only with reduced sodium current and ACh administration or under heart failure remodelling conditions, which also led to re-entry.

Furthermore, simulations with administration of 25 nM of ACh showed that the fibroblasts greatly reduced the S N R T compared to the heterogeneous tissues (−23% for UF and −29% for HF), which in turn allowed earlier recovery than in the uniform case (−36%, Figure 6). In particular, the S N R T values in the UF and HF setups with ACh administration were much closer to the values without ACh (1,479 ± 223 ms vs 1,124 ± 9 ms and 1,371 ± 165 ms vs 1,007 ± 163 ms) than that with the H setup (1,923 ± 325 ms vs 1,047 ± 18 ms). Our finding of SNRT = 3,020 ms with the U setup, which shows nodal dysfunction (corrected indirect S N R T > 525 ms), is in line with the ex vivo experimental findings on dogs (Lou et al., 2013, 2014) and humans (Li et al., 2020). Modelling work on humans (Li et al., 2020) congruently showed atrial pauses of > 5 s due to exit blocks at the SEP exits upon ACh application when the fast sodium current was blocked by 20%. In the rabbit model by Karpaev et al. (2018), the durations of the pauses were proportional to the number of fibroblasts. This discrepancy may be ascribed to the different basal heart rate. No pauses were reported in the model by Zhao et al. (2023). During pacing, none of the models showed filtering of the paced stimulus through the SEPs, as previously observed in optical mapping experiments in human hearts at frequencies lower than or equal to 2 Hz (Li et al., 2017).

In the uniform setup, the models showed C L s equal to the single-cell values (814 ms); this is shorter than values reported in earlier computational works (vs 883 ± 1.7 ms in Amsaleg et al. (2022), 999 ms in Li et al. (2020), and 930 ms in Zhao et al. (2023)) but still within the range of human heart rates (74 vs 60 − 100 bpm).

The occurrences of exit blocks in the baseline H setups explain the large atrial σ C L and average values close to those of the uniform setup. Interestingly, the simultaneous presence of heterogeneity and fibroblasts resulted in lower C L and σ C L . Unlike the observations in the preliminary human single-cell simulations reported by our group (Ricci et al., 2021), when the SAN cells were hyperpolarised by the atrium, the duration for which their membrane voltages were below the resting potential of the fibroblasts (−50 mV) increased. Thus, the fibroblasts acted as current sources over a larger portion of the diastolic depolarisation ( D D ) phase, accelerating it and strengthening atrial depolarisation. In addition, owing to their depolarised resting potentials, the fibroblasts softened the intrinsic rate differences due to parameter randomisation among the cells by pushing them toward similar beating frequencies (Supplementary Figure S2) and provided enough current to avoid exit blocks from the SEPs, further supporting the hypothesis of protective actions of the fibroblasts from hyperpolarisation.

By itself, the presence of fibroblasts prolongs the C L minimally (Figure 3); this small effect is consistent with the findings of Karpaev et al. (2018) but contrasts in amplitude with the larger prolongation found by Oren and Clancy (2010), with the experiments of Csepe et al. (2015) and Fahrenbach et al. (2007), as well as with our previous results in isolated human SAN models (Ricci et al., 2021). However, this discrepancy may be explained by the voltage difference between the resting potential of the fibroblasts and M D P of the SAN cells when they are also coupled to the atrium, according to the mechanism described above.

Interestingly, when ACh is applied, the combination of heterogeneity and fibroblasts prevents the occurrence of exit blocks in the H model and shortens the C L in the U setup (1,448 ms, +78%, slightly longer than that reported by Li et al. (2020) but in line with the findings of Zhao et al. (2023) who obtained +28% prolongation with 10 nM of ACh), thereby avoiding bradycardia (991 ± 76 ms, equal to 61 bpm) and reducing C L variability (Figure 5). During ISO administration, the HF setup also showed a tendency to limit the tachycardia (109 vs 113 bpm) achieved in the H setup (Figure 7), in agreement with the lower average rates obtained with the fibroblasts alone (90 bpm). Therefore, from this perspective, the fibroblasts provide robustness to atrial driving but reduce heart rate flexibility at the same time.

In the full model, a heart rate of 50 bpm was obtained; this slight bradycardia was the consequence of the reduced R gap of both the SAN cells and fibroblasts needed to achieve a physiological atrial activation time. Furthermore, when considering the simultaneous simulated infusion of ACh and ISO, the model shows the predominance of cholinergic over adrenergic effects, in line with the experimental observation that the intrinsic rate (e.g. in the denervated heart) is higher than the in vivo rate owing to the absence of the large basal parasympathetic tone (Opthof, 2000; Peters et al., 2020).

Depending on both the randomised ionic properties and the fibroblast distribution, clusters with different intrinsic properties (e.g. beating rate) are formed inside the tissue and become dominant. The activation sequence originates according to the locations of these clusters (Supplementary Figure S5 and Supplementary Movie S6) and more distally from the atrium compared to the uniform model.

Thus, the physiological levels of heterogeneity ( σ = 0.2 for the log-normal distribution (Sobie, 2009; Campana et al., 2022; Maltsev et al., (2022)) and fibroblasts (Csepe et al., 2016, 2017; Li et al., 2023) can overcome the atrial sink effect to shift the leading pacemaker (LPM) location at least at the coupling values adopted. Interventions that decrease the strength of the source or increase the sink, such as reducing R gap in the full model, administering 25 nM of ACh, or blocking I f , shift the LPM away from the SEP exits. Consequently, the electrical excitation follows a more regular pattern of activation, propagating from the left part of the SAN toward the atrium (Figure 10).

Interestingly, maintaining a fixed distribution of the diffused fibroblasts while changing the heterogeneity distribution always resulted in atrial activation from the same SEP (not shown), suggesting that the fibroblasts may guide propagation. This is further strengthened by symmetrical simulations in which the distribution of cellular heterogeneity was fixed and that of the fibroblasts was changed, resulting in atrial activation via different SEPs.

A steep gradient in the gap-junctional coupling was found to be necessary to achieve atrial driving. This suggests that the most important parameter for achieving robust atrial excitation is intercellular coupling, as highlighted in previous works (Kharche et al., 2017; Inada et al., 2014; Amsaleg et al., 2022). Indeed, the slope, half value, and transition shapes inside and outside the SEPs were seen to deeply affect propagation, with only a fine balance among these parameters allowing propagation in the model, as reported in the Methods.

The number of SEPs does not influence the average C L or the S F , suggesting that the designed coupling gradient allows the SAN to be a sufficiently strong source that is unaffected by increased sink (Fabbri et al., 2017b). With regard to the SEP widths, the simulations show that there are lower limits on the numbers of cells comprising the exit pathways in both the K (Koivumäki et al., 2011) and MBS (Mazhar et al., 2023) models for excitation of the crista terminalis, as found in the detailed 3D model by Amsaleg et al. (2022). However, Li et al. (2020) managed to achieve SAN–atrial propagation with one SEP as narrow as three cells, even in the context of simulated heart failure (20% I f block; I N a blocks in the SAN, non-conductive fibrosis in the SAN and SEP; 5% I N a block in the atrium). This difference could possibly be explained by the fact that the authors scaled the maximal conductance of the sodium current in the SEP by a factor of five, thereby increasing the ability of the pathway to capture the atrium.

As the SEP widths increased, slight increases in the conduction velocity (and consequent reductions in the conduction times) were obtained in both the K and MBS models (Table 1). However, the C L did not appear to be affected by the SEP width, contrary to the findings of the 3D model by Amsaleg et al. (2022). This indeed showed that enlarging the SEPs resulted in loss of driving due to excessive hyperpolarisation arising from a larger sink. This observation was in agreement with a previous study that reported a C L prolongation up to SAN suppression with wider SEPs, while narrower SEPs allowed pacing and driving over a larger range of conductivities (Zyantekorov et al., 2018). This discrepancy may be ascribed to different aspects. First, in a 2D setting like the one considered in this work, the interface between the SAN and crista terminalis grows linearly with increasing SEP width but quadratically in 3D. Second, Zyantekorov et al. (2018) considered SEPs that were open on both sides of the SAN. Third, the K model has a membrane capacitance of 50 pF, which is less than the those used in Amsaleg et al. (2022) (100 pF) and Zyantekorov et al. (2018) (81 pF). All these factors could lower the sink effect in the present model, limiting C L modulation by the SEP width.

In the last set of simulations, the MBS model was less easily excitable than the K model; the reason for this is that the present work only involved the 29 equations describing the APs of the atrial cells, while neglecting the ones devoted to modelling the mechanoelectric feedback and calmodulin-kinase activity. Although the AP features were still in accordance with experimental values (Supplementary Table S1), this set of equations was not optimised to describe the electrophysiology of single atrial cells. For example, the maximal upstroke velocity was lower than in the original MBS model ( d V / d t max = 177 vs 199 mV ms⁻¹). Thus, even when the maximal conductance of the fast sodium current was doubled to reproduce the experimental values of the conduction velocity in the crista terminalis, the atrial tissue comprising the MBS model cells struggled to depolarise.

Although further refinements are needed to use the electrophysiological equations of the MBS model in a standalone manner, these simulations provide some useful information on atrial driving. First, atrial cells interspersed in the exit pathways (so-called mosaic model) actually help propagation compared to a flat interface. Although this setup was not needed with the K model, its addition to the gradients in the gap-junctional coupling could represent a more robust and hence more likely anatomical architecture of the SAN (Karpaev et al., 2018; Amsaleg et al., 2022). Therefore, it was included in the full model. Second if the SAN heterogeneity increases the amount of charge delivered to the atrial cells (as discussed above), this is insufficient to overcome the gaps that can be addressed only by coupling or anatomical aspects (such as the gradient and mosaic configurations). Indeed, without the addition of the mosaic setup, only one SAN heterogeneous distribution was shown to drive atrial tissues comprising the MBS model (and to an unphysiologically low rate owing to exit blocks).

The present model has several shortcomings. First, to the best of our knowledge, no SAN-specific fibroblast model exists for humans or any other species. In this work, we chose to use the model by Morgan et al. (2016) as it is a more similar phenotype (atrial) to other human fibroblast models (Sachse et al., 2008; Maleckar et al., 2009; Liu et al., 2021).

Second, the fibroblasts were distributed diffusely according to a concentration gradient, with higher amounts distally from the crista terminalis and lower amounts when approaching the ends of the SEPs (Li et al., 2023). However, histological sections of the SANs of different species of mammals showed that the fibroblasts actually formed structural networks supporting and protecting different clusters of SAN myocytes (Camelliti et al., 2004; Sánchez-Quintana et al., 2002; Smithers et al., 2021).

Third, the 2D geometry adopted herein represents a rough simplification of the complex structure of the SAN (Chandler et al., 2011; Fedorov et al., 2010; Csepe et al., 2016) and is simpler than other recent models that also consider fibre direction (Amsaleg et al., 2022; Zhao et al., 2023). Considering that the anatomy has a clear impact on function in an elaborate structure like the SAN (Nikolaidou et al., 2012), our stylised model geometry affects the results quantitatively. For example, conically shaped SEPs may increase the driving strength and limit re-entry occurrence by providing less favourable source–sink balance for the excitations entering the SAN. Moreover, we did not consider the presence of spatially dependent distributions of ACh receptors (Lang et al., 2011) or localisation of different intrinsic pacemaking properties in the specialised SAN clusters (Brennan et al., 2020; Maltsev et al., 2023), which were aspects proposed to explain the shifts in the leading pacemaker location following autonomic modulation.

However, the main drawback of the proposed model is the slow conduction within the exit pathways. Considering a SAN CL of 67 μm, the average conduction velocity in the SAN was ∼ 1 cm s⁻¹ (lower than the experimental value of 5 cm s⁻¹ reported by Fedorov et al. (2010)), which dropped to < 0.5 cm s⁻¹ in the SEPs. These velocities correspond to conduction times exceeding 300 ms in the SEPs (Table 1). Therefore, on average, atrial activation was achieved more than 500 ms after the stimulus originated in the SAN. this value is much higher than in the other models ( < 100 ms reported by Amsaleg et al. (2022) and Li et al. (2020)) as well as outside the experimental range (82 ± 17 ms as reported by Fedorov et al. (2010)). The slow sequence of activation sets the stage for re-entrant activity, as in the H configuration under the control conditions and in the H and HF cases after external pacing (Supplementary Movies S1, S2). Amsaleg et al. (2022) reported velocities as low as 1 cm s⁻¹ at the beginning of the gap-junctional gradient in the SEP, but these sharply increased along the SEP. When we used the mosaic architecture with the MBS model, the limited width of the 2D SEP did not allow enough cells for smooth transition to an atrial cell phenotype; thus, the CV remained very low for most of the SEP length before reaching 100 cm s⁻¹ just before exiting the SEP.

Finally, as a first approach, all the conditions tested in this work were only with regard to the SAN. This is reasonable considering the presence of fibroblasts and vagal tone, which are known to be minimal in the right atrium under physiological conditions. However, this approach neglects the heterogeneity in the ionic properties of the atrial tissue and its responses to sympathetic stimulation, which likely play roles in modulating the ability of the SAN to drive the atrium.