Maturation of the Autonomic Nervous System in Premature Infants: Estimating Development Based on Heart-Rate Variability Analysis

This study aims at investigating the development of premature infants' autonomic nervous system (ANS) based on a quantitative analysis of the heart-rate variability (HRV) with a variety of novel features. Additionally, the role of heart-rate drops, known as bradycardias, has been studied in relation to both clinical and novel sympathovagal indices. ECG data were measured for at least 3 h in 25 preterm infants (gestational age ≤32 weeks) for a total number of 74 recordings. The post-menstrual age (PMA) of each patient was estimated from the RR interval time-series by means of multivariate linear-mixed effects regression. The tachograms were segmented based on bradycardias in periods after, between and during bradycardias. For each of those epochs, a set of temporal, spectral and fractal indices were included in the regression model. The best performing model has R² = 0.75 and mean absolute error MAE = 1.56 weeks. Three main novelties can be reported. First, the obtained maturation models based on HRV have comparable performance to other development models. Second, the selected features for age estimation show a predominance of power and fractal features in the very-low- and low-frequency bands in explaining the infants' sympathovagal development from 27 PMA weeks until 40 PMA weeks. Third, bradycardias might disrupt the relationship between common temporal indices of the tachogram and the age of the infant and the interpretation of sympathovagal indices. This approach might provide a novel overview of post-natal autonomic maturation and an alternative development index to other electrophysiological data analysis.

1. Introduction

Premature infants represent 10% of the neonatal population and are at higher risk for developmental disorders that can lead to adverse outcome (Aylward, 2014). The investigation of maturation via multiple physiological biomarkers is part of the clinical practice to prevent lower cognitive, motor, or language outcomes later on in life (Franke et al., 2012; Koolen et al., 2016). A common probe to inspect the development of the neurovegetative functions or Autonomic Nervous System (ANS) is the heart-rate fluctuation, simply known as Heart-rate variability (HRV).

The guidelines of the adult HRV task force clearly specify the association between the different frequency tones of the tachograms and the stimulation of the ANS branches (Camm et al., 1996). The stimulation of the sympathetic branch is normally represented by the low-frequency band (LF, [0.04−0.15] Hz) of the HRV, while the high-frequency band (HF, [0.15−0.4] Hz) reflects the parasympathetic branch. The sympathovagal balance can be expressed by the power ratio of the two frequency bands $(\frac{LF}{HF})$, while the very-low-frequency band (VLF, [0−0.04] Hz) is usually associated to thermal and hormonal regulation. On the contrary, the fetal and preterm HRV frequency bands are still the subject of an intensive discussion in the literature. The early exposure to the ex-utero environment induces an aberrant sympathetic response and delays autonomic maturation (Smith et al., 2013; Javorka et al., 2017). The association between the common HRV frequency bands and the sympathovagal regulation is far less documented in infants and fetuses (Doret et al., 2015). Other factors are known to play a role in the definition of the oscillations of the heart rate, such as intermittent breathing cycles with high respiratory frequency and the actual delay in maturation of the autonomic nervous system. Therefore, David et al. (2007), Hoyer et al. (2013), and Doret et al. (2015) suggested that new ways to investigate the sympathovagal balance should be examined. Since the fetal heart-rate is characterized by a strong slow-wave baseline, David et al. redefine the frequency bands for fetuses as follows: VLF = [0.02−0.08] Hz, LF = [0.08−0.2] Hz, HF = [0.2−3] Hz. While adults normally present an HRV spectrum with two clear peaks at HF and LF (Camm et al., 1996), infants and fetuses have a 1/f spectrum up to 0.1 Hz (Karin et al., 1993). Consequently, the full description of the preterm ANS has to consider all the possible frequency bands (VLF,LF,HF) (Clairambault et al., 1992; Curzi-Dascalova, 1994; Mazursky et al., 1998; Longin et al., 2006). This could explain why the $\frac{LF}{HF}$ ratio can give contradictory results: Krueger et al. (2010) did not find any specific change in this ratio in a longitudinal study with preterm patients, while Longin et al. (2006) found a decrease in $\frac{LF}{HF}$ from preterm to term age. The rapid development and the unclear definition of the sympathovagal frequency bands might not give a simple interpretation of $\frac{LF}{HF}$ as it is for adults. Surprisingly, infants show greater changes in the absolute power of the three main bands VLF, LF and HF than relative power (Longin et al., 2006). Hoyer et al. (2013) argued that predominant principles of autonomic development are not only an increase in heart-rate variability but also an increase in the complexity and pattern formation. Consequently, HRV indices can be chosen to reflect these principles in order to describe the sympathovagal balance maturation. Pattern formation can be described by tachogram skewness and the new ratio $\frac{VLF}{LF}$, while the increasing complexity is characterized by an increasing HRV entropy. It should also be stressed that the computation of power ratios, such as $\frac{LF}{HF}$, requires stationarity, which can be questioned in the case of infants heart-rate time series. Therefore, Abry et al. (2010) and Doret et al. (2015) proposed fractal analysis as an alternative method to investigate the sympathovagal balance in fetal heart-rate. It focuses on quantities, such as oscillations or increments at different scales to tackle the absence of stationarity and determines specific relations between the fractal exponents (such as the Hurst Exponent) and the $\frac{LF}{HF}$ ratios. However, those methods were never applied to premature infants.

One example of non-stationarity is the presence of bradycardias. These are normally heart-rate drops below 70% of the heart-rate average, which last at least for 4 s and may be associated with apneas or hypoxias (Poets et al., 1993). These drops can alter oxygen saturation and blood flow, putting organs at risk of damage (Paolillo and Picone, 2013). Apneic spells that occur with bradycardias or hypoxic events are most likely to affect brain homeostasis. In addition, those physiological instabilities are the probable consequence of the immature respiratory system (Porges, 1995; Atkinson and Fenton, 2009). However, it has also been shown different HR reactions can be triggered by hyperoxia and hypoxias via chemoreception in term infants (Søvik et al., 2001). Additionally, Poets et al. specifically highlighted that bradycardias can occur independently from apneic or gas events (Poets et al., 1993). Bradycardias can be considered a consequence of heart-rate dysregulation, which can disrupt the state-space and the probability density function of the tachogram (Gee et al., 2017). Any proper model that tries to describe the development of the infants' ANS has to include not simply the slow variation of the basal heart-rate, but the sudden drops of the tachogram, independently from any other conditioning factor. Those non-stationary events possibly affect the most common HRV temporal or spectral features used in clinical practice, such as the standard deviation and the $\frac{LF}{HF}$ ratio. They are commonly used for the assessment of development outcome, sleep or pathologies diagnosis (Javorka et al., 2017), and any disruption of these features can bias conclusions made by the medical community. For example, bradycardias can forcefully increase the variability of the tachogram or its regularity.

In order to address the shortcomings using the studies outlined above and the lack of autonomic growth charts for premature infants, a new framework to describe autonomic maturation in healthy preterm babies has been provided. This research can be divided into two main strands. First, both spectral analysis and multifractal analysis have been employed to investigate the neurovegetative development of the sympathovagal balance and its complexity and track maturation. Second, the impact of bradycardias on both clinical and novel ANS maturation indices has been investigated. This study tries to provide a complete overview of autonomic maturation in premature neonates, including non-stationary events, such as bradycardias. The final clinical objective is to provide novel maturation charts for the premature autonomic nervous system in the first weeks of life and correct the effect of heart-rate events on common clinical HRV indices. Those normative charts might be used as references to investigate early-life and ex-utero factors that can deviate from normal premature development and define suitable therapies in the neonatal intensive care unit.

2. Methods

2.1. Dataset

The dataset consists of electrocardiograms (ECG) of 25 preterm infants, which were recorded at the Neonatal Intensive Care Unit of the University Hospital of Leuven. It was collected in a multimodal setting for another research study related to brain development and a sleep-stage analysis (Koolen et al., 2016; Dereymaeker et al., 2017). Inclusion criteria were as follows: a normal neurodevelopmental outcome at 9 and 24 months corrected age (Bayley Scales of Infant Development-II, mental and motor score >85), no severe brain lesions, assessed by ultrasound and not taking any sedative or antiepileptic drugs during the EEG registration. The sampling ECG frequency was 250 or 500 Hz and the average length of the recording was 4 h 44 min. An overview of the dataset is reported in Table 1, while a complete description is reported in Supplementary Tables 1, 2. The latter shows the heterogenous interperiod sessions among recordings, which indicate that the measurements were not scheduled at the same PMA for each patient. In addition, the Supplementary Materials 1, 2 show the ECG sampling frequency for each of the recording.

TABLE 1

Table 1. Demographics of the 25 patients: average duration of the tachogram in minutes (Duration_Rec), average duration of the annotated bradycardias in s (Duration_WB), average number of the annotated bradycardias (Number_WB), average RR amplitude during the bradycardia in ms (RR_WB), post-menstrual age in weeks (PMA) and gestational age in weeks (GA).

2.2. Pre-processing

The HRV represents the instantaneous fluctuations of heart rate and is usually expressed by the tachogram which visualizes the variations of the time interval between two consecutive R-peaks (RR intervals, RR_i). In order to compute a RR_i time series, the R peaks of the ECG have been detected via the Matlab toolbox by Moeyersons et al. (2019), which is based on enveloping procedure. This graphical user-interface also allows for correction and deletion in case of erroneous R-peaks. In case of a single missing R-peak, the value was replaced by using the following formula:

$\begin{array}{ll}{\widehat{R}}_t=\frac{R_{t-1}+R_{t+1}}{2},& (1)\end{array}$

where ${\widehat{R}}_t$ is the estimated position of the missed R-peak, while R_t−1 and R_t+1 are the location of the previous and following R-peak. In the case of two or more missing R-peaks due to ECG flat lines or muscle artifacts, which made the QRS detection impossible, the contaminated parts of the signal were discarded. In case that less 20 min of noise-free signal remained, the signal was discarded. The length of each recording (in min) is reported in Supplementary Tables 1, 2. The progressive number of the recording ID shows that some of them were fully discarded. The full overview shows that all included recordings had at least 50 min of available data (Duration_Rec column) resulting in a total of 74 recordings.

Besides the preprocessing of artifacts and before the feature extraction, we also dealt with the sudden drops of heart-rate, known as bradycardias. Although those phenomena are completely natural in the developing infant, they can suddenly increase the frequency content of the RR_i series. Therefore, traditional linear spectral and temporal analysis might not be suitable since the instantaneous variance and mean of the heart-rate can vary over time, as explained in detail by Gee et al. (2017). According to the same study, the heart-rate activity that precedes sudden drops might differ from the drops itself and other bradycardia-free periods. Consequently, bradycardias have been detected in the current studies before any further processing. Based on the definitions of apnea of prematurity and bradycardias by Paolillo and Picone (2013), the bradycardia spells were detected as sudden $\overline{R{R}_i}$ increases above $\theta =1.5*\overline{R{R}_i}$ that persist for more than 4 s, where $\overline{R{R}_i}$ is the median tachogram of the entire recording. We defined the onset of the bradycardia as the moment that the tachogram exceeds θ. Conversely, the offset was defined as the moment that the amplitude decreases below the same threshold. Subsequently, three different windowing strategies were applied:

1. Post-bradycardia (PB) windowing: the 10-min period that starts 10 s after the bradycardia offset was considered a candidate for features extraction. This window did not include the bradycardia itself.

2. Between-bradycardias (BB) windowing: all non-overlapping 10-min windows contained between bradycardic events were considered as candidate epochs for features' extraction. The first viable window was at least 10 min from the bradycardia offset in order to guarantee that the signal was stabilized.

3. Within-bradycardia (WB) windowing: a 10-min window was considered from the bradycardia onset. This windowing should involve both the information related to the heart-rate drop and the recovery period. Based on the definition of the PB, the PB and WB windowing schemes overlap almost fully for the period after the bradycardia offset, but WB is the only scheme fully containing the bradycardic event.

A visual description of the windowing scheme is reported in Figure 1. The gray dashed boxes highlight the three types of windows (WB, BB, PB) that can be determined in a single trace, while the dot-dash box shows typical bradycardia events. The average duration and amplitude of a bradycardia event are reported in Table 1, which also shows the average number of bradycardias in the entire dataset. A full overview recording by recording is reported in the Supplementary Tables 1, 2, which display the number of bradycardias as well as the average intensity and the average duration of heart-rate drops patient by patient. The Supplementary Material shows that some of the recordings did not have any heart-rate drop according to the reported definition. Therefore, the windowing scheme based on bradycardias was not applicable. In this specific case, the design choice was a segmentation in non-overlapping 10-min windows and assign the results of feature extraction to post-bradycardia windowing scheme (see Figure 2).

FIGURE 1

Figure 1. Visual representation of the three windowing schemes applied in this study: post-bradycardia scheme (PB), between-bradycardia scheme (BB) and within-bradycardia scheme (WB). The selected windows in each trace are indicated with gray dashed boxes, while the dot-dashed boxes show examples of annotated bradycardia. In case of BB windowing, a period T greater that 10 min is present between the end of the bradycardia and the first available window.

FIGURE 2

Figure 2. The block diagram shows the main steps of the study. For each RR signal, artifact preprocessing is performed and associated resampling of the tachogram. The signal is split in different windows according to the scheme of Figure 1. For each of these epochs, temporal, spectral and fractal features undergo a grand-median process if there is more than one epoch per scheme. The three datasets are then used to estimate the age of the recording in a linear mixed effects (LME) regression.

2.3. Feature Extraction

In each of the windows defined according to the PB, BB, and WB schemes, a set of temporal, spectral and fractal features were derived to describe the autonomic nervous system of the premature infants and its relationship with development. These features were chosen based on the principles of variability increase, complexity increase and pattern formation by Hoyer et al. (2013). An overview of the different attributes is reported in Table 2.

TABLE 2

Table 2. Overview of all computed features.

2.3.1. Temporal Indices

Based on the most common guidelines related to HRV processing (Camm et al., 1996; Javorka et al., 2017), the first- and the second-order moments of the RR_i, i.e., the mean of the tachogram (μ_RR) and the standard deviation (σ_RR), were computed in order to assess the level of the variability.

2.3.2. Spectral Analysis

The sympathovagal activity is normally assessed by the computation of the spectral power in the different HRV frequency bands (Camm et al., 1996). Unlike adults, premature infants have a higher mean heart rate with very slow oscillation around it (Clairambault et al., 1992; David et al., 2007). Therefore, the frequency bands of the premature patients were defined as follows: VLF = [0, 0.08] Hz, the low-frequency LF = [0.08, 0.2] Hz and high frequency HF = [0.2, 3.0] Hz. Additionally, the RR time series of the premature infant can be non-stationary due to a series of events, like bradycardias or other heart-rate dysregulation. Therefore, the power spectral density was computed with time-frequency (TF) methodologies, which allows us to investigate the principle of pattern formation as discussed by Hoyer et al. (2013). Namely, the HRV power spectral density (PSD) was estimated with three specific approaches: the Welch's periodogram, the quadratic smoothed pseudo Wigner-Ville distribution (SPWD) (Orini et al., 2012) and the continuous wavelet transform (CWT) (David et al., 2007).

Given a fixed window size, Welch's algorithm estimates multiple periodograms in overlapping subwindows and averages them. The Welch's window length was set at 3 min and the overlap at 50%. Based on the suggestions of the HRV guidelines (Camm et al., 1996), the window length was set to investigate the very-low-frequency band.

One can also estimate the instantaneous autospectrum S_RR(t, f). Based on the CWT of the tachogram

$\begin{array}{ll}{W}_{RR}(t,s)={\displaystyle {\int }_{-\infty }^{+\infty }R}R(\tau ){\psi }^{*}\left(\frac{t-\tau }{s}\right)d\tau ,& (2)\end{array}$

S_RR(t, f) can be computed as the scalogram of the wavelet transform of the signal as follows:

$\begin{array}{ll}{S}_{RR}(t,f)=|W_{RR}(t,f){|}^2,& (3)\end{array}$

where ψ is the mother wavelet (Analytic Morlet), while s stands for the scale of the wavelet transform and, in the general, s ≈ f⁻¹. However, the S_RR(t, f) based on CWT risks to be distorted by interference terms which can be present with linear time-frequency approaches. Therefore, Orini et al. (2012) proposed to estimate the instantaneous autospectrum via a quadratic time-frequency distribution, such as SPWD. Then S_RR is then estimated as follows:

$\begin{array}{ll}{S}_{RR}(t,f)={\displaystyle \int }{{\displaystyle \int }}_{-\infty }^{+\infty }{\Phi }_{RR}(\tau ,\nu )A_{RR}(\tau ,\nu )e^{j2\pi (t\nu -\tau f)}d\tau d\nu ,& (4)\end{array}$

where A_RR(τ, ν) is the ambiguity function, which is defined as the Fourier Transform of the time-dependent auto-correlation of RR(t) as follows

$\begin{array}{ll}{A}_{RR}(\tau ,\nu )={\displaystyle {\int }_{-\infty }^{+\infty }}RR(t+\tau /2)R{R}^{*}(t-\tau /2)e^{-j2\pi t\nu }dt,& (5)\end{array}$

The smoothing of the time-frequency cross-coupling in (4) is done via the exponential kernel in the ambiguity domain defined as

$\begin{array}{ll}{\Phi }_{RR}(\tau ,\nu )=\exp \left\{-\pi {\left[{\left(\frac{\nu }{{\nu }_0}\right)}^2+{\left(\frac{\tau }{{\tau }_0}\right)}^2\right]}^{2\lambda }\right\},& (6)\end{array}$

Following the study by Widjaja et al. (2013), ν₀, τ₀, λ were set to 0.050, 0.046, and 0.3, respectively, leading to a kernel function with a TF resolution of [Δt, Δf] = [10.9 s, 0.03 Hz]. Both Welch's approach and the CWT were computed with MATLAB subroutines, while the implementation details and software download for the SPWD are reported in Orini et al. (2012).

Based on the given methodologies, the instantaneous power in the β = [f₁, f₂] band of interest can be obtained as

$\begin{array}{ll}{P}_{\beta }(t)={\displaystyle {\int }_{f1}^{f2}}{S}_{RR}(t,f)df.& (7)\end{array}$

In particular, one can compute the absolute power in three bands VLF, LF and HF as the reported integral in (7). Besides, the ratios $\frac{VLF}{LF}$ and $\frac{LF}{HF}$ were also computed alongside two indices to represent the normalized LF power: $\frac{LF}{VLF+LF}$ and $\frac{LF}{HF+LF}$. In case of Welch's algorithm, there is no dependency from the time variable t. On the contrary, CWT and SPWD generate a time series for each selected frequency band, as highlighted by (3), (4), and (7). In order to obtain one single value for each window, the median of this time-series was taken into account. The set of spectral features derived for each methodology is reported in the central block of Table 2.

2.3.3. Multifractal Analysis

Since spectral analysis requires stationarity of data and the very definition of the tachogram series frequency bands have been questioned, the HRV was also analyzed according to the fractal or multifractal paradigm. The multifractal analysis aimed to describe the principle of complexity increase discussed by Hoyer et al. (2013). As shown in Doret et al. (2015), the infant's tachogram is a fractal or scale-free signal, which presents a power-law decay spectrum as follows:

$\begin{array}{ll}{S}_{RR}(f)=\left|C\right|f^{-2(H_{exp}-1)}& (8)\end{array}$

where H_exp is known as the Hurst exponent and controls the decay of the power function. H is also a representative parameter for fractal time series and there can be more than one exponent for each signal. A signal with one single exponent is commonly known as monofractal, while a signal with multiple exponents h is known as multifractal (Ivanov et al., 1999). Small values of h represent sharp and transient regularity or singularity, while large values represent smooth changes (Leonarduzzi et al., 2010).

An efficient method to determine the amount of exponents or singularities h is the multifractal formalism based on the wavelet transform modulus maxima. The discrete wavelet transform (DWT) decomposes the signal RR(t) into elementary time-frequency components based on different scales a. Large scales describe smooth and low frequency oscillations, while small scales describe the sharp transitions in the signal. According to Popivanov et al. (2006) and Wendt et al. (2007), a partition function $Z(a,q)=Z_L(2^j,q)$ can be estimated using the wavelet leader L_f(j, k), as follows:

$\begin{array}{ll}Z(a,q)=Z_L(2^j,q)=\frac{1}{n_k}{\displaystyle \sum _{k=1}^{n_k}}|L_f(j,k){|}^q~2^{j\tau (q)},& (9)\end{array}$

where L_f(j, k) represents a specific type of wavelet transform, where the maximum of wavelet coefficients is considered in a narrow time neighborhood. More details can be found in Wendt et al. (2007) and Abry et al. (2010).

One can immediately notice the similarity between (8) and (9), especially between the Hurst exponents and τ(q). For certain values of q, the scaling exponent τ(q) (SE) has a specific meaning: for positive, respectively negative q, Z(a, q) reflects scaling of large, respectively short fluctuations. In general, for each q, the partition function exhibits a power-law decay characteristic, such as the power spectrum of 1/f noise (8). The scaling exponents τ(q) (SE) associated with this decay can be obtained by computing the slope of Z vs. the scales in a log-log diagram from a certain scale $a_1=2^{j_1}$ to a certain scale $a_2=2^{j_2}$. The log-transform clearly shows the advantage to define scales as power quantities. In case of a monofractal signal, τ(q) is a linear function of q and H_exp, which is τ(q) = qH_exp−1 (as also shown in 8). In case of a multifractal signal, the τ(q) is a non-linear function of the local exponents h and its fractal dimensions D(h), known also as the singularity spectrum (SS) (Popivanov et al., 2006).

The fractal paradigm fully describes the properties of the signal by means of the singularity spectrum D(h), obtained as the Legendre transform of τ(q) (Abry et al., 2010). This function represents the embedding dimensions in the function of the different singularities of the signals and two main D(h) attributes are normally derived: the maximum and the width of the D(h) distribution, known also as the parameters C₁ and C₂. They can be computed as cumulants or coefficients of the Taylor expansion of τ(q) and they are used to represent the main Hurst exponent of the multifractal signal (H_exp or C₁) and the “variety” of fractals in the time series (C₂).

The multifractality parameters (H_exp, C₂) were computed in the entire non-overlapping window according to three schemes discussed in section 2.2. Specifically, the multifractal features were derived using the Wavelet p-Leader and Bootstrap based MultiFractal analysis (PLBMF) MATLAB toolbox, described in Wendt et al. (2007). This toolbox can be downloaded from https://www.irit.fr/~Herwig.Wendt/software.html. A fundamental design choice is the scale range $\left[2^{j_1},2^{j_2}\right]$ from which the exponent τ(q) is estimated from (9) (Wendt et al., 2007; Abry et al., 2010). In case of HRV, the exponents [j₁, j₂] are normally set equal to [3,12]. Given the fact that the scale can be written as $a=2^j=(f_s/2)/f$ with f_s as sampling frequency of the signal, it follows that the range [j₁, j₂] = [3,12] approximately represents the frequency band ≈ [0.375, 0.001] Hz with f_s = 6 Hz. In case that [j₁,j₂]=[5,12], the chosen scale range approximately represents the frequency band ≈ [0.094, 0.001] Hz. It is clear the first range considers part of the HF band, while the latter solely focuses on the combination of LF and VLF. A window size of 10 min does not allow an investigation of oscillations lower than up to 0.001 Hz. However, since the VLF band of infants is limited to 0.01 Hz, we can state that the scale range [j₁,j₂] = [5,12] fully describes slow-waves of infants HRV. Since the chosen scale range might influence the multifractal attributes, both ranges were tested in this study to investigate which frequency bands mostly reflects the sympathovagal balance. In fact, the main Hurst exponent H_exp or C₁ parameter is able to influence the ratio $\frac{LF}{HF}$. Based on (8), one can rewrite the spectral ratios as follows:

$\begin{array}{ll}\frac{LF}{HF}=\frac{{\int }_{f_m}^{f_I}{S}_{RR}(f)df}{{\int }_{f_I}^{f_M}{S}_{RR}(f)df}=\frac{(f_I^{2-2H_{exp}}-f_m^{2-2H_{exp}})}{(f_M^{2-2H_{exp}}-f_I^{2-2H_{exp}})}& (10)\end{array}$

where [f_m, f_I],[f_I, f_M] represent the frequency bands of LF and HF. Taking into account that the Hurst exponent and the $\frac{LF}{HF}$ are related and taking also into account that the chosen [j₁, j₂] decides which frequency bands the multifractal parameters are related to, the scales investigation of the fractal properties can shed a light which bands mainly reflect the sympathovagal activity. The set of fractal features derived for each methodology is reported in the last block of Table 2.

2.3.4. Algorithmic Pipeline and Statistical Analysis

The processing pipeline of the current study is reported in Figure 2. For each HR time series, the signal was split according to the PB, BB, and WB windowing scheme reported in Figure 1 and all the features reported in Table 2 were extracted. Besides artifact removal, a fundamental preprocessing step is the resampling of the tachogram. The behavior of non-linear features can depend on the sample rate, as also shown by the definition of the scales and their range for the multifractal parameters (section 2.3.3). Based on the findings by Bolea et al. (2016), the following sampling frequencies were tested for fractal indices: [6, 8, 12] Hz. In contrast, the sampling frequencies for the spectral and temporal indices was set to 6 Hz in order to include the higher respiratory frequency of premature infants (Camm et al., 1996; Javorka et al., 2017). The data were resampled with linear interpolation. The results are reported for only one sampling frequency in the main manuscript; however, a more complete discussion of the use of different sampling frequencies and the associated results are included in the Supplementary Material.

As described in section 2.3, the tuning and design parameters for the spectral and fractal analysis were chosen in accordance with the absence of stationarity and the persistent slow-wave baseline of the premature HRV signal. The necessity to investigate long-range fluctuations and a recovery period after events, such as bradycardias justify the segmentation in 10 min. Normally, time-frequency approaches use windows longer than 600 s to describe evolution in HRV spectrum (Orini et al., 2012; Widjaja et al., 2013) and the fractal indices also require windows of this size to fully investigate changes in regularity (Abry et al., 2010; Doret et al., 2015). Additionally, the BB and WB schemes can generate a set of windows and therefore an array of features based on the number of bradycardias present in each recording (on average, seven bradycardias per recording, as reported by Table 1). In order to obtain one representative value for each recording in each windowing scheme, the median of this array of attributes over the different windows was computed for each recording, as highlighted by the grand-median block in Figure 2.

After the features' extraction process and the grand-median step, three datasets were then obtained according to the three different windowing schemes. The number of features extracted for each dataset was then 27 in total: 21 for the spectral attributes, two for the temporal ones and four for fractal indices, as shown in Table 2.

In order to investigate the ANS maturation, the HRV features were used to estimate the PMA of the patient, as shown by the last block of the diagram in Figure 2. Since the PMA is known for each recording, a linear mixed effects (LME) regression model was developed for each dataset with PMA as response variable (Lavanga et al., 2018). The actual regression consisted of two steps. First, the features were selected via the least absolute shrinkage and selection operator (LASSO) due to the high number of features, after that the absolute power features were log-transformed. Specifically, the LASSO was repeated for 20 iterations on the entire dataset and the features which were selected in more than 40% of the total number of iterations (eight iterations out of 20) were included in the regression model (Lavanga et al., 2018). Second, a linear mixed-effect regression model was built with the selected subsets with multiple random splits of the data. In particular, the dataset was split into 70% training set and 30% test set for 20 iterations and the model was developed on the train set and tested for test set for each iteration. The performance was then assessed as mean absolute error MAE on the test set, as well as explained variance R², both on train and test set ($R_{train}^2$, $R_{test}^2$). We also reported the p_value of the F-statistics for each iteration. The results were reported as median(IQR) (where IQR stands for InterQuartile Range) over the 20 iterations. A linear mixed effect model requires the definition of a grouping variable that introduces the random effect, and the patient ID was taken as a grouping variable since a set of one or more recordings belong to a patient [as discussed in a previous study (Lavanga et al., 2018)]. Furthermore, the LME regression with the LASSO procedure was not simply examined for the entire subset of features, but also for the three subset feature groups: temporal, spectral and fractal attributes. In case of temporal features regression, the LASSO step was not performed.

On top of that, the trend for the ANS features throughout the patients' development was also reported as median(IQR) in three age groups (PMA ≤ 31 weeks, PMA ∈(31−36] weeks, and PMA > 36 weeks) as well as Pearson correlation coefficient with PMA.

3. Results

The overview of the dataset is reported in Table 1, which shows certain traits of the annotated bradycardias. The mean length is around 18 s. On average, there are seven bradycardias per recording and the mean RR_i during bradycardias is ~631 ms. The overview shows that the infants have an average PMA of 34 weeks and 35 recordings are collected in the range (32–36] weeks. A total of 22 recordings is collected in the first days of life, while 17 recordings were included from the weeks close to discharge. A detailed overview of each recording is reported in Supplementary Tables 1, 2.

Figure 3 and Table 3 report the trends in the three different windowing schemes for the following features: the mean μ_RR and standard deviation of the HRV σ_RR, the absolute power in the LF band (and its logarithmic transform), the relative $\frac{LF}{LF+VLF}$ power, the Hurst Exponent in the range [j₁,j₂]=[5,12] H_exp[j₁,j₂=5,12] and the width of the singularity spectrum in the same range C₂[j₁,j₂=5,12]. The overview of all features for all different windowing schemes (PB, BB, and WB) are reported in the Supplementary Tables 3–5. Figure 3 reports the results for the windowing scheme for the within-bradycardia epochs on the left column, while the results for the between-bradycardia epochs are shown on the right column. The power in the LF band and the relative LF power ($\frac{LF}{LF+VLF}$) increase with increasing PMA in both scenarios: in particular, Pearson correlations ρ_xy are 69% (72% with logarithm transform) and 64% for bradycardia epochs, respectively, while ρ_xy are 71% (71% with logarithm transform) and 48%, respectively, for between-bradycardia windows. Concerning the post-bradycardia period, the Table 3 shows a Pearson correlation of 69% for the power in LF band and 57% for $\frac{LF}{LF+VLF}$. Results are here reported for the wavelet approach, but the other spectral methodologies exhibit similar trends (see Supplementary Tables 3–5). In addition, the Hurst exponent (derived as the c₁ of the singularity spectrum) decreases with development (ρ_xy are −45% in the bradycardias scenario, −47% in post-bradycardia scenario, and −50% in the between-bradycardia scenario), while the width of the singularity spectrum (c₂ parameter) increases with increasing PMA (ρ_xy are 54% in the bradycardia scenario, 45% in the post-bradycardia scenario and 43% in the between-bradycardia scenario). The greatest contrast was found with the variability of the heart-rate, σ_RR. While the standard deviation increases with infants' maturation in the between-bradycardia epochs, the σ_RR does not increase with age within the bradycardic event. Moreover, it is higher in the bradycardia epochs than in the between-bradycardia scenario (ρ_bradycardia = −4% vs. ρ_between = 64% with p_v = 0.77 vs. p_v ≤ 0.01).

FIGURE 3

Figure 3. (A–J) The figure shows the linear-mixed effect regressions between the post-menstrual age and the following HRV features: the standard deviation of the tachogram σ_RR, the absolute and the relative power in the LF band $({log}_{10}(LF),\frac{LF}{LF+VLF})$, the Hurst exponent H_{exp, [j1,j2=5,12]} and the parameter C₂. The sampling frequencies for the fractal indices is f_s = 8 Hz. The left column—magenta circles report the results for the bradycardia epochs, while the right column—indigo diamonds the results for the between-bradycardia epochs. ρ is the correlation coefficient with the associated significance p_value.

TABLE 3

Table 3. The main temporal, spectral and fractal features are reported for the three windowing schemes (PB, BB, and WB).

Table 4 shows the regression results for the linear mixed-effect models, while Table 5 reports the features selected by LASSO. Those two tables report the results for the three different windowing schemes (PB, BB, WB) in three different blocks, while the rows report the results for the different feature groups (temporal, spectral and fractal attributes) and sampling frequencies. The different columns, respectively report the explained variance in the training set ($R_{train}^2$), the mean absolute error (MAE) and the explained variance in the test set ($R_{test}^2$). The best performance is reached for the combination of all features in the PB scheme ($R_{train}^2=0.75$, MAE = 1.83 weeks, $R_{test}^2=0.57$) as well as between bradycardias ($R_{train}^2=0.68$, MAE = 1.56 weeks, $R_{test}^2=0.59$). During the bradycardia event (WB), the best performance is achieved with the spectral features ($R_{train}^2=0.73$, MAE = 1.9 weeks, $R_{test}^2=0.62$). Table 5 shows that the selected features are the absolute spectral power in LF and VLF together with C₂ parameter in the range [j₁,j₂] = [5,12] for the first two schemes. For the WB scheme, the selected feature is simply the power in LF band.

TABLE 4

Table 4. Linear mixed-effect model performances.

TABLE 5

Table 5. LASSO selected features for the linear mixed-effect model.

Figure 4 shows that the relationship of (10) between the H_exp and the ratios $\frac{VLF}{LF}$ and $\frac{LF}{HF}$. The first row shows the relationship between $\frac{VLF}{LF}$ and H_{exp, [j1,j2=5,12]} in the three windowing schemes: WB in magenta circles, PB in light-blue squares and BB in indigo diamonds. The Pearson correlation coefficients are 21, 49, and 43%, respectively. The second row shows the relationship between $\frac{LF}{HF}$ and H_{exp, [j1, j2 = 3, 12]} in the same three schemes. The Pearson correlation coefficients are 18, 20, and 36%, respectively.

FIGURE 4

Figure 4. (A–F) The figure shows results for the linear-mixed effect regression that models the relationship between H_{exp, [j1,j2=5,12]} and $\frac{VLF}{LF}$ in the first row and between H_{exp, [j1, j2 = 3, 12]} and $\frac{LF}{HF}$ in the second row. The three columns, respectively represents bradycardia epochs (magenta circle data points), post-bradycardia epochs (light-blue squares data points), and between-bradycardia epochs (indigo diamonds data points). ρ is the correlation coefficient of the regression and p_value represents the significance of the correlation.

4. Discussion

This study provides an overview of the autonomic nervous system maturation in preterm infants and aims to estimate the post-menstrual age of the infants based on the HRV. Since the neonatal tachogram is a signal characterized by lack of stationarity and strong slow-wave baseline (Abry et al., 2010; Hoyer et al., 2013; Doret et al., 2015), the current study investigated the maturation of sympathetic and parasympathetic branches with the combination of temporal, spectral and fractal indices. Three main novel findings can be reported. First, Table 4 shows that the maturation of infants can be assessed with different spectral and fractal HRV indices with comparable performances to other maturation models for fetal and preterm development by Hoyer et al. (2013), De Wel et al. (2017), and Lavanga et al. (2017, 2018). Second, Figure 4 reports that the spectral ratio $\frac{VLF}{LF}$ and the Hurst exponent in the range [j₁,j₂]=[5,12] are more correlated than the $\frac{LF}{HF}$ and the Hurst exponent [j₁,j₂]=[3,12]. This might indicate that neonates do not have a sympathovagal balance that relies on the typical interplay between LF and HF (Abry et al., 2010; Doret et al., 2015). Third, the bradycardias can impact HRV maturational features, especially the most common temporal indices that are used in clinical practice (Javorka et al., 2017), such as the $\frac{LF}{HF}$ ratio and the standard deviation (Tables 3, 4). Additionally, the relationship between spectral ratios and H_exp is strongly reduced in the WB scheme, as stressed by Figure 4.

The different age models that were derived in this study can outperform or can be compared to the other developmental models reported in the literature (Hoyer et al., 2013; Lavanga et al., 2018). Specifically, Table 4 highlights the capacity of spectral features to outperform all other features in the PMA estimation in all three windowing schemes. Furthermore, the LF power is consistently selected by LASSO for all the different f_s and with any type of windowing scheme. These results are not simply in line with a decrease of $\frac{VLF}{LF}$ by Hoyer et al. (2013), but they are also supported by other clinical findings. Namely, an increase of the short-term variability of the tachogram was found during the first days of life (De Souza Filho et al., 2019) and the absolute LF power can discriminate preterm and full-term infants with 84% accuracy (Mulkey et al., 2018). However, the highest performances in the PB and BB schemes are achieved when the fractal and spectral features are combined, as highlighted in bold in Table 4 and suggested by the concomitant increase of entropy and short-term variability of HRV found by De Souza Filho et al. (2019). Interestingly, the highest performances are also achieved when the between bradycardias epochs are considered (MAE = 1.56 weeks), which further reveals a bias effect of bradycardias in the description of autonomic maturation. In line with other studies (Clairambault et al., 1992; Curzi-Dascalova, 1994; Longin et al., 2006; Hoyer et al., 2013), we found that the tachogram mean μ_RR and its standard deviation σ_RR increase with maturation together with the absolute power in all investigated frequency bands. These findings also confirm the results by Mulkey et al. (2018), who found a greater discrimination ability of the absolute power than relative indices in the classification between preterm and full-term neonates. The lack of stationarity and the 1/f spectrum behavior makes it difficult to describe the autonomic maturation without all the frequency bands in place or the fractal indices H_exp and c₂, which are also involved in the estimation of the development (Table 5).

These findings seem to suggest that the ratio $\frac{LF}{HF}$ is not the most suitable index for the sympathovagal balance and the common HRV frequency bands are suitable for infants. As anticipated, David et al. noticed that the fetal heart-rate has such enhanced slow-wave baseline, which increases the power in the VLF band such that both David et al. (2007) and Hoyer et al. (2013) used the ratio $\frac{VLF}{LF}$ as a possible index to describe the sympatho-vagal interplay. This approach is confirmed by the results in Figure 4. As discussed by Abry et al. (2010) and Doret et al. (2015), the spectral ratio is linked to H_exp via (10). The panels suggest that the ANS modulation and its fractal regularity lie in the lower-frequency bands since the H_exp is more correlated with $\frac{VLF}{LF}$. The Pearson correlation coefficients ρ_xy between $\frac{VLF}{LF}$ and H are respectively 21, 49, 43% according to different windowing schemes. These are significantly higher than ρ_xy between $\frac{LF}{HF}$ and H_exp (18, 20, and 36%). It is important to notice that the H_exp matches the spectral ratio if its estimation range [j₁, j₂] matches the frequency bands with most of the exponential decay in the PSD. In this study, [j₁, j₂ = 5, 12] covers specifically the lower frequency bands. These frequency band's importance is confirmed by the features selected by LASSO (Table 5). In line with Doret et al. (2015), the current findings clearly suggest a redefinition of $\frac{LF}{HF}$ with an extension of frequency bands from the most common adults' range, e.g., LF = [0.02−0.15] Hz and HF = [0.15−1.3] Hz. They also highlight the greater prominence of the slower oscillations in the description of premature ANS.

However, the results also highlight the disruptive role of bradycardias in maturation analysis. As anticipated, the best regression results are achieved in the between bradycardia epochs (Table 4), and the relationship between the spectral ratios and the H_exp is disrupted with WB windowing (panels with magenta circle in Figure 4). Most importantly, the relationship between the temporal features and maturation is lost, as highlighted by the poor R² (Table 4 and panels with magenta circles in Figure 3). In addition, Gee et al. (2017) observed that the LF power, the variance and the regularity of the heart-rate increase before bradycardias. The results in Figure 3 and Supplementary Tables 4, 5 support this increase in variance and regularity (as can be easily noticed by the y-axis of σ_RR or any other features of the left column in Figure 3). This finding clearly implies that the exclusion of bradycardias is fundamental whenever using the standard deviation and the mean of the tachogram to assess the maturation of ANS. Figure 4 also shows that bradycardias annihilate the expected relationship between the spectral ratios and the Hurst exponent. This is further proof that bradycardias disrupt the vagal tone (Porges, 1995), which can distort the PSD power-law in Equation (8). However, Table 4 shows that the autonomic age models within bradycardia can maintain comparable performance to the other two windowing schemes thanks to the spectral features. In particular, Table 5 confirms that the most selected power attribute is the power in LF = [0.08 − 0.2] Hz, as further proof of the central role of this range in the bradycardic event and ANS maturation (David et al., 2007; Hoyer et al., 2013; Gee et al., 2017). The role of bradycardias might be further investigated via proper conditioning with respiration or SpO₂ signals. Unfortunately, saturation and respiration data were not available in this dataset and ECG derived respiration does not properly estimate the breathing activity due to the small rib cage of the infants and possible skin stripping (Pereira et al., 2017). As already mentioned, bradycardias can also arise independently from hypoxic or apneic events (Poets et al., 1993). Therefore, this study solely looked at the HR instabilities, but future analysis might comprehend a full cardiovascular assessment to describe the ANS maturation and take into account the influence of the recovery period from the bradycardia spike.

It is also important to highlight some limitations of the current study. One may object to the exclusion of proper sleep-staging in the current analysis, as normally done by Curzi-Dascalova (1994). However, the specific focus on the bradycardia effect strongly limits the number of windows available. On top of that, bradycardias are events normally associated with active sleep (Porges, 1995) and most of the annotated bradycardias in this study were found during states that were not associated with quiet sleep. Similarly, one may also find the number of patients limited, but it was caused by the difficulties in the follow-up. All the included patients had normal developmental outcome at 2 years and the development assessment process is normally characterized by large drop-outs. Concerning the methodology, the different spectral methods (Welch, Wavelet, and Wigner-Ville) show very similar spectral trends, but LASSO more frequently tends to select time-frequency distribution features (Wavelet and Wigner-Ville, Supplementary Tables 4, 5). Although there are studies that claim the superiority of the quadratic time-frequency methods (Orini et al., 2012; Widjaja et al., 2013), the current findings show the wavelet approach would suffice for the spectral analysis. Concerning other sampling frequencies, negligible differences were found and a full discussion is reported in the Supplementary Material. It should also be mentioned that the current study was designed to provide growth charts based on the three principles by Hoyer et al. (2013): increase of variability, pattern formation and increase of complexity. We decided to replace the multiscale entropy with the multifractality due to the non-stationarity of the HRV signals and the relationship between spectral and fractal features (especially the spectral ratio in 10). However, the state space and the increase of complexity could also be monitored by entropy measures, such as the multiscale entropy or the asymmetries of HRV (Porta et al., 2008). In order to have a complete overview of the autonomic changes, future studies should not only analyze the specific frequency bands and the fractal properties of the signal, but they should include changes in the probability densities and the conditional entropy of the signal (Porta et al., 2015, 2017). This might not only provide a universal framework to describe the development of the autonomic nervous system in infants but also a further assessment of the bradycardia impact on the state-space of the tachogram. As also mentioned earlier, a full extension of this analysis should also include respiration signals and arterial blood pressure to provide a complete overview of the cardiovascular regulation of the premature infant (Porta et al., 2006; Montalto et al., 2014).

In a nutshell, the HRV analysis might be a useful tool for development monitoring, but two important factors have to be taken into account. First, the neonatal HRV is characterized by a very-low-frequency tone which requires a redefinition of the different frequency bands to the autonomic stimulation. Second, bradycardias have a disruptive role in the assessment of maturation.

5. Conclusion

The present study investigated the maturation of the preterm autonomic nervous system by means of temporal, spectral and fractal features of HRV. Three main findings can be reported. First, infants' maturation can be described by means of multifractal and spectral analysis, which show an increasing trend of LF power as well as a decreasing trend of fractal regularity with increasing post-menstrual age. The best obtained performances (R² = 0.68 and MAE = 1.56 weeks) are obtained as combination of fractal and spectral features and are comparable to other developmental models reported by different authors (Hoyer et al., 2013; De Wel et al., 2017; Lavanga et al., 2017, 2018). Second, this predominance of LF and VLF bands as well as the lower scales for the multifractal indices suggest that the sympathovagal balance of neonates might not simply be related to the ratio $\frac{LF}{HF}$, but the entire HRV band and the regularity of the tachogram should be included to have a better understanding of the ANS maturation. Third, bradycardias might forcefully increase the variance of the heart rate and disrupt the relationship between autonomic indices and age, especially for commonly used metrics in clinical practice. The PMA estimation models based on novel HRV indices provide a more comprehensive understanding of post-natal autonomic maturation. They can be also considered an alternative automated maturity index to other electrophysiological data analysis for the NICU. This research might be a first step to design personalized therapies or preventive care to preserve infants' development.

Data Availability Statement

The datasets presented in this article are not readily available because the clinical metadata of patients (such as hospital ID, age, recording time-stamps, pain scales) are subject to the European data-privacy policy. Requests to access the datasets should be directed to bWxhdmFuZ2EmI3gwMDA0MDtlc2F0Lmt1bGV1dmVuLmJl. The authors will try to provide an anonymized version of the dataset in compliance with the privacy policy of the University Hospitals of Leuven, which is the owner of the data.

Ethics Statement

The studies involving human participants were reviewed and approved by Ethical Committee of University Hospitals Leuven, Belgium. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin. Written informed consent was obtained from the individual(s), and minor(s)' legal guardian/next of kin, for the publication of any potentially identifiable images or data included in this article.

Author Contributions

ML wrote the article and conducted the data analysis. EH and JM supported the data analysis. AC and SV supervised the data analysis. KJ and GN conducted the data collection. EH, JM, BB, KJ, EO, GN, SV, and AC reviewed and corrected the manuscript. All authors contributed to the article and approved the submitted version.

Funding

The research was supported by Bijzonder Onderzoeksfonds KU Leuven (BOF) : C24/15/036 The effect of perinatal stress on the later outcome in preterm babies, EU: H2020 MSCA-ITN-2018: INtegrating Functional Assessment measures for Neonatal Safeguard (INFANS), funded by the European Commission under Grant Agreement #813483. This research received funding from the Flemish Government (AI Research Program). SV and ML were affiliated to Leuven.AI-KU Leuven institute for AI, B-3000, Leuven, Belgium. Funders were not involved in the study design, collection, analysis, interpretation of data, the writing of this article or the decision to submit it for publication.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

ML was a SB Ph.D. fellow at Fonds voor Wetenschappelijk Onderzoek (FWO), Vlaanderen, supported by the Flemish government.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2020.581250/full#supplementary-material

References

Abry, P., Wendt, H., Jaffard, S., Helgason, H., Goncalves, P., Pereira, E., et al. (2010). “Methodology for multifractal analysis of heart rate variability: from LF/HF ratio to wavelet leaders,” in 2010 Annual International Conference of the IEEE Engineering in Medicine and Biology (Buenos Aires: IEEE), 106–109. doi: 10.1109/IEMBS.2010.5626124

CrossRef Full Text | Google Scholar

Atkinson, E., and Fenton, A. C. (2009). Management of apnoea and bradycardia in neonates. Paediatr. Child Health 19, 550–554. doi: 10.1016/j.paed.2009.06.008

CrossRef Full Text | Google Scholar

Aylward, G. P. (2014). Neurodevelopmental outcomes of infants born prematurely. J. Dev. Behav. Pediatr. 35, 394–407. doi: 10.1097/01.DBP.0000452240.39511.d4

CrossRef Full Text | Google Scholar

Bolea, J., Pueyo, E., Orini, M., and Bailón, R. (2016). Influence of heart rate in non-linear HRV indices as a sampling rate effect evaluated on supine and standing. Front. Physiol. 7:501. doi: 10.3389/fphys.2016.00501

PubMed Abstract | CrossRef Full Text | Google Scholar

Camm, A., Malik, M., Bigger, J., and Günter, B. (1996). Task force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. Heart rate variability: standards of measurement, physiological interpretation and clinical use. Circulation 93, 1043–1065. doi: 10.1161/01.CIR.93.5.1043

CrossRef Full Text | Google Scholar

Clairambault, J., Curzi-dascalovab, L., Kauffmann, F., and Mkdiguea, C. (1992). Heart rate variability in normal sleeping and preterm neonates full-term. Early Hum. Dev. 28, 169–183. doi: 10.1016/0378-3782(92)90111-S

PubMed Abstract | CrossRef Full Text | Google Scholar

Curzi-Dascalova, L. (1994). Development of sleep and autonomic nervous system contorl in premature and full-term newborns. Archiv. Pediatr. 2, 255–262. doi: 10.1016/0929-693X(96)81138-9

PubMed Abstract | CrossRef Full Text | Google Scholar

David, M., Hirsch, M., Karin, J., Toledo, E., and Akselrod, S. (2007). An estimate of fetal autonomic state by time-frequency analysis of fetal heart rate variability. J. Appl. Physiol. 102, 1057–1064. doi: 10.1152/japplphysiol.00114.2006

CrossRef Full Text | Google Scholar

De Souza Filho, L. F. M., De Oliveira, J. C. M., Ribeiro, M. K. A., Moura, M. C., Fernandes, N. D., De Sousa, R. D., et al. (2019). Evaluation of the autonomic nervous system by analysis of heart rate variability in the preterm infants. BMC Cardiovasc. Disord. 19:198. doi: 10.1186/s12872-019-1166-4

PubMed Abstract | CrossRef Full Text | Google Scholar

De Wel, O., Lavanga, M., Dorado, A., Jansen, K., Dereymaeker, A., Naulaers, G., et al. (2017). Complexity analysis of neonatal EEG using multiscale entropy: applications in brain maturation and sleep stage classification. Entropy 19:516. doi: 10.3390/e19100516

CrossRef Full Text | Google Scholar

Dereymaeker, A., Pillay, K., Vervisch, J., Van Huffel, S., Naulaers, G., Jansen, K., et al. (2017). An automated quiet sleep detection approach in preterm infants as a gateway to assess brain maturation. Int. J. Neural Syst. 27:1750023. doi: 10.1142/S012906571750023X

PubMed Abstract | CrossRef Full Text | Google Scholar

Doret, M., Spilka, J., Chudáček, V., Gonçalves, P., Abry, P., and Stanley, H. (2015). Fractal analysis and hurst parameter for intrapartum fetal heart rate variability analysis: a versatile alternative to frequency bands and LF/HF ratio. PLoS ONE 10:e0136661. doi: 10.1371/journal.pone.0136661

PubMed Abstract | CrossRef Full Text | Google Scholar

Franke, K., Luders, E., May, A., Wilke, M., and Gaser, C. (2012). Brain maturation: predicting individual BrainAGE in children and adolescents using structural MRI. Neuroimage 63, 1305–1312. doi: 10.1016/j.neuroimage.2012.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Gee, A. H., Barbieri, R., Paydarfar, D., and Indic, P. (2017). Predicting bradycardia in preterm infants using point process analysis of heart rate. IEEE Trans. Biomed. Eng. 64, 2300–2308. doi: 10.1109/TBME.2016.2632746

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoyer, D., Tetschke, F., Jaekel, S., Nowack, S., Witte, O. W., Schleußner, E., et al. (2013). Fetal functional brain age assessed from universal developmental indices obtained from neuro-vegetative activity patterns. PLoS ONE 8:e74431. doi: 10.1371/journal.pone.0074431

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C., Amaral, L. a., Goldberger, a. L., Havlin, S., Rosenblum, M. G., Struzik, Z. R., et al. (1999). Multifractality in human heartbeat dynamics. Nature 399, 461–465. doi: 10.1038/20924

PubMed Abstract | CrossRef Full Text | Google Scholar

Javorka, K., Lehotska, Z., Kozar, M., Uhrikova, Z., Kolarovszki, B., Javorka, M., et al. (2017). Heart rate variability in newborns. Physiol. Res. 66, 203–214. doi: 10.33549/physiolres.933676

CrossRef Full Text | Google Scholar

Karin, J., Hirsch, M., and Akselrod, S. (1993). An estimate of fetal autonomic state by spectral analysis of fetal heart rate fluctuations. Pediatr. Res. 34, 134–138. doi: 10.1203/00006450-199308000-00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Koolen, N., Dereymaeker, A., Räsänen, O., Jansen, K., Vervisch, J., Matic, V., et al. (2016). Early development of synchrony in cortical activations in the human. Neuroscience 322, 298–307. doi: 10.1016/j.neuroscience.2016.02.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Krueger, C., van Oostrom, J. H., and Shuster, J. (2010). A longitudinal description of heart rate variability in 28–34-week-old preterm infants. Biol. Res. Nurs. 11, 261–268. doi: 10.1177/1099800409341175

PubMed Abstract | CrossRef Full Text | Google Scholar

Lavanga, M., De Wel, O., Caicedo, A., Jansen, K., Dereymaeker, A., Naulaers, G., et al. (2017). Monitoring effective connectivity in the preterm brain: a graph approach to study maturation. Complexity 2017:9078541. doi: 10.1155/2017/9078541

CrossRef Full Text | Google Scholar

Lavanga, M., Wel, O. D., Caicedo, A., Jansen, K., Dereymaeker, A., Naulaers, G., et al. (2018). A brain-age model for preterm infants based on functional connectivity. Physiol. Meas. 39:044006. doi: 10.1088/1361-6579/aabac4

PubMed Abstract | CrossRef Full Text | Google Scholar

Leonarduzzi, R. F., Schlotthauer, G., and Torres, M. E. (2010). “Wavelet leader based multifractal analysis of heart rate variability during myocardial ischaemia,” in 2010 Annual International Conference of the IEEE Engineering in Medicine and Biology (Buenos Aires), 110–113. doi: 10.1109/IEMBS.2010.5626091

PubMed Abstract | CrossRef Full Text | Google Scholar

Longin, E., Gerstner, T., Schaible, T., Lenz, T., and König, S. (2006). Maturation of the autonomic nervous system: differences in heart rate variability in premature vs. term infants. J. Perinatal Med. 34, 303–308. doi: 10.1515/JPM.2006.058

PubMed Abstract | CrossRef Full Text | Google Scholar

Mazursky, J. E., Birkett, C. L., Bedell, K. A., Ben-Haim, S. A., and Segar, J. L. (1998). Development of baroreflex influences on heart rate variability in preterm infants. Early Hum. Dev. 53, 37–52. doi: 10.1016/S0378-3782(98)00038-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Moeyersons, J., Amoni, M., Huffel, S. V., Willems, R., and Varon, C. (2019). R-DECO: an open-source Matlab based graphical user interface for the detection and correction of R-peaks. bioRxiv 560706. doi: 10.7717/peerj-cs.226

CrossRef Full Text | Google Scholar

Montalto, A., Faes, L., Marinazzo, D., Porta, A., and Vicente, R. (2014). MuTE: a MATLAB toolbox to compare established and novel estimators of the multivariate transfer entropy. PLoS ONE 9:e109462. doi: 10.1371/journal.pone.0109462

PubMed Abstract | CrossRef Full Text | Google Scholar

Mulkey, S. B., Kota, S., Swisher, C. B., Hitchings, L., Metzler, M., Wang, Y., et al. (2018). Autonomic nervous system depression at term in neurologically normal premature infants. Early Hum. Dev. 123, 11–16. doi: 10.1016/j.earlhumdev.2018.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Orini, M., Bailon, R., Mainardi, L. T., Laguna, P., and Flandrin, P. (2012). Characterization of dynamic interactions between cardiovascular signals by time-frequency coherence. IEEE Trans. Biomed. Eng. 59, 663–673. doi: 10.1109/TBME.2011.2171959

PubMed Abstract | CrossRef Full Text | Google Scholar

Paolillo, P., and Picone, S. (2013). Apnea of prematurity. J. Pediatr. Neonatal Individ. Med. 2:e020213. doi: 10.7363/020213

CrossRef Full Text | Google Scholar

Pereira, C. B., Heimann, K., Venema, B., Blazek, V., Czaplik, M., and Leonhardt, S. (2017). “Estimation of respiratory rate from thermal videos of preterm infants,” in 2017 39th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) (Seogwipo: Institute of Electrical and Electronics Engineers Inc.), 3818–3821. doi: 10.1109/EMBC.2017.8037689

CrossRef Full Text | Google Scholar

Poets, C. F., Stebbens, V. A., Samuels, M. P., and Southall, D. P. (1993). The relationship between bradycardia, apnea, and hypoxemia in preterm infants. Pediatr. Res. 34, 144–147. doi: 10.1203/00006450-199308000-00007

PubMed Abstract | CrossRef Full Text | Google Scholar

Popivanov, D., Stomonyakov, V., Minchev, Z., Jivkova, S., Dojnov, P., Jivkov, S., et al. (2006). Multifractality of decomposed EEG during imaginary and real visual-motor tracking. Biol. Cybern. 94, 149–156. doi: 10.1007/s00422-005-0037-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Porges, S. W. (1995). Orienting in a defensive world: Mammalian modifications of our evolutionary heritage. A polyvagal theory. Psychophysiology 32, 301–318. doi: 10.1111/j.1469-8986.1995.tb01213.x

CrossRef Full Text | Google Scholar

Porta, A., Bari, V., Marchi, A., De Maria, B., Cysarz, D., Van Leeuwen, P., et al. (2015). Complexity analyses show two distinct types of nonlinear dynamics in short heart period variability recordings. Front. Physiol. 6:71. doi: 10.3389/fphys.2015.00071

PubMed Abstract | CrossRef Full Text | Google Scholar

Porta, A., Baselli, G., and Cerutti, S. (2006). Implicit and explicit model-based signal processing for the analysis of short-term cardiovascular interactions. Proc. IEEE 94, 805–818. doi: 10.1109/JPROC.2006.871774

CrossRef Full Text | Google Scholar

Porta, A., Casali, K. R., Casali, A. G., Gnecchi-Ruscone, T., Tobaldini, E., Montano, N., et al. (2008). Temporal asymmetries of short-term heart period variability are linked to autonomic regulation. Am. J. Physiol. Regul. Integr. Compar. Physiol. 295, R550–R557. doi: 10.1152/ajpregu.00129.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Porta, A., De Maria, B., Bari, V., Marchi, A., and Faes, L. (2017). Are nonlinear model-free conditional entropy approaches for the assessment of cardiac control complexity superior to the linear model-based one? IEEE Trans. Biomed. Eng. 64, 1287–1296. doi: 10.1109/TBME.2016.2600160

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, S. L., Haley, S., Slater, H., and Moyer-Mileur, L. J. (2013). Heart rate variability during caregiving and sleep after massage therapy in preterm infants. Early Hum. Dev. 89, 525–529. doi: 10.1016/j.earlhumdev.2013.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Søvik, S., Lossius, K., and Walløe, L. (2001). Heart rate response to transient chemoreceptor stimulation in term infants is modified by exposure to maternal smoking. Pediatr. Res. 49, 558–565. doi: 10.1203/00006450-200104000-00019

PubMed Abstract | CrossRef Full Text | Google Scholar

Wendt, H., Abry, P., and Jaffard, S. (2007). Bootstrap for multifractal analysis. IEEE Signal Process. Mag. 24, 38–48. doi: 10.1109/MSP.2007.4286563

CrossRef Full Text | Google Scholar

Widjaja, D., Orini, M., Vlemincx, E., and Huffel, S. V. (2013). Cardiorespiratory dynamic response to mental stress: a multivariate time-frequency analysis. Comput. Math. Methods Med. 2013:451857. doi: 10.1155/2013/451857

PubMed Abstract | CrossRef Full Text | Google Scholar