1. Introduction

Understanding the functional and physiological changes in the human brain across a lifespan is an important aim in neuroscience. Brain measurements (i.e., biomarkers) that accurately describe the brain trajectory throughout the lifespan could promote holistic views of the brain, brain developments, and many brain diseases [1], especially in the early and late stages of life. Having such trajectory measures would be like obtaining a “growth chart” for the brain. By comparing with charts, we can assess the maturity of an infant brain, as well as the degenerative stage of an elderly brain. Moreover, we might be able to predict the progress of certain brain diseases based on the assessment of individual trajectories. It has been suggested that the infant brain is not a miniature version of the adult brain, but a continuously self-organizing system [2]. Short-range and long-range connectivity among multiple regions forms functional networks in this self-organizing system. For instance, the autism spectrum disorder (ASD) is a developmental brain disease in children. Establishing typical and atypical trajectories of brain development could improve our understanding and treatment of ASD [3]. The same is true for mild cognitive impairment and Alzheimer’s disease in older adults [4]. In short, a lifespan study has important implications for disease investigation because it provides a baseline for evaluating brain impairment in age-related neuropsychiatric disorders [1].

Identifying a proper imaging modality to track functional and physiological changes in the human brain across a lifespan is a major challenge. Many neuroimaging tools have been used to investigate brain development and degeneration, such as functional magnetic resonance imaging (fMRI), positron emission tomography (PET), and electroencephalography (EEG) [5–7]. However, fMRI and PET are costly and ill-suited for bedside or portable measurement, and EEG is sensitive to motion artifacts and has low spatial resolution. Also, due to the spatial confinement, fMRI and PET are not suitable for functional brain imaging studies on children (especially infants), the elderly, or people with claustrophobia [8].

In recent years, functional near-infrared spectroscopy (fNIRS), as a non-invasive hemodynamic-based neuroimaging technique, has been extensively employed to investigate brain function, brain development, and brain disease. Its advantages over fMRI and PET include good portability, higher temporal resolution, and lower sensitivity to motion artifacts [9,10]. Additionally, its spatial resolution is superior to that of EEG. fNIRS measures relative changes in concentrations of oxy-hemoglobin (Δ[HbO]) and deoxy-hemoglobin (Δ[Hb]) secondary to neuronal activity. Since the fNIRS has less restriction on subjects and higher motion-tolerance, it is an ideal functional brain imaging tool for studying infants, children, and the elderly, in daily settings [11,12].

Another challenge is identifying effective biomarkers that reflect the trajectory [13]. An important measure of the brain is complexity, which is broadly thought to be represented by the differentiation or diversity of the neural activity [14]. Several methods have been proposed to quantify brain complexity, mostly in the time-domain. These include permutation entropy [15] and multiscale entropy [4]. In the spatial domain, Scharthner et al. have [16] proposed amplitude coalition entropy and synchrony coalition entropy to quantify spatiotemporal variability in neural activity. Additionally, Puglia et al. [17] used the Lempel-Ziv algorithm to measure the spatial complexity of EEG during general anesthesia. They found that spatiotemporal Lempel-Ziv complexity increased with age at the baseline state in the 8 to 16 years age group. These intriguing findings demonstrate the importance and potential of spatial complexity in understanding the brain. However, all of these studies were based on EEG measurement. Given that fNIRS is still a novel modality in neuroscience, few groups have applied spatial complexity in fNIRS studies. However, fNIRS has much better spatial resolution (than EEG), and is sensitive to a different signal (i.e., blood signal). Thus, it would be useful to develop spatial complexity methods for fNIRS studies. To our knowledge, only one method derived from fNIRS signals, namely normalized global spatial complexity (NGSC), has been proposed and used to investigate ASD [18]. This led to the discovery that an ASD group’s NGSC value is significantly higher than that of a group with typical development. To further explicate the importance of spatial complexity, and to extend its application, in this study, we developed another promising spatial complexity method for fNIRS studies, and compare it to the NGSC method.

The time dependence between neural signals of spatially different brain regions is an important way to quantify their functional connectivity [19]. In fNIRS and fMRI studies, a time delay measure for low-frequency hemodynamic oscillations was employed to investigate the propagation of the cerebral blood flow (CBF) [20,21]. In the present study, we propose spatial time-delay entropy (STDE) to analyze the varying complexity of hemodynamic oscillations during brain development and degeneration. Traditional methods to assess the time delay have been based on Pearson’s correlation coefficient—a time-averaged metric which cannot characterize brain signals’ dynamic characteristics. In our study, we employ cross-power spectral density (CPSD) to calculate the cross-spectrum phase of two time series in a limited frequency band. The time delay map is achieved by calculating the phase lag at the peak coherence frequency. In order to investigate the variability of complexity throughout a lifespan, we analyze the sleeping or resting state fNIRS signals at the following age groups: early preterm, term, 3-4 months old, children (6-11 years old), adults and healthy seniors. We also calculated the NGSC for comparison. We hypothesized that the spatial complexity increased with brain development, and decreased with brain degeneration [4,22].

This paper is organized as follows. Section 2 describes the participants’ information and data recording, fNIRS preprocessing, spatial and temporal complexity measurements, and the statistical methods employed. In Section 3, we analyze the noise influence on STDE and the spatial complexity changes across age groups, and compare them with the NGSC and sample entropy (SampEn) measures. Finally, Section 4 is the discussion and conclusion.

2. Material and methods

2.1 Participants

We enrolled 162 healthy participants, ranging in age from infancy to healthy seniors. The fNIRS data were collected from two sites [23] and are described as follows.

1) The infant fNIRS dataset

We collected three age groups’ infant fNIRS data from Keio University, Tokyo, Japan: (1) 38 early preterm infants (21 boys and 17 girls, gestational age (GA) > 22 weeks); (2) 23 term infants (13 boys and 10 girls, GA $\ge $37 weeks); and (3) 35 3-4-month-olds infants (18 boys and 17 girls, GA >36 weeks). All subjects’ parent(s) signed an informed consent from the Keio University Hospital ethics committee before the measurements (No. 20090189). The experiments were conducted in a dimly lit room while the neonates and infants were sleeping (thus, not resting state). The sleeping state was judged by the frequency of their motor activity and rapid eye movement. We employed an Apgar score and Enjoji’s developmental test to evaluate the health states of both neonates (early preterm and term) and 3-4-month-old infants. The results of the tests showed that the condition and development of both neonates and infants were within the normal range. Detailed information on all subjects is presented in Table 1.

Table 1. Demographic characteristics of the subjects

View Table | View all tables in this article

We used an ETG-4000 NIRS system (Hitachi Medical Corporation) to measure the hemodynamic oscillations of Δ[HbO] and Δ[Hb] (millimolar millimeter). The layouts of the fNIRS probes used on the early preterm and term neonates are presented in Fig. S1(a)-(c). A set of a 22-channel array (3 sources and 5 detectors) covers the frontal area, and two 12-channel arrays (3 sources and 3 detectors) cover the left and right temporal areas. On the brain model, the X’s and the O’s are the positions of the light-emitting sources and the detectors. The exact channel arrays are shown in Fig. S1(d), where the red dots are the sources and the blue dots are the detectors. The channels are located between sources and detectors. The layouts of the fNIRS probes used on the 3-4-month-old infants are shown in Fig. S1(e)-(g). The two sets have 44 source-detector channels with 3 sources and 5 detectors arrays to cover the frontal, left and right temporal regions. The locations corresponding to these channels are shown in Fig. S1(h). In all infants’ fNIRS recordings, the distance between the source and the detector was approximately 2 cm.

The fNIRS system’s sampling rate is 10 Hz and the near-infrared lights’ wavelengths are 695 and 830 nm. We used the NIRS-SPM toolbox (resources were obtained from www.nitrc.org) to convert the detected light intensity information into Δ[Hb] and Δ[HbO] [24]. The differential path-length factor (DPF) was set as 6.1718 and 5.5374 for the two wavelengths (695 and 830 nm).

2) The children, adults, and the healthy seniors’ fNIRS dataset

We collected the fNIRS data from children to healthy seniors from Beijing Normal University, Beijing, China. There were three age groups: (1) children (n = 21, aged 6-11, 10 males, 11 females), (2) adults (n = 18, aged 19-27, 12 males, 6 females), and (3) healthy seniors (n = 27, aged 58-77, 15 males, 12 females). All subjects were right-handed based on the Edinburgh Handedness Inventory assessment. The subjects (or their parents for minors) signed an informed consent form prior to the experiment. The Institutional Review Board at Beijing Normal University, China approved all the experiments. Previous studies have reported detailed information on children [25] and adults [26]. All participants were healthy with normal weight. To exclude any participants with cognitive impairment, we conducted the standard neuropsychological tests—the Chinese version of the Mini-Mental State Examination (MMSE), the Beijing version of the Montreal Cognitive Assessment (MoCA), and for the healthy seniors we also conducted the Clinical Dementia Rating (CDR).

The detailed fNIRS experimental setup has been described by Cai et al. [27] . In our study, we used the CW6 fNIRS system (TechEn Inc., MA, U.S.A) to record the hemodynamic response in resting state. Then, we placed 12 light sources and 24 detectors constituting 46 measurement channels on a stretchable cap covering most of the head. The source-detector distance was 3.2 cm with wavelengths of 690 nm and 830 nm, respectively. The system’s sampling rate was 50 Hz. The layouts are shown in Fig. S1(i)-(k) and the corresponding locations of those channels are shown in Fig. S1(l). They cover the cortical frontal, parietal, temporal, and occipital regions. The experiment took around 11 minutes and was conducted in a “resting state”, defined as eyes closed, relaxed, and awake state [25].

We used the Homer2 toolbox to convert the light intensity into Δ[Hb] and Δ[HbO] (resources were obtained from https://www.nitrc.org/projects/homer2/) [9], and Gratzer’s extinction factors for the calculation, with DPF set to 6.0 for all participants [6].

2.2 fNIRS signal preprocessing

Signal quality is crucial to data analysis. As in our previous study [23], we identified any fNIRS signal that lacks a cardiac component (within a range of 1-2 Hz) in the spectrum as the noise signal, and removed them from further analysis. Then, we eliminated the motion artifact using the kurtosis-based Wavelet Filtering method [28]. After that, we downsampled the fNIRS signals to 2 Hz for further analysis. Finally, we extracted the fNIRS signals (i.e. Δ[HbO] and Δ[Hb]) in the 0.01-0.05 Hz, 0.05-0.1 Hz, and 0.01-0.1 Hz frequency bands with a zero-phase digital band-pass filter (MATLAB function filtfilt.m, butterworth order = 3) for spatial complexity analysis.

2.3 Spatial time delay entropy

Hemodynamic oscillation time dependence between regions reflects the hemodynamics changes in various stages [20,21]. In this study, we proposed a spatial entropy—spatial time delay entropy—to measure the overall hemodynamics variation with the brain development.

The STDE algorithm can be described as follows: First, consider two Δ[HbO] time series (see Fig. 1(a)) or two Δ[Hb] time series with 120 s epochs, and filter the data to a given frequency band (see Fig. 1(b)). In this study, we analyzed the fNIRS signals in the 0.01-0.1 Hz frequency band. Then, we used the CPSD method to calculate the cross-spectrum phase of two time series. The magnitude of the CPSD represents the coherence between two signals, as shown in Fig. 1(c). The phase of the cross spectrum can be achieved by the angle function (MATLAB: angle.m). The frequency point with the highest coherence value corresponds to the phase lag with two signals. Therefore, the time delay of two time series was calculated by $timedelay = \frac{{{\theta _{\max }}}}{{2\pi {f_{\max }}}}$, where ${f_{\max }}$ is the frequency point with the highest CPSD value of the two signals. The ${\theta _{\max }}$ is the corresponding phase under the frequency ${f_{\max }}$ . The time delay change corresponding to the frequency is presented in Fig. 1(d). Next, we calculated all the time delay values between all paired-two fNIRS channels. In this study, there were 46 source-detector channels, and in theory, we could have obtained 46*45/2 = 1,035 absolute values for each epoch, as shown in Fig. 1(e).

 

Fig. 1. The STDE calculation process based on two-channel Δ[HbO] signals. (a) The raw Δ[HbO] signals recorded from an early-preterm participant. (b) Filtered Δ[HbO] signals in the 0.01-0.1 Hz frequency band. (c) The CPSD of two Δ[HbO] signals. The red circle indicates the peak frequency. (d) Time delay values corresponding to the frequency changes. The time delay corresponding to the peak frequency is marked with two orthogonal red lines. (e) Time delay matrix of all paired-channels in an early-preterm participant. (f) Time delay histogram distributed in [0,20].

Download Full Size | PDF

Then, we calculated a histogram of all the time delay values with one subject, with the number of bins set at 20 (see Fig. 1(f)). The entropy of the time delay distribution can be calculated by $E ={-} \sum\limits_{j = 1}^b {{p_j}\log {p_j}}$, where ${p_j}$ is the proportion of bins and b is the number of the bins. We termed the normalized entropy STDE, and calculated it by ${E_{STDE}} = E/\log \frac{1}{N}$, where N is the number of bins.

To calculate the time delay and entropy, it is important to determine a maximum time delay threshold for all age groups. It is well-known that a correlation coefficient greater than 0.3 is meaningful for investigating the functional connection in fNIRS and fMRI studies. Therefore, we selected paired channels with correlation coefficients exceeding 0.3 for calculation [29]. Approximate normal distributions for all time delay values with cross coefficients greater than 0.3 for early preterm, term, 3-4 months, children, adults and healthy seniors are presented in Fig. 2(a)-(f). The fig’ x-axes are the time delay values after power conversion by the boxcox function (Matlab: boxcox.m). Based on these distributions, we used a significance threshold of p < 0.05 for each measure to avoid spurious results. The significant thresholds for early preterm, term, 3-4 months, children, adults and healthy seniors were 2.27, 2.24, 2.08, 2.18, 2.01 and 2.00, respectively. The corresponding time delays for these age groups were 11.36 s, 13.91 s, 14.65 s, 18.10 s, 12.07 s and 9.38 s. For the convenience of calculation, we analyzed only time delay values less than 20 s. All time delay values with cross coefficients greater than 0.3 for early preterm, term, 3-4 months, children, adults and healthy seniors are presented in Fig. S2(a)-(f).

 

Fig. 2. Approximate normal distribution plot after power conversion of all time delay values for early preterm, term, 3-4 months, children, adults and healthy seniors.

Download Full Size | PDF

2.4 Normalized global spatial complexity

Jia et al. proposed the NGSC method [18] in which they employed the Shannon entropy of a feature spectrum of its covariance matrix to quantify spatial complexity by using fNIRS signals. The NGSC algorithm is described below:

Firstly, we conducted the principal component analysis (PCA) on the Δ[HbO] or Δ[Hb] signals, which yielded m principal components and a spectrum of eigenvalues, where m is the number of channels across the scalp. Then, to assess the relative contribution of each principal component’s relative contribution to the total variance, we normalized the eigenvalues of principal components. The normalized eigenvalue of the i-th principal component was calculated as

The lowest NGSC value was 0, meaning that the hemodynamic signals of all fNIRS channels consisted of exactly one principal component or spatial mode. This indicates a minimum brain spatial complexity. Meanwhile, the highest NGSC value was 1, indicating that the total data variance was uniformly distributed across all m principal components, i.e., maximum brain spatial complexity.

2.5 Sample entropy

Richman et al. have proposed SampEn [30], an important measure for describing the unpredictability and randomness of a non-linear time series. The detailed calculation is described below:

Firstly, we consider a signal with N data points from a time series $\{ x(n)\} = x(1),x(2),\ldots ,x(N)$ to define SampEn. We set the tolerance coefficient (r) and the embedding dimension (m) to construct the template vector ${X_m}$ and define the generalized formula of ${X_m}$ as: ${X_m}(k) = \{{x(k),x(k + 1),\ldots ,x(k + 1 - m)} \}$.

Then, we compute how many distances were less than $R = (R = r \times SD)$ (standard deviation of the original data)). The distance function is defined as:

The total number of templates matches under the embedding dimension m is:

Finally, we define SampEn as:

In this study, we set the parameters as m = 2 and r = 0.2SD, based on a previous study [31]. SampEn measures the complexity of one time series. For each channel, we can obtain a SampEn value. Thus, we used the mean SampEn value for all channels to represent the time domain complexity, for comparison.

2.6 Time delay calculations based on simulated signals

Time delay can be calculated with the cross-correlation method. However, in this study, we used CPSD to achieve the time delay value. To evaluate performance in the time delay calculation of CPSD methods, we used two simulated signals by combining two sinusoid waves, at 0.2 Hz and 0.5 Hz, under a sampling rate of 10 Hz (see Fig. 3(a)). Compared to signal ${y_2}$, the signal ${y_1}$ phase lagged $\frac{\pi }{2}$ at 0.2 Hz and phase lead $\frac{\pi }{3}$ at 0.5 Hz. The data length of the two simulated signals was 300 s. The formulas are:

 

Fig. 3. Comparison of time delay calculation methods. (a) simulated time-series of ${y_1}$ and ${y_2}$. (b) CPSD of two simulated time-series. (c) The phase of two simulated time-series corresponds with the frequency changes. The red circles and dashed lines+ indicate the peak frequency and corresponding phase. (d) Time delay values corresponding with the frequency changes. The time delay corresponding to the peak frequency is marked with two orthogonal red lines and a red circle.

Download Full Size | PDF

The signals of ${y_1}$ and ${y_2}$ have CPSD peaks at 0.2 Hz and 0.5 Hz, respectively (see Fig. 3(b)). The CPSD peak at 0.2 Hz is greater than the peak at 0.5 Hz.

There are two steps to calculating the cross-spectrum phase. In the first step, we used the MATLAB function cpsd.m to derive a vector of frequencies and corresponding cross power spectral density which is a complex value. In the second step, we used the MATLAB function angle.m to derive the cross-spectrum phase. We can achieve the exact phase lag in a specific frequency based on this procedure. The ${\Phi _{0.2}}$ and ${\Phi _{0.5}}$ are the cross-spectrum phase at 0.2Hz and 0.5Hz, respectively. Based on above calculation, the cross spectrum phases of ${y_1}$ and ${y_2}$ are $\frac{\pi }{2}$ at 0.2 Hz, and $- \frac{\pi }{3}$ at 0.5 Hz, respectively. This means that compared to time-series ${y_2}$, the time-series ${y_1}$ phase lagged $\frac{\pi }{2}$ at 0.2 Hz, and phase advanced $\frac{\pi }{3}$ at 0.5 Hz (see Fig. 3(c)).

The time delay can be calculated by $timedelay = \frac{{{\Phi _f}}}{{2\pi f}}$, where f is the frequency point and ${\Phi _f}$ is the corresponding cross-spectrum phase at this frequency point. Therefore, the time delay values of ${y_1}$ and ${y_2}$ are 1.25 s at 0.2 Hz and -0.33 s at 0.5 Hz (see Fig. 3(d)).

We regarded the time delay determined by the dominant frequency as the time delay between the two fNIRS signals. Therefore, we used a time delay corresponding to the maximum CPSD as the delay value between the two signals.

Further, the CPSD (calculating the time delay) was performed in a noisy environment. As shown in Fig. S3, we added random noise of varying intensities (from 1 dB to 25 dB) to the simulated signals. The results showed that the time delay and phase of the corresponding frequencies had not changed significantly. In our simulation, the only deviation of the time delay was observed at 0.2 Hz, and only by 1.6% (from 1.25 s to 1.27 s) (see Fig. S3). Therefore, the CPSD is resistant to the noise introduced in calculating the time delay.

2.7 Impacts of random noise on STDE and NGSC calculations

In order to investigate the noise impact on the STDE and NGSC calculation measures, we added random noise at varying intensities into the raw optical density signals of varying numbers of channels. Then, we used the NIRScout 816 imaging system (NIRx Medical Technologies LLC) for data collection. Based on our previous study [23], we knew that the raw optical density signals had distinct heart rate information, and therefore we assumed them to be “clean” signals. For convenience, we selected 21 channels’ fNIRS signals from normal men which had at least 10 minutes of resting-state data for simulation. We selected the noise signal by using only those fNIRS signals which had been recorded in a dark environment. Next, we placed the detectors in a dark environment away from the subject so that we could get a pure noise signal containing no hemoglobin signals. Then, we added the noise signal into the raw optical density signals with the noise coefficient from 0 to 1. To obtain the changes in STDE and NGSC at varying noise levels, we set the step to 0.001 within a range of 0 to 1. Then, we calculated the signal-to-noise ratio (SNR) for various levels of noise by dividing the raw optical density energy by the noise energy. Because STDE and NGSC are spatial complexity measurements, the number of contaminated channels is also an important influence. We set the step to 1 within a range of 0 to 21 (signal channels), and added noise signals of various levels (the SNR step to 1 within the range 0 to 40 dB) for simulation. Finally, we calculated STDE and NGSC under various SNR levels and numbers of contaminated channels.

2.8 Statistical analysis

The aim of this study was to (1) evaluate STDE’s performance in tracking brain development and degeneration, and (2) compare it with measures of NGSC. We used a box plot to show the sample distributions in participant-based statistical analysis. For each of these indices, we used a Lilliefors test to determine whether the values had normal distributions. Then, we used a Kruskal-Wallis test and a multiple comparison test to estimate the statistical significance of the differences between the indices of different age groups. We performed a post-hoc test with Bonferroni correction to adjust the p-values. Finally, any adjusted p-values less than 0.05 were considered as indicating a statistically significant difference. In all figures and tables, we used ‘*’, ‘**’, and ‘***’ to represent p < 0.05, p < 0.01, and p < 0.001, respectively.

3. Results

3.1 Effects of STDE and NGSC with the influence of noise

In this study, we analyzed noise’s influence on spatial complexity measures. Because the calculation incorporated all fNIRS channels which cover the brain, we varied the number of contaminated channels, as well as the noise intensity in the simulation. Figure 4 shows the effects of STDE and NGSC under the influence of noise at varying intensities and with different numbers of contaminated channels. Figure 4(a) presents the raw optical density signals (785 nm and 830 nm) derived from the fNIRS system, and Fig. 4(e) and (i) are the optical density signals contaminated by random noise with different intensities. The relevant power spectra of raw optical density and optical density signal with noise are shown in Fig. 4 (b), (f) and (j), respectively. We converted the optical data into hemodynamic oscillations of Δ[Hb] and Δ[HbO] with the Homer2 toolbox which uses Gratzer’s extinction factor [32]. The corresponding Δ[HbO] and Δ[Hb] signals are shown in Fig. 4(c), (g) and (k). The filtered Δ[HbO] and Δ[Hb] signals in the 0.01-0.1 Hz frequency band are shown in Fig. 4(d), (h) and (l). As shown in figure 4, the Δ[HbO] and Δ[Hb] waveforms were distorted with the influence of noise.

 

Fig. 4. The effects of the STDE and NGSC for Δ[Hb] and Δ[HbO] with the influence of noise. (a, e, i) Light intensity signals for raw, SNR = 25 dB, and SNR = 16 dB. (b, f, j) Spectrum power from the 785 nm light intensity from the left Fig, respectively. (c, g, k) Generated Δ[Hb] and Δ[HbO] signals calculated by the modified Beer-Lambert law from the raw signals in (a, e, i), respectively. (d, h and l) Filtered Δ[Hb] and Δ[HbO] in the 0.01-0.1 Hz frequency band derived from the left Fig. (m) STDE_HbO changes with SNR increase and contaminated noise channel increase. (n) STDE_Hb changes with SNR increase and contaminated noise channel increase. (o) NGSC_HbO changes with SNR increase and contaminated noise channel increase. (p) NGSC_Hb changes with SNR increase and contaminated noise channel increase.

Download Full Size | PDF

Next, we evaluated noise’s influence on the indices (i.e., STDE and NGSC) by varying amounts of contaminated fNIRS channels (from 0 to 21, step 1) and varying noise levels (SNR from 0 to 40 dB, step 1 dB). If the indices change more than 10% as a result of increasing noise levels and increasing number of contaminated channels, we consider the indices are “unreliable”. Because STDE and NGSC can be calculated based on the hemodynamic oscillations of either Δ[HbO] or Δ[Hb], we analyzed the STDE and NGSC changes in both signals. We termed the STDE based on Δ[HbO] and Δ[Hb] as STDE_HbO and STDE_Hb, respectively; and the NGSC based on Δ[HbO] and Δ[Hb] as NGSC_HbO and NGSC_Hb, respectively.

Figure 4(m) and (n) present the STDE_HbO and STDE_Hb with varying noise levels and number of contaminated fNIRS channels, respectively. As shown in Fig. 4(m), the STDE_HbO is relatively stable when SNR is greater than 25dB (a higher SNR value indicates superior signal quality) and when there are less than 10 contaminated channels. The STDE_HbO abruptly increased (more than 10% relative to ‘baseline’) when there were more than 10 contaminated channels and SNR was less than 25 dB. The STDE_HbO became unreliable when the SNR was less than 25 dB and noise-contaminated channels exceeded 47.6% of the total channel. Compared with STDE_HbO, the STDE_Hb decreased when there were more than 5 contaminated channels and SNR was less than 15 dB, indicating that the STDE_Hb was more prone to noise.

Figures 4(o) and 4(p) present changes in NGSC based on Δ[HbO] or Δ[Hb] with varying noise levels and numbers of contaminated channels. As shown, the changes in NGSC_HbO indices are similar to those of NGSC_Hb, which dropped when there were more than 10 contaminated channels and SNR was less than 25 dB. However, compared with the STDE_HbO, the NGSC_Hbo and NGSC_Hb values decreased with the increase in contaminated channels and the decrease in SNR. We used the offset of the metrics to compare the performances of the STDE and NGSC under noisy conditions (i.e., anti-noise performance). Under the same noise levels, the larger index change (i.e., offset) means that the index is more prone to noise (i.e., poor “anti-noise performance”). As a result, the biggest offsets of STDE_HbO, STDE_Hb, NGSC_Hbo and NGSC_Hb were 44.31%, 56.87%, 90.08%, and 96.14%, respectively (with SNR = 0∼40 dB; contaminated channel number = 0∼21). This indicates that among these four spatial complexity measures, STDE_HbO performs the best at noise resistance. Also, the simulations showed that STDE_HbO, STDE_Hb, NGSC_Hb and NGSC_Hb are reliable when the SNR exceeds 25 dB and the noise-contaminated channels are less than 47.6% of the total channel.

3.2 Spatial complexity changes paralleling brain development and degeneration

In this study, we proposed a new measurement to quantify brain development and degeneration by quantifying the spatial complexity of hemodynamic oscillations. In Fig. 5, we calculated the time delay matrixes of six participants from six age groups—early preterm, term, 3-4 months old, children, adults, and healthy seniors. We also analyzed the Δ[HbO] in the 0.01-0.1 Hz frequency band. In these time delay matrixes, the red dots indicate longer time delay between two hemodynamics oscillations, while the blue dots indicate a short time delay. The histograms are the distribution of the time delay values in the upper panel. As shown, the time delay distribution widens with age within the first 3-4 months of life. This indicates that the time delay between all channels diversifies with age in the infant group. Compared with adults, the distribution of the time delay values narrows in healthy seniors. This indicates that the variation in the time-delay dependency between all channels decreases in healthy seniors. In addition to the Δ[HbO], we also analyzed the STDE based on Δ[Hb]. Fig. S4 shows the time delay matrixes and time delay value distributions for Δ[Hb] in the 0.01-0.1 Hz frequency band in different age groups. As shown, the distribution for the infant of 3-4 months old widens more than those of early-preterm and term. However, we observe no difference between the healthy seniors and the adults.

 

Fig. 5. Time delay map and time delay histogram for Δ[HbO] for various subjects across age groups in the 0.01-0.1 Hz frequency band. (a)-(f) are the time delay map and time delay histogram for the subjects with age of PNA = 17, 5 and 132 days. The color of the dots in the time delay map represents the time delay value; the red represent a time delay approach of 20, and the blue represent a time delay approach of 0. The histograms show the time delay distributed in [0,20]. (g)-(l) are the time delay map and time delay histogram for subjects aged 6, 27 and 70 years.

Download Full Size | PDF

To further investigate the complexity change with brain development and degeneration at the group level, we statistically analyzed the STDE, NGSC and SampEn measurements in each age group. We calculated the significances of the spatial entropy indices among various age groups with the Kruskal-Wallis test and multiple comparison tests. Figure 6(a) and (b) show boxplots for STDE_HbO, STDE_Hb and SampEn in each age group in the 0.01-0.1 Hz frequency band. As shown, the STDE_HbO in the 3-4 months old age group exceeded those of the early-preterm and term age groups. Among all age groups, children have the highest values. Compared with the children and adults, the healthy seniors’ values decreased. The STDE_Hb in the 3–4-month-old age group exceeded that of the early-preterm and term age groups (see Fig. 6(b) and Table 2). However, there were no significant differences between the age groups of children, adults, and healthy seniors. The box plots for STDE_HbO and STDE_Hb in the 0.01-0.05 Hz and 0.05-0.1 Hz frequency bands are shown in Fig. S5(a)-(b) and S6(a)-(b), respectively. As shown, the STDE_HbO and STDE_Hb can distinguish the age groups of early-preterm, term, and 3–4-month-old infants in 0.01-0.05 Hz. Of note, the STDE_HbO has superior performance in distinguishing the children, adults, and healthy senior age groups. The mean ± SD of STDE_HbO and STDE_Hb in all three frequency bands are presented in Table 2, S1 and S2. All statistically significant difference values are presented in Table S3-S8. All these findings were statistically significant at a level of p < 0.05, or lower.

 

Fig. 6. Boxplots of STDE and NGSC for 0.01-0.1Hz. (a) is a boxplot of STDE_HbO for early-preterm, term, 3-4 month, children, adults and healthy seniors within 0.01-0.1Hz. (b) is the boxplot for STDE_Hb for early-preterm, term, 3-4 month, children, adults and healthy seniors within 0.01-0.1Hz. (c) is the boxplot for NGSC_HbO for early-preterm, term, 3-4 month, children, adults and healthy seniors within 0.01-0.1Hz. (d) is the boxplot for NGSC_Hb for early-preterm, term, 3-4 month, children, adults and healthy seniors within 0.01-0.1Hz. (e) is the boxplot for SampEn_HbO for early-preterm, term, 3-4 month, children, adults and healthy seniors within 0.01-0.1Hz. (f) is the boxplot for SampEn_Hb for early-preterm, term, 3-4 month, children, adults and healthy seniors within 0.01-0.1Hz.

Download Full Size | PDF

Table 2. STDE, NGSC and SampleEntropy range at 0.01-0.1Hz band (mean${\pm} $SD)

View Table | View all tables in this article

For the NGSC measure, the NGSC_HbO increased with age until adulthood and then decreased, as shown in Fig. 6(c) and 6(d). However, there were no significant differences between the term and 3-4 months infant age groups (p > 0.05). The NGSC_Hb outperformed the NGSC_HbO at distinguishing the infant age groups (early preterm, term, and 3-4 months old). The NGSC_HbO and NGSC_Hb indices increased from children to adults, and then decreased in the healthy seniors. These trends are different from the STDE, which has the highest value in children. The statistical analysis of the NGSC indices in the 0.01-0.05 Hz and 0.05-0.1 Hz frequency bands are presented in Fig. S5(c)-(d) and S6(c)-(d), respectively. NGSC_Hb has superior performance in distinguishing the developmental states in infants (early-preterm, term, and 3-4 months old), as well as in distinguishing the degeneration from adults to healthy seniors. The statistically significant difference values are presented in Table S9-S14. The mean ± SD of NGSC_HbO and NGSC_Hb in all three frequency bands are presented in Table 2, S1 and S2.

For the SampEn measure, the SampEn_HbO and SampEn_Hb values in the early preterm age group were higher (but not statistically significant) than those of the term and 3–4-month-old age groups, as shown in Figs. 6(e) and 6(f). Also, there were no significant differences between the early pre-term, term, and 3-4 month infant age groups (p > 0.05) for the SampEn measure. This means that SampEn could not distinguish the developmental changes in infants. SampEn_HbO in 0.01-0.1 Hz had the highest value in adults, compared to other groups (p < 0.01). SampEn_Hb in 0.01-0.1 Hz also had the highest value in adults, higher than other groups (except children) (p < 0.001). The statistical analysis of the SampEn indices in the 0.01-0.05 Hz and 0.05-0.1 Hz frequency bands are presented in Figs. S5(e)-(f) and S6(e)-(f), respectively. The mean ± SD of SampEn_HbO and SampEn_Hb in all three frequency bands are presented in Tables 2, S1 and S2. The statistically significant differences are presented in Tables S15-S20.

4. Discussion

Tracking brain development and degeneration with brain functional imaging tools has been a critical aim of neuroscience in recent years. Some fNIRS studies have explored functional changes in the brain, however, over a short period [33,34]. In this study, we have assessed brain development over a long lifespan (neonates to healthy elderly people), and more importantly, we validated that STDE, a new spatial complexity measurement, as an effective fNIRS biomarker, can be applied to different age groups (without too much variation), and used to track brain development and degeneration.

Our main findings are as follows. First, complexity is an important, but understudied parameter for revealing the mystery of the brain [35]. The method we proposed (STDE) could be used to investigate the spatial complexity of the hemodynamic oscillations, which has rarely been studied (especially in fNIRS). Our study has not only explored the spatial complexity in fNIRS, but has also produced the first complexity-based trajectory for healthy subjects over a long lifespan. Second, we have verified STDE’s performance under various noise conditions, and have validated its effectiveness in distinguishing brain states in various age groups. The STDE_HbO had the best noise resistance among four of the spatial complexity indices. Third, we have compared the STDE with another spatial complexity measure (NGSC) and a time-domain entropy (SampEn). We evaluated them based on their unique strengths in reflecting various aspects of brain growth. The findings in this paper will pave the way for subsequent use of the methods/parameters as a biomarker for brain development assessment.

4.1 STDE can track cerebral vascular changes during brain development and degeneration

Brain development and degeneration throughout the lifespan is a complex process entailing physiological, neural, and functional network changes. The physiological changes include development and degeneration in blood volume/flow/speed and oxygenation [36]. The neural alternations include the changes in synaptogenesis in the cortex, the dendritic connections, and the white-matter connectivity [37]. Meanwhile, development and degeneration of the brain functional networks includes the changes in cortex connections [38], brain networks [2], and functional complexities [39]. None of these structural or functional changes are linear with aging.

From the perspective of physiological development and degeneration in the brain, a previous study [40] has suggested that brain development causes the following vascular changes in the infant state: (1) increased cerebral capillary and venous blood flow; (2) increased cerebral blood volume for the increased capillary density; and (3) increased mean arterial pressure limited by cerebral auto-regulation. On the contrary, the physiological changes in the aging brain include declining metabolic rates and cerebrovascular function, as well as neuron degeneration, and these processes are interrelated [41]. It is also known that the aging brain may induce hypoperfusion and blood-brain barrier dysfunction. Cerebrovascular aging entails (1) microvascular rarefaction, (2) oxidative stress in inflammation, (3) arterial stiffness, (4) narrowing of the vascular lumen, and (5) endothelial replicative senescence [41]. In short, the cerebrovascular system is intertwined with the process of brain development, maturation, and aging.

Based on our previous research on fMRI and fNIRS, we found that time delays between the low-frequency oscillations of fMRI/fNIRS signals are related to the propagation of the cerebral blood [20,21]. The range of the time delays from the whole head represents the cerebral blood transit time [20]. The STDE measures the diversity of hemodynamic spatial time-delay in the cerebral cortex. Therefore, we used the STDE to quantify cerebrovascular-related development and degeneration, including changes in cerebral blood flow, volume, and oxygenation. In our study, the STDE values increased from preterm and term to 3-4 months old, and had the highest values in children (6-11 years old) both for Δ[HbO] and Δ[Hb] in the fNIRS frequency band of 0.01-0.1 Hz. Then, they decreased in adults and healthy elderly people (see Figs. 6(a)-(b)).

Similar trajectories have been found in CBF from other imaging studies. For example, in a PET [42] and a transcranial Doppler study (TSD) [43] on CBF from participants ranging from 10 days to 16 years old and 3 to 18 years old, respectively, CBF was found to increase at an early age. It then peaked around age 7 before starting to decrease, and stabilize at age 15 [43]. In another study [44], the researchers employed intracranial 4D flow MRI to analyze CBF in children (7 months - 17 years old) and adults (19-61 years old). They measured the peak velocities in several major arteries of the brain and observed a similar phenomenon of age-related CBF fluctuations before age 20. These observations demonstrate the close relationship between STDE and hemodynamic parameters, such as CBF.

4.2 Relationship between STDE and NGSC

In this study, we compared the performance of NGSC’s spatial complexity in tracking the brain development and degeneration with STDE. Although they are all spatial complexity measures, the theoretical basis of the NGSC approach differs from that of STDE. The NGSC calculations were based on PCA and normalized entropy, which is thought to estimate the global functional connectivity between brain signals. While STDE measures the time-dependence relationship between two hemodynamic oscillations and reflects the cerebrovascular change between two regions [20,21].

Jia et al. [18] used NGSC to analyze differences in global functional connectivity levels between children with ASD and age- and gender-matched children with normal development. It has been suggested that the NGSC of neural signals changes inversely with global functional connectivity (FC). They found that the NGSC was greater in ASD children than that in healthy controls, which implied that the ASD children have weaker FC. Cai et al. investigated brain functional development from early childhood (7-8 years old) to adulthood, based on resting-state fNIRS signals [27]. They used graph theory to construct brain functional network analysis and a metric of global efficiency. Their findings showed that young children had higher FC (i.e., low spatial complexity) than adults. In the present study, we found that NGSC intensified from childhood to adulthood, which matches cognitive development and is consistent with the existing studies. Therefore, NGSC is good at measuring changes in brain FC at the global level.

Compared with NGSC, we found that the trends in STDE indices from infancy to childhood and adulthood are more similar to the development of cerebrovascular parameters (i.e., CBF) than to that of FC.

For all calculations, we analyzed Δ[Hb] and Δ[HbO] in the 0.01-0.1 Hz frequency band. The results are consistent with existing studies [30]. The origins of these blood-related low frequency oscillations are still unclear. The sources of these oscillations include neural (resting state neural activity) and vascular (fluctuations in blood oxygenation and flow arising from cardiac and blood pressure changes, respiratory effects, etc.) components [45]. The results of our study suggest that STDE and NGSC characterize the physiological and neural components, respectively, thus indicating that there are both neuronal and physiological components in the low frequency signals of Δ[Hb] and Δ[HbO].

In summary, this study supplies a tool for assessing development and aging, from infancy to old age. STDE and NGSC quantify brain development from the perspective of spatial complexity, which is more suitable for analyzing non-linear and non-stational brain oscillations. The spatial complexity can characterize the brain development better than the time-domain complexity measure (i.e., SampEn). Furthermore, these two approaches have potential value in evaluating cerebrovascular deficits related to brain injuries (such as stroke and cerebral ischemia), as well as developmental and degenerative brain diseases (such as ASD and Alzheimer’s disease), in a ‘resting-state’.

4.3 Other potential fNIRS biomarkers

In addition to the complexity measures, many theories have been employed to investigate the functional changes in spatial-temporal scales [36] such as FC [46,47] and graph theory [22]. FC analysis is an important tool for investigating the coupling and cooperation between brain regions during tasks or resting states [46,47]. Cortical regions are thought to be functionally connected when these regions activate simultaneously [48]. Several studies have also used the FC to investigate cognitive decline [49], brain development [50] and degeneration [51]. The human brain is a complex network with non-trivial global and local topological features [52]. Using graph-theory approaches, fNIRS-derived brain networks can quantify the functional and topological characteristics of network organizations within the brain [53]. Graph-theory based metrics have also been used to investigate brain growth trajectory [54,55]. These studies suggest that graph-theory approaches are powerful tools for characterizing brain status. Since the FC and graph theory are sensitive to the NIRS probe configurations and coverage, it was hard to calculate them in this study. We hope to compare all of these parameters in future lifespan studies.

5. Limitations

This study has the following limitations. First, the fNIRS probe covers a smaller cortical region in infants (preterm, term, 3-4 months) than in the other age groups. From the perspective of brain development research, it is always valuable to investigate the whole brain [56]. However, due to the limitations of the fNIRS device and the study population, whole-brain coverage might not be possible. Based on our simulation and results, to obtain robust results, we suggest 1) keeping at least 20 fNIRS channels, and 2) covering at least 50% of the cortical region. Second, we collected the infants’ fNIRS data during sleep as in similar infant studies [57,58], while collecting fNIRS data during resting states for other age groups. Therefore, the conclusions of this study could only be drawn based on these two conditions. Third, the age gaps between children, adults, and healthy elderly people were too wide; there were no participants in the age range from 1 to 6 years old. However, other studies have suggested that profound brain growth transpires during this period [59,60]. Therefore, meta-analysis should be conducted with more datasets and smaller age gaps in the future. Finally, all the fNIRS channels in this study are single-distance channels. It has been documented that fNIRS is prone to “noise” from the superficial layer. People use multi-distance optodes and independent component analysis to remove the superficial noise [61]. The noise is more profound in adults than in infants (due to less superficial noise). In many fNIRS studies (including ours), the number of fNIRS optodes is limited, thus restraining the multi-distance measurements. Based on our results, a source-detector distance of 3 cm is sufficient for measuring a reliable trajectory. In the future, it would be useful to conduct and compare multi-distance and signal-distance STDE measurements.

6. Conclusion

In this study, we have presented a novel spatial complexity to assess brain development and degeneration. In summary, the spatial complexity of brain near-infrared signals increases as the brain develops and diminishes in old age. The STDE is sensitive to the vascular and metabolic changes associated with brain development in different age groups. Meanwhile, the NGSC is sensitive to the functional changes during brain development. This study has demonstrated that STDE and NGSC derived from fNIRS are potential biomarkers for assessing brain development across the lifespan.