Prediction of intraventricular haemorrhage in preterm infants using time series analysis of blood pressure and respiratory signals

Despite the decline in mortality rates of extremely preterm infants, intraventricular haemorrhage (IVH) remains common in survivors. The need for resuscitation and cardiorespiratory management, particularly within the first 24 hours of life, are important factors in the incidence and timing of IVH. Variability analyses of heart rate and blood pressure data has demonstrated potential approaches to predictive monitoring. In this study, we investigated the early identification of infants at a high risk of developing IVH, using time series analysis of blood pressure and respiratory data. We also explore approaches to improving model performance, such as the inclusion of multiple variables and signal pre-processing to enhance the results from detrended fluctuation analysis. Of the models we evaluated, the highest area under receiver-operator characteristic curve (5th, 95th percentile) achieved was 0.921 (0.82, 1.00) by mean diastolic blood pressure and the long-term scaling exponent of pulse interval (PI α2), exhibiting a sensitivity of >90% at a specificity of 75%. Following evaluation in a larger population, our approach may be useful in predictive monitoring to identify infants at high risk of developing IVH, offering caregivers more time to adjust intensive care treatment.

Intraventricular Haemorrhage (IVH) remains a serious threat to survival for preterm infants and neurodevelopmental outcomes 1 . Despite advances in modern neonatal care such as antenatal steroids, artificial surfactant treatment and the use of neuroprotective agents such as magnesium sulphate given to mothers in labour, rates of IVH, particularly high grade, remain unchanged. Prematurity, respiratory-distress syndrome and mechanical ventilation are among the factors that may predispose infants to IVH. Recent studies have also suggested an association between IVH and cerebral pressure passivity, that is, where changes in cerebral blood flow correspond to changes in blood pressure 2 . The need for resuscitation and cardiorespiratory management of preterm infants within the first 24 hours of life play an important role in the development and timing of IVH 3,4 , where the majority of these cases can be detected at their full extent by the end of the first postnatal week 5 . The potential to identify infants at high risk of developing IVH is thus, particularly important.
Retrospective studies of premature infants after the diagnosis of IVH have highlighted altered autonomic functions which are reflected by heart rate variability analysis 6,7 . In particular, one study showed that these differences could be detected using electrocardiogram data from the first 24 hours of life 8 . Variability of beat-to-beat systolic blood pressure and mean arterial pressure has also been shown to offer useful information in distinguishing infants who later developed IVH from those who did not 9 . Such distinctions were demonstrated using detrended fluctuation analysis (DFA), a non-linear time domain technique that is able to quantify long-range power law correlations in a given time series. Its application is characterised by a scaling exponent (α) which can be calculated over different time scales and indicates the corresponding degree of correlation 10,11 .
More recent work in this area by Fairchild et al. has demonstrated associations between a heart rate characteristic index and adverse neurodevelopmental outcomes or white matter damage 12 . Models for early prediction of IVH have explored either clinical risk factors, as in the case of Luque et al. 13 with an AUC of 0.79, or employed time series analysis techniques and physiological signals as in Tuzcu et al. who reported a sensitivity of 70% and specificity of 79% for their model using heart rate variability 8 .
The objective of this study was twofold; firstly, to explore means of improving prediction of IVH from DFA through pre-processing, and secondly, to evaluate the potential of multivariable or multimodal models in the prediction of IVH. The latter objective focused specifically on combinations of blood pressure and respiratory features and inherently involved an evaluation of robustness as applied in a clinical context. The features evaluated comprised of the mean (μ) as well as short-and long-term scaling exponents derived from DFA (α 1 and α 2 , respectively), extracted from five different time series. These were: mean arterial pressure (MAP), systolic blood pressure (SBP), diastolic blood pressure (DBP) and pulse interval (PI) series as derived from the arterial blood pressure data, as well as the interbreath intervals (IBI) from the respiratory air flow data.

Results
Study Population. The study cohort consisted of 27 low birth-weight (< 1500 g) infants, 8 of which subsequently developed IVH. We examined the differences in physiological characteristics and other metadata between the two groups, as summarised in Table 1. Although the IVH and non-IVH groups did not exhibit significant differences in the collected metadata, certain clinical characteristics, namely, the mean DBP values were observed to be significantly different for the two groups after non-linear trend removal (p < 0.05). This established a foundation for fitting the univariate logistic regression models, though evaluation of their overall performance requires reference to leave-one-out cross-validation (LOOCV) results.

Effect of Detrending.
In the initial stages of model fitting, we observed that in certain instances, the arterial blood pressure was subject to non-linear drift. We examined the effect of detrending the arterial blood pressure signal and noted that prediction performance of the fitted logistic regression models could be improved. This detrending affected the linear blood pressure features in particular, which also propagated to the beat-to-beat blood pressure values and thus the scaling exponents derived from DFA (α 1 and α 2 , the short-and long-term scaling exponents, respectively). The Mann-Whitney U-test comparisons of the non-detrended features are shown in Table 2, where the effect of detrending was characterised by the changes in AUC and p values from the two-sided Mann-Whitney U-test. For example, the AUC scores of the mean DBP model improved from 0.757 to 0.807 subsequent to detrending. A similar increase from 0.757 to 0.771 was also observed for the univariate model of the long-term scaling exponent of DBP (DBP α 2 ), motivating the inclusion of this pre-processing step for the subsequent analyses. The histograms of mean DBP and DBP α 1 are also shown in Fig. 1.

Univariate Predictors of IVH.
We fitted univariate logistic regression models using various linear and DFA features, taking the mean across all qualifying time windows of data for each subject. We evaluated the AUC, the 95% confidence interval (5th, 95th percentile) according to the Delong method for determining standard error 14 , as well as the positive likelihood ratio and threshold corresponding to a specificity of 75%. These results are summarised in Table 3. Overall, the short-term fractal exponents (α 1 ) derived from the MAP, SBP and DBP signals as well as the mean DBP yielded the highest AUCs. Respiratory variables were not found to be strongly predictive.

Multivariable Predictors of IVH.
It was also observed that model performance could also be improved through combination of the predictors extracted for the univariate models, as shown in Table 4. Predictor combinations that were significantly correlated were excluded to mitigate effects of collinearity. Many of the multivariable models exhibited higher AUC scores than univariate models, with the highest being the combination of mean DBP and PI α 2 . Similar results were obtained when evaluating all qualifying windows and with the inclusion of gestational age. Note however that the comparisons with the best performing univariate model (i.e. mean DBP, Table 3) were not statistically significant (p > 0.05) due to high degree of overlap. The ROC curves for a number of these models are displayed in Fig. 2, where the non-linear detrending process to obtain the DBP c μ feature is illustrated in Fig. 3.
Leave-One-Out Cross-Validation. The models were further evaluated using LOOCV, where the probability estimates of each testing sample were used to construct a receiver-operator characteristic (ROC) curve 15,16 , as summarised in the latter column of Table 4. Delong comparison of these LOOCV ROC curves with the corresponding LOOCV mean DBP model did not exhibit statistically-significant results (p > 0.05), with the exception of the mean DBP and SBP α 1 combination. Sensitivity Analysis. The models used to obtain these results were based on the mean feature(s) calculated from all qualifying 10 min windows for each subject, with a 30 sec overlap. To evaluate how robust these results were, we examined the effect of including all individual windows rather than the mean feature per subject, of using non-overlapping windows and of including gestational age. The first two cases involved use of a mixed-model allowing for repeated measures while the latter involved the addition of gestational age as a predictor in the existing multivariable models. In all three cases, we found that mean DBP and PI α 2 remained the most predictive combination for IVH (where AUC = 0.88, 0.86 and 0.85 for the three cases, respectively, compared to the AUC = 0.92 reported in Table 4).

Discussion
Summary of findings. This study evaluated the blood pressure and respiratory features in discerning infants with IVH from those without. The highest AUC achieved was 0.921 (95% CI 0.82, 1.00) by the model fitted with PI α 2 and mean DBP. The results of cross-validation also supported this, with an AUC LOOCV of 0.821 (0.66, 0.99). This model exhibited a sensitivity of > 90% at a specificity of 75% which is greater than that reported for the heart rate variability-based model from Tuzcu et al. (70% sensitivity, 79% specificity) 8 . This latter cohort was of a similar  size (n = 24), though it was limited to very low birthweight infants (< 1000 g) as opposed to our low birthweight cohort, potentially contributing to the difference in IVH representation observed (41.7% compared to 29.6%). Although the univariate model fitted with mean DBP exhibited an AUC of 0.807 in the initial analysis, results from LOOCV cautioned its use as a sole predictor, with an AUC of 0.607 and a non-significant 95% confidence interval of (0.38, 0.88).

Effect of detrending.
Out of all the factors we examined, detrending as part of the pre-processing phase of analysis resulted in the greatest improvement in prediction of IVH. It rendered both mean DBP and DBP α 2 significantly different (p < 0.05) between the IVH and non-IVH groups, where the AUC scores for two univariate models fitted with these features increased from 0.757 to 0.807 and 0.771, respectively (Tables 2 and 3). These improvements suggest that long-term drift and/or baseline wander of the blood pressure signals, among other aspects of signal quality, may confound the results from DFA as well as linear parameters. The significant contribution of mean DBP is of particular importance, given previous assertions that linear features alone were not an informative proxy of cerebral perfusion pressure 17 and reports that such features were not significantly different between the IVH and non-IVH groups 9 . This relationship to DBP may also have been a reflection of a widened pulse pressure, seen in symptomatic patent ductus arteriosus. It is necessary to note however that this feature was not explicitly evaluated in the previous studies. Aside from the detrending, this analysis differed in other aspects including the quality control measures for window selection and the evaluation process; the features extracted represent the mean of all qualifying 10 min windows across the recording, rather than a single segment. In the clinical application of DFA to monitored signals, we strongly recommend examination of signal data to determine whether overall detrending is necessary prior to analysis.
Use of bivariate models. Subsequent to detrending, a further increase in AUC achieved through fitting of bivariate models, where the combination of mean DBP with MAP α 1 and PI α 2 , for example, exhibited respective scores of 0.871 and 0.921. This would suggest that relevant, non-redundant information may be captured using linear and DFA-based approaches in the prediction of IVH. The short-term scaling exponents for the beat-to-beat MAP and SBP, along with mean DBP were shown to be relatively robust markers in prediction of IVH for the studied cohort. Although all of the studied infants had triggered ventilation modes (synchronised intermittent positive pressure ventilation), breathing termination was not employed and so the potential for adverse patient ventilator interaction was possible. Previous work examining the impact of patient-ventilator asynchrony indicates the significant potential for IVH with this phenomenon 18,19 . These results align with those previously reported 9 and the altered vagal nerve activity in infants with IVH 20 . From further evaluation of the bivariate models, it was clear that the initial analysis did not necessarily translate to robust and consistent performance in leave-one-out cross-validation. The model with the highest AUC in the initial analysis achieved an AUC LOOCV of 0.821 (0.66, 0.99), though it did not exhibit a statistically significant improvement on the univariate  reference model (mean DBP) 14 . It was interesting to note the inclusion of pulse interval-based features in the highest-scoring model in both the initial analysis and cross-validation, given its use as an estimate of heart rate variability and the reported high correlation between the two 21,22 . The accuracy of this estimation however, has not been clarified, particularly in the neonatal context, though electrocardiogram-based heart rate variability has been found to offer useful information in distinguishing infants with and without IVH 8 .
Addition of respiratory signals. In this study, we found that the addition of respiratory signals did not considerably improve model performance. The fractal dynamics of respiration have been applied in the context of preterm infants 23,24 , though not with respect to IVH. As for the models fitted with interbreath interval (IBI) features, mechanical ventilation may have contributed to their observed lack of prediction capacity (p > 0.05), despite the relevance of respiration mechanics in the development of IVH 18 . It is also possible that the IBI-based features were not suited to characterising patient-ventilator asynchrony.
Clinical significance and application. Hypercarbia, high ventilator pressure and patency of the ductus arteriosus are among the factors and events that may influence the fluctuation of blood pressure of preterm infants in the neonatal intensive care unit 25 . Infants who later developed IVH exhibited lower mean DBP and a higher DBP α 2 (p < 0.05) across the entire recording in this study. Recent studies have reported a range of observations pertaining to blood pressure and IVH, with the main focus on characterising cerebral perfusion. These include reports of IVH being associated with the elevated diastolic closing margin 17 and significant deviation above a defined optimal MAP value in infants who later developed IVH 26 . This approach may be applied to a clinical context in a manner similar to that shown in Fig. 4, offering examples of both correct and incorrect classification of IVH from the studied cohort. A threshold may be defined according to the dashed line in each of the cases (a) to (d), where calculated probabilities exceeding this threshold  Table 4. Multivariable logistic regression models. Features included mean (μ), short-and long-term scaling exponents (α 1 and α 2 , respectively) for mean arterial (MAP), systolic (SBP) and diastolic (DBP) blood pressure, as well as pulse interval (PI) time series. LR denotes the positive likelihood ratios, the 95% confidence intervals (CI) are reported for the AUC. p values are derived from the Delong comparison 14 with the non-detrended mean MAP model. The corresponding AUCs from leave-one-out cross-validation (AUC LOOCV ) are also reported, where *denotes a statistically-significant (p < 0.05) difference from the Delong comparison of the LOOCV mean DBP model 14 . could flag infants at high-risk of developing IVH. Further model evaluation requires validation on a larger and more balanced cohort to estimate the prediction error and support its potential application in a clinical context.

Limitations.
We acknowledge that this study was limited by the size of the dataset (n = 27) as well as representation of IVH (29.6%), slightly lower than the referenced 35-45% of incidence reported in neonatal care facilities 27 . Model evaluation was also limited by the low number of IVH cases (n = 8), though our LOOCV and sensitivity analyses showed the main findings to be consistent. Another limiting factor was the signal quality of the recordings which was managed by implementing quality control measures as part of the feature extraction process.

Conclusion
In conclusion, this study found mean DBP and short-term scaling exponents from beat-to-beat MAP, DBP and SBP to be useful markers in the prediction of IVH in preterm infants. Non-linear trend removal and the inclusion of additional features such as the short-term scaling exponent (α 1 ) of MAP was able to improve model performance. Of the models evaluated, the one that performed consistently in both the initial analysis and cross-validation was fitted with mean DBP and PI α 2 . In a clinical context, such an approach to signal processing and predictive monitoring could be applied, where a running 10 min window could continuously evaluate the relevant features from qualifying segments of data. Following evaluation in a larger population, these features may be helpful in identifying infants at high-risk of developing IVH, offering caregivers more time to adjust intensive care treatment.

Methods
Data Collection. Physiological data was collected from the infants within 1-3 hours of birth as part of a prospective clinical investigation at a large tertiary neonatal intensive care unit in Sydney, Australia. The study was approved by the Sydney West Area Health Service Human Research and Ethics and conducted according to the World Medical Association Declaration of Helsinki. Informed parental consent was obtained in all cases. Inclusion criteria for the cohort comprised low birthweight (< 1500 g), gestational age (< 30 weeks) and an absence of significant congenital anomalies. Of the 46 infants enrolled, 27 infants had arterial blood pressure and air flow wave recordings with sufficiently long, artefact-free segments. The average (SD) length of recording was 156 (34) mins. Intra-arterial blood pressure was measured using an umbilical or peripheral arterial catheter, following single-point calibration to atmospheric pressure, collected using a bedside patient monitor (Philips Agilent Systems, Philip Healthcare, North Ryde, Australia), while the raw air flow wave was acquired from a ventilator (Babylog 8000, Drägerwerk, Lübeck, Germany). Both signals were sampled at 1 kHz and recorded by an analog data acquisition system (ADInstruments, Sydney, Australia). Cranial ultrasounds were performed at 2, 12, 24 and 36 hours then daily for the first week. The presence and grade of IVH was determined according to the Papile system 27 .
Signal Processing and Data Analysis. Signal processing and feature extraction was completed in Python (Python Software Foundation, version 2.7. https://www.python.org/). Each of the arterial blood pressure and air flow signals were down-sampled to 125 Hz prior to analysis for computational efficiency. This frequency was sufficient for peak detection in both respiratory and blood pressure signals. From the downsampled signals, the following time series were extracted; the beat-to-beat MAP, SBP, DBP and PI, as derived from arterial blood pressure, as well as IBIs derived from air flow data. The signal quality constraints of the air flow data limited extraction of other respiratory features such as peak flow.
Only the arterial blood pressure signal was found to exhibit significant drift, defined by non-linear trends in the diastolic and systolic blood pressure ranges. Thus, the detrending was applied solely to this signal, as shown in Fig. 3. Such a correction would also minimally impact the derivation of the IBI-based features from the air flow signal. The overall trend of each signal was determined using a large-window median filter (window width = 1000 ms) on a further downsampled signal and the mean-centred trend was subsequently removed from the original arterial blood pressure signal. This approach was similar to the baseline wander removal that has been applied widely prior to feature extraction from the electrocardiogram signal 28 . An example of this detrending process is shown in Fig. 3.
The features used in IVH prediction were extracted from a running 10 min window of arterial blood pressure and air flow data, shifted in 30 sec increments across the total recording length. This approach was adopted to simulate the application of these techniques in a clinical setting, where windows which fulfilled the quality criteria were included for feature extraction. This criteria comprised defined ranges for the allowable number of beats and breaths in a given window (40-250 beats per minute and > 20 breaths per minute), a maximum limit for an absence of detected beats (15 sec) as well as an absence of large spikes in the arterial blood pressure signal (range . Arterial blood pressure data and the predicted probability for IVH using the highest scoring model, mean DBP and PI α 2 for correct classification of (a) IVH and (b) non-IVH, as well as incorrect classification of (c) IVH and (d) non-IVH. The threshold for classifying IVH, designated by the dashed line was defined at 90% specificity and 85% sensitivity. Red and blue markers represent windows that exceeded and did not exceed the threshold, respectively.
Scientific RepoRts | 7:46538 | DOI: 10.1038/srep46538 of beat-to-beat SBP < 30 mmHg). For each respective time series, outliers were removed by imposing a maximum change of 150% from the previous data point and also a maximum loss of 30% for each window. Features including the mean (μ), short-and long-term scaling exponents (α 1 , α 2 , respectively) from DFA of the five time series (MAP, SBP, DBP, PI and IBI) were subsequently extracted.
Developed by Peng and co-workers 10 , DFA is able to quantify long-range power law correlations and accommodate for confounding non-stationarities often seen in real-world signals. It does this through the detrending, that is, linear trend removal, step prior to calculating the root-mean squared fluctuation as defined in equation 1. where y(k) is any given time series, y n (k) the local linear trend for a given segment, and N the number of data points in the series for a given round of analysis. The application of DFA is further explained by Thamrin et al. 11 .
The scaling exponent α is calculated from the gradient of the line fitted to the Log-Log relationship between n and F(n). In this case, the short-term scaling exponent was defined over 4-15 beats, as aligned with similar observations of this data 9 and similarly defined for heart rate variability analysis of preterm infants 8 . The long-term scaling exponent was determined across 15-50 beats.
Model Fitting and Evaluation. Statistical analysis was completed using R 3.3.1 software 29 . Logistic regression models were used to fit the mean of extracted predictors across all qualifying 10 min windows, while the AUC was used to assess accuracy in predicting IVH. Fitted models were evaluated using leave-one-out cross-validation, where the predicted probability of each test sample was subsequently compiled and used to generate a ROC curve for performance comparison.