Autonomic development in preterm infants is associated with morbidity of prematurity.

BACKGROUND: Previous studies have described an association between preterm birth and maturation of the autonomic nervous system (ANS); however, this may be impacted by multiple factors, including prematurity-related complications. Our aim was to evaluate for the effect of prematurity-related morbidity on ANS development in preterm infants in the NICU.
METHODS: We compared time and frequency domains of heart rate variability (HRV) as a measure of ANS tone in 56 preterm infants from 2 NICUs (28 from each). One cohort was from a high-morbidity regional referral NICU, the other from a community-based inborn NICU with low prematurity-related morbidity. Propensity score matching was used to balance the groups by a 1:1 nearest neighbor design. ANS tone was analyzed.
RESULTS: The two cohorts showed parallel maturational trajectory of the alpha 1 time-domain metric, with the cohort from the high-morbidity NICU having lower autonomic tone. The maturational trajectories between the two cohorts differed in all other time-domain metrics (alpha 2, RMS1, RMS2). There was no difference between groups by frequency-domain metrics.
CONCLUSIONS: Prematurity-associated morbidities correlate with autonomic development in premature infants and may have a greater impact on the extrauterine maturation of this system than birth gestational age. IMPACT: Autonomic nervous system development measured by time-domain metrics of heart rate variability correlate with morbidities associated with premature birth. This study builds upon our previously published work that showed that development of autonomic tone was not impacted by gestational age at birth. This study adds to our understanding of autonomic nervous system development in a preterm extrauterine environment. Our study suggests that gestational age at birth may have less impact on autonomic nervous system development than previously thought.

Introduction

An appropriately mature autonomic nervous system (ANS) is key for successful transition from the intrauterine to extrauterine environment at birth. Abnormal ANS development has been linked to neurodevelopmental, neuropsychiatric, and cardiovascular outcomes later in life.[1-4] The ANS undergoes a critical period of maturation during the latter half of gestation and during the early neonatal period. Previous studies have described an association between preterm birth and subsequent extrauterine maturation on the developing ANS.[5-8] However, there are multiple potential factors beyond birth gestational age (GA) that may influence ANS development in the extrauterine environment. We recently studied a relatively low-morbidity population of premature infants and found no significant association between duration of preterm extrauterine maturation and metrics of ANS development from birth to term-corrected age.[9] We proposed that birth GA may have less impact on ANS development than previously reported.

Prior studies investigating ANS maturation in preterm infants have usually involved infants with higher levels of prematurity-related complications.[6-8] Consequently, ANS development in preterm infants may be related more to medical comorbidities than to birth GA itself. To address this question, we sought to compare ANS development in preterm infants from two neonatal intensive care units (NICUs) with different levels of prematurity-related morbidities. In this study, we set out to determine the degree to which complications of prematurity impact ANS development by comparing two populations of preterm infants with differing levels of complications of prematurity.

Methods

Participants

We compared a panel of time- and frequency-domain ANS metrics in preterm infants (birth GA ≤ 34 weeks) from two different level-4 NICUs. The lower morbidity (LOW) NICU at Inova Fairfax Hospital is a large (108 bed) suburban community unit caring for inborn infants from a large obstetric program with about 10,000 deliveries per year. In a longitudinal prospective study of ANS maturation between May 2017 and May 2019, preterm infants were enrolled within 96 h of birth. Infants with a birth GA of 23 1/7 to 33 6/7 weeks, requiring NICU admission, without a suspected genetic condition at the time of enrollment, and expected to survive to NICU discharge at the time of enrollment were included in the study. The higher morbidity cohort (HIGH) cohort was recruited from the Children’s National Hospital NICU, an urban quaternary referral center, between June 2017 to November 2018, and included preterm infants with the same eligibility criteria. These subjects were retrospectively enrolled and their medical records and archived continuous physiological signals were analyzed. All infants in the Children’s National cohort (HIGH) were outborn and admitted before 72 h of age for medical or surgical complications of prematurity.

For the purposes of this study, we recorded significant postnatal medical complications at any time during the subjects’ NICU admission, including any of the following: culture-positive infection occurring at <7 days and ≥7 days of life (infection), necrotizing enterocolitis (NEC) diagnosis requiring medical or surgical intervention, presence of a patent ductus arteriosus (PDA) requiring medical or surgical intervention, brain injury on head ultrasound including periventricular leukomalacia and intraventricular hemorrhage grade 3 or 4 (neuro), diagnosis of bronchopulmonary dysplasia or chronic lung disease, retinopathy of prematurity at stage ≥1, and need for mechanical ventilation/intubation (mechanical ventilation). We also added up the total number of above conditions for each subject in the study for a summative “comorbidity score.” Data on medical comorbidities was abstracted from the medical record. We added the total number of the above conditions for each subject in the study to calculate a “comorbidity score.” The studies were approved by the Children’s National Hospital and INOVA Hospital Institutional Review Boards.

Physiologic signal collection

In the INOVA cohort (LOW), we collected up to 96 h of continuous electrocardiogram (ECG) recording directly from the infants’ bedside monitor immediately following enrollment and arrival in the NICU using custom software developed in either MATLAB 2018 (Mathworks, Inc., Natick, MA) or Labview (National Instruments, Austin, TX). The sampling rate of the ECG was 1000 Hz. After this initial prolonged ECG data collection period, we recorded up to 7 h of ECG data during the daytime at weekly intervals until NICU discharge.

In the Children’s National cohort (HIGH), ECG data were retrospectively retrieved from the Phillips Data Warehouse system (DWS) (IntelliVue MP-70, Philips, MA). We retrieved ECG data at weekly intervals from the Children’s National NICU admission until 42 weeks corrected GA during daytime hours, with the first collection occurring within the first 72 h of life in order to match the timing of the initial and weekly ECG data collection periods in the INOVA (LOW) cohort. The sampling rate of the ECG retrieved from the DWS is 250 Hz. All data were processed offline using computer programs developed in MATLAB.

Heart rate variability (HRV) analysis

We bandpass filtered the ECG signals between 0.5 and 60 Hz, identified R-waves, and calculated beat-to-beat interval (RRi).[10,11] We partitioned RRi into 10-min end-to-end epochs and calculated HRV metrics for RRi in every epoch. We estimated a power spectrum of RRi in every 10-min epoch using the Welch periodogram approach with a frequency resolution of 0.0167 Hz. We calculated low frequency (LF) and high frequency (HF) as the median of the logarithm of the power in 0.05–0.25 and 0.3–1 Hz, respectively. We also calculated normalized low frequency (nLF) and normalized high frequency (nHF) as the ratio of the sum of the powers in low- and high-frequency bands to the total power. Total power was defined as the sum of the powers in the frequency band between 0.05 and 2 Hz.

To calculate the time-domain metrics, we performed the detrended fluctuation analysis (DFA), which is a modified root mean square (RMS) analysis. In short, DFA involves the following five steps. First, the mean value of the RRi is subtracted from the RRi series and the cumulative sum, called the profile function, is calculated. Second, the profile function is partitioned into end-to-end windows containing “s” number of beats and fit to the profile in each window using a fourth-degree polynomial. Third, the local fluctuation function is calculated as the RMS of the difference between the profile and best fits. Fourth, the fluctuation function is calculated as the median of the local fluctuation function from all the windows. Last, steps 1–3 are repeated for different window sizes.[12] In our calculation, we varied the window size from 6 beats to one-fourth of the number of RRi samples in the 10-min epoch. Using the DFA fluctuation function, we calculated time-domain metrics as follows.[13,14] RMS1 (s) and RMS2 (s) were calculated as the maximum value of the fluctuation function for “s” between 15 and 50 beats and 100–150 beats, respectively. We also calculated alpha 1 and alpha 2 as the slope of the fluctuation and “s” in double logarithmic plot from the regions of 15–30 beats (short-term scale) and 35–150 beats (long term scale/ultralow frequency), respectively. The alpha metrics characterize the autocorrelation in the RRi, while the RMS (s) metrics characterize the variability in the RRi.[15] We calculated eight HRV metrics (LF, HF, nLF, nHF, alpha 1, alpha 2, RMS1, RMS2) in every 10-min epoch and averaged every metric over the entire recording period and obtained one representative value for the metric per study. In our earlier study, we showed that HF power calculated from DWS is compromised and hence we excluded this metric from further analysis.

For infants on the invasive ventilator support, RRi contains respiration-induced (artificial) oscillations in the short time scale. These oscillations affect the estimation of the alpha 1 estimation. To attenuate the artificial respiratory oscillations, we reconstructed the respiratory waves using the amplitude changes in the ECG signal. Using the reconstructed wave, we attenuated the respiratory components in the RRi using a frequency-dependent subtraction approach. We calculated alpha 1 using the respiration attenuated RRi.

Propensity score matching (PSM)

To balance the HIGH and LOW morbidity cohorts, PSM was used. The probability that a subject was in the HIGH cohort given a set of clinically important variables such as sex, maternal age, birth GA, birth weight Z-score, birth head circumference Z-score, maximum APGAR score, and length of stay was estimated using an unconditional multivariate logistic regression model. We then implemented one-to-one (1:1) nearest neighbor matching design.

Statistical analysis

Following matching, we performed separate latent growth modeling using a repeated-measures random-effects mixed model approach to determine the relationship between postnatal age (PNA) days on the selected variables. PNA days were treated as random effects. Models were assessed for a second-order polynomial function, and for consistency across models, given that only one term (alpha 1) was quadratic, a decision was made to treat PNA days as a linear function. Several various covariance matrices were evaluated using Akaike’s Information Criterion with an unstructured covariance matrix determined to be best fitting.

Two-sided significance testing was used. All analyses were conducted using SAS (ver. 9.3, Cary, NC) with statistical significance considered for p < 0.05. No adjustment was made for multiple comparisons. Descriptive and comparative statistics were performed for patient baseline characteristics and treatment outcomes. Data are presented as median (Q1–Q3) or frequency and percentage. Group comparisons were made using either t tests or Fisher’s exact test, where appropriate.

Results

We evaluated autonomic metrics over time from birth to NICU discharge or 42 weeks corrected GA in 56 infants (28 in each of the HIGH and LOW cohorts), propensity score-matched to minimize the effects of treatment site. The clinical and demographic characteristics of the infants are reported in Table 1. There was no significant difference in demographic characteristics, as expected, due to propensity matching. LOW infants were more likely to have been delivered by cesarean section (p = 0.04). Prenatal and postnatal conditions of the pregnant women and infants are presented in Table 2. There was no significant difference in the rates of multiple gestation, fetal growth restriction/intrauterine growth restriction, pre-eclampsia, or pregnancy-induced hypertension between the cohorts. However, LOW infants were more likely to have received antenatal steroids than HIGH infants (p = 0.03). HIGH infants had higher rates of PDA (p = 0.002), culture-positive infections (p = 0.053), and mechanical ventilation (p = 0.002) and had higher number of comorbidities per infant than LOW infants (average no. of comorbidities 2.6 ± 1.8 vs. 1.3 ± 1.4, p = 0.003).

Table 1

Demographic data.

Demographic data	HIGH (n = 28)	LOW (n = 28)	p value
Maternal age at birth (years)	31.5 (29.0 to 33.5)	31.5 (28.5 to 33.5)	0.95
Gestational age at birth (weeks)	28.8 (27.3 to 30.4)	28.6 (26.2 to 30.2)	0.58
Birth weight (g)	1140.0 (800.0 to 1482.5)	1090.00 (817.5 to 1485.0)	0.67
Z-score	−0.10 (−0.69 to 0.55)	0.05 (−0.41 to 0.39)
Head circumference (cm)	26.0 (23.3 to 28.0)	26.5 (23.5 to 29.0)	1.0
Z-score	−0.39 (−1.06 to 1.00)	0.61 (0.01 to 0.86)
Apgar 1 min	5.0 (3.0 to 7.5)	4.5 (3.5 to 7.0)
Apgar 5 min	7.5 (6.0 to 9.0)	8.0 (7.0 to 9.0)
Sex
Female	15 (54)	11 (39)	0.28
Male	13 (46)	17 (61)
Delivery type
Vaginal	8 (29)	2 (7)	0.04
C-section	20 (71)	26 (93)

Categorical data represented as N (%) and continuous data are reported as median (Q1–Q3) or as mean (SD); p values are calculated by Wilcoxon-Rank-Sum Test and T Test.

Bold values indicate statistical significance p < 0.05.

Table 2

Prenatal and postnatal comorbidities.

	HIGH (n = 28)	LOW (n = 28)	p value
Prenatal
Multiple gestation	5 (18)	2 (7)	0.42
FGR/IUGR	5 (18)	2 (7)	0.42
Pre-eclampsia	5 (18)	8 (29)	0.53
Pregnancy-induced HTN	2 (7)	0 (0.0)	0.75
Steroids	17 (61)	25 (89)	0.03
Postnatal
Length of stay	70.5 (56.0–90.0)	74.6 (42.0–106.6)	0.56
Neuro	7 (25)	2 (7)	0.14
NEC	2 (7)	1 (4)	1.0
ROP stage ≥1	9 (32)	7 (25)	0.77
BPD/CLD	14 (50)	13 (46)	1.0
PDA	11 (39)	1 (7)	0.002
ID	6 (21)	0 (0.0)	0.053
<7 days of life	1 (17)	0 (0.0)	—
≥7 days of life	5 (83)	0 (0.0)	0.43
Mechanical ventilation	23 (82)	11 (39)	0.002
Average no. of comorbidities	2.6 ± 1.8	1.3 ± 1.4	0.003

Categorical data represented as N (%) and continuous data are reported as median (Q1–Q3) or as mean (SD); p values are calculated by Fisher’s Exact Test.

Bold values indicate statistical significance p < 0.05.

Both HIGH and LOW cohorts showed changes over time for all autonomic metrics (Figs. 1 and 2 and Table 3). There was an absolute difference in the alpha 1 measure between HIGH and LOW cohorts, shown by nearly parallel maturational trajectories in Fig. 1 and in a significant cohort effect in Table 3. No other HRV metric showed a similar difference between cohorts. All other time-domain metrics (alpha 2, RMS1, RMS2) showed significant differences in the group-by-time interaction term (Table 3), reflecting different maturational trajectories between the groups. RMS1 and RMS2 absolute values were similar between cohorts soon after delivery, but the curves for the two groups diverged, with the LOW cohort showing more robust development than the HIGH cohort over time (Fig. 1). Alpha 2 was lower in the HIGH cohort compared with the LOW cohort at early time points but was higher in the HIGH cohort close to term equivalent age (Fig. 1). There was no difference between the frequency-domain metrics by cohort (LF, nLF, nHF; Table 3 and Fig. 2).

Fig. 1

Trajectory of development of time-domain metrics of heart rate variability.

Children’s cohort (HIGH) in blue, INOVA cohort (LOW) in red. PNA = postnatal age in days. a Alpha 1 development. b Alpha 2 development. c RMS1 development. d RMS2 development.

Fig. 2

Trajectory of development of frequency-domain metrics of heart rate variability.

Children’s National cohort (HIGH) in blue, INOVA cohort (LOW) in red. a Normalized low frequency (nLF) development. b Low frequency (LF) power development. c Normalized high frequency (nHF) development.

Table 3

Time, cohort, and cohort by time interaction effects.

	Alpha 1		Alpha 2		RMS1		RMS2		nHF		nLF		LF
	F	p	F	p	F	p	F	p	F	p	F	p	F	p
Time effect	76.0	<0.001	122.6	<0.001	139.3	<0.001	165.2	<0.001	32.8	<0.001	37.6	<0.001	154.4	<0.001
Cohort effect	7.4	0.03	3.4	0.226	0.4	0.56	0.1	0.75	2.2	0.14	0.7	0.41	0.03	0.85
Cohort × time interaction	0.6	0.05	29.4	<0.001	14.0	<0.001	6.9	0.01	2.8	0.09	1.0	0.32	2.5	0.11

F value from the F-test for overall significance for time effect, cohort effect, and cohort × time interaction term. p values are calculated by Fisher’s exact test.

Bold values indicate statistical significance p < 0.05.

Discussion

