ECG-based beat-to-beat assessment of AV node conduction properties during AF

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 ($\pm$ std) of 169$\pm$14 ms for the refractory period in the fast pathway; 131$\pm$13 ms for the conduction delay in the fast pathway; 67$\pm$10 ms for the refractory period in the slow pathway; and 178$\pm$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.

1 Introduction

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 $\mathit{\varphi }=$ [$R^{FP}$, $R^{SP}$, $D^{FP}$, and $D^{SP}$] – 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 $\mathit{\varphi }$.

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 $\mathit{\varphi }$ 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 $\mathit{\varphi }$ 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 ($\mathit{\varphi }$ 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 $\mathit{\varphi }$ 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 $\mathit{\varphi }$ (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 $\pm$ 10 years, 60% male, and an average signal duration of 58 $\pm$ 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 $\pm$ 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 $\mathit{\alpha }$, is generated by a Gaussian random walk described by a mean $({\mu }^{\alpha })$ and standard deviation $({\sigma }^{\alpha })$. Both ${\mu }^{\alpha }$ and ${\sigma }^{\alpha }$ 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 $({\mu }^f)$ and standard deviation $({\sigma }^f)$ of the inverse f-wave frequency are calculated, as well as the $SQI$. Subsequently, ${\mu }^{\alpha }$ is drawn from $\mathcal{N}({\mu }^f,\text{max}{(\mathrm{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, ${\sigma }^{\alpha }$ is set to $4{\sigma }^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).

Figure 1

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

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 $(R_i(n))$ and conduction delay $(D_i(n))$ for node $i$ are updated for each incoming impulse $n$ according to Equations 1–3,

$R_i\left(n\right)=R_{\mathit{min}}+\mathrm{\Delta }R\left(1-e^{-{\tilde{t}}_i\left(n\right)/{\tau }_R}\right)(1)$

$D_i\left(n\right)=D_{\mathit{min}}+\mathrm{\Delta }D{e}^{-{\tilde{t}}_i\left(n\right)/{\tau }_D},(2)$

${\tilde{t}}_i\left(n\right)=t_i\left(n\right)-\left(t_i\left(n-1\right)+R_i\left(n-1\right)\right),(3)$

where ${\tilde{t}}_i(n)$ is the diastolic interval preceding impulse $n$ and $t_i(n)$ is the arrival time of impulse $n$ at node $i$. When ${\tilde{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 $\mathit{\theta }=[R_{\mathit{min}}^{FP},\mathrm{\Delta }{R}^{FP},{\tau }_R^{FP},R_{\mathit{min}}^{SP},\mathrm{\Delta }{R}^{SP},{\tau }_R^{SP},D_{\mathit{min}}^{FP},\mathrm{\Delta }{D}^{FP},{\tau }_D^{FP},D_{\mathit{min}}^{SP},\mathrm{\Delta }{D}^{SP},$${\tau }_D^{SP}]$. 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 $t_i(n)$ in a queue $(\mathit{q}(t))$ of impulses is propagated next (or blocked depending on ${\tilde{t}}_i(n)$). Three types of data are needed to run the model. First, the parameter vector $\mathit{\theta }$ is necessary, corresponding to the properties of the AV node. Second, the queue $\mathit{q}(t)$ is required, where each impulse is represented by a tuple containing arrival time $t_i(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 $\mathit{q}(t)$ with corresponding node index for the first nodes in the slow and fast pathway. Finally, a vector containing the repolarization times $(\mathit{R}\mathit{T}(t))$ for each of the 21 nodes in the model is needed, corresponding to the end of the diastolic interval $({\tilde{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 $\mathit{\theta }$ are assumed to be fixed between two heartbeats and the corresponding AV node conduction properties $\mathit{\varphi }$ are estimated in the following way: for a given pathway, e.g., FP, $R^{FP}$ and $D^{FP}$ are estimated as the medians of all of $R_i(n)$ and $D_i(n)$ (Equations 1, 2) in the time interval, with $D^{FP}$ 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 $\mathit{\varphi }$ 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 ${\widehat{\mathit{\theta }}}_{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=\mathrm{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 ${\widehat{\mathit{\theta }}}_{1,j}$ with the model, using the current repolarization times $({\mathit{R}\mathit{T}}_j(0)=0)$ and the current queue $({\mathit{q}}_j(0))$ filled by the AA series extracted from the EGM, until each particle has generated a ventricular activation at time $t={\widehat{V}}_{k,j}$, i.e., a simulated heartbeat. The resulting ${\mathit{q}}_j({\widehat{V}}_{k,j})$ and ${\mathit{R}\mathit{T}}_j({\widehat{V}}_{k,j})$ at the time of each new heartbeat ${\widehat{V}}_{k,j}$ is saved. After all $N$ particles in the filter have been used to simulate heartbeats, each ${\widehat{V}}_{k,j}$ is used together with the time of the true heartbeat $(V_k)$ to calculate the weight $w_{k,j}$ of particle $j$. The weight is related to the probability that ${\widehat{V}}_{k,j}-V_k$ was drawn from a normal distribution with zero mean and standard deviation ${\sigma }_w$, i.e.,

$w_{k,j}=\mathcal{N}\left(\left({\widehat{V}}_{k,j}-V_k\right)|0,{\sigma }_w^2\right),(4)$

where ${\sigma }_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 $({\widehat{\mathit{\theta }}}_{k+1,j})$ with corresponding ${\mathit{q}}_j({\widehat{V}}_{k,j})$ and ${\mathit{R}\mathit{T}}_j({\widehat{V}}_{k,j})$ are drawn with replacement from ${\widehat{\mathit{\theta }}}_{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 $\mathcal{N}(0,\mathrm{\Sigma })$. Details on $\mathrm{\Sigma }$ 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.

Algorithm 1

Algorithm 1. EGM-particle filter.

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=\mathrm{40,000}$ particles are first independently drawn from the same twelve-dimensional uniform distribution as the EGM-PF during initialization, after which $N_{\mathit{ECG}}=25$ copies of each particle are created $({\widehat{\mathit{\theta }}}_{k,j})$. Additionally, normally distributed noise drawn from $\mathcal{N}(0,\mathrm{\Sigma })$ is added to each particle, with previously defined $\mathrm{\Sigma }$ (Supplementary Material S2). This creates $N\cdot N_{\mathit{ECG}}$ unique particles, identical to the number of particles in the EGM-PF.

Each unique particle is evaluated with a different AA series $({\mathit{\alpha }}_{k,j})$, with each ${\mathit{\alpha }}_{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 ${\widehat{\mathit{\theta }}}_{1,j}$ until the first ventricular activation has been simulated $({\widehat{V}}_{1,j})$.

For the following time steps $k>1$, the AA series $(\mathit{\alpha })$ 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 ${\widehat{V}}_{k-1,j}$, denoted $a_{k-1,j}$, so that ${\mathit{\alpha }}_{k,j}=\mathit{\alpha }+a_{k-1,j}$. Moreover, each particle is evaluated by running the model until two ventricular activations have been simulated (${\widehat{V}}_{k-1,j}^{\prime }$ and ${\widehat{V}}_{k,j}$). Before the time ${\widehat{V}}_{k-1,j}$, the parameter vector ${\widehat{\mathit{\theta }}}_{k-1,j}$ is used to evaluate the model, whereas ${\widehat{\mathit{\theta }}}_{k,j}$ is used after ${\widehat{V}}_{k-1,j}$, as illustrated in Figure 2. Around 98% of the time, ${\widehat{V}}_{k-1,j}^{\prime }$ is equal to ${\widehat{V}}_{k-1,j}$, and particles where ${\widehat{V}}_{k-1,j}^{\prime }$ differ from ${\widehat{V}}_{k-1,j}$ are excluded. Consequently, when running the particle filter, ${\widehat{V}}_{k-1,j}^{\prime }$ is equivalent to ${\widehat{V}}_{k-1,j}$. Therefore, impulses in ${\mathit{\alpha }}_{k-1,j}$ that arrive after $a_{k-1,j}$ do not affect ${\widehat{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 ${\widehat{V}}_{k,j}$. To ensure this, re-running the previous time step in this manner was performed.

Figure 2

Figure 2. An AA series ${\mathit{\alpha }}_{k,j}$ and corresponding time series of simulated ventricular activations (${\widehat{V}}_{k,j}^{\prime }$ and ${\widehat{V}}_{k,j}$), where $a_{k-1,j}$ leads to ${\widehat{V}}_{k,j}^{\prime }$. 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.

After ${\widehat{V}}_{k,j}$ has been generated by the model, values for ${\widehat{V}}_{k,j}$, ${\mathit{q}}_j({\widehat{V}}_{k,j})$, and ${\mathit{R}\mathit{T}}_j({\widehat{V}}_{k,j})$ are saved for each particle. As for the EGM-PF, each ${\widehat{V}}_{k,j}$ is used together with a corresponding measured value $V_k$ in Equation 4 to calculate the weight $w_{k,j}$.

In the resampling phase, $N$ new particles with corresponding ${\mathit{q}}_j({\mathit{a}}_{k,j})$ and ${\mathit{R}\mathit{T}}_j({\mathit{a}}_{k,j})$ are drawn with replacement from ${\widehat{\mathit{\theta }}}_{k,j}$ based on their weights, before $N_{\mathit{ECG}}=25$ copies of each particle are created. Normally distributed noise drawn from $\mathcal{N}(0,\mathrm{\Sigma })$ is added to each of the copied particles, which functions as the propagation phase, thereby creating ${\widehat{\mathit{\theta }}}_{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.

Algorithm 2

Algorithm 2. ECG-particle filter. Differences from the EGM-PF are marked with ‘*‘.

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 $w_{K,j}$ and denoted $\mathit{x}(K)$, where $\mathit{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 ${\widehat{\mathit{\theta }}}_{k,j}$ originates from the selected particle at the time step of the previous iteration ${\widehat{\mathit{\theta }}}_{k+1,\mathit{x}(k+1)}$, with previously defined $\mathrm{\Sigma }$ (see Supplementary Material S2), according to Equation 5.

${\widehat{w}}_{k,j}=w_{k,j}\mathcal{N}\left({\widehat{\mathit{\theta }}}_{k+1,\mathit{x}\left(k+1\right)}|{\widehat{\mathit{\theta }}}_{k,j},\mathrm{\Sigma }\right).(5)$

A new particle $j$ is selected with probability proportional to ${\widehat{w}}_{k,j}$ and assigned to $\mathit{x}(k)$. After completion, the vector $\mathit{x}$ contains one trajectory of indices corresponding to parameters ${\widehat{\mathit{\theta }}}_{k,x(k)}$ sampled from the smoothing probability density function. Running the smoothing algorithm $M=\mathrm{20,000}$ times generates $M$ trajectories $({\mathit{x}}_m)$ of $\widehat{\mathit{\theta }}$, 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 $\widehat{\mathit{\varphi }}$ (see Section 2.3), the $M$ trajectories ${\mathit{x}}_m(k)$ also yield $M$ trajectories of $\widehat{\mathit{\varphi }}$ 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 ${\tilde{\mathit{\varphi }}}_m(k)$.

Algorithm 3

Algorithm 3. Smoothing algorithm.

2.5 Evaluation of particle filters

To evaluate the accuracy of the proposed methodology, estimates obtained using the EGM-PF $({\tilde{\mathit{\varphi }}}_m^{\mathit{EGM}}(k))$ and ECG-PF $({\tilde{\mathit{\varphi }}}_m^{\mathit{ECG}}(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 ${\mathit{\theta }}^{*}(k)$ with associated AA series and f-waves, needed for the EGM-PF and ECG-PF, respectively. Each parameter trend ${\mathit{\theta }}^{*}(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 $({\mathit{\varphi }}^{*}(k))$.

To compute the most probable value at each time step $k$ for each AV node property, the mode of the estimated trends ${\tilde{\mathit{\varphi }}}_m^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_m^{\mathit{ECG}}(k)$ is computed and compared to the mode of the ground truth ${\mathit{\varphi }}^{*}(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 ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(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 $l^1$-norm $(\overline{l^1})$, between ${\mathit{\varphi }}^{*}(k)$ and ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(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, $\overline{l^1}$ is also normalized to a percentage by dividing by the range between the highest and lowest values for each property in the simulations ${\mathit{\varphi }}^{*}(k)$, denoted $\mathit{r}$ (values found in Supplementary Material S1). In addition, ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{EGM}}(k)$, ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(k)$, and ${\mathit{\varphi }}^{*}(k)$ are averaged over each 1-min simulation, and the $l^1$-norm between the averages, denoted $\overline{l_{\mathit{min}}^1}$, 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 ${\mathit{\varphi }}^{*}(k)$, denoted $C{R}_{95}(k)$, is calculated by finding the values at the 2.5th and 97.5th percentiles of ${\tilde{\mathit{\varphi }}}_m^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_m^{\mathit{ECG}}(k)$ and evaluating how often ${\mathit{\varphi }}^{*}(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 ${\tilde{\mathit{\varphi }}}_m^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_m^{\mathit{ECG}}(k)$ are estimated using the two particle filters and the smoothing algorithm before ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(k)$ are calculated as described in Section 2.5. The concordance between ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(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, ${\tilde{\mathit{\varphi }}}_m^{\mathit{ECG}}(k)$ is estimated using the ECG-PF and smoothing algorithm, before ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(k)$ is calculated as described in Section 2.5. Further, for each patient, ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(k)$ is averaged over each tilt phase to obtain ${\bar{\mathit{\varphi }}}^{\mathit{Supine}}$, ${\bar{\mathit{\varphi }}}^{\mathit{HDT}}$, and ${\bar{\mathit{\varphi }}}^{\mathit{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 $\overline{l^1}$, $\overline{l_{\mathit{min}}^1}$, $\overline{l^1}/\mathit{r}$, and $\overline{l_{\mathit{min}}^1}/\mathit{r}$ obtained using ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(k)$, respectively, as well as the paired difference between ${\bar{\mathit{\varphi }}}^{\mathit{Supine}}$ and ${\bar{\mathit{\varphi }}}^{\mathit{HDT}}$ and between ${\bar{\mathit{\varphi }}}^{\mathit{Supine}}$ and ${\bar{\mathit{\varphi }}}^{\mathit{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 ${\tilde{\mathit{\varphi }}}_m^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_m^{\mathit{ECG}}(k)$ was on average 3.7 and 17.5 min per minute of data, respectively.

3.1 Evaluation of particle filters

An example of ${\tilde{\mathit{\varphi }}}_m^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_m^{\mathit{ECG}}(k)$, estimated from simulated data, with corresponding ground truth ${\mathit{\varphi }}^{*}(k)$, is displayed in Figure 3. The shown example displays typical patterns, such as the estimates of $R^{SP}$ having a more narrow credibility region compared to the estimates of $R^{FP}$, $D^{FP}$, and, $D^{SP}$, as well as more accurate tracking of changes for $D^{SP}$ estimates (Figure 3d) obtained from EGM recordings compared to ECG recordings. Moreover, neither PF-EGM nor PF-ECG managed to capture fast changes in ${\mathit{\varphi }}^{*}(k)$ for $R^{FP}$ and $D^{FP}$ (Figures 3a,c).

Figure 3

Figure 3. AV node estimates ${\tilde{\mathit{\varphi }}}_m^{\mathit{EGM}}(k)$ (blue) and ${\tilde{\mathit{\varphi }}}_m^{\mathit{ECG}}(k)$ (red) obtained based on simulated data and corresponding ground truth ${\mathit{\varphi }}^{*}(k)$ (black): (a) $R^{FP}$, (b) $R^{SP}$, (c) $D^{FP}$, (d) $D^{SP}$. For comparison with Table 2, EGM and ECG $\overline{l^1}$ are 183 and 160 ms in $R^{FP}$, 46 and 58 ms in $R^{FP}$, 114 and 127 ms in $R^{FP}$, and 76 and 138 ms in $R^{FP}$. The modes ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(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.

The results from the analysis of all simulated data are summarized in Table 1. As expected, the estimation error $\overline{l^1}$ in $R^{SP}$, $D^{FP}$, and $D^{SP}$ is significantly lower for EGM compared to ECG. However, the opposite is true for $R^{FP}$. The estimate of $R^{SP}$ has the closest match to the ground truth independently of the type of recording used, with an EGM $\overline{l^1}/\mathit{r}$ of 4.48% and ECG $\overline{l^1}/\mathit{r}$ of 6.33%. Further, the estimation error of $D^{SP}$ increases the most between EGM and ECG recordings, as seen in an almost doubling of $\overline{l^1}/\mathit{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 $\overline{l_{\mathit{min}}^1}$ compared to $\overline{l^1}$ for ECG and EGM for all AV node conduction properties. Additionally, $C{R}_{95}$ is typically slightly above 95% for all AV node conduction properties, indicating that the uncertainty bounds produced are conservative if not exact.

Table 1

Table 1. The mean $\pm$ standard deviation of $\overline{l^1}$ and $\overline{l_{\mathit{min}}^1}$ between simulated patients, using the EGM-PF and ECG-PF.

3.2 Analysis of iafdb data

An example of ${\tilde{\mathit{\varphi }}}_m^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_m^{\mathit{ECG}}(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 $D^{SP}$ (Figure 4d). Notably, the difference between ${\tilde{\mathit{\varphi }}}_m^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_m^{\mathit{ECG}}(k)$ is more prominent for the first heartbeats, as seen clearly in the $R^{FP}$ (Figure 4a). This is likely an artifact from the particle filter, which tends to give a higher uncertainty for early time steps.

Figure 4

Figure 4. The AV node estimates ${\tilde{\mathit{\varphi }}}_m^{\mathit{EGM}}(k)$(blue) and ${\tilde{\mathit{\varphi }}}_m^{\mathit{ECG}}(k)$(red) obtained for Patient 6 in the iafdb database: (a) $R^{FP}$, (b) $R^{SP}$, (c) $D^{FP}$, (d) $D^{SP}$. The modes ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{EGM}}(k)$and ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(k)$are shown as solid lines and the 80% credibility region as shaded background.

The results of ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(k)$ for all patient is summarized in Table 2, and the Bland–Altman plots of ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(k)$ are presented in Figure 5. On average, $R^{SP}$, $D^{FP}$, and $D^{SP}$ estimated from ECG are slightly lower than the corresponding estimates obtained from EGM, whereas the opposite is true for $R^{FP}$. However, the average difference is less than 5% for all AV node conduction properties. Differences between ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(k)$ do not seem to be patient-specific.

Table 2

Table 2. The mean $\pm$ standard deviation of ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{EGM}}(k)$ and ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(k)$ computed over all time steps $k$ for all patients in the iafdb.

Figure 5

Figure 5. Bland-Altman plot comparing the concordance between ECG-PF and EGM-PF estimates of (a) $R^{FP}$, (b) $R^{SP}$, (c) $D^{FP}$, and (d) $D^{SP}$obtained from Patient 1 (blue “o”), Patient 2 (red “+”), Patient 3 (yellow “*”), Patient 4 (black “x”), and Patient 6 (green “$\square$”).

3.3 Analysis of tilt test data

An example of ${\tilde{\mathit{\varphi }}}_m^{\mathit{ECG}}(k)$ estimated from one patient during the tilt test is shown in Figure 6. Notably, transient changes in $R^{SP}$ and $D^{SP}$ can be seen in response to HDT and HUT.

Figure 6

Figure 6. The estimated modes ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(k)$ (lines) and the 80% credibility region (shaded background) for one patient in the tilt test study. Panel (a) show $R^{FP}$, (b) $R^{SP}$, (c) $D^{FP}$, and (d) $D^{SP}$. Dashed vertical lines indicate the time of tilting, starting in supine position, following five minutes in HDT, before HUT.

The resulting average estimates for each patient and each phase, ${\bar{\mathit{\varphi }}}^{\mathit{Supine}}$, ${\bar{\mathit{\varphi }}}^{\mathit{HDT}}$, and ${\bar{\mathit{\varphi }}}^{\mathit{HUT}}$ are summarized in Table 3. A significant decrease between ${\bar{\mathit{\varphi }}}^{\mathit{Supine}}$ and ${\bar{\mathit{\varphi }}}^{\mathit{HUT}}$ for $R^{FP}$ ($p<0.05$, 16 of 21 patients), $R^{SP}$ ($p<0.001$, 18 of 21 patients), and $D^{FP}$ ($p<0.01$, 16 of 21 patients) could be seen, with no significant difference in $D^{SP}$. Moreover, no significant difference between ${\bar{\mathit{\varphi }}}^{\mathit{Supine}}$ and ${\bar{\mathit{\varphi }}}^{\mathit{HDT}}$ was seen.

Table 3

Table 3. The population mean $\pm$ standard deviation of the average ${\tilde{\mathit{\varphi }}}_{\mathit{Mode}}^{\mathit{ECG}}(k)$ for each tilt position for the patient in the tilt test.

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 $\mathit{\varphi }$ 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 $\mathit{\varphi }$ 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 $\varphi$ with a 10-min resolution using an approximate Bayesian computation approach (Karlsson et al., 2024). In that study, $\varphi$ could be estimated with a mean absolute error of 111 ms for $R^{FP}$, 12 ms for $R^{SP}$, 90 ms for $D^{FP}$, and 110 ms for $D^{SP}$. 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 $\overline{l^1}$ and ECG $\overline{l_{\mathit{min}}^1}$ in Table 1, it is evident that the beat-to-beat estimates in the present study have a higher error (ECG $\overline{l^1}$) for $R^{FP}$, $D^{FP}$, and $D^{SP}$, and a drastically higher error for $R^{SP}$, illustrating a trade-off between temporal resolution and estimation accuracy. Conversely, the error of the 1-min averaged estimates (ECG $\overline{l_{\mathit{min}}^1}$) 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 $R^{SP}$ (6%) and the highest for $D^{SP}$ (16%), as seen in Table 1 (ECG $\overline{l^1}/\mathit{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 $R^{FP}$ and $R^{SP}$ 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, $R^{SP}$ 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 $R^{FP}$ is generally not possible. Instead, the 1-min average, with an error of 92 ms, may be sufficient to detect transient changes in $R^{FP}$. 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 $R^{SP}$. Furthermore, the results in Table 1 suggest that variations in the $D^{FP}$ and $D^{SP}$ 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 $D^{SP}$, 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 $R^{FP}$ 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 $R^{FP}$ $(p<0.05)$, $R^{SP}$ $(p<0.001)$, and $D^{FP}$ $(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 $\overline{l_{\mathit{min}}^1}$) 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 $R^{SP}$ and $D^{SP}$ 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 cHlvdHIucGxhdG9ub3ZAbWVkLmx1LnNl.

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

Acknowledgments

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

References

Abdollahpur M., Engström G., Platonov P. G., Sandberg F. (2022). A subspace projection approach to quantify respiratory variations in the f-wave frequency trend. Front. Physiology 13, 976925. doi:10.3389/fphys.2022.976925

PubMed Abstract | CrossRef Full Text | Google Scholar

Al-Khatib S. M., Allen LaPointe N. M., Chatterjee R., Crowley M. J., Dupre M. E., Kong D. F., et al. (2014). Rate-and rhythm-control therapies in patients with atrial fibrillation: a systematic review. Ann. Intern. Med. 160, 760–773. doi:10.7326/M13-1467

PubMed Abstract | CrossRef Full Text | Google Scholar

Andrew N. E., Thrift A. G., Cadilhac D. A. (2013). The prevalence, impact and economic implications of atrial fibrillation in stroke: what progress has been made? Neuroepidemiology 40, 227–239. doi:10.1159/000343667

PubMed Abstract | CrossRef Full Text | Google Scholar

Aponte-Becerra L., Novak P. (2021). Tilt test: a review. J. Clin. Neurophysiology 38, 279–286. doi:10.1097/WNP.0000000000000625

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamin E. J., Muntner P., Alonso A., Bittencourt M. S., Callaway C. W., Carson A. P., et al. (2019). Heart disease and stroke statistics-2019 update a report from the American heart association. Circulation 139, e56–e528. doi:10.1161/CIR.0000000000000659

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanck Z., Dhala A. A., Sra J., Deshpande S. S., Anderson A. J., Akhtar M., et al. (1995). Characterization of atrioventricular nodal behavior and ventricular response during atrial fibrillation before and after a selective slow-pathway ablation. Circulation 91, 1086–1094. doi:10.1161/01.cir.91.4.1086

PubMed Abstract | CrossRef Full Text | Google Scholar

Botteron G. W., Smith J. 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. doi:10.1109/10.387197

PubMed Abstract | CrossRef Full Text | Google Scholar

Bridge C. P., Ioannou C., Noble J. A. (2017). Automated annotation and quantitative description of ultrasound videos of the fetal heart. Med. image Anal. 36, 147–161. doi:10.1016/j.media.2016.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Chopin N., Papaspiliopoulos O. (2020). An introduction to sequential monte carlo, 4. Springer.

Google Scholar

Climent A. M., Guillem M. S., Husser D., Castells F., Millet J., Bollmann A. (2010). Role of the atrial rate as a factor modulating ventricular response during atrial fibrillation. Pacing Clin. Electrophysiol. 33, 1510–1517. doi:10.1111/j.1540-8159.2010.02837.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Climent A. M., Guillem M. S., Zhang Y., Millet J., Mazgalev T. (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

Denes P., Wu D., Dhingra R. C., Chuquimia R., Rosen K. M. (1973). Demonstration of dual av 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

Goldberger A. L., Amaral L. A., Glass L., Hausdorff J. M., Ivanov P. C., Mark R. G., et al. (2000). Physiobank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals. Circulation 101, e215–e220. doi:10.1161/01.cir.101.23.e215

PubMed Abstract | CrossRef Full Text | Google Scholar

Henriksson M., Corino V. D., Sörnmo L., Sandberg F. (2015). A statistical atrioventricular node model accounting for pathway switching during atrial fibrillation. IEEE Trans. Biomed. Eng. 63, 1842–1849. doi:10.1109/TBME.2015.2503562

PubMed Abstract | CrossRef Full Text | Google Scholar

Henriksson M., Petrėnas A., Marozas V., Sandberg F., Sörnmo L. (2018). Model-based assessment of f-wave signal quality in patients with atrial fibrillation. IEEE Trans. Biomed. Eng. 65, 2600–2611. doi:10.1109/TBME.2018.2810508

PubMed Abstract | CrossRef Full Text | Google Scholar

Hindricks G., Potpara T., Dagres N., Arbelo E., Bax J. J., Blomström-Lundqvist C., et al. (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. doi:10.1093/eurheartj/ehaa612

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Jørgensen P., Schäfer C., Guerra P. G., Talajic M., Nattel S., Glass L. (2002). A mathematical model of human atrioventricular nodal function incorporating concealed conduction. Bull. Math. Biol. 64, 1083–1099. doi:10.1006/bulm.2002.0313

PubMed Abstract | CrossRef Full Text | Google Scholar

Karlsson M., Platonov P. G., Ulimoen S. R., Sandberg F., Wallman M. (2024). Model-based estimation of av-nodal refractory period and conduction delay trends from ECG. Front. Physiology 14, 1287365. doi:10.3389/fphys.2023.1287365

PubMed Abstract | CrossRef Full Text | Google Scholar

Karlsson M., Sandberg F., Ulimoen S. R., Wallman M. (2021). Non-invasive characterization of human AV-nodal conduction delay and refractory period during atrial fibrillation. Front. Physiol. 12, 1849. doi:10.3389/fphys.2021.728955

PubMed Abstract | CrossRef Full Text | Google Scholar

Karlsson M., Wallman M., Platonov P. G., Ulimoen S. R., Sandberg F. (2022). ECG based assessment of circadian variation in AV-nodal conduction during AF – influence of rate control drugs. Front. Physiology 2015. doi:10.3389/fphys.2022.976526

PubMed Abstract | CrossRef Full Text | Google Scholar

Kurian T., Ambrosi C., Hucker W., Fedorov V. V., Efimov I. R. (2010). Anatomy and electrophysiology of the human av node. Pacing Clin. Electrophysiol. 33, 754–762. doi:10.1111/j.1540-8159.2010.02699.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee J., McManus D. D., Bourrell P., Sörnmo L., Chon K. H. (2013). Atrial flutter and atrial tachycardia detection using bayesian approach with high resolution time–frequency spectrum from ECG recordings. Biomed. Signal Process. Control 8, 992–999. doi:10.1016/j.bspc.2013.04.002

CrossRef Full Text | Google Scholar

Leffler C. T., Saul J. P., Cohen R. 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. doi:10.1111/j.1540-8167.1994.tb01109.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Mangin L., Vinet A., Pagé P., Glass L. (2005). Effects of antiarrhythmic drug therapy on atrioventricular nodal function during atrial fibrillation in humans. EP Eur. 7, S71–S82. doi:10.1016/j.eupc.2005.03.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Masè M., Marini M., Disertori M., Ravelli F. (2015). Dynamics of av coupling during human atrial fibrillation: role of atrial rate. Am. J. Physiology-Heart Circulatory Physiology 309, H198–H205. doi:10.1152/ajpheart.00726.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Nathan V., Jafari R. (2017). Particle filtering and sensor fusion for robust heart rate monitoring using wearable sensors. IEEE J. Biomed. Health Inf. 22, 1834–1846. doi:10.1109/JBHI.2017.2783758

PubMed Abstract | CrossRef Full Text | Google Scholar

Ng J., Sehgal V., Ng J. K., Gordon D., Goldberger J. J. (2013). Iterative method to detect atrial activations and measure cycle length from electrograms during atrial fibrillation. IEEE Trans. Biomed. Eng. 61, 273–278. doi:10.1109/TBME.2013.2290003

PubMed Abstract | CrossRef Full Text | Google Scholar

Östenson S., Corino V. D., Carlsson J., Platonov P. G. (2017). Autonomic influence on atrial fibrillatory process: head-up and head-down tilting. Ann. Noninvasive Electrocardiol. 22, e12405. doi:10.1111/anec.12405

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryzhii M., Ryzhii E. (2023). A compact multi-functional model of the rabbit atrioventricular node with dual pathways. Front. Physiology 14, 1126648. doi:10.3389/fphys.2023.1126648

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaffer F., Ginsberg J. P. (2017). An overview of heart rate variability metrics and norms. Front. Public Health 5, 258. doi:10.3389/fpubh.2017.00258

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen M. J., Zipes D. P. (2014). Role of the autonomic nervous system in modulating cardiac arrhythmias. Circulation Res. 114, 1004–1021. doi:10.1161/CIRCRESAHA.113.302549

PubMed Abstract | CrossRef Full Text | Google Scholar

Shkurovich S., Sahakian A. V., Swiryn S. (1998). Detection of atrial activity from high-voltage leads of implantable ventricular defibrillators using a cancellation technique. IEEE Trans. Biomed. Eng. 45, 229–234. doi:10.1109/10.661270

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar