Cerebral hemodynamic responses in preterm-born neonates to visual stimulation: classification according to subgroups and analysis of frontotemporal–occipital functional connectivity

Abstract. How neurovascular coupling develops in preterm-born neonates has been largely neglected in scientific research. We measured visually evoked (flicker light) hemodynamic responses (HRs) in preterm-born neonates (n=25, gestational age: 31.71±3.37 weeks, postnatal age: 25.48±23.94 days) at the visual cortex (VC) and left frontotemporal lobe (FTL) using functional near-infrared spectroscopy (fNIRS) neuroimaging. We found that the HR characteristics show a large intersubject variability but could be classified into three groups according to the changes of oxyhemoglobin concentration at the VC [(A) increase, (B) decrease, or (C) inconclusive]. In groups A and B, the HRs at the left FTL were correlated with those at the VC, indicating the presence of a frontotemporal–occipital functional connectivity. Neonates in group A had a higher weight at measurement compared to those in group B, and had the lowest baseline total hemoglobin concentration and hematocrit compared to group C. To the best of our knowledge, this is the first fNIRS study showing (1) that the HRs of preterm-born neonates need to be classified into subgroups, (2) that the subgroups differed in terms of weight at measurement, and (3) that HRs can be observed also at the FTL during visual stimulation. These findings add insights into how neurovascular coupling develops in preterm-born neonates.


Introduction
Brain activity is associated with physiological changes, which alter the optical properties of tissue. These optical changes can be detected by functional near-infrared spectroscopy (fNIRS) neuroimaging. 1 This technique is well suited for the investigation of fragile newborn infants because it enables noninvasive measurements, is portable, and delivers critical physiological information about cerebral tissue oxygenation and hemodynamics. Neuroimaging with fNIRS has the potential to provide insights into the origin of neonatal brain lesions, as well as into functional development of both the normal and the abnormal brain. 2,3 The ability of fNIRS to access neurovascular coupling associated with neuronal activity offers a useful opportunity to investigate, for example, the neural basis of visual cognitive development in infants. This is of special interest in preterm-born infants, as many of them develop visual cognitive impairment in later childhood. 4 This visual impairment is thought to result from injury to the posterior visual pathway involving periventricular white matter lesions. 5 The visual impairment can develop even in the absence of major neuromotor impairments. 6 The visual function deficits seen in children born prematurely may be related to the networks involving the cortical dorsal stream and its connection to the parietal, frontal, and hippocampal areas. 7 Even if preterm-born infants do not present with major brain injuries, higher brain dysfunctions may appear during development. There are reasons to conclude that earlier exposure to the extrauterine environment influences the developmental trajectory of the brain. [8][9][10] From a medical point of view, higher brain dysfunction, such as aberration in the visual pathway, should be identified in the early stages of development to enable early neuroprotective intervention. The assessment of the hemodynamic response (HR) to cognitive processing could reveal some aspects of physiological traits of these higher brain functions.
A growing number of functional magnetic resonance imaging (fMRI) studies have shown that preterm-born neonates and infants have brain structures with aberrant volumes, morphologies, and networks at term equivalent age. 11 Given the fact that preterm birth before 37 weeks of gestation accounts for more than 10% of all live births, there is a clear need for a practical and safe method to detect the status of brain development in preterm-born neonates. 12 Neuroimaging with fNIRS has been shown to be a practical and powerful tool to examine the status of hemodynamic physiology in infants as well as in adults. 2,13,14 Neuronal activity is associated with changes in tissue metabolism (neurometabolic coupling) and tissue hemodynamics (neurovascular coupling). These changes lead in turn to characteristic changes in the oxygenated (oxyhemoglobin, O 2 Hb) and deoxygenated (deoxyhemoglobin, HHb) states of hemoglobin: a distinct increase in the concentration in [O 2 Hb] and a slight decrease in [HHb]. 1 The change in hemodynamics (blood flow/volume) follows a specific pattern (functional hyperemia) known as the hemodynamic response function (HRF). 15,16 Neuroimaging with fNIRS relies on the measurement of changes in [O 2 Hb] and [HHb], which can be linked to neural activity (provided that other factors influencing the measured fNIRS signals can be excluded 17 ).
In previous fNIRS studies by our group, we demonstrated the ability to measure evoked hemodynamic changes in the occipital cortex after visual stimulation in term-born neonates within their first days of life. [18][19][20] However, the development of the HRF pattern in the human brain remained unclear; some fNIRS studies reported an inverted HRF pattern (i.e., a decrease in blood flow or [O 2 Hb]) in young infants in response to perceptual stimuli, whereas some studies did not. 2,21,22 This has been a controversial issue, and investigations about this aspect in preterm-born neonates are particularly sparse, partly because the investigation of HRFs of these neonates using fNIRS is a relatively new subject. [22][23][24] Notably, a recent fNIRS study 23 reported developmental changes in phase differences of spontaneous low-frequency oscillations of [O 2 Hb] and [HHb] in pretermborn infants. Although such developmental changes in neurovascular regulation are very intriguing, these results come only from resting-state (rs) measurements.
Relatively little and controversial information exists on the HR to visual stimulation in the visual cortex (VC) of the preterm brain. The development of the dynamic properties of hemodynamics and metabolism due to task-or stimulus-evoked changes in brain activity in preterm-born neonates has yet to be investigated in more detail. Moreover, functional brain connectivity has been only rarely accessed with fNIRS in preterm-born neonates.
Therefore, this study was undertaken specifically to elucidate the functional developmental process in the preterm brain with reference to the effects of visual stimulation. In our study, we specifically examined two aspects of the brain response in premature infants: (1) the HR pattern in the VC and (2) the frontotemporal-occipital connectivity during visual stimulation and characterization of the child-dependent characteristics of it.

Participants
The study was performed at the Department of Neonatology, University Hospital Zurich. The study protocol was accepted by the ethical committee of Zurich (KEK 2010-0102/2) and Swissmedic (2010-MD-0019). Clinically stable preterm-born neonates (n ¼ 36) breathing spontaneously on room air were enrolled in our study. Infants with congenital malformations were not recruited for the study. Parental consent was obtained in all cases prior to enrolment.
Of the 36 preterm-born neonates measured, 11 data sets (i.e., 30.5%) had to be excluded from further analysis due to insufficient data quality: 4 of the 11 data sets contained too frequent severe movement artifacts and in 7 data sets high noise was present caused by technical interference.
The specific exclusion of data had to be done manually based on our experience with fNIRS data from the past since the complexity of the signal quality assessment could not be automated sufficiently. The final data set used for the analysis therefore comprised measurements of 25 preterm-born neonates. The demographic data as well as baseline physiological data for this population are shown in Table 1.

Study Protocol
Preterm-born neonates selected for the measurements were brought into a quiet examination room near the neonatal intensive care unit. After setting up the equipment, the physiological measurements were started as soon as the neonates were in a calm state (no crying, no large movements, being asleep). Before the fNIRS measurement of the stimulus-evoked HRs, baseline physiological values were measured with fNIRS, as described in a previous validation study of the new NIRS device. 25 The total measurement lasted about 15 min. The measurements were conducted by two or three persons, enabling one person to respond quickly to the behavior of the neonate while the other(s) focused on performing the measurements. Behavioral states of the newborn infant were closely monitored at the bedside, according to the criteria of Prechtl. 26

Visual Stimulation
Visual stimulation was elicited by flicker light stimulation with a portable handheld 7″ TFT LCD monitor (667GL-70NP/H/Y/S, Lilliput Ltd.; brightness: 450 cd∕m 2 , contrast: 500:1). The light Table 1 Demographic data and baseline physiological data for the whole population of analyzed neonates (n ¼ 25). The numerical data are given as mean AE standard deviation (SD), median (first, third quartile).

Parameter
All (n ¼ 25) It implements a self-calibrating multidistance (SCMD) measurement technique, i.e., a combination of the spatially resolved spectroscopy 28,29 and the self-calibration method. 30 We term this technique the SCMD method. SCMD reduces artifacts caused by (1) differences of light coupling between the sensors and the detectors, (2) influences from blood perfusion changes in the superficial part of the tissue, and (3) movement of the tissue with respect to the sensor causing movement artifacts. 31 As depicted in Fig. 1(e), the source-detector separations (SDS) of light channels of the optode have different lengths (15, 20, 30, and 35 mm), making it possible to perform depth-resolved NIRS measurements. For a more detailed listing of the sensor's features, see Ref. 25.
The raw optical data recorded by the NIRS devices were converted into relative concentration values using the SCMD method. This involved calculating the slope of intensity values versus distance. Based on this and the known absorption spectrum of [O 2 Hb] and [HHb], the relative (to a baseline value) concentrations of [O 2 Hb] and [HHb] were calculated (denoted as Δ½O 2 Hb and Δ½HHbÞ. In addition, Δ½O 2 Hb and Δ½HHb were also determined using the traditional single-distance approach (for all available source-detector combination and thus SDS values), employing the modified Beer-Lambert law (MBLL). 1 From the Δ½O 2 Hb and Δ½HHb signals, two additional signals were calculated: absolute tissue oxygen saturation (StO 2 ) and relative total hemoglobin concentration (Δ½tHb ¼ Δ½O 2 Hb þ Δ½HHb). The device also measured the absolute hemoglobin concentration [tHb], assuming no individual differences between neonates in the optical light-scattering properties of head tissue. Since the validity of this assumption is not known, [tHb] measured by OxyPrem v1.3 has to be used in a proper way taking into account this aspect.
Pulse oximetry (Sensmart Model X-100 Universal Oximetry System, 6100CN cloth sensor-neonatal pulse oximetry sensor; Nonin Medical, Inc., Plymouth, Minnesota) was used to measure the peripheral arterial oxygen hemoglobin saturation (SpO 2 ) during the whole experiment. The SpO 2 sensor was placed at the right hand of the neonate [see Fig. 1 All data were available with a sampling frequency of 5 Hz.
1. Signals were visually inspected for the presence of technical noise in form of strong periodic fluctuations in the region of about 1 Hz. This type of noise had been observed already by a previous analysis of measurements with our NIRS device. 25 In about 20% of the data sets, such a noise component was observed (a sinusoidal signal with a frequency of about 1 Hz and an amplitude modulation with a frequency of about 0.05 Hz). Data sets with this strong interference signal had to be completely removed from the final data sets used for the next signal-processing steps. The origin of the interference signal could not be determined exactly. The authors suspect an influence of technical electromagnetic fields present in the measurement room, which could affect the measurement device under specific conditions.
2. Movement artifacts in the data (mainly evident as spikes and baseline shifts in the data) were identified and corrected using a semiautomatic approach based on moving-standard deviation and spline interpolation developed by our group. 32 For the present study, further development of this method (replacing the spline interpolation with the shape-preserving piecewise cubic interpolation) has been used. Each data set (n ¼ 25) available had to be processed by this method due to the inability to avoid movements of the neonates during the measurement.
3. To remove high-frequency noise, improve the signalto-noise ratio (SNR) and to enhance the spectral band with the physiologically relevant information, all NIRS data were band-pass filtered (lower cut-off frequency: 0.0167 Hz ¼ 60 s, upper cut-off frequency: 0.2 Hz ¼ 5 s, phase-corrected finite-impulse response filter of order 500; the high order of the filter was chosen according to a recent investigation pointing out the necessity to use a high-order filter for fNIRS signal analysis 33 ) and afterward low-pass filtered with a Savitzky-Golay filter to further reduce noise (fifth order, window length = 4 s). The selection of the filter type and filter parameter values was done empirically to obtain the optimal SNR. The filter specifications were based on the fact that the stimulation was 20-s long and the interstimulus interval (ISI) was 35.8 AE 3.4 s (mean AE SD). In addition, it proved to be important to set the lower cut-off frequency as high as possible to remove strong low-frequency physiological oscillations and to obtain an optimal SNR. In practice, a lower cut-off of 0.0167 Hz was chosen, Neurophotonics 045005-3 Oct-Dec 2019 • Vol. 6(4) even though it is higher than the first harmonic of one of the 15 repetition, because this removes lowfrequency interference in the data of all 15 repetitions and hence leads to a higher overall SNR.
4. Owing to the often strong movements of the neonates or problems with the attachment of the sensors to the skin, several data sets processed with steps 1 and 2 still contained artifacts. Consequently, all data points were  removed when they exceeded a specific threshold value (T) given as T ¼ α × median [moving standard deviation (MSD)], with MSD of the whole time-series (with a window length of 20 samples i.e., 4 s) and a scaling factor that was empirically chosen (e.g., α ¼ 2.5) for each data set to ensure that no data from artifacts were present in the cleaned data. This signalprocessing step was essential to ensure that only data of sufficient quality were used in the further data analysis steps.

Block averaging
To detect stimulus-evoked HRs in the NIRS data, block averages were calculated for each signal, measurement location, and neonate. To this end, the median over the 15 single trials was calculated as well as the standard error of the median (SEM), according to SEM ¼ ð0.741 þ IQRÞ∕sqrtðnÞ, where IQR is the interquartile range and n is the number of trials (n ¼ 15).

Statistical analysis
Whether the stimulus-evoked hemodynamic changes of the fNIRS signals from each neonate were statistically significant or not was evaluated by the following approach: the maximum and minimum values during the stimulation (20 s) of the median-based block average were determined (x max ; x min ). All values from all 15 trials that centered around x max and x min (AE2 samples, total length of vector: 5 samples) were then collected to test if they have a distribution that is different from one, with a median value of zero (Wilcoxon signed-rank test). This approach was used since we noticed that time to peak of the hemodynamic responses of the neonates showed a significant intersubject-variability. The method we used in our study (based on detecting the x max and x min values of the individual shape of the neonate's hemodynamic response) ensured to take this intersubject variability into account for the statistical analysis. For each data set (each neonate), 16 p-values were then obtained [4 (Δ½O 2 Hb, Δ½HHb, Δ½tHb, StO 2 ) × 2 (x max ; x min ) × 2 (left FTL, VC)]. A false discovery rate correction was then applied to the p-values in order to correct the multiple comparison situation. The corrected p-values were then used for further evaluation of the statistical significance of the single-subject HRs (first-level analysis). Group averages of the HRs were also calculated and their statistical significance was also evaluated using the same approach (second-level analysis). The statistical analysis of the HRs was performed in Matlab. Further statistical analysis has been conducted with JASP (v.0.9.2, The JASP Team, Amsterdam, The Netherlands).
For the statistical analysis of possible differences in demographic and baseline physiological parameters between groups of neonates, a Kruskal-Wallis test (one-way ANOVA on ranks) has been used. Post-hoc testing was performed with Dunn's post-hoc test, correcting for the multiple comparison situation with the Holm-Bonferroni method.
Possible relationships between the magnitude of HRs and demographic and baseline physiological parameters have been analyzed with linear correlation analysis.
Figures were created using Matlab, R (version 3.5.2, R Foundation for Statistical Computing, Vienna, Austria), Adobe Illustrator (version CS6, Adobe systems software, Ireland), and Adobe Photoshop (version CS6, Adobe systems software, Ireland).

Stimulus-Evoked Hemodynamic Responses
Can Be Measured in Preterm-Born Neonates Both at the Visual Cortex and the Frontotemporal Lobe As already described in Sec. 2, data sets of 25 neonates contained data of sufficient quality for a reliable statistical analysis. Stimulus-evoked hemodynamic changes of different magnitude and waveform were observed at the VC as well as at the left FTL. A significant HR (positive or negative) at the VC was detected in 16 of the 25 neonates (64%).

Stimulus-Evoked Hemodynamic Responses Show a Large Intersubject Variability
We noticed high intersubject variability in the HRs in Δ½O 2 Hb, Δ½HHb, Δ½tHb, and StO 2 with respect to the magnitude and waveform. The specific form of the HRF was thus quite different between the neonates. Figure 2 shows the evoked changes measured in six neonates exemplarily (the subgrouping used in this figure is explained in the next section).

Magnitude and Waveform of the Stimulus-Evoked Hemodynamic Responses Can Be Classified into Three Groups
Looking at the HRs of the individual neonates, we came to the conclusion that in general three pattern types can be observed for the time during the stimulation: a decrease, increase, or an inconclusive change. Furthermore, we noticed that the HR of Δ½O 2 Hb at the VC were mostly pronounced and had clear features that could be used to classify them. Based on these insights, we classified the 25 data sets according to the sign and statistical significance of Δ½O 2 Hb changes at the VC. This resulted in three groups: (1) statistically significant increase of Δ½O 2 Hb at the VC (group A), (2) statistically significant decrease of Δ½O 2 Hb at the VC (group B), and (3) no statistically significant change of Δ½O 2 Hb at the VC (group C). Based on this assignment, group averages of the single subject-specific HRs were calculated (second-level analysis) ( Fig. 3) and the statistical significance of the group averages was determined. The validity of this classification was evaluated by testing if a machine-learning approach comes to the same conclusion. Therefore, the block-averages of Δ½O 2 Hb of all 25 neonates were analyzed by k-means clustering (assigned number of clusters: k ¼ 3). The analysis revealed that groups A and B were classified by k-means clustering in a very similar way to the semimanual approach (Fig. 4). Overall, the classification accuracy (evaluated by the Rand index R) was 0.69.
Individual HRs, classified into the three groups, can be seen in Fig. 2.
We observed high similarity in the changes in Δ½O 2 Hb and Δ½tHb for all measurement locations (left FTL, VL) and groups (A to C). The shapes of the changes in Δ½O 2 Hb and StO 2 at the VC for groups A and B differed to a specific degree, especially in the poststimulation recovery phase.
The changes in Δ½HHb at the VC and left FTL in groups A and B were not statistically significant, compared to the baseline due to the high intersubject variability of the Δ½HHb responses.
Neurophotonics 045005-5 Oct-Dec 2019 • Vol. 6(4) In group C, the shapes of the HRs were different compared to those in groups A and B, confirming our approach of classifying the HRs into three groups. Significant stimulus-evoked changes in Δ½O 2 Hb, Δ½HHb, Δ½tHb, and StO 2 were observed in group C, but the magnitude of the changes was in general significantly weaker than for groups A and B.

Subgrouping of Hemodynamic Responses
Is Associated with the Weight and Height of the Neonate The subgroup effect raised the question whether the demographic and baseline physiological data of the neonates were also different for these groups. Table 2 lists the data according to the subgroups. The statistical analysis (results see Table 3) revealed that the weight of the neonate at measurement and the Hct were statistically significantly different for the three groups [weight: According to the post-hoc analysis, those neonates with a positive HR at the VC (group A) were heavier at measurement than those of groups B and C (A versus B: p ¼ 0.048, A versus C: p ¼ 0.024) and had a larger body length than those of group B (p ¼ 0.026) [see also Fig. 5(a)]. In addition, the Hct and hemoglobin concentration (cHb) values of the neonates were statistically significantly different (p < 0.05) between groups A and C; neonates with high Hct and cHb values had a higher chance of being classified into group C (no statistically significant HR detected). No other demographic or baseline physiological data were different for the three groups. An overview of all data distributions depending on the groups is shown in Fig. 5. We also investigated if the subgroup effect appears when creating scatterplots showing the Hct versus weight at measurement, cHb versus weight at measurement, as well as gestational age (GA) versus postmenstrual age (PMA) [ Fig. 5(b)]. Data from group A can be seen to cluster in regions of high body weight at measurement and low Hct as well as low cHb. That group A had a tendency to have lower Hct and cHb values can also be seen in the respective subplots in Fig. 5(a).
To further investigate a possible dependence between the sign of the HR and the individual Hct and cHb values, we plotted the magnitude of the HR at the VC (determined as the median of the maximum value of each single HR during the time of visual stimulation) with respect to several baseline parameters (measured before the visual stimulation), i.e., total hemoglobin measured at the VC with NIRS ([tHb]), total hemoglobin  Fig. 6. Linear correlation analysis revealed no statistically significant correlation between the magnitude of the HRs and the seven parameters.

Hemodynamic Responses Show a Frontotemporal-Occipital Functional Connectivity
We were surprised to find that the shape of the HR of Δ½O 2 Hb at the left FTL was dependent on that of the VC: in group A at the VC the canonical HR shape was evident (increase in Δ½O 2 Hb during the stimulation with a decrease after reaching the peak) [ Fig. 3(b)], whereas at the left FTL, Δ½O 2 Hb showed a double oscillation with a fast increase having a maximum at about 15 s after stimulation onset and a minimum at about 30 s [ Fig. 3(a)]. Fascinatingly, the patterns in the left FTL was reversed comparing group A with group B-in accordance with the shape of the HR at the VC, which also reversed the sign comparing group A with group B. We, therefore, concluded the presence of a form of frontotemporal-occipital functional connectivity (fto-FC). The fto-FC is also evident in the Δ½tHb data but not in StO 2 and Δ½HHb.
For a better visualization of the fto-FC, we plotted the groupaveraged stimulus-evoked time course of changes in Δ½O 2 Hb, Δ½HHb, Δ½tHb, and StO 2 at the left FTL with respect to those at the VC (Fig. 7). This allows easy identification of nonlinear correlations between time series that have different oscillation frequencies; a circle will result if one of the time series has a frequency that is twice that of the frequency of another. Such circles clearly are visible in Fig. 7, indicating the presence of a specific fto-FC. In addition, color-coding the time information in the scatter plots shows visually that the circles of group A have a right turn (counterclockwise rotation), whereas those of group B show a left one (clockwise rotation). This specific feature agrees with the shape changes (see in Fig. 3).

Hemodynamic Responses Are Present in the Cerebral and Extracerebral Tissue Layer
Having the ability to measure with a multidistance NIRS oximeter allowed us to investigate possible differences in the shape of the HR depending on the position measured in the tissue. In addition to the fNIRS signals determined by the multidistance self-calibrating measurement technique, the fNIRS signals were also calculated with the traditional single-distance MBLL approach involving an SDS of either 1.5 or 3.5 cm. All three methods have a different sensitivity to tissue regions: Method 1 (SCMD) is most sensitive to deeper tissue layers and reduces the impact of hemodynamic changes happening in superficial ones. Method 2 (MBLL with SDS ¼ 1; 5 cm) is most sensitive to extracerebral hemodynamics in the skin. Method 3 (MBLL with SDS ¼ 3.5 cm) probes hemodynamics in deeper layers as well as in the extracerebral compartment. Figure 8 shows the HRs from five randomly chosen neonates at both measurement locations calculated with all three methods. Three things can be noticed: (1) there are hemodynamic changes happening in the extracerebral part of the head; (2) the hemodynamic changes in the extracerebral compartment are generally larger than those in the cerebral one in case of the left FTL measurements, but the opposite is generally true for the VC; and (3) the shape of the HR (in Δ½O 2 Hb) is similar for all three methods, but method 1 (SCMD) is able to remove disturbing extracerebral hemodynamic changes, resulting in the most reliable approach to detect hemodynamic changes in brain tissue due to neurovascular coupling.

Detection Success Rate of Stimulus-Evoked Hemodynamic Responses in Preterm-Born Neonates
In our study, we had a sensitivity rate of 64% of detecting HRs at the VC. This value is slightly higher than we found in previous studies in term-born infants (50% 18 and 61% 20 ). These two and other fNIRS studies involving neonates generally used only group averages (second-level analysis) to calculate the HRs and did not perform a single-subject analysis. 3 The information about their sensitivity for single-subject detection of HR is missing. One exception is the study of Zimmermann et al. 34 In their study using auditory stimulation, the authors reported a sensitivity of 100% in a small sample (n ¼ 6) of preterm-born neonates. However, an increase in Δ½O 2 Hb has only been found in 7 of 18 cases, whereas the remaining 11 measurements showed a decreased Δ½HHb during the auditory stimulation. An absent HR could be due to multiple factors such as (1) technical ones (e.g., lack of sensitivity of the NIRS devices to detect small changes in hemodynamics with a good SNR, measuring at the wrong position where no neuronal-caused hemodynamic change happened, artifacts in the measured signals due to changes in the sensor attachment or influences on the electronics from environmental factors, suboptimal experimental paradigm, different fNIRS signal measurement methodologies, and signal-processing methods used) or (2) physiological ones (e.g., immature neurovascular coupling, different states of anatomical/vascular development, different states of consciousness). 1,3 Therefore, we made sure to optimize the technical aspects to obtain and analyze fNIRS signals of good quality and to perform a sophisticated data analysis; however, we had to discard 11 data sets from the final analysis due to insufficient signal quality. Since for the final data analysis only fNIRS signals of sufficient quality were used, we were able to exclude technical reasons as the cause of our inability to measure HRs in 36% of our neonates. Physiological factors appear, therefore, to be the more obvious reason.

Hemodynamic Responses in Preterm-Born
Neonates: Why a Subgroup Analysis Is Always Needed One of the main findings of our study is that the shapes of the HRs can be classified into three subgroups: increase, decrease, or no change in [O 2 Hb] at the VC. Interestingly, the changes in Δ½O 2 Hb observed in both brain regions (VC and the left FTL) showed a slightly opposite direction. In group A, we found a clear increase in Δ½O 2 Hb, whereas in the left FTL, a slight increase in Δ½O 2 Hb was followed by a decrease in Δ½O 2 Hb.
In group B, a clear decrease in Δ½O 2 Hb at the VC was evident, whereas in the left FTL, first a decrease followed by an increase in Δ½O 2 Hb was found. This highlights once more the question of why we see such differences in the hemodynamic changes in our group of preterm-born neonates, and whether it will generally be possible to define a typical HR for different regions in the cortex of the developing brain. Previously published fNIRS studies investigating HRs of preterm-born neonates [35][36][37][38][39][40][41][42][43] did not perform a subgroup analysis.  Instead, a group average was calculated. Since it is most likely that in these studies also neonates showed HRs of different signs, the calculated group averages (ignoring the individual sign of the HRs) should be interpreted with care.

What Is the Reason for the Differences in Hemodynamic Responses in Preterm-Born Neonates?
The measurements in our study were performed on the neonates during their spontaneous sleep; the sleeping stages were visually determined by the behavioral state of the newborn according to Prechtl. 26 Previous studies showed that visual stimulation during sleep produced increases in Δ½O 2 Hb at the occipital cortex in newborns 18,44-46 and preterm-born neonates investigated at termcorrected age. 47 Only the changes in Δ½HHb showed a large variability in these studies. In a recent study by Taga et al., 24 the same visual stimulus produced different responses in different age groups within the first year of life [see Fig. 6(i)]. While the youngest group of neonates (2 to 3 months of age) showed a pattern of increase in Δ½O 2 Hb and decrease in Δ½HHb, similar to the other studies in newborns, the oldest age group (10 to 11 months of age) showed an opposite change with a decrease in Δ½O 2 Hb and an increase in Δ½HHb in the VC. The authors concluded from their findings that visual processing undergoes development of inhibitory mechanisms during sleep and these phenomena primarily reflect neural development rather than vascular development. 24 In our study, we found both types of responses, such as the one seen in term newborns and up to 3 months of age, but we also observed changes such as those normally present in older infants aged 10 to 11 months 24 [see Fig. 6(i)]. Therefore, we speculate that some of the neonates who showed also a decrease in Δ½O 2 Hb have already developed some inhibitory mechanisms as seen in older ones. This may be explained by earlier exposure to the extrauterine environment in the neonatal intensive care unit. Thus, it is unlikely that the developmental state is the cause of the variation observed. In addition, spontaneous sleep has an effect on the shape and latency of visually evoked potentials in infants, 48,49 and sleep stages influence the cerebral cHb. 50 Brain maturity in terms of both angiogenesis and synaptogenesis in preterm-born and term-born neonates is a possible factor influencing the HR pattern. A recent fNIRS study examining the phase relationship of spontaneous oscillations of Δ½O 2 Hb and Δ½HHb 23 further revealed rapid changes in neurovascular regulation from neonates with a PMA of 34 to term, which was similarly seen in response to phonetic changes of speech in preterm-born and term-born infants. 22 Watanabe et al. 23 showed that the phase relationship of spontaneous oscillations of Δ½O 2 Hb and Δ½HHb in very preterm-born infants was lower than term-born infants and had a slower development compared with late preterm-born and term-born infants at later PNA, suggesting that (1) early preterm-born infants may have an advanced development of their circulatory process after birth and that (2) the structural and/or functional aspects of the development of the neurovascular system may be altered during the first month of life.
So far, the physiological mechanism underlying the variable response of preterm-born neonates remains unclear. However, it could possibly be explained by individual immature neurovascular and/or metabolic functional ability in the cerebral cortex. Specifically, immature functioning of the synaptic structure with less myelination might result in inefficient energy use requiring more [O 2 Hb]. 51 Furthermore, due to the immaturity of the arteriole vessels and capillaries, blood flow may not respond sufficiently to brain activity in preterm-born neonates. We measured heart rate and peripheral SpO 2 but found no significant changes during the stimulations. We did not measure blood pressure, which is unfortunate since systemic blood pressure changes may occur during the stimuli and can confound fNIRS signals. 17,52 The measurement of continuous blood pressure in neonates is difficult to realize and not routinely performed.
Another crucial factor explaining variability of HR patterns is the targeted brain region. For our measurements, we placed one sensor 1 to 2 cm above the inion, where the VC is suspected to be according to the international 10-20 system modified for neonates 53 and the second one on the left FTL (Fig. 1). Looking at the HR in the VC, we found in group A an expected increase and in group B a decrease in Δ½O 2 Hb, whereas in group C no significant changes in the individual responses were detected. A finding is that the left FTL also respond to visual stimulus in premature infants during spontaneous sleep. This has so far, to the best of our knowledge only been investigated once in a group of term-born infants. 46 The difference we found between these two brain regions was that in the VC the hemodynamic change was more robust with an increase (group A) or decrease (group B) in Δ½O 2 Hb, whereas in the FTL the answer was more biphasic with slight increase followed by a decrease in Δ½O 2 Hb in group A and in group B just the opposite direction (Fig. 3). A variation in the shape of the HR between different cortical regions is reported in the literature about brain development. 54 One explanation might be that brain maturation is not homogenous between cortical regions so that HR may vary from one brain area to another and based on the exact age at measurement. In summary of the recent review by Issard and Gervain, 54 the temporal cortex seems to present canonical responses earlier than in the occipital and frontal cortices.
The occipital cortex showed a canonical response at birth, but inverted later in infancy. The frontal cortex showed more variable responses, depending on the stimulus complexity and age of participants (>3 months). However, the prefrontal HR to visual stimulation in our study implies that the frontotemporal cortex may already contribute to some function in the premature brain. In addition to the different hemodynamic changes in the VC and FTL, we also calculated the functional connectivity (FC) between these two areas of the premature brain (see Sec. 4.7).

Can a Typical Stimulus-Evoked Hemodynamic Response Be Defined for Preterm-Born
Neonates?
The HR pattern in infants is still a controversial and critical issue in the neuroimaging literature. Information about this aspect in premature infants is particularly sparse. In contrast to the Δ½O 2 Hb increase typically observed in adults, young infants sometimes show an inverted (i.e., decreased) pattern of Δ½O 2 Hb or a negative fMRI blood-oxygenation-level-dependent (BOLD) response, corresponding to an increase in Δ½HHb. 55 However, results have been inconsistent and other studies reported typical Δ½O 2 Hb and Δ½HHb changes even in neonates. 46 In a recently published review investigating the pattern of the HR in fNIRS studies in term newborns less than 1 month old, the authors came to the conclusion that the typical response in newborns is not well established. 2 The majority of papers identified for this review reported an increase in Δ½O 2 Hb, with only some of them reporting a decrease in Δ½HHb. Furthermore, they did not observe a clear association between the reported response and stimulus type. In contrast, a further review including infants above the newborn period reported variable HRs with a canonical response of an increase in Δ½O 2 Hb and a decrease in Δ½HHb or an inverted response of a decrease in Δ½O 2 Hb and increase in Δ½HHb, with some studies also reporting changes in the same direction. 54 At present, however, direct comparison of optical signals between premature infants, neonates, infants at older ages, or adults is not simple because of the possibly different tissue locations measured with fNIRS (related to the age-dependent head anatomy), the choice of cortex areas studied, age, task demand, and individual responsiveness. 21 The comparison of fNIRS study results is also difficult since the measurement technique used and the specific signal processing applied to the data can lead to different results. 31,56,57,58 From our results in preterm-born neonates, we came to the conclusion that a typical HR, i.e., a canonical HRF, cannot be clearly defined so far. Thus, there is a desire to identify a typical response in the premature brain such that future studies can enable detection of normal and abnormal developmental patterns in the developing brain in this fragile patient group. Otherwise, it might be possible that the ambiguity in Δ½O 2 Hb response observed (i.e., either an increase, decrease, or no change in Δ½O 2 HbÞ reflects a normal pattern as these differences reveal the maturation of the developing premature brain on an individual basis. In this case, subgrouping is required as we performed in our study.

Why Does the Shape of the Hemodynamic
Response Depend on the Weight and Height of the Preterm-Born Neonate?
One surprising finding of our study is that the groups of neonates (A, B, and C) differ in their weight at measurement. This might be explained by the fact that neonates with higher weight are also more likely to be at an advanced state of brain development, despite no significant difference in head circumference. This may explain why group A showed HRs similar to those of adults. PMA per se might, therefore, not be the primary factor explaining the variable shapes of the HRs. This finding is in contrast to other studies. 22,23 In our results, weight at measurement was more important. Unfortunately, previous fNIRS studies of preterm-born neonates did not analyze a possible association between the weight and the HR shapes/signs. Attributing causality to weight and height requires more research.
The length and weight are most likely only parameters that correlate with the state of the neonate with respect to his or her (in utero, ex utero) development and physiological constitution. Future research is needed to investigate in detail the physiological reason for this correlation we discovered. Instead of treating height and weight as confounding factors for the analysis, these parameters give valuable information with regard to the subgrouping and thus possible relations to the neurodevelopmental state.

Hematocrit and the Total Hemoglobin
Concentration: Relevant to Explain the Variability of Hemodynamic Responses in Preterm-Born Neonates?
Correlations of the HRs with Hct and Hb were also observed in our study but only with respect to distinguishing neonates with or without statistically significant HRs at the VC. A dependence on the fNIRS-measured baseline [tHb] concentration as well as other baseline physiological parameters was not evident. This is contrary to the findings of Zimmermann et al. 34 In their study, a relationship between the shape of visually evoked HRs in preterm-born neonates and their [tHb] values was observed (higher , whereas in our case, we used the SCMD method that delivered [tHb] values that are influenced not only by the actual Hb concentration but also by the amount of light scattering in the tissue, which can vary between neonates and which could not be directly measured using the NIRS device. It is well known that the transition from fetal to adult Hb starts at birth, regardless of GA and that Hct decreases faster in more premature infants. 59 Previous assessments of cerebral blood flow and oxygen metabolism in premature infants using fNIRS also showed that only tissue StO 2 undergoes a decrease as a function of PNA until around 6 to 8 weeks, and it was speculated that this observation might be caused by a decrease in Hct. [60][61][62] In contrast, Watanabe et al. 23 suggested that a decrease in cerebral tissue StO 2 is not induced by the decrease in cHb or Hct but by the change in the oxygen affinity of hemoglobin.

Frontotemporal-Occipital Functional
Connectivity: What Is the Underlying Physiological Process?
The FC of the cortex in infants during sleep has been assessed based on average-time correlations among fNIRS signals from multiple locations of the cortex. 63 The developmental changes involved in the organization of cortical networks showed regional dependency and dynamic changes in the course of development. In the temporal, parietal, and occipital regions, FC increased between homologous regions in the two hemispheres and within hemispheres, and in the frontal regions, it decreased progressively. 63 However, it seems that FC among the cortical regions is generally higher during sleep than during wakefulness, especially for visual processing. 24 The visual system is one important functional network for humans because it processes visual information from the outside world and is consistently found to develop first. 64 In addition, the visual processing is extremely important for the maturation processes of other functional networks, such as attention networks, since attention requires accurate and prompt visual information processing. 11 The resting-state fMRI studies revealed that the maturing order of different brain functional networks seems to follow gradients, i.e., from posterior to anterior, from inferior to superior, from medial to lateral areas, as well as from primary functional systems to high-level cognition-related networks. 65 Of all the primary brain functional networks, the visual area seems to develop earlier and mature faster. 66 Intranetwork and internetwork FC reorganization could both increase during early brain development but at a different pace or slope. 67 For a better visualization of the fto-FC in our study, we plotted the stimulus-evoked time course of changes in Δ½O 2 Hb, Δ½tHb, and Δ½HHb at the left FTL with respect to those at the VC (Fig. 7). This illustrates nonlinear correlations between time series that have different oscillation frequencies.
These circle structures were clearly visible for Δ½O 2 Hb and Δ½tHb but not for Δ½HHb. There seems to already be some development of processing of visual stimulation in the premature brain. It will be interesting to investigate with future studies whether these findings are important development precursors for a later development of cognitive processes, such as social or other communicative tasks.  70 It cannot completely be excluded that stimulus-evoked systemic blood pressure increases also affected our measured HRs, although an advanced SCMD NIRS measurement technique has been used. Such stimulus-evoked systemic blood pressure increases might also be partially responsible for the ftc-FC observed in our study. Another factor is changes in local oxygenation and hemodynamics at the temporal side of the head elicited by mouth movements (causing muscle activity) that can influence the fNIRS measurements. 71 However, the muscle activity-related artifacts seem to be not significant for our study since the neonates were sleeping and were generally not showing significant mouth movements during the functional measurements.

Strengths and Limitations of the Study
Our study has a number of important strengths that give confidence in its findings and interpretations. (1) We used a NIRS device that enables reliable measurements of cerebral oxygenation and hemodynamics with the SCMD technique. 25 (2) A bright visual stimulation has been used in our study making sure that the neonates were exposed to a sufficiently intense stimulation to elicit HRs. between the three groups of neonates and specific demographic and baseline physiological parameters; we recognized that the weight of the neonate at measurement is a relevant factor explaining the variability of HR shapes.
Our study also has a few limitations. (1) Since our NIRS measurement device measured in the CW mode, the measurement of changes in scattering in the tissue was not possible. Since changes in scattering may happen during a functional task as proven by us, 72 a specific part of the variability of the measured HRs in our study might be related to this aspect.
(2) Our study analyzed data from 25 neonates. A larger number would have helped in the correlation analysis performed to investigate possible associations between the types of HRs and demographic as well as baseline physiological parameters.
(3) It cannot be ruled out completely that the measured HRs also contain influences from non-neuronal evoked changes in tissue oxygenation and hemodynamics, i.e., changes due to systemic physiology happening in the extracerebral and/or in the cerebral tissue. 1

Summary of the Main Findings
The main findings of our study were as follows. (1) HRs due to visual stimulation were observed in about one-third of the preterm-born neonates analyzed. (2) The magnitude and shape of the HRs varied from neonate to neonate (large intersubject variability). (3) The characteristics of the HR at the VC was used to subgroup the neonates into three groups [positive (A) or negative (B) HR, or absence (C) of a HR]. (4) The subgrouping of the neonates according to the HR characteristics was correlated with differences in the weight of the neonate at measurement. Differences between the groups were also evident for the height at measurement, Hct and cHb. (5) HRs at the left FTL were observed for groups A and B. (6) The HRs at the left FTL were correlated with those at the VC, indicating the presence of an fto-FC. (7) The hemodynamic changes measured were depthdependent and the SCMD-NIRS measurement technique proved to be useful and necessary to recover hemodynamic changes happening in the deeper layers of the head and thus being more associated with neurovascular coupling.

Implications of Our Study for Future Functional Near-Infrared Spectroscopy Studies in Preterm-Born Neonates
Our study shows that future fNIRS studies investigating stimulus-evoked responses in preterm-born neonates should perform a first-level analysis (single-subject analysis of HRs) as well as a second-level analysis (group average) based on subgrouping the neonates according to their HR characteristics. We are not aware of any other fNIRS studies in neonates that performed a subgroup analysis based on the shapes of the hemodynamic responses. We hope that our publication will motivate fNIRS researchers to perform such a subgrouping analysis in their future fNIRS studies, especially in full-term and preterm neonates. This aspect could be also relevant for fNIRS studies in adults. Furthermore, it is recommended to use a multidistance NIRS device approach that enables a depth-sensitive measurement. Ideally, absolute values of the absorption and scattering properties in tissue should be measured (using frequencydomain or time-domain NIRS). Covering the whole head with light sensors and detectors would be ideal since this would enable NIRS imaging (tomography) to be performed with the possibility to determine precisely the depth-resolved spatiotemporal changes of tissue oxygenation and hemodynamics. 74,75 Extending the fNIRS measured signals is another possibility for future fNIRS studies with preterm-born neonates. Especially measuring changes in the cytochrome-C-oxidase 76,77,78 with broadband NIRS (B-NIRS) may give additional information, in particular regarding the changes in cerebral metabolism elicited by the functional stimulation. In addition, a direct noninvasive measurement of cerebral blood flow with diffuse correlation spectroscopy (DCS) would give additional insights. B-NIRS and DCS can also be combined, as recently demonstrated. 79 We also recommend measuring systemic physiological changes in parallel, based on the systemic-physiology-augmented fNIRS neuroimaging approach recently introduced by us. 80,81 Especially the measurement of stimulus-evoked systemic blood pressure changes might be of particular interest. Also, the measurement of skin conductance changes in parallel to fNIRS neuroimaging measurement (as shown to be useful already by studies of our group 80,82 ) would add to additional insights into possible confounding physiological processes happening during the functional experiment.
In addition, an extended set of demographic and baseline physiological data should be assessed, including the weight at measurement and absolute hemoglobin (tissue and systemic) values. We also recommend extending the measured population to be more inclusive to neonates having a different race, ethnicity, and ancestry 83 than the normally measured population (Caucasian, bright hairs) in order to take into account that Western society is getting more diverse and multicultural. It should be also investigated whether the sex of the neonate has an impact on the HRs measured.
Future fNIRS studies should also use different stimulation modalities (visual, auditory, tactile, olfactory, etc.) as single stimuli or in combination in order to assess how the neonatal brain reacts with HRs. Varying the stimulus intensity and duration would add to the insights of these investigations. In addition, longitudinal fNIRS measurement in preterm-born neonates (i.e., measuring stimulus-evoked HRs at different PNAs in single neonates) would deliver valuable information about the development of neurovascular coupling and dynamic cerebral autoregulation.
Finally, machine learning and predictive modeling are exciting techniques that might be useful for classifying and quantifying the recorded HR of preterm-born neonates. For example, measured HRs in combination with evoked changes in systemic physiology, as well as baseline demographic data, could be processed by machine learning to predict medically relevant parameters such as a brain maturation score or a risk index for neurodevelopmental problems. OxyPrem AG had no influence on how the data were analyzed by F.S. specialization in computational intelligence from the University of Salerno (UNISA, Italy) in March 2018. He carried out his master's thesis on the reduction of false alarm rates in neonatal intensive care units using machine learning techniques, at the BORL, University Hospital Zurich. His research interests include biomedical signal processing, neuroinformatics, and machine learning.
Dirk Bassler is a professor of neonatology and director of the Department of Neonatology, USZ. He achieved a master's degree in health research methodology at the Department of Clinical Epidemiology and Biostatistics of the McMaster University in Ontario, Canada. He has extensive scientific experience in clinical research, among others as principal investigator in international randomized control trials. He is an expert in studying methodology and on compiling systematic reviews and meta-analyses.
Martin Wolf is a professor of biomedical optics at the University of Zürich, Zürich, Switzerland. He received his PhD from ETH Zurich. He heads the BORL, which specializes in developing techniques to measure and quantitatively image oxygenation of brain, muscle, tumor, and other tissues. His aim is to translate these techniques to clinical application for the benefit of adult patients and preterm infants.
Felix Scholkmann received his PhD from the University of Zürich, Zürich, Switzerland, in 2014. A postdoc and research associate at the BORL of the USZ, University of Zurich, and a research associate at the University of Bern, his research focuses on biomedical signal processing, biophotonics/neurophotonics (development and application of NIRS), neuroscience, integrative physiology, and biophysics.