In this study, we found that ANS developmental trajectory in preterm-born infants may relate more to prematurity-associated medical comorbidity than to birth GA itself. In a preterm infant cohort with a low level of medical morbidity, we recently demonstrated that duration of extrauterine development did not seem to significantly impact ANS developmental trajectory from birth to NICU discharge.[9] Our finding is in contrast to previous other studies in which significant effects of prematurity on ANS development were reported.[7-9] These previous reports included preterm populations with higher levels of medical complexity.[7-9] The data from our current study supports the hypothesis that the medical complexity, in addition to birth GA, has an impact on ANS development in preterm infants.

The ANS undergoes a prolonged period of development that spans the fetal period and continues after birth.[16,17] The sympathetic division of the ANS begins to develop early in gestation shortly after the unmyelinated vagus (which has little function until after birth despite early structural development[18]) and continues to steadily mature throughout the fetal period. In the third trimester, the myelinated vagus nerve develops and begins to exert parasympathetic influence,[18,19] resulting in a steep increase in parasympathetic tone near term.[1,20,21] Within the central ANS, brainstem and hypothalamic centers develop first and develop multiple, complex connections with supratentorial structures, including the insula, the anterior thalamus, the anterior cingulate gyrus, the hippocampus, the hypothalamus, and the amygdala, later in gestation and during the early neonatal period.[5,22] Overall, sympathetic innervation develops and begins exerting influence earliest, with an increase in parasympathetic influence over time as higher cortical processes begin to integrate autonomic control during late gestation and the early neonatal period.[23,24]

HRV, the variation in beat-to-beat (R-R) intervals between heart beats, is the most widely used measure of autonomic function in the neonate. HRV metrics can be analyzed by both time and frequency-domain approaches. Short-term variability in time-domain metrics reflects parasympathetic function, while long-term variability is influenced by both sympathetic and parasympathetic function.[25] Sympathetic and parasympathetic tone and changes in heart rate related to baroreflex function are characterized by LF HRV,[26,27] while HF HRV reflects parasympathetic tone.[27]

In our study, there were significant differences in rates of culture-positive infection, PDA requiring treatment, need for mechanical ventilation, and overall number of comorbidities between the two cohorts. Infectious and inflammatory conditions such as neonatal sepsis and NEC are known to be associated with short-term depression and alterations of autonomic tone, which may precede symptoms of infection and inflammation by hours to days.[28-33] Preterm infants with early-onset sepsis are more likely to have lower ANS tone and abnormal neurologic outcomes.[4] Autonomic tone has previously been shown to differ in infants with specific types of medical comorbidities, including structural heart disease,[34] PDA,[35] and hypoxic–ischemic encephalopathy (HIE).[28,36] In term infants with HIE and stroke, brain injury has been associated with alterations in ANS function in the neonatal period.[37,38] This study adds to the current body of knowledge by showing that more medically complex preterm infants as a group have ANS maturation different from that of preterm infants with fewer medical comorbidities.

Alpha 1 is a variable that reflects quality of regulation of the ANS and is the only HRV metric studied for which there was an absolute difference between cohorts, with nearly parallel trajectories of development. In this study, alpha 1 started and remained lower in the HIGH cohort compared with the LOW cohort, suggesting lower autonomic regulation. Whether this difference is due to the underlying medical complexity within the Children’s National cohort (HIGH) or is an inherent preceding factor that increases risk for medical complications in preterm infants that would necessitate transfer to a quaternary referral NICU cannot be determined from this study. This is a question worthy of further investigation. The role of prenatal factors should also be examined. Our LOW and HIGH cohorts differed only in the administration of antenatal steroids, and we cannot say from this study whether this had any influence on ANS development or on alpha 1. Studies focused on prenatal influence on ANS function and fetal HRV may enable determination of whether infant differences precede birth and the development of early postnatal medical complications.

Alpha 2 decreased over time in the LOW cohort (as we have shown previously[9]) but increased over time in the HIGH cohort. This finding may reflect sympathetic over-compensation, perhaps due to greater prevalence of painful or stressful events in the HIGH cohort. RMS1 and RMS2 metrics were very similar between cohorts early after birth, but trajectories diverged over time, suggesting that maturation is slowed in the HIGH cohort. Given that our two cohorts are matched for major demographic criteria but differ in their levels of medical complexity, these data suggest that levels of medical complexity have a significant impact on ANS development as measured by time-domain metrics. However, given the lack of difference in frequency-domain metrics, there are likely multiple influences on ANS development in the preterm infant with medical complexity being one factor.

While our study focused on postnatal comorbidities, we would be remiss not to mention potential prenatal differences between the cohorts that could have contributed to alterations in ANS development. To this end, we found no difference between prenatal factors studied, except the rate of antenatal steroids. We hypothesize that pregnant mothers in the LOW cohort, who are more likely to be of higher socioeconomic status, may have utilized prenatal care to a greater extent than the HIGH cohort, resulting in earlier diagnosis of risk for preterm labor and greater rates of administration of steroids. However, due to the retrospective nature of this study, we cannot determine with certainty why this difference was present. We found similar time-domain metrics between groups at similar PNAs that diverged over time in the NICU. While prenatal factors are critical for early ANS development, the altered trajectory of postnatally measured RMS1 and RMS2 metrics of ANS development occurred over time in the NICU and differed between the cohorts. Interestingly, alpha 1 differed between cohorts at the time of birth. Future studies are needed to separate the influence of prenatal vs. postnatal factors on different metrics of HRV and their contribution to ANS developmental trajectory.

Our study has a number of strengths, including the longitudinal HRV recordings in two preterm infant cohorts from different NICUs, but also has limitations. First, we used preterm populations with different levels of medical comorbidity from two separate NICU sites. This was necessary in order to have an appropriate number of cases in each group for comparison, and we attempted to minimize the impact of the possible differences in care/hospital population through propensity matching. Nevertheless, there may be differences between NICU care and populations that we are not able to address in this study analysis. We chose to report the data as PNA rather than postmenstrual age in this study, in order to show ANS maturation for infants who had been matched for GA at birth. The median GA at birth in our two cohorts was slightly higher than in previous studies.[8] We do not have follow-up HRV measurements in these cohorts following NICU discharge and therefore cannot comment on long-term differences in ANS development between cohorts.

Because of the higher acuity at Children’s National (HIGH), more invasive ventilation is required, and we speculate that more handling procedures may have been required than for the relatively “well” preterm infants in the INOVA cohort (LOW). Practices such as minimal early medical handling, less invasive ventilation, early skin-to-skin contact, and prone positioning are known to impact preterm infant development and ANS maturation.[39-41] Neonates with depressed autonomic function at 24 h of age demonstrate inappropriate responses to stimulation, impaired cerebral autoregulation during painful procedures, and vagal activation during endotracheal suctioning, potentially complicating NICU care,[42] and we cannot rule out that depressed autonomic tone predisposed the Children’s National cohort (HIGH) toward more medical complexity or toward the need for transfer to this referral NICU. Likewise, the cohort at INOVA (LOW) was inborn while the infants at Children’s National (HIGH) were all outborn. The mothers of LOW infants may have had more access and interaction with their infants than mothers of HIGH infants, who may have had to travel long distances to visit their infants in the NICU. It is not possible for us to separate the effects of these factors from the impact of the medical morbidities themselves, and likely the changes we see are multifactorial. Infants with impaired autonomic function are at increased risk for maladaptive responses during intensive care and may be more likely to have events after discharge.[42,43] The LOW cohort had a significantly higher rate of cesarean-section delivery compared with the HIGH cohort. While the impact of delivery mode on HRV has not been well studied in preterm infants, our group has previously shown that delivery mode does not seem to have a significant impact on HRV metrics in full-term infants.[15]

Why it is that we showed significant differences in time-domain but not frequency-domain metrics of HRV between cohorts is not clear at this point, but several possible explanations warrant consideration. Time-domain metrics, including alpha 1, alpha 2, RMS1, and RMS2, may more accurately reflect ANS function under non-stationary conditions, such as critical illness. This is a question that warrants further investigation. Future, prospective studies should include real-time monitoring of preterm infants in the NICU setting and include data regarding handling procedures, skin-to-skin time, time with mother’s voice, and other interventions known to alter HRV, such as pacifier use.[40,44-48] In summary, we found that cohorts of preterm infants with differing levels of medical complexity show differences in their ANS development, specifically for time-domain metrics of HRV. The causal nature and direction of this association is the subject of ongoing investigation.