Heart rate characteristic based modelling of atrial fibrillatory rate using implanted cardiac monitor data

Objective. The objective of the present study is to investigate the feasibility of using heart rate characteristics to estimate atrial fibrillatory rate (AFR) in a cohort of atrial fibrillation (AF) patients continuously monitored with an implantable cardiac monitor. We will use a mixed model approach to investigate population effect and patient specific effects of heart rate characteristics on AFR, and will correct for the effect of previous ablations, episode duration, and onset date and time. Approach. The f-wave signals, from which AFR is estimated, were extracted using a QRST cancellation process of the AF episodes in a cohort of 99 patients (67% male; 57 ± 12 years) monitored for 9.2(0.2–24.3) months as median(min-max). The AFR from 2453 f-wave signals included in the analysis was estimated using a model-based approach. The association between AFR and heart rate characteristics, prior ablations, and episode-related features were modelled using fixed-effect and mixed-effect modelling approaches. Main results. The mixed-effect models had a better fit to the data than fixed-effect models showing h.c. of determination (R2 = 0.49 versus R2 = 0.04) when relating the variations of AFR to the heart rate features. However, when correcting for the other factors, the mixed-effect model showed the best fit (R2 = 0.04). AFR was found to be significantly affected by previous catheter ablations (p < 0.05), episode duration (p < 0.05), and irregularity of the RR interval series (p < 0.05). Significance. Mixed-effect models are more suitable for AFR modelling. AFR was shown to be faster in episodes with longer duration, less organized RR intervals and after several ablation procedures.


Introduction
As the general population ages, the incidence of atrial fibrillation (AF) rises (Murray et al 2016). Nowadays, AF has become the most common arrhythmia encountered in clinical practices leading to increased mortality. However, the underlying mechanisms are still under investigation and appropriate patient selection for treatment still remains a challenge (Rottner et al 2020, Saglietto et al 2021. The opportunity to find additional ways of characterizing atrial electromechanical and anatomical properties and ways of predicting subsequent outcome after therapeutic intervention would favour timely therapy selection. More specifically, AF is related to a compromised atrial function caused by a fast and irregular atrial depolarization which can be characterized from the f-waves in the ECG. There have been many parameters introduced in the literature to characterize f-waves including their amplitude, frequency, morphology Sörnmo 2001, Corino et al 2007) and, atrial organization (Alcaraz and José 2007).
The f-wave frequency, often referred to as the atrial fibrillatory rate (AFR) is an AF characteristic which has been subject to considerable clinical attention. It has been shown to be a useful tool for monitoring drug effects (Platonov et al 2014) as well as for predicting outcomes from clinical procedures such as successful AF cardioversion (Bollmann et al 2003) and early AF recurrence (Bollmann et al 2008).
The correlation of AFR and several well-known HRV features describing the variability and irregularity of RR intervals during AF has been previously explored. While the analysis showed that variability parameters were independent from AFR, the irregularity parameters were significantly correlated with AFR (Corino et al 2013). However, the extrapolation of their findings for all clinical types of AF remains to be determined as the population used in the study included only patients with underlying congestive heart failure. In addition, the analysis did not account for variations on AFR due to circadian cycles, episode duration or previous ablations.
Circadian variations in the AF frequency within a 24 h period have been studied using Holter recordings in pursuit of understanding the underlying mechanisms of AF. It was concluded that AFR was significantly lower during night-time than during daytime (Bollmann et al 2000, Cosson et al 2000, Meurling et al 2001. These studies, however, have a drawback as AFR was computed from sparse measurements with several hours in between estimates. In addition, in one of the studies, two different sets of patients were identified: one which showed an increase (minority) while the other showed a decrease (majority) in nocturnal AFR (Bollmann et al 2000). A later study used more advanced signal processing techniques and obtained a more robust estiamates and a continuous AF frequency trend (Sandberg et al 2010). This study found that circadian variations were present in most of the patients with long-standing persistent AF analysed (13/18). However, the former and the latter studies have the drawback of having their insight on atrial electrophysiological characteristics during AF constrained to 24 h long Holter registrations, limited datasets of up to 30 patients and were all based on persistent and chronic AF patients.
Other studies have shown a positive correlation between episode duration and AFR (R = 0.53, p < 0.05) (Bollmann et al 1999). However, this study was conducted in a small dataset with only 31 episodes from 11 paroxysmal AF patients.
Implantable cardiac monitors (ICMs) with highly sensitive AF detection algorithms showing detection rates up to 96% (Hindricks et al 2010) are widely used in the clinical practice for the diagnosis and monitoring of AF. Clinical applications of these monitors include AF diagnosis after catheter ablation (Kapa et al 2013, Verma et al 2013 and cryptogenic stroke (Etgen et al 2013). These devices offer the advantage of long-term monitoring which can lead to a more detailed characterization of the AFR behaviour. These long monitoring periods spanning several months, thus including several episodes per patient, allow an analysis of the joint effects of HRV derived features and, circadian variation, episode duration and previous ablations on AFR.
The estimation of AFR from RR series data would enable wearable based assessment of AFR, e.g. wristband PPG, which would lead to a better characterization of the patient's condition. The aim of this study is to model variations in AFR using RR series features, by correcting for the effect of time of the day of episode onset, episode duration and previous ablations on AFR, in a cohort of AF patients continuously monitored with an ICM. The AFR was estimated using a model-based approach for f-wave analysis (section 2.2) that allowed for a robust estimation of AFR from ICM data. In a previous study regarding AFR and HRV, the analysis was conducted only on patients with underlying congestive heart failure and did not account for the presence of confounding factors (Corino et al 2013). We used a simple fixed-effect (FE) modelling approach using RR series features and compared it to two more complex mixed-effect (ME) modelling approaches: one to study both the population and patient specific effects of RR series in AFR and another that allowed correction for the effect of episode duration, previous ablations, and possible circadian variations. In addition, ME modelling will account for the heterogeneity within AF patients.

Patient population
The present study population consisted of a subset of the cohorts of patients with symptomatic AF enrolled in the Reveal LINQ Usability study, a multi-centre single-arm clinical study (ClinicalTrials.gov Identifier: NCT01965899) (Sanders et al 2016), and a Slovak tertiary-care arrhythmia centre database (Bou Ezzeddine et al 2015). In both cohorts, the patients provided written informed consent, and the study protocols were reviewed and approved by the Human Research Ethics Committee of each participating institution in accordance with the Declaration of Helsinki. 99 patients (67% male, 33% female; 57 ± 12 years) had pre-and post-ablation monitoring and were included in the analysis. Nineteen patients had a previous failed ablation before the ablation procedure considered in this study. The clinical baseline characteristics of the analysed patients are shown in table 1.
The devices used in the Usability and the Slovakia studies were the Reveal LINQ and Reveal XT (Medtronic Inc., Minneapolis, MN), respectively which were implanted within the fourth intercostal space (V2-V3 electrode orientation) near the apex of the heart. The feasibility of extracting atrial activation from ICM data has been previously explored . Both devices record the signal with sampling frequency of 256 Hz and then sense and detect the rhythm by analysing RR interval patterns to compute an AF evidence score every 2 min. Due to memory restrictions, the devices store a single-lead ECG signal of 2 min of the AF episode detected as well as the ventricular sense, i.e. the positions of the R-peaks. In addition, the devices store the detected episode onset date and time, and the total duration of the episode.
The 99 patients included in the study had the ICM implanted and were followed-up for 9.2 (0.2-24.3) months as median (min-max). The ablation procedure was performed 5.8 (1.0-14.4) months after the implant. The ablation procedures were either pulmonary vein isolation (PVI) only (76 patients, 77%) or PVI plus extra lesions, which included roof and mitral lines, and ablation of complex fractionated atrial electrograms (23 patients, 23%). AF recurrence was defined as an AF episode detected by the ICM after a 3 month blanking period following catheter ablation and only those episodes outside the blanking period were considered in the analysis. The blanking period is based on reports describing how early recurrences could be caused by post-ablation inflammation or short-term autonomic imbalance rather than ablation failure (Hindricks et al 2021). In the analysed cohort, 31 (31%) had AF recurrence, 38 (38%) had no AF recurrence, and 30 patients (30%) left the study before the ending of the 3 months blanking period so there is no available information of their recurrence status. To evaluate the variations of AFR, the episodes occurring during the full monitoring period, except the 3 months blanking period, of the 99 patients included in the study were considered.

Estimation of AFR
The atrial activity, from which the AFR is computed, is extracted from ECG signals using Cardiolund ECG Parser (www.cardiolund.com).
Once the atrial activity (x(n)) having a length of N samples is obtained, a harmonic f-wave model (Henriksson et al 2018) is used to estimate the local f-wave frequency. In this model, the f-waves are modelled by a complex signal s(n, θ), defined as the sum of two harmonically related, complex exponentials with fundamental frequency f ( ) where A m and f m define the amplitude and phase, respectively, of the m-th exponential (first and second), f s is the sampling frequency and θ, the parameter vector.
The complex valued analytic representation x a (n) of the observed f-wave signal x(n), is assumed to be composed of ( ) q s n, and white, complex Gaussian noise e is a 2 × 1 vector containing the amplitude and phase information: and ( ) w Z 0 is anŃ 2 Vandermonde matrix defined as: The model is evaluated in 5 s windows by locally fitting the model in K 0.5 s overlapping sub-segments and estimatingq using a maximum likelihood approach. For each subsegment, the local frequency estimateŵ k 0, is determined by: The local frequency estimates are then averaged over the 5 s windows and the AFR of the segment is determined as:ˆˆ( The model's fit is then evaluated using the model errorˆ( ) ( ) (ˆ) q = e n x n s n; a to estimate the signal quality index (SQI) The SQI is confined between [0, 1] with larger values associated to a better fit. A fixed threshold is used to indicate whether f-waves have sufficient quality for the analysis. For this study, the SQI was estimated every 5 s windows and the threshold was set to 0.3 as suggested by Henriksson et al (2018). Figure 1 shows the f-wave extraction process for the AFR estimation.
For further details on the estimation of AFR and the SQI, the reader is referred to Henriksson et al (2018). For each episode, the local f-wave frequency calculated for each 5 s segment is then averaged over the estimates of sufficient signal quality, thus obtaining the AFR. Figure 1. Illustration of f-wave extraction for AFR estimation. From top to bottom: ECG signal extracted from ICM, QRST-cancelled signal (x(n)), signal quality index (solid line) with threshold for acceptable signal quality (dashed line) and, estimated frequency (ˆ) f and estimated f-wave signal ( (ˆ) q s n; ) with signal segment with signal quality below the acceptable threshold marked as light grey.
The mean AFR of at least one acceptable 5 s segment was considered to be representative of the first minutes of the AF episode under the assumption that the AFR was stable within 2 min of AF. The stability of AFR within the episodes was studied by selecting the episodes where more than 80% of the episode had acceptable levels of SQI and iteratively computing the mean AFR for decreasing percentages of the signal. In each iteration, an additional random 5 s segment was left out and the relative absolute error between the mean AFR of the reduced signal segments and the mean AFR of the whole acceptable signal was evaluated.
2.3. Modelling FE models are statistical models which only contain fixed effects. In contrast, ME models contain both fixed effects and random effects, and are useful when dealing with data involving multiple sources of random error such as repeated measures within subjects (Stroup 2013).
In general terms, for the ME model in this analysis, we consider N patients, with the index i representing the ith patient (   i N 1 ). Each patient has n i measurements y , ij with the index j representing the jth episode (   j n 1 i ). P then is defined as the total number of episodes included in the analysis so that There are M random effects considered, and the patient-specific random effect is represented by ). The y ij of each episode j in each patient i is assumed to follow a Gaussian distribution: is the conditional expectation of the observations, in this case the AFR of each episode y , ij given by random vector b , i containing the patient specific random effects [ and s 2 is the dispersion parameter of the distribution.
Linear ME models can be represented as: where y is the set of y , ij X is the design matrix for fixed effects, i.e. not dependent on the patient, Z is the design matrix for random effects, b is a vector of fixed effect parameters and e contains the random errors associated with the ijth observation. The design matrix ( ) is aṔ M matrix which contains the values for each of the M features describing the episodes considered in the analysis. Similarly, the design matrix for the random effects ( ⁎( ) ) contains the values z bm ij which represent the value of feature m describing the jth episode of the ith patient.
The maximum pseudo likelihood method is employed for the parameter estimation. For ME models with G N 0, , being R and G, variance matrices for distributions of | y b and b respectively, the log-likelihood equation can be written as (Stroup 2013): The estimator is then given by: see (Stroup 2013) for details on the maximum pseudo likelihood method. FE models can be represented as in equation (10) but with the random effects coefficient b being equal to 0. Hence, the log-likelihood equation can be written as: and the estimator is given by: In this study, a linear FE model is used to evaluate the effects of changes in RR series characteristics on AFR. In addition, the results obtained will be compared to two ME models which will evaluate both the fixed and random effects of circadian variations, previous ablations, episode duration and RR series characteristics on AFR.
Following the previous formulations, y ij corresponds to the AFR ij of the ijth episode, X and Z contain the features comprised by time, number of ablations, duration of episode and heart rate characteristic features, and b and b contain the fixed and random effects of the model.
The FE and ME models (described bellow) were fitted and evaluated using Matlab R2022a (The Mathworks Inc., Natick, Massachusetts). The Akaike Information Criterion (AIC) (Akaike 1974) and the deviance residual, i.e. an index of model fit, where a model with a higher deviance provides a poorer model fit to the data than a model with a lower deviance (Zeng et al 2014), are used to select the model which fits best the dataset. In addition, to check the goodness of fit of the models, the fitted values of AFR are compared to the observed values of AFR, and the coefficients of determination (R 2 ) (Wright 1921) are computed.
2.3.1. Fixed-effect model of AFR From the ventricular sense of the AF episode stored in the ICM, heart rate characteristic features are computed to represent the variability and irregularity of the RR intervals in each episode. For this study, the parameters extracted were the mean RR intervals (mean interval between ventricular senses in milliseconds), the mean squared differences of successive RR intervals (RMSSD) expressed in milliseconds and calculated as: where RR i is the ith RR interval, + RR i 1 the successive interval, and N is the total number of RR intervals, and the sample entropy (SampEn). SampEn estimates the irregularity of the RR series and was computed as described in (Richman and Moorman 2000). To compute SampEn, the length of the vectors to be compared in the calculation (m) and the noise rejection level (r) must be selected. In this study, the chosen parameters were = m 2 and s = r 0.2 , where s is the standard deviation of the RR intervals. The FE model assumes fixed (population) effects of RR mean, RR variability and RR irregularity on AFR (Mean RR, RMSSD and SampEn) hence an estimate of AFR (ÂFR ij ), for patient i and episode j, is given by: where b 0 is the intercept estimate for the fixed effect. In a similar way, b , 1 b 2 and b 3 represent the fixed effects of the Mean RR, RMSSD and SampEn respectively, in Model FE.

Mixed-effect model of AFR
In order to account for the heterogeneity between patients, an ME model is used. The ME model assumes both the fixed and random effects of changes in RR series on AFR by introducing the patient-specific random effects to equation (18): where b i 0 represents the intercept estimate for the patient-specific random effect, and b , i 1 b i 2 and b i 3 represent the patient-specific random effect of Mean RR, RMSSD and SampEn respectively.

Mixed-effect model of AFR and correct for confounding factors
There are several variables that could affect AFR and need to be accounted for. In this study, the effect of circadian variations quantified by the episode onset time, the effect of multiple ablations, the effect of the long monitoring periods quantified by the time since ICM implant, and the effect of episode duration will be used to correct the ME model.
The circadian variations of AFR were modelled by considering the time of the day of the onset of each AF episode. In order to relate AFR with circadian variation which has a cyclical nature, the time onset parameter was transformed into a sinusoid where −1 represents the middle of the day (12:00) and 1 represents the middle of the night (00:00) using the following expression: where Time 24h represents the 24 h based AF onset and Time the modified AF onset used in the analysis. Furthermore, the patients of this cohort have undergone one or two ablation procedures in their lifetimes. The 'nAblations' parameter represents the number of ablations undergone by the patient at the time of the episode onset. This parameter can either be 0 for patients with no previous ablations in episodes occurring before their first ablation, 1 for patients after their first ablation or before their second ablation, or 2 for patients with a previously failed ablation and after their second ablation. This parameter will correct for the effect on AFR of episodes occurring pre-and post-ablation as well as for patients with previously failed ablation procedures with compromised atria.
In addition, the patients of this cohort have been followed up for long monitoring periods, so the feature 'DaysSinceImplant' represents the period between the episode onset date and the date the patients were implanted with the ICM. This parameter will correct for the effects of the long monitoring periods on AFR; the change in DaysSinceImplant provides a common timescale for the episodes included in the analysis.
Lastly, in addition to the onset time and date, the Reveal LINQ also stores the total duration of the episode in minutes, i.e. 'Duration'.
The complete ME' model assumes the fixed and random effects of changes in the RR series, and corrects for circadian variation, number of ablations and episode duration: where b ¢ , 4 b ¢ , 5 b ¢ 6 and b ¢ 7 represent the fixed effects of the circadian variations, number of ablations, period between implant date and onset and episode duration respectively, and ¢ b , to the fixed and patient-specific random effects introduced in equations (18) and (19) but for the ¢ ME model.

Model evaluation
To evaluate the model's performance with unseen data, the episodes are partitioned into training and test sets. As the ME and ME' models are patient-specific, the training set was comprised of the earliest 80% of the episodes from each patient while the test set was comprised of the latest 20% of the episodes.
For comparison the FE model, that is patient-independent, was also evaluated using a training set comprised of the episodes from 80% of the patients and a test set comprised of the episodes from the remaining 20% patients.
In both cases, the accuracy was defined as the proportion of episodes where the estimated AFR deviated less than a certain absolute value from the true AFR.

Statistical analysis
Continuous data are presented as mean ± standard deviation if the null hypothesis H 0 of the Kolmogorov-Smirnov test (H : 0 data is normally distributed) was not rejected. Otherwise, continuous data are presented as median (min-max). Categorical data are presented as absolute frequency (relative frequency in percentage). The null hypothesis was rejected when p < 0.05, then set as the level of significance. The relationship between the features was studied using the patient-average of the features and Pearson's linear pair-wise correlation between them was evaluated. This average was defined for each patient as the mean feature extracted from their episodes for those normally distributed features (such as AFR, MeanRR and SampEn) and the median feature extracted from their episodes for those not normally distributed features (such as Time and Duration). The normality of the feature was determined by analysing the complete set of episodes. In addition, to study the multicollinearity among features, Belsley's collinearity test was performed on the complete set of episodes. To check for collinearity, the number of near dependencies was first determined as the number of dependencies with conditional index ( ) h p for p dependencies over the threshold h = 30. * The features will be considered to have collinearity if their variance-decomposition proportion within the near dependencies exceeds the threshold p = 0.5 * as suggested in Belsley (1991). The statistical analysis was performed using Matlab R2022a (The Mathworks Inc., Natick, Massachusetts).

Results
During the monitoring period, the ICMs stored 3739 episodes out of which 2908 (77%) episodes contained at least one 5 s segment of sufficient signal quality for estimation of AFR. Only one segment of sufficient signal quality was required as AFR was assumed to be stable during the first 2 min of the AF episode. The stability analysis showed a maximum relative error of 7.2 (4.4)% for the 24 episodes where more than 80% of the episode had acceptable levels of SQI. Figure 2 shows the evolution of the relative error (%) for varying percentages of signal.
Out of the 2908 included in the analysis, 1796 (62%) occurred before ablation, 657 (23%) occurred after ablation and 455 (15%) occurred during the 3 months blanking period and were excluded from the analysis.
Each patient had a median (min-max) of 20.5 (2-114) episodes with 16 (0-77) episodes pre-ablation, and for the patients with AF recurrence, 10 (1-61) episodes post-ablation. Figure 3(A) shows the histogram of the AFR extracted from the episodes with at least one acceptable 5 s segment in the whole population and (B) the distribution of AFR for each patient.
The full dataset of AF episodes has a median AFR (min-max) of ( ) -5.3 4.0 9.7 Hz. A high proportion (85%) of the AFR extracted from the episodes considered in the analysis were comprised between 4 and 6 Hz. A considerable variation in AFR between episodes is observed for some patients, whereas AFR remains similar between episodes for others.
The scatter plots of the patient-average AFR, modified time of AF onset, duration of AF episode (in log scale), and Mean RR, RMSSD and SampEn are shown in figure 4.
The patient-average AFR shows a mild significant linear correlation with the patient-average SampEn (R = 0.27, p < 0.05). Regarding the three HR characteristic features, the patient-average mean RR interval had a strong correlation with patient-average RMSSD (R = 0.677, p < 0.001) and a mild correlation with patientaverage SampEn (R = 0.312, p < 0.05). The patient-average time also shows a mild correlation with the partientaverage mean RR interval (R = 0.203, p < 0.05). In addition, the frequency distribution histograms of the patient-averages of the features studied are also aligned diagonally on the subpanels of figure 4. This frequency distribution histograms show that the average AF episode duration for each patient is predominantly short episodes with 54 patients (55%) having an average episode duration shorter than 20 min. Overall, the duration of the episodes had a median of 14 min with the minimum duration being the detection threshold, i.e. 2 min, and the maximum duration spanning 62.7 d.     The results of the different models are summarized in table 2, and the fitted values of AFR were plotted against the true values and illustrated in figure 5.
The FE model, which only included the fixed effect of changes in RR series had the lowest complexity but had also the largest AIC value. When also considering the random effects in the ME model, the values of AIC and the deviance decreased. Lastly, when accounting for the circadian variation, the duration of the episode, the number of ablations before the episode and correcting for the monitoring period as in the ¢ ME model, the deviance and AIC value is minimum. This suggests that the ¢ ME model outweighs the other models also when accounting for the increased model complexity. The coefficients of determination (R 2 ) also corroborate that the ¢ ME model has the better fit as its coefficient is the highest ( = R 0.56 2 ). The patient-wise root mean square error (RMSE) was computed for the ME' model and a median (minmax) RMSE of 0.35 (0.04-1.19) was found. This shows that for a subset of patients, the model is capable of estimating the AFR. In contrast, the other models had a patient-wise RMSE of 0.55 (0.24-1.74) in the case of the fixed model (FE) and 0.43 (0.02-1.39) when introducing the random effects (ME), which indicate a worse fit to the data. In addition, the Bland-Altman test (Altman and Bland 1983) is illustrated in figure 6, from which the limits of agreement can be seen to lie between 0.82 and −0.82 Hz (mean difference of 0 Hz and standard deviation of 0.42 Hz).
The estimation and the results for the fixed effect coefficients for the ¢ ME model are summarized in table 3. The parameter associated to the intercept (b ¢ = < 4.979, p 0.001 0 ) is comparable to the average AFR, see figures 3 and 4.
The irregularity in RR intervals has a significant effect on the AFR (b¢ 3 = 0.105, p < 0.05) with higher AFR for higher irregularity. The number of prior ablations also has a significant effect on the AFR (b ¢ = < p 0.168, 0.05 5 ) with higher AFR after multiple ablations, and so does the episode duration (b ¢ =´-1.182 10 , 7 5 p < 0.05) with AFR being higher for longer episodes. The models performance in unseen data, evaluated as described in section 2.4, is illustrated in figure 7. In (A), the training set consists of the first 80% of episodes for each patient and the test set consist of the remaining 20% of the episodes. From the available 2884 episodes, 2346 (81.3%) were included in the training set and the remaining 538 (18.7%) were included in the test set. The accuracy of the models was then computed for varying absolute differences between the estimated AFR and the true AFR. As before, the results show that ME models fit the training set better than the FE model. However, both models overfit the data and the accuracies for the test set are lower with fewer than 50% of the estimated episodes within 0.3 Hz of the true AFR. In (B), the patients were divided into training and test sets. As FE models do not rely on patient-specific estimates, the episodes of these patients were stratified into one of the two sets. From the available 98 patients, 78 (79.5%) were selected for the training set including 2349 episodes and the remaining 20 (20.5%) were selected for the test set including 547 episodes. In this case the accuracy of the training set and the test set were similar, however, the accuracy of both is relatively low.

Discussion
The rapidly increasing use of continuous monitoring devices for patients diagnosed with AF (Lee and Mittal 2018) offers the chance of a more detailed characterization and a better understanding of the patients' condition in order to select the most appropriate therapy but above all, the right timing in order to avoid disease deterioration. Continuous assessment of AFR extracted from ECG strips is a non-trivial matter. However, the estimation of AFR from RR series data would enable wearable based assessment of AFR, e.g. wristband PPG. The aim of the study was to model variations in AFR based on changes of RR series characteristics (mean RR interval, RMSSD and sample entropy). After a first attempt using FE modelling approach (FE model), the results were compared to more complex ME modelling approaches: considering both population and specific effects of RR series (ME model) and allowing correction for the effect of episode duration, previous ablations, and possible circadian variations ( ¢ ME model). In addition, ME modelling would account for the heterogeneity within AF patients. In order to apply the different models, the AFR stability within 2 min of AF episodes was evaluated. In the 24 episodes with enough data to run the analysis we found that for different percentages of signal, the error between the mean AFR of the signal and the mean AFR on the remaining segments was 7.2 (4.4)%. With this result in mind, AFR was considered stable and the mean AFR calculated on a single segment (5 s) was considered representative of the whole signal.
We conclude from this initial analysis that there are no strong linear correlations between AFR and RR series characteristics in the present dataset and that the relationship between RR characteristics and AFR is more complex. It should be noted that the aim of the present study is to investigate the feasibility of estimating AFR from RR series data and not to develop an AFR estimator. The small size of the present dataset prevents us from using more complex modelling approaches due to the risk of overfitting. In the study of the multicollinearity among the variables, the dataset exhibited no multicollinearity in the Belsley's collinearity test. In addition, when evaluating the model performance with unseen data, the results show that while the ME models are able to fit the data better than the FE model, they overfit the data and produce low accuracies in the test set. This points out the difficulties and unreliability of estimating AFR using such simple regression models such as ME or FE and calls for further investigation with larger datasets and more complex models to be able to reach clinically acceptable accuracies.
The correlation between AFR and HRV features describing the variability and irregularity of RR intervals has been explored before in a cohort of patients with underlying congestive heart failure (Corino et al 2013). However, to the best of our knowledge, this is the first study to assess the variation of AFR in AF patients without congestive heart failure which were continuously monitored over several months. It is also the first study to model AFR variations using ME models which account for the heterogeneity of the patient population and confounding factors.
The model's parameters showed a significant effect of RR irregularity quantified by SampEn (b ¢ = 3 0.105, p < 0.05), number of ablations (b ¢ = < p 0.168, 0.05 5 ), and episode duration (b ¢ =´-1.182 10 , 7 5 p < 0.05) on AFR. Due to the heterogeneity between patients, the ME model developed with correction for the confounding factors ( ¢ ME model) was able to better fit the data compared to the FE model, the FE modelling approach, (R 2 = 0.56 versus R 2 = 0.04).
The fact that average SampEn is correlated with average AFR ( = < R p 0.27, 0.05) and that it has a significant effect on AFR when evaluating the ¢ ME model, is in line with a previous study (Corino et al 2013) which evaluated the relationship between AFR and HRV derived features. In the aforementioned study, the regularity statistic quantifying the unpredictability of fluctuations in a time series used was the approximate entropy (ApEn) and was shown to have a Pearson's correlation of = < R p 0.26, 0.05, indicating that the higher the AFR, the less organized the RR series. In our study, the effect the regularity of the RR series has in modelling AFR was further evaluated in the ¢ ME model showing a significant effect (b ¢ 3 = 0.105, p < 0.05). Previous studies have shown that the RR irregularity during AF change in response to changes in autonomic tone induced by drugs , Cygankiewicz et al 2015 and tilt-test (Patel et al 2018). However, it should be noted that increased RR irregularity may also be a direct effect of changes in atrial electrical activity as quantified by an increased AFR or variations in the atrioventricular node conduction (Plappert et al 2021). Out of the three HRV parameters considered, only SampEn was correlated with average AFR and had significant effect on AFR when considering the ME model. However, SampEn alone does not convey as much information on the RR characteristics as when combined with the mean RR and RMSSD. A comparison between the ¢ ME model and a model excluding the mean RR and RMSSD showed that the full model had a better fit to the data (R 2 = 0.56, AIC = 4298 versus R 2 = 0.51, AIC = 4408).
Out of the 99 patients included, 19 (19%) had a previous ablation. Evidence of negative correlation between the percentage of fibrotic tissue in the left atria with the fibrillatory frequency have been reported (Swartz et al 2009). However, our study suggests that the number of ablations has a positive significant effect on AFR (b ¢ = < 0.168, p 0.05 5 ) with patients having gone through a higher number of ablations, having higher AFR.
Though further investigation with mapping and electrogram recordings is needed, this result suggests that creating lesions in the atria reduces macro-reentrant pathways but at the same time, promotes micro-reentrant circuits and/or faster, smaller rotors which would increase the fibrillatory rate of the atria.
Trying to understand AF behaviour, the link between AF episode duration and AFR had been previously studied showing a positive correlation between them ( = < R p 0.53, 0.05) (Bollmann et al 1999). The study showed that having a higher AFR at the start of the AF episode could predict longer episodes. However, this study was conducted in a small dataset with only 31 episodes from 11 paroxysmal AF patients and the correlation was describing the whole database without explicitly considering intra-and inter-patient effects. The results of the present study confirm Bollmann's results as the episode duration had a significant effect on the AFR in the ¢ ME model (b ¢ =´-1.182 10 , 7 5 < p 0.05). Finding of longer episodes of AF having higher AFR could be explained by an increase of sympathetic drive with longer duration of arrhythmia (Carnagarin et al 2019). The Reveal LINQ detects AF episodes of 2 min duration and longer, and hence AF episodes of 30 s to 2 min remain undetected. However, 'Duration' is still considered a necessary variable in the model due to the large spread of episode durations on the data: with the minimum duration being the detection threshold, i.e. 2 min, and the maximum duration spanning 62.7 d. With these large duration ranges, even if the episodes are evaluated in bins of 2 min, the effect of long versus short durations on AFR can be evaluated accurately.
In several studies, circadian dynamics of AFR were observed in patients using Holter monitoring (Bollmann et al 2000, Cosson et al 2000, Meurling et al 2001. In these studies, AFR showed a decrease at night and an increase during the morning hours, with a peak during the afternoon. In one of those studies, two different sets of patients were identified: one which showed an increase (minority) while the other showed a decrease (majority) in nocturnal AFR (Bollmann et al 2000). However, the insight in circadian behaviour of atrial electrophysiological characteristics during AF was constrained to 24 h long Holter registrations and limited datasets of up to 30 patients. The present study is based on a larger study population (99 patients) and longer monitoring periods (0.2−24.3 months), and no significant fixed effect of the onset time was observed when modelling the AFR. Circadian variation in the AFR is caused by autonomic modulation in the atrial electrical activity and could potentially be used to guide clinical strategies such as time of the day medication should be administered, or which patient would benefit the most from a catheter ablation procedure.
Currently, there is still uncertainty about the optimal selection for catheter ablation and overall, the challenges to define clear patients' phenotypes for appropriate AF management are still prominent (Karamichalakis et al 2015, Heijman et al 2018. Our study suggests that continuous assessment of AFR has the potential to estimate the impact of therapies and therefore, to help the stratification of patients towards AF ablation especially when it can help the physician decide whether or not to recommend ablation therapy and persist in drug treatment. To assess the importance of including patient specific dependencies on RR series characteristics in the model, the results were compared to an ME model that only considered b . 0 This model had a better performance compared to the FE model ( = R 0.37; 2 = AIC 4727 versus = R 0.05; 2 = AIC 5699) but had a worse fit to the data than the ME model ( = R 0.49; 2 = AIC 4424). This shows that the patients still have a high heterogeneity not being addressed by the model but some of it can be found in their RR characteristics.
Detailed results for the estimated random effects of the ME' model can be found in the supplementary material. For most patients, the random effects coefficients were non-significant. However, for a subset of patients, the HR variables had negative significant effect on the AFR (mean RR: 6, RMSSD: 5, and SampEn: 3 patients) while for another subset, they had a positive significant effect (mean RR: 8, RMSSD: 6, and SampEn: 8 patients). The rest of the random effect coefficients had a negative significant effect on the AFR of a subset of patients (Time: 5, number of Ablations: 2, Days from Implant: 5, and Duration: 1 patient) while for another subset of patients, they had a positive significant effect (Time: 3, number of Ablations: 2, Days from Implant: 4, and Duration: 4 patients). For the rest of the patients, there were no significant random effect coefficients. This shows the heterogeneity in the dataset and proves the rationale of needing ME models to analyze such cohort. This retrospective study combining 2 different cohorts with limited clinical baseline data, is selected from a patient population implanted with the Reveal LINQ based on clinical indications including suspected AF, AF ablation monitoring or AF management. Although this unique database offers the advantage of long monitoring periods and a complete monitorization of patients suffering AF, some limitations should be noted. The patient population included in the study is relatively young (57 ± 12 years) compared to the general AF population and with a low degree of cardiovascular risk. Due to the retrospective nature of the study, the medication administered to each patient during the monitoring period and the ablation technique in those patients with a previous failed ablation were not available. In a similar manner, the gender of the patient was not available and its influence on AFR could not be assessed. In any case, possible influence of medication, scar tissue from previous ablations and gender on AFR is modelled as a patient specific random effect. Furthermore, due to the limited memory of the device, ECG data was only stored for a subset of the detected episodes; 57 840 episodes were detected by the device out of which ECG data from 3739 (6.5%) were stored. The number of episodes detected but not stored is not only linked to the number of episodes suffered by the patient but also to the frequency of visits to the hospital as the data was downloaded and saved each time the patient had a check-up. Overall, the patients had a median (min-max) of 130 (1-4564) episodes detected pre-ablation out of which 43.4 (1.9-100)% were stored. Post-ablation, the device detected 1 (0-4245) episodes out of which 60.2 (2.9-100)% were stored. The aim of this study involves analyzing the relationship between AFR and HRV features during the initial 2 min of the episodes while also accounting for various confounding factors. However, these first 2 min may not be representative of the whole episode; results from a previous study suggest that AFR may accelerate during the first 3-4 h . Therefore, the relationship between AFR and the HRV features might be different in other instances of the episode.
The extraction of atrial activity has been a difficult problem since the atrial and ventricular activities overlap spectrally and there is no uniform standard. There are two main approaches that exploit the property that atrial and ventricular activities arise from different bioelectrical sources: principal component analysis (PCA) (Raine et al 2004) and independent component analysis (ICA) (Rieta et al 2004). However, these approaches derive a global atrial signal with contributions from all the leads while the QRST cancellation technique used in this study (Stridh and Sörnmo 2001), extracts the atrial signal in a specific lead and strives at not modifying the original f-wave morphology. Despite their differences, in a comparative study between these three algorithms applied to the same surface ECG recordings, the dominant rate of the atrial signals extracted by each of the algorithms was within 1Hz (Langley et al 2002). The performance of the QRST cancellation technique was explored in both AF simulated signals and ECG data with AF and it was showed to be well-suited for atrial activation extraction (Stridh and Sörnmo 2001). In addition, there have been direct comparisons between endocardial signals and surface ECG signals subjected to the QRST cancellation technique used, summarized in (Bollmann et al 2006). The feasibility of extracting atrial activation from ICM data has been previously explored  where patients with paroxysmal AF received an implantable loop recorder (same one used in this study) and the AFR was extracted from episodes longer than 1 h also using the AFR tracker software provided by CardioLund. Conversely, in this study AFR was extracted from any episode longer than 2 min detected by the ICM and from both paroxysmal and persistent patients. The harmonic f-wave model used for AFR estimation was originally developed for surface ECG, and the characteristics of the f-waves in the ECG recorded by the ICM may differ from these since the electrodes are placed next to the apex of the heart (with V2-V3 electrode orientation). However, our results indicate that the model fit was sufficient in most cases; 2908 of the 3739 detected AF episodes in our 99-patient cohort had sufficient signal quality to be analysed.

Conclusion
Fixed and mixed effects modelling approaches were used to investigate the effect of changes in RR series characteristics corrected for episode onset and duration, previous ablations, and onset date, on variations in AFR in a study population of 99 patients monitored for 9.2 (0.2-24.3) months as median (min-max). The ME modelling approach was shown to be superior to the FE modelling approach due to the heterogeneity of the patient population and the presence of confounding factors. The fixed effects extracted from the ¢ ME model showed that AFR is slightly higher in episodes with less organized RR series and of longer duration and is affected by catheter ablations. The use of ME models combined with long term monitoring of patients offers the chance of continuously estimating the AFR from RR series and episode-based characteristic and will lead to a more detailed characterization and a better understanding of the patients' condition which could potentially aid the clinicians in their decision-making process. However, the results obtained point out the difficulties and unreliability of estimating AFR using such simple regression models. Further research with larger datasets that would allow for more complex models is warranted.