The institutional review board of the Charité – Universitätsmedizin Berlin approved the study (EA2/104/21). Due to the retrospective nature of the study, the need for patient consent was waived.

Continuously and routinely monitored vital parameters were extracted from the electronic health systems and harmonized for further analysis.

Data were exported for all patients admitted to the clinic between January 1, 2013, and February 1, 2022. The data set was then refined on the basis of the following inclusion criteria, and all calculations were later performed within the refined data set. To include all patients with continued clinical suspicions of withdrawal symptoms but to exclude those with one-time-only suspicions or accidental documentation, we performed the export by selecting patients with at least 3 documented FSs, since a pharmaceutical intervention was usually not initiated on the basis of a single scoring result. We cross-referenced this export selection with patients classified as having NAS in accordance with ICD-10 (International Statistical Classification of Diseases, Tenth Revision) criteria for quality control purposes. We categorized patients into subgroups of primary NAS and iNAS based on opioid medication, history of surgery, and time after birth before documenting the first FS for each patient.

Primary NAS was coded when at least 1 FS was documented before any opioid medication was administered, any surgery was performed, and the patient had not yet approached postpartum day 8. iNAS was coded when any opioid medication was administered before the patient’s first FS regardless of postnatal age. Patients with documented FS who did not meet any of these criteria were excluded. To design a sensitivity analysis, the analytic code was also applied to total study population without exclusion due to unclear NAS classification.

Each variable was checked for availability within the hospital information system (SAP/Cerner) as well as the PDMS (COPRA) used in the neonatal wards (levels 1‐3).

Medications were not named consistently; hence, their names had to be preprocessed manually. We exported all unique medication entries from the PDMS and categorized them manually. The complete list of medication categories is provided in Multimedia Appendix 1.

To evaluate patients’ demographics, we recorded their sex, gestational age at birth, birth weight, mode of delivery, number of documented FSs, and whether a time frame for individual baseline calculation was available and, if so, whether data were available within this time frame for the abovementioned vital parameters (including heart rate, respiratory rate, peripheral oxygen saturation, and mean blood pressure).

We calculated means for all variables listed below within specified time periods (t–1 and t–2, see the Time Periods section): heart rate, respiratory rate, peripheral oxygen saturation, and blood pressure.

Additionally, we generated an individual baseline for each patient by calculating the mean for each variable in a period of up to 5 days before documenting the first FS, which we defined as the relevant beginning of withdrawal (Figure 1). When calculating the individual baseline, we excluded spacer periods immediately post partum to minimize effects from postnatal adaptation and those immediately before documenting the first FS to minimize the effects of early-onset withdrawal. We set this spacer period to 1 day for patients with primary NAS and 3 days for those with iNAS (Figure 2). The sensitivity analysis was carried out with both spacer periods for the whole collective. We introduced the difference between this individual baseline and the mean of the respective vital parameter within the specified period before documenting any FS as a new variable to use as an alternative to the mean of the vital parameters and henceforth referred to this variable as the “baseline-controlled mean” (Figures 1 and 2).

Body temperature was assessed in variable patterns, and we calculated the mean and baseline-controlled mean body temperature within 1 day before documenting each FS and used these data for both time periods (t–1 and t–2).

Furthermore, we included the following variables. (1) The pharmacodynamics of buprenorphine—the first-line pharmacotherapy for NAS—is highly complex, and data on transferability to neonates are limited [19,20]. Hence, we did not attempt to estimate pharmacodynamics in neonates but considered the time between the last documented opioid medication and each FS instead. As these hours since medication can only be recorded for patients who have been administered any opioids before documenting the respective FS, either because of iNAS or because of treatment of any type of NAS, we used the date of birth as the date of the last opioid medication for infants with primary NAS if no opioid medication was documented more recently. (2) The last body weight measure before documenting each FS was considered as a percentage of the individual’s birth weight. (3) The current gestational age at each FS documentation was recorded.

All graphs and models were created for 2 time periods. With the goal of exploring options for predicting withdrawal symptoms, we focused on time periods before each FS. The first time period, t–1, was set to 1 hour±10 minutes before each FS, resulting in a 20-minute period from 10 minutes before until 10 minutes after the time point of 1 hour prior of each FS documentation. The second time period was set in the same manner to 2 hours±10 minutes before the respective FS, resulting in an earlier time period, t–2. The time periods are visualized in Figure 2.

Data analysis was carried out using RStudio (version 2022.07.1+554) and R (version 4.2.1; 2022-06-23 ucrt) [21,22] using the following packages and their dependencies in addition to the function included in R, RStudio, and the ::base-package R during data extraction and harmonization: cli [23], data.table [24], dplyr [25], lubridate [26], tibble [27], and tidyverse [28].

We used the consort package to generate Figure 3 [29]. Table 1 was created using the tableone package; significance was tested using chi-square tests for categorical variables, Wilcoxon tests for skewed variables and t tests for normally distributed metric variables [30]. Skewness was assessed using the summary-function from tableone in accordance with tableone documentation [30]. We generated graphs with ggplot2 [31] and fitted our mixed effects models using lme4 [32]. Goodness of fit parameters of the mixed effects models was calculated using MuMIn [33]. Table output from RStudio was facilitated using flextable [34]. All code is has been published previously [35].

We visualized data availability in a clustered bar plot, reporting the number of data points per variable available for each patient within each period (Figure 4). By computing 2 statistical measures for each of the 4 vital parameter variables in each of the 2 time periods, we obtained a total of 16 graphs (see Multimedia Appendix 1). To visualize data distribution, we generated heat map plots. The color of each tile is set by the number of data points weighted by the number of measurements generating the data point so that data points based on more observations contribute more to the color scale. Hence, 2 data points based on 1 observation each and 1 data point based on 2 observations both result in the same color of the respective tile.

We developed several mixed effects models. The models were fitted to descriptively analyze the relationship between the FS and vital parameters as well as NAS type, the time elapsed since the last opioid medication, gestational age, and the percentage of birth weight reached by the most recent body weight measurement. We controlled for interindividual differences by including the patient identifier as the random effect. We generally selected one set of vital parameter variables mixing neither means or individual baseline-controlled means nor the time periods of different vital parameters within a single mixed effects model. The models were not fitted to predict the FS; therefore, neither cross-validation nor bootstrapping were applicable. The regression equation for the full model was as follows, where “mean” could be substituted with “baseline-controlled mean” and “t–1” with “t–2” throughout the equation:

We evaluated the goodness of fit again for models with a simplified set of variables, excluding variables with small effects (estimates of <0.05) and high degree of missingness (>30% missing values). Goodness of fit was determined using the Akaike information criteria (AIC), Bayesian information criteria (BIC), and R²_m and R²_c and is listed in Multimedia Appendix 1.

For the sensitivity analysis, we excluded the NAS-type variable from the full model. The model’s results and goodness-of-fit data from the sensitivity analysis can be found in Multimedia Appendix 2.

Patient demographic data are provided in Table 1, and patient allocation is demonstrated in Figure 3. Out of 205 neonates with NAS, 78 (38%) had been ICD-10–coded. Abstinence after intrauterine exposure (P96.1) was coded in 17.4% (n=4) of patients who were classified as having primary NAS and 5.5% (n=10) of patients who were classified as having iNAS. Abstinence after therapeutic exposition (P96.2) was coded in 21.7% (n=5) and 34.6% (n=63) of the respective patient groups. We obtained 7050 FS data points and calculated the baseline-controlled mean for up to 166 (81%) infants (Table 1). Two reasons prevented us from doing so in the other cases: no definable time period (as illustrated in Figure 2) or no data points within the time period to calculate upon.

As shown in Figure 4, heart rate measurements showed the highest data density and mean blood pressure measurements showed the lowest. Heart rate measurements were the only variable for which at least 5 data points within a time period were commonly available, the only other variable being mean arterial blood pressure in very rare instances.

The plots show a widely distributed pattern for mean and baseline-controlled mean heart rate (Figure 5) and mean arterial blood pressure (Figure 6), each plotted against the FS. The tile color and different ranges shown in the color legends illustrate differences in the density of data available for heart rate and mean arterial blood pressure measurements; this corresponds to the data density described above and is shown in Figure 4. While there was no direct relationship between the FS and the respective parameters visible, the graphs for the baseline-controlled version of each parameter showed a narrower spectrum on the x-axis. In particular, the baseline-controlled blood pressure shows a discernible trend of greater-than-zero values, indicating a rise in blood pressure in neonates in comparison with that before withdrawal assessment. However, this rise in baseline-controlled blood pressure does not clearly increase with an increase in the FS.

Estimates for vital parameters stayed either positive or negative across all fitted models and varied only in the corresponding SEs, CIs, and P values. Among all models, the model containing all variables, using the later time period of t–1 and the baseline-controlled mean as the statistical measure, yielded the highest R²_c and R²_m and second-lowest AIC and BIC (Table 2). This model was fitted on 357 observations obtained from 84 infants.

An increasing individual baseline-controlled heart rate (estimate 0.04, 95% CI 0.02 to 0.07) and an increasing arterial blood pressure (estimate 0.05, 95% CI 0.01-0.08) correlated significantly with an increased FS.

Furthermore, a decreasing individual baseline-controlled peripheral oxygen saturation (estimate –0.09, 95% CI −0.19 to 0.01), decreasing gestational age (estimate −0.06, 95% CI −0.18 to 0.06) as well as increasing individual baseline-controlled respiratory rate (estimate 0.00, 95% CI −0.02 to 0.03), increasing hours since the last medication (estimate 0.00, 95% CI 0.00-0.00), increasing baseline-controlled body temperature (estimate 0.33, 95% CI −0.63 to 1.28), increasing percentage of birth weight (estimate 0.00, 95% CI −0.01 to 0.00), and the status of iatrogenic (vs primary) NAS (estimate −0.27, 95% CI −3.41 to 2.92) were associated with an increasing FS. However, while the inclusion of these variables improved the goodness of fit of the model, the estimates showed SEs and 95% CIs too large to be considered significant.

The irregularity and scarcity of blood pressure measurements in neonatal standard care resulted in missing values for the individual baseline-controlled mean arterial blood pressure and thereby reduced the amount of complete data on which the model could be fitted. Due to the large number of missing values (mean blood pressure unavailable for >90% of FS), imputation was not applicable. When blood pressure was not included—thereby reducing the complexity of the model but increasing the number of complete observations—all measures for goodness of fit decreased (worst model with blood pressure [R²_m=0.08, AIC=2345]; best model without blood pressure [R²_m=0.07, AIC=25,992]). On excluding the NAS type from our sensitivity analysis, these results were confirmed. Including previously excluded patients and using spacer periods of the same length for all patients during baseline calculation resulted in similar results (Multimedia Appendix 2).

The abovementioned estimated effects are small and may seem clinically irrelevant. However, due to the nomenclature of mixed effects models, the estimates refer to 1-unit changes of the respective variable. This implies that a heart rate increase of 10 beats per minute from the individual baseline would coincide with an FS increase of 0.4 and a 10-mm Hg increase (increase of 0.6) in the mean arterial blood pressure assuming that all other values remained constant.