Neonatal apnea and hypopnea prediction in infants with Robin sequence with neural additive models for time series

Our Contribution

We built on recent work in neural additive models [16], which form their prediction by summing over the output of several neural networks. The additive nature allows users to inspect the magnitude of the additive contributions for each prediction to gauge the importance of the underlying feature or signal modality. This inherent interpretability makes neural additive models a good choice for clinical prediction tasks like individual apnea or hypopnea prediction. Since our data came in the form of time series, we replaced classically-used neural network architectures by architectures which are tailored to time series. We show that our model achieved a high level of prediction, with an average area under curve of 0.80 for the task of apnea and hypopnea prediction for infants with Robin sequence. Furthermore, our model allowed us to gauge the importance of different signal modalities for a performed prediction as well as to localise discriminative features within the signals.

Neonatal dataset and problem setup

We performed our analyses on a dataset of n = 19 infants with Robin sequence who had been admitted to the Tübingen neonatology department between May 2020 and April 2021. Gestational age at birth was 39 (32–41) [median (range)] weeks and birth weight 3,390 (1,320–4,380) g. The recording methods are described in detail in [17]. At the time of conducting the respiratory polygraphy, infants were 17 (1–73) days old and weighed 3,392 (2,642–4,380) g. In total, 185 hours of respiratory polygraphy data were recorded, containing 122 hours of total sleep time. Adverse events, namely obstructive, mixed and central apneas, and hypopneas, as well as infant movement were annotated via visual inspection by a domain expert after the recording. Infants experienced 27 (3–112) obstructive, mixed and central apnea events per hour, and 28 (1–58) obstructive and mixed hypopnea events per hour, with events lasting for 3.8 (3.9–5.7) and 4.9 (2.8–7.3) seconds per event, respectively. After removing parts of the signals that were annotated as either adverse events or movements, signals were divided into 30 second time windows. The target time windows were taken with an offset of 15 seconds, whereas the control time windows were required to be at least three minutes away from both the start and end of any adverse event. This procedure resulted in 214 (75-606) time windows per infant to be used for classification. For most patients there were more control than target time windows. The imbalance (percentage of target time windows) was 36% (4%-66%). The overall setup is depicted in Fig. 1 and is described in more detail in the Method Section.

• Download figure
• Open in new tab

Figure 1: Overview of the setup and core results.

a) 19 overnight respiratory polygraph recordings of infants with Robin sequence were collected. Six signals were used for the analysis: Nasal pressure (NP), the sum of thoracic and abdominal respiratory effort (T+A), photoplethysmogram (PPG), heart rate (HR), SpO₂, and PCO₂ levels. b) Adverse events (apneas and hypopneas, indicated as pink vertical boxes) as well as infant movement (blue boxes) were annotated manually. From these annotations, 30 second target time windows shortly preceding the adverse events (red), as well as control time windows (green) were extracted. c) Example of an extracted target time window. d) To discriminate between target and control time windows, a Neural Additive Model (NAM) [16] was trained and tested with patient-based leave-one-out cross validation. e) Classification performance of the NAM: Pointwise average over the ROC curves of individual patients. Uncertainty corresponds to the pointwise standard deviation. f) Area under ROC values for the individual patient-based test sets sorted in increasing order. Random performance is indicated in dark green. g) The average test performance as a function of an increasing prediction horizon. Uncertainty corresponds to the standard deviation over the performance of different patients.

Performance of neonatal adverse event prediction

We constructed a time-series based Neural Additive Model (NAM) [16] to perform the classification into target (pre-adverse event) and control time windows based on seven signal modalities: Nasal pressure, thoracic and abdominal respiratory effort, heart rate, SpO₂, and PCO₂ levels. Because of their strong correlation, thoracic and abdominal respiratory effort were summed to obtain a single signal, resulting in six modalities. We measured the performance of our classifier by performing patient-based leave-one-out cross validation. The average area under the receiver operating characteristic curve (AuROC) over all leave-one-out test sets was 0.80 with a standard deviation of 0.093 (Fig. 1). For some patients, we achieved a classification performance of over 0.9 AuROC. We conducted permutation tests [18] to assess that our obtained test performances are significant. For all but one patient (p_per = 0.121), the performance was significantly better than random (p_per < 0.001).

We also investigated how much performance drops as the prediction horizon between target time window and adverse event is increased or decreased. Fig. 1f shows the average AuROC for completely retrained classifiers as a function of the prediction horizon. As expected, for an increasing prediction horizon, the classification performance decreased and reached near chance level after about 150 seconds. On the other hand, for no offset between annotated events and target time windows, the performance further increased.

Neural Additive Model for signal modality importance

Our neural additive model allowed us to investigate how predictions are formed by inspecting the corresponding additive contributions. A key question was to understand which of the six physiological signal modalities carry predictive information. Identifying which signals carry the most predictive information could enable both neonatologists and possible future automated prediction systems to focus on fewer signal modalities, thus reducing strain on the infants. The NAM classifier allowed us to investigate the importance of different signal modalities by analysing the additive contribution of each signal modality to the overall classification score. This analysis was possible because, by construction, the overall classification score is a sum of the individual additive contributions from each signal modality. The higher an individual additive contribution is compared to the other contributions, the more influence it has on the overall classification. To gauge overall importance of the six signal modalities used for classification, we analysed the additive contributions pooled over all patient test sets (Fig. 2a). To analyse differences between patients, we also computed the standard deviations over the additive contributions of individual patients (Fig. 2b). Furthermore, we performed a single modality analysis, where the subnetwork of each signal modality is trained independently. To analyse the significance of potential differences between modalities, we repeated our training procedure five times and computed Wilcoxon signed-rank (p_wil) tests over the n = 19 patients. See the Appendix for more details about the statistical tests.

• Download figure
• Open in new tab

Figure 2: Importance analysis of the signal modalities.

For the NAM, the importance of a signal modality can be measured by analysing its associated additive contributions. a) Histograms of additive contributions pooled over all patient-based test sets. b) Box plots of the distribution over the standard deviation of additive contributions of all individual patients. c) Average AuROC scores of the individually trained single modality networks together with the standard deviation. All analyses show that the nasal pressure signal contains the most predictive information on average.

• Download figure
• Open in new tab

Figure 3: Exemplary activation maps together with the associated additive contributions.

a) Correctly classified target (pre-adverse-event) time window. Light blue indicates negative activation and a classification towards no adverse event. Pink indicates positive activation and a classification towards an upcoming adverse event. The classifier detected different types of irregularities in the breathing signal as well as variations in heart rate and SpO₂ levels. b) Correctly classified control time window c) Incorrectly classified control time window (false positive). d) Incorrectly classified target time window (false negative).

The large standard deviations over both the pooled and patient based additive contributions indicated that the nasal pressure signal contained the most predictive information (Fig. 2a/b). The other breathing-related signal, the sum of thoracic and abdominal respiratory effort, was comparatively less informative (p_wil < 0.001). This result was confirmed by the performances of the single modality networks: The network for nasal pressure achieved the highest average AuROC, which was significantly better than the performance achieved by the network for thoracic and abdominal respiratory effort (p_wil = 0.032) (Fig. 2c). The importance of both breathing-related signals was in accordance with the clinical perspective as infants with Robin sequence mainly suffer from obstructive apnea and hypopnea. Furthermore, both the NAM and single modality analysis showed that, apart from the breathing-related signals, the other modalities PPG, heart rate, PCO₂ and especially SpO₂ also carried predictive information. The SpO₂ single modality performance was on par with the nasal pressure (no significant difference: p_wil = 0.241). The informativeness of SpO₂ was expected clinically: apneas and hypopneas often appear in clusters and hence decreased SpO₂ levels may indicate immediately preceding adverse events which increase the probability of another upcoming adverse event. A similar effect could be observed for the heart rate as apneas and hypopneas are often followed by bradycardia. Furthermore, the PPG signal is very sensitive to small movements, and thus captured pre-apnea or pre-hypopnea arousal.

Although some of the signal modalities contained less information than others, additively combining their information increased overall performance: The NAM achieved significantly better average AuROC than each of the trained single modality networks (p_wil = 0.009).

Activation maps for local interpretability

While analysing the additive contributions of the NAM revealed which of the signal modalities contained the most predictive information, it could not answer what parts or features of the respective signal modalities were discriminative for an upcoming event. In other words, we wanted to answer the question: What does the classifier “look at” within the signal modalities? As we input whole time series into our classifier, we lost global interpretability, that is, the ability to visualise the full regression function, which is usually provided in generalised additive models. We could, however, achieve local interpretability. Local interpretability refers to the the explanation for the prediction made with respect to an individual data point. To this end, the fully convolutional subnetworks of the NAM allowed us to use convolutional activation maps to investigate which part of the time window activates the convolutional filters.

In our case, each time window classified by the NAM resulted in a total of six activation maps and additive contributions, one for each signal modality. Additive contributions were simply computed as the average activation across time. Positive activations indicated that the corresponding segment of the signal was discriminative towards a target time window. In contrast, negative activations indicated that a segment of the signal was discriminative towards a control time window. Based on a visual analysis, the classifier focused on irregularities in the breathing signal both for the nasal pressure and the sum of the thoracic and abdominal respiratory efforts. The detected irregularities in the nasal pressure could include general change in respiratory rate and amplitude or a “mini-apnea” (flat nasal pressure signal not meeting annotation criteria). On top of these features, the classifier used arousal in the plethysmograph, reduced SpO₂ levels or SpO₂ desaturations as well as variations in the heart rate to make a classification towards a pre-adverse-event time window.

Limitations

While this study presents a first step in the direction of individual neonatal apnea and hypopnea prediction, there are several obstacles that hinder the application of the presented model in a real-world clinical setting. The ultimate goal would be to create a system that monitors an infant, predicts adverse events, and automatically triggers an intervention (that is, a closed loop apnea management system).

The presented model is not yet well suited for this task. For some patients, the prediction performance was not sufficient to be of practical use. Especially for systems that perform automated interventions, high prediction accuracy will be necessary as mistakes in the clinical setting are costly. Note that we evaluated performance by training and validating on 18 of the 19 infants and testing on the one left out. Following a domain adaptation approach, where the classifier adapts to the current infant in an online manner, could potentially mitigate the low performance for some of the neonates. Furthermore, additional measures would need to be taken to ensure the robustness of the model to unseen data, outliers and malfunctioning recording equipment to guarantee safe deployment.

There is also room for improvement regarding the interpretability of the classifier. As presented, the classifier did not discriminate between different predictive features within the time series that occur before an adverse event. A model that can provide more context about these predictive features will further increase the trust in its ability to provide correct predictions.

More generally, the current NAM classifier ignores possible interactions between time series. This shortcoming, common to all generalized additive models (GAMs), can lower the model expressiveness and cause concurvity issues that negatively affect the interpretability and stability of the model. There is work on including two-feature interaction in tabular GAMs [19], but it is currently unclear how to extend this to the time-series based NAM framework without losing some of the inherent interpretability.

Conclusion

Despite its limitations, our neural additive model presents a first step towards automated neonatal apnea and hypopnea prediction. When performing reliably, such a model would have the potential to optimise care for vulnerable neonates with upper airway obstruction.

Data collection

We based our study on a physiological dataset from infants with Robin sequence. This group is well suited for automated prediction because of their homogeneous pathophysiology, that is, apneas and hypopneas are mainly caused by a narrow upper airway. Between May 2020 and April 2021, 19 infants with Robin sequence underwent whole night recordings in the Department of Neonatology at Tübingen University Hospital, using standard digital cardiorespiratory polysomnography (Remlogic, Natus Medical Incorporated, California, USA). Modalities recorded included nasal airflow obtained via a nasal pressure transducer and recorded at 200Hz, as well as thoracic and abdominal respiratory efforts via respiratory inductance plethysmography recorded at 50Hz. Furthermore, the heart was monitored via an EKG and PPG recorded at 200Hz and 100Hz respectively. The beat-to-beat heart rate was automatically derived from the EKG signal based on the RR-intervals. Finally, the blood SpO₂ obtained from pulse oximetry and transcutaneous PCO₂ readings were recorded. Manual annotations of adverse events were derived from the recordings [17]. Five different types of adverse events were annotated by a domain expert according to the criteria of the 2020 American Academy of Sleep Medicine (AASM) guidelines. These were central apnea, obstructive apnea, and mixed apnea as well as central and obstructive hypopnea. Additionally, intermittent hypoxia events, defined as a fall in SpO₂ levels by more than three percentage points within a duration of five seconds (desaturations) and body movements were annotated. The movement annotation was based on a low-resolution infrared video frame. Details about the annotation criteria can be found in the Appendix.

General Setup

We framed the problem of apnea and hypopnea prediction as a binary time-series classification problem. All types of apnea and hypopnea were combined into one single type of adverse event. We made no attempt to differentiate between different types of apneas and hypopneas. Intermittent hypoxia events were not considered as adverse events. Parts of the signal that were annotated as either an adverse event or movement were removed from the signal. The remaining signal was then divided into 30 second non-overlapping time windows. If a time window preceded an adverse event with an offset of 15 seconds, it was labeled as a target time window. This offset was well above the maximal duration of the annotated apneas and hypopneas, and thus ensured that our prediction was not based on imprecise annotation time stamps. If a time window was at least three minutes away from both the start and end of an adverse event, it was labeled as a control time window (Fig. 1).

As mentioned in the Results Section, the number of extracted time windows ranged between 75 and 606 time windows with a median number of 214. For all but one patient, the time window generation resulted in more control time windows than pre-adverse-event time windows. The imbalance ranged between 4% and 66% with a median imbalance of 36%. To avoid training issues due to class imbalance, we undersampled the control time windows on a patient basis during training to balance control and target time windows. The classifier performance was evaluated on the full (imbalanced) test sets. To assess the average performance of our classifier, we performed patient-based leave-one-out cross validation. That is, we used all generated time windows from a single patient as a test set and trained on the union of all time windows of the remaining patients. This procedure was repeated for every patient and average performance scores and standard deviations are reported. For every repetition, we selected the hyperparameters with nested cross validation [20]: All time windows from every patient apart from the test patient time windows were used as a validation set. After exploring different hyperparameter configurations, we picked the hyperparameters that resulted in the highest average validation AuROC for training. More details are given in the Section “Architecture and training details”.

We focused on the importance of six signal modalities: Nasal pressure, the sum of thoracic and abdominal respiratory effort, heart rate (which is derived from the ECG), PPG, SpO₂, and PCO₂ signals. Inter-breath intervals based on thoracic respiratory effort and heart rate have “classically” been used to predict apnea and hypopnea [6]. Because thoracic and abdominal respiratory effort are highly correlated signals, we treated them as a single signal by summing both signals traces pointwise. Creating the sum of these two modalities is a well-established practice in the analysis of apnea and hypopnea episodes [21, 22]. Additionally, nasal pressure presented another interesting breathing-related signal that potentially carries different information than the “classically” used thoracic and abdominal respiratory effort. Moreover, we included the PPG signal and the two blood related parameters SpO₂ and PCO₂.

Preprocessing

The three breathing-related signals nasal pressure, thoracic and abdominal respiratory effort, sampled at 200Hz or 50 Hz respectively, were downsampled to 5 Hz after applying an anti-aliasing 8th order Chebyshev filter. The same preprocessing was applied to the PPG signal. Both the SpO₂ and PCO₂ signals were downsampled from 2Hz to 1Hz. The derived heart rate was left unchanged at 1Hz. We then standardised all breathing-related signals and the PPG signal on a time window basis. The summed signal of the thoracic and abdominal respiratory effort was obtained by summing both standardised signals pointwise and standardising the result again. The heart rate as well as SpO₂ and PCO₂ were range-normalised from 50 to 240bpm, 60-100%, and 30-70 mmHg to a range of -1 to 1. Our rationale for standardising the breathing-related and PPG modalities but range-normalising the heart rate, SpO₂ and PCO₂ signals, was that signal values from first group of modalities can vary strongly between patients. This is not true for the latter group of modalities. Apart from downsampling and standardising or normalising, we performed no additional preprocessing or filtering, allowing the capture of systematic noise in the signals, which might be predictive for an upcoming adverse event.

Periodic breathing

The choice of 30 second time windows was not arbitrary. Since we did not perform padding for shorter time windows, the fixed length ensured a minimal distance between annotated adverse events. This fixed length was necessary to avoid the very short inter-apnea episodes that occur during periodic breathing and that are easily distinguishable from normal breathing [23, 6]. With the choice of 30 seconds, we were above the definition of Kelly and Shannon [24], who define periodic breathing as “three or more episodes of central apnea lasting at least 4 seconds, each separated by no more than 20 seconds of normal breathing”.

Interpretability

In medical adverse event prediction, there are usually several signal modalities available that might be used to predict adverse events. Furthermore, data is often noisy and heterogeneous. It is therefore important to have some level of interpretability to avoid unforeseen shortcuts and failure modes in the prediction model. Our case of neonatal apnea and hypopnea prediction was no different in this regard. Since there was little prior work on possible features for individual apnea and hypopnea prediction [10], we performed automatic feature extraction with deep neural networks. Despite using deep neural networks, we required our classifier to be inherently interpretable. There were two questions we wanted our classifier to be able to answer:

• Question 1: Which physiological signal modalities of the multivariate polygraphy signal carry the strongest information to make a prediction?

• Question 2: What features within the signal modality traces themselves are discriminative for the prediction?

The first question could be answered by inspecting the additive contributions of a prediction made by a neural additive model and the second by generating activation maps. Both approaches were easily unified into a single architecture, as we describe in the following.

Neural Additive Models (Question 1)

Generalised Additive Models (GAMs) have historically been used to analyse tabular data especially in the medical domain. They combine expressiveness with built in interpretability [25, 19, 26]. For tabular data, GAMs are a linear combination of non-linear functions f_i, which each take one feature of the data as input. Linear and logistic regression are two special cases of GAMs where the f_i are constrained to be linear functions. The basic GAM for binary classification is given by

Here denotes the estimated probability that the input corresponds to a positive sample and σ(x) = (1 + e^x)⁻¹ denotes the sigmoid function. GAMs were traditionally fitted with splines or decision trees, but very recent work parameterised the functions f_i by simple multi-layer perceptrons giving rise to the so called Neural Additive Models (NAM) [16]. GAMs, and as a special case NAMs, are inherently interpretable because it is possible to visualise the individual functions f_i as a function of their corresponding feature. This allows users to globally interpret the model. Moreover, given a specific data sample x_i, it is possible to analyse the individual importance of each feature by comparing the contribution, that is, the individual summands α_if_i(x_i) to the overall classification score.

Unlike classical GAMs, NAMs are not limited to tabular data. It is possible to choose f_i to be any type of neural network architecture and consequently input whole time series instead of tabular values. There exists little work that considers the GAM framework in this more general non-tabular setting. See for example [27] for an application in genomics. When using more general networks to input whole time series x_i instead of tabular features, the “global” interpretability of f_i is lost. Because x_i is now high dimensional, it is impossible to visualise f_i as a simple function of x_i. We could, however, still read out the additive contribution of each network α_if_i(x_i) and therefore still quantify the importance of an individual time series modality to the overall classification. This analysis could either be done for a single specific input or over a whole set of inputs to gauge the overall signal modality importance.

Activation Maps (Question 2)

To alleviate the loss of the global interpretability, we chose deep fully convolutional networks (FCNs) [28] as our functions f_i. By parameterising the modality functions f_i as FCNs, we could achieve local interpretability by computing activation maps [29]. FCNs have been used extensively in time series classification [30]. They consist of several layers of intertwined operations of convolutional filters and non-linear activation functions, as well as a final output layer, which computes a global average over the convolutional filter activations across the temporal dimension.

The computed activation maps are inherently linked to the inductive bias of the FCN: The FCN detects temporally local patterns in the time series on the scale of the size of the used convolutional filters. As a result, activation maps are temporally local “explanations” that reflect the filter activations across the temporal dimension and thus allowed us to investigate which parts of the signal are discriminative for the classification (Question 2). Due to the linearity of the final global average operation, activation maps were easily integrated into the NAM framework: The additive contributions for each individual signal modality were just taken to be the average filter activation across time.

Architecture and training details

We implemented the presented NAM, as well as all downstream analysis in Python using the Pytorch [31] and the Scikit-learn [32] library. The individual subnetworks consisted of three convolutional layers. Kernel sizes were selected with nested cross-validation and ranged between 5 and 17. We used zero padding to preserve the time dimension of each input signal. Like the kernel size, the number of convolutional hidden channels was also selected with nested cross validation, which resulted in 20 hidden channels per layer. After each convolutional layer, we applied batch normalisation [33] and the rectified linear units (ReLUs) [34]. Before the global average pooling layer, network activations were linearly combined per time point to produce one single time series activation map for each signal modality.

The Adam optimiser [35] was used to train the classifier with a standard binary cross entropy loss. We trained the networks for 10 epochs with a learning rate of 0.0001 and a weight decay of 0.001. These choices resulted in optimal performance in every validation fold. The previously mentioned network hyper-parameters were also stable across validation folds. The stability of hyperparameters across subjects was a consequence of our leave-one-out approach, where training sets differ only slightly between validation folds.