A compact multi-functional model of the rabbit atrioventricular node with dual pathways

The atrioventricular node (AVN) is considered a “black box”, and the functioning of its dual pathways remains controversial and not fully understood. In contrast to numerous clinical studies, there are only a few mathematical models of the node. In this paper, we present a compact, computationally lightweight multi-functional rabbit AVN model based on the Aliev-Panfilov two-variable cardiac cell model. The one-dimensional AVN model includes fast (FP) and slow (SP) pathways, primary pacemaking in the sinoatrial node, and subsidiary pacemaking in the SP. To obtain the direction-dependent conduction properties of the AVN, together with gradients of intercellular coupling and cell refractoriness, we implemented the asymmetry of coupling between model cells. We hypothesized that the asymmetry can reflect some effects related to the complexity of the real 3D structure of AVN. In addition, the model is accompanied by a visualization of electrical conduction in the AVN, revealing the interaction between SP and FP in the form of ladder diagrams. The AVN model demonstrates broad functionality, including normal sinus rhythm, AVN automaticity, filtering of high-rate atrial rhythms during atrial fibrillation and atrial flutter with Wenckebach periodicity, direction-dependent properties, and realistic anterograde and retrograde conduction curves in the control case and the cases of FP and SP ablation. To show the validity of the proposed model, we compare the simulation results with the available experimental data. Despite its simplicity, the proposed model can be used both as a stand-alone module and as a part of complex three-dimensional atrial or whole heart simulation systems, and can help to understand some puzzling functions of AVN.

1 Introduction

The atrioventricular node (AVN) plays a key role in the cardiac electrical conduction system. It is located between the atria and ventricles, namely, in the base of the right atrium (Figure 1A). The AVN is the only site responsible for transmitting impulses originating in the sinoatrial node (SN) to the ventricles of the heart, coordinating the relationship between contraction periods of the heart chambers, and introducing a variable delay, allowing effective pumping of the blood in a wide range of cardiac rhythms (George et al., 2017; Billette and Tadros, 2019). Figure 1B demonstrates the detailed anatomical location of the rabbit AVN and surrounding tissues (de Carvalho et al., 1959).

FIGURE 1

FIGURE 1. (A) Conduction system of the mammalian heart. SN—sinus node, AVN—atrioventricular node, RA and LA—right and left atria, PB—penetrating bundle, HB—His bundle, TV—tricuspid valve, RV and LV—right and left ventricles. (B) Anatomic drawing of the rabbit atrioventricular node and surrounding tissues with excitation wave latency (from (de Carvalho et al., 1959)). CT—crista terminalis, CS—coronary sinus, TZ—transitional zone, and INE—inferior nodal extension. The numbers in milliseconds are given with respect to the sinus node pacemaker. The fast and slow pathways are indicated by red and blue lines, respectively.

The rabbit AVN contains two built-in functional pathways conducting impulses from the atria to His bundle (HB)—so-called fast pathway (FP) with a longer effective refractory period, and slow pathway (SP) with a shorter one (Figure 1B). The FP also has a smaller conduction latency than the SP (Moe et al., 1956; Medkour et al., 1998; Nikolski et al., 2003). The dual pathway AVN physiology provides protection of ventricles from atrial fibrillation (AFB) and supraventricular tachycardia by filtering excessive atrial electrical impulses (Billette and Tadros, 2019). On the other hand, the interplay of the FP and SP in the AVN can be a potential substrate for atrioventricular nodal re-entrant tachycardia (AVNRT) (Nikolski et al., 2003; Katritsis and Josephson, 2013). The AVN also acts as a subsidiary pacemaker and provides atrial and ventricular pacing in the case of SN failure. According to experimental studies, in such a case, the AVN pacemaker activity in the rabbit heart originates in the SP between the coronary sinus and the tricuspid valve in the inferior nodal extension (Dobrzynski et al., 2003).

Despite numerous research efforts, many aspects of the AVN electrophysiological behavior remain unclear and debated (George et al., 2017). In addition to animal experiments and clinical studies on humans, mathematical modeling at different levels of detalization and complexity may offer insight into the AVN “black box”. Currently, very few AVN models or atrial models incorporating AVN with dual conduction pathways exist (Inada et al., 2009; Climent et al., 2011; Masè et al., 2012; Podziemski and Żebrowski, 2013; Li et al., 2014; Ryzhii and Ryzhii, 2019). The rabbit cardiac biophysically-detailed multicellular ion-channel model (Inada et al., 2009) is quite complex and requires reconfiguration, e.g., changing the number of model cells for simulation of different AVN phenomena. Moreover, some reproducibility issues were reported recently (Schölzel et al., 2021). The model computes normal action potential propagation from the SN through the atrium to AVN, AVNRT, the effect of the FP and SP ablation, I_Ca,L block effect, the AVN filtering function during AFB, and AVN pacemaking.

There were also attempts to develop AVN models based on more computationally effective simplified non-linear two-variable models. In the work (Podziemski and Żebrowski, 2013) the AVN model based on hundreds of Lienard-transformed modified Van der Pol (van der Pol and van der Mark, 1928) and FitzHugh-Nagumo (FitzHugh, 1961; Nagumo et al., 1962) equations, reproduces normal sinus rhythm and AVNRT. The three-dimensional anatomically-detailed model of the rabbit right atrium (Li et al., 2014) includes the AVN with dual pathways and utilizes the cellular automaton method and modified FitzHugh-Nagumo two-variable models. The model simulates normal sinus rhythm and nodal behavior during AFB. In work (Ai et al., 2018) the human real-time hybrid automaton model was proposed, which can emulate sinus rhythm and the AVNRT. In the works (Climent et al., 2011; Masè et al., 2012) the authors used a simple, functional mathematical model of the rabbit AVN with the inclusion of dual pathway physiology. This model consists of exponential mathematical equations and is based on rabbit heart preparation data, and is able to predict AVN conduction time and the interaction between the FP and SP wavefronts during regular and irregular atrial rhythms. The functions of the models include simulation of nodal recovery, concealed conduction, Wenckebach phenomena, AFB and atrial flutter (AFL), and FP and SP ablations. The network model of human AVN (Wallman and Sandberg, 2018) consists of FP and SP interconnected at one end, which includes functions of the HB and Purkinje network. Both pathways are composed of interacting nodes, and the refractory periods and conduction delays of the latter are determined by the stimulation history of each node and described by exponential functions (Lian et al., 2006). The model is able to represent RR interval series extracted from ECG data (both in the forms of histograms, Poincaré plots, and autocorrelation), and the latest version also incorporates autonomic tone (Plappert et al., 2022).

In our previous work (Ryzhii and Ryzhii, 2019), we considered a simplified human AVN model with dual pathways, which consists of about 40 cells including excitable cells and pacemaking cells described by Aliev-Panfilov (A-P) (Aliev and Panfilov, 1996) and modified Van der Pol models. The model reproduces normal sinus rhythm, AVN automaticity, reentry, and filtering function. However, we found that various variants of the modified Van der Pol equation are not suitable for the description of cells of subsidiary pacemakers.

In this work, we present a compact mathematical model of the rabbit AVN based on a single non-linear two-variable A-P model for the description of both quiescent excitable and pacemaking cells, which allows representing various parts of cardiac conduction system with different tissue types in a uniform and convenient way. The model consists of only 33 cells and requires relatively low computational resources. Adjusting to available experimental data for the rabbit heart made it possible to demonstrate a wide range of behavioral phenomena of AVN. This includes proper bidirectional conduction in the FP and SP both in the normal case and after the FP and SP ablations reflected in the calculated recovery curves, AVNRT, the AVN automaticity, filtering function during AFB and AFL, and Wenckebach periodicity. To achieve correct reproduction of the retrograde conduction characteristics, we assumed the presence of an asymmetry of the coupling between the AVN cells and implemented it in the model. In addition, we included visualization demonstrating the inner processes in the AVN “black box”, in particular, the laddergrams (Lewis ladder diagrams) (Johnson and Denes, 2008) with the exact timing of each model cell in both FP and SP pathways during the excitation wave propagation. Finally, the comparison with published experimental findings provides validation of the simulated results.

2 Materials and methods

2.1 Structure and equations of the model

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.

FIGURE 2

FIGURE 2. Schematic representation of the rabbit cardiac conduction model with dual pathways. SN - sinus node cell, PS1–PS3 - peripheral sinus node cells, AM1–AM3 - right atrial muscle cells, FP1–FP8 - fast pathway cells, SP1–SP10 - slow pathway cells, PB - penetrating bundle cell, HB1–HB6 - His bundle cells. Arrows Antero and Retro denote the points for anterograde and retrograde S1S2 stimuli application, arrow AF indicates the point of stimuli application for random pacing (atrial fibrillation) and regular pacing (atrial flutter), and AM* is an additional intermediate atrial muscle cell. The gray-shaded cells indicate primary and subsidiary pacemaker cells.

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:

$\frac{\partial u}{\partial t}=c_t\left[ku\left(u-a_1\right)\left(1-u\right)-vu\right]+I^{\mathit{coupl}}+I^{\mathit{stim}},(1)$

$\frac{\partial v}{\partial t}=c_t\epsilon \left[-v-ku\left(u-a_2-1\right)\right],(2)$

$\epsilon ={\epsilon }_0+\frac{v{\mu }_1}{u+{\mu }_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^{\mathit{coupl}}=d_{i-1}\left(u_{i-1}-u_i\right)+d_i\left(u_{i+1}-u_i\right),(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^{\mathit{coupl}}=d_{i-1}\left(u_{i-1}-{\alpha }_{i-1}{u}_i\right)+d_i\left(-u_i+{\alpha }_i{u}_{i+1}\right),(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).

2.2 Parameters

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.

FIGURE 3

FIGURE 3. Distributions of refractory period for uncoupled model cells, coupling coefficient d, and coefficient of coupling asymmetry α. Opaque and filled circles denote excitable and pacemaking cells, respectively. For d and α plots, the data in SP with indexes 6.5 and 15.5 correspond to “verical” coupling between the rows of the AVN model.

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).

2.3 External pacing

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.

2.4 Visualization of laddergrams

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).

2.5 Numerical methods

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).

3 Results and their validation

3.1 Remarks on figures of excitation propagation

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.

FIGURE 4

FIGURE 4. Primary and subsidiary pacemaking. Here and after top scheme illustrates propagation conduction sequences of interacting fast (FP, red arrows) and slow (SP, blue arrows) pathways. Top panel: the laddergram (Lewis ladder diagram) with the exact timing of each model cell excitation. Smaller circles denote quiescent excitable model cells, and larger circles - natural pacemaker cells. Red and blue traces correspond to the propagation through the FP and SP, respectively. The trace color in the HB part corresponds to the leading pathway. Dark red arrowhead markers on the right side and the adjacent lines denote the points of latency time measurements. Middle panel: the sequence of action potentials passing through the FP between the SN and HB. Bottom panel: the sequence of action potentials passing through the SP between the SN and HB. Action potentials in red correspond to the latency time measurement points, and in black color - to the pacemaking cells. Time scale is common for all three panels. (A) Rabbit’s normal sinus rhythm of 166 bpm (361 ms interval). The numbers on the traces indicate excitation wave latency with respect to the SN in milliseconds. (B) In the case of a silent SN pacemaker, the spontaneous activity with 100 bpm rate (600 ms interval) originates in the SP.

FIGURE 5

FIGURE 5. Atrial pacing in the control case with S1S2 interval of 140 ms and 130 ms. (A) At 140 ms interval, the excitation wave arrives via the FP and SP to the PB almost simultaneously with a slight advance of the FP. (B) The conduction in the FP is stopped due to a functional conduction block. In the top panel, arrowhead markers indicate points of stimulation. Upper (between the markers) and lower numbers correspond to the S1S2 and H1H2 intervals, respectively. Short stimulation impulses are shown at the top of the middle panel.

FIGURE 6

FIGURE 6. (A) Atrial pacing in the control case with S1S2 interval of 110 ms. Perpetual SP–FP AVNRT. (B) Atrial pacing for the control case with S1S2 interval of 94 ms. Full nodal conduction block.

FIGURE 7

FIGURE 7. SP and FP ablations. (A) Atrial pacing with S1S2 interval of 110 ms after SP ablation. See Figure 6A for comparison. (B) Atrial pacing with S1S2 interval of 140 ms after FP ablation. See Figure 5A for comparison.

FIGURE 8

FIGURE 8. Retrograde AVN conduction, HB pacing with 175 ms S1S2 interval. (A) Control case. (B) After FP ablation. (C) After SP ablation. Upper and lower numbers correspond to A1A2 and S1S2 (H1H2) intervals. Short stimulation impulses are shown at the bottom of the middle panel.

3.2 Normal sinus rhythm

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).

3.3 AVN automaticity

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.

3.4 Atrial pacing

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.

3.5 Effects of FP and SP ablations on anterograde conduction

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).

FIGURE 9

FIGURE 9. Comparison of experimental (Reid et al., 2003) (left panels) and simulated (right panels) anterograde and retrograde recovery curves in the control case and after FP and SP ablations. (A) Experimental recovery curves in the control case and after FP ablation and simulated recovery curves. (B) Experimental recovery curves in the control case and after SP ablation and simulated recovery curves. Dashed lines correspond to the simulated retrograde conduction in the control case without coupling asymmetry.

3.6 Retrograde conduction and FP and SP ablations

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.

3.7 AVN conduction curves

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).

FIGURE 10

FIGURE 10. Comparison of experimental (left panels) and simulated (right panels) conduction curves. (A) Experimental (Tadros and Billette, 2012) and simulated AVN refractory, His-atrial, and recovery curves. (B) Experimental (Xu et al., 2006) and simulated refractory curves before and after SP ablation. The simulated control refractory curves (H1H2) in the right panels are the same.

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.

3.8 Filtering of atrial rhythm during atrial fibrillation

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).

FIGURE 11

FIGURE 11. Filtering of atrial rhythm during atrial fibrillation stimulated by random impulse sequence in 75–125 ms range applied to AM* model cell. The short stimulation impulses and resulting AM* cell action potentials are shown in the top part of the middle panel.

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).

FIGURE 12

FIGURE 12. Impact of ablation on AVN filtering function. (A) Laddergrams of atrial fibrillation stimulated by the same random impulse sequence as in Figure 11 after SP ablation (top panel) and after FP ablation (bottom panel). (B) Distribution of consecutive HH intervals for intact AVN and after ablations. Random atrial pacing in the 75–125 ms range was applied for 350 s. Control case (left panel), after SP ablation (middle panel), and after FP ablation (right panel). Mean HH intervals and HB rates for each case are shown at the bottom.

3.9 Wenckebach periodicity during atrial flutter

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).

FIGURE 13

FIGURE 13. Simulated 5:4 Wenckebach pattern (atrial/His rates 480/384 bpm) at regular atrial pacing (atrial flutter) with the interval of 125.2 ms. (A) Laddergram and action potential propagation through the FP and SP. (B) AH conduction time of consecutive beats in the Wenckebach cycle (top panel), and HH interval (bottom panel). Each Wenckebach cycle begins with FP conduction (red squares), AH delays progressively prolonged, and HH intervals progressively shortened.

4 Discussion

4.1 Major achievements

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.

4.2 Comparison with previous AVN 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.

4.3 Limitations

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).

5 Conclusion

In this work, we have developed a compact multi-functional model of a rabbit atrioventricular node with dual pathways. The one-dimensional model is relatively simple, contains a small number of model cells, and is computationally efficient, allowing for fast simulations. Using the developed model, we successfully reproduced several electrophysiological phenomena accompanied by visualization of laddergrams. We believe the proposed AVN model has a high potential for utilization in research and education both as an independent module and as a part of complex three-dimensional whole-heart simulation systems and real-time testbed systems. The proposed model can also be used for quick preliminary testing of hypotheses before accurate and detailed simulations with computationally expensive and complex cellular models.

Data availability statement

The datasets presented in this study can be found in the article/Supplementary Material. The MATLAB source code is available in Figshare repository (https://doi.org/10.6084/m9.figshare.22015289.v1).

Author contributions

MR wrote the program code and performed computer simulations. ER analyzed the results. All authors contributed to the conception and design of the study and manuscript preparation and revisions, read and approved the submitted version.

Funding

This work was supported by Grant No. 20K12046, JSPS KAKENHI.

Acknowledgments

The authors would like to thank Katrina Armstrong for her assistance with the manuscript.

Conflict of interest

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2023.1126648/full#supplementary-material

References

Ai, W., Patel, N. D., Roop, P. S., Malik, A., Andalam, S., Yip, E., et al. (2018). A parametric computational model of the action potential of pacemaker cells. IEEE Trans. Biomed. Eng. 65, 123–130. doi:10.1109/TBME.2017.2695537

PubMed Abstract | CrossRef Full Text | Google Scholar

Ai, W., Patel, N. D., Roop, P. S., Malik, A., and Trew, M. L. (2020). Cardiac electrical modeling for closed-loop validation of implantable devices. IEEE Trans. Biomed. Eng. 67, 536–544. doi:10.1109/TBME.2019.2917212

PubMed Abstract | CrossRef Full Text | Google Scholar

Aliev, R. R., and Panfilov, A. V. (1996). A simple two-variable model of cardiac excitation. Chaos Solit. Fractals 7, 293–301. doi:10.1016/0960-0779(95)00089-5

CrossRef Full Text | Google Scholar

Biasi, N., Seghetti, P., Mercati, M., and Tognetti, A. (2022). A reaction-diffusion heart model for the closed-loop evaluation of heart-pacemaker interaction. IEEE Access 10, 121249–121260. doi:10.1109/access.2022.3222830

CrossRef Full Text | Google Scholar

Billette, J. (1987). Atrioventricular nodal activation during periodic premature stimulation of the atrium. Am. J. Physiol. 252, H163–H177. doi:10.1152/ajpheart.1987.252.1.H163

PubMed Abstract | CrossRef Full Text | Google Scholar

Billette, J., and Tadros, R. (2019). An integrated overview of AV node physiology. Pacing Clin. Electrophysiol. 42, 805–820. doi:10.1111/pace.13734

PubMed Abstract | CrossRef Full Text | Google Scholar

Carmeliet, E. (2019). Conduction in cardiac tissue. Historical reflections. Physiol. Rep. 7, e13860. doi:10.14814/phy2.13860

PubMed Abstract | CrossRef Full Text | Google Scholar

Climent, A. M., Guillem, M. S., Zhang, Y., Millet, J., and Mazgalev, T. N. (2011). Functional mathematical model of dual pathway AV nodal conduction. Am. J. Physiol. Heart Circ. Physiol. 300, H1393–H1401. doi:10.1152/ajpheart.01175.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

de Carvalho, A. P., de Mello, W. C., and Hoffman, B. F. (1959). Electrophysiological evidence for specialized fiber types in rabbit atrium. Am. J. Physiol. Leg. Content 196, 483–488. doi:10.1152/ajplegacy.1959.196.3.483

CrossRef Full Text | Google Scholar

Denes, P., Wu, D., Dhingra, R. C., Chuquimia, R., and Rosen, K. M. (1973). Demonstration of dual A-V nodal pathways in patients with paroxysmal supraventricular tachycardia. Circulation 48, 549–555. doi:10.1161/01.cir.48.3.549

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobrzynski, H., Nikolski, V. P., Sambelashvili, A. T., Greener, I. D., Yamamoto, M., Boyett, M. R., et al. (2003). Site of origin and molecular substrate of atrioventricular junctional rhythm in the rabbit heart. Circ. Res. 93, 1102–1110. doi:10.1161/01.RES.0000101913.95604.B9

PubMed Abstract | CrossRef Full Text | Google Scholar

Feld, G. K., Fleck, R. P., Fujimura, O., Prothro, D. L., Bahnson, T. D., and Ibarra, M. (1994). Control of rapid ventricular response by radiofrequency catheter modification of the atrioventricular node in patients with medically refractory atrial fibrillation. Circulation 90, 2299–2307. doi:10.1161/01.cir.90.5.2299

PubMed Abstract | CrossRef Full Text | Google Scholar

FitzHugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1, 445–466. doi:10.1016/s0006-3495(61)86902-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Garrigue, S., Tchou, P. J., and Mazgalev, T. N. (1999). Role of the differential bombardment of atrial inputs to the atrioventricular node as a factor influencing ventricular rate during high atrial rate. Cardiovasc. Res. 44, 344–355. doi:10.1016/s0008-6363(99)00201-1

PubMed Abstract | CrossRef Full Text | Google Scholar

George, S. A., Faye, N. R., Murillo-Berlioz, A., Lee, K. B., Trachiotis, G. D., and Efimov, I. R. (2017). At the atrioventricular crossroads: Dual pathway electrophysiology in the atrioventricular node and its underlying heterogeneities. Arrhythmia Electrophysiol. Rev. 6, 179–185. doi:10.15420/aer.2017.30.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Gillette, K., Gsell, M., Strocchi, M., Grandits, T., Neic, A., Manninger, M., et al. (2022). A personalized real-time virtual model of whole heart electrophysiology. Front. Physiol. 13, 907190. doi:10.3389/fphys.2022.907190

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansom, S. P., Golian, M., and Green, M. S. (2021). The Wenckebach phenomenon. Curr. Cardiol. Rev. 17, 10–16. doi:10.2174/1573403X16666200719022142

PubMed Abstract | CrossRef Full Text | Google Scholar

Hong, K. L., Borges, J., and Glover, B. (2020). Catheter ablation for the management of atrial fibrillation: Current technical perspectives. Open Heart 7, e001207. doi:10.1136/openhrt-2019–001207

PubMed Abstract | CrossRef Full Text | Google Scholar

Hucker, W. J., Fedorov, V. V., Foyil, K. V., Moazami, N., and Efimov, I. R. (2008). Images in cardiovascular medicine. Optical mapping of the human atrioventricular junction. Circulation 117, 1474–1477. doi:10.1161/CIRCULATIONAHA.107.733147

PubMed Abstract | CrossRef Full Text | Google Scholar

Hwang, H. J., Ng, F. S., and Efimov, I. R. (2014). “Mechanisms of atrioventricular nodal excitability and propagation,” in Cardiac electrophysiology: From cell to bedside. Sixth Edition (Elsevier), 275–285. doi:10.1016/B978-1-4557-2856-5.00028-5

CrossRef Full Text | Google Scholar

Inada, S., Hancox, J. C., Zhang, H., and Boyett, M. R. (2009). One-dimensional mathematical model of the atrioventricular node including atrio-nodal, nodal, and nodal-His cells. Biophys. J. 97, 2117–2127. doi:10.1016/j.bpj.2009.06.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Inada, S., Shibata, N., Iwata, M., Haraguchi, R., Ashihara, T., Ikeda, T., et al. (2017). Simulation of ventricular rate control during atrial fibrillation using ionic channel blockers. J. Arrhythmia 33, 302–309. doi:10.1016/j.joa.2016.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, N., and Denes, P. (2008). The ladder diagram (a 100+ year history). Am. J. Cardiol. 101, 1801–1804. doi:10.1016/j.amjcard.2008.02.085

PubMed Abstract | CrossRef Full Text | Google Scholar

Katritsis, D. G., and Josephson, M. E. (2013). Classification of electrophysiological types of atrioventricular nodal re-entrant tachycardia: A reappraisal. EP Eur. 15, 1231–1240. doi:10.1093/europace/eut100

PubMed Abstract | CrossRef Full Text | Google Scholar

Katritsis, D. G., and Morady, F. (2022). “Basic intervals and atrial and ventricular conduction curves,” in Clinical cardiac electrophysiology. Editors D. G. Katritsis, and F. Morady (Elsevier), 54–72. doi:10.1016/B978-0-323-79338-4.00015-7

CrossRef Full Text | Google Scholar

Kawashima, T., and Sato, F. (2021). First in situ 3D visualization of the human cardiac conduction system and its transformation associated with heart contour and inclination. Sci. Rep. 11, 8636. doi:10.1038/s41598-021-88109-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerr, C. R., Prystowsky, E. N., Browning, D. J., and Strauss, H. C. (1980). Characterization of refractoriness in the sinus node of the rabbit. Circ. Res. 47, 742–756. doi:10.1161/01.res.47.5.742

PubMed Abstract | CrossRef Full Text | Google Scholar

Kléber, A. G., and Rudy, Y. (2004). Basic mechanisms of cardiac impulse propagation and associated arrhythmias. Physiol. Rev. 84, 431–488. doi:10.1152/physrev.00025.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Inada, S., Schneider, J. E., Zhang, H., Dobrzynski, H., and Boyett, M. R. (2014). Three-dimensional computer model of the right atrium including the sinoatrial and atrioventricular nodes predicts classical nodal behaviours. PLoS ONE 9, e112547. doi:10.1371/journal.pone.0112547

PubMed Abstract | CrossRef Full Text | Google Scholar

Lian, J., Müssig, D., and Lang, V. (2006). Computer modeling of ventricular rhythm during atrial fibrillation and ventricular pacing. IEEE Trans. Biomed. Eng. 53, 1512–1520. doi:10.1109/TBME.2006.876627

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, L. J., Billette, J., Medkour, D., Reid, M. C., Tremblay, M., and Khalife, K. (2001). Properties and substrate of slow pathway exposed with a compact node targeted fast pathway ablation in rabbit atrioventricular node. J. Cardiovasc. Electrophysiol. 12, 479–486. doi:10.1046/j.1540-8167.2001.00479.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Masè, M., Glass, L., Disertori, M., and Ravelli, F. (2012). Nodal recovery, dual pathway physiology, and concealed conduction determine complex AV dynamics in human atrial tachyarrhythmias. Am. J. Physiol. Heart Circ. Physiol. 303, H1219–H1228. doi:10.1152/ajpheart.00228.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Mazgalev, T., Mowrey, K., Efimov, I., Fahy, G. J., Wagoner, D. V., Cheng, Y., et al. (1997). Mechanism of atrioventricular nodal facilitation in rabbit heart: Role of proximal AV node. Am. J. Physiol. Heart Circ. Physiol. 273, H1658–H1668. doi:10.1152/ajpheart.1997.273.4.H1658

PubMed Abstract | CrossRef Full Text | Google Scholar

Medkour, D., Becker, A. E., Khalife, K., and Billette, J. (1998). Anatomic and functional characteristics of a slow posterior AV nodal pathway: Role in dual-pathway physiology and reentry. Circulation 98, 164–174. doi:10.1161/01.cir.98.2.164

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitsuoka, T., Mazgalev, T., Dreifus, L. S., and Michelson, E. L. (1987). Differential vagal effects on antegrade vs. retrograde atrioventricular conduction. Am. J. Physiol. 253, H1059–H1068. doi:10.1152/ajpheart.1987.253.5.H1059

PubMed Abstract | CrossRef Full Text | Google Scholar

Moe, G. K., Preston, J. B., and Burlington, H. (1956). Physiologic evidence for a dual A-V transmission system. Circ. Res. 4, 357–375. doi:10.1161/01.res.4.4.357

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagumo, J., Arimoto, S., and Yoshizawa, S. (1962). An active pulse transmission line simulating nerve axon. Proc. IRE 50, 2061–2070. doi:10.1109/jrproc.1962.288235

CrossRef Full Text | Google Scholar

Narula, O. S. (1974). Retrograde pre-excitation. Comparison of antegrade and retrograde conduction intervals in man. Circulation 50, 1129–1143. doi:10.1161/01.cir.50.6.1129

PubMed Abstract | CrossRef Full Text | Google Scholar

Nikolski, V. P., Jones, S. A., Lancaster, M. K., Boyett, M. R., and Efimov, I. R. (2003). Cx43 and dual-pathway electrophysiology of the atrioventricular node and atrioventricular nodal reentry. Circ. Res. 92, 469–475. doi:10.1161/01.RES.0000059304.97120.2F

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Hara, G., Gendreau, R., Billette, J., Amellal, F., Nayebpour, M., Talajic, M., et al. (1993). Rate-dependent functional properties of retrograde atrioventricular nodal conduction in experimental animals. Am. J. Physiol. 265, H1257–H1264. doi:10.1152/ajpheart.1993.265.4.H1257

PubMed Abstract | CrossRef Full Text | Google Scholar

Palacios-Prado, N., Huetteroth, W., and Pereda, A. E. (2014). Hemichannel composition and electrical synaptic transmission: Molecular diversity and its implications for electrical rectification. Front. Cell. Neurosci. 8, art. 324. doi:10.3389/fncel.2014.00324

PubMed Abstract | CrossRef Full Text | Google Scholar

Parsons, S. P., and Huizinga, J. D. (2020). A myogenic motor pattern in mice lacking myenteric interstitial cells of Cajal explained by a second coupled oscillator network. Am. J. Physiol. Gastrointest. Liver. Physiol. 318, G243. doi:10.1152/ajpgi.00311.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Plank, G., Loewe, A., Neic, A., Augustin, C., Huang, Y. L., Gsell, M. A. F., et al. (2021). The openCARP simulation environment for cardiac electrophysiology. Comput. Methods Programs Biomed. 208, 106223. doi:10.1016/j.cmpb

CrossRef Full Text | Google Scholar

Plappert, F., Wallman, M., Abdollahpur, M., Platonov, P. G., Östenson, S., and Sandberg, F. (2022). An atrioventricular node model incorporating autonomic tone. Front. Physiol. 13, 976468. doi:10.3389/fphys.2022.976468

PubMed Abstract | CrossRef Full Text | Google Scholar

Podziemski, P., and Żebrowski, J. J. (2013). Liénard-type models for the simulation of the action potential of cardiac nodal cells. Phys. D. 261, 52–61. doi:10.1016/j.physd.2013.06.007

CrossRef Full Text | Google Scholar

Reid, M. C., Billette, J., Khalife, K., and Tadros, R. (2003). Role of compact node and posterior extension in direction-dependent changes in atrioventricular nodal function in rabbit. J. Cardiovasc. Electrophysiol. 14, 1342–1350. doi:10.1046/j.1540-8167.2003.03382.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryzhii, E., and Ryzhii, M. (2014). A heterogeneous coupled oscillator model for simulation of ECG signals. Comput. Methods Programs Biomed. 117, 40–49. doi:10.1016/j.cmpb.2014.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryzhii, M., and Ryzhii, E. (2019). Optimization of dual pathway AV nodal conduction model. J. Phys. Conf. Ser. 1372, 012078. doi:10.1088/1742-6596/1372/1/012078

CrossRef Full Text | Google Scholar

Ryzhii, M., and Ryzhii, E. (2022). Pacemaking function of two simplified cell models. PLoS ONE 17, e0257935. doi:10.1371/journal.pone.0257935

PubMed Abstract | CrossRef Full Text | Google Scholar

Schölzel, C., Blesius, V., Ernst, G., Goesmann, A., and Dominik, A. (2021). Countering reproducibility issues in mathematical models with software engineering techniques: A case study using a one-dimensional mathematical model of the atrioventricular node. PLoS ONE 16, e0254749. doi:10.1371/journal.pone.0254749

PubMed Abstract | CrossRef Full Text | Google Scholar

Shampine, L. F., and Reichelt, M. W. (1997). The MATLAB ODE suite. SIAM J. Sci. Comput. 18, 1–22. doi:10.1137/s1064827594276424

CrossRef Full Text | Google Scholar

Tadros, R., and Billette, J. (2012). Rate-dependent AV nodal function: Closely bound conduction and refractory properties. J. Cardiovasc. Electrophysiol. 23, 302–308. doi:10.1111/j.1540-8167.2011.02180.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Talajic, M., Papadatos, D., Villemaire, C., Glass, L., and Nattel, S. (1991). A unified model of atrioventricular nodal conduction predicts dynamic changes in Wenckebach periodicity. Circ. Res. 68, 1280–1293. doi:10.1161/01.res.68.5.1280

PubMed Abstract | CrossRef Full Text | Google Scholar

Temple, I. P., Inada, S., Dobrzynski, H., and Boyett, M. R. (2013). Connexins and the atrioventricular node. Heart rhythm. 10, 297–304. doi:10.1016/j.hrthm.2012.10.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Trayanova, N. A. (2011). Whole-heart modeling: Applications to cardiac electrophysiology and electromechanics. Circ. Res. 108, 113–128. doi:10.1161/CIRCRESAHA.110.223610

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Pol, B., and van der Mark, J. (1928). LXXII. The heartbeat considered as a relaxation oscillation, and an electrical model of the heart. Lond. Edinb. Dublin Physiological Mag. J. Sci. 6, 763–775. doi:10.1080/14786441108564652

CrossRef Full Text | Google Scholar

Wallman, M., and Sandberg, F. (2018). Characterisation of human AV-nodal properties using a network model. Med. Biol. Eng. Comput. 56, 247–259. doi:10.1007/s11517-017-1684-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, B., Billette, J., and Lavallée, M. (2006). Concealed conduction in nodal dual pathways: Depressed conduction, prolonged refractoriness, or reset excitability cycle? Heart rhythm. 3, 212–221. doi:10.1016/j.hrthm.2005.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, W. G., Yue, B., Aoyama, H., Kim, N. K., Cameron, J. A., Chen, H., et al. (2017). Junctional delay, frequency, and direction-dependent uncoupling of human heterotypic Cx45/Cx43 gap junction channels. J. Mol. Cell. Cardiol. 111, 17–26. doi:10.1016/j.yjmcc.2017.07.117

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Bharati, S., Mowrey, K. A., and Mazgalev, T. N. (2003). His electrogram alternans reveal dual atrioventricular nodal pathway conduction during atrial fibrillation: The role of slow-pathway modification. Circulation 107, 1059–1065. doi:10.1161/01.cir.0000051464.52601.f4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Bharati, S., Sulayman, R., Mowrey, K. A., Tchou, P. J., and Mazgalev, T. N. (2004). Atrioventricular nodal fast pathway modification: Mechanism for lack of ventricular rate slowing in atrial fibrillation. Cardiovasc. Res. 61, 45–55. doi:10.1016/j.cardiores.2003.10.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., and Mazgalev, T. N. (2011). AV nodal dual pathway electrophysiology and Wenckebach periodicity. J. Cardiovasc. Electrophysiol. 22, 1256–1262. doi:10.1111/j.1540-8167.2011.02068.x

PubMed Abstract | CrossRef Full Text | Google Scholar