The Dissimilar Impact in Atrial Substrate Modificationof Left and Right Pulmonary Veins Isolation after Catheter Ablation of Paroxysmal Atrial Fibrillation

Since the discovery of pulmonary veins (PVs) as foci of atrial fibrillation (AF), the commonest cardiac arrhythmia, investigation revolves around PVs catheter ablation (CA) results. Notwithstanding, CA process itself is rather neglected. We aim to decompose crucial CA steps: coronary sinus (CS) catheterization and the impact of left and right PVs isolation (LPVI, RPVI), separately. We recruited 40 paroxysmal AF patients undergoing first-time CA and obtained five-minute lead II and bipolar CS recordings during sinus rhythm (SR) before CA (B), after LPVI (L) and after RPVI (R). Among others, duration, amplitude and atrial-rate variability (ARV) were calculated for P-waves and CS local activation waves (LAWs). LAWs features were compared among CS channels for reliability analysis. P-waves and LAWs features were compared after each ablation step (B, L, R). CS channels: amplitude and area were different between distal/medial (p≤0.0014) and distal/mid-proximal channels (p≤0.0025). Medial and distal showed the most and least coherent values, respectively. Correlation was higher in proximal (≥93%) than distal (≤91%) areas. P-waves: duration was significantly shortened after LPVI (after L: p=0.0012, −13.30%). LAWs: insignificant variations. ARV modification was more prominent in LAWs (L: >+73.12%, p≤0.0480, R: <−33.94%, p≤0.0642). Medial/mid-proximal channels are recommended during SR. CS LAWs are not significantly affected by CA but they describe more precisely CA-induced ARV modifications. LPVI provokes the highest impact in paroxysmal AF CA, significantly modifying P-wave duration.


Introduction
Atrial fibrillation (AF) is the prevailing cardiac arrhythmia in the western world. Prolonged lifespan and the connection with a plenty of other comorbidities contribute to the ever-growing AF incidence. Health and economic burden caused by AF alert the need for thorough investigation on its pathophysiology [1]. AF springs principally from pulmonary veins (PVs) [2] and propagates through cardiac structures [3]. The main mechanism assisting the AF propagation is structural remodeling and fibrosis is especially contributing to the alteration of the cardiac anatomy, causing conduction heterogeneity, hence favoring the AF perpetuation [1,4]. Although conduction heterogeneity is more prominent during AF, the anatomical substrate can still be present for both atria even substrate modification or which channels of CS catheter are the most appropriate for the analysis are two vital issues that remain unexplored. During CS cannulation, the most proximal pair of electrodes (9)(10) is placed close to RA and the most distal pair (1)(2) close to LA [44,47]. Notwithstanding, CS catheterization may be rather challenging due to variable CS anatomy and shape, aggrravated by myocardial contraction or the existence of CS dilation, factors that can lead to unstable recordings, especially from the distal tip of the catheter [50][51][52][53]. Additionally, anatomical alterations of the, adjacent to CS extremes, mitral annulus across the cardiac cycle in SR may affect furthermore the stability of recordings acquired from distal and proximal electrodes of the CS catheter [54,55]. Considering the aforementioned factors, information recorded across the CS catheter electrodes could vary significantly and the choice of the appropriate channel recruited for the analysis should be made with extreme caution.
The present work aims to elucidate the aforeposed issues regarding the CA procedure, in order to arise the understanding on the mechanisms of important CA steps and their interaction with the CA result. In the first place, the ability of CS channels to describe with the highest precision possible the AF dynamics during SR is assessed and the most and least recommended CS channels are defined. Afterwards, the relevance of CS in substrate modification evaluation due to CA is investigated via analysis of features traditionally applied to ECG recordings and cross-referenced by P-waves and HRV analysis of the ECGs. Finally, the evolution of P-waves and CS LAWs after isolation of either sides of PVs is tracked in order to define the PV side that has the highest impact in atrial substrate modification due to CA.
The manuscript is organized as follows. Section 2 briefly describes the database recruited for the analysis, as well as the preprocessing and analysis steps. Section 3 presents the results, which are further interpreted in Section 4. Main findings are stated in Section 5.

Database
Initial database consisted of 61 paroxysmal AF patients without any previous CA sessions. Twenty-one patients were discarded due to extremely low amplitude or presence of noise and artifacts in the extreme channels of the CS recordings, probably due to dynamical change of mitral annulus anatomy during SR and vigorous myocardial contraction causing the movement of the CS cathether. The final database consisted of the remaining 40 patients.
Recordings from a standard 12-lead electrocardiogram (ECG) and a decapolar CS catheter with sampling frequency at 1 kHz were acquired by a Labsystem™ PRO EP recording system (Boston Scientific, Marlborough, MA, USA). Five-minute continuous segments before RFCA initiation (step B), after LPVI (step L) and after RPVI, which coincides with the end of the RFCA procedure ( step R), were chosen. The evolution of each ablation step is illustrated in Figure 1a. Step B corresponds to 0% of the procedure, step L corresponds to 100% of LPVI, 0% of RPVI and 50% of the total ablation procedure, while step R corresponds to 100% of the overall procedure. The effect of ablation of each PV side is shown in Figure 1b. A statistically significant difference in features between steps B and L indicates that LPVI is critical in atrial substrate modification, while a difference between steps L and R or between B and R but without the difference between steps B and L being statistically significant, proves a significant effect of RPVI. Surface analysis was limited to lead II, as P-waves are more prominent in this lead [56].
For the study of the reliability of CS channels in preserving the AF dynamics, final database consisted of a total of 58 step B or R recordings from the 40 patients of the database. CS catheter consisted of the following channels of bipolar signals: distal (D), mid-distal (MD), medial (M), mid-proximal (MP) and proximal (P), with D channel being the closest to LA and P channel to the RA. Firstly, a multichannel comparison allowed us to define the channels that can record the AF dynamics to the most reliable degree. Selection of the analysis channel for the CS study was then performed among the most robust channels, with unique favorable criteria the high signal amplitude and low baseline fluctuation. Specific attention was paid so that the same channel would be employed for all CA steps for the same patient. Steps of CA procedure for which recordings were extracted and analyzed. In step B, no ablation has been performed yet (0%). In step L, LPVI has been completed (100%) and hence, we are in the middle of the procedure (50%).

RPVI
Step R corresponds to RPVI and to the end of the procedure. Therefore, each step is completed (100% Regarding the CA procedure, all patients underwent circumferential RFCA of PVs, guided by 3-D electroanatomical mapping during SR. RFCA was initiated by performing a crown surrounding left PVs (step L), followed by a crown surround right PVs (step R). Non-inducibility of AF was confirmed by pacing in all patients and was the endpoint of the procedure.

Preprocessing
For ECG recordings, powerline interference and high frequency muscle noise were removed by a wavelet-based denoising method [57] followed by a bidirectional low-pass filter with cut-off requency at 70 Hz [58]. Baseline wander was also removed [58].
EGM denoising and mean removal were the first preprocessing steps for CS recordings analysis [59]. Although presence of ventricular activity is not dominant in atrial bipolar signals, far-field activity in line with the R-peak of the ECG recordings has been observed in some cases. Removal of ventricular activity was performed with an adaptive cancellation method [60].
Next, ectopic beats correction was performed. Ectopic beats are premature atrial or ventricular contractions that affect the HRV. In our analysis, ectopic beats in ECGs, if present, did not exceed 4% of total beats. Their correction included the detection and cancellation of the premature complexes and their replacement by a new beat via linear interpolation [61]. Among various ectopic replacement methods, linear interpolation was chosen due to its better performance for time-domain HRV features [62].
Finally, detection and delineation of atrial activations was carried out. P-waves were firstly detected by an adaptive search window prior to the R-peak [63] and then delineated [64]. Local activation waves (LAWs) of CS were detected with an algorithm based on an alternative Botteron's technique [65]. Delineation was performed by firstly smoothing the LAW with a five-point moving average filter [66]. Delineation of both ECG and intracardiac recordings was visually inspected and corrected, if needed, by an expert.

Main Analysis
Once preprocessed, the duration, amplitude, root mean square (RMS) value, area, number of deflections and inflections (NODI) and slope rate were calculated for P-waves and LAWs, as shown in Figure 2. Final values of these features were calculated by signalaveraging. A brief description of these characteristics is provided as follows. Further details are described elsewhere [66].

1.
Duration: Distance between the onset and offset of each activation.

2.
Amplitude: Amplitude of positive and negative maximum of each activation were considered as positive (PosAmp) and negative amplitude (NegAmp), respectively. Peakto-peak amplitude (PPAmp) was the distance between positive and negative maximum points. As P-waves are positive in lead II, only maximum amplitude was calculated for ECG analysis.

3.
RMS: Let X n be a time-series, so that X n = {x 1 , x 2 , . . . , x n }. RMS value is the quadratic mean of the function that defines the time-series. In our case, this function is defined by either the P-wave or LAW waveform.

4.
Area: Area is calculated as the integration of the signal over the time interval. Trapezoidal method allows this integration, by splitting each signal into smaller and easier to calculate trapezoids. Final area is defined by the cumulative sum of these trapezoids. As LAWs contain both positive (PosAr) and negative (NegAr) parts, this method was separately applied to each one of them.

5.
NODI: Deflections and inflections were calculated from the points that cross two auxiliary baselines, at ±25% of the signal amplitude. This metric was only calculated for LAWs, as P-waves do not show multiple major deflections and inflections. 6.
Slope rate: The rhythm of increasing or decreasing slope was calculated at sample points equal to i% of the activation duration, with i = 5, 10, 20. Slope rate at the maximum point was also computed. The equation calculating these slope rates was the following: where Amp(t i ) is the amplitude at the i% of the activation duration, Amp(t onset ) is the amplitude at the onset of the activation, t i is the sample point at the i% of the activation duration and t onset is the sample point corresponding to the onset of the activation.
Afterwards, features calculated across each recording were analyzed: morphology variability (MV), dispersion and time-domain HRV features.

1.
MV: A reference signal was firstly created by the 20 most similar activations of the channel under analysis and then correlated with each and every activation, using an adaptive signed correlation index (ASCI) with 12% tolerance [67]. MV was then defined as the percentage of signals that correlated <95% with the reference signal.

2.
Dispersion: Traditionally, for the calculation of dispersion, more than one ECG lead is employed and the difference between maximum and minimum activation duration across channels is computed. Alternatively, in our case, lead II was just extracted and P-wave dispersion analysis was performed in this channel because atrial activity presents the highest amplitude. Dispersion was then defined as the difference between the 25th and 75th percentiles of the atrial activations duration of each recording. This way, the effect of signal delineation accuracy is minimized and an extremely long or short activation caused by various factors will not affect significantly the results [68]. Dispersion of EGM recordings was calculated in the same way.

3.
Time-domain HRV features: HRV analysis is normally performed on R-R intervals, thus describing ventricular response. As in the present work we are focused on atrial analysis, we modified the techniques by substituting R-R peaks by P-wave to P-wave for ECG and LAW to LAW for EGM recordings. As these features describe the atrial response, thus neglecting the effect of the atrioventricular node and other cardiac structures, they will be referred in the remainder of this document as atrial rate variability (ARV) features. Standard deviation of normal-to-normal beat interval (SDNN), variance of normal to normal beat interval (VARNN) and RMS of successive interbeat differences (RMSSD) were calculated for each recording.

Heart Rate Adjustment
Time-domain features of P-wave analysis are affected by variable HR [69]. More specifically, as HR increases, intervals between fiducial points of ECGs shorten. Therefore, HR adjustment is proposed in order to moderate this effect. For this purpose, additionally to the original analysis, a simple HR-adjustment factor is performed. As sampling frequency is 1 kHz, a 60 beat-per-minute recording would show one activation every 1000 sample points. However, as HR is diverse and often deviant from these values, the adjustment factor for the i th activation was set as where IBI i is the interbeat interval between the i th and the (i − 1) th activations. Duration and area were normalized by this factor, while slope rates were inversely scaled by it. HR-adjusted values will be shown as HRA(y), where y = duration, area or S i .

Statistical Analysis
Normality and homoscedasticity were tested with Saphiro-Wilk and Levene tests, respectively [70,71]. According to the results, non-parametric tests were employed for the comparison between populations.
Reliability analysis of the CS channels with respect to AF dynamics was performed on Duration, Amplitude, RMS, Area and NODI. In the first place, a multichannel comparison via a Kruskal-Wallis (KW) test was employed [72]. Comparison in pairs of two channels was performed with a Mann-Whitney U-test (MWU) [73] with Bonferroni correction. Median values were also calculated at each channel and any significant differences between each one and the remaining channels was explored by as well using a MWU with Bonferroni correction. Afterwards, a reference signal for each channel was calculated as described from MV analysis in Section 2.3. Then, the correlation between the morphology of each channel's reference signal in pairs of two was calculated for each recording, using an ASCI with 12% tolerance.
For P-waves, LAWs and ARV analysis, comparison between ablation steps is performed with KW and post-hoc tests to define which step is crucial are performed with MWU with Bonferroni correction and median and interquartile range calculations. Analysis is performed for P-waves and CS LAWs separately. Percentage of variation (POV) of features was specified for each recording and between two ablation steps as where r i is the recording of the ith patient, V 2 is the value of each feature at the posterior step and V 1 is the value of the same feature at the prior step. How POV was modified between ablation steps B-L and L-R for both P-waves and LAWs was tested with MWU. Finally, HR was measured at each ablation step and compared among all steps with KW. Table 1 shows the results of the multichannel comparison as well as the comparison in pairs of channels, for the selected features. Amplitude and Area show different values among the CS catheter channels. Paired analysis showed that differences are mostly located between D and M or D and MP channels. A trend has also been observed between values of MD and M channels. Due to Bonferroni correction, α is 0.005, as 10 paired comparisons have been performed. Comparison between D-MD, MD-P and M-MP channels are not illustrated, as they did not show any statistically significant differences (p > 0.03, 0.38 and 0.35, respectively).   Figure 3 shows the median values obtained for each analyzed feature from the CS. As can be seen, the distal channel showed the longest LAWs Duration and lowest LAWs Amplitude and Area values for most of the features. Medial channel contained the shortest LAWs Duration and highest LAWs Amplitude values. Area was higher in mid-proximal channel, followed by medial channel. Mid-proximal channel showed overall similar behavior as medial channel, with rather high Amplitude and short LAWs Duration values.  The analysis of channels that varied significantly with respect to the others at each feature is shown in Table 2. Statistical comparison allows the detection of the channels that show the most and least deviant values and can corroborate the results shown in Figure 3. None of the channels showed a statistically different LAWs Duration comparing to the remaining ones. Nevertheless, a trend was observed for distal and medial channels. Distal and medial channels showed additionally statistically different values with respect to the remaining channels regarding LAWs Amplitude and Area features. Mid-proximal channel also showed statistically significant difference than the other channels for positive Amplitude and a trend for the remaining Amplitude and Area features. Combining the aforementioned observations with the median values of the LAWs features for each CS channel presented in Figure 3, we can conclude that Ampitude and Area are statistically smaller in distal channel and higher in medial channel of the catheter, while duration tends to be longer in distal and shorter in medial channels. Mid-proximal channel shows a trend for high Amplitude and Area values as well.

Analysis of CS Features between Channels
Finally, how LAWs morphology of each channel correlated with the morphology of LAWs at each of the remaining channels can be appreciated from Figure 4. Channels of proximal area (MP, P) show higher correlations between their LAWs morphology compared to correlations from distal area (D, MD). Additionally, medial channel showed stronger correlation with proximal (93-95%) than distal area (86-91%). Although all adjacent channels showed relatively strong correlations, the highest values were observed in proximal area.

Analysis of Features from P-waves and LAWs
The HR measurements did not reveal any statistical difference between HR at different ablation steps. However, a decrease in HR was observed in step L, as shown in Table 3. Table 3. Heart-rate at each ablation step and comparison between three steps (KW). As result indicated a non-significant comparison, no MWU has been performed. KW: Kruskal-Wallis; MWU: Mann-Whitney U-test; B: before CA; L: after LPVI; R: after RPVI.  Table 4 shows the median and interquartile range for the calculated features in Pwaves at each ablation step. Multiple comparison among ablation steps is then shown in the fifth column (KW) and comparison of each feature between ablation steps in pairs of two using Bonferroni correction can be observed in the last three columns (MWU). Duration varied statistically among channels. When analysis in pairs was conducted, this variation was detected between steps B-L and B-R. As L-R comparison did not show any significant variation, the B-R significance is probably due to the B-L significance. Hence, step L is considered the critical step for the reduction in Duration, which falls from 120 ms in the beginning of the procedure to 104 ms, almost the final value observed after the end of CA of PVs. In the same context, amplitude showed a downward trend after step L. Despite the key role of step L in the modification of Duration, when HR adjustment was performed, this step did not show any effect in HRA(Duration) and neither did step R. However, modification of HRA(Duration) showed a trend when values in the beginning and the end of the procedure were measured (B-R comparison), possibly due to a cumulative effect of CA in AF substrate which can be better appreciated quantitatively when isolation is totally performed. Although Area modification was not originally found to be statistically significant at each of the ablation steps, HR adjustment revealed for it a trend in B-L comparison, falling from 26.10 mV×ms to 19.40 mV × ms. Finally, it can also be observed that measurements after step L showed, in a non-significant level, lower Amplitude and Area values and higher ARV, MV values and Slope rate values than steps B and R.
Regarding LAWs, comparison of each feature among all channels revealed a significant variation of RMSSD. Apart from this observation, none of the features varied statistically at any of the ablation steps. However, some trends have been observed. The respective results are shown in Table 5. LAWs showed a trend for shortening between steps B-R, as shown by both Duration and HRA(Duration) results. MV showed a trend for amplification after step L (B-L comparison) which was almost statistically significant. All ARV features employed in this study showed an increasing trend after step L (B-L comparison) and a decreasing trend after step R (L-R comparison). Note that threshold for MWU is α = 0.0167 due to Bonferroni correction. POV of all features between each two ablation steps are illustrated in Table 6 for both P-waves and LAWs. Comparison of POV that corresponds to successive step transitions B-L and L-R can also be seen in the last two columns. In general, POV in P-waves seem to be more prominent than respective POV values in LAWs, especially in the B-L comparison. Duration in P-waves got reduced by −13.30% after step L, while in LAWs by −5.49%.
Step R did not additionally modify at average the P-waves at all (0.00% POV in L-R comparison). For LAWs, step R did reduce LAW Duration by an additional −1.93%. Comparison between B-L and L-R alterations in P-wave Duration presented a significant difference (p < 0.0001), which was not observed in LAWs, indicating a different effect of LPVI and RPVI in P-wave but not in CS LAWs Duration alteration. Table 6. POV for between every two ablation steps for P-waves and LAWs and comparison between POV of successive step transitions B-L and L-R. Statistically significant results are shown in (*). POV: percentage of variation. HR-adjustment did not have a different effect in B-L comparison of P-waves with respect to the same comparison for non-normalized Duration variation (−13.73% vs. −13.30%, respectively). However, HR-adjustment revealed a slightly higher, albeit not statistically significant, effect of step R in HRA(Duration), showing an incrementation of +1.60%. For LAWs, HRA(Duration) slightly mitigated the effect of step L (−3.89% for HR-adjustment) and potentiated the effect of step R (−4.46%). These results were not statistically significant either. Amplitude and HRA(PosAr) values showed the same kind of variations for P-waves and LAWs, with comparison between B-L and L-R transitions showing a trend for P-waves and low statistical power for LAWs.

Features P-Waves LAWs P-Waves LAWs P-Waves
For the remaining features, MV showed a rather strong magnification after step L by +144.90% in LAWs, whereas the corresponding step in P-waves showed only a +15.00% of MV magnification. After step R, MV dropped by −5.93% in LAWs and by −6.11% in Pwaves. Due to very intense variations in MV of LAWs, POV evolution between transitions from B-L and L-R varied significantly (p = 0.0176). As already observed in Tables 4 and 5, ARV features showed a notable incrementation after step L in both P-waves and LAWs. This is even more prominent in POV analysis, where ARV got increased by up to +225.90% in LAWs and up to +64.86% in P-waves (in VARNN in both cases). However, ARV after the end of the procedure (step R) got decreased by up to −55.91% in LAWs and −42.64% in P-waves (in VARNN in both cases). These very high variations in LAWs analysis led the B-L and L-R comparison to be statistically significant for all ARV features. As P-waves POV analysis also showed considerable variations but to a lesser extent, only POV in RMSSD varied statistically between B-L and L-R. Figure 5 then shows the POV box and whisker plots for the features that showed any kind of significant alteration or a trend, so that alterations in POV can be better illustrated. From y−axis can be additionally observed that POV in LAWs shows more scattered values that span along a wider range than P-waves analysis. This is noticeable at most of the boxplot pairs of Figure 5, but can be especially seen in the subplot of MV, where POV in P-waves is in the range of 0 to 1200% while in LAWs in the range of 0 to 12,000%. [%] [%] [%] [%] [%] [%] [%] [%] [%]

Discussion
This multi-approach study had three main objectives. First of all, to define the CS channels that can record with the highest precision and robustness the AF dynamics during SR. The analysis revealed the existence of variability among CS channels especially in Duration and Amplitude features. Differences were mostly found betweeen distal and medial and between distal and mid-proximal, with a trend between mid-distal and medial channels. A combined interpretation of the analysis of medians and the one-vs-all analysis indicates that distal channel showed the longest Duration, whilst the shortest Duration has been recorded by the medial channel. Regarding Amplitude and Area values, these were smaller in distal channel and larger in medial and mid-proximal channels. Proximal area showed the strongest morphological correlations between its channels. On the contrary, correlations in distal area or between distal and proximal channels were weaker. Hence, being the least susceptible channels to exogenous factors during SR, medial and mid-proximal channels are recommended while distal and mid-distal channels are not suggested.
Various studies corroborate these conclusions. CS EGM fractionation analysis was related with AF recurrence during SR in proximal and medial but not in distal channel in a recent study [47]. Fractionation of proximal CS EGMs indicated AF patients in another study employing recordings during AF [74]. Another work found that AF cycle length analysis in distal channel failed to predict AF termination after CA. In the same study, only mid-proximal channel predicted freedom from AF recurrence [75].
The second objective of the present study was to investigate if CS recordings can describe adequately the substrate modification due to CA, as observed by P-wave analysis. Parallel P-waves and LAWs analysis has been conducted for this purpose. The former represent the entire atria while the latter provide very specific yet crucial information on CS function. Variations were observed to a higher degree in most of the features in P-waves than LAWs. MV and ARV features were modified to an exceptionally higher extent in LAWs than P-waves. As CS recordings are closer to the tissue under ablation than surface ECG recordings, variability caused by RF energy deliverance may be illustrated with higher precision by LAWs [38,39]. Variation was more consistent in P-waves than LAWs, while the latter showed higher Dispersion in values across all features.
The last but not least purpose of this work was the evaluation of additional recordings acquired during CA in order to understand the role that the ablation of each PV side plays to the modification of the studied features and, as a consequence, to the atrial substrate alteration. A significant P-wave shortening was observed after CA of both PV sides, in line with a plethora of previous studies mainly attributed to fibrotic areas causing conduction delays [20][21][22][23]. Interestingly enough, this reduction was observed right after LPVI, with RPVI not showing any additional effect to this feature. Duration before CA was 120 ms, dropping down to 104 ms after LPVI and showing a minor increase after RPVI, to 106.5 ms. P-wave Amplitude also tended to show a lower value after LPVI which was slightly increased after RPVI, but remained overall smaller than the pre-ablative measurements.
HRV attenuation after CA is considered an indicator of CA success [38][39][40]. In line with previous studies, ARV in the present study showed a non-significant reduction after the end of the procedure. Nevertheless, recordings obtained after LPVI showed a trend for amplification of ARV values by up to +64.86%. Previous works studying the effect of RF energy in rabbits and students in lying position found that RF exposure can cause HRV incrementation and HR attenuation [76,77]. These findings explain the aforementioned results of the present study. Additionally, HR was indeed found to decrease after LPVI in the present work. Finally, HR-adjustment preserved the variation that was observed in Duration after LPVI, although losing statistical power. It also incremented the different effect that RPVI had on Duration, showing a slight non-statistical increase of +1.60% after RPVI. The reason for this slight incrementation is the HR acceleration after the end of CA procedure with respect to recordings after LPVI, where RF energy deliverance was still going on. Higher HR leads to generally narrower P-waves, the size of which is retrieved after HR-adjustment. An additional factor that may have a minor effect on this incrementation is a possible deviation of 1-2 ms in P-wave delineation precision, which due to its size (<+2.00% of duration values range) is considered acceptable.
Amplitude in P-waves showed a trend for reduction after LPVI. Although final Amplitude was non-statistically reduced with respect to the beginning of the procedure, as in previous studies [27], RPVI slightly but non-significantly increased Amplitude values. Contrastly, Slope rate was non-significantly increased after LPVI but decreased after RPVI in most of the studied time instances. It is highly possible that RF exposure also has had an effect on Amplitude and Slope rate features, explaining these variations.
Regarding LAWs, variations of most features show weaker statistical power and lower POV at each ablation step with respect to P-waves analysis. As mentioned afore, MV and ARV varied more prominently in LAWs than P-waves. MV showed a high incrementation in the order of +144.9% after LPVI and a decrease of −5.93% after RPVI. The reasons for these dramatic changes are not clear. Exposure to RF energy, not only affecting ARV but also MV may be an explanation. Recently a study found MV in lead V1 P-waves of paroxysmal AF undergoing CA of PVs to decrease after the procedure, as a sign of a successful ablation [78]. Apart from the fact that these results come from P-waves analysis, which in our case show a slight attenuation overall, MV in the present analysis was extracted by a template of the 20 most similar activations and not by considering all activations of one recording. ARV was also increased after LPVI in LAWs and to a higher extent than P-waves (up to +225.9%), possibly due to proximity to the tissue under ablation as explained previously.
To our knowledge, this is the first complete comparative study to perform simultaneous P-waves and LAWs analysis on recordings obtained not only before and after but also during CA of PVs. Recently, a relevant study calculated the organization of ECG recordings before, during and after CA of PVs and did not find any step to affect significantly the organization indices under calculation [79]. As many differences exist with respect to our study, a comparison would not be straightforward. In the first place, individuals studied were persistent AF patients while the present study employed exclusively paroxysmal AF patients. As persistent AF shows more complicated atrial substrate and often presents AF drivers outside of PVs, efficiency of CA of PVs is notably lower with respect to success rates in paroxysmal AF patients [19,80]. Secondly, the procedure consisted of CA of PVs, CA of CFAEs and linear CA of LA. Recordings during the procedure were the recordings after CA of PVs and before CA of CFAEs and linear CA of LA. Hence, the intermediate stage of their analyis would be the final stage of ours and no information is provided about the contribution of left or right PVs. Finally, features employed are different from the features employed in the present study.
The key aspects of CA of PVs for paroxysmal AF patients investigated in the present study improve significantly the understanding of the AF mechanisms during SR and contribute to the knowledge on how these mechanisms respond to each step of CA. Moreover, the CA procedure itself is reconsidered and the most reliable means to analyze CS EGMs are explored. Overall, a more detailed perspective of the CA procedure and the effect of RF exposure to atrial tissue is obtained.

Conclusions
LPVI is the critical part of CA of PVs for paroxysmal AF patients, altering significantly the P-wave duration. RF exposure tends to cause temporary ARV incrementation, which is reversed right after the end of the CA procedure. The effect of CA of PVs on CS is less straightforward and takes place to a lesser extent. Thus, other atrial structures may be more indicative of the ablation outcome and should be assessed as alternative references.
It should be noted, however, that ARV modifications regarding RF energy are more prominently observed in CS LAWs, possibly due to the vicinity with the tissue under RF exposure. Hence, the employment of CS recordings may be beneficial for the study of ARV alterations during and after CA of PVs.
Finally, studies interested in employing CS analysis are encouraged to extract and investigate medial or mid-proximal channels, as they were found to be the most robust, showing the highest coherence between LAWs morphologies. Distal and mid-distal channels, on the other hand, should be avoided as they were prone to variable morphology and less clear activations.

Informed Consent Statement:
Written informed consent was granted from all the subjects participating in the present research. All acquired data were anonymized before processing.

Data Availability Statement:
The data supporting reported results and presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors have no association with commercial entities that could be viewed as having an interest in the general area of the submitted manuscript. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.