Next Article in Journal
Fast Compression of MCMC Output
Previous Article in Journal
Anisotropic Strange Star in 5D Einstein-Gauss-Bonnet Gravity
Previous Article in Special Issue
Complexity of Body Movements during Sleep in Children with Autism Spectrum Disorder
 
 
Erratum published on 20 October 2021, see Entropy 2021, 23(11), 1375.
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Bispectral Analysis of Heart Rate Variability to Characterize and Help Diagnose Pediatric Sleep Apnea

by
Adrián Martín-Montero
1,*,
Gonzalo C. Gutiérrez-Tobal
1,2,
David Gozal
3,
Verónica Barroso-García
1,2,
Daniel Álvarez
1,2,
Félix del Campo
1,2,4,
Leila Kheirandish-Gozal
3 and
Roberto Hornero
1,2
1
Biomedical Engineering Group, University of Valladolid, 47002 Valladolid, Spain
2
CIBER-BBN, Centro de Investigación Biomédica en Red en Bioingeniería, Biomateriales y Nanomedicina, 28029 Madrid, Spain
3
Department of Child Health and the Child Health Research Institute, The University of Missouri School of Medicine, Columbia, MO 65212, USA
4
Sleep-Ventilation Unit, Pneumology Service, Río Hortega University Hospital, 47012 Valladolid, Spain
*
Author to whom correspondence should be addressed.
Entropy 2021, 23(8), 1016; https://doi.org/10.3390/e23081016
Submission received: 5 July 2021 / Revised: 28 July 2021 / Accepted: 3 August 2021 / Published: 6 August 2021
(This article belongs to the Special Issue Entropy and Sleep Disorders II)

Abstract

:
Pediatric obstructive sleep apnea (OSA) is a breathing disorder that alters heart rate variability (HRV) dynamics during sleep. HRV in children is commonly assessed through conventional spectral analysis. However, bispectral analysis provides both linearity and stationarity information and has not been applied to the assessment of HRV in pediatric OSA. Here, this work aimed to assess HRV using bispectral analysis in children with OSA for signal characterization and diagnostic purposes in two large pediatric databases (0–13 years). The first database (training set) was composed of 981 overnight ECG recordings obtained during polysomnography. The second database (test set) was a subset of the Childhood Adenotonsillectomy Trial database (757 children). We characterized three bispectral regions based on the classic HRV frequency ranges (very low frequency: 0–0.04 Hz; low frequency: 0.04–0.15 Hz; and high frequency: 0.15–0.40 Hz), as well as three OSA-specific frequency ranges obtained in recent studies (BW1: 0.001–0.005 Hz; BW2: 0.028–0.074 Hz; BWRes: a subject-adaptive respiratory region). In each region, up to 14 bispectral features were computed. The fast correlation-based filter was applied to the features obtained from the classic and OSA-specific regions, showing complementary information regarding OSA alterations in HRV. This information was then used to train multi-layer perceptron (MLP) neural networks aimed at automatically detecting pediatric OSA using three clinically defined severity classifiers. Both classic and OSA-specific MLP models showed high and similar accuracy (Acc) and areas under the receiver operating characteristic curve (AUCs) for moderate (classic regions: Acc = 81.0%, AUC = 0.774; OSA-specific regions: Acc = 81.0%, AUC = 0.791) and severe (classic regions: Acc = 91.7%, AUC = 0.847; OSA-specific regions: Acc = 89.3%, AUC = 0.841) OSA levels. Thus, the current findings highlight the usefulness of bispectral analysis on HRV to characterize and diagnose pediatric OSA.

1. Introduction

Obstructive sleep apnea (OSA) is a common respiratory disorder affecting up to 5% of the general pediatric population [1]. OSA is characterized by the occurrence of either total upper airway obstruction (apnea events) and/or events of significant airflow reduction (hypopnea events) during sleep, leading to decreased blood oxygen saturation and/or sleep fragmentation [2,3]. As a result, pediatric OSA has been associated with several adverse cognitive and cardiovascular consequences, such as learning deficits and reduced academic performance [4,5], as well as alterations in autonomic regulation and vasomotor tone, increases in systemic and pulmonary vascular blood pressure, nocturnal cardiac strain and both left and right ventricular hypertrophy [1,3]. The increased cognitive and cardiovascular risks obviously threaten the long-term cardiovascular health [1,3,6] and academic potential of children [7], such that early detection and treatment of pediatric OSA are essential.
Nocturnal laboratory-based polysomnography (PSG) is considered the standard diagnosis technique for pediatric OSA. This test not only allows for the detection of the presence or absence of pediatric OSA but also enables estimates of OSA severity [8,9]. During a PSG, a set of leads is placed on the child’s body to register biological signals, including an electrocardiogram (ECG), electroencephalography, oximetry (SpO2) or airflow, among others [8]. Then, medical experts evaluate these signals following well-established rules to extract indices of respiratory disturbance [10]. Among those indices, the common choice to illustrate and report OSA severity is the Apnea–Hypopnea Index (AHI), which reflects the number of total apneic and hypopneic events per hour (e/h) of sleep [8,10]. Despite the usefulness of PSG to diagnose pediatric OSA, the procedure needs a specialized laboratory and is resource-intensive and time-consuming. Moreover, the high number of sensors connected to the child along with having to sleep outside of home makes PSG also especially uncomfortable and inconvenient for patients and caretakers [11,12]. These drawbacks have motivated the search for alternatives to diagnose pediatric OSA and to study its consequences while reducing the number of signals required for the diagnosis [11,12,13]. A recently published systematic review [14] analyzed and conducted a meta-analysis on machine learning techniques employed to automatically diagnose pediatric OSA. Among the shortcomings identified in the literature, the authors highlighted that most of the studies were exclusively based on the analysis of SpO2 signals. Thus, there exists a lack of studies performing machine learning techniques using other physiological signals as alternative approaches to the gold standard PSG [14].
Evaluation of the adverse cardiovascular implications of pediatric OSA has highlighted autonomic dysfunction as the main reason for cardiac alterations, since the autonomic nervous system (ANS) plays a major role in cardiovascular system regulation [15,16,17]. In this sense, analysis of heart rate variability (HRV) has gained obvious relevance as it reflects the ANS state [17,18]. HRV non-invasively assesses the variations in the heart rate over time, meaning it can be used to analyze cardiovascular modulation due to OSA [18,19,20,21]. As a result, a variety of studies have evaluated HRV in the context of pediatric OSA [16,17,22,23,24,25,26,27,28]. The vast majority of these studies focused on HRV analysis in the temporal or frequency domain. When using these traditional analysis approaches, nonlinearity and non-Gaussianity are ignored [29]. Although HRV signals are usually nonlinear and non-Gaussian in nature [29,30], under some specific conditions, HRV can show dynamics where these nonlinearities are not always present [31]. Nevertheless, in the context of pediatric OSA, the nonlinear dynamics of HRV signals can be increased during sleep [21,32], as well as by cardiac alterations due to OSA [20,21,30,32]. Furthermore, in the present study, the presence of HRV nonlinear dynamics in the pediatric OSA context has been demonstrated (see Appendix C). Despite this, only two of the previous studies on pediatric OSA performed non-lineal analysis of HRV signals using Poincaré plots [17,33]. Notwithstanding, application of bispectral analysis of HRV would likely constitute a further advance, as it allows for reflection of not only nonlinear behaviors but also non-Gaussian and non-stationary events, all the while retaining phase and amplitude information, and being more immune to noise [29]. The unique advantages conferred by the bispectrum analysis properties have led to its application in HRV for some purposes such as evaluation of cardiac state [34,35,36], congestive heart failure [37], major depression [38] or even OSA diagnosis in adults [30].
OSA was initially described in adults, and its cardiac implications and effects on HRV have been extensively studied in recent decades [39,40,41]. Nevertheless, pediatric OSA presents differentiating etiological, diagnostic and treatment considerations [5]. It has led to an increasing interest in the study of pediatric OSA to gain insights into the mechanism underlying pediatric OSA pathogenesis in recent years [3,5]. The advances in the study of pediatric OSA have demonstrated that the long-term effects of the disease in the child population can be avoided with a timely treatment [11]. Accordingly, more stringent rules were defined to diagnose pediatric OSA [10]. In this sense, while an adult is considered as an OSA patient when presenting an AHI = 5 e/h, the presence of an AHI = 1 e/h in children is enough to diagnose pediatric OSA. Similarly, an adult patient with an AHI = 10 e/h is considered as mild OSA, whilst the same AHI in children is considered as severe OSA [42]. Furthermore, an apneic event is scored in adults when it lasts over 10 s, while in children, 6 s is enough to score a respiratory event as apneic [10]. These distinctions between children and adults also produce a marked difference in the HRV alterations due to OSA and consequently affect the HRV bispectrum in a distinct way, evidencing the necessity of independent HRV pediatric study. However, to date, HRV bispectral analysis has never been evaluated in the pediatric OSA context.
Of note, previous bispectral analysis study in adults focused on OSA [30], as well as some of the aforementioned studies [34,35], considered HRV bispectral analysis that was restricted to the non-redundant bispectral region. This region, as a result of the bispectrum symmetric properties, is known to completely define the overall information contained in the bispectrum [29]. Other studies analyzed particular bispectral regions defined based on the classic frequency ranges of HRV analysis [36,37,38]. However, in a study using regular spectral analysis [27], we recently showed that there exist OSA-specific frequency ranges that allow for better characterization of the alterations occurring in HRV in the context of pediatric OSA.
Based on these considerations, we hypothesized that bispectral analysis of HRV contains novel and useful information to characterize and diagnose pediatric OSA. Consequently, the main objectives of this study were to perform, for the first time in the field of pediatric OSA, a characterization and evaluation of bispectral regions bounded with classic and OSA-specific HRV frequency ranges. Therefore, the main novelty of this study is the analysis of HRV bispectral regions defined by classic and OSA-specific frequency ranges to characterize and diagnose pediatric OSA. Furthermore, we propose a novel bispectral parameter, which is analyzed here for the first time.

2. Subjects and Signals under Study

The population included in this study was the same as that in a previous work [27] and was composed of 1738 children between 0 and 13 years of age. This cohort consists of two large databases. The first one was generated in the Pediatric Sleep Unit at the University of Chicago (UofC) Medicine Comer Children’s Hospital (Chicago, IL, USA) and involved 981 children referred for a nocturnal laboratory-based PSG due to clinical symptoms suggestive of OSA. The second database was a subset of 757 PSGs performed during the multicenter Childhood Adenotonsillectomy Trial (CHAT) study [43,44].
In regard to the UofC dataset, the Ethics Committee of the UofC Medicine Comer Children’s Hospital approved the protocol (#11-0268-AM017, #09-115-B-AM031 and #IRB14–1241), and the study was conducted in accordance with the Declaration of Helsinki. Additionally, before the acquisition of the ECGs during the PSG, informed consent from all children caretakers was acquired. For comparison purposes, the training–test division carried out was the same as that in our preceding work [27], meaning the 981 children from the UofC sample served as the training set.
Regarding the CHAT set (clinical trial identifier: NCT00560859), information about the rationale, study design and methodological aspects can be perused in the published literature [43,44]. Initially, there were 779 pediatric PSG recordings, but 22 subjects were removed from the study because they did not fulfill the signal pre-processing protocol detailed below. Accordingly, the remaining 757 children were retained as the test set.
The PSG studies from both databases were evaluated by pediatric sleep experts from the different centers who diagnosed the children following the scoring rules established by the American Academy of Sleep Medicine [10]. Subsequently, the AHI was obtained and used to determine OSA severity. In accordance with previous pediatric OSA studies [11,27,42,45,46], we defined four OSA severity groups based on three common AHI cutoff points (1, 5 and 10 e/h). Thus, the severity groups in this study included: no-OSA (AHI < 1 e/h), mild OSA (1 ≤ AHI < 5 e/h), moderate OSA (5 ≤ AHI < 10 e/h) and severe OSA (AHI ≥ 10 e/h) groups. Clinical and demographic data of the children under study are shown in Table 1.
In order to conduct the HRV analysis, ECGs from both datasets were equally pre-processed [27]. First, we removed the initial and last 15 min of the ECG recordings to prevent the signals from periods of early and final artifacts [27]. After this removal, we assessed that, at least, 3 h of sleep recording was available for each child [46,47]. The next step consisted in R peak detection. To carry out this process, the algorithm developed by Benítez et al. [48] was used, which is based on the Hilbert transform and has been evaluated in previous studies [27,49]. Then, HRV signals were derived by measuring the time between consecutive R peaks, i.e., RR intervals [19]. An artifact rejection step was then performed to ensure that all the intervals considered were physiologically plausible (N-N intervals). To this effect, we rejected those RR intervals that did not satisfy the following criteria [21]: (i) RR interval duration was between 0.33 and 1.5 s, and (ii) difference from the previous RR interval was higher than 0.66 s. After the artifact rejection, we confirmed that the duration of signals still surpassed 3 h of sleep. At this point, the remaining HRV samples were not equally distributed among the time; therefore, all the signals were resampled to a constant sampling frequency of 3.41 Hz [21,49], allowing us to perform higher-order spectra (HOS) analyses.

3. Methods

The methodology followed can be divided into three stages. First, we conducted a feature extraction stage to characterize each of the bispectral regions considered in the study. Second, we performed an automatic feature selection stage based on the fast correlation-based filter (FCBF) algorithm [50] considering two cases, one for the regions defined by classical HRV frequency ranges, and the other for the regions defined by OSA-related frequency ranges. Finally, a classification stage was conducted based on multi-layer perceptron (MLP) neural networks for the three AHI cutoffs used for partitioning OSA severity groups (1 e/h, 5 e/h and 10 e/h) and for each feature subset resulting from the feature selection stage.

3.1. Bispectrum Estimation

The estimation of power spectra has been one of the main tools for the analysis of biological signals for decades [29]. This technique contains information of the autocorrelation sequence, which is enough to characterize Gaussian signals, but information regarding the phase relationship among frequency components, as well as phase coupling, is lost during the process [29,34]. HRV signals, as with many other biological signals, are essentially nonlinear, non-Gaussian and non-stationary [29]. Therefore, power spectrum analysis may not be able to completely characterize changes in HRV series [34]. HOS analysis, meanwhile, contains both amplitude and phase information and can be used to characterize Gaussianity, stationarity and linearity deviations [29].
HOS are spectral representations of higher-order cumulants of a random process [29]. In particular, the bispectrum refers to the HOS for the third-order cumulant, reflecting spectral decomposition of the signal skewness over the frequency [29]. Bispectrum computation is based in the 2D Fourier transform of the third-order cumulant, and it can be defined in terms of the Fourier transform as [29,34]
B ( f 1 , f 2 ) = X ( f 1 ) · X ( f 2 ) · X * ( f 1 + f 2 ) ,             f 1 ,   f 2 = 0 , ,   f N  
where X(f) is the Fourier transform of a signal, f1 and f2 are the frequency indices, and fN is the Nyquist frequency. The resultant matrix reflects the phase coupling degree between frequency components for each pair f1,f2 [29].
As the bispectrum preserves both amplitude and phase information, it allows for assessment of interactions between signal patterns [29]. Likewise, bispectral analysis is used to evaluate changes in signal Gaussianity, where bispectral values = 0 indicate that signal components are Gaussian, and deviations from the Gaussianity of the components are otherwise detected [29]. Furthermore, bispectral analysis detects linearity deviations through phase coupling between its frequency components [29]. Phase coupling between three harmonic components at the f1, f2 and f3 frequencies and phase angles φ1, φ2 and φ3 is described as f3 = f1 + f2 and φ3 = φ1 + φ2 [29]. Thus, if phase coupling exists, it means that there are nonlinear interactions between harmonic components [29].
In this work, the HRV bispectrum was estimated employing a Hamming window of 210 samples with 50% overlapping and an FFT of 211 samples. After bispectrum matrix computation, a normalization was applied by dividing each element of the matrix by the sum of all matrix elements as [29,46]
B N ( f 1 , f 2 ) = B ( f 1 , f 2 ) B P ,             f 1 ,   f 2 = 0 , ,   f N
where BP is the total bispectral power, measured as the sum of all magnitudes of the bispectral matrix. This normalization was applied in order to ensure that all elements of the matrix were bounded between 0 and 1, reducing subject inter-variability due to physiological conditions other than OSA [46].

3.2. Determination of Bispectral Regions

The bispectral matrix, by definition, presents symmetric properties that render evaluation of a triangular non-redundant area sufficient for a full bispectrum characterization [29,34]. This region is commonly named the region of interest (ROI) and satisfies 0 ≤ f1f2f1 + f2fN [29,34]. As it was mentioned in the Introduction section, previous studies analyzing HRV bispectra focused their analysis along the whole ROI [30,34,35]. The analysis of HRV in the frequency domain, however, is commonly performed along the classic HRV spectral bands, i.e., the very low frequency (VLF, 0–0.04 Hz), low frequency (LF, 0.04–0.15 Hz) and high frequency (HF, 0.15–0.4 Hz) bands [18]. Past studies defined sub-band regions inside the bispectral ROI, bound by those frequencies [36,37,38], which we have termed classic bispectral regions. Furthermore, in a previous study applying a common spectral analysis, three pediatric OSA-related spectral ranges for HRV analysis were identified, which outperformed the classic spectral bands for pediatric OSA characterization and diagnosis [27]: BW1 (0.001–0.005 Hz), BW2: (0.028–0.074 Hz) and BWRes (0.04 Hz around HF peak). A detailed explanation of the process that led us to obtain those frequency ranges can be found in Appendix A. Following a similar reasoning to those studies that analyzed classic bispectral regions, three OSA-specific bispectral regions can be defined as bound by those OSA-related frequency ranges.
Therefore, in this study, six sub-band regions were assessed, three based on the classic HRV frequency ranges (VLF, LF and HF bispectral regions), and three based on the HRV OSA-related frequency ranges (BW1, BW2 and BWRes bispectral regions). To provide an overview, Figure 1 shows the averaged bispectrum magnitude in the training set for the four disease severity groups considered in the range 0–0.4 Hz. It can be observed that the no-OSA bispectral power is mainly concentrated below 0.02 Hz, and it spreads over a higher range of frequencies as the severity increases. The 3D bispectral region representations, averaged for each severity group, are shown in Appendix A Figure A1, Figure A2, Figure A3, Figure A4, Figure A5 and Figure A6.

3.3. Feature Extraction Stage

In order to characterize the bispectral regions under study, we computed features based on four different approaches: bispectral region amplitude, bispectral region entropy, bispectral region moment and weighted center of bispectrum (WCOB) region features. Furthermore, we introduced a new bispectral feature in this study.
As explained in Appendix A, BWRes is an adaptive frequency band based on the maximum value inside the HF range. It means that the frequency range is different for each subject, as it depends on the location of this peak. Therefore, contrary to the rest of the regions, BWRes might not be centered in the main diagonal of the bispectral matrix. Before feature extraction, we confirmed that this occurred in most of the subjects considered. As some of the features included in this study were computed over the diagonal of each region, the physiological meaning of the diagonal elements was lost (f1f2); therefore, these features were not computed over the BWRes region in this study.

3.3.1. Bispectral Region Amplitude Features

  • Maximum amplitude (Bmax), measured as the maximum magnitude value inside each of the regions considered [46]:
    B m a x = max ( | B N ( f 1 ,   f 2 ) | ) ,             f 1 ,   f 2     ,
    where represents one of the six regions considered.
  • Minimum amplitude (Bmin), measured as the minimum magnitude value inside each of the regions considered [46]:
    B m i n = min ( | B N ( f 1 ,   f 2 ) | ) ,             f 1 ,   f 2    
  • Total bispectral power (Btotal), measured as the sum of all magnitudes inside each of the regions considered [46]:
    B t o t a l = f 1 ,   f 2   | B N ( f 1 ,   f 2 ) | .
This parameter allows measuring deviations from Gaussianity [46].
Following a similar tendency to the spectral approach [27], as severity increases, higher values of bispectral amplitude features are expected in regions related to OSA, such as BW2 or LF. Consequently, lower values with OSA in regions related to respiration, such as HF or BWRes, are also expected.

3.3.2. Bispectral Entropy Features

  • Normalized bispectral entropy (BE1), normalized squared bispectral entropy (BE2) and normalized cubed bispectral entropy (BE3). These parameters, based on Shannon’s entropy, quantify the irregularity of the bispectral distribution in each region and are computed as [29,34]
    B E i = j   p j · log ( p j ) ,           i = 1 , 2 , 3
    where p is the magnitude distribution of the region elements:
    p j = | B N ( f 1 ,   f 2 ) | i f 1 ,   f 2   | B N ( f 1 ,   f 2 ) | i ,             i = 1 , 2 , 3
The values of the bispectral entropies increase with the randomness of a process, meaning changes in the HRV irregularity as a result of OSA [30] would be captured by the bispectral entropies of the regions.
  • Phase entropy (PE), which quantifies the phase regularity of the region [29]. PE, as with the bispectral entropies, is higher as the randomness of a process increases, meaning it would be zero for a harmonic, periodic and predictable process [34]. PE computation is performed applying Shannon’s entropy to the normalized distribution of the region phase angles [29,46]:
    P E = n   p ( Ψ n ) · log ( p ( Ψ n ) )
    where
    p ( Ψ n ) = 1 L f 1 ,   f 2   I n d ( φ [ B N ( f 1 ,   f 2 ) ]     Ψ n ) ,
    Ψ n = { φ | π + 2 π n N   φ < π + 2 π   ( n + 1 ) N } ,           n = 0 , 1 , , N 1
    where I n d ( · ) is the indicator function (equal to 1 if φ is within the range of histogram bins Ψ n ), φ is the phase angle of the region, and N is the bin number of the histogram.

3.3.3. Bispectral Region Moment Features

  • The sum of the logarithmic magnitude values of the region (H1), sum of the logarithmic magnitude values of the diagonal of the region (H2) and first- and second-order spectral moments of the magnitude values of the diagonal elements of the region (H3 and H4, respectively). These features were included as they allow characterizing the nonlinearity of the regions and are computed as follows [46]:
    H 1 = f 1 ,   f 2   log ( | B N ( f 1 ,   f 2 ) | ) .
    H 2 = f k ,     Γ d i a g log ( | B N ( f k ,   f k ) | ) .
    H 3 = f k ,     Γ d i a g k · log ( | B N ( f k ,   f k ) | ) .
    H 4 = f k ,     Γ d i a g ( k H 3 ) 2 · log ( | B N ( f k ,   f k ) | )
Those children suffering from OSA would be expected to present an increased bispectral power concentration in the region defined by frequency ranges related to the occurrence of apneic events (BW2). This would mean an increase in the phase coupling between the frequency components of this region and, accordingly, higher nonlinearity [30,46]. Therefore, OSA children would be expected to present higher values of bispectral region moment features in regions related to OSA, and lower values in the respiratory-related regions.

3.3.4. Bispectral WCOB Features

  • WCOB allows reflecting the interaction of different frequency components through the assignment of a weight to each bispectral point of the region [46]. The weighted center of each region is composed of two vectors, f1m and f2m, which indicate the coupling focus of the region as a summary of the frequency interaction [46]. Those components of WCOB are computed as [46]
    f 1 m = f 1 ,   f 2   f 1 · B N ( f 1 ,   f 2 ) f 1 ,   f 2   B N ( f 1 ,   f 2 ) ,
    f 2 m = f 1 ,   f 2   f 2 · B N ( f 1 ,   f 2 ) f 1 ,   f 2   B N ( f 1 ,   f 2 )
WCOB parameters are associated with bispectral peak values, with decreases in f1m and f2m values implying an activity shift towards lower frequencies [46]. Hence, as OSA children are expected to present higher activity in bispectral regions related to apneic events, their WCOB would be centered around these regions.

3.3.5. Relative Power of the Diagonal, a Novel Bispectral Feature

  • The relative power of the diagonal (RPDiag), computed as the sum of the bispectral amplitudes of the diagonal elements of the region, after a normalization applied over the whole diagonal. This novel parameter evaluates the relative bispectral magnitude value inside the diagonal of the region with respect to the complete bispectral diagonal magnitude:
    R P D i a g = f k     Γ d i a g | D i a g N ( f k ) | ,
    where Γ d i a g represents the diagonal elements of one of the regions considered except BWRes, and D i a g N is the normalized bispectral diagonal after the normalization performed such that
    D i a g N ( f k ) = D i a g ( f k ) D P ,           f k = 0 , ,   f N
    where Diag is the diagonal of the bispectral matrix, and DP is the diagonal power, measured as the sum of all amplitudes of the Diag elements.
The diagonal elements of a region are a particular case of the bispectral matrix when f1 = f2; therefore, this parameter, as well as H2, is intended to measure the phase coupling between the harmonic components of HRV signals, such that f3 = 2f1 and φ3 = 2φ1 [30,46].
RPDiag and H2 present two important differences. First, a normalization over the whole diagonal is applied in this novel feature. As a result of this normalization, the sum of all bispectral amplitudes of the diagonal elements is equal to 1, meaning RPDiag evaluates the proportion of the total diagonal bispectral power contained in the region. Then, as the normalization is scaling the values of the diagonal elements, we use a linear scale to compute the sum of the relative power instead of the logarithmic scale applied in H2. The rationale of this parameter lies in the normalization considering all of the frequency range. When OSA occurs, there is an alteration in the synchronization of the heart rhythm [30], leading to a redistribution of HRV activity to frequency components associated with the occurrence of apneic events. Thus, the normalization applied here considers not only the bispectral power in the diagonal of the region evaluated but also that in other diagonal elements. This influence of the redistribution to other frequency ranges is lost when applying a logarithmic scale.
Table 2 summarizes the bispectral features computed in every region under study.

3.4. Feature Selection Stage

Once each region was characterized with the features detailed in Table 2, we constructed two different optimal feature subsets (classic and OSA-specific bands) via the FCBF algorithm. This method, which has been previously demonstrated to be useful in pediatric OSA diagnosis [42,45,46,51], allows creating a non-redundant and relevant feature set based on the symmetrical uncertainty [50].
We performed the selection stage over 1000 bootstrap replicates from the training dataset in order to obtain generalizable and non-dependent subsets [52]. Those features selected more than 500 times were chosen to form the optimal subsets [42,45,51].

3.5. Classification Stage

Similar to previous pediatric OSA diagnosis studies, we conducted the classification stage using MLP neural networks [42,46,51]. These neural networks are typically formed by an input layer, a hidden layer and an output layer, each one composed of a different number of neurons, called perceptrons [53]. Each perceptron of an MLP network layer is connected to all the perceptrons from the next layer, with a weight associated with this connection [53]. The number of perceptrons in the first layer is equal to the number of input features. The number of perceptrons in the output layer depends on the objective of the network. We performed binary classification for the three severity thresholds (1, 5 and 10 e/h), as in our previous work [27]. This implies three different MLP neural networks for each feature subset, with one perceptron in the output layer providing the posterior probability of belonging to the severity group considered at each case [53]. The number of perceptrons in the hidden layer (NH) is a parameter to be optimized. To deal with overfitting, we also introduced a regularization parameter (λ) in the tuning of the network weights, which were randomly initialized [53].
The optimization of the hidden layer design parameters (NH and λ) was performed, again, by means of 1000 bootstrap replicates from the training dataset, but different from the replicates employed in the feature selection stage. We computed Cohen’s kappa (k) for each NH/λ combination and selected those values where k was maximum [42,46,51].
Thus, six MLP neural networks were optimized, one for each feature subset and severity threshold.

3.6. Statistical Analysis

In the training set, the features included in this study did not pass normality and homoscedasticity tests. For this reason, we assessed statistically significant differences (p-value < 0.01 after applying Bonferroni correction) between bispectral features from OSA severity groups through the non-parametric Kruskal–Wallis test. To provide a visual representation of these differences, along with the distribution followed by the features in each severity group, boxplots for each selected feature were also constructed.
Subsequently, we conducted a correlation analysis. To this effect, relationships between the features selected and some polysomnographic indices were evaluated using Spearman’s partial correlation coefficient (ρS), controlling the possible effect of age. The polysomnographic indices, related to OSA, as well as sleep structure and quality, were the same as those in [27]: AHI, Obstructive AHI (OAHI), obstructive apnea index (OAI), oxygen desaturation index (ODI), wake after sleep onset (WASO), number of awakenings during total sleep time (#Awakenings), percentage of total sleep spent in sleep stages N1, N2, N3 and rapid eye movement (%N1, %N2, %N3 and %REM, respectively) and total arousal index per hour of sleep (TAI). Correlation analysis was performed on the test set.
Finally, after the optimization of the MLP neural networks in the training set, the diagnostic performances of each individual selected feature and optimized model were evaluated in the test set in terms of sensitivity (Se), specificity (Sp), accuracy (Acc) and area under the receiver operating characteristic curve (AUC).

4. Results

4.1. Feature Selection in the Training Set

We conducted two feature selection processes. As Table 2 shows, the number of input features for the classic bispectral region set was 42, while in the case of the OSA-specific region set, there were up to 38 features. The feature selection through FCBF allowed reducing the amount of redundant information while assessing feature complementarity. Figure 2 shows the number of times that each feature was selected by the algorithm in the 1000 bootstrap replicates for both cases. In the case of the classic regions (Figure 2a), it can be observed that the optimum subset (BISPClassic) was formed by three features, one of each region: VLF_f2m, LF_BE2 and HF_PE. Regarding the OSA-specific region set (Figure 2b), the optimum subset (BISPSpecific) was composed of four features selected more than 500 times: BW2_RPDiag, BW2_BE1, BWRes_Bmin and BWRes_BE3. None of the BW1 region features considered were selected over 500 times.

4.2. Descriptive Analysis of the Features Selected

Figure 3 and Figure 4 show the boxplot distributions of the four OSA severity groups for the features selected in both the BISPclassic and BISPSpecific subsets, respectively. The p-value resulting from the Kruskal–Wallis test is also depicted in these figures. It can be appreciated that, in the BISPClassic subset, VLF_f2m experienced an increase with OSA severity, while a decrease in LF_BE2 and HF_PE occurred as the OSA severity increases. For the features included in BISPSpecific, there was a clear rise in the BW2_RPDiag values, along with a slight increase in BWRes_BE3 with OSA severity. In contrast, the BW2_BE1 values experienced a decrease with the disease. BWRes_Bmin was the only parameter showing an unclear tendency among the severity groups, which led it to be the only one that did not show statistically significant differences. The remaining six parameters showed statistically significant differences among the four OSA severity groups (p-value < 0.01 after Bonferroni correction).

4.3. MLP Network Optimization and Training

After extraction of the optimum feature subsets, six MLP models were optimized: three models with BISPclassic features as input (MLP1Classic, MLP5Classic and MLP10Classic, with AHI = 1, 5 and 10 e/h as thresholds for binary classification, respectively), and three models with BISPSpecific as input features (MLP1Specific, MLP5Specific and MLP10Specific, with AHI = 1, 5 and 10 e/h as a threshold for binary classification, respectively). For each model, NH varied from 2 to 20 in steps of 1, and from 22 to 50 in steps of 2. Similarly, λ varied from 0.5 to 10 in steps of 0.5. Each NH/λ pair resulted in an averaged k through 1000 bootstrap replicates of the training set; therefore, we selected the NH/λ combination with the higher averaged k. NH = 2 and λ = 5 were the optimized design parameters selected in four out of six MLP models: MLP1Classic, MLP5Classic, MLP5Specific and MLP10Specific. For the MLP10Classic model, the optimized design parameters were NH = 34 and λ = 5. Finally, the optimized parameters in MLP1Specific were NH = 38 and λ = 5.

4.4. Correlation Analysis in the Test Set

The results of the correlation study are shown in Table 3. Although the values of |ρS| obtained were generally low, some of the correlations evaluated were statistically significant (p-value < 0.01 after Bonferroni correction). VLF_f2m and BW2_RPDiag showed similar behaviors, with a statistically significant positive ρS with the four respiratory indices (AHI, OAHI, OAI and ODI) and TAI, as well as a negative ρS with %REM. In the opposite way, BW2_BE1 obtained a negative ρS with the four respiratory indices and TAI. LF_BE2 also reached a statistically significant negative ρS with AHI, OAHI and TAI. BW2_RPDiag reached the highest absolute correlation values among almost all of these statistically significant correlations, only being equaled by the ρS reached between VLF_f2m and OAHI. None of the selected BWRes features nor HF_PE obtained statistically significant correlations with any of the polysomnographic indices considered in this study.

4.5. Diagnostic Ability Assessments

Table 4 shows the diagnostic performance obtained in the test set by each feature individually, along with the classification results reached by each MLP model optimized for the three severity thresholds.
Regarding the individual performance, in the 1 e/h threshold, BW2_RPDiag obtained the highest results in terms of Acc and AUC. For the 5 e/h threshold, again, BW2_RPDiag showed a higher Acc than the other features selected, being only slightly surpassed by VLF_f2m in terms of AUC. When considering 10 e/h as severity cutoff for binary classification, LF_BE2 was the feature showing the higher Acc, but with a lower AUC and a more unbalanced Se/Sp pair than BW2_RPDiag.
All the individual diagnostic yields were outperformed by the MLP models. The MLP models formed by the OSA-specific region features obtained the highest results of this study in terms of Acc and AUC in the 1 and 5 e/h thresholds. The MLP10Classic model was the only one that achieved a higher diagnostic performance than the OSA-specific models in terms of Acc and AUC at the cost of a strongly unbalanced Se/Sp pair and a very low Se value (43.5%).

5. Discussion

In this work, bispectral analysis of HRV signals, a process performed for the first time in pediatric OSA research, was conducted herein. Features from bispectral HRV regions based on the classic and OSA-specific frequency ranges were extracted to assess their usefulness in the characterization and diagnosis of pediatric OSA. In both types of regions separately, the selected features showed their complementarity, and the models constructed achieved a high diagnostic performance. The OSA-specific region models generally outperformed the classic ones, highlighting the importance of these novel bispectral features in the study of pediatric OSA.

5.1. Physiological Interpretation of the Features Selected

The averaged normalized bispectrum in the training set shown in Figure 2 serves as a summary of the degree of phase coupling in the frequency range 0–0.4 Hz, covering all the regions considered. The bispectral power is mainly concentrated under 0.02 Hz in the no-OSA subjects, and then it spreads with severity to higher frequencies. This coupling focus is markedly lower in the severe OSA group, as it can be observed in the amplification depicted in each upper right corner of the representations from Figure 2. It reflects an increase in the linearity and Gaussianity of HRV signals at very low frequencies due to apneic events [29]. This shift to higher frequencies is also reflected in the selected feature VLF_f2m. The values of VLF_f2m increase with OSA severity, as depicted in Figure 3a. This means that the focus of coupling in the VLF range is displaced by apneic events, generating higher HRV activity at higher frequencies [54]. This feature makes sense, as the higher frequencies of the VLF range (0–0.04 Hz) overlap with the lower frequencies of BW2 (0.028–0.074 Hz), which has been defined as the frequency range related to the duration of apneic events [27].
Inside this BW2 apneic-related region, BW2_RPDiag was one of the two features selected. This parameter, as with VLF_f2m, showed an increasing tendency with OSA severity (Figure 4a). As explained in the Methods section, RPDiag is intended to measure phase coupling between the harmonic components (f3 = 2f1 and φ3 = 2φ1) of the HRV. The increment in BW2_RPDiag with OSA severity would indicate increasing nonlinear interactions between those harmonics of OSA-affected children. Consequently, there is an increment in less random/more periodic harmonics in HRV signals inside the BW2 region due to apneic events. The increase in VLF_f2m and BW2_RPDiag with OSA severity is supported by the correlation study, showing statistically significant correlations with all respiratory indices, as well as with TAI. Furthermore, BW2_RPDiag generally obtained the highest individual diagnostic performance. Taken together, these facts highlight the importance of analyzing the BW2 bispectral region in the field of pediatric OSA, especially characterized through our new proposed parameter RPDiag.
The entropy features also showed their usefulness to characterize the bispectral regions of the HRV. Three entropies of the bispectral amplitude distribution were selected: BE1 inside the BW2 region, BE2 inside the LF region and BE3 inside BWRes. These parameters measure the irregularity of the HRV from the bispectral distribution in each region, with the inclusion of quadratic and cubic components scaling the differences in the bispectral amplitude [29]. Bispectral distributions averaged for each severity group in these three regions are shown in Figure A2 (LF region), A4 (BW2 region) and A5 (BWRes region). In the case of BW2_BE1, a decrease with severity was observed, reflecting a reduction in irregularity with apneic events in the HRV components linked to the frequencies of this region. As a result of apneic events, it can be appreciated that the bispectral amplitude in BW2 is more concentrated at low frequencies in the severe group (Figure A4d) and starts to distribute more randomly to other frequencies as OSA decreases. This may be due to the aforementioned increment in the less random harmonics in this range due to OSA, leading to the reduction with severity experimented in BW2_BE1 (Figure 4b). In the case of LF_BE2, this parameter also decreases with the severity of the disease, again reflecting a reduction in irregularity with OSA in the HRV associated with this frequency region. Figure A2 points to the fact that the bispectral power distribution of no-OSA subjects is more dispersed over the whole LF region and starts to concentrate at lower frequencies as OSA severity increases. Interestingly, this parameter only showed negative statistically significant correlation values when considering respiratory indices that include hypopneas (AHI and OAHI). It seems that, as apneas are less frequent than hypopneas in the database under study, the only effect of apneas is not enough to decrease the bispectral HRV irregularity associated with the LF frequency range. Similarly, collective apneic effects are better captured with the inclusion of the quadratic amplitude, suggesting that BE2 is more accurate in detecting alterations in the bispectral distribution of LF as a result of all apneic events. Regarding BWRes_BE3, Figure 4d shows that no-OSA, mild OSA and moderate OSA subjects presented lower values than the severe OSA group. As it can be seen in Figure A6, there is a higher bispectral power concentration around the respiratory peak for the first three severity groups, and a lower concentration for severe OSA, whose distribution spreads over other frequencies. The lower coupling around the respiratory peak in the severe group makes sense as OSA results in a redistribution in the bispectral power to frequency ranges related to apneic events, such as the BW2 region. These are milder differences than those observed in BW2 and LF; therefore, an increase in HRV irregularity due to OSA appears to be better captured through BE3 from BWRes.
In addition to the bispectral entropies, PE was also selected in the HF range. This parameter showed a decreasing tendency with OSA severity (Figure 3c). As a reduction in PE indicates that a process becomes less random [34], this result suggests that OSA alterations lead to a reduction in the irregular behavior of the HRV phase along the HF region. The fact that every entropy measure included in this study was selected, at least, in one region highlights the importance of the entropy features when analyzing the bispectral HRV distribution in pediatric OSA.
BWRes_Bmin was the remaining feature selected. The normalization applied over each bispectral matrix allowed Bmin to estimate the minimum coupling within this region [46]. Despite the unclear tendency and the absence of differences obtained in this parameter (Figure 4c), BWRes_Bmin was selected by the algorithm. This implies that BWRes_Bmin contains information that is complementary to the other features selected in the OSA-specific region feature subset.
Although BW1 showed its potential utility in the spectral analysis [27], none of the features included in this region were selected by the FCBF algorithm a number of times above the threshold established. This suggests that, when analyzing pediatric OSA alterations in the HRV bispectrum, the assessment of the BW2 and BWRes regions would be enough to characterize OSA effects. Similarly, the features related to bispectral moments were not selected either, probably due to the redundancy when introducing RPDiag, which seems to be more accurate in the characterization of apneic alterations.

5.2. Diagnostic Performance of the Bispectral Models

In this study, the information extracted from the bispectral HRV regions obtained an overall high diagnostic performance both in the classic and OSA-specific bispectral region models in the test set. The results obtained in the MLP models outperform the individual diagnostic yield, highlighting the utility of the FCBF algorithm and MLP neural networks to assess complementarity between features and to diagnose pediatric OSA, respectively. When comparing classic against OSA-specific region models, the latter generally obtained a higher performance. In the 1 e/h severity cutoff, despite the unbalanced Se/Sp pair, MLP1Specific obtained a higher Acc (63.4% versus 54.7% from the MLP1Classic model) and a higher AUC (0.627 versus 0.600). In the 5e/h threshold, the Acc obtained by MLP5Specific and MLP5Classic was 81.0% in both cases, but with the specific model, showing a more balanced Se/Sp, a higher AUC was found (0.774 versus 0.791). Finally, in the 10e/h severity cutoff, MLP10Classic obtained the highest Acc and AUC from the study. However, the Se/Sp pair was strongly unbalanced, with a very low Se value of 43.5%. Nevertheless, the MLP10Specific model, at the cost of a very slight reduction in terms of Acc and AUC (89.3% versus 91.7%, and 0.847 versus 0.841, respectively), resulted in a more balanced Se/Sp pair (66.7%/91.6%). These results reinforce the conclusion from our previous work [27] about the importance of analyzing OSA-specific HRV frequency ranges whenever pediatric OSA is under study.

5.3. Comparison with Previous Work

As far as we know, this is the first study where HRV bispectral analysis of pediatric OSA was conducted. In this sense, a direct comparison of classification results with other research studies using an HRV bispectral approach is not possible. However, HRV analysis to diagnose pediatric OSA has been previously performed. A research group carried out three classification studies [26,55,56], where they derived HRV features from decreases in the amplitude fluctuations in the photoplethysmography signal. In their studies, they included a population of 21 children (10 children with OSA and 11 controls) and obtained an Acc ranging from 73.3% to 80%, Se from 62.5% to 87.5% and Sp from 71.45% to 85.7%, when considering HRV features only. The highest results obtained in the present work outperform these diagnostic results in terms of Acc and Sp, with Se in the same range. Nonetheless, the difference in the child population, along with the different criteria followed to establish OSA, makes it difficult to perform a more comprehensive comparison. The study developed by Cohen and de Chazal [57] was the first study conducting an automated classification of children with OSA using only the HRV signal. Notwithstanding, this classification was based on detecting apneic events, rather than global classification of each subject, meaning a comparison with the diagnostic results of the present study is not possible [57].
Thus, our previous study [27] is the only one that establishes a benchmark with which we can compare the diagnostic performance reached in the present study. Furthermore, the comparison of these results is direct as the database, the distribution of training/test sets and the criterion of OSA diagnosis were the same. In the previous work, we computed the individual diagnostic performances of each relative power (RP) from the power spectral density in the VLF, LF, HF, BW1, BW2 and BWRes bands, and the LF/HF ratio too. Then, we constructed two linear discriminant analysis (LDA) models to assess the joint diagnostic yield of the relative powers in the classic bands versus the OSA-specific bands. Following this methodology, the best classification outcomes obtained in terms of Acc and AUC using 1 e/h as the threshold for binary classification were 59.2% Acc (RPBW1 individually) and 0.592 AUC (LDA band of interest model). With 5 e/h as the severity threshold, the highest results were 76.6% Acc (RPBW2) and 0.688 AUC (LDA band of interest model). Finally, in the 10 e/h severity cutoff, the highest results reached were 82.8% Acc and 0.796 AUC (both using the LDA band of interest model). It can be observed from Table 4 that, in the present work, the best results achieved with the MLP optimized models surpassed, by far, those achieved previously, and even some of the individual bispectral features eventually outperformed several of these results. It can be argued that the improvement in the diagnostic outcomes may be mediated by the increased complexity of the classification algorithm. To deal with this issue, and also for a fair comparison, we have included in Appendix B the results obtained using the same classification methodology as in [27]. Table A1 shows the classification results obtained including the features selected from each approach using an LDA classifier. The results obtained with both the LDAClassic and LDASpecific models also outperform the diagnostic performance obtained in the previous work, being only surpassed in terms of Acc in the 1 e/h threshold. Thus, the diagnostic utility of the features extracted from the bispectral HRV region analysis is clearly demonstrated, with the MLP models reaching the highest diagnostic performance in the literature when using HRV features exclusively to generate an automated classification of pediatric subjects into the presence or absence of OSA and to estimate the severity grouping. Moreover, the higher diagnostic performance reached by the bispectral analysis highlights the usefulness of this analysis in the pediatric OSA context, which seems to be more accurate than traditional frequency analysis to evaluate HRV alterations. The presence of HRV nonlinear dynamics demonstrated in Appendix C and captured through the HRV bispectrum estimation may be behind this improvement over the traditional techniques.
Finally, although the results of the diagnostic performance from studies using physiological data from different sources must be carefully compared, it is interesting to contrast our results with the meta-analysis of a recently published systematic review [14]. In this work, Gutiérrez-Tobal et al. gathered the pooled Se and Sp results from nineteen studies of machine learning methods to diagnose pediatric OSA that fulfilled their eligibility criteria [14]. A meta-analysis was performed for the same OSA severity thresholds that were employed here, which obtained an Se of 84.9%, 71.4% and 65.2% for the 1, 5 and 10 e/h cutoffs, respectively. Similarly, the meta-analysis obtained an Sp of 49.9%, 83.2% and 93.1% for the 1, 5 and 10 e/h cutoffs, respectively. Table 4 shows that the MLP models presented here for both classic and specific approaches obtained a diagnostic performance in the same range for Sp in the 5 and 10 e/h cutoffs, and also MLP10Specific obtained similar Se results. It is worth noting that, among the studies included in the meta-analysis, none of them considered HRV signals. Furthermore, the sample size from seventeen out of nineteen studies included in the systematic review was markedly smaller than the databases analyzed here, only being surpassed by the cohort included in the studies of Hornero et al. [42] and Vaquerizo-Villar et al. [58]. Therefore, these considerations along with the similar diagnostic performance reached in the 5 and 10 e/h severity thresholds reinforce the support for the use of HRV bispectral analysis as a potential alternative to overnight PSG for pediatric OSA diagnosis.

5.4. Limitations and Outlook

Despite the high diagnostic performance obtained with the methodology followed in this study, some limitations deserve mention. First, although the sample size included is markedly large, accounting for 1738 overnight HRV recordings, some imbalance between severity groups is apparent. In this sense, an increase in the population included in future studies, trying to balance the severity groups’ distribution, would be desirable to raise the robustness and generalizability of our results.
Additionally, none of the features extracted from the region bound by BW1 frequencies (0.001–0.005 Hz) were selected. This frequency range has been linked to sleep fragmentation due to OSA, and its usefulness in spectral analysis has been demonstrated [27]. Thus, a combination of features from this region with features derived from different approaches would be performed in future studies to assess if BW1 OSA alterations are really reflected in bispectral analysis, and also if they contain non-redundant information on the features presented here.
A previous study performing bispectral HRV analysis in the adult OSA context only analyzed the non-redundant bispectral region [30]. In this sense, despite the well-known differences between children and adults in the OSA context, the high diagnostic performance achieved here in analyzing classic and OSA-specific bispectral regions serves as a motivation to search for adult OSA-specific frequency ranges. Consequently, a replication of the methodology carried out here to extract OSA-specific regions in adults and its analysis in the frequency and bispectral domains is one of our future research aims.
Additionally, the highly satisfactory results achieved illustrate the utility of bispectral HRV analysis to characterize and diagnose pediatric OSA, outperforming the previous spectral analysis diagnosis yield [27]. Therefore, the binary classification performed serves as a first step, and further explorations of more complex predictive models, as well as estimation of the AHI (regression), instead of binary classification, are also some of our future research aims.
Finally, despite the high diagnostic performance achieved here through bispectral analysis and the MLP models constructed, the promising results obtained by deep learning techniques in healthcare issues in recent years highlight the potential utility of these methods to automate the diagnosis of pediatric OSA [14]. Accordingly, in trying to increase the diagnostic performance, the inclusion of deep learning methods in the pediatric OSA context is a future research need.

6. Conclusions

To the best of our awareness, this is the first work where bispectral HRV analysis has been conducted to characterize and diagnose pediatric OSA. Our methodology allowed us to obtain two feature subsets, one containing information regarding bispectral regions based on classic HRV frequency ranges, and the other one with OSA-specific bispectral regions. Those subsets were formed by features containing complementary information about alterations in the non-Gaussianity, nonlinearity and irregularity behavior of the HRV due to OSA. Among the features selected, a novel bispectral measure presented here, RP_Diag, showed its utility, generally achieving the highest individual diagnostic performance as well as the highest correlations with polysomnographic indices. Furthermore, the MLP models outperformed the previous results of diagnostic performance based on spectral analysis, with the MLP1Specific, MLP5Specific and MLP10Classic models achieving the highest diagnostic yield from the study for each severity cutoff. These results highlight the usefulness of bispectral HRV analysis in the pediatric OSA context, especially when analyzing bispectral regions bounded by OSA-specific frequency ranges. Thus, we conclude that information extracted from HRV bispectra allows for the characterization and diagnosis of pediatric OSA, leading us to propose this approach as a potential alternative to PSG.

Author Contributions

Conceptualization, A.M.-M., G.C.G.-T. and R.H.; methodology, A.M.-M., G.C.G.-T. and D.Á.; software, A.M.-M., G.C.G.-T. and V.B.-G.; validation, A.M.-M., G.C.G.-T. and D.Á.; formal analysis, A.M.-M., G.C.G.-T. and D.Á.; investigation, A.M.-M.; resources, L.K.-G., F.d.C. and R.H.; data curation, A.M.-M. and V.B.-G.; writing—original draft preparation, A.M.-M. and G.C.G.-T.; writing—review and editing, A.M.-M., G.C.G.-T., L.K.-G., V.B.-G., D.Á., F.d.C., D.G. and R.H.; supervision, R.H.; project administration, F.d.C., D.G. and R.H.; funding acquisition, D.G. and R.H. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by ‘Ministerio de Ciencia, Innovación y Universidades-Agencia Estatal de Investigación’ and ‘European Regional Development Fund (FEDER)’ under projects DPI2017-84280-R and RTC-2017-6516-1, by ‘European Commission’ and ‘FEDER’ under project ‘Análisis y correlación entre la epigenética y la actividad cerebral para evaluar el riesgo de migraña crónica y episódica en mujeres’ (‘Cooperation Programme Interreg V-A Spain-Portugal POCTEP 2014–2020′) and by ‘CIBER en Bioingeniería, Biomateriales y Nanomedicina (CIBER-BBN)’ through ‘Instituto de Salud Carlos III’ co-funded with FEDER funds, as well as under the project SleepyHeart from 2020 valorization call. A.M.-M. is in receipt of an ‘Ayudas para contratos predoctorales para la Formación de Doctores’ grant from the Ministerio de Ciencia, Innovación y Universidades (PRE2018-085219). V.B.-G. is in receipt of an ‘Ayuda para financiar la contratación predoctoral de personal investigador’ grant from ‘Consejería de Educación de la Junta de Castilla y León’ cofunded by ‘European Social Fund’. D.A. is supported by a ‘Ramón y Cajal’ grant (RYC2019-028566-I) by the ‘Ministerio de Ciencia e Innovación–Agencia Estatal de Investigación’ co-funded by ESF. L.K.-G. and D.G. are supported by National Institutes of Health (NIH) grant HL130984, the Leda J. Sears Foundation and by a Tier 2 grant from the University of Missouri. D.G. is also supported by NIH grants HL140548, and AG061824.

Institutional Review Board Statement

This study was performed following the guidelines of the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of the UofC Medicine Comer Children’s Hospital (#11-0268-AM017, #09-115-B-AM031 and #IRB14–1241).

Informed Consent Statement

Informed consent from caretakers of all children included in the study was acquired.

Data Availability Statement

The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the privacy of individuals that participated in the study.

Acknowledgments

The Childhood Adenotonsillectomy Trial (CHAT) was supported by the National Institutes of Health (HL083075, HL083129, UL1-RR-024134, UL1 RR024989). The National Sleep Research Resource was supported by the National Heart, Lung, and Blood Institute (R24 HL114473, 75N92019R002).

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Determination of HRV OSA-Specific Frequency Ranges and Averaged Bispectral Regions in the Training Set

The classic HRV spectral frequency ranges present some limitations when analyzing alterations in the autonomic nervous system as a result of pediatric OSA. This is the main motivation that led to us to the search of pediatric HRV OSA-specific frequency ranges in a recent work [27]. To perform the extraction of those spectral bands, we first established four severity groups based on the AHI. The HF frequency range, also known as the respiratory frequency band, is strongly modulated by age. For this reason, the search was separated in two different analyses to extract bands of interest based on the power spectral density. First, we performed an analysis over the 0–0.15 Hz frequency range, as it should not be altered by age [27]. The seek was performed through the evaluation of statistical differences (Mann–Whitney U test), frequency by frequency, between the amplitude of power spectral densities along this range for each severity group. Those frequency bands where, at least, two of the statistical tests showed statistically significant differences (p-value < 0.01 after Bonferroni correction) were selected as OSA-specific frequency ranges. This methodology allowed us to find two bands of interest, termed BW1 (0.001–0.005 Hz) and BW2 (0.028–0.074 Hz).
The second analysis performed was oriented to consider age modulation in the respiratory range. In this sense, age plays an important role in the location of the respiratory peak within HF (0.15–0.40 Hz) [59,60], meaning the assessment was an adaptive individual analysis. We first located the individual maximum amplitude inside the HF range. Therefore, as it was explained in the Discussion section in [27], a frequency range of 0.04 Hz around this respiratory peak was enough to assess pediatric OSA alterations in this adaptive individual frequency band, which we termed BWRes. A complete and detailed explanation of the whole methodology to extract those HRV OSA-specific frequency bands of interest can be found in the Methods section from our previous study [27].
Thus, now that the extraction of HRV OSA-specific frequency ranges has been described, below, we include the averaged bispectral distribution for each severity group, bounded by both classic and OSA-specific frequencies. Figure A1, Figure A2, Figure A3, Figure A4, Figure A5 and Figure A6 show the 3D representation for the VLF, LF, HF, BW1, BW2 and BWRes regions, respectively. The x and y axes from the BWRes region cannot be represented in terms of frequency, as it changes for each subject. That is why they appear in terms of the number of samples for this region.
Figure A1. Averaged bispectrum magnitude in the very low frequency region (0–0.04 Hz) for the four severity groups in the training set. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Figure A1. Averaged bispectrum magnitude in the very low frequency region (0–0.04 Hz) for the four severity groups in the training set. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Entropy 23 01016 g0a1
Figure A2. Averaged bispectrum magnitude in the low-frequency region (0.04–0.15 Hz) for the four severity groups in the training set. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Figure A2. Averaged bispectrum magnitude in the low-frequency region (0.04–0.15 Hz) for the four severity groups in the training set. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Entropy 23 01016 g0a2
Figure A3. Averaged bispectrum magnitude in the high-frequency region (0.15–0.40 Hz) for the four severity groups in the training set. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Figure A3. Averaged bispectrum magnitude in the high-frequency region (0.15–0.40 Hz) for the four severity groups in the training set. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Entropy 23 01016 g0a3
Figure A4. Averaged bispectrum magnitude in the BW1 frequency region (0.001–0.005 Hz) for the four severity groups in the training set. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Figure A4. Averaged bispectrum magnitude in the BW1 frequency region (0.001–0.005 Hz) for the four severity groups in the training set. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Entropy 23 01016 g0a4
Figure A5. Averaged bispectrum magnitude in the BW2 frequency region (0.028–0.074 Hz) for the four severity groups in the training set. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Figure A5. Averaged bispectrum magnitude in the BW2 frequency region (0.028–0.074 Hz) for the four severity groups in the training set. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Entropy 23 01016 g0a5
Figure A6. Averaged bispectrum magnitude in the BWRes frequency region (0.04 Hz around the respiratory peak inside HF) for the four severity groups in the training set. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Figure A6. Averaged bispectrum magnitude in the BWRes frequency region (0.04 Hz around the respiratory peak inside HF) for the four severity groups in the training set. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Entropy 23 01016 g0a6

Appendix B. Diagnostic Performance of Bispectral Region Models with a Linear Discriminant Analysis Classifier

To perform the main automated classification analysis of this study, we optimized the MLP models for each severity threshold (results available in Table 4). However, for comparison purposes, we have included in this appendix the diagnostic performance achieved by each optimized feature subset following the same methodology as in our previous work [27]. To this effect, the diagnostic performance of the BISPClassic and BISPSpecific feature subsets was evaluated through two models based on LDA (LDAClassic and LDASpecific models, respectively). Both classifiers were trained in the training set for each severity threshold, and the diagnostic performance was evaluated in the test set. Table A1 shows the results achieved by each model, along with the results from our previous work [27]. It can be observed that, as it is commented on in the Discussion section, the new diagnostic performance surpassed our previous classification results in terms of Acc and AUC in almost all the parameters. Only the individual Acc in 1 e/h for RPBW1 outperformed the LDA results presented here, but with a lower AUC.
Table A1. Diagnostic performance in the test set achieved following a spectral and a bispectral approach. The individual performances correspond to each relative power extracted in our previous work. Joint performance was assessed by constructing LDA classifiers.
Table A1. Diagnostic performance in the test set achieved following a spectral and a bispectral approach. The individual performances correspond to each relative power extracted in our previous work. Joint performance was assessed by constructing LDA classifiers.
Feature/ModelAHI Threshold = 1 e/hAHI Threshold = 5 e/hAHI Threshold = 10 e/h
SeSpAccAUCSeSpAccAUCSeSpAccAUC
RPVLF68.931.656.30.51833.065.060.20.45640.664.262.10.495
Previous work
(frequency analysis approach)
RPLF43.562.950.10.55752.758.457.60.59059.458.458.50.666
RPHF35.571.947.80.52339.368.163.80.54043.576.773.70.605
LF/HF37.770.348.70.54045.566.863.70.56749.370.868.80.643
RPBW166.345.359.20.55965.254.055.60.62169.652.353.90.624
RPBW232.778.148.10.59145.582.076.60.67058.078.276.40.751
RPBWRes45.556.649.30.53244.664.061.20.57149.364.062.60.628
LDA Classic Bands25.781.344.50.55946.472.268.40.63350.775.373.10.685
LDA Bands of Interest42.572.352.60.59250.080.976.40.68863.884.782.80.796
Present work
(bispectral analysis approach)
LDAClassic30.181.347.40.60153.685.380.60.77966.789.787.60.847
LDASpecific37.977.351.30.61563.482.879.90.79271.085.984.50.842
RP: relative power; VLF: very low frequency; LF: low frequency; HF: high frequency; LDA: linear discriminant analysis; AHI: Apnea–Hypopnea Index. The highest accuracy and AUC in each severity threshold are highlighted in bold.

Appendix C. Surrogate Data Approach

Appendix C.1. Testing for Nonlinearities

Despite the previous literature and evidence supporting the utility of nonlinear HRV analysis in the pediatric OSA context commented on in the Introduction section, testing the presence of nonlinearities in our database can increase the rationale in the selection of bispectral analysis over other traditional methods. To this effect, a surrogate data approach can be performed, which allows for these types of test [61]. The surrogate data generation techniques to test nonlinearities construct surrogate data from the original signal, preserving the power spectrum and data distribution while removing nonlinear properties [61,62]. After seeking for different surrogate data generation techniques, we decided to implement the iterative amplitude-adjusted Fourier transform method (IAAFT) [61,63], which is one of the most popular solutions [61].
Testing for nonlinearities with the IAAFT method first implies constructing a set of surrogate data. We performed this test over four randomly selected subjects from the training set, each one belonging to an OSA severity group (no-OSA: AHI < 1; mild OSA: 1 ≤ AHI < 5; moderate OSA: 5 ≤ AHI < 10; severe OSA: AHI ≥ 10). Thus, we generated 100 realizations of surrogate data from the original selected HRVs. The null hypothesis established was that the data represent a linear process, with all surrogate data being consistent with this hypothesis [61,63]. The next step was to compute a nonlinear measure for both original and all surrogate data and rank them in numerical order. We selected the Lempel-Ziv complexity (LZC) as a nonlinear measure that evaluates the generation rate of patterns present in a data sequence [64]. If the null hypothesis is true, the nonlinear measure evaluated in the original data should be indistinguishable from the measure computed for the surrogates [63]. However, for a two-sided test, if the nonlinear parameter is lower or higher than the rest of the measures computed in the surrogates, the null hypothesis can be rejected with a p-value < 0.02, and the existence of statistically significant nonlinearities in the original signal can be ensured. We confirmed that it occurs in the four HRV signals included, which means that the presence of nonlinearities in HRV signals under our study conditions has been demonstrated.

Appendix C.2. Bispectrum with Surrogate

The surrogate data method can also be employed to check for the significance of bispectral peaks. Although peaks derived from frequency coupling alone should not appear in the bispectrum computation, sometimes they appear mixed with true peaks resulting from the frequency and phase coupling [62]. To solve this problem, the bispectrum with surrogate (BWS) computation method can be utilized [62]. Computing a bispectrum with and without the BWS technique and comparing these results allow assessing if the main peaks present in our study can be considered as significant.
The BWS technique is based on surrogate data. Thus, following the approach established in the subsection above, we performed BWS computation over the same four HRVs from the different OSA severity groups. Once the 100 realizations of surrogate data for each HRV were computed, we estimated the bispectrum for the original HRV and also for each surrogate. Then, over the surrogate set, we computed the mean and standard deviation of the bispectral magnitude, obtaining 100 means and standard deviations [62]. Based on these statistics, we constructed a 95% statistical threshold as the mean plus two times the standard deviation for each surrogate and selected the maximum threshold value as the threshold to apply over the original bispectrum estimation [62]. Any bispectral peak present in the original bispectrum estimation above the threshold computed was considered as significant, meaning BWS estimations were constructed canceling the bispectral magnitude values below the threshold. The next figures (Figure A7, Figure A8, Figure A9 and Figure A10) show the original bispectrum estimation along with the BWS estimation represented in 2D and 3D for each severity group.
Figure A7. Bispectrum magnitude estimation with and without surrogate method for a no-OSA subject. (a) 2D bispectrum estimation; (b) 2D BWS estimation; (c) 3D bispectrum estimation; (d) 3D BWS estimation.
Figure A7. Bispectrum magnitude estimation with and without surrogate method for a no-OSA subject. (a) 2D bispectrum estimation; (b) 2D BWS estimation; (c) 3D bispectrum estimation; (d) 3D BWS estimation.
Entropy 23 01016 g0a7
Figure A8. Bispectrum magnitude estimation with and without surrogate method for a mild OSA subject. (a) 2D bispectrum estimation; (b) 2D BWS estimation; (c) 3D bispectrum estimation; (d) 3D BWS estimation.
Figure A8. Bispectrum magnitude estimation with and without surrogate method for a mild OSA subject. (a) 2D bispectrum estimation; (b) 2D BWS estimation; (c) 3D bispectrum estimation; (d) 3D BWS estimation.
Entropy 23 01016 g0a8
Figure A9. Bispectrum magnitude estimation with and without surrogate method for a moderate OSA subject. (a) 2D bispectrum estimation; (b) 2D BWS estimation; (c) 3D bispectrum estimation; (d) 3D BWS estimation.
Figure A9. Bispectrum magnitude estimation with and without surrogate method for a moderate OSA subject. (a) 2D bispectrum estimation; (b) 2D BWS estimation; (c) 3D bispectrum estimation; (d) 3D BWS estimation.
Entropy 23 01016 g0a9
Figure A10. Bispectrum magnitude estimation with and without surrogate method for a severe OSA subject. (a) 2D bispectrum estimation; (b) 2D BWS estimation; (c) 3D bispectrum estimation; (d) 3D BWS estimation.
Figure A10. Bispectrum magnitude estimation with and without surrogate method for a severe OSA subject. (a) 2D bispectrum estimation; (b) 2D BWS estimation; (c) 3D bispectrum estimation; (d) 3D BWS estimation.
Entropy 23 01016 g0a10
The first outcome that can be appreciated is that the original bispectrum estimations (Subfigures a and c from each figure) present the main bispectral peak as well as some smaller peaks to a greater or lesser degree. After applying the BWS method (Subfigures b and d), most of these spurious peaks disappeared, with only the main peaks remaining, which are considered as significant [62]. In the Discussion section, we identified and discussed three phase coupling peaks: the VLF peak, the HF respiratory peak and a shift to the BW2 region range in the presence of OSA. For the no-OSA and mild OSA children (Figure A7 and Figure A8), the bispectral power is mainly concentrated around the VLF range, and a shift coupling is presented to a slight degree in a moderate OSA subject (Figure A9) and markedly in a severe OSA subject (Figure A10). For severe OSA, in fact, the main bispectral peak moved to the OSA-specific BW2 region. These results are in accordance with the conclusions extracted in the Discussion section. Another remaining coupling peak after BWS correction that deserves mention is the respiratory peak (typically inside the HF region). This peak can be appreciated in Figure A7 and Figure A8. However, it is absent in the moderate and severe OSA children, supporting our original conclusions about the coupling reduction around the respiratory peak and the redistribution of the bispectral power with OSA severity, especially in the severe group. Therefore, the results obtained through the BWS method support our original conclusions, increasing the robustness of our results.

References

  1. Marcus, C.L.; Chapman, D.; Ward, S.D.; McColley, S.A.; Herrerias, C.T.; Stilwell, P.C.; Howenstine, M.; Light, M.J.; Schaeffer, D.A.; Wagener, J.S.; et al. Clinical practice guideline: Diagnosis and management of childhood. Pediatrics 2002, 109, 704–712. [Google Scholar]
  2. Society, A.T. Standards and indications for cardiopulmonary sleep studies in children. Am. J. Respir. Crit. Care Med. 1996, 153, 866–878. [Google Scholar] [CrossRef]
  3. Tauman, R.; Gozal, D. Obstructive sleep apnea syndrome in children. Expert Rev. Respir. Med. 2011, 5, 425–440. [Google Scholar] [CrossRef] [PubMed]
  4. Gozal, D. Sleep-Disordered Breathing and School Performance in Children. Pediatrics 1998, 102, 616–620. [Google Scholar] [CrossRef]
  5. Hunter, S.J.; Gozal, D.; Smith, D.L.; Philby, M.F.; Kaylegian, J.; Kheirandish-Gozal, L. Effect of sleep-disordered breathing severity on cognitive performance measures in a large community cohort of young school-aged children. Am. J. Respir. Crit. Care Med. 2016, 194, 739–747. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Kwok, K.L.; Ng, D.K.K.; Cheung, Y.F. BP and Arterial Distensibility in Children With Primary Snoring. Chest 2003, 123, 1561–1566. [Google Scholar] [CrossRef]
  7. Gozal, D.; Pope, D.W., Jr. Snoring During Early Childhood and Academic Performance at Ages Thirteen to Fourteen Years. Pediatrics 2001, 107, 1394–1399. [Google Scholar] [CrossRef] [Green Version]
  8. Iber, C.; Ancoli-Israel, S.; Chesson, A.L.; Quan, S.F. The AASM Manual for the Scoring of Sleep and Associated Events: Rules, Terminology and Technical Specifications, 1st ed.; American Academy Sleep Medicine: Westchester, IL, USA, 2007. [Google Scholar]
  9. Marcus, C.L.; Brooks, L.J.; Ward, S.D.; Draper, K.A.; Gozal, D.; Halbower, A.C.; Jones, J.; Lehmann, C.; Schechter, M.S.; Sheldon, S.; et al. Diagnosis and management of childhood obstructive sleep apnea syndrome. Pediatrics 2012, 130, e714–e755. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Berry, R.B.; Budhiraja, R.; Gottlieb, D.J.; Gozal, D.; Iber, C.; Kapur, V.K.; Marcus, C.L.; Mehra, R.; Parthasarathy, S.; Quan, S.F.; et al. Rules for Scoring Respiratory Events in Sleep: Update of the 2007 AASM Manual for the Scoring of Sleep and Associated Events. J. Clin. Sleep Med. 2012, 08, 597–619. [Google Scholar] [CrossRef] [Green Version]
  11. Tan, H.-L.; Gozal, D.; Ramirez, H.M.; Bandla, H.P.R.; Kheirandish-Gozal, L. Overnight polysomnography versus respiratory polygraphy in the diagnosis of pediatric obstructive sleep apnea. Sleep 2014, 37, 255–260. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Alonso-Álvarez, M.L.; Terán-Santos, J.; Ordax Carbajo, E.; Cordero-Guevara, J.A.; Navazo-Egüia, A.I.; Kheirandish-Gozal, L.; Gozal, D. Reliability of Home Respiratory Polygraphy for the Diagnosis of Sleep Apnea in Children. Chest 2015, 147, 1020–1028. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Chiner, E.; Cánovas, C.; Molina, V.; Sancho-Chust, J.N.; Vañes, S.; Pastor, E.; Martinez-Garcia, M.A. Home Respiratory Polygraphy is Useful in the Diagnosis of Childhood Obstructive Sleep Apnea Syndrome. J. Clin. Med. 2020, 9, 2067. [Google Scholar] [CrossRef] [PubMed]
  14. Gutiérrez-Tobal, G.C.; Álvarez, D.; Kheirandish-Gozal, L.; Campo, F.; Gozal, D.; Hornero, R. Reliability of machine learning to diagnose pediatric obstructive sleep apnea: Systematic review and meta-analysis. Pediatr. Pulmonol. 2021, 25423. [Google Scholar] [CrossRef]
  15. Vlahandonis, A.; Walter, L.M.; Horne, R.S.C. Does treatment of SDB in children improve cardiovascular outcome? Sleep Med. Rev. 2013, 17, 75–85. [Google Scholar] [CrossRef]
  16. Horne, R.S.C.; Shandler, G.; Tamanyan, K.; Weichard, A.; Odoi, A.; Biggs, S.N.; Davey, M.J.; Nixon, G.M.; Walter, L.M. The impact of sleep disordered breathing on cardiovascular health in overweight children. Sleep Med. 2018, 41, 58–68. [Google Scholar] [CrossRef]
  17. Vitelli, O.; Del Pozzo, M.; Baccari, G.; Rabasco, J.; Pietropaoli, N.; Barreto, M.; Villa, M.P. Autonomic imbalance during apneic episodes in pediatric obstructive sleep apnea. Clin. Neurophysiol. 2016, 127, 551–555. [Google Scholar] [CrossRef]
  18. Malik, M.; Bigger, J.T.; Camm, A.J.; Kleiger, R.E.; Malliani, A.; Moss, A.J.; Schwartz, P.J. Heart rate variability. Standards of measurement, physiological interpretation, and clinical use. Eur. Heart J. 1996, 17, 354–381. [Google Scholar] [CrossRef] [Green Version]
  19. Acharya, U.R.; Joseph, K.P.; Kannathal, N.; Lim, C.M.; Suri, J.S. Heart rate variability: A review. Med. Biol. Eng. Comput. 2006, 44, 1031–1051. [Google Scholar] [CrossRef]
  20. Guilleminault, C.; Winkle, R.; Connolly, S.; Melvin, K.; Tilkian, A. Cyclical Variation of the Heart Rate in Sleep Apnoea Syndrome. Lancet 1984, 323, 126–131. [Google Scholar] [CrossRef]
  21. Penzel, T.; Kantelhardt, J.W.; Grote, L.; Peter, J.H.; Bunde, A. Comparison of Detrended fluctuation Analysis and Spectral Analysis of Heart Rate Variability in Sleep and Sleep Apnea. IEEE Trans. Biomed. Eng. 2003, 50, 1143–1151. [Google Scholar] [CrossRef] [Green Version]
  22. Liao, D.; Li, X.; Vgontzas, A.N.; Liu, J.; Rodriguez-Colon, S.; Calhoun, S.; Bixler, E.O. Sleep-disordered breathing in children is associated with impairment of sleep stage-specific shift of cardiac autonomic modulation. J. Sleep Res. 2010, 19, 358–365. [Google Scholar] [CrossRef] [Green Version]
  23. Walter, L.M.; Nixon, G.M.; Davey, M.J.; Anderson, V.; Walker, A.M.; Horne, R.S.C. Autonomic dysfunction in children with sleep disordered breathing. Sleep Breath. 2013, 17, 605–613. [Google Scholar] [CrossRef]
  24. Shouldice, R.B.; O’Brien, L.M.; O’Brien, C.; De Chazal, P.; Gozal, D.; Heneghan, C. Detection of Obstructive Sleep Apnea in Pediatric Subjects using Surface Lead Electrocardiogram Features. Sleep 2004, 27, 784–792. [Google Scholar] [CrossRef] [Green Version]
  25. Kwok, K.L.; Yung, T.C.; Ng, D.K.; Chan, C.H.; Lau, W.F.; Fu, Y.M. Heart Rate Variability in Childhood Obstructive Sleep Apnea. Pediatr. Pulmonol. 2011, 46, 205–210. [Google Scholar] [CrossRef]
  26. Gil, E.; Mendez, M.; Vergara, J.M.; Cerutti, S.; Bianchi, A.M.; Laguna, P. Discrimination of sleep-apnea-related decreases in the amplitude fluctuations of ppg signal in children by HRV analysis. IEEE Trans. Biomed. Eng. 2009, 56, 1005–1014. [Google Scholar] [CrossRef]
  27. Martín-Montero, A.; Gutiérrez-Tobal, G.C.; Kheirandish-Gozal, L.; Jiménez-García, J.; Álvarez, D.; del Campo, F.; Gozal, D.; Hornero, R. Heart rate variability spectrum characteristics in children with sleep apnea. Pediatr. Res. 2021, 89, 1771–1779. [Google Scholar] [CrossRef] [PubMed]
  28. Baharav, A.; Kotagal, S.; Rubin, B.K.; Pratt, J.; Akselrod, S. Autonomic cardiovascular control in children with obstructive sleep apnea. Clin. Auton. Res. 1999, 9, 345–351. [Google Scholar] [CrossRef] [PubMed]
  29. Chua, K.C.; Chandran, V.; Acharya, U.R.; Lim, C.M. Application of higher order statistics/spectra in biomedical signals-A review. Med. Eng. Phys. 2010, 32, 679–689. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  30. Atri, R.; Mohebbi, M. Obstructive sleep apnea detection using spectrum and bispectrum analysis of single-lead ECG signal. Physiol. Meas. 2015, 36, 1963–1980. [Google Scholar] [CrossRef]
  31. Porta, A.; Bari, V.; Marchi, A.; De Maria, B.; Cysarz, D.; Van Leeuwen, P.; Takahashi, A.C.M.; Catai, A.M.; Gnecchi-Ruscone, T. Complexity analyses show two distinct types of nonlinear dynamics in short heart period variability recordings. Front. Physiol. 2015, 6, 71. [Google Scholar] [CrossRef] [Green Version]
  32. Martín-González, S.; Navarro-Mesa, J.L.; Juliá-Serdá, G.; Ramírez-Ávila, G.M.; Ravelo-García, A.G. Improving the understanding of sleep apnea characterization using Recurrence Quantification Analysis by defining overall acceptable values for the dimensionality of the system, the delay, and the distance threshold. PLoS ONE 2018, 13, e0194462. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Aljadeff, G.; Gozal, D.; Schechtman, V.L.; Burrell, B.; Harper, R.M.; Davidson Ward, S.L. Heart Rate Variability in Children With Obstructive Sleep Apnea. Sleep 1997, 20, 151–157. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Chua, K.C.; Chandran, V.; Acharya, U.R.; Lim, C.M. Cardiac state diagnosis using higher order spectra of heart rate variability. J. Med. Eng. Technol. 2008, 32, 145–155. [Google Scholar] [CrossRef]
  35. Chua, K.C. Cardiac Health Diagnosis Using Higher Order Spectra and Support Vector Machine. Open Med. Inform. J. 2009, 3, 1–8. [Google Scholar] [CrossRef] [Green Version]
  36. Saliu, S.; Birand, A.; Kudaiberdieva, G. Bispectral analysis of heart rate variability signal. In Proceedings of the 11th European Signal Processing Conference, Toulouse, France, 3–6 September 2002; IEEE: Piscataway, NJ, USA, 2015; Volume 2002. [Google Scholar]
  37. Yu, S.-N.; Lee, M.-Y. Bispectral analysis and genetic algorithm for congestive heart failure recognition based on heart rate variability. Comput. Biol. Med. 2012, 42, 816–825. [Google Scholar] [CrossRef] [PubMed]
  38. Garcia, R.G.; Valenza, G.; Tomaz, C.A.; Barbieri, R. Instantaneous bispectral analysis of heartbeat dynamics for the assessment of major depression. In Proceedings of the 2015 Computing in Cardiology Conference (CinC), Nice, France, 6–9 September 2015; IEEE: Piscataway, NJ, USA, 2016; Volume 42, pp. 781–784. [Google Scholar]
  39. Shao, S.; Wang, T.; Song, C.; Chen, X.; Cui, E.; Zhao, H. Obstructive Sleep Apnea Recognition Based on Multi-Bands Spectral Entropy Analysis of Short-Time Heart Rate Variability. Entropy 2019, 21, 812. [Google Scholar] [CrossRef] [Green Version]
  40. Zheng, L.; Pan, W.; Li, Y.; Luo, D.; Wang, Q.; Liu, G. Use of Mutual Information and Transfer Entropy to Assess Interaction between Parasympathetic and Sympathetic Activities of Nervous System from HRV. Entropy 2017, 19, 489. [Google Scholar] [CrossRef] [Green Version]
  41. Liu, D.; Yang, X.; Wang, G.; Ma, J.; Liu, Y.; Peng, C.-K.; Zhang, J.; Fang, J. HHT based cardiopulmonary coupling analysis for sleep apnea detection. Sleep Med. 2012, 13, 503–509. [Google Scholar] [CrossRef]
  42. Hornero, R.; Kheirandish-Gozal, L.; Gutiérrez-Tobal, G.C.; Philby, M.F.; Alonso-Álvarez, M.L.; Alvarez, D.; Dayyat, E.A.; Xu, Z.; Huang, Y.S.; Kakazu, M.T.; et al. Nocturnal oximetry-based evaluation of habitually snoring children. Am. J. Respir. Crit. Care Med. 2017, 196, 1591–1598. [Google Scholar] [CrossRef]
  43. Redline, S.; Amin, R.; Beebe, D.; Chervin, R.D.; Garetz, S.L.; Giordani, B.; Marcus, C.L.; Moore, R.H.; Rosen, C.L.; Arens, R.; et al. The Childhood Adenotonsillectomy Trial (CHAT): Rationale, Design, and Challenges of a Randomized Controlled Trial Evaluating a Standard Surgical Procedure in a Pediatric Population. Sleep 2011, 34, 1509–1517. [Google Scholar] [CrossRef] [Green Version]
  44. Marcus, C.L.; Moore, R.H.; Rosen, C.L.; Giordani, B.; Garetz, S.L.; Taylor, H.G.; Mitchell, R.B.; Amin, R.; Katz, E.S.; Arens, R.; et al. A randomized trial of adenotonsillectomy for childhood sleep apnea. N. Engl. J. Med. 2013, 368, 2366–2376. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Barroso-García, V.; Gutiérrez-Tobal, G.C.; Kheirandish-Gozal, L.; Álvarez, D.; Vaquerizo-Villar, F.; Núñez, P.; del Campo, F.; Gozal, D.; Hornero, R. Usefulness of recurrence plots from airflow recordings to aid in paediatric sleep apnoea diagnosis. Comput. Methods Programs Biomed. 2020, 183, 105083. [Google Scholar] [CrossRef]
  46. Barroso-García, V.; Gutiérrez-Tobal, G.C.; Kheirandish-Gozal, L.; Vaquerizo-Villar, F.; Álvarez, D.; del Campo, F.; Gozal, D.; Hornero, R. Bispectral analysis of overnight airflow to improve the pediatric sleep apnea diagnosis. Comput. Biol. Med. 2021, 129, 104167. [Google Scholar] [CrossRef] [PubMed]
  47. Barroso-García, V.; Gutiérrez-Tobal, G.C.; Kheirandish-Gozal, L.; Álvarez, D.; Vaquerizo-Villar, F.; Crespo, A.; del Campo, F.; Gozal, D.; Hornero, R. Irregularity and variability analysis of airflow recordings to facilitate the diagnosis of paediatric sleep apnoea-hypopnoea syndrome. Entropy 2017, 19, 447. [Google Scholar] [CrossRef] [Green Version]
  48. Benitez, D.; Gaydecki, P.A.; Zaidi, A.; Fitzpatrick, A.P. The use of the Hilbert transform in ECG signal analysis. Comput. Biol. Med. 2001, 31, 399–406. [Google Scholar] [CrossRef]
  49. Gutiérrez-Tobal, G.C.; Álvarez, D.; Gomez-Pilar, J.; Del Campo, F.; Hornero, R. Assessment of time and frequency domain entropies to detect sleep apnoea in heart rate variability recordings from men and women. Entropy 2015, 17, 123–141. [Google Scholar] [CrossRef] [Green Version]
  50. Yu, L.; Liu, H. Efficient Feature Selection via Analysis of Relevance and Redundancy. J. Mach. Learn. Res. 2004, 5, 1205–1224. [Google Scholar]
  51. Vaquerizo-Villar, F.; Álvarez, D.; Kheirandish-Gozal, L.; Gutiérrez-Tobal, G.C.; Barroso-García, V.; Crespo, A.; del Campo, F.; Gozal, D.; Hornero, R. Utility of bispectrum in the screening of pediatric sleep apnea-hypopnea syndrome using oximetry recordings. Comput. Methods Programs Biomed. 2018, 156, 141–149. [Google Scholar] [CrossRef]
  52. Guyon, I.; Elisseeff, A.; Kaelbling, L.P. An introduction to variable and feature selection. J. Mach. Learn. Res. 2003, 3, 1157–1182. [Google Scholar]
  53. Bishop, C. Neural Networks for Pattern Recognition; Oxford University Press: Oxford, UK, 1995. [Google Scholar]
  54. Wang, R.; Wang, J.; Li, S.; Yu, H.; Deng, B.; Wei, X. Multiple feature extraction and classification of electroencephalograph signal for Alzheimers’ with spectrum and bispectrum. Chaos Interdiscip. J. Nonlinear Sci. 2015, 25, 013110. [Google Scholar] [CrossRef] [PubMed]
  55. Gil, E.; Bailón, R.; Vergara, J.M.; Laguna, P. PTT variability for discrimination of sleep apnea related decreases in the amplitude fluctuations of PPG signal in children. IEEE Trans. Biomed. Eng. 2010, 57, 1079–1088. [Google Scholar] [CrossRef]
  56. Lázaro, J.; Gil, E.; Vergara, J.M.; Laguna, P. Pulse rate variability analysis for discrimination of sleep-apnea-related decreases in the amplitude fluctuations of pulse photoplethysmographic signal in children. IEEE J. Biomed. Health Inform. 2014, 18, 240–246. [Google Scholar] [CrossRef] [PubMed]
  57. Cohen, G.; de Chazal, P. Automated detection of sleep apnea in infants: A multi-modal approach. Comput. Biol. Med. 2015, 63, 118–123. [Google Scholar] [CrossRef]
  58. Vaquerizo-Villar, F.; Alvarez, D.; Kheirandish-Gozal, L.; Gutierrez-Tobal, G.C.; Barroso-Garcia, V.; Santamaria-Vazquez, E.; Del Campo, F.; Gozal, D.; Hornero, R. A convolutional neural network architecture to enhance oximetry ability to diagnose pediatric obstructive sleep apnea. IEEE J. Biomed. Health Inform. 2021. Early Access. [Google Scholar] [CrossRef] [PubMed]
  59. Milagro, J.; Gracia, J.; Seppa, V.-P.; Karjalainen, J.; Paassilta, M.; Orini, M.; Bailon, R.; Gil, E.; Viik, J. Noninvasive Cardiorespiratory Signals Analysis for Asthma Evolution Monitoring in Preschool Children. IEEE Trans. Biomed. Eng. 2019, 67, 1863–1871. [Google Scholar] [CrossRef] [PubMed]
  60. Milagro, J.; Gil, E.; Lazaro, J.; Seppa, V.P.; Pekka Malmberg, L.; Pelkonen, A.S.; Kotaniemi-Syrjanen, A.; Makela, M.J.; Viik, J.; Bailon, R. Nocturnal Heart Rate Variability Spectrum Characterization in Preschool Children with Asthmatic Symptoms. IEEE J. Biomed. Heal. Inform. 2018, 22, 1332–1340. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  61. Porta, A.; Faes, L. Wiener–Granger Causality in Network Physiology With Applications to Cardiovascular Control and Neuroscience. Proc. IEEE 2016, 104, 282–309. [Google Scholar] [CrossRef]
  62. Siu, K.L.; Ahn, J.M.; Ju, K.; Lee, M.; Shin, K.; Chon, K.H. Statistical Approach to Quantify the Presence of Phase Coupling Using the Bispectrum. IEEE Trans. Biomed. Eng. 2008, 55, 1512–1520. [Google Scholar] [CrossRef]
  63. Theiler, J.; Eubank, S.; Longtin, A.; Galdrikian, B.; Doyne Farmer, J. Testing for nonlinearity in time series: The method of surrogate data. Phys. D Nonlinear Phenom. 1992, 58, 77–94. [Google Scholar] [CrossRef] [Green Version]
  64. Maestri, R.; Pinna, G.D.; Porta, A.; Balocchi, R.; Sassi, R.; Signorini, M.G.; Dudziak, M.; Raczak, G. Assessing nonlinear properties of heart rate variability from short-term recordings: Are these measurements reliable? Physiol. Meas. 2007, 28, 1067–1077. [Google Scholar] [CrossRef]
Figure 1. Averaged bispectrum magnitude in the range 0–0.4 Hz in the training set for the four severity groups. To improve the visualization of the coupling focus at very low frequencies, an amplification between 0 and 0.05 Hz is depicted in the upper right corner for each figure. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Figure 1. Averaged bispectrum magnitude in the range 0–0.4 Hz in the training set for the four severity groups. To improve the visualization of the coupling focus at very low frequencies, an amplification between 0 and 0.05 Hz is depicted in the upper right corner for each figure. (a) No-OSA; (b) mild OSA; (c) moderate OSA; (d) severe OSA.
Entropy 23 01016 g001
Figure 2. Features selected in each feature set after applying FCBF algorithm over 1000 bootstrap replicates of the training set. (a) Features selected in the classic region set; (b) features selected in the OSA-specific region set.
Figure 2. Features selected in each feature set after applying FCBF algorithm over 1000 bootstrap replicates of the training set. (a) Features selected in the classic region set; (b) features selected in the OSA-specific region set.
Entropy 23 01016 g002
Figure 3. Boxplot distribution of the features selected in the bispectral classic region feature subset for the four OSA severity groups in the training set. The p-value obtained with the Kruskal–Wallis test is shown in each subplot. (a) VLF_f2m boxplots and p-value; (b) LF_BE2 boxplots and p-value; (c) HF_PE boxplots and p-value.
Figure 3. Boxplot distribution of the features selected in the bispectral classic region feature subset for the four OSA severity groups in the training set. The p-value obtained with the Kruskal–Wallis test is shown in each subplot. (a) VLF_f2m boxplots and p-value; (b) LF_BE2 boxplots and p-value; (c) HF_PE boxplots and p-value.
Entropy 23 01016 g003
Figure 4. Boxplot distribution of the features selected in the bispectral OSA-specific region feature subset for the four OSA severity groups in the training set. The p-value obtained with the Kruskal–Wallis test is shown in each subplot. (a) BW2_RPDiag boxplots and p-value; (b) BW2_BE1 boxplots and p-value; (c) BWRes_Bmin boxplots and p-value; (d) BWRes_BE3 boxplots and p-value.
Figure 4. Boxplot distribution of the features selected in the bispectral OSA-specific region feature subset for the four OSA severity groups in the training set. The p-value obtained with the Kruskal–Wallis test is shown in each subplot. (a) BW2_RPDiag boxplots and p-value; (b) BW2_BE1 boxplots and p-value; (c) BWRes_Bmin boxplots and p-value; (d) BWRes_BE3 boxplots and p-value.
Entropy 23 01016 g004
Table 1. Clinical and demographic data from children included in this study.
Table 1. Clinical and demographic data from children included in this study.
AllTraining Set (UofC)Test Set (CHAT)
Subjects (n)1738981757
Age (years)6.4 [3.3]6.0 [6.0]7.0 [2.4]
Males (n)962 (55.35%)602 (61.37%)360 (47.95%)
BMI (kg/m2)17.63 [5.37]18.02 [5.86]17.28 [4.64]
AHI (e/h)2.23 [5.27]3.8 [7.76]1.46 [2.07]
AHI ≥ 1 (e/h)1309 (75.31%)808 (82.36%)501 (66.18%)
AHI ≥ 5 (e/h)519 (29.86%)407 (41.49%)112 (14.80%)
AHI ≥ 10 (e/h)298 (17.15%)229 (23.34%)69 (9.11%)
Data are shown as median [interquartile range] or n (percentage). BMI: body mass index; AHI: Apnea–Hypopnea Index; UofC: University of Chicago; CHAT: Childhood Adenotonsillectomy Trial.
Table 2. Summary of the bispectral features initially computed in each region. Features related to the diagonal of the region were excluded in the BWRes region.
Table 2. Summary of the bispectral features initially computed in each region. Features related to the diagonal of the region were excluded in the BWRes region.
Classic Region Feature SetOSA-Specific Region Feature Set
FeaturesVLFLFHFBW1BW2BWRes
RPDiagVLF_RPDiagLF_RPDiagHF_RPDiagBW1_RPDiagBW2_RPDiag-
BmaxVLF_BmaxLF_BmaxHF_BmaxBW1_BmaxBW2_BmaxBWRes_Bmax
BminVLF_BminLF_BminHF_BminBW1_BminBW2_BminBWRes_Bmin
BTotalVLF_BTotalLF_BTotalHF_BTotalBW1_BTotalBW2_BTotalBWRes_BTotal
BE1VLF_BE1LF_BE1HF_BE1BW1_BE1BW2_BE1BWRes_BE1
BE2VLF_BE2LF_BE2HF_BE2BW1_BE2BW2_BE2BWRes_BE2
BE3VLF_BE3LF_BE3HF_BE3BW1_BE3BW2_BE3BWRes_BE3
PEVLF_PELF_PEHF_PEBW1_PEBW2_PEBWRes_PE
H1VLF_H1LF_H1HF_H1BW1_H1BW2_H1BWRes_H1
H2VLF_H2LF_H2HF_H2BW1_H2BW2_H2-
H3VLF_H3LF_H3HF_H3BW1_H3BW2_H3-
H4VLF_H4LF_H4HF_H4BW1_H4BW2_H4-
f1mVLF_f1mLF_f1mHF_f1mBW1_f1mBW2_f1mBWRes_f1m
f2mVLF_f2mLF_f2mHF_f2mBW1_f2mBW2_f2mBWRes_f2m
Table 3. Results of the partial correlation study in the test set between features selected for each subset and the polysomnographic indices considered.
Table 3. Results of the partial correlation study in the test set between features selected for each subset and the polysomnographic indices considered.
BISPClassic Features
PSG IndexVLF_f2mLF_BE2HF_PE
ρSp-ValueρSp-ValueρSp-Value
AHI0.274<<0.01−0.185<<0.01−0.1120.002 *
OAHI0.261<<0.01−0.149<<0.01−0.0970.008
OAI0.167<<0.01−0.1050.004 *−0.0640.079
ODI0.215<<0.01−0.1230.001 *−0.0540.138
#Awakenings−0.0750.039−0.0270.461−0.0200.586
WASO−0.0030.9290.0650.076−0.0220.538
%N10.0890.014−0.0710.052−0.0300.404
%N2−0.0340.3570.0990.007 *0.0130.715
%N30.0340.355−0.0250.497−0.0440.23
%REM−0.1250.001−0.0520.1540.0590.108
TAI0.213<<0.01−0.158<<0.01−0.1150.002 *
BISPSpecific Features
PSG IndexBW2_RPDiagBW2_BE1BWRes_BminBWRes_BE3
ρSp-ValueρSp-ValueρSp-ValueρSp-Value
AHI0.308<<0.01−0.180<<0.010.0540.1360.0450.214
OAHI0.261<<0.01−0.180<<0.010.0980.007 *0.0280.435
OAI0.177<<0.01−0.173<<0.010.0710.0510.0580.112
ODI0.247<<0.01−0.1390.0010.0190.610.0720.047
#Awakenings−0.0330.372−0.0010.994−0.0060.8760.0350.331
WASO0.0710.050.0780.031−0.0180.6220.0560.126
%N10.1070.003 *−0.0610.0930.0230.5270.0280.441
%N2−0.0610.0910.0080.8370.0480.1840.0590.104
%N30.0530.1470.0080.817−0.0750.04−0.0920.011
%REM−0.1390.0010.0480.192−0.0070.8550.0130.722
TAI0.225<<0.01−0.144<<0.010.0680.0620.0680.06
PSG: polysomnographic; AHI: Apnea–Hypopnea Index; OAHI: Obstructive AHI; OAI: obstructive apnea index; ODI: oxygen desaturation index; #Awakenings: number of awakenings during total sleep time; WASO: wake after sleep onset; %N1: percentage of sleep spent in N1; %N2: percentage of sleep spent in N2; %N3: percentage of sleep spent in N3; %REM: percentage of sleep spent in REM; TAI: total arousal index. Those p-values below 10−4 appear as << 0.01. * Non-significant after Bonferroni correction. Statistically significant correlations (p-value < 0.01 after Bonferroni correction) appear in bold.
Table 4. Diagnostic performance achieved in the test set by each feature selected and each MLP optimized model for the binary classification in each severity threshold. Results are shown in terms of sensitivity (Se %), specificity (Sp %), accuracy (Acc %) and AUC.
Table 4. Diagnostic performance achieved in the test set by each feature selected and each MLP optimized model for the binary classification in each severity threshold. Results are shown in terms of sensitivity (Se %), specificity (Sp %), accuracy (Acc %) and AUC.
Threshold: AHI = 1 e/h
Feature/ModelSeSpAccAUC
VLF_f2m44.572.353.90.605
LF_BE242.172.752.40.581
HF_PE42.963.349.80.550
BW2_RPDiag50.964.855.60.629
BW2_BE147.159.451.30.559
BWRes_Bmin40.557.446.20.482
BWRes_BE341.557.446.90.513
MLP1Classic52.359.454.70.600
MLP1Specific76.338.363.40.627
Threshold: AHI = 5 e/h
Feature/ModelSeSpAccAUC
VLF_f2m62.572.270.80.749
LF_BE256.374.471.70.670
HF_PE45.572.168.20.628
BW2_RPDiag60.777.775.20.747
BW2_BE156.370.168.00.671
BWRes_Bmin58.945.347.30.567
BWRes_BE347.358.456.80.569
MLP5Classic50.986.281.00.774
MLP5Specific62.584.281.00.791
Threshold: AHI = 10 e/h
Feature/ModelSeSpAccAUC
VLF_f2m63.876.775.60.784
LF_BE258.081.579.40.740
HF_PE53.672.170.40.663
BW2_RPDiag68.176.075.30.789
BW2_BE147.876.073.40.692
BWRes_Bmin56.550.651.10.557
BWRes_BE355.159.459.00.614
MLP10Classic43.596.591.70.847
MLP10Specific66.791.689.30.841
VLF: very low frequency; LF: low frequency; HF: high frequency; MLP: multi-layer perceptron; AHI: Apnea–Hypopnea Index. The highest ACC and AUC for each severity threshold are highlighted in bold.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Martín-Montero, A.; Gutiérrez-Tobal, G.C.; Gozal, D.; Barroso-García, V.; Álvarez, D.; del Campo, F.; Kheirandish-Gozal, L.; Hornero, R. Bispectral Analysis of Heart Rate Variability to Characterize and Help Diagnose Pediatric Sleep Apnea. Entropy 2021, 23, 1016. https://doi.org/10.3390/e23081016

AMA Style

Martín-Montero A, Gutiérrez-Tobal GC, Gozal D, Barroso-García V, Álvarez D, del Campo F, Kheirandish-Gozal L, Hornero R. Bispectral Analysis of Heart Rate Variability to Characterize and Help Diagnose Pediatric Sleep Apnea. Entropy. 2021; 23(8):1016. https://doi.org/10.3390/e23081016

Chicago/Turabian Style

Martín-Montero, Adrián, Gonzalo C. Gutiérrez-Tobal, David Gozal, Verónica Barroso-García, Daniel Álvarez, Félix del Campo, Leila Kheirandish-Gozal, and Roberto Hornero. 2021. "Bispectral Analysis of Heart Rate Variability to Characterize and Help Diagnose Pediatric Sleep Apnea" Entropy 23, no. 8: 1016. https://doi.org/10.3390/e23081016

APA Style

Martín-Montero, A., Gutiérrez-Tobal, G. C., Gozal, D., Barroso-García, V., Álvarez, D., del Campo, F., Kheirandish-Gozal, L., & Hornero, R. (2021). Bispectral Analysis of Heart Rate Variability to Characterize and Help Diagnose Pediatric Sleep Apnea. Entropy, 23(8), 1016. https://doi.org/10.3390/e23081016

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop