Model-based estimation of AV-nodal refractory period and conduction delay trends from ECG

Introduction: Atrial fibrillation (AF) is the most common arrhythmia, associated with significant burdens to patients and the healthcare system. The atrioventricular (AV) node plays a vital role in regulating heart rate during AF by filtering electrical impulses from the atria. However, it is often insufficient in regards to maintaining a healthy heart rate, thus the AV node properties are modified using rate-control drugs. Moreover, treatment selection during permanent AF is currently done empirically. Quantifying individual differences in diurnal and short-term variability of AV-nodal function could aid in personalized treatment selection. Methods: This study presents a novel methodology for estimating the refractory period (RP) and conduction delay (CD) trends, and their uncertainty in the two pathways of the AV node during 24 h using non-invasive data. This was achieved by utilizing a network model together with a problem-specific genetic algorithm and an approximate Bayesian computation algorithm. Diurnal variability in the estimated RP and CD was quantified by the difference between the daytime and nighttime estimates, and short-term variability was quantified by the Kolmogorov-Smirnov distance between adjacent 10-min segments in the 24-h trends. Additionally, the predictive value of the derived parameter trends regarding drug outcome was investigated using several machine learning tools. Results: Holter electrocardiograms from 51 patients with permanent AF during baseline were analyzed, and the predictive power of variations in RP and CD on the resulting heart rate reduction after treatment with four rate control drugs was investigated. Diurnal variability yielded no correlation to treatment outcome, and no prediction of drug outcome was possible using the machine learning tools. However, a correlation between the short-term variability for the RP and CD in the fast pathway and resulting heart rate reduction during treatment with metoprolol (ρ = 0.48, p < 0.005 in RP, ρ = 0.35, p < 0.05 in CD) were found. Discussion: The proposed methodology enables non-invasive estimation of the AV node properties during 24 h, which—indicated by the correlation between the short-term variability and heart rate reduction—may have the potential to assist in treatment selection.


I. INTRODUCTION
Atrial fibrillation (AF) is the most common sustained cardiac arrhythmia and a significant burden for patients and the healthcare system [2].The prevalence of AF is currently estimated to be between 2 and 4% worldwide [3].In addition, the number of AF cases in the European Union is estimated to increase by 89% between 2016 and 2060 [4].Atrial fibrillation is characterized by disorganized electrical activity in the atria, leading to rapid and irregular contraction, and is associated with an increased risk of mortality, predominantly due to heart failure or stroke [5].
The atrioventricular (AV) node acts as the only electrical connection between the atria and ventricles and partly protects the ventricles from the rapid and irregular electrical activity in the atria during AF.It can be functionally divided into two pathways, the fast pathway (FP) and the slow pathway (SP), interconnected at the Bundle of His [6].The AV node either blocks an incoming impulse, based on its refractory period (RP), or sends it through with a delay, based on its conduction delay (CD).The AV node is thus the most essential part in regulating the heart rate during AF, and the RP and CD are the two most important properties of the AV node, deciding its filtering capability.
The AV node during permanent AF is in many cases insufficient in regards to maintaining a healthy heart rate.Therefore, the AV node properties are often modified by treatment with rate control drugs, with β-blockers and calcium channel blockers recommended as first-line treatment [2].Common β-blockers for AF treatment are metoprolol and carvedilol, which block the β1 receptors in the heart in order to reduce the effect of the sympathetic nervous system on the heart [7].Common calcium channel blockers are verapamil and diltiazem, which prevent the L-type calcium channels in the cardiac myocytes from opening in order to reduce conduction in the AV node [8].However, due to the significant and poorly understood individual variability, the choice of drug is currently made empirically for each patient [2].This could lead to a prolonged time until successful treatment, and possibly result in a suboptimal final choice of drug.Since the two recommended first-line treatments have different physiological effects on the AV node, assessing the patient-specific properties of the AV node has the potential to assist in treatment selection.Specifically, we hypothesize that β-blockers would exhibit an increased effect (more reduced heart rate) when variations in the AV node properties are prominent since β-blockers reduce the effect of the sympathetic nervous system.
The AV-node has previously been studied using several mathematical models based on invasive data from humans and animals [9,10,11,12,13,14,15,16].However, in order for a model to be clinically applicable on an individual level, the model parameters should ideally be identifiable from noninvasive data, such as the ECG.A statistical model of the AV node with dual pathway physiology using the RR interval series and the atrial fibrillatory rate (AFR) for model fitting has been proposed [17,18,19].However, the model lumps RP and CD together, limiting its interpretability.
We have previously proposed a network model of the AV node [20] together with a framework for continuously estimating twelve model parameters describing the RP and CD in the two pathways from 24-hour Holter ECG [1].Although promising, the characterization of the AV node was still limited by the number of model parameters and their intrinsic complex dependencies, where a large change in the model parameters could result in a very small change in the RP or CD, thus, making their interpretation a non-trivial task.For a modeling approach to gain acceptance in a clinical context, the outcome should be readily interpretable by medical professionals; a fact that has become especially relevant with the increasing use of advanced modeling and machine learning techniques [21,22].Additionally, in [1], a version of Sobol's method was applied to quantify uncertainty in the parameter estimates.However, these uncertainty estimates were not directly interpretable as probabilities and could thus only be used as a relative measure between the model parameters, between patients, or between different times of the day.When the extent of the uncertainty is unknown, uncertain estimates have the potential to mislead decision-making processes or further analysis of the trends.A proper quantification of the uncertainty is thus advantageous in order to fully understand the estimates.
In the present study, we propose a novel methodology for estimating the RP and CD of both pathways of the AV node and the associated uncertainty continuously over 24 hours.The methodology comprises a genetic algorithm (GA) for initial model parameter estimation and an approximate Bayesian computation (ABC) algorithm to refine the estimates, together with a simulation approach to map model parameters to RP and CD in order to increase interpretability.In addition to refining the estimates, the ABC algorithm provides samples from the Bayesian posterior distribution of the AV node properties, hereafter denoted the posterior, enabling proper quantification of the uncertainty of the estimated properties.We employ these novel tools in an exploratory manner to analyze Holter ECGs from 51 patients during baseline in combination with their respective drug responses to identify potential markers for differences in drug response.Specifically, we analyze the correlation between diurnal and short-term variability and drug outcomes, as well as train several machine learning models to predict drug outcomes.

II. MATERIALS AND METHODS
The overall method for assessing the RP and CD of the two pathways in the AV node for each patient (pat) can be divided into four stages, as shown in Figure 1.Firstly, 24-hour Holter ECGs are processed to extract RR interval series and AFR trends, divided into ten-minute segments (s) with a 50% overlap, as described in Sections II-A and II-B.Secondly, the parameters for the network model of the AV node, described in Section II-C, are fitted to the RR interval series and AFR in each segment using a problem-specific dynamic GA as described in Section II-D1.The GA-derived estimates are subsequently used as inputs to an ABC algorithm to refine the estimates and estimate the posterior of the model parameters, as described in Section II-D2.These model parameter estimates are finally used to simulate data with the model while tracking the RP and CD used for the two pathways, as described in Section II-D3.This results in a distribution of the RP and CD in the FP and the SP for each ten-minute segment.Finally, the possibility to predict treatment outcomes using the estimated distributions is evaluated, as described in Section II-E.

A. ECG Data
Data from the Rate Control in Atrial Fibrillation (RATAF) study, a randomized, investigator-blind, crossover study, approved by the regional ethics committee and the Norwegian Medicines Agency and conducted in accordance with the Helsinki Declaration, is analyzed in this study [23].Specifically, 24-hour ambulatory ECGs from 60 patients (mean age 71 ± 9 years, 18 women) with permanent AF, no heart failure, or symptomatic ischemic heart disease, recorded during baseline, are used for the estimation of patient-specific AV node properties.In addition to the baseline ECG, the relative change in the 24hour average heart rate (∆HR) for treatment with the two Fig. 1: A schematic overview of the methodology, from ECG to estimations of the RP and CD.Previous study refers to [1].calcium channel blockers verapamil and diltiazem and the two β-blockers metoprolol and carvedilol are used to evaluate the therapeutic implications of the estimated AV node properties.The calculation of ∆HR is based on the RR interval series extracted from the ECG, as explained in Section II-B.

B. ECG Processing
The RR interval series is extracted from the ECG for each patient and divided into ten-minute segments with a 50% overlap (RR(pat, s)), where RR intervals following and preceding QRS-complexes with deviating morphology are excluded from the series [24].Segments with excessive noise can lead to a large number of undetected beats and thus an unrealistically low heart rate.Hence, each ten-minute segment is divided into minute-long non-overlapping intervals, and the whole ten-minute segment is excluded from further analysis if any one-minute interval has fewer than 20 detected beats.Patients with RR interval series with a total duration shorter than 12 h are excluded from further analysis.The RR interval series corresponding to the four rate control drugs are calculated equivalently.
Spatiotemporal QRST cancellation is employed to extract the f-waves from the ECG [25].Subsequently, the fundamental frequency of the extracted f-waves is tracked using a hidden Markov model-based method to extract an AFR trend for each patient with a resolution of one minute [26].For time segments where the AFR could not be obtained due to excessive noise, but the RR interval series could, the AFR is set to the closest observed AFR value.

C. Network Model of the AV Node
Our network model of the AV node, introduced in [20], describes the AV node as two pathways (the SP and the FP) comprising 10 nodes each.These two pathways are connected by a coupling node, as illustrated in Figure 2.Each pathway node corresponds physiologically to a localized section of the respective pathway, and the coupling node corresponds to the Purkinje fibers and Bundle of His.
Atrial impulses are modeled by a Poisson process with mean arrival rate λ.The impulses are assumed to reach the first nodes of SP and FP simultaneously.Each network node can be either in a refractory state or in a non-refractory state.A node in its refractory state will block incoming impulses, and a node in its non-refractory state will transmit an incoming impulse to all adjacent nodes with an added conduction delay before entering its refractory state.The RP (R i (n)) and CD (D i (n)) for node i are updated for each incoming impulse n according to Equations 1, 2, and 3, where, ti (n) is the diastolic interval preceding impulse n and t i (n) is the arrival time of impulse n at node i.When ti (n) < 0, the node is in its refractory state and will block incoming impulses.All parameters are fixed for each pathway, resulting in three model parameters for the RP in the FP (R F P min , ∆R F P , τ F P R ); three model parameters for the CD in the FP (D F P min , ∆D F P , τ F P D ); three model parameters for the RP in the SP (R SP min , ∆R SP , τ SP R ); three model parameters for the CD in the SP (D SP min , ∆D SP , τ SP D ).These twelve model parameters constitute the mode parameter vector θ.In addition, the RP in the coupling node is fixed to the mean of the ten shortest RR intervals in the data, and its CD is fixed at 60 ms [20].

D. Parameter Estimation
For each ten-minute segment, the mean arrival rate for the Poisson process λ is estimated as the mean of the AFR trend ( λ(pat, s)), and the model parameters θ(pat, s) are estimated using a GA together with an ABC algorithm.
An error function (ϵ) based on the Poincaré plot, i.e., a scatter plot of successive pairs of RR intervals, is used to quantify the difference between RR(pat, s) and a simulated RR interval series ( RR).The successive pairs of RR intervals for RR(pat, s) and RR are placed in two-dimensional bins covering the interval between 250 and 1800 ms in steps of Fig. 2: 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 [20].50 ms, resulting in K = 961 bins, which we refer to as the Poincaré histogram.The error function, based on the work presented in [20], is computed according to Equation 4, where x k and xk are the numbers of RR intervals in the k-th bin of RR(pat, s) and RR, respectively.Additionally, t norm acts as a normalizing constant and is calculated as the duration of RR divided by the duration of RR(pat, s). 1) Genetic Algorithm: A problem-specific dynamic GA based on the work presented in [1] is used to get an initial estimate of θ(pat, s) in every segment.This results in an estimate denoted as θGA m (pat, s), where m denotes the m-th fittest individual in the population after completion of the GA, i.e. the individual with the m-th lowest ϵ.The hyper-parameters in the algorithm are tuned during the optimization using the difference between the Poincaré histograms in pairs of consecutive segments (∆P ) [1].This difference is calculated using Equation 4 with x k and xk as the number of RR intervals in each bin of the current segment and the following one, respectively.
The GA uses a population of 300 individuals, where each individual is a model parameter vector θ.The algorithm uses tournament selection, a two-point crossover, and creep mutation.To avoid premature convergence and to increase performance, immigration through replacement of the leastfit individuals in the population is performed, following the work in [1].Furthermore, ∆P is used to determine the number of generations that the GA runs before moving to the next data segment, between two and seven.The initialization of individuals is done using latin hypercube sampling within the ranges given in Table I.These values also act as boundaries for the model parameters in the GA.For further details about the algorithm, see [1].
2) Approximate Bayesian Computation: To estimate the posterior p(θ|RR(pat, s), λ(pat, s)), an approximate Bayesian computation population Monte Carlo sampling (ABC PMC) algorithm is used [27].The pseudo-code for the problem-specific ABC PMC is shown in Algorithm 1.The ABC PMC uses a set of N p = 100 particles to estimate the posterior in each RR segment independently, which are updated iteratively for eight iterations (j).Each particle corresponds to a model parameter vector, denoted θABC v,j , where v corresponds to the v-th particle for the j-th iteration.The algorithm is sped up by utilizing the results from the GA to create the initial population.To construct the initial population, twenty particles are drawn from five different normal distributions, N ( , where the covariance matrix Σ init = Cov( θGA 1:25 ) and 1 : 25 denotes [1, 2, ..., 25] for convenience.During each iteration, each particle has a probability of being chosen based on an assigned weight, computed according to Equation 5 [28] where w v,j is the weight for the v-th particle in the j-th iteration and is the probability of θABC k,j−1 given the normal distribution with mean θABC v,j and covariance Σ j−1 , where Σ j = 2Cov( θABC 1:Np,j ).Furthermore, the chosen particle (θ * ) is perturbed to create a proposal particle (θ * * ) using a transition kernel set as N (0, Σ j ) [28].The model is used to simulate data using θ * * to calculate an associated proposal error (ϵ * * ) according to Equation 4. If ϵ * * is lower than a set threshold (T j ), θ * * is accepted and used in the next iteration of the algorithm; if not, a new particle is chosen and perpetuated to create a new proposal particle.Note that the boundaries for the ABC PMC algorithm are more inclusive compared to the GA to accommodate the full width of the estimated posteriors, as shown in Table I.A proposal particle outside the boundaries is always rejected.The next iteration starts when N p new proposal particles have been accepted, and w v,j , T j , and Σ j are then updated.The threshold changes based on the results from the GA; where T 1 = θGA 10 (pat, s), T 2 = θGA 8 (pat, s), T 3 = θGA 5 (pat, s), T 4 = θGA 3 (pat, s), and T 5:8 = θGA 1 (pat, s).Hence, after the eighth iteration, the ϵ for all particles is lower than the ϵ for the fittest individual found by the GA.Thus, the final population is assumed to be N p samples from p(θ|RR(pat, s), λ(pat, s)).
3) Parameter Reduction: The posterior estimate of the parameter vector θ(pat, s) is obtained using the resulting N p samples ( θABC 1:Np,8 (pat, s)) from the ABC PMC algorithm.Each θABC 1:Np,8 (pat, s) is utilized within the model together with the associated λ(pat, s) to simulate a ten-minute RR interval series.For each simulation, R i (n) and D i (n) are stored for each activation n in each pathway node i and used as the sample distribution of the RP and CD for the SP and the FP, respectively.The samples from these four distributions, denoted as Φ(pat, s) = [R F P (pat, s), R SP (pat, s), D F P (pat, s), D SP (pat, s)], serves as a translation from the twelve model parameters θ to four more interpretable AV node properties Φ, taking into account not only the model parameters but also the mean AFR associated with the current RR-segment.To quantify these distributions, their corresponding empirical probability density functions are computed using the MATLAB function ksdensity (MATLAB R2022b) with default bandwidth.From the empirical probability density functions, the maxima are obtained, denoted φmax (pat, s) = [R F P max (pat, s), R SP max (pat, s), D F P max (pat, s), D SP max (pat, s)].In addition, the 5th percentile and the 95th percentile are obtained from Φ(pat, s), denoted φ5 (pat, s) = [R F P 5 (pat, s), R SP 5 (pat, s), D F P 5 (pat, s), D SP 5 (pat, s)], and φ95 (pat, s) = [R F P 95 (pat, s), R SP 95 (pat, s), D F P 95 (pat, s), D SP 95 (pat, s)], respectively.Furthermore, the number of impulses traveling through the FP and SP (N F P and N SP , respectively) is stored, and the ratio is denoted as The patient-specific diurnal variability (∆DV ) in the AV node properties is quantified by the average value of φmax during daytime (9:00 A.M. to 9:00 P.M) divided by the average value of φmax during nighttime (2:00 A.M. to 6:00 A.M).In addition, the patient-specific short-term variability in the AV node properties is quantified by the average Kolmogorov-Smirnov distance (∆KS) between consecutive segments of Φ during the full 24-hour (8:00 A.M to 8:00 A.M).The Kolmogorov-Smirnov distance represents the maximal sep-aration between the empirical cumulative distribution functions between consecutive segments [29].

E. Prediction of Treatment Outcome
The predictive power of the estimates Φ, φ5 , φ95 , φmax , and SP ratio in relation to ∆HR for the different rate control drugs is evaluated in three ways; by analyzing the correlation between the diurnal and short-time variability and ∆HR; by training a feature-based regression model on statistical properties of the trends to predict ∆HR; and by training a convolutional neural network on the trends to predict ∆HR.
To quantify the correlation between diurnal and short-term variability in the AV node properties and ∆HR after treatment with the four rate control drugs, Spearman's rank correlation is used.Due to the exploratory nature of the study, no hypothesis test is performed and hence no correction of p-values is applied [30,31].
Three different feature-based regression models (linear regression, random forest [32], and k-nearest neighbor [33]) are trained on 66 statistical properties of the trends.These statistical properties are; the mean ± std of the four AV node properties φmax during daytime (8 properties), during nighttime (8 properties), and the full 24-hour (8 properties); the mean ± std of the 90% credibility region -calculated as the difference between φ5 and φ95 -during daytime (8 properties),  5) end for Update the covariance for the transition kernel Σ j ← 2Cov( θABC 1:Np,j ) end for nighttime (8 properties), and the full 24-hour (8 properties); the mean ± std of the SP ratio during daytime (2 properties), nighttime (2 properties), and the full 24-hour (2 properties); ∆DV in the four AV node properties (4 properties); the shortterm variability in the four AV node properties (4 properties); as well as the age, gender, weight, and height of the patient.
Deep learning approaches have achieved the current state-of-the-art performance for time-series classification and regression [34].Hence, the prediction of ∆HR for the different rate control drug is evaluated using the time series for φ5 , φ95 , φmax , SP ratio , AFR, and the RR interval series as an input to three convolutional neural networks with different architectures, based on only fully connected layers [35], the ResNet architecture [35], and the Inception architecture [36], respectively.To incorporate the age, gender, weight, and height of the patients, the last fully connected layer of the networks is modified to also include these properties as input neurons.The networks were trained using the tsai library [37], with the Adam solver [38] and the Huber loss [39].Leave-one-out cross-validation is used, so that the network is trained on data from all but one patient and tested on the left-out patient.The average mean square error (MSE) of the predicted and true ∆HR for the whole population is calculated and compared between approaches.

III. RESULTS
As described in Section II-A, this study is based on a population of 60 patients.However, due to excessive noise, some patients are excluded from analysis, as described in Section II-B, resulting in a total of 51 patients.In addition, excessive noise in the ECG during treatment with the four rate control drugs leads to missing values for ∆HR for some patients.Thus, of the remaining 51 patients at baseline, two lack data for verapamil, three lack data for diltiazem, two lack data for metoprolol, none lack data for carvedilol, and one lacks data for both verapamil and metoprolol.The mean ± standard deviation of ∆HR in the population are 19% ± 23% for verapamil; 24% ± 18% for diltiazem, 17% ± 18% for metoprolol; and 11% ± 6% for carvedilol.

A. Parameter Trends
The 24-hour trends of φmax (pat, s), φ5 (pat, s), and φ95 (pat, s) for two patients, denoted A and B, are presented in Figure 3 and 4. Figure 3 shows a low short-term variability in the RP and CD in both pathways for patient A (∆KS = [0.27,0.19, 0.24, 0.33] for R F P , R SP , D F P , and, D SP ), whereas patient B in Figure 4 has a larger shortterm variability (∆KS = [0.41,0.55, 0.40, 0.40]).Conduction mainly occurs through the SP in both patients, as indicated by an SP ratio over 0.5, which results in a wider credibility region in the R F P compared to the R SP .However, for patient B, there are segments where the FP is more prevalent, e.g. between 5 PM and 6 PM.In these segments, the RP and CD have a very low variability indicating a stationary behavior of the AV node.A notable shift in RP occurs at 8 AM for patient A, probably as a response to waking up from sleep, resulting in a clear change in autonomic regulation.No notable difference between the average R F P , R SP , and D F P during daytime and during nighttime could be seen for patient A, with a slight difference in D SP (∆DV = [0.80,0.81, 0.99, 1.39]).For patient B, only D F P showed a notable difference (∆DV = [0.81,0.92, 2.60, 1.19]).
Similar observations can be made for the whole population, as displayed in Table II, which includes the mean and standard deviation of φmax (pat, s), the 95% credibility region, and ∆KS, during daytime, nighttime, and during 24 hours, as well as ∆DV , for the RP and CD in the FP and the SP for all patients.For convenience, the total CD, calculated by multiplying the CD for one node by ten, is listed.From Table II, it is evident that the RP on average is higher and the CD is lower during nighttime compared to daytime, probably linked Fig. 3: The estimated RF (top) and CD (middle) for φmax (pat, s) (dotted) as well as φ5 (pat, s) and φ95 (pat, s) (filled) for the FP (blue) and SP (red), as well as the SP ratio (bottom) are shown for patient A, marked with a black circle in Figure 6.
to the lower heart rate during sleep and/or circadian autonomic variations.Figure 5 illustrates the population average trends of φmax (pat, s), φ5 (pat, s), and φ95 (pat, s).To reduce the influence of outliers, only segments containing data from over 20% of the population are shown.A distinct separation between RP and CD of the two pathways exists, indicating different functionality.Additionally, the credibility region for the R F P is larger compared to the R SP .Moreover, the credibility region for D F P , in proportion to its mean value, is larger than that of D SP .The differences in credibility regions between FP and SP reflect the SP ratio , which is 0.78 ± 0.11 (mean ± std) during the day, 0.79 ± 0.12 during the night, and 0.78 ± 0.10 during the full 24-hour, indicating that the SP is dominant on average.

B. Prediction of Treatment Outcome
Spearman's rank correlation between the patient-specific ∆DV , as described in Section II-E, and ∆HR showed no clear correlation (p < 0.05) for any combination of drug and AV node property.Hence, no relationship between diurnal variability and drug outcome was found.
The Spearman correlation between the patient-specific shorttime variability, quantified by ∆KS, and ∆HR showed no clear correlation (p < 0.05) for the RP and CD in the SP.A moderate correlation was however found between ∆KS and ∆HR for R F P in the β-blocker metoprolol (ρ = 0.47, p = 0.0011) and for D F P in metoprolol (ρ = 0.35, p = 0.017).Figure 6 shows the individual ∆KS plotted against ∆HR and their linear relation for all four drugs, with the left panel showing R F P and the right panel showing D F P .Interestingly, a similar relation between ∆KS, and ∆HR is not present in the other β-blocker carvedilol.
The ability to predict ∆HR using machine learning approaches is evaluated by the average MSE between the predicted and true ∆HR for the four drugs using the leaveone-out validation method.The average MSE is benchmarked against the population variance of ∆HR for the four drugs.Hence, if the average MSE is larger than the population variance at 0.0071%, the population mean yields a more accurate predictor.Using the feature-based regression models, as described in Section II-E, resulted in an average MSE of 0.0073% for the linear regression, an average MSE of 0.0074% for the random forest, and an average MSE of 0.074% for the k-nearest neighbor.In addition, using the convolutional neural network resulted in an average MSE of 0.0073% for the fully connected architecture, an average MSE of 0.0079% for the ResNet architecture, and an average MSE of 0.0074% for the Inception architecture.Overall, all the machine-learning approaches resulted in an average MSE higher than 0.0071% and thus in a poor fit to new-seen data.

IV. DISCUSSION
A mathematical model with an associated framework for patientspecific estimation and proper uncertainty quantification of the RP and CD in the FP and SP of the AV node using only non-invasive data has been proposed.
Individual estimation of trends and variability in AV node properties using non-invasive data has the potential to increase the patient-specific understanding of the AV node during AF, which in turn can be used to enhance informatics approaches for the next generation of personalized medicine.The two most dominant properties of the AV node, the RP and CD, together with the ratio of impulses conducted through the different pathways, have the potential to increase the understanding of the AV node and its function during AF.
Due to the physiological differences between the effect of β-blockers and calcium channel blockers, where β-blockers reduce the effect of the sympathetic nervous system, our hypothesis was that β-blockers could exhibit an increased effect when variations in the AV node properties are prominent since this would indicate a larger influence of the autonomic nervous system.The population-averaged trends (Figure 5) Fig. 4: The estimated RF (top) and CD (middle) for φmax (pat, s) (dotted) as well as φ5 (pat, s) and φ95 (pat, s) (filled) for the FP (blue) and SP (red), together with the SP ratio (bottom) are shown for patient B, marked with a red circle in Figure 6.
show an increase in RP and a slight decrease in CD during nighttime compared to daytime, suggesting that the decreased sympathetic activity during nighttime affects the RP and CD.However, no correlation was found between diurnal variations in AV properties and reduction in heart rate during treatment with β-blockers.
However, a potential association between the shorttime variability and the treatment outcome with metoprolol was found.The findings depicted in Figure 6 demonstrate a moderate correlation between ∆KS and the change in heart rate (∆HR) in the RP and CD for the FP for metoprolol, but not for any other drugs.The lack of correlation between ∆HR after treatment with carvedilol (also a β-blocker) and ∆KS could potentially be attributed to its modest overall effect observed in the RATAF study, likely stemming from its rapid elimination as acknowledged in [40].For the possible relationship between short-time variation in the RP and CD in the FP between metoprolol and treatment outcome suggested by this analysis, additional studies are needed to confirm the results.
The possibility to predict ∆HR for the different rate control drugs used in this study was evaluated using three featuredbased regression models and three different architectures of a convolutional neural network (Section III-B).With a resulting average MSE higher than the variance of ∆HR for the population, it appears impossible to predict ∆HR with any certainty in the present data set.Either there is not enough information relevant for predicting the heart rate reduction after drug treatment in the AV node property trends -possibly due to the 10-minute resolution, limiting the information about autonomic regulation -or the data set size of 51 patients is too low given the inter-individual variability present in the measurements.
Prior iterations of the model and framework focused on estimating the model parameter trends rather than the patient-specific property trends of the AV node [1].This approach imposed limitations on the interpretability of the results, since the interpretation of the model parameters in terms of common cardiology terminology such as RP and CD is not straightforward.In contrast, the current work introduces a novel methodology that enables the estimation of the RP and CD for each ECG segment individually, facilitating a more comprehensible and interpretable analysis.The ability to derive such estimates is vital as it allows for effective communication of the analysis results.Furthermore, this advancement in methodology opens up new avenues for gaining a deeper understanding of the AV node and its diurnal and short-term variations.
The estimation of the posterior by obtaining a range of plausible values, as opposed to relying on a point estimate of the AV node properties, offers notable advantages.For example, the credibility region for R F P in Figure 4 is very broad during most segments at nighttime, reflecting a high uncertainty.In scenarios where the extent of the uncertainty is unknown, these uncertain estimates have the potential to influence decisionmaking processes or further analysis of the trends.As a result, the usefulness and reliability of these estimates may be decreased, emphasizing the need for a posterior estimation approach.
It has previously been shown that a GA can estimate the time-varying network model parameters during 24 hours [1].However, in order to estimate the posterior distribution, the ABC approach was necessary.The ABC approach has in recent years been used for the personalization of the electrophysiological properties in cardiac models [41].Although ABC approaches are generally computationally expensive [27], starting in a promising area of the model parameter space, derived from the GA results, reduced the computation time by a factor of around 50 (data not shown).The GA was also used to decide on a reasonable threshold level for the ABC PMC algorithm, which is not straightforward since imperfections in the model make certain RR series more challenging to replicate Fig. 5: The average RF (top) and CD (middle) for φmax (pat, s) (dotted) as well as φ5 (pat, s) and φ95 (pat, s) (filled) for the FP (blue) and SP (red), together with the mean (black, dotted) and standard deviation (black, filled) of the SP ratio (bottom).than others, resulting in a higher average ϵ.Hence, an ϵ value corresponding to a good fit for one RR interval series could correspond to a poor fit for another, making thresholds very data-dependent.Using the GA to find the threshold levels ensures a reasonable threshold level specified for each data segment.

A. Study Limitations and Future Perspectives
The estimated RP and CD have not been validated against intracardiac measurements, since obtaining such measurements during AF -if at all possible -would be very difficult and time-consuming.The average RP and CD for the two pathways can however be compared with invasive electrophysiological measurements of the AV node from two patients with paroxysmal supraventricular tachycardia and evidence of dual AV nodal conduction found in the literature [42].The two patients had an RP in the FP of 820 ms and 495 ms; an RP in the SP of 540 ms and 414 ms; a CD in the FP of 125 ms and 150 ms; and a CD in the SP of 500 ms and 300 ms.Comparing these values to the daytime estimates seen in Table II, it is evident that the measured values for the RP and CD in both pathways are within the range of our estimated values.It should be noted that the measured functional RP values come from an S1-S2 protocol during sinus rhythm, thus the comparison is not trivial.The functional RP is the smallest AA-interval preceding a conducted impulse.It is however still dependent on the previous pacing frequency, which is not well-defined during AF.Nevertheless, since AF leads to high frequencies, the RP should be reasonably close to the functional RP.In addition, the estimated CD from our model and framework shown in Table II corresponds to the peak of the probability density function of all CDs in each pathway multiplied by 10.Hence, it differs slightly from the measured CD, since it also captures CDs corresponding to impulses that are blocked within the node.
In this study, short-time variability was estimated as the difference between adjacent 10-minute intervals.However, limiting the short-time variability to ten minutes also limits the information about the autonomic nervous system -which is known to operate on a higher resolution -to a ten-minute resolution.Hence, improving the time resolution of the analysis has the possibility to increase the information extracted by the model and framework, which could improve the results.Furthermore, to extract even more information about the impact of the autonomic nervous system on the AV node, an extension of the model has been proposed in [43].A similar framework to the one presented in this work could be employed for that model to estimate model parameters and simulate the RP and CD.This could further refine the estimates and thus the information about the AV node.
Moreover, analyzing the RP and CD trends for all the patients, a high inter-individual variability with a wide range of diurnal and short-time variability could be seen, likely due to the inherent individual differences.This, in combination with the relatively low number of patients (51), indicates that the results in this paper should be verified in a larger study.

V. CONCLUSION
We have proposed a novel framework for estimating patientspecific 24-hour trends of the RP and CD in the FP and SP of the AV node by mapping estimated model parameters.These estimates include the full posterior of the RP and CD and could be estimated using only non-invasive data.Additionally, a correlation between short-term variability in both the RP and CD for the FP and drug-induced changes to the heart rate was found.The individual estimates of AV node properties offer patient-specific trends in RP and CD, which may have the potential to assist in treatment selection.Fig. 6: Scatter plot of the 24-hour ∆HR and ∆KS for the R F P (left) and D F P (right) for the four drugs, with patient A (as shown in Figure 3) marked with black and patient B (as shown in Figure 4) marked with red.

TABLE I :
Parameter ranges for the GA and the ABC PMC algorithm.