Front. Physiol.Frontiers in PhysiologyFront. Physiol.1664-042XFrontiers Media S.A.162440310.3389/fphys.2025.1624403PhysiologyOriginal ResearchECG-based beat-to-beat assessment of AV node conduction properties during AFKarlsson et al.10.3389/fphys.2025.1624403KarlssonMattias12PlappertFelix2PlatonovPyotr G.3ÖstensonSten4WallmanMikael1†SandbergFrida2*†1Department of Systems and Data Analysis, Fraunhofer-Chalmers Centre, Gothenburg, Sweden2Department of Biomedical Engineering, Lund University, Lund, Sweden3Department of Cardiology, Clinical Sciences, Lund University, Lund, Sweden4Department of Internal Medicine and Department of Clinical Physiology, Central Hospital Kristianstad, Kristianstad, Sweden
Edited by:Mark Potse, UMR5251 Institut de mathématiques de Bordeaux (IMB), France
Reviewed by:Michela Masè, University of Trento, Italy
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Introduction
The refractory period and conduction delay of the atrioventricular (AV) node play a crucial role in regulating the heart rate during atrial fibrillation (AF). Beat-to-beat variations in these properties are known to be induced by the autonomic nervous system (ANS) but have previously not been assessable during AF. Assessing these could provide novel information for improved diagnosis, prognosis, and treatment on an individual basis.
Methods
To estimate AV nodal conduction properties with beat-to-beat resolution, we propose a methodology comprising a network model of the AV node, a particle filter, and a smoothing algorithm. The methodology was evaluated using simulated data and using synchronized electrogram (EGM) and ECG recordings from five patients in the intracardiac atrial fibrillation database. The methodology’s ability to quantify ANS-induced changes in AV node conduction properties was evaluated by analyzing ECG data from 21 patients in AF undergoing a tilt test protocol.
Results
The estimated refractory period and conduction delay matched the simulated ground truth based on ECG recordings with a mean absolute error (± std) of 169±14 ms for the refractory period in the fast pathway; 131±13 ms for the conduction delay in the fast pathway; 67±10 ms for the refractory period in the slow pathway; and 178±28 ms for the conduction delay in the slow pathway. These errors decreased when using simulated ground truth based on EGM recordings. Moreover, a decrease in conduction delay and refractory period in response to head-up tilt was seen during the tilt test protocol, as expected under sympathetic activation.
Discussion
These results suggest that beat-to-beat estimation of AV nodal conduction properties during AF from ECG is feasible, with different levels of uncertainty, and that the estimated properties agree with expected AV nodal modulation.
Atrial fibrillation (AF), characterized by disorganized electrical activity in the atria, is the most common sustained cardiac arrhythmia with an estimated prevalence between 2% and 4% globally (Benjamin et al., 2019). The disorganized electrical activity in the atria leads to rapid and irregular contraction of the atria and ventricles, resulting in an increased risk of mortality, predominantly due to heart failure or stroke (Andrew et al., 2013). Despite extensive research on AF, very little robust evidence exists to inform the best type and intensity of rate control treatment on an individual level (Hindricks et al., 2020; Al-Khatib et al., 2014).
The atrioventricular (AV) node normally functions as the sole electrical connection between the atria and ventricles. During AF, the AV node plays a crucial role in protecting the ventricles from the rapid and irregular impulses originating in the atria. This function is accomplished through two distinct pathways; the fast pathway (FP) and the slow pathway (SP), which converge at the Bundle of His (Kurian et al., 2010). The FP conducts impulses faster than SP but has a longer refractory period. Depending on the refractoriness of its pathways, the AV node can either block an incoming impulse or send it through with a conduction delay. Therefore, the refractory period and conduction delay of the two pathways – here denoted ϕ= [RFP, RSP, DFP, and DSP] – are critical determinants of its filtering capability. The AV node thus serves an essential role in regulating the heart rate during AF, and can functionally be characterized by its properties ϕ.
The autonomic nervous system (ANS) has been shown to contribute to the initiation and maintenance of AF (Shen and Zipes, 2014), suggesting that inter-patient variability in ANS activity might influence individual responses to AF treatment. During normal sinus rhythm, the ANS affects the heart rate primarily through changes to the sinus node automaticity, which can be quantified using heart rate variability (Shaffer and Ginsberg, 2017). However, during AF, the disorganized electrical activity in the atria overrides the organized electrical signals from the sinus node, preventing it from regulating the heart rate. Instead, the ANS affects the heart rate primarily through changes to the atrial fibrillatory rate and AV node conduction properties. Therefore, heart rate variability is not applicable as a tool for quantifying ANS modulation during AF. As an alternative, modulation of AV nodal function could be used to quantify the ANS function during AF. Since the AV node function mainly depends on the refractory period and conduction delay of the two pathways, estimating beat-to-beat changes to these properties might give insights into the ANS function.
Assessing the AV-nodal function under AF is a complex task, since its behavior is influenced by multiple factors such as atrial impulses, autonomic modulation, as well as its intrinsic dynamics and structure. Thus, standard signal processing tools are insufficient, and a model-based analysis is required. Several mathematical models of the AV node have previously been proposed with varying complexity tailored to different intents, including but not limited to Jørgensen et al. (2002), Mangin et al. (2005), Inada et al. (2009), Climent et al. (2011), Masè et al. (2015), Henriksson et al. (2015), Ryzhii and Ryzhii (2023). For clinical application on an individual level, a model should ideally have parameters identifiable from non-invasive recordings. To the best of our knowledge, the only model incorporating the refractory period and conduction delay of both pathways while simultaneously allowing for identification of model parameters based on non-invasive recordings is our previously proposed model (Karlsson et al., 2021). Using this model, the individual 24-h trends of the AV node properties ϕ have previously been estimated with a temporal resolution of 10 min using an error function based on the Poincaré plot of the RR interval series (Karlsson et al., 2022; 2024). However, the ANS is known to modulate AV node conduction with beat-to-beat resolution (Leffler et al., 1994). Thus, beat-to-beat resolution of the AV node conduction property estimates would be preferable for studying the ANS.
Because the Poincaré plot relies on statistical information gathered over a sequence of multiple heartbeats, it is of limited use for beat-to-beat analysis. To increase temporal resolution, we propose a particle filter to estimate ϕ with beat-to-beat resolution using our previously proposed model of the AV node (Karlsson et al., 2021). Particle filters approximate the solution to the filtering problem – estimating the current state of a system (ϕ in our case) based on past and current observations. Due to their ability to leverage information from previous time points effectively, particle filters are suitable for beat-to-beat estimation. Particle filters have proven especially powerful for nonlinear and non-Gaussian problems, and have previously been used for, e.g., atrial flutter detection (Lee et al., 2013), to robustly track heart rate (Nathan and Jafari, 2017), and to automatically annotate ultrasound videos of the fetal heart (Bridge et al., 2017). Moreover, by combining the resulting estimates from a particle filter with a smoothing algorithm, estimates of the current state of a system based on past, current, and future observations can be obtained (Chopin and Papaspiliopoulos, 2020).
This study aims to present and evaluate two particle filter-based frameworks for beat-to-beat assessment of AV node conduction properties, based on intracardiac electrogram (EGM) and ECG data, respectively. The evaluation is done in three steps. Step one is to evaluate the estimation accuracy for the EGM and ECG-based methods on simulated data. Step two is to compare the estimates obtained from synchronized ECG and EGM measurements to estimates derived from ECG measurements only. Finally, step three is to analyze the dynamics of the AV node properties during a tilt test protocol using ECG recordings from 21 patients in order to evaluate the method’s ability to quantify expected changes in AV node conduction properties.
2 Materials and methods
The data used in this study are described in Section 2.1, and are followed by a description of the signal processing used to derive an atrial activation time series (AA series) and a ventricular activation time series (RR series) from the EGM and ECG recordings, in Section 2.2. Furthermore, the network model of the AV node is described in Section 2.3, and the computation of the posterior distribution of ϕ using a particle filter is described in Sections 2.4.1, 2.4.2. Finally, the posterior estimates are refined using a smoothing algorithm, as described in Section 2.4.3.
2.1 Datasets
Two previously obtained datasets are used in this study. The publicly available intracardiac atrial fibrillation database (iafdb) provides synchronized EGM and ECG measurements during AF and is used to assess coherence between EGM and ECG-based estimates of ϕ (Goldberger et al., 2000). Additionally, ECG recordings from a previously conducted tilt test study are used to study modulation of the AV node properties in response to changes in ANS activity (Östenson et al., 2017).
2.1.1 Intracardiac atrial fibrillation database
The iafdb data consist of EGM recordings from four separate regions of the right atrium with synchronized three-lead ECG from eight patients with AF or flutter, sampled at 1,000 Hz (Goldberger et al., 2000). The recordings at the tip of the tricuspid valve annulus are used in this study, due to its proximity to the AV-node entrance. Five recordings contain solely AF and were selected for analysis, with an average patient age of 73 ± 10 years, 60% male, and an average signal duration of 58 ± 7 s. In addition, recordings with the catheter resting against the atrial free wall are used to create realistic simulated data, described in Supplementary Material S1.
2.1.2 Tilt test study
The tilt test study includes ECG recordings from 40 patients with persistent AF (Östenson et al., 2017). For the current study, data with sufficient quality from 21 patients were used (average age of 67 ± 7 years, and 67% male). Eleven patients were excluded due to missing ECG data, five due to the inability of the CardioLund software to detect R peaks (see Section 2.2.1) and three due to failure to extract f-waves (see Section 2.2.2). The tilt test protocol involved standard 12-lead ECG recordings taken between 1 and 3 PM in a quiet room. Participants transitioned from supine position after approximately 5 min to head-down tilt (HDT) position (−30°) for approximately 5 minutes before finally a head-up tilt (HUT) position (+60°) for approximately 5 minutes.
2.2 Signal processing
The frameworks presented in this study for assessing the AV node conduction properties with beat-to-beat resolution rely on simultaneous analysis of the RR and AA series. These are obtained using different signal processing methods depending on whether synchronized EGM and ECG recordings are available, or only ECG recordings, as described below.
2.2.1 RR and AA series from synchronized EGM and ECG
The synchronized ECG and EGM recordings from the iafdb were used to derive the RR and AA series. The RR series is extracted using R-peak detection performed by the CardioLund ECG parser (www.cardiolund.com). The AA series is extracted from the EGM recordings using an iterative method (Ng et al., 2013) following average beat subtraction-based ventricular far-field cancellation and standard pre-processing (Shkurovich et al., 1998; Botteron and Smith, 1995).
2.2.2 RR and AA series from ECG
Using solely ECG, the RR series is again extracted from the R-peak detection performed by the CardioLund ECG parser. However, the AA series cannot be extracted from the ECG. Instead, multiple AA series are generated for each RR interval based on the f-wave characteristics of the corresponding ECG segment. Each AA series, denoted α, is generated by a Gaussian random walk described by a mean (μα) and standard deviation (σα). Both μα and σα are estimated based on the f-wave signal extracted from the ECG by applying QRST-cancellation using the CardioLund ECG parser, before a harmonic model (Henriksson et al., 2018) is fitted to the f-waves to estimate the f-wave frequency and a signal quality index (SQI), sampled at 50 Hz, as described by Abdollahpur et al. (2022). For each RR interval, the mean (μf) and standard deviation (σf) of the inverse f-wave frequency are calculated, as well as the SQI. Subsequently, μα is drawn from N(μf,max(0,0.3−SQI)4), where a SQI greater than 0.3 is deemed sufficient based on previous studies (Henriksson et al., 2018), and the exponent 4 is chosen to get a quadratic decrease on the variance for SQI below 0.3. Further, σα is set to 4σf, where the factor 4 is chosen empirically.
2.3 Network model of the AV node
Our previously introduced network model of the AV node (Karlsson et al., 2021) describes it as two pathways (FP and SP). Each pathway comprises 10 nodes, interconnected with a coupling node at the end, which can transmit impulses to the ventricles (see Figure 1). Each node corresponds anatomically to a localized section of its respective pathway, while the coupling node represents the Purkinje fibers and Bundle of His (Kurian et al., 2010).
A schematic representation of the network model where the yellow node represents the coupling node, the red nodes the SP, the green nodes the FP, and arrows the direction for impulse conduction. For readability, only a subset of the 21 nodes is shown (Karlsson et al., 2021).
Diagram illustrating a cardiac conduction model. Input as AA interval series leads from Atria through AV node's slow and fast pathways. Paths converge at a coupling node, proceeding to output as RR interval series in the Ventricles.
The AA series extracted from data (see Section 2.2) arrives at the first nodes of the FP and the SP simultaneously. Each node can be refractory (blocking impulses) or non-refractory (transmitting impulses). Transmitted impulses arrive at adjacent nodes with an added conduction delay, and nodes immediately become refractory after transmitting an impulse. The refractory period (Ri(n)) and conduction delay (Di(n)) for node i are updated for each incoming impulse n according to Equations 1–3,Rin=Rmin+ΔR1−e−t̃in/τRDin=Dmin+ΔDe−t̃in/τD,t̃in=tin−tin−1+Rin−1,where t̃i(n) is the diastolic interval preceding impulse n and ti(n) is the arrival time of impulse n at node i. When t̃i(n)<0, the node is in its refractory state and will block incoming impulses. Each pathway has three parameters for the refractory period and three for the conduction delay, totaling 12 model parameters θ=[RminFP,ΔRFP,τRFP,RminSP,ΔRSP,τRSP,DminFP,ΔDFP,τDFP,DminSP,ΔDSP,τDSP]. The coupling node’s refractory period is fixed to the shortest RR interval in the data used for estimation minus 50 ms, and its conduction delay is fixed at 60 ms (Karlsson et al., 2021).
The model is evaluated using a modified version of Dijkstra’s algorithm (Wallman and Sandberg, 2018). Impulses are propagated through the network in an event-based fashion where the impulse with the lowest ti(n) in a queue (q(t)) of impulses is propagated next (or blocked depending on t̃i(n)). Three types of data are needed to run the model. First, the parameter vector θ is necessary, corresponding to the properties of the AV node. Second, the queue q(t) is required, where each impulse is represented by a tuple containing arrival time ti(n) and node index i. These impulses can arrive from the atria or from a transmitting node within the model. The atrial activation times from an AA series are placed in q(t) with corresponding node index for the first nodes in the slow and fast pathway. Finally, a vector containing the repolarization times (RT(t)) for each of the 21 nodes in the model is needed, corresponding to the end of the diastolic interval (t̃i(n)). The model code and a basic user example can be found at https://github.com/FraunhoferChalmersCentre/AV-node-model.
The AV node model parameters θ are assumed to be fixed between two heartbeats and the corresponding AV node conduction properties ϕ are estimated in the following way: for a given pathway, e.g., FP, RFP and DFP are estimated as the medians of all of Ri(n) and Di(n) (Equations 1, 2) in the time interval, with DFP multiplied by 10 to account for the cumulative delay of the whole pathway. Values for SP are computed analogously.
2.4 Parameter estimation
To estimate ϕ with beat-to-beat resolution, a particle filter is first used to solve the filtering problem – estimating the current state of the system based on current and past observations – before a smoothing algorithm is applied to the estimated states to solve the smoothing problem – estimating the current state of the system based on current, past, and future observations. Moreover, two versions of the particle filter were developed, one designed for RR and AA series extracted from ECG and EGM recordings (EGM-PF), and one designed for AA series generated based on the f-wave frequencies and RR series extracted from the ECG (ECG-PF).
2.4.1 EGM-particle filter
A basic particle filter can be described by its four phases: initialization, weighting, resampling, and propagation. These phases all affect the particles in the particle filter. In this work, each particle corresponds to a model parameter vector θ̂k,j, where k denotes the RR interval index (also referred to as time step) and j is the particle index. The EGM-PF is initialized by drawing N=1,000,000 particles independently from a twelve-dimensional uniform distribution (ranges found in Supplementary Material S2), where particles with an SP refractory period greater than the FP refractory period or an SP conduction delay less than the FP conduction delay are excluded. Initialization is followed by a weighting phase. This starts by evaluating all θ̂1,j with the model, using the current repolarization times (RTj(0)=0) and the current queue (qj(0)) filled by the AA series extracted from the EGM, until each particle has generated a ventricular activation at time t=V̂k,j, i.e., a simulated heartbeat. The resulting qj(V̂k,j) and RTj(V̂k,j) at the time of each new heartbeat V̂k,j is saved. After all N particles in the filter have been used to simulate heartbeats, each V̂k,j is used together with the time of the true heartbeat (Vk) to calculate the weight wk,j of particle j. The weight is related to the probability that V̂k,j−Vk was drawn from a normal distribution with zero mean and standard deviation σw, i.e.,wk,j=NV̂k,j−Vk|0,σw2,where σw was set to 30 ms to account for uncertainties in R wave detection. After all j weights have been calculated, they are further normalized by the sum of all weights. The weighting phase is followed by a resampling phase, where new particles (θ̂k+1,j) with corresponding qj(V̂k,j) and RTj(V̂k,j) are drawn with replacement from θ̂k,j with probability proportional to their weights, thereby approximating the posterior distribution at time step k. In the subsequent propagation phase, each particle is propagated one time step forward by adding normally distributed noise drawn from N(0,Σ). Details on Σ and the propagation phase are found in Supplementary Material S2. The propagation phase is followed by a new weighting phase. The resampling, propagation, and weighting are repeated sequentially for each time step, from k=2 to k=K, where K denotes the last time step. The pseudo-code for the EGM-PF is shown in Algorithm 1.
EGM-particle filter.
Initialization (k = 1):
forj=1 to Ndo
Sample θ̂1,j∼U with exclusion criteria described in Section 2.4.1.
Weighting:
Simulate V̂1,j by running the model with θ̂1,j, qj(0), and RTj(0)=0.
Save qj(V̂1,j) and RTj(V̂1,j) at ventricular activation time V̂1,j.
Compute w1,j=N((V̂1,j−V1)|0,σw) (Equation 4).
end
Normalize w1,j←w1,j/∑jw1,j.
fork=2 to Kdo
forj=1 to Ndo
Resampling:
Generate θ̂k,j by resampling θ̂k−1,j using wk−1,j.
Propagation:
Sample θ̂k,j∼N(θ̂k−1,j,Σ).
Weighting:
Simulate V̂k,j by running the model with θ̂k,j, qj(V̂k−1,j), and RTj(V̂k−1,j).
Save qj(V̂k,j) and RTj(V̂k,j) at ventricular activation time V̂k,j.
Compute wk,j=N((V̂k,j−Vk)|0,σw) (Equation 4).
end
Normalize wk,j←wk,j/∑jwk,j.
end
2.4.2 ECG-particle filter
Similar to the EGM-PF, the ECG-PF is described by its four phases: initialization, weighting, resampling, and propagation. In contrast to the EGM-PF, the atrial activity is not fully known from the ECG. To evaluate several possible AA series, N=40,000 particles are first independently drawn from the same twelve-dimensional uniform distribution as the EGM-PF during initialization, after which NECG=25 copies of each particle are created (θ̂k,j). Additionally, normally distributed noise drawn from N(0,Σ) is added to each particle, with previously defined Σ (Supplementary Material S2). This creates N⋅NECG unique particles, identical to the number of particles in the EGM-PF.
Each unique particle is evaluated with a different AA series (αk,j), with each αk,j generated by a Gaussian random walk, as described in Section 2.2.2. For the first time step, each particle is evaluated by running the model with θ̂1,j until the first ventricular activation has been simulated (V̂1,j).
For the following time steps k>1, the AA series (α) is generated based on a Gaussian random walk (see Section 2.2.2), with an added impulse at time zero. However, the arrival time of the first impulse is set to the arrival time of the atrial impulse leading to V̂k−1,j, denoted ak−1,j, so that αk,j=α+ak−1,j. Moreover, each particle is evaluated by running the model until two ventricular activations have been simulated (V̂k−1,j′ and V̂k,j). Before the time V̂k−1,j, the parameter vector θ̂k−1,j is used to evaluate the model, whereas θ̂k,j is used after V̂k−1,j, as illustrated in Figure 2. Around 98% of the time, V̂k−1,j′ is equal to V̂k−1,j, and particles where V̂k−1,j′ differ from V̂k−1,j are excluded. Consequently, when running the particle filter, V̂k−1,j′ is equivalent to V̂k−1,j. Therefore, impulses in αk−1,j that arrive after ak−1,j do not affect V̂k−1,j in the particle filter and thus do not affect the probability of being selected for the previous resampling. Hence, these impulses should not affect V̂k,j. To ensure this, re-running the previous time step in this manner was performed.
An AA series αk,j and corresponding time series of simulated ventricular activations (V̂k,j′ and V̂k,j), where ak−1,j leads to V̂k,j′. Note that it is not necessarily the first AA impulse after a ventricular activation that leads to the next ventricular activation since impulses may be blocked.
Diagram showing a sequence from \( a_{k-1,j} \) to \( a_{k,j} \), marked by circles and crosses on a line labeled \( \alpha_{k,j} \). Dashed lines indicate transitions. Below, \( \hat{R}R \) is shown with two sections marked \( \hat{V}'_{k-1,j} \) and \( \hat{V}_{k,j} \), labeled "Simulating using \( \theta_{k-1,j} \)" and "Simulating using \( \theta_{k,j} \)" respectively.
After V̂k,j has been generated by the model, values for V̂k,j, qj(V̂k,j), and RTj(V̂k,j) are saved for each particle. As for the EGM-PF, each V̂k,j is used together with a corresponding measured value Vk in Equation 4 to calculate the weight wk,j.
In the resampling phase, N new particles with corresponding qj(ak,j) and RTj(ak,j) are drawn with replacement from θ̂k,j based on their weights, before NECG=25 copies of each particle are created. Normally distributed noise drawn from N(0,Σ) is added to each of the copied particles, which functions as the propagation phase, thereby creating θ̂k+1,j.
The propagation phase is followed by a new weighting phase before the resampling, propagation, and weighting are repeated sequentially for each time step, from k=2 to k=K. The pseudo-code for the ECG-PF is shown in Algorithm 2. In addition, Matlab code with a usage example can be found at https://github.com/FraunhoferChalmersCentre/AV-node-model.
ECG-particle filter. Differences from the EGM-PF are marked with ‘*‘.
Initialization (k = 1):
forj=1 to Ndo
Sample Nθ̂∼U with exclusion criteria described in Section 2.4.1.
* Copy θ̂NECG times to generate θ̂1,j
* Sample θ̂1,j∼N(θ̂1,j,Σ)
Weighting:
* Simulate V̂1,j by running the model with θ̂1,j, qj(0)=α1,j, and RTj(0)=0.
Save qj(a1,j) and RTj(a1,j) at ventricular activation time V̂1,j.
Compute w1,j=N((V̂1,j−V1)|0,σw) (Equation 4).
end
Normalize w1,j←w1,j/∑jw1,j.
fork=2 to Kdo
forj=1 to Ndo
Resampling:
Generate Nθ̂ by resampling θ̂k−1,j using wk−1,j.
* Copy θ̂NECG times to generate θ̂k,j
Propagation:
Sample θ̂k,j∼N(θ̂k,j,Σ).
Weighting:
* Simulate V̂k,j by running the model as described in Section 2.4.2 with θ̂k−1,j, θ̂k,j, qj(ak−1,j), RTj(ak−1,j), and αk,j generated as described in Section 2.2.2.
Compute wk,j=N((V̂k,j−Vk)|0,σw) (Equation 4).
end
Normalize wk,j←wk,j/∑jwk,j.
end
2.4.3 Smoothing algorithm
The combined particle filter and smoothing algorithm utilized in this work is commonly referred to as the forward filtering backward sampling algorithm (Chopin and Papaspiliopoulos, 2020). The smoothing algorithm is applied after either the EGM-PF or the ECG-PF, and functions as the backward sampling step.
Starting at the last time step k=K, one of the j particles is selected with probability proportional to wK,j and denoted x(K), where x is a vector of indices. The algorithm continues iteratively from time step k=K−1 to k=1. The weights are updated based on the likelihood that θ̂k,j originates from the selected particle at the time step of the previous iteration θ̂k+1,x(k+1), with previously defined Σ (see Supplementary Material S2), according to Equation 5.ŵk,j=wk,jNθ̂k+1,xk+1|θ̂k,j,Σ.
A new particle j is selected with probability proportional to ŵk,j and assigned to x(k). After completion, the vector x contains one trajectory of indices corresponding to parameters θ̂k,x(k) sampled from the smoothing probability density function. Running the smoothing algorithm M=20,000 times generates M trajectories (xm) of θ̂, all sampled from the smoothing probability density function. The pseudo-code for the smoothing algorithm is shown in Algorithm 3. Since each simulated heartbeat is associated with a ϕ̂ (see Section 2.3), the M trajectories xm(k) also yield M trajectories of ϕ̂ sampled from the smoothing probability density function. These are used as estimates of the posterior distribution of the AV node conduction delays and refractory periods and denoted ϕ̃m(k).
To evaluate the accuracy of the proposed methodology, estimates obtained using the EGM-PF (ϕ̃mEGM(k)) and ECG-PF (ϕ̃mECG(k)), respectively, are compared with the corresponding ground truth from the simulations. Recordings from the iafdb (see Section 2.1.1) are used to create realistic simulated data, as described in Supplementary Section S1. This results in 50 model parameter trends θ*(k) with associated AA series and f-waves, needed for the EGM-PF and ECG-PF, respectively. Each parameter trend θ*(k) with corresponding AA series was used to simulate an RR series of 1 min duration with associated ground truth AV node conduction property trends (ϕ*(k)).
To compute the most probable value at each time step k for each AV node property, the mode of the estimated trends ϕ̃mEGM(k) and ϕ̃mECG(k) is computed and compared to the mode of the ground truth ϕ*(k). The mode is obtained independently for each time step k by sorting the M values into histogram bins with 5 ms width before identifying the center of the histogram bin with the highest count, resulting in ϕ̃ModeEGM(k) and ϕ̃ModeECG(k). The bin width of 5 ms was chosen to strike a balance between temporal resolution and robustness to noise in the estimated distributions. The average l1-norm (l1̄), between ϕ*(k) and ϕ̃ModeEGM(k) and ϕ̃ModeECG(k), respectively, for each AV node property in each simulation is used to quantify how well the most probable estimate aligns with the ground truth. To facilitate comparison between the AV node conduction properties, l1̄ is also normalized to a percentage by dividing by the range between the highest and lowest values for each property in the simulations ϕ*(k), denoted r (values found in Supplementary Material S1). In addition, ϕ̃ModeEGM(k), ϕ̃ModeECG(k), and ϕ*(k) are averaged over each 1-min simulation, and the l1-norm between the averages, denoted lmin1̄, is used to evaluate how well the proposed methodology performs with a temporal resolution of 1 min.
Furthermore, the percentage of heartbeats for which the 95% credibility region covered ϕ*(k), denoted CR95(k), is calculated by finding the values at the 2.5th and 97.5th percentiles of ϕ̃mEGM(k) and ϕ̃mECG(k) and evaluating how often ϕ*(k) lies in between. In theory, this should converge towards 95% as the number of samples grows.
2.6 Analysis of iafdb data
EGM recordings from the tip of the tricuspid valve annulus with synchronized ECG recordings, as described in Section 2.1.1, are analyzed to compare the concordance between the resulting estimates from the EGM-PF and the ECG-PF. The ϕ̃mEGM(k) and ϕ̃mECG(k) are estimated using the two particle filters and the smoothing algorithm before ϕ̃ModeEGM(k) and ϕ̃ModeECG(k) are calculated as described in Section 2.5. The concordance between ϕ̃ModeEGM(k) and ϕ̃ModeECG(k) is visualized using the Bland–Altman plot.
2.7 Analysis of tilt test data
ECG recordings from a tilt test protocol, as described in Section 2.1.2, are analyzed to evaluate the method’s ability to quantify expected changes in AV node characteristics. First, ϕ̃mECG(k) is estimated using the ECG-PF and smoothing algorithm, before ϕ̃ModeECG(k) is calculated as described in Section 2.5. Further, for each patient, ϕ̃ModeECG(k) is averaged over each tilt phase to obtain ϕ̄Supine, ϕ̄HDT, and ϕ̄HUT, which are used to analyze changes between phases.
2.8 Statistical analysis
The paired one-sided Wilcoxon signed rank test is used to quantify significant increase or decrease in this study, since the data do not generally follow a normal distribution according to the Shapiro-Wilk test (p<0.05). This includes the paired significance test between l1̄, lmin1̄, l1̄/r, and lmin1̄/r obtained using ϕ̃ModeEGM(k) and ϕ̃ModeECG(k), respectively, as well as the paired difference between ϕ̄Supine and ϕ̄HDT and between ϕ̄Supine and ϕ̄HUT.
3 Results
To recapitulate, three datasets were used in this study. The simulated data (see Section 2.5) are used to evaluate the estimation accuracy for the EGM and ECG-based methods. The iafdb recordings (see Section 2.1.1) are used to compare the estimates obtained from synchronized EGM and ECG recordings. The tilt recordings (see Section 2.1.2) are used to evaluate the method’s ability to quantify expected changes in AV node characteristics.
The computation time (performed on a desktop computer with a Ryzen 9 5900X CPU, using the twelve cores in parallel) to obtain posterior distributions ϕ̃mEGM(k) and ϕ̃mECG(k) was on average 3.7 and 17.5 min per minute of data, respectively.
3.1 Evaluation of particle filters
An example of ϕ̃mEGM(k) and ϕ̃mECG(k), estimated from simulated data, with corresponding ground truth ϕ*(k), is displayed in Figure 3. The shown example displays typical patterns, such as the estimates of RSP having a more narrow credibility region compared to the estimates of RFP, DFP, and, DSP, as well as more accurate tracking of changes for DSP estimates (Figure 3d) obtained from EGM recordings compared to ECG recordings. Moreover, neither PF-EGM nor PF-ECG managed to capture fast changes in ϕ*(k) for RFP and DFP (Figures 3a,c).
AV node estimates ϕ̃mEGM(k) (blue) and ϕ̃mECG(k) (red) obtained based on simulated data and corresponding ground truth ϕ*(k) (black): (a)RFP, (b)RSP, (c)DFP, (d)DSP. For comparison with Table 2, EGM and ECG l1̄ are 183 and 160 ms in RFP, 46 and 58 ms in RFP, 114 and 127 ms in RFP, and 76 and 138 ms in RFP. The modes ϕ̃ModeEGM(k) and ϕ̃ModeECG(k) are shown as solid lines and the 80% credibility region as shaded background. The 80% credibility region is used over the 95% for clarity of visualization.
Four line graphs labeled a, b, c, and d compare performances over time steps. Graph a shows \(R^{FP}\); b shows \(R^{SP}\); c shows \(D^{FP}\); and d shows \(D^{SP}\). Lines represent different models: blue for EGM, red for ECG, and black for a reference model, each with shaded confidence intervals. Values range from 0 to 1500 ms on the y-axis.
The results from the analysis of all simulated data are summarized in Table 1. As expected, the estimation error l1̄ in RSP, DFP, and DSP is significantly lower for EGM compared to ECG. However, the opposite is true for RFP. The estimate of RSP has the closest match to the ground truth independently of the type of recording used, with an EGM l1̄/r of 4.48% and ECG l1̄/r of 6.33%. Further, the estimation error of DSP increases the most between EGM and ECG recordings, as seen in an almost doubling of l1̄/r, from 8.97% to 16%. Not surprisingly, the accuracy increases when compressing the 1-min trend into a single estimate, as seen by the decrease in lmin1̄ compared to l1̄ for ECG and EGM for all AV node conduction properties. Additionally, CR95 is typically slightly above 95% for all AV node conduction properties, indicating that the uncertainty bounds produced are conservative if not exact.
The mean ± standard deviation of l1̄ and lmin1̄ between simulated patients, using the EGM-PF and ECG-PF.
Quantification type
RFP
RSP
DFP
DSP
EGM l1̄ (ms)
174 ± 16.5
51.4 ± 11.5
105 ± 8.7
99.5 ± 23.6
ECG l1̄ (ms)
169 ± 13.7a
66.5 ± 9.78b
131 ± 13.3b
178 ± 28.1b
EGM lmin1̄ (ms)
141 ± 35.4
16.2 ± 18.4
24.4 ± 19.5
44.3 ± 33.4
ECG lmin1̄ (ms)
92.3 ± 52.3a
23.6 ± 17.5b
89.6 ± 33.8b
152 ± 65.5b
EGM l1̄/r (%)
12 ± 1.14
4.89 ± 1.1
11 ± 0.916
8.97 ± 2.13
ECG l1̄/r (%)
11.7 ± 0.943a
6.33 ± 0.931b
13.8 ± 1.4b
16 ± 2.53b
EGM lmin1̄/r (%)
9.69 ± 2.44
1.54 ± 1.75
2.57 ± 2.05
3.99 ± 3.01
ECG lmin1̄/r (%)
6.36 ± 3.60a
2.24 ± 1.67b
9.43 ± 3.56b
13.7 ± 5.90b
EGM CR95 (%)
99.5 ± 1.1
99.4 ± 1.39
98.9 ± 1.67
98.8 ± 2.38
ECG CR95 (%)
99.8 ± 0.543
96.4 ± 3.72
98.5 ± 1.97
93 ± 6.48
Indicates a significant decrease for ECG compared to EGM (p<0.05).
Indicates a significant increase for ECG compared to EGM (p<0.05).
3.2 Analysis of iafdb data
An example of ϕ̃mEGM(k) and ϕ̃mECG(k) estimated from synchronized EGM and ECG recordings from one patient is shown in Figure 4. Similar to the results from the analysis of simulated data, the largest difference between using EGM and ECG recordings is seen in DSP (Figure 4d). Notably, the difference between ϕ̃mEGM(k) and ϕ̃mECG(k) is more prominent for the first heartbeats, as seen clearly in the RFP (Figure 4a). This is likely an artifact from the particle filter, which tends to give a higher uncertainty for early time steps.
The AV node estimates ϕ̃mEGM(k)(blue) and ϕ̃mECG(k)(red) obtained for Patient 6 in the iafdb database: (a)RFP, (b)RSP, (c)DFP, (d)DSP. The modes ϕ̃ModeEGM(k)and ϕ̃ModeECG(k)are shown as solid lines and the 80% credibility region as shaded background.
Four line graphs labeled a, b, c, and d show measurements over time steps ranging from zero to fifty-five. Each graph plots two lines: blue for EGM and red for ECG, with shaded areas representing eighty percent confidence intervals. Units are in milliseconds. Graph a measures \(R_{FP}\), b shows \(R_{SP}\), c depicts \(D_{FP}\), and d illustrates \(D_{SP}\).
The results of ϕ̃ModeEGM(k) and ϕ̃ModeECG(k) for all patient is summarized in Table 2, and the Bland–Altman plots of ϕ̃ModeEGM(k) and ϕ̃ModeECG(k) are presented in Figure 5. On average, RSP, DFP, and DSP estimated from ECG are slightly lower than the corresponding estimates obtained from EGM, whereas the opposite is true for RFP. However, the average difference is less than 5% for all AV node conduction properties. Differences between ϕ̃ModeEGM(k) and ϕ̃ModeECG(k) do not seem to be patient-specific.
The mean ± standard deviation of ϕ̃ModeEGM(k) and ϕ̃ModeECG(k) computed over all time steps k for all patients in the iafdb.
Patient
Data
RFP (ms)
RSP (ms)
DFP (ms)
DSP (ms)
1
EGM
917 ± 141
535 ± 104
257 ± 67.8
593 ± 135
1
ECG
993 ± 76.3
492 ± 81.6
231 ± 42.5
627 ± 78.9
2
EGM
818 ± 149
392 ± 69.9
240 ± 77.5
600 ± 132
2
ECG
926 ± 74.8
355 ± 61.3
221 ± 34.4
574 ± 68.2
3
EGM
1,015 ± 108
663 ± 59.7
283 ± 120
761 ± 164
3
ECG
1,037 ± 101
656 ± 79
219 ± 68.8
655 ± 104
4
EGM
1,240 ± 102
1,058 ± 113
355 ± 92.3
507 ± 81.4
4
ECG
1,217 ± 72.2
1,004 ± 37.7
259 ± 93.7
478 ± 67.7
6
EGM
952 ± 251
576 ± 183
253 ± 82
538 ± 159
6
ECG
1,035 ± 89.6
559 ± 161
221 ± 50.8
485 ± 72.9
Bland-Altman plot comparing the concordance between ECG-PF and EGM-PF estimates of (a)RFP, (b)RSP, (c)DFP, and (d)DSPobtained from Patient 1 (blue “o”), Patient 2 (red “+”), Patient 3 (yellow “*”), Patient 4 (black “x”), and Patient 6 (green “□”).
Four Bland-Altman plots labeled a to d compare differences in measurements between EGM and ECG for five patients. Data points are represented by circles, plus signs, asterisks, crosses, and squares, each indicating a different patient. Horizontal lines indicate averages and limits of agreement. Averages are noted as -4.6%, 3%, 4.2%, and 2.8% for plots a, b, c, and d respectively. Plots depict variability across different measurement ranges, with corresponding X and Y axes labeled with mathematical symbols.3.3 Analysis of tilt test data
An example of ϕ̃mECG(k) estimated from one patient during the tilt test is shown in Figure 6. Notably, transient changes in RSP and DSP can be seen in response to HDT and HUT.
The estimated modes ϕ̃ModeECG(k) (lines) and the 80% credibility region (shaded background) for one patient in the tilt test study. Panel (a) show RFP, (b)RSP, (c)DFP, and (d)DSP. Dashed vertical lines indicate the time of tilting, starting in supine position, following five minutes in HDT, before HUT.
Four graphs labeled a, b, c, and d display various measurements in milliseconds over time in minutes, with vertical dashed lines indicating key moments. Graph a) shows \( R_{FP} \), b) shows \( R_{SP} \), c) shows \( D_{FP} \), and d) shows \( D_{SP} \). All graphs include an orange line representing ECG data with a shaded area for ECG 80% confidence.
The resulting average estimates for each patient and each phase, ϕ̄Supine, ϕ̄HDT, and ϕ̄HUT are summarized in Table 3. A significant decrease between ϕ̄Supine and ϕ̄HUT for RFP (p<0.05, 16 of 21 patients), RSP (p<0.001, 18 of 21 patients), and DFP (p<0.01, 16 of 21 patients) could be seen, with no significant difference in DSP. Moreover, no significant difference between ϕ̄Supine and ϕ̄HDT was seen.
The population mean ± standard deviation of the average ϕ̃ModeECG(k) for each tilt position for the patient in the tilt test.
Tilt position
RFP
RSP
DFP
DSP
ϕ̄Supine (ms)
926.2 ± 26.2
349.4 ± 68.8
203.5 ± 13.7
523.8 ± 11.3
ϕ̄HDT (ms)
925.5 ± 24.1
338.7 ± 64.8
202.0 ± 12.4
528.8 ± 12.8
ϕ̄HUT (ms)
921.5 ± 25.3a
327.8 ± 63.4a
199.4 ± 14.2a
525.8 ± 11.7
Indicate a significant decrease (p<0.05) compared with supine position.
4 Discussion
This study proposes a method for estimating AV node conduction properties with beat-to-beat resolution during AF utilizing a network model of the AV node and a particle filter together with a smoothing algorithm. The method was evaluated using simulated data to analyze how well the particle filters estimated the AV node conduction properties for each heartbeat. Additionally, synchronized EGM and ECG recordings were analyzed to study the concordance between estimates obtained from the two data types. Finally, ECG recordings obtained during a tilt test protocol were analyzed to evaluate the method’s ability to quantify expected changes in AV node characteristics.
The proposed framework builds on a physiological model where ϕ must be estimated for each heartbeat based on noisy observations. In this setting, optimization-based techniques such as a genetic algorithms could be used to estimate ϕ beat by beat, but these would only yield point estimates and not the posterior. While Markov chain Monte Carlo methods, such as Metropolis-Hasting, could approximate the full posterior, they are too slow for dynamic beat-to-beat inference. Alternatively, recursive Bayesian methods such as the Kalman filter are computationally efficient but rely on strong assumptions, such as linear dynamics and Gaussian noise, which are not applicable here. Particle filtering and smoothing, by contrast, allow for efficient approximation of the full posterior in a nonlinear, non-Gaussian, and time-varying system, making them well-suited for this application.
The estimation accuracy of the present study can be compared to results from our previous study, where the same AV node model was used to estimate ϕ with a 10-min resolution using an approximate Bayesian computation approach (Karlsson et al., 2024). In that study, ϕ could be estimated with a mean absolute error of 111 ms for RFP, 12 ms for RSP, 90 ms for DFP, and 110 ms for DSP. The present study, by contrast, introduces beat-to-beat tracking of AV nodal dynamics, which represents a major advancement in temporal resolution and allows for the analysis of short-term variations that were not accessible in the earlier framework.
Comparing these results to the ECG l1̄ and ECG lmin1̄ in Table 1, it is evident that the beat-to-beat estimates in the present study have a higher error (ECG l1̄) for RFP, DFP, and DSP, and a drastically higher error for RSP, illustrating a trade-off between temporal resolution and estimation accuracy. Conversely, the error of the 1-min averaged estimates (ECG lmin1̄) is comparable to the error Karlsson et al. (2024). The computation time in the present study (17.5 min per minute of ECG) was multiple magnitudes greater than the 40 s per minute required using the framework of the previous study. Hence, the proposed methodology enables a more detailed analysis of the AV node dynamics at the expense of increased computational cost.
The accuracy of the ECG-based estimates differs between AV node conduction properties, with the lowest error obtained for RSP (6%) and the highest for DSP (16%), as seen in Table 1 (ECG l1̄/r). Comparing the estimated refractory period trends to clinical data is difficult, since the AV node refractory period is assessed using pacing protocols resulting in a point estimate rather than a trend. In addition, pacing protocols assess the refractory period under controlled conditions, not during the irregular conditions of AF, as is the case in this study. Nevertheless, the results in Table 1 suggest that the model and framework cannot identify variations in the RFP and RSP smaller than 169 and 67 ms, respectively.
The differences in estimation accuracy across AV node properties arise from both signal characteristics and model dynamics. Fast pathway properties are estimated with higher error for the simulated data since fewer beats conduct through it, reducing its effect on the output as well as the amount of data available for inference. Moreover, RSP is estimated with similar accuracy using ECG-PF and the EGM-PF, as seen in Table 1, likely because it strongly influences the timing of RR intervals and can therefore be inferred reliably even from ECG recordings. In contrast, the conduction delay is more accurately estimated from EGM recordings, since they provide direct information on atrial activation that is not visible in ECGs.
Only a few sources exist for functional and effective refractory periods in the SP and FP of the AV node in humans. In the study by Denes et al., the functional refractory period for two patients with paroxysmal supraventricular tachycardia and evidence of dual AV nodal conduction was 820 and 495 ms in the FP and 540 and 414 ms in the SP (Denes et al., 1973). In the study by Blanck et al., the effective refractory period in 18 patients with inducible sustained atrioventricular nodal reentrant tachycardia varied from 230 to 440 ms in the FP, and from 180 to 420 ms in the SP (Blanck et al., 1995). Assuming that intra-individual variations in the AV node refractory period are smaller than the difference between individuals, comparing the reported results to the results in Table 1 indicates that identifying beat-to-beat variations in RFP is generally not possible. Instead, the 1-min average, with an error of 92 ms, may be sufficient to detect transient changes in RFP. Moreover, the beat-to-beat error of 67 ms suggests that the model and framework may be sufficient to identify beat-to-beat variations in the RSP. Furthermore, the results in Table 1 suggest that variations in the DFP and DSP smaller than 178 ms cannot be identified. Beat-to-beat variations in AV conduction delay of such large magnitude are generally not observed during normal sinus rhythm (Denes et al., 1973). However, it cannot be excluded that larger variations in conduction delay may occur during AF.
The iafdb was used to study the concordance between estimates obtained from synchronized ECG and EGM recordings. As seen in Figure 4, there are some differences when using EGM recordings compared to ECG recordings. This is clearest in DSP, where estimates obtained using EGM recordings have higher beat-to-beat variability than those obtained using ECG recordings. Thus, the ability to track fast changes is somewhat reduced when only using non-invasive recordings. A difference in convergence time also exists, most clearly visible for RFP in Figure 4. Such convergence behavior was seen for all patients in the iafdb database. Nevertheless, the estimates for all four AV node properties are similar using PF-EGM and PF-ECG, with no average difference larger than 5%.
The AV node conduction properties estimated based on the ECG recordings from the tilt test study, summarized in Table 3, show a significant decrease from supine to HUT position for RFP(p<0.05), RSP(p<0.001), and DFP(p<0.01). Changing the position from supine to HUT is associated with increased sympathetic activity (Aponte-Becerra and Novak, 2021), and hence a decrease in AV nodal conduction delay and refractory period, as previously discussed by Plappert et al. (2022). Thus, the results summarized in Table 3 are in line with what can be expected for HUT, suggesting that the proposed methodology enables the detection of tilt-induced changes in AV node conduction properties. Although the decrease is significant and in line with what is expected for the studied population, the differences are smaller than the corresponding 1-min averaged error (ECG lmin1̄) seen in Table 1, limiting the possibility to detect the changes on an individual level. In addition, the resulting RR interval is influenced not only by the AV node properties, but also by the electrical activity in the atria (Masè et al., 2015; Climent et al., 2010). For the tilt test dataset specifically, a significant increase compared to baseline in the atrial fibrillatory rate was observed during HDT (Östenson et al., 2017). Hence, variations in the RR interval across tilt phases cannot be attributed solely to changes in AV node properties, complicating the analysis.
Moreover, it is less clear how changing the position from supine to HDT affects the sympathetic and parasympathetic activity, and the results from our analysis do not show any significant changes in AV node conduction properties between supine and HDT. Further, as suggested by the estimated trends of RSP and DSP in Figure 6, the proposed methodology may enable the detection of transient changes in the AV node conduction properties.
4.1 Study limitations and future perspectives
The network model accounts for several important properties of the AV node conduction during AF, however, it is by no means a perfect replica. For example, it does not include ventricular escape rhythm and the network topology used in this work excludes some uncommon AV node structures such as multiple slow pathways. Nevertheless, these simplifications are essential to developing a model with a manageable number of parameters. Moreover, given that these are uncommon structures, these limitations are not likely to affect the results.
The estimated AV node properties have only been validated using ground truth data generated from the same AV node model. However, obtaining the exact refractory period and conduction delay in both pathways from patients suffering from AF – if at all possible – would be very difficult and time-consuming. In addition, the number of patient recordings analyzed in this study is limited, which restricts the ability to generalize the findings across the wider AF population.
While a computation time of 17.5 min per minute of analyzed data prevents real-time clinical application of this framework (e.g., using it as part of a monitoring system), this could be addressed in several ways. Since each particle in the particle filter is computed separately, parallellizing the method across CPU cores would reduce processing time. In addition, the number of particles is directly proportional to computation time. In this study, 1,000,000 particles were used to explore the potential of the framework to study beat-to-beat variations. Many applications, such as comparing tilt positions in the tilt test data, do not require beat-to-beat resolution. In those cases, using fewer particles and basing the analysis on 1-min averages could be sufficient. However, to study the transition period between tilt positions, beat-to-beat resolution is necessary, which is only achievable with a time delay using the current implementation of the method. In general, the number of particles, and thereby the computation time and accuracy, needs to be chosen to match the requirements of the clinical application.
The developed methodology constitutes an important step towards non-invasive characterization of autonomic activity in AF patients. This could potentially inform personalized therapeutic strategies and refine prognostic assessments. However, further investigations are required to establish which static and dynamic properties derived from this method exhibit the strongest correlation with clinically significant outcomes. It would thus be interesting to apply the proposed framework to ECG signals obtained during common AF treatments, such as different rate control drugs.
5 Conclusion
We have proposed a novel framework for estimating patient-specific AV node properties with beat-to-beat resolution and conservative uncertainty estimates utilizing a mathematical model combined with a particle filter and smoothing algorithm. By using synchronized EGM and ECG recordings for the parameter estimation, the loss of estimation accuracy using non-invasive recordings could be studied, suggesting a sufficient accuracy for capturing beat-to-beat changes in the refractory period of the SP.
We illustrate the potential of the proposed methodology by analyzing a tilt-test dataset. Results suggest that changes in AV node conduction properties can be assessed from ECG using the proposed method.
Data availability statement
The estimated AV node properties and the simulated data supporting the conclusions for this article will be available from MK upon request. The code for the model together with a user example can be found at https://github.com/FraunhoferChalmersCentre/AV-node-model. The iafdb data are publicly available on https://physionet.org/content/iafdb/1.0.0/. The tilt test protocol data is owned by the Department of Cardiology, Clinical Sciences, Lund University, Sweden, and access requests should be directed to pyotr.platonov@med.lu.se.
Ethics statement
The studies involving humans were approved by Regionala Etikprövningsnämnden in Lund. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and institutional requirements.
Author contributions
MK: Methodology, Investigation, Conceptualization, Software, Visualization, Validation, Writing – original draft, Formal Analysis, Writing – review and editing, Project administration. FP: Data curation, Writing – review and editing. PP: Supervision, Writing – review and editing, Funding acquisition. SÖ: Data curation, Investigation, Writing – review and editing. MW: Methodology, Conceptualization, Resources, Writing – review and editing, Validation, Investigation, Project administration, Supervision, Formal Analysis, Software, Funding acquisition. FS: Supervision, Conceptualization, Investigation, Writing – review and editing, Methodology, Funding acquisition, Project administration, Validation, Formal Analysis, Resources.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported by the Swedish Foundation for Strategic Research (Grant FID18-0023), the Swedish Research Council (Grant VR2019-04272), and the Crafoord Foundation (Grant 20200605).
The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at Tetralith partially funded by the Swedish Research Council through grant agreement no. 2022-06725.
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.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
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.2025.1624403/full#supplementary-material
ReferencesAbdollahpurM.EngströmG.PlatonovP. G.SandbergF. (2022). A subspace projection approach to quantify respiratory variations in the f-wave frequency trend. Front. Physiology13, 976925. 10.3389/fphys.2022.976925Al-KhatibS. M.Allen LaPointeN. M.ChatterjeeR.CrowleyM. J.DupreM. E.KongD. F. (2014). Rate-and rhythm-control therapies in patients with atrial fibrillation: a systematic review. Ann. Intern. Med.160, 760–773. 10.7326/M13-1467AndrewN. E.ThriftA. G.CadilhacD. A. (2013). The prevalence, impact and economic implications of atrial fibrillation in stroke: what progress has been made?Neuroepidemiology40, 227–239. 10.1159/000343667Aponte-BecerraL.NovakP. (2021). Tilt test: a review. J. Clin. Neurophysiology38, 279–286. 10.1097/WNP.0000000000000625BenjaminE. J.MuntnerP.AlonsoA.BittencourtM. S.CallawayC. W.CarsonA. P. (2019). Heart disease and stroke statistics-2019 update a report from the American heart association. Circulation139, e56–e528. 10.1161/CIR.0000000000000659BlanckZ.DhalaA. A.SraJ.DeshpandeS. S.AndersonA. J.AkhtarM. (1995). Characterization of atrioventricular nodal behavior and ventricular response during atrial fibrillation before and after a selective slow-pathway ablation. Circulation91, 1086–1094. 10.1161/01.cir.91.4.1086BotteronG. W.SmithJ. M. (1995). A technique for measurement of the extent of spatial organization of atrial activation during atrial fibrillation in the intact human heart. IEEE Trans. Biomed. Eng.42, 579–586. 10.1109/10.387197BridgeC. P.IoannouC.NobleJ. A. (2017). Automated annotation and quantitative description of ultrasound videos of the fetal heart. Med. image Anal.36, 147–161. 10.1016/j.media.2016.11.006ChopinN.PapaspiliopoulosO. (2020). An introduction to sequential monte carlo, 4. Springer.ClimentA. M.GuillemM. S.HusserD.CastellsF.MilletJ.BollmannA. (2010). Role of the atrial rate as a factor modulating ventricular response during atrial fibrillation. Pacing Clin. Electrophysiol.33, 1510–1517. 10.1111/j.1540-8159.2010.02837.xClimentA. M.GuillemM. S.ZhangY.MilletJ.MazgalevT. (2011). Functional mathematical model of dual pathway AV nodal conduction. Am. J. Physiol. Heart Circ. Physiol.300, H1393–H1401. 10.1152/ajpheart.01175.2010DenesP.WuD.DhingraR. C.ChuquimiaR.RosenK. M. (1973). Demonstration of dual av nodal pathways in patients with paroxysmal supraventricular tachycardia. Circulation48, 549–555. 10.1161/01.cir.48.3.549GoldbergerA. L.AmaralL. A.GlassL.HausdorffJ. M.IvanovP. C.MarkR. G. (2000). Physiobank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals. Circulation101, e215–e220. 10.1161/01.cir.101.23.e215HenrikssonM.CorinoV. D.SörnmoL.SandbergF. (2015). A statistical atrioventricular node model accounting for pathway switching during atrial fibrillation. IEEE Trans. Biomed. Eng.63, 1842–1849. 10.1109/TBME.2015.2503562HenrikssonM.PetrėnasA.MarozasV.SandbergF.SörnmoL. (2018). Model-based assessment of f-wave signal quality in patients with atrial fibrillation. IEEE Trans. Biomed. Eng.65, 2600–2611. 10.1109/TBME.2018.2810508HindricksG.PotparaT.DagresN.ArbeloE.BaxJ. J.Blomström-LundqvistC. (2020). 2020 ESC guidelines for the diagnosis and management of atrial fibrillation developed in collaboration with the European association of cardio-thoracic surgery (EACTS). Am. J. Physiol. Heart Circ. Physiol. 10.1093/eurheartj/ehaa612InadaS.HancoxJ.ZhangH.BoyettM. (2009). One-dimensional mathematical model of the atrioventricular node including atrio-nodal, nodal, and nodal-his cells. Biophysical J.97, 2117–2127. 10.1016/j.bpj.2009.06.056JørgensenP.SchäferC.GuerraP. G.TalajicM.NattelS.GlassL. (2002). A mathematical model of human atrioventricular nodal function incorporating concealed conduction. Bull. Math. Biol.64, 1083–1099. 10.1006/bulm.2002.0313KarlssonM.PlatonovP. G.UlimoenS. R.SandbergF.WallmanM. (2024). Model-based estimation of av-nodal refractory period and conduction delay trends from ECG. Front. Physiology14, 1287365. 10.3389/fphys.2023.1287365KarlssonM.SandbergF.UlimoenS. R.WallmanM. (2021). Non-invasive characterization of human AV-nodal conduction delay and refractory period during atrial fibrillation. Front. Physiol.12, 1849. 10.3389/fphys.2021.728955KarlssonM.WallmanM.PlatonovP. G.UlimoenS. R.SandbergF. (2022). ECG based assessment of circadian variation in AV-nodal conduction during AF – influence of rate control drugs. Front. Physiology2015. 10.3389/fphys.2022.976526KurianT.AmbrosiC.HuckerW.FedorovV. V.EfimovI. R. (2010). Anatomy and electrophysiology of the human av node. Pacing Clin. Electrophysiol.33, 754–762. 10.1111/j.1540-8159.2010.02699.xLeeJ.McManusD. D.BourrellP.SörnmoL.ChonK. H. (2013). Atrial flutter and atrial tachycardia detection using bayesian approach with high resolution time–frequency spectrum from ECG recordings. Biomed. Signal Process. Control8, 992–999. 10.1016/j.bspc.2013.04.002LefflerC. T.SaulJ. P.CohenR. J. (1994). Rate-related and autonomic effects on atrioventricular conduction assessed through beat-to-beat pr interval and cycle length variability. J. Cardiovasc. Electrophysiol.5, 2–15. 10.1111/j.1540-8167.1994.tb01109.xManginL.VinetA.PagéP.GlassL. (2005). Effects of antiarrhythmic drug therapy on atrioventricular nodal function during atrial fibrillation in humans. EP Eur.7, S71–S82. 10.1016/j.eupc.2005.03.016MasèM.MariniM.DisertoriM.RavelliF. (2015). Dynamics of av coupling during human atrial fibrillation: role of atrial rate. Am. J. Physiology-Heart Circulatory Physiology309, H198–H205. 10.1152/ajpheart.00726.2014NathanV.JafariR. (2017). Particle filtering and sensor fusion for robust heart rate monitoring using wearable sensors. IEEE J. Biomed. Health Inf.22, 1834–1846. 10.1109/JBHI.2017.2783758NgJ.SehgalV.NgJ. K.GordonD.GoldbergerJ. J. (2013). Iterative method to detect atrial activations and measure cycle length from electrograms during atrial fibrillation. IEEE Trans. Biomed. Eng.61, 273–278. 10.1109/TBME.2013.2290003ÖstensonS.CorinoV. D.CarlssonJ.PlatonovP. G. (2017). Autonomic influence on atrial fibrillatory process: head-up and head-down tilting. Ann. Noninvasive Electrocardiol.22, e12405. 10.1111/anec.12405PlappertF.WallmanM.AbdollahpurM.PlatonovP. G.ÖstensonS.SandbergF. (2022). An atrioventricular node model incorporating autonomic tone. Front. Physiology13, 1814. 10.3389/fphys.2022.976468RyzhiiM.RyzhiiE. (2023). A compact multi-functional model of the rabbit atrioventricular node with dual pathways. Front. Physiology14, 1126648. 10.3389/fphys.2023.1126648ShafferF.GinsbergJ. P. (2017). An overview of heart rate variability metrics and norms. Front. Public Health5, 258. 10.3389/fpubh.2017.00258ShenM. J.ZipesD. P. (2014). Role of the autonomic nervous system in modulating cardiac arrhythmias. Circulation Res.114, 1004–1021. 10.1161/CIRCRESAHA.113.302549ShkurovichS.SahakianA. V.SwirynS. (1998). Detection of atrial activity from high-voltage leads of implantable ventricular defibrillators using a cancellation technique. IEEE Trans. Biomed. Eng.45, 229–234. 10.1109/10.661270WallmanM.SandbergF. (2018). Characterisation of human AV-nodal properties using a network model. Med. Biol. Eng.56, 247–259. 10.1007/s11517-017-1684-0