A schematic view of the one-dimensional rabbit cardiac conduction system model incorporating AVN with dual pathways is shown in Figure 2. The model consists of a two-row matrix, where the upper row includes the SN and peripheral sinus node (PS1–PS3) cells, atrial muscle (AM1–AM3), fast pathway (FP1–FP8), penetrating bundle (PB), and HB (HB1-HB6) cells, while the lower row includes only ten slow pathway (SP1–SP10) cells and additional intermediate atrial muscle cell (AM*). The latter is added to provide a smooth transformation of short rectangular-shaped excitation impulses into typical atrial action potential shapes in the case of AFB (with random pacing) and AFL (regular pacing) simulations. The arrows Antero and Retro denote the points of application of external S1 and S2 stimuli for anterograde and retrograde propagation study with S1S2 stimulation protocol, and the arrow AF indicates entrance for “internal” stimuli in the cases of AFB and AFL simulations. The premature S1S2 stimulation (Katritsis and Morady, 2022) refers to pacing that is performed at a constant basic cycle length (shorter than normal sinus rhythm) for 8–10 S1 impulses, followed by the introduction of an extra stimulus S2, and with stepwise reduction of the S1S2 interval until AVN refractoriness (total conduction block) is reached. Therefore, S1 and S2 mark the basic and test stimulus, respectively, and the same subscripts apply to atrial (A) and HB (H) stimulation or response.

Each model cell represents a group of real cardiac cells and is described by two-variable A-P model. The model for monodomain description of the cardiac cell consists of the following set of reaction-diffusion type non-linear ordinary differential equations: ∂ u ∂ t = c t k u u − a 1 1 − u − v u + I coupl + I stim , (1) ∂ v ∂ t = c t ε − v − k u u − a 2 − 1 , (2) ε = ε 0 + v μ 1 u + μ 2 , Where 0 ≤ u ≤ 1 and v are the dimensionless normalized transmembrane potential and slow recovery variables, respectively. Parameter k controls the magnitude of transmembrane current, and parameters μ ₁, μ ₂, a ₁, and a ₂ are adjusted to reproduce characteristics of cardiac tissue, ɛ sets the time scale of the recovery process determining the restitution properties of the action potential, c _t is the time scaling coefficient introducing physical time in seconds into the system. In the general case, a ₁ = a ₂ > 0. Realistic transmembrane voltage values of the AVN cells in millivolts can be obtained with the rescaling of the u value by 120u − 80 mV for the AM, and 100u − 70 mV for HB (Inada et al., 2009). Since the A-P model was originally developed for the description of canine ventricular myocyte, the parameters were adjusted for the rabbit AVN in terms of action potential duration and refractoriness (see Supplementary Tables S1, S2). I ^stim is the external stimulation current, I ^coupl = ∇ ⋅ (D∇u) is the coupling (diffusion) current, where ∇ is a spatial gradient operator defined within the model tissue geometry, and D is a tensor of diffusion coefficients characterizing electrotonic interactions between neighboring cells via gap junctional coupling conductance. For the considered one-dimensional system, the coupling term in the discrete form reads as follows I i coupl = d i − 1 u i − 1 − u i + d i u i + 1 − u i , (3) where i is the model cell index, and d is the coupling (diffusion) coefficient in s⁻¹ assigned for each pair of cells normalized on dimensionless distance. As will be discussed below, to enhance the direction-dependent conduction properties of the AVN model, we implemented a coupling asymmetry with the modified coupling term: I i coupl = d i − 1 u i − 1 − α i − 1 u i + d i − u i + α i u i + 1 , (4) where α > 0 is the coefficient of asymmetry applied to the right-hand side cell in each coupled pair of model cells. The values α = 1 correspond to symmetric coupling, α < 1 to the accelerated anterograde and slowed retrograde conduction, and α > 1 to the accelerated retrograde and slowed anterograde conduction.

The cells playing the role of natural primary (SN) and subsidiary (SP7 through SP10, FP8, and PB) pacemakers are gray-shaded on the scheme (Figure 2). They are described by the same A-P model with a negative coefficient a ₁. The oscillating (pacemaking) properties of the A-P model were considered recently (Ryzhii and Ryzhii, 2022). Implementing the oscillating A-P model cells simplifies the AVN model description, providing easy control of intrinsic oscillation frequency and reliable overdrive suppression. That substitutes a significant advantage of the cardiac A-P model over the FitzHugh-Nagumo and modified Van der Pol models and their combinations, used in the modeling of the simplified cardiac conduction system (Podziemski and Żebrowski, 2013; Ryzhii and Ryzhii, 2014), due to the absence of the quiescent state in the modified Van der Pol model and the absence of well-defined threshold potential in both FitzHugh-Nagumo and modified Van der Pol models. Also, close values of the minimal diastolic potential of both non-pacemaking and pacemaking A-P model cells allow the driving of a relatively big number of non-pacemaking cells by a single pacemaking A-P model cell (Ryzhii and Ryzhii, 2022).

As a starting point for the search of optimal model cell parameters, we took the experimental data from (Reid et al., 2003; Billette and Tadros, 2019): minimal and maximal AVN conduction time for normal case and after the FP and SP ablations for anterograde conduction and retrograde conduction, the nodal effective refractory period (ERPN) for control case and after the SP ablation, and the exponental-like shapes of conduction curves for both anterograde and retrograde directions.

Taking into account prominent difference in antrograde and retrograde conduction curves observed experimentally, we assumed that the characteristics of FP and SP should include non-linear conduction properties (some degree of rectification, e.g., diode-like conduction) due to spatial gradients of refractory period, local intercellular coupling, coefficient of coupling asymmetry, and action potential duration. For the sake of simplicity, in the course of the model development, we paid attention to the differences between the FP and SP local conduction delays due to spatial distribution of coupling coefficient but not the macroscopic conduction velocities. Faster conduction velocity in the FP is not confirmed in experiments (Hucker et al., 2008). Anatomically shorter FP distance compared with the SP is considered a primary cause of a shorter conduction delay in the FP (Hwang et al., 2014). Thus, though in our phenomenological model, the SP is longer than the FP in terms of the number of model cells (10 vs. 8), this difference in the model cell numbers may not correspond to the real anatomical structure of the rabbit AVN. Proper pathway conduction delays were obtained by setting spatial distributions of the coupling coefficient d and coupling asymmetry α values, and refractory periods of the model cells (see Figure 3 and Supplementary Tables S1, S2). In this regard, in our model, the cells represent groups of real cells and are distinguished not by biophysical types (de Carvalho et al., 1959; Billette, 1987; Inada et al., 2009) but by a spatial gradual distribution of the parameters. Figure 3 shows the distribution of refractoriness in the uncoupled model cells, local coupling coefficient d, and coefficient of coupling asymmetry α. The refractory period values are in line with the available experimental data: 166 ± 30 ms for the rabbit SN and 68 ± 11 ms for atrial muscle (Kerr et al., 1980), 91 ± 10 ms for the inferior nodal extension (Reid et al., 2003), and 141 ± 15 ms for the transitional zone (Reid et al., 2003), which is an anatomical substrate of the FP. The profile of d approximately corresponds to the distribution of coupling conductance in the ion-channel multicellular model (Inada et al., 2009). Action potential profiles of the uncoupled model cells in the FP and SP can be found in Supplementary Figure S1.

The FP and SP ablations were simulated by setting diffusion coefficient d = 0 in the middle of the pathways (between FP5 and FP6, and SP6 and SP7 cells).

The S1S2 stimulation protocol consists of ten basic S1 impulses followed by a premature S2 impulse. We created the AVN conduction curves applying constant S1S1 basic cycles with 360 ms intervals (shorter than spontaneous sinus rhythm of 361 ms) followed by the test premature S1S2 cycle with stepwise reduction starting from the basic cycle length until the full AVN functional refractory block occurs (360–90 ms) (Billette and Tadros, 2019). The current stimulation impulses were of 1 ms duration and 1.3× greater than the threshold value (I ^stim = 280.0). The AFB and AFL were simulated by applying similar stimulation impulses through the intermediate cell AM* to eliminate the influence of the pacing impulse shape. The correct shape of AM* action potentials allowed more natural stimulation of the atrial part of the conduction system model. For the simulation of the AFB, a random sequence of the impulses with the same amplitude and duration was generated using MATLAB pseudorandom number generator with uniform distribution in the 75–125 ms range.

An original algorithm was developed to visualize the laddergrams with the exact activation timing of each model cell. It allows tracing the wavefront propagation for each excitation sequence and in each pathway separately (with red color for the FP, and blue color for the SP) based on the activation time in each cell. The activation times in each model cell were recorded at the level of 0.25 during the action potential upstroke.

For anterograde propagation, the algorithm determines the leading pathway which provides the excitation to reach the PB first, thus, the traces coming down to the HB are marked in the same color. This is very useful since it allows visual differentiation between conduction through the SP and FP obtained in the simulations and easy comparison with the experimental method (Zhang et al., 2003; Climent et al., 2011), in which pronounced change in the signal amplitude in the inferior HB domain occurs when conduction switches from one pathway to another, and is especially important for visualization of chaotic rhythms (Zhang et al., 2003).

The solution of the set of Eqs. (1) and (2) was performed with MATLAB (R2022a) standard ODE solver ode23 with relative tolerance 10^–7 and absolute tolerance 10^–10. The ode23 is a single-step solver and implements the explicit Runge-Kutta (2,3) pair (Shampine and Reichelt, 1997). Since ode23 uses an adaptive time stepping, for the reliable periodic and random pacing, MATLAB mod function was used to force the solver to accept input impulses at predefined time moments. The solver output with the variable time steps was evaluated with the finite discretization of 0.1 ms for further analysis. Open-end boundary conditions were used at both SN and HB6 ends of the model.

The MATLAB source code for this study is available in Figshare repository (https://doi.org/10.6084/m9.figshare.22015289.v1).

Figures 4–8 demonstrate different cases of simulated excitation propagation in the AVN. Each figure consists of a scheme and three panels. The scheme illustrates the propagation conduction sequence of interacting pathways in each case. In Figures 5–8 the scheme corresponds to the response on S2 stimulus. The top panel shows the laddergram with exact timing of each model cell excitation. The trace color in the HB part corresponds to the leading pathway. The middle and bottom panels demonstrate the sequences of action potentials passing between the SN and HB through the FP and SP, respectively.

First, normal sinus rhythm with a cycle length of 361 ms (166 bpm) which is typical for rabbit heart (363 ± 21 ms, (Garrigue et al., 1999)) was simulated, and the results are shown in Figure 4A. The numbers in the top panel indicate excitation wave latency with respect to the SN and the sinus rhythm in milliseconds. The excitation wave initiated at the SN propagates along the peripheral sinus node and AM cells and reaches the last atrium point (AM3) at 43.7 ms. Then the wave propagates through both the FP and SP. The wave in the FP arrives first at the PB cell at 92.4 ms, causing retrograde propagation through the SP. Both SP anterograde and retrograde excitation fronts meet in the inferior nodal extension at 99.8 ms and cancel each other. The leading FP conduction wave travels further down from the PB to the last HB cell of the model and arrives there at 139.6 ms. The obtained latency values are very close to the experimental measurements (de Carvalho et al., 1959) shown in Figure 1B for the right atrium of the rabbit heart and computer simulations with complex 1D ion-channel based rabbit AVN model (Inada et al., 2009). A similar scheme of AVN conduction via the FP and SP was also observed in the case of atrial pacing at long S1S2 intervals (see Figure 5A below).

In addition to its essential role in controlling the electrical impulse transmission between the atria and ventricles, the AVN can also serve as a subsidiary pacemaker when the SN fails (Hwang et al., 2014). In rabbits, the AVN leading pacemaker site is defined in the inferior nodal extension, generating electrical excitation with an average cycle length of about 600 ms (627 ± 111 ms, (Dobrzynski et al., 2003)). In this case, the excitation spreads in two directions - upward to the atria and downward to the HB. Figure 4B shows the laddergram of AVN automaticity and the corresponding propagating action potentials in the absence of SN activity. The SP is the leading pathway for both anterograde and retrograde directions.

The most important properties of the AVN behavior are reflected in the so-called AV nodal conduction or recovery curve (Talajic et al., 1991). The characteristic steps of the AVN anterograde conduction curve calculation during atrial pacing are presented in Figures 5, 6. In the range of long S1S2 intervals, the conduction scheme is similar to that demonstrated in Figure 4A with the FP leading.

At the S1S2 interval shortened to approximately 140 ms, a transitional stage occurs when the FP and SP excitation waves reach the PB almost simultaneously with minimal difference in latencies (Figure 5A). In such a case, while our model determines the leading pathway mathematically, from a physiological point of view both pathways provide conduction down to the HB. Simultaneous input of both pathways was experimentally demonstrated in HB electrograms with intermediate height (Zhang et al., 2003).

For the S1S2 intervals shorter than 135 ms, only the SP becomes functional and allows passage of the excitation from the SN to the HB (Figure 5B).

By reducing the interval below 123 ms, we observed the onset of perpetual AVN reentrant tachycardia (AVNRT) (Nikolski et al., 2003; Katritsis and Josephson, 2013), when excitation repeatedly circulates within the AVN ring (Figure 6A). The existence of AVNRT in the rabbit AVN at short S1S2 intervals was reported in previous experiments (Lin et al., 2001). An excitation wave propagating anterogradely along the SP, retrogradely travels back up the FP, and then activates the atrial muscle, after which subsequent reentry cycles occur. This phenomenon takes place because the FP has a significantly longer refractory period than the SP (Figure 3). When excitation travels via the SP, the upper part of the FP is blocked, but the lower FP part recovers its excitability and conducts in a retrograde direction (Billette and Tadros, 2019). According to the AVNRT classification (Katritsis and Josephson, 2013), the observed AVNRT belongs to the typical slow-fast type. AVNRT persisted with further shortening of S1S2 intervals until full functional block was achieved. In our simulations, the full functional block occurred at S1S2 intervals shorter than 95 ms when there was no conduction through both pathways (Figure 6B).

In all the above cases, the excitation wave from the SN pacemaker meets and annihilates with the retrograde wave from the atrial muscle caused by premature S1 impulses, and the SN activity is suppressed by S2 impulses.

Animation of S1S2 protocol simulation for anterograde atrial pacing is shown in Supplementary Video S1.

Catheter ablation of the heart is used to destroy small areas in one of the AVN pathways, thereby stopping the propagation of unusual or undesirable electrical signals that travel through the AVN and cause irregular arrhythmias (Hong et al., 2020). Figure 7 shows the simulation results of the atrial pacing after the SP and FP ablations. In the former case (Figure 7A), when S1S2 impulses are delivered at 110 ms intervals, the conduction of excitation is totally canceled since, at such short pacing intervals, functional conduction block in the FP takes place. Therefore, due to the broken AVN ring, SP ablation eliminates AVNRT shown in Figure 6A (Katritsis and Josephson, 2013).

FP ablation (Figure 7B) leads to increased conduction time (Billette and Tadros, 2019). Although the response to the S2 impulse looks similar to that shown in Figure 5B, where FP conduction is absent due to the functional conduction block associated with FP refractoriness, the increase in the conduction time is due to a different conduction propagation scenario caused by the preceding S1 impulse. Because ablation destroys pathway tissue, this conduction scenario extends over the entire range of S1S2 intervals (see complete anterograde conduction curves in Figure 9).

Retrograde conduction refers to conduction from the lower part of the cardiac conduction system to the atria through the AVN. Experimental and clinical studies have demonstrated the ability of the AV node to conduct impulses in the anterograde and retrograde directions (Narula, 1974; Mitsuoka et al., 1987; O’Hara et al., 1993; Billette and Tadros, 2019).

We simulated retrograde conduction applying S1S2 (H1H2) stimulation to the HB (at HB1 model cell, Retro arrow in Figure 2) and plotted pertinent recovery curves (see AVN conduction curves Subsection). Figure 8 demonstrates retrograde conduction in the control case (Figure 8A), after FP ablation (Figure 8B), and after SP ablation (Figure 8C) for S1S2 interval of 175 ms.

In the control and SP ablation cases, excitation propagates retrogradely via the FP, creating almost the same A1A2 intervals and with close H1A1 and H2A2 conduction times. In the former case, similar to the case of atrial pacing with long intervals, retrograde and anterograde waves in SP meet and annihilate, but in the upper part of the AVN model (Figure 8A). In the FP ablation case, the A1A2 interval was shorter (Figure 8B), while conduction passing through SP leads to longer H1A1 and H2A2 latencies (see complete retrograde conduction curves in Figure 9).

Animation of S1S2 protocol simulation for retrograde His bundle pacing in the control case is shown in Supplementary Video S2.

Applying the S1S2 stimulation protocol to the atrial and HB cells, we constructed various conduction curves for control (pre-ablation) and post-ablation cases. The recovery curve represents the atrial-His conduction time A2H2 for anterograde conduction or His-atrial conduction time H2A2 for retrograde conduction plotted either versus the S1S2 interval (A1A2 for atrial pacing and H1H2 for HB pacing) or versus the recovery interval (H1A2 for atrial pacing and A1H2 for HB pacing) (Mazgalev et al., 1997).

According to the experimental findings, the AVN conduction curves demonstrate pronounced direction-dependent behavior (Reid et al., 2003; Billette and Tadros, 2019). Figure 9 shows the comparison of experimentally obtained (left panels) (modified from (Reid et al., 2003)) and simulated (right panels) anterograde and retrograde recovery curves in the control and ablation cases plotted versus S1S2 intervals. The modification from the dependence on recovery time H1A2 to S1S2 (A1A2) interval was performed by shifting each conduction curve on its minimal conduction time (A1H1) value according to formula A1A2 = H1A2 + A1H1 (Mazgalev et al., 1997; Tadros and Billette, 2012).

Simulated anterograde recovery curve (Figure 9A, right panel) in the control case demonstrates a typical exponential-like rise with reducing pacing interval (Reid et al., 2003; Climent et al., 2011; Masè et al., 2012) with a pronounced bend and tilting at the point where conduction switches from the FP to SP with smooth change of the conduction curve from flat to steep without gap, which is typical for rabbit AVN (Reid et al., 2003; Climent et al., 2011; Billette and Tadros, 2019). At long S1S2 intervals, the SP conduction is overcome by the FP conduction (flat portion of the control curve, see also Figure 5A), while at short S1S2 intervals, excitation propagates through the SP only (Figures 5B, 6A) due to longer refractory period of the FP (Figure 3). In the S1S2 range of 95–122 ms, the perpetual AVNRT was observed (Lin et al., 2001), as mentioned above (Figure 6A).

The ERPN is set as the longest S1S2 pacing interval that results in neither HB nor atrial response in the cases of atrial pacing (anterograde propagation) and HB pacing (retrograde propagation), respectively. FP ablation markedly increased conduction time in both anterograde (Lin et al., 2001; Reid et al., 2003) and retrograde (Reid et al., 2003) conduction cases. The ERPN was 94 ms for anterograde conduction and did not change after FP ablation (91 ± 10 ms in (Reid et al., 2003)). For the retrograde conduction, the ERPN was 163 ms and also did not change after FP ablation (143 ± 16 ms before and 149 ± 14 ms after FP ablation in (Reid et al., 2003)).

During our initial simulations, the value of the upward shift of the retrograde control curve was relatively small (dashed line in Figure 9), while the experiments (Reid et al., 2003) demonstrated much slower retrograde conduction in the control case. The introduction of the asymmetry to the intercellular coupling in the FP, PB, and upper part of HB allowed additional slowing of the retrograde control conduction, better reflecting the AVN direction-dependent conduction properties (Billette and Tadros, 2019).

SP ablation blocked its conduction, while FP conduction remained unaffected. The ablation led to the prolongation of ERPN to 132 ms (141 ± 15 ms in (Reid et al., 2003)) and the disappearance of the steep part of the anterograde conduction curve in the short range of S1–S2 intervals (Figure 9B). The ablation also eliminated the AVNRT, as observed in clinical practice (Katritsis and Josephson, 2013) (compare Figures 6A, 7A).

In the retrograde case, SP ablation did not affect retrograde conduction and SP post-ablation curve virtually coincides with the control curve (Figure 9B). This, together with conduction scenarios shown in Figures 8A,C, confirms retrograde conduction through the FP in the control case. The ERPN was also prolonged to approximately the same value (163 ms) as in the retrograde pre-ablation case, matching with the experimental results (151 ± 12 ms in (Reid et al., 2003)).

Figure 10 shows three main anterograde AVN functional curves (refractory curve H1H2, His-atrial curve H1A2, and recovery curve A2H2) obtained in experiments (left panels) and our simulations (right panels). The AVN functional refractory period (FRPN) is determined as the shortest interval between two consecutive impulses propagated from the atria to the HB (H1H2). In our simulations, FRPN was 180.5 ms, which is in agreement with experiments (175 ± 10 ms in (Reid et al., 2003).

Each point of the refractory curve H1H2 obtained in our simulations equals the sum of A2H2 value on the recovery curve and H1A2 value on the His-atrial curve at the corresponding S1S2 interval (Figure 10A, right), which is in line with the previously proposed quantitative relationship between the three function curves (Tadros and Billette, 2012) (Figure 10A, left). Figure 10B demonstrates another comparison of experimental (Xu et al., 2006) and simulated AVN refractory curves before and after SP ablation. SP ablation decreased the left steeply rising part of the refractory curve compared to the control curve, resulting in an increase in ERPN, but not FRPN.

Note that the experimental conduction curves shown in Figures 9, 10 were obtained for different rabbit preparations and slightly different experimental setups (Reid et al., 2003; Xu et al., 2006; Tadros and Billette, 2012). In particular, the AVN conduction measurements before and after FP and SP ablations were performed separately for each group of rabbits, thus making the experimental conduction curves rather different in Figures 9, 10. On the other hand, our model parameters were adjusted for the rabbit data corresponding to Figure 9A, thus our results show a greater quantitative difference with the experimental conduction curves demonstrated in Figures 9B, 10.

It has been demonstrated that dual pathways are involved in AVN conduction during AFB by limiting the number of atrial impulses transmitted to the ventricles (Zhang et al., 2003). To show the AVN filtering function and reveal the interaction between the pathways during AFB we set random atrial pacing with uniform distribution in the 75–125 ms range (Figure 11). In Figure 11, different conduction sequences for each atrial impulse are present, including complete passage of the excitation to HB through either pathway with the domination of the SP conduction over the FP (which may depend on mean random atrial rate), and simultaneous conduction block in both pathways. One can also see many episodes of retrograde propagation annihilating with an anterograde impulse resulting from a subsequent atrial stimulus. Besides, the sinus rhythm was completely suppressed by faster atrial pacing. This filtering behavior during AFB corresponds to the results observed in experiments (Zhang et al., 2003; Zhang et al., 2004) and in simulations with complex 1D ion-channel based model (Inada et al., 2009; Inada et al., 2017).

As an additional validation of the AVN model filtering function, we simulated the effect of SP and FP ablations on the AFB. Figure 12A shows laddergrams of SP ablation (top panel) and FP ablation (bottom panel) cases simulated with the same random atrial impulse pattern, as shown in Figure 11. The HB rate was visually lower for the SP ablation case. This is also confirmed by the statistical distributions obtained for all three cases simulated with random atrial stimulation during 350 s (Figure 12B). The reduction of the mean HB rate from 339 bpm in the control case to 292 bpm in the SP ablation case and to 318 bpm in the FP ablation case and prolongation of the mean HH interval from 177 ms to 205 ms and 189 ms, respectively, revealed SP ablation as a relatively more effective way of ventricular rate control than FP ablation, which is in agreement with experimental (Zhang et al., 2004) and clinical (Feld et al., 1994) findings, and with mathematical modeling study (Climent et al., 2011).

Dual pathway electrophysiology was shown to be also directly related to the occurrence of the Wenckebach phenomenon (Climent et al., 2011; Zhang and Mazgalev, 2011). Figure 13 demonstrates an example of regular atrial pacing mimicking AFL that yields 5:4 Wenckebach periodicity corresponding to four out of five atrial impulses conducted to the HB, obtained at a regular atrial pacing interval of 125.2 ms (atrial/His rates 480/384 bpm). Figure 13A shows laddergram and action potential propagation through the FP and SP, while Figure 13B demonstrates AH conduction times (top panel) and HH intervals (bottom panel) of consecutive beats in the Wenckebach cycle. This “typical” Wenckebach pattern is characterized by the progressive increase of AVN conduction times (AH) and progressive shortening of HH intervals before the conduction block occurs (Hansom et al., 2021). Also, exactly in line with the experimentally obtained HB electrograms, the first atrial stimulus after the conduction block passed through the FP (depicted by red squares in Figure 13B), followed by a transition to the SP until its ultimate block (Climent et al., 2011; Zhang and Mazgalev, 2011).

In this work, we have developed a compact one-dimensional mathematical model of rabbit AVN extended structure, which includes SN and HB. The AVN model is based on the Aliev-Panfilov non-linear system of two-variable differential equations, describing both quiescent excitable and pacemaking cells. Despite its simplicity and a small number of model cells (33 elements), the model is multi-functional and does not require structural changes for reproduction of a wide variety of AVN dual-pathway electrophysiological phenomena. In particular, the model can correctly reproduce normal sinus rhythm, AVN subsidiary pacemaking, AVNRT, filtering AVN function during AFB and AFL with Wenckebach patterns, and situations after FP and SP ablations.

For the first time, we simulated the retrograde conduction in the AVN and implemented the asymmetry of the intercellular coupling in the FP, PB, and upper part of HB for accurate reproduction of AVN direction-dependent conduction properties (Billette and Tadros, 2019). Several conduction curves (Figures 9, 10) reflecting the general behavior of AVN were presented along with experimental data to demonstrate their good agreement.

Direction-dependent conduction can originate from specific 3D structural or functional characteristics (Carmeliet, 2019), such as the source-sink (current-to-load) mismatch (Kléber and Rudy, 2004), or expression of connexins in heterotypic gap junctions (Temple et al., 2013; Ye et al., 2017). The acceptable direction-dependent conduction behavior was achieved in the initial simulations with the spatial dispersion of action potential duration and gradients of local refractoriness and coupling coefficient (see Figure 3). However, the upward shift of the retrograde control curve (dashed line in Figure 9) was insufficient to reproduce the much slower retrograde conduction in the control case shown in experiments (Reid et al., 2003). This led us to conclude that an additional factor is responsible for slowing down retrograde conduction in the FP. We hypothesized that the asymmetry of the intercellular coupling (Eq. (4)) plays the role of this factor of slow retrograde conduction. The gap junctions with asymmetric or rectifying electrical transmission properties are well known in neuroscience (Palacios-Prado et al., 2014). Rectifying coupling was also implemented in the simulation of the small intestine of mice (Parsons and Huizinga, 2020). The introduction of the asymmetry to the intercellular coupling in the FP, PB, and upper part of HB allowed additional slowing of the retrograde control conduction, better reflecting the AVN direction-dependent conduction properties (Reid et al., 2003; Billette and Tadros, 2019).

Our model is accompanied by a visualization tool, which discloses processes taking place inside the AVN “black box” by creating laddergrams and revealing the leading pathway for anterograde propagation. With the help of the laddergrams, it became possible to compare FP and SP functional state before and after ablation, including hidden SP activation over a broad range of atrial and HB pacing intervals, as well as the involvement of dual AVN pathways in the Wenckebach phenomena and during AFB.

We believe the proposed AVN model is suitable for cardiac research and educational purposes. The spatial dispersion of local refractoriness, coupling coefficient, and coupling asymmetry shown in Figure 3 can serve as a template for implementing ionic cellular models with different levels of complexity in the AVN modeling. Our model can be adjusted for other animal species and humans, taking into account specific features, for example, a characteristic gap (“jump”) in the conduction curve of human AVN (Denes et al., 1973). The model can be integrated into the whole-heart models (Trayanova, 2011) and simulation frameworks such as openCARP (Plank et al., 2021) for near-realtime or even real-time simulations (Gillette et al., 2022). Another possible application of the model is utilization in the real-time testbed systems for implantable cardiac devices (Ai et al., 2020; Biasi et al., 2022). The proposed model can be also used for quick preliminary testing of hypotheses before accurate simulations with computationally expensive and biophysically detailed cellular models.

Similar to the complex multi-cellular ion-channel model (Inada et al., 2009) and functional mathematical model (Climent et al., 2011), our model is based on experimental rabbit data. Due to the small number of elements and much simpler non-linear equation systems, our model allows faster computation and easier tuning than more biophysically detailed model (Inada et al., 2009), keeping the model structure unchanged. On the other hand, functional mathematical models like (Climent et al., 2011) benefit from effortless fitting of exponential equations to particular rabbit cardiac preparation. The presented model is capable of reproducing almost all electrophysiological properties of AVN found in both types of models with the additional feature of retrograde conduction with realistic ERPN in control and FP and SP post-ablation cases.

The concept of the presented model is based on monodomain formulation of electrical propagation in a 1D cable structure where FP and SP are electrically isolated from each other. It is possible that some macroscopic (anatomical) and microscopic (cytological and extracytologic matrix) properties (Kawashima and Sato, 2021) related to the complexity of real 3D AVN structure were not taken into account in the proposed AVN model. This can also explain the quantitative difference between experimental and simulated anterograde conduction curves in the range of short S1S2 pacing cycles and the retrograde conduction cases.

Normal conduction case shown in Figure 4A assumes that retrograde conduction exists in the SP. In Figure 4A, one can see that upon the excitation wave in the FP reaches the PB, the retrograde wave goes up through the SP, and the annihilation of the anterograde and retrograde waves in the SP occurs at about 100 ms. However, it was shown in the experiments (Reid et al., 2003) that in some rabbit heart preparations, retrograde conduction totally disappeared after the FP ablation, which means the absence of retrograde conduction in the intact SP. In such a case, normal conduction demonstrated in Figure 4A would look different for the SP. The SP trace will go straight down and end up at the PB due to the refractory state of the latter at about 110 ms. Similarly, in the case of spontaneous activity originating in the SP (Figure 4B), the retrograde part of the excitation wave in the SP will be absent, and retrograde conduction will be provided by the FP only.

Due to the simplicity of the two-variable A-P model used in this study, the proposed AVN model does not allow simulation of subtle effects on cardiac cell ion channels, such as the impact of drugs or change in ion concentrations as demonstrated in the more complex multicellular ion-channel based models (Inada et al., 2017).