This was a two-center, prospective, observational cohort study conducted in academic tertiary and quaternary level neonatal intensive care units (NICU) at Prentice Women’s Hospital and Lurie Children’s Hospital in Chicago, Illinois. Infants born <32 weeks’ gestation between September 1, 2022, and December 31, 2023, were screened. Infants with chromosomal aneuploidy and critical congenital heart disease (apart from PDA and PFO) were excluded from the study. Data pertaining to the patients’ baseline characteristics, mode and duration of respiratory support, and morbidities associated with prematurity were recorded prospectively until hospital discharge.

The primary outcome was the post-menstrual age (PMA) in weeks required to achieve liberation of respiratory support (LRS), defined as room air or <1 L/min (FiO2 = 100%). Secondary clinical factors of interest were tracheostomy placement or death prior to 6-months of age, moderate or severe bronchopulmonary dysplasia, treatment of retinopathy of prematurity (laser or anti-vascular endothelial growth factor (VEGF) injections), grade III-IV intraventricular hemorrhage, and pulmonary hypertension at 36 weeks’ PMA as defined by echocardiogram. Echocardiographic determination of pulmonary hypertension was per the attending cardiologist’s discretion and graded as mild, moderate, or severe. BPD was defined as grade 1–3 [15].

Lung ultrasounds were completed weekly beginning at 2 weeks of age by three trained sonographers (CM, RW, SM). Scans were completed with a commercial transportable ultrasound device (SonoSite SII) using a high frequency linear probe (13–6 MHz). The scanning protocol included five lung zones in each lung (anterior superior, anterior inferior, lateral, posterior superior, and posterior inferior: as shown in Fig. 1 in Zong et al.) [16]. Anterior and lateral zones were assessed in the supine position, while posterior zones were assessed in the lateral decubitus position. Two, 3–5 s cine clips were taken of each lung zone and stored for later analyses. The patient’s position was recorded (supine, prone or lateral decubitus) one hour prior to the scan. Aeration in each lung zone was scored from 0 to 3 points (total score ranging from 0 to 30), using a validated extended lung ultrasound scoring protocol [17]. The sum of the total score serves as an index for the severity of lung disease. The LUSSc was assigned as follows: 0 indicates A-pattern (defined by the presence of only A lines); 1, B-pattern (defined by the presence of ≥ 3 well-spaced B-lines within the frame); 2, severe B pattern (defined as the presence of coalescent B-lines with or without consolidation < 1 cm in depth from the pleural line); and 3, coalescent B-lines with extended sub-pleural consolidation (Fig. 1). Lung ultrasound findings have been shown to correlate closely with extravascular lung water measured by transpulmonary thermodilution [18]. Per multiple prior studies, A-lines represent reflections of the pleural line that are visible when ultrasound waves diffuse through a non-obstructed air-filled lung. B-lines occur when there is increased lung density filling the interstitium and/or alveolar space [16,17,18,19,20]. 

A 0 points; normally aerated lung with A line pattern (A-line denoted by arrows) that represents reflections of sonographic waves off the pleura in a non-obstructed air-filled lung. B 1 point; mild interstitial syndrome with non-confluent B-line pattern (B-line denoted by arrow heads) representing transmission of sonographic waves through lung parenchyma with increased density filling the pulmonary interstitial and/or alveolar spaces. C 2 points; severe interstitial syndrome with confluent B-line pattern with or without evidence of consolidation < 1 cm in depth from the pleural line. D 3 points; extensive lung consolidation (Denoted by asterisk) with coalescent B-lines with sub-pleural consolidation ≥ 1 cm in depth from the pleural line extending into the parenchyma.

Lungs scans were scored by two physician reviewers credentialed in point of care ultrasound acquisition and interpretation. One physician reviewer who did not scan these infants was blinded (RH) to clinical course/outcomes; the other reviewer scanned infants, and thus, was aware of infants’ clinical trajectory (SM). The interrater intraclass correlation coefficient was quantified.

Description of LUSSc by PMA was the first goal to quantify how these sonographic markers of pulmonary pathology evolve during the NICU hospitalization. Given the hypothesized potential difference in quantifying the LUSSc, we reported the trajectory of LUSSc as modeled by our two reviewers of each scan. The interrater intraclass correlation coefficient was calculated to compare the agreement between the 2 reviewers’ scores.

Then, to understand the presence and/or magnitude of the relationship between LUSSc and age to achieve LRS, we utilized the maximum LUSSc obtained between 30-32 weeks’ PMA, a full 4 weeks before the diagnosis of BPD is established. We chose this time point after reviewing our descriptive data reported in Fig. 4. This maximum score was used as the main exposure on the outcome of age-at LRS (in weeks PMA). Candidate confounding variables were incorporated as those associated with the exposure and outcome and/or those hypothesized to be related to the severity of BPD (e.g., SGA < 10th centile, sex, position of infant during the scan, and respiratory support during scan).

We conducted quantile regression to examine the association between [exposure variable, e.g., total lung ultrasound score (LUSSc)] and the outcome variable ([Time to LRS, expressed in PMA]) in unadjusted association. Quantile regression estimates conditional quantiles (e.g., the median) of the outcome distribution rather than the mean. We used this approach because of our small sample size and because the outcome distribution of time to LRS is non-parametrically distributed and skewed. We estimated median regression models (τ = 0.5) and reported regression coefficients with 95% confidence intervals. Our limited sample size precluded estimations of quantiles other than the reported median (e.g., 25th%ile, 75th%ile). Of note, gestational age was not considered in these models as it is embedded in the calculation of age (PMA).

In secondary analyses, to describe the determinants of LUSSc and their inpatient trajectories, we the collected clinical variables (e.g., sex, gestational age, SGA<10th centile, grade 2–3 BPD) and quantified the unadjusted associations between these variables and LUSSc [21]. Those variables that were significantly related to LUSSc (p < 0.2) in unadjusted analyses were incorporated into the multivariable linear regression model via backward selection. For these computations, we employed a time-series regression analyses with autocorrelations to examine how each infant’s LUSSc changed over time with respect to their clinical/demographic characteristics. We report the model incorporating clinical characteristics associated with the weekly change in LUSSc.

This study was reviewed by the IRB at Ann & Robert H. Lurie Children’s Hospital of Chicago (IRB 2021-4037). Informed parental/guardian consent was obtained prior to participation. Exams were deferred to a later date if a patient was deemed too clinically unstable to tolerate the exam by the treating medical team. Redcap with two-factor authentication was used for confidential storage of patient clinical characteristics and STATA (v17.1) was used for statistical analyses. All hypothesis tests were 2-sided, and p < 0.05 was statistically significant.