To facilitate comprehensive data collection and analysis for all enrolled infants, we utilized the Integrated Data Repository (IDR) in collaboration with the University of Florida Clinical and Translational Science Institute (UF CTSI). The IDR aggregates clinical data from the electronic health record (EHR), incorporating both structured data (e.g., demographics, clinical variables, laboratory results) and unstructured data (e.g., clinical notes, bedside physiologic monitoring). All retrieved data underwent systematic cleaning and preprocessing using R Studio to ensure data integrity and consistency prior to analysis. This thorough data management process allowed the integration of the physiologic variables into machine learning models, creating a complete and accurate physiologic profile for each infant.

In addition to continuous monitoring with the cardiorespiratory physiologic monitor, regional cerebral oxygenation was measured using INVOS™ (Medtronic, formerly Covidien, Minneapolis, MN) Cerebral/Somatic Oximetry Infant-Neonatal Sensors. These sensors utilize near-infrared spectroscopy (NIRS) to monitor regional mixed venous-arterial saturations. Sensors were placed on either the left or right side of the neonate's forehead, and data were continuously collected throughout the hypothermia phase, during the rewarming period, and immediately post-rewarming. This comprehensive monitoring approach provided a detailed dataset for analysis, capturing both systemic and cerebral physiologic parameters.

MRI were performed at either 3–5 (n = 18) days of age following rewarming or 7–12 days of age (n = 7, not stable for imaging at 3–5 days of age). One subject had an MRI performed at day 1 of life and another subject at day 20. Neonates were imaged on a Siemens Magnetom Verio 3 T scanner (Siemens, Malvern, PA) at UF Health Gainesville. A single blinded subspecialty board-certified neuroradiologists with over 10 years of experience in neonatal imaging interpreted all the MRI images using the Weeke scoring system (21). The Weeke scoring system evaluates brain injury across three regions: deep grey matter, white matter/cortex, and the cerebellum, with an additional subscore assessing the presence of intraventricular hemorrhage (IVH), subdural hemorrhage (SDH), and cerebral sinovenous thrombosis (CSVT). Each anatomical region is systematically scored based on the extent and distribution of injury. The deep grey matter subscore, with a maximum of 23 points, assesses the thalamus, basal ganglia, posterior limb of the internal capsule (PLIC), brainstem, perirolandic cortex, and hippocampus. Injury in these areas is scored as 0 (no injury), 1 (focal injury affecting <50%), or 2 (extensive injury affecting ≥50%). Additionally, injuries are noted as either unilateral (score of 1) or bilateral (score of 2). The white matter/cortex subscore, with a maximum of 21 points, evaluates damage to the cerebral cortex, cerebral white matter, optic radiations, corpus callosum, punctate white matter lesions (PWML), and parenchymal hemorrhage. Injury is scored as 0 (no injury), 1 (focal injury involving one lobe), or 2 (extensive injury affecting multiple lobes). The cerebellum subscore, with a maximum of 8 points, examines lesions or hemorrhages within the cerebellum. A score of 0 is given for no injury, 1 for focal lesions smaller than 0.5 cm, and 2 for extensive lesions larger than 0.5 cm or multiple lesions. In addition to these anatomical regions, the Weeke score incorporates an “Additional” subscore assessing the presence of IVH, SDH, and CSVT. For each of these conditions, the scoring is binary: 0 for absence and 1 for presence, contributing to a maximum additional score of 3. The total score is calculated by summing the scores from the grey matter, white matter/cortex, cerebellum, and additional categories, with a maximum score of 55. If 1H-MRS data are available, abnormalities in the basal ganglia and thalamus, such as reduced N-acetyl aspartate (NAA) or elevated lactate peaks, are incorporated into the grey matter subscore, increasing the total score to a maximum of 57.

The dataset consisted of physiological biomarkers from 138 subjects monitored over 72 h during therapeutic hypothermia. The recorded biomarkers included heart rate (HR), systemic oxygen saturation (SpO₂), cerebral oximetry (rcSO2), arterial systolic blood pressure (ArtSBP), arterial diastolic blood pressure (ArtDBP), and mean arterial pressure (MAP). Some measurements were missing throughout the data collection process. To address these gaps, missing values were imputed using the average of the five nearest neighboring timestamps. Each biomarker feature was standardized using Z-score normalization, transforming the data into a standard normal distribution with a mean of zero and a standard deviation of one.

Our dataset had small, randomly distributed gaps due to the retrospective nature of this study and the inclusion of out-born infants, which resulted in occasional missing measurements. Most abnormal physiological events were concentrated early in the monitoring period, whereas the majority of the data reflected stable physiological patterns. We specifically chose the five-nearest-neighbor averaging imputation method due to the following advantages: (1) Local Adaptability: Unlike simpler methods such as forward-fill or global linear interpolation, averaging the five closest temporal neighbors (including timestamps before and after a missing point) preserves local temporal context and prevents oversmoothing; (2) Outlier Resilience: Utilizing multiple neighboring points reduces the likelihood of propagating single-point measurement errors or outliers, a risk inherent in methods like forward-fill; (3) Trend Preservation: Local averaging is capable of accurately capturing subtle, short-term trends and variations in the signal, maintaining physiological realism more effectively than linear interpolation, especially when abrupt physiological changes occur.

We carefully considered alternative methods, including multiple imputation (MI), model-based imputation techniques (such as mixed-effects models), random forest-based approaches, and deep learning methods and evaluated the pros and cons. Multiple Imputation and Model-based Approaches: Methods like multivariate linear mixed-effects models are robust for multilevel or longitudinal data, leveraging hierarchical structures to mitigate bias. However, these methods require assumptions about data distributions and sufficient variability in hierarchical structures, which may not align well with our minute-level physiological data. Random Forest Imputation: Although flexible and capable of capturing nonlinear relationships, random forest-based methods may introduce biases when handling skewed physiological variables or complex interactions inherent to clinical biomarkers. Deep Learning Methods: These models excel with extensive datasets, where substantial training data are available. Given our limited dataset size and small gap intervals, deep learning models are not feasible, as they typically require large training samples and may not reliably capture the short-duration gaps accurately.

After evaluating these methods' strengths and limitations relative to our dataset's characteristics and clinical context, we concluded that the five-nearest-neighbor averaging approach is most suitable for accurately and reliably filling small gaps in our time-sensitive physiological data (22).

We built a Long-Short Term Memory (LSTM) model to predict no-to-mild vs. moderate-to-severe brain injury on MRI, based on the time sequence biomarker features. LSTM is a specialized type of recurrent neural network (RNN) designed to overcome the limitations of classic RNNs, particularly the vanishing gradient problem, by preserving information over long temporal spans through its gated memory cell architecture. In this study, we chose LSTM over the alternative architectures such as Gated Recurrent Unit (GRU) due to its enhanced ability to capture long-range temporal dependencies, which are essential for early diagnosis. Although GRU offers faster training and a simpler structure, LSTM provides finer control over memory retention—critical in our application, where diagnostic accuracy depends on the length of observation (e.g., 24 h vs. 48 h) (23).

There are three major components in the LSTM model: (1) the forget gate (Ft) which will decide whether we should keep the information from the previous timestamp or forget it. Then, a sigmoid function is applied over it. That will make ft a number between 0 and 1; (2) the input gate (It) which will decide which of the values from the inputs is to be used to change the memory. The sigmoid function determines whether or not to allow 0 or 1 values through. In addition, using the tanh function, you can assign weights to the data, determining their importance on a scale of −1 to 1, and 3) the output gate (Ot) which will generate final output based on the input and memory of the block. When the sigmoid function is used, it determines whether the 0 or 1 value should be allowed through. Given h hidden units, and the batch size is n, and the dimension of inputs is d. We will introduce the function of each gate in the following:Ft=σ(XtWxf+Ht−1Whf+bf)It=σ(XtWxi+Ht−1Whi+bi)Ot=σ(XtWxo+Ht−1Who+bo)Where,

Xt∈Rn×d: input to the current timestamp t.

Wxf,Wxi,Wxo∈Rd×h: weight associated with the input.

Ht−1∈Rn×h: the hidden state of the previous timestamp t−1.

Whf,Whi,Who∈Rh×h: the weight matrix associated with hidden state.

bf,bi,bo∈R1×h: the bias.

The LSTM model consists of a single LSTM layer with 64 hidden units for temporal feature extraction, followed by fully connected layers serving as the classifier. Given the limited cohort size relative to the model complexity and training duration, there is an inherent risk of overfitting, where the model may memorize training data instead of learning generalizable patterns. To mitigate this, we employed several overfitting prevention strategies. First, we reduced model complexity by using only one LSTM layer in feature extraction layer. Second, we introduced L2 regularization (weight decay = 1e-2) in the optimizer to penalize large weights. Third, we incorporated early stopping via validation monitoring with patience of 10 epochs. We split the dataset into training, validation, and test sets with a 70%/15%/15% stratified sampling strategy within each class, ensuring balanced representation across sets. Given that the duration of the time-series input plays a critical role in diagnostic performance, we further explored the trade-off between input length and predictive accuracy. Specifically, we compiled two datasets using 24-h and 48-h observation windows, and trained LSTM models on each to compare performance. Models were trained for up to 1,500 epochs using stochastic gradient descent (SGD) with a learning rate of 0.01 and momentum of 0.8. Each setup was trained 10 times and the average accuracy and standard deviation on the test were reported to ensure a robust evaluation. The implementation was developed in PyTorch and executed on a system running Ubuntu 22.04 with an NVIDIA RTX 4080 GPU.

A total of 138 subjects were included in this cohort analysis (39% outborn). Seventy-three neonates (52.9%) exhibited moderate to severe injury on MRI, while 65 neonates (47.1%) had no or mild injury. Both groups were similar in several key characteristics. The majority of neonates were male (60% in the moderate/severe group vs. 54% in the no/mild group). Gestational age was comparable between the groups (38.2 ± 2 weeks for both), as was birth weight (3,212 ± 740 g for the moderate/severe group vs. 3,219 ± 752 g for the no/mild group). Additionally, the incidence of cesarean delivery was similar (56% in the moderate/severe group vs. 65% in the no/mild group), as were umbilical cord gas pH and base deficit values.

Notably, a sentinel event was reported in 43% of the neonates in the no/mild group compared to 33% in the moderate/severe group. One-, five-, and ten-minute Apgar scores did not differ significantly between the groups. However, neonates in the moderate/severe group had a higher incidence of seizures (49% vs. 14%) and a greater proportion with Sarnat scores of III on the initial neurologic examination (37% vs. 8%) compared to those in the no/mild group (p < 0.05). For a detailed description of the neonatal characteristics, please refer to Table 1.

Neonates in the moderate/severe group required respiratory support during hypothermia in 96% of cases, including mechanical ventilation in 84% and CPAP in 11%. In comparison, 85% of neonates in the no/mild group required respiratory support, with 68% receiving mechanical ventilation and 8% on CPAP (p < 0.05). Similarly, neonates in the moderate/severe group required ionotropic support during hypothermia in 60% of cases, compared to 34% in the no/mild group (p < 0.05).

The median time to MRI was 4 days (IQR: 4–6) for the moderate/severe group and 4 days (IQR: 4–6) for the no/mild group. The median length of hospital stay was 15 days for both groups, with an IQR of 10–27 days for the moderate/severe group and 9–32 days for the no/mild group (Table 1).