Multivariate classification of Brugada syndrome patients based on autonomic response to exercise testing

Ventricular arrhythmias in Brugada syndrome (BS) typically occur at rest and especially during sleep, suggesting that changes in the autonomic modulation may play an important role in arrhythmogenesis. The autonomic response to exercise and subsequent recovery was evaluated on 105 patients diagnosed with BS (twenty-four were symptomatic), by means of a time-frequency heart rate variability (HRV) analysis, so as to propose a novel predictive model capable of distinguishing symptomatic and asymptomatic BS populations. During incremental exercise, symptomatic patients showed higher HFnu values, probably related to an increased parasympathetic modulation, with respect to asymptomatic subjects. In addition, those extracted HRV features best distinguishing between populations were selected using a two-step feature selection approach, so as to build a linear discriminant analysis (LDA) classifier. The final features subset included one third of the total amount of extracted autonomic markers, mostly acquired during incremental exercise and active recovery, thus evidencing the relevance of these test segments in BS patients classification. The derived predictive model showed an improved performance with respect to previous works in the field (AUC = 0.92 ± 0.01; Se = 0.91 ± 0.06; Sp = 0.90 ± 0.05). Therefore, based on these findings, some of the analyzed HRV markers and the proposed model could be useful for risk stratification in Brugada syndrome.


Introduction
Brugada syndrome (BS) is an inherited disease presenting a typical pattern on the electrocardiogram (ECG), characterized by a distinct ST-segment elevation in right precordial leads, associated with a high risk for unexpected sudden cardiac death (SCD), secondary to ventricular fibrillation (VF) in absence of any apparent structural cardiopathy [1,2]. Since its initial description in 1992 as a new cardiac syndrome [3], BS has raised a great interest due to its high incidence, especially in far eastern countries, and its association with sudden death in young a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 adults and, less frequently, in infants and children. It has been estimated that BS is responsible for 4-12% of the total amount of SCD and for 20% of SCD in patients with structurally normal hearts [4,5].
Although several methods have been evaluated for the prediction of VF occurrence, studies based on the largest clinical series including BS patients only proved two consistent and reliable predictors of major cardiac events: documented symptoms and spontaneous type 1 Brugada-like ECG pattern [5,6]. As a consequence, risk stratification in order to determine the best treatment approach for these patients still remains challenging, especially for asymptomatic individuals without documented VF episodes.
Ventricular arrhythmias (VA) in BS typically occur at rest and especially during sleep, suggesting that parasympathetic activity may play a relevant role in the arrhythmogenesis of the disease [7,8]. Moreover, a sympathetic autonomic dysfunction on BS patients has been reported in previous works on cardiac autonomic nervous system (ANS) analysis based on positron emission tomography [8][9][10][11]. Thus, changes in the autonomic modulation captured by heart rate variability (HRV) analysis may provide useful information for the prediction of VA in these patients. Indeed, the autonomic function has already been studied in BS, but most previously reported autonomic markers are based on long-term measurements, being timeconsuming and leading to contradictory results [12][13][14][15][16][17][18]. The evaluation of the autonomic response can be better characterized by stimulating the ANS in a controlled and repeatable fashion, by applying standardized maneuvers such as physical stress testing, pharmacological stimulations or the head-up tilt test. Physical exercise induces an increase in sympathetic activity and a parasympathetic withdrawal, resulting in higher heart rates (HR). Conversely, postexercise cardio-deceleration is mediated by a progressive rise in vagal activity [19], as well as a continued sympathetic recession [20]. Indeed, previous works have reported the potential of exercise testing to predict VA in patients suffering from BS [21][22][23][24][25]. Nevertheless, we are not aware of any study having analyzed the temporal progression of the autonomic response under conditions of exercise and recovery in this population.
In this work, a time-frequency approach was applied on 105 BS patients at different levels of risk (symptomatic and asymptomatic) so as to characterize the temporal evolution of several HRV features in response to exercise. Then, a multivariate approach based on a step-based machine learning method was implemented to identify those extracted HRV features best distinguishing between populations. Based on the hypothesis that changes in the autonomic function could improve prognosis interpretation in BS, the main objective of the study was to build a multivariate classifier capable of identifying patients at high risk.

Study population
The standard 12-lead ECG recordings from 118 consecutive patients diagnosed with Brugada syndrome who took part in a physical stress test were collected during a prospective, multicentric study conducted between 2009 and 2013 in the Cardiology department of the Rennes University Hospital (CHU de Rennes), in France. Participants were enrolled in 6 French hospitals located in Rennes, Saint Pierre  In accordance with the current guidelines [1,2], BS was diagnosed when a coved ST-segment elevation (! 0.2 mV) was recorded in at least one right precordial lead (V1 and/or V2) located in the 2 nd , 3 rd or 4 th intercostal space, in the presence or absence of sodium-channelblocking agent.
In order to characterize populations with different levels of risk, patients were classified as symptomatic and asymptomatic, based on their medical history. Twenty-four patients presented documented symptoms of ventricular origin: syncope (50%), cardiac arrest (41.7%), dizziness (12.5%) and, less frequently, palpitations and nocturnal convulsions (4.2%). The remaining 81 patients were considered as asymptomatic.
Participant ages ranged from 19 to 74 years old (45.17 ± 13.62 years old) and 76.2% were males. ICDs had been implanted in 18 of 81 (22.2%) asymptomatic patients, based on a positive EPS (Electrophysiological Study) test, whereas all symptomatic patients had ICDs implanted. Among 76 patients (19 were symptomatic) in whom genetic analysis was performed, an SCN5A mutation was found in 27 (35.5%). Table 1 summarizes the clinical characteristics of patients included in the study.
The cardiac response to exercise is influenced by the complex interaction of many factors including age, gender, physical conditioning, sympathetic drive, baroreceptor reflexes and venous return [26]. Nevertheless, since no significant differences in age, gender and SCN5A-mutation presence between symptomatic and asymptomatic groups were noted (p-value>0.05), similar baseline characteristics were assumed between populations.

Data acquisition and test
Patients underwent a triangular exercise test recommended by the American Heart Association [27] where the load was increased until it reached the 80% of the theoretical maximum heart rate of each patient, defined by the formula MHR = 220 − age [28]. The test was performed in a cyclo ergometer (Ergoline 900 Egamed, Piestany, Slovakia) and divided in the following phases: • Exercise phase: • Warm-up phase: for men, initial load of 50 watts (W); for women, initial load of 30 W, both for 2 minutes.
• Incremental exercise phase: for men, initial load of 80 W for 2 minutes and then incrementing 20 W every 2 minutes; for women, initial load of 50 W increasing load 20 W every 2 minutes.
• Recovery phase: • Active recovery: for men, fixed load of 50 W; for women, fixed load of 30 W, both for 3 minutes. • Passive recovery: total cessation of effort for 3 minutes.
The standard 12-lead ECG recordings sampled at 1000 Hz from each patient were collected and analyzed by the central board (Rennes University Hospital). Fig 1 illustrates the global methodology proposed in this paper in order to differentiate symptomatic and asymptomatic BS patients. This methodology is based on a general machine learning approach built from the following four main steps:

Proposed classification methodology
A. Feature extraction. The standard 12-lead ECG signals acquired during exercise testing were analyzed to detect each QRS complex and extract the RR and ECG-Derived Respiration (EDR) series of each patient. A time-frequency (TF) method based on the smoothed M , respectively, where N b refers to the 160 observations after class balancing (79 symptomatic and 81 asymptomatic samples). C) Feature selection, which starts by randomly defining patient subsets for training, (N tr , 75% of patients, 59 symptomatic and 60 asymptomatic) and testing (N te , the rest of patients, 20 symptomatic and 21 asymptomatic), followed by the estimation of a minimal feature dimension M w < M, that maximizes classification performance, using filtering and wrapper methods. D) The final step is dedicated to classification and performance evaluation. pseudo Wigner-Ville distribution (SPWVD) of RR series that adapts frequency bands to respiratory information resulting from EDR signals was applied to estimate the evolution of different spectral HRV markers. Estimated features include the LF, LF nu , HF, HF nu and LF HF mean values at different time periods of the exercise test. The output of this step is matrix R N M , which contains the calculated raw features, with M = 60 HRV markers for the N = 105 patients available on the whole database.
B. Feature conditioning. In order to handle the impact of markers measured at different scales, all features of R N M were standardized, leading to matrix F N M . Then, in order to reduce the effect of imbalanced classes, 55 synthetic symptomatic patients were generated and included in the analysis by a class balancing approach, resulting in matrix F N b M , with N b = 160 observations. C. Feature selection. After balancing all standardized HRV features, F N b M was divided in a training subset F Ntr M (Ntr = 119, 75% randomly selected patients) and the remaining testing subset F Nte M (Nte = 41). Then, a two-step feature selection process including a filter and a wrapper method was applied in order to capture the most relevant HRV features. In Fig  1C,  A more detailed description of each step is presented in the following sections.

Feature extraction
RR series extraction. From the standard 12-lead ECG recordings of each patient, RRinterval and R-peak amplitude series were extracted by using a noise-robust wavelet-based algorithm for QRS complex detection and subsequent R-peak location [29]. After performing manual corrections when necessary, a cubic-spline interpolation was applied to RR-interval time series, to obtain uniformly sampled data at a rate of 4 Hz. A representative example of RR series observed during each phase of the exercise test is shown in Fig 2. Time-frequency HRV analysis. According to the Task Force on HRV [30], classic spectral HRV indices require stationary data to provide accurate estimates of ANS modulation. As illustrated by an increasing heart rate (decreasing RR series) in proportion to exercise workload in Fig 2, given that signals on a physical stress test are typically non-stationary, spectral characteristics associated with HRV were analyzed using a time-frequency (TF) approach.
First, in order to remove the very low frequency component, RR series were high-pass filtered at 0.03 Hz with a 4 th order Butterworth filter applied in both forward and backward directions so as to remove phase distortion. Then, a smoothed pseudo Wigner-Ville distribution (SPWVD) transform from the Time-Frequency toolbox [31] was employed since it has proved its usefulness for the analysis of cardiovascular signals [32].
The Wigner Ville distribution is a quadratic time-frequency method defined as the Fourier transform of the instantaneous autocorrelation function [33]. However, since it is affected by significant interference terms, the SPWVD introduces a smoothing kernel function C(τ, υ), defined in Costa et al [34], that attenuates interferences while maintaining a suitable time-frequency resolution. Being A RR (τ, υ) the Ambiguity Function of the RR series, x RR (t), the SPWVD is defined as: Kernel parameters were adjusted to υ 0 = 0.06 and τ 0 = 0.03, obtaining temporal and spectral resolutions of 16.7 seconds and 0.033 Hz, respectively. Among all the analyzed combinations, this one led to the most efficient interference terms cancellation for the lowest TF filtering. Then, HRV was measured as the total power in LF and HF bands (noted as LFb and HFb), obtained from the SPWVD: Assuming that sympathetic activity always lies within the standard LF band, this band was fixed between 0.04 and 0.15 Hz for the whole stress test. However, the total power in the HF band captures parasympathetic activity, closely related to respiratory sinus arrhythmia (RSA). Since respiratory frequency during exertion is not restricted to the classic HF band (0.15-0.4 Hz) and can increase up to 0.7 Hz, HRV analysis within the standard frequency band would lead to unreliable measures of the parasympathetic activity. In order to overcome this limitation, we defined a time-varying HF band, based on an estimation of the respiratory activity from the ECG signal, by applying an ECG-Derived Respiration (EDR) method [35].
Time-varying respiratory frequency estimation. The applied EDR method estimates respiratory information from the amplitude modulation of R-wave peaks [36]. Cubic-spline interpolation at a rate of 4 Hz was also applied to the obtained series. Then, a band-pass 4 th order Butterworth filter between 0.15 and 0.7 Hz was applied in both forward and backward directions to remove frequencies out of the respiratory range. The same SPWVD transform used for RR series was then applied to EDR filtered signals to estimate the instantaneous respiratory frequency.
The simplest estimation method consists in finding the frequency presenting the largest peak in the spectrum at each time instantf ðtÞ. However, in order to avoid spurious peak detections, for each time instant t k , the search interval was limited to frequencies between 2δ Hz, centered around a reference frequency This reference frequency was defined as an exponential average of previous estimates: where β is the forgetting factor. As in [37], a value of β = 0.7 was used, based on real respiratory patterns during exercise testing and δ = 0.01, since respiratory frequency variations are not supposed to be faster than 0.01 Hz per 0.25 s. Moreover, to reduce the risk of spurious frequency detections in the initialization of the reference frequency, the first instantaneous respiratory frequency f r (t 0 ) was selected within the standard HF band (0.15 − 0.4 Hz).
Once the estimated respiratory frequency series f r (t) was obtained, the time-varying HF band for HRV analysis was defined as HFb(t) = [f r (t) − 0.125, f r (t) + 0.125] Hz, with t covering the whole test.
Final HRV features extraction. Unlike classic autonomic indices, SPWVD leads to timefrequency HRV estimates that are indeed time series that vary during the exercise testing. These markers, accounting for the sympathetic and parasympathetic influences of the ANS on heart rate, were normalized and expressed as percentages of the total power (TP), defined as the sum of both spectral bands (TP(t) = LF(t) + HF(t)), leading to the time series LF nu (t) and HF nu (t): Since each test differed in the incremental exercise phase duration and the shortest case lasted less than 5 minutes, in order to compare the same time periods between groups of patients, only the first 3 minutes of incremental exertion (EX1, EX2 and EX3), as well as the last minute before peak effort (PE), were assessed. In addition, the entire warm-up (WU1 and WU2) and both active (AR1, AR2 and AR3) and passive recovery (PR1, PR2 and PR3) phases were compared between symptomatic and asymptomatic patients. Fig 2 displays the analyzed periods along the different phases of the exercise test on a representative example of RR-interval series, indicating the peak effort instant. From each HRV series, 12 intra-patient 1-minute means were calculated, leading to 60 features per patient that, after being compared between populations, were included in the multivariate classification approach described herein. Data in S1 Table include the resulting HRV features for the whole clinical series.
Although the autonomic response to exercise significantly depends on test conditions [38], since no statistically significant differences in workload at peak effort were observed between populations (symptomatic: 175.4 ± 56.6 W; asymptomatic: 175.1 ± 55.3 W; p-value: 0.957), similar exercise intensities in both groups of patients were assumed.
Statistical comparison. Comparisons between symptomatic and asymptomatic patients at each analyzed minute of the physical stress test were evaluated by Mann-Whitney U nonparametric tests. In order to compare the last minute of exertion and recovery, all patients had to be synchronized with respect to the peak effort instant. The analysis was made using the commercially available software MatLab (Mathworks Inc., MI, USA) and setting the level of significance at p < 0.05.

Feature conditioning
All features extracted from HRV analysis were considered as candidates for the construction of a model classifying symptomatic and asymptomatic BS patients. The initial feature subset R N M was composed of N = 105 observations (24 symptomatic and 81 asymptomatic patients) by M = 60 features (5 HRV markers for 12 analyzed minutes of test). In order to equalize the contribution of all features to multivariate analysis, each raw HRV marker j for each patient i was standardized as follows: where μ j is the mean and σ j the standard deviation of a specific feature j, taking into account the data from all patients i = 1, Á Á Á, N. The new standardized dataset is defined in matrix F N M . Then, to attenuate the impact of imbalanced classes, synthetic symptomatic samples were generated by applying the ADASYN approach [39]. Since this method randomly chooses examples from the minority class to generate new samples, the algorithm was applied 50 times and the mean from all realizations was kept as the final balanced dataset F N b M , where N b = 160 observations. Fig 3 illustrates the feature conditioning process.

Feature selection
To reduce the number of attributes included in the model so as to decrease its computational cost, we applied a two-step feature selection approach that identified the most relevant features in distinguishing between symptomatic and asymptomatic patients. As previously mentioned, Fig 1 specifies the methodology followed for feature extraction and posterior feature selection, only applied to a randomly selected sample of 75% of the feature database (F Ntr M : training subset). The remaining 25% (F Nte M : testing subset) was then used for model validation.
ReliefF filter method. The first step in the feature selection process was a simple filter method based on the ReliefF algorithm [40]. Since this approach ignores the effects of attributes on classification, it can rapidly remove some irrelevant and redundant features.
The algorithm estimates feature weights W according to their capability of distinguishing between data from different classes, here symptomatic and asymptomatic patients, based on a k-Nearest Neighbors (k-NN) approach. The extension of the original Relief method used in this study, ReliefF, not only deals with multiclass problems but it is also more robust with incomplete and noisy data.
Hence, the algorithm assigns a relevance weight ranged from −1 to 1 to each feature, with large positive weights allocated to significant attributes. However, it should be noted that, since this method is based on a k-NN approach, feature weights usually depend on k. For small values of k, the estimates can be unreliable for noisy data; while for k values comparable with  the number of observations, the algorithm can fail to find significant attributes. Thus, ReliefF was computed for k = 10, Á Á Á, 19 and W was defined as the average of all weights.
Moreover, since HRV feature values significantly vary over patients, a bootstrap technique was applied [41]. The algorithm was run 50 times on different randomly chosen subsets including 60% of the training dataset F Ntr M (35 symptomatic and 36 asymptomatic observations), here represented as F Ntr i M . Then, the relevance of each feature was obtained as the median of the 50 realizations analyzed. Fig 4 displays the methodology followed to select M f as the 75% most relevant attributes (45 best ranked out of 60 features) that were kept for further analysis.
LDA-based wrapper method. For final feature selection, a second step was applied on the reduced subset of attributes resulting from the previous stage. It consisted in a wrapper algorithm with both forward and backward search strategies (floating method), using a Linear Discriminant Analysis (LDA) classifier as a black box. Since this approach is based on classification performance, the final subset is only optimized for this particular classifier [42]. Fig 5 represents the wrapper feature selection process in more detail. As in the previous step, it was repeated 50 times on different randomly chosen subsets of training data F Ntr i M f . Those features appearing more than L times, among the 50 realizations, formed the final subset F Ntr M w . The value of L was optimized, based on performance metrics, so as to find the best selection of features M w leading to the finest classifier distinguishing between symptomatic and asymptomatic BS patients.

LDA-based classifier
After selecting the best set of features, an LDA classifier [43] was implemented using a 5-fold cross-validation approach in order to reduce classification error. This technique divides the entire training subset F Ntr M w into 5 blocks where each classifier is firstly trained on 4 portions and then tested on the 5 th block. This is performed for the 5 different possible combinations of blocks for training/testing so the outputs of each solution are then averaged.
In addition, to estimate the mean performance variability of the classifier when applied to testing data, 5-fold cross-validation was run 10 times on differently divided subsets of training data. Fig 6 illustrates the steps followed to train and test the LDA classifier. As specified in previous sections, 75% of the data F Ntr M w were used for training, and the remaining 25% F Nte M w for testing.

Performance evaluation
Model performance evaluation was based on the resulting confusion matrix, which specifies the number of true positives (TP), true negatives (TN), false positives (FP) and false negatives (FN), when comparing true and predicted labels and considering symptomatic patients as positives.
First, the AUC or area under the ROC (Receiver Operating Characteristic) curve was computed to quantify the classifier performance. This measure was also assessed in order to optimize the threshold L that led to the best performing classifier after wrapper feature selection.
Moreover, classical sensitivity (Se = TP/(TP + FN)) and specificity (Sp = TN/(TN + FP)) measures, associated with the optimal operating point in the ROC curve, were calculated to quantify the classifier capability of correctly detecting symptomatic and asymptomatic patients, respectively.

Single-patient representation
Given that all patients presented similar tendencies in RR and HRV series, the following representative example illustrates changes induced by exercise testing on those time series involved in HRV features extraction.
The upper panel of Fig 7 displays an example of EDR series. Below, the SPWVD spectral power of respiration is shown, along with its estimated instantaneous respiratory frequency f r (t), represented with a dashed red line. Note that, in this example, as the patient approaches the peak effort (dashed vertical line), respiratory frequency exceeds the standard HF band upper limit of 0.4 Hz.
Based on this respiratory information, the time-varying patient-specific HF band was identified. Fig 8 displays the RR series and its associated SPWVD spectral power for the same patient, where the LF and respiration-centered HF bands are represented in dashed white lines (second panel).
Although certain stationarity can be observed during the first two minutes corresponding to the warm-up phase, the RR series displays a significant non-stationarity represented by a progressive decrease during incremental exercise and a continuous increase at recovery. As previously reported by [44], the cardiovascular response properly balanced the performed physical activity intensity. Moreover, it should be noted the non-negligible power centered in the HF band as a result of respiration. Since HF exceeds the standard band, the use of classic spectral limits would have led to unreliable measures of the parasympathetic activity in this patient.
The third panel shows the time series LF(t) and HF(t) extracted from TF analysis, then normalized and expressed as percentages of the total power. Finally, the last panel displays the obtained HF nu (t) series, where the 1-minute window to calculate the mean value of HF nu during the first minute of active recovery (HF AR1 nu ) is indicated. The results show a progressive increase in HF nu (t) during exercise, as well as a decrease at recovery. Although not represented in Fig 8, due to the normalization step applied to LF(t) and HF(t), the complementary effect was observed for LF nu (t).

Inter-group comparison
HRV features extracted during exercise and recovery were compared between symptomatic and asymptomatic patients. Apart from warm-up and recovery phases, since exercise duration was highly variable among subjects, only the first 3 minutes and the last minute of incremental exertion were analyzed.
During the second minute of incremental exercise, statistically significant differences were found in mean normalized HF (HF EX2 nu , p = 0.041) and thus in LF EX2 nu . Symptomatic patients showed an increased HF nu , and a reduced LF nu , in this time period with respect to asymptomatic patients. However, no significant differences between groups were noted after exertion, during active or passive recovery. Table 2 summarizes the mean, standard deviation and associated p-values for LF, LF nu , HF, HF nu and LF HF obtained during the warm-up phase, the first 3 minutes of incremental exercise and the last minute before peak effort.

Classification
After feature conditioning, feature selection was performed by a two-step approach. Filter selection was followed by the repeated application of a wrapper method to the selected features. The final subset contained those features appearing more than a specific number of times L among realizations, optimized based on the performance metric AUC. Fig 9 displays the mean and standard deviation of the AUC associated with each value of L.
Based on these results, the final classifier placed the threshold at L = 18, where a maximum AUC using the minimum number of features was found. When only those parameters   appearing more than 18 times were kept, the final subset contained 22 features (listed in Table 3), leading to an AUC = 0.92 ± 0.01, Se = 0.91 ± 0.06 and Sp = 0.90 ± 0.05. Fig 10 displays the mean ROC curves resulting from each 5-fold cross-validation, as well as the global ROC curve and its optimal operating point, for the proposed classifier.

Discussion
In this study, the autonomic response to exercise testing was analyzed in 105 BS patients. Although many previous studies have already assessed the autonomic function in BS, only a few have exploited the potential of exertion to predict cardiac events [22][23][24][25]. Moreover, since symptoms have been related to VF occurrence on this population [5,6], our aim was to compare, according to symptomatic status, the time-varying HRV changes induced by effort and subsequent recovery. To our knowledge, this work presents the first time-frequency HRV analysis under conditions of physical activity in BS. Furthermore, due to the complexity of distinguishing between symptomatic and asymptomatic patients by means of univariate analyses, a multivariate classifier based on the combination of the extracted HRV features was proposed.
Most reported approaches assessing the autonomic response of BS patients are based on measurements obtained from 24-hour recordings that lead to controversial conclusions. In a clinical series of 17 patients with BS, of whom 10 were asymptomatic subjects with Brugada ECG, and 45 controls, Krittayaphong et al [12] concluded that BS patients presented a lower HRV and a lower vagal tone at night compared to controls, as well as lower diurnal and higher overnight heart rates compared to asymptomatic subjects and controls. Likewise, Hermida et al [13] reported a significantly lower HRV at night in 21 symptomatic, with respect to 26 asymptomatic BS patients. Pierre et al [14] also asserted in a clinical series of 46 BS patients and 46 controls that HRV in the first group was significantly lower with respect to healthy subjects. Tokuyama et al [15] results showed a significantly lower HRV in BS patients when analyzing a series of 12 symptomatic, 17 asymptomatic and 16 healthy individuals. The results also reflected a significant reduction in both sympathetic and parasympathetic tones in BS patients, as well as a decreased circadian variation of the autonomic function over 24 hours, with respect to controls. Kostopoulou et al [17] examined autonomic disorders in 20 patients with BS and 20 age-matched controls. In that case, HRV analysis did not reveal any significant difference between groups, but a high susceptibility to vasovagal syncope was observed in BS patients, possibly being a disease-related symptom. Nakazawa et al [18] analyzed, using a 24-hour continuous ECG monitoring, the autonomic properties of 27 BS patients (10 of them had a history of VF and 17 did not) and of 26 healthy subjects, finding higher vagal and reduced sympathetic tones in symptomatic BS patients. Likewise, in a recent work from our group where the 24-hour Holter recordings from 118 BS patients were analyzed, symptomatic subjects showed an increased parasympathetic activity during both daytime and nighttime [16] Subramanian et al [23] proved the usefulness of some electrocardiographic markers extracted during exercise testing for risk stratification in BS. Moreover, our group has recently reported significant differences in heart rate complexity between symptomatic and asymptomatic patients, during periods of recovery after exertion [24]. Nevertheless, Amin et al [22] published the first work measuring the autonomic function of 50 BS patients and 35 controls during exercise, finding a higher parasympathetic reactivation during early recovery in patients with prior VF events. Likewise, Makimoto et al [21] analyzed the autonomic function of 93 BS patients and 102 controls during recovery from treadmill exercise testing. They studied parasympathetic reactivation by computing the Heart Rate Recovery (HRR) after peak exercise, concluding that a higher vagal activity was related to the occurrence of cardiac events in BS.
In our study, as illustrated in Fig 8, all patients displayed a progressive increase in the mean normalized HF (HF nu ) during incremental exercise, as well as an HF nu decrease at recovery. Although data on direct sympathetic nerve recording and plasma catecholamines measurements have reported that a decreased parasympathetic and an increased sympathetic activity play a major role in the autonomic response to exercise, many studies on cardiac autonomic function based on HRV analysis have failed to represent this response, even in healthy subjects [38]. Indeed, the LF component does not provide an index of sympathetic tone but rather reflects a complex interplay among many factors including the sympathetic and parasympathetic contributions to ANS. Similarly, just as parasympathetic neural activity influences LF values, sympathetic activation also modulates the HF component [45]. Moreover, since an increasing power can be noted at the HF band center corresponding to the respiratory frequency (second panel in Fig 8), the gradual increase observed in HF nu during incremental exercise may be significantly influenced by respiration. Thus, LF and HF indices should not be analyzed as accurate representations of, respectively, the sympathetic and parasympathetic tones. They should be interpreted as estimates of the autonomic function that may capture relevant tendencies in HR modulation, potentially useful for the detection of differences between BS patients at different levels of risk.
According to the inter-group comparison of HRV markers, statistically significant differences were observed during the second minute of incremental exercise in HF nu , and thus LF nu , between groups. Since no significant differences were found in terms of spectral power at the respiratory frequency, HF nu differences between populations might mostly be due to vagal activity. Thus, as previously reported [15,16,18], symptomatic patients seem to experience an increased parasympathetic modulation with respect to asymptomatic patients, supporting the idea that higher vagal responses could be related to a worse prognosis in BS. These results may be explained by the dysfunction on presynaptic norepinephrine (NE) recycling and the reduction in the concentration of NE at the synaptic cleft found on previous works based on positron emission tomography on BS patients [8][9][10][11].
The lack of significant results in univariate analysis reveals the difficulty of distinguishing between symptomatic and asymptomatic groups. Therefore, a multivariate approach following a step-based machine learning method was designed in order to improve classification performance. The proposed solution significantly reduced the final subset of features included in the predictive LDA-based model to one third of the total amount of HRV features, leading to a mean AUC of 92.1%. First, a filter feature selection method discarded the least relevant and most redundant features, holding the 75% of the initial features subset, to which the LDAbased wrapper algorithm was applied. On the one hand, the results after filtering show that all autonomic markers during the last minutes of incremental exercise and recovery were kept, evidencing the relevance of these test segments in classifying BS patients. On the other hand, although HF EX2 nu led to significant results in univariate analysis, the applied filter method identified this marker as a redundant feature and only kept LF EX2 nu for further analysis. Since classification performance significantly depends on the number of chosen features after wrapper feature selection, the algorithm was optimized to obtain the best AUC using the minimum number of features. Thus, when selecting only those features appearing more than 18 times after wrapper application, an optimal classifier containing 22 parameters was implemented. Among the final subset of features, only one LF HF marker was kept, acquired during the second minute of active recovery. The remaining parameters equally belonged to LF and LF nu or HF and HF nu measures. Regarding test phases, only one feature came from the warm-up phase (HF WU1 nu ) and other 5 markers were acquired during the passive recovery stage. Most parameters were measured during incremental exercise and active recovery, and more specifically during the last minutes of both phases.
Recent studies have also proposed prediction models for VA risk stratification in BS patients using non-invasive parameters [46]. However, our model outperformed previous approaches, evidencing the interest of analyzing HRV features during exercise testing to better understand VF risk in this population. Indeed, the results from our previous study, where the classification potential of these markers was already presented [47], were also improved by enlarging the clinical series under study and applying a more robust feature selection approach.
This study presents some relevant limitations that should be noted. The clinical value of autonomic features can only be proved if a close relationship between HRV markers and ventricular events is established. Since no VF was induced during the test, HRV variations in symptomatic patients cannot be directly related to this phenomenon. Moreover, a synthetic oversampling approach was applied in order to overcome complications found when learning from imbalanced datasets. Thus, the obtained results should be validated by enlarging the studied population. Finally, since some asymptomatic patients may develop symptoms in the future and thus present high-risk patterns during the analyzed recordings, a more suitable clinical database for risk stratification should include follow-up information. Thereby, autonomic changes could have been related to the probability of developing symptoms rather than to the identification of a high-risk group having already suffered these symptoms.
Nonetheless, previous studies have shown the need of new autonomic markers with higher predictive values, such as those here presented, to better stratify risk in patients suffering from Brugada syndrome. According to international guidelines [1], ICD implantation is recommended in BS patients being survivors of a cardiac arrest and/or having documented spontaneous sustained ventricular tachycardias (class I) and can be useful in patients with a spontaneous diagnostic type 1 Brugada-like ECG pattern having a history of syncope caused by ventricular arrhythmias (class IIa). However, the decision of implanting an ICD on asymptomatic subjects is still contentious, even if they represent around the 60% of diagnosed patients. Thus, the proposed model is presented as a potential instrument to better identify those asymptomatic BS patients at high risk who may benefit from an ICD implantation. Moreover, the proposed model might also be used for processing HRV data acquired from ICDs on implanted BS patients, in order to control their risk of VF occurrence during followup.

Conclusions
In this study, the autonomic function of 105 BS patients who underwent a standardized physical stress test was analyzed so as to characterize symptomatic and asymptomatic populations. Based on the hypothesis that changes in the ANS induced by exercise testing could improve prognosis interpretation, a classifier capable of identifying patients at high risk was then designed.
First, the extracted time-varying HRV features were compared between populations. Statistically significant differences were found in LF nu and HF nu during incremental exercise, suggesting that symptomatic patients seem to experience an increased vagal function with respect to asymptomatic BS patients.
Then, a predictive model based on a two-step feature selection strategy identified the most discriminant HRV features to distinguish symptomatic and asymptomatic patients. Despite the difficulty in finding differences between these populations, classification results show the potential of autonomic markers when identifying symptoms in BS.
Although the present study presents some limitations and is based on a relatively small population of BS patients, the results indicate important trends of clinical relevance that could be useful for risk stratification in asymptomatic patients for whom the decision to implant a cardioverter defibrillator is complex and controversial.
Supporting information S1