Skip to main content
Advertisement
  • Loading metrics

Oscillator decomposition of infant fNIRS data

  • Takeru Matsuda ,

    Roles Conceptualization, Formal analysis, Funding acquisition, Investigation, Software, Writing – original draft, Writing – review & editing

    takeru.matsuda@riken.jp

    Affiliation RIKEN Center for Brain Science, RIKEN, Wako, Japan

  • Fumitaka Homae,

    Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Writing – original draft, Writing – review & editing

    Affiliations Department of Language Sciences, Tokyo Metropolitan University, Tokyo, Japan, Research Center for Language, Brain and Genetics, Tokyo Metropolitan University, Tokyo, Japan

  • Hama Watanabe,

    Roles Data curation

    Affiliation Graduate School of Education, The University of Tokyo, Tokyo, Japan

  • Gentaro Taga,

    Roles Conceptualization, Data curation, Funding acquisition, Writing – original draft, Writing – review & editing

    Affiliation Graduate School of Education, The University of Tokyo, Tokyo, Japan

  • Fumiyasu Komaki

    Roles Conceptualization, Funding acquisition, Supervision, Writing – review & editing

    Affiliations RIKEN Center for Brain Science, RIKEN, Wako, Japan, Graduate School of Information Science and Technology, The University of Tokyo, Tokyo, Japan, International Research Center for Neurointelligence (IRCN), The University of Tokyo, Tokyo, Japan

Abstract

The functional near-infrared spectroscopy (fNIRS) can detect hemodynamic responses in the brain and the data consist of bivariate time series of oxygenated hemoglobin (oxy-Hb) and deoxygenated hemoglobin (deoxy-Hb) on each channel. In this study, we investigate oscillatory changes in infant fNIRS signals by using the oscillator decompisition method (OSC-DECOMP), which is a statistical method for extracting oscillators from time series data based on Gaussian linear state space models. OSC-DECOMP provides a natural decomposition of fNIRS data into oscillation components in a data-driven manner and does not require the arbitrary selection of band-pass filters. We analyzed 18-ch fNIRS data (3 minutes) acquired from 21 sleeping 3-month-old infants. Five to seven oscillators were extracted on most channels, and their frequency distribution had three peaks in the vicinity of 0.01-0.1 Hz, 1.6-2.4 Hz and 3.6-4.4 Hz. The first peak was considered to reflect hemodynamic changes in response to the brain activity, and the phase difference between oxy-Hb and deoxy-Hb for the associated oscillators was at approximately 230 degrees. The second peak was attributed to cardiac pulse waves and mirroring noise. Although these oscillators have close frequencies, OSC-DECOMP can separate them through estimating their different projection patterns on oxy-Hb and deoxy-Hb. The third peak was regarded as the harmonic of the second peak. By comparing the Akaike Information Criterion (AIC) of two state space models, we determined that the time series of oxy-Hb and deoxy-Hb on each channel originate from common oscillatory activity. We also utilized the result of OSC-DECOMP to investigate the frequency-specific functional connectivity. Whereas the brain oscillator exhibited functional connectivity, the pulse waves and mirroring noise oscillators showed spatially homogeneous and independent changes. OSC-DECOMP is a promising tool for data-driven extraction of oscillation components from biological time series data.

Author summary

The functional near-infrared spectroscopy (fNIRS) can detect hemodynamic responses in the brain and the data consist of bivariate time series of oxygenated hemoglobin (oxy-Hb) and deoxygenated hemoglobin (deoxy-Hb) on each channel. In this study, we investigate oscillatory changes in fNIRS signals of infants by using a statistical method for extracting oscillators from time series data, which we call the oscillator decompisition method (OSC-DECOMP). OSC-DECOMP determines the number, frequencies and powers of oscillators as well as their projection patterns in a data-driven manner and does not require arbitrary selection of band-pass filters. Three types of oscillators (brain activity, pulse wave, mirroring noise) were found and each oscillator showed a characteristic spatial synchrony. Model comparison with the Akaike Information Criterion (AIC) demonstrated that the time series of oxy-Hb and deoxy-Hb on each channel originate from common oscillatory activity. We believe that OSC-DECOMP will become a promising tool of neural oscillation analysis for not only fNIRS but also LFP, EEG, MEG and fMRI.

Introduction

Spontaneous changes in cerebral oxygenation reflect hemodynamic response to spontaneous neural activity. The functional Near-Infrared Spectroscopy (fNIRS) has been used to detect the low frequency oscillations of cerebral oxygenation with frequency range between 0.01 and 0.1 Hz in adults [14], elders [5] and infants [6]. Temporal correlation of the oscillatory signals measured at multiple locations revealed the functional connectivity of the cortex in adults [79] and infants [10]. However, there are three outstanding problems. (1) In addition to neurogenic changes, hemodynamic oscillations include cardiovascular oscillations such as cardiac (-1 Hz), respiratory (-0.3 Hz), Mayer wave (-0.1 Hz), and vasomotion (-0.1 Hz) [1, 3, 4, 11, 12]. Conventionally, such oscillation components in fNIRS are extracted by applying band-pass filters to the raw data. (2) Although signals with a wide frequency band (0.01–0.1 Hz) have been used to analyze functional connectivity, subdivision of the frequency range revealed that distinct functional connectivity is dependent on a specific frequency range [13]. It is not clear which frequency range is involved in the functional connectivity in different behavioral states such as wake and sleep, and different populations such as infants and the elderly. (3) The fNIRS can detect both oxygenated and deoxygenated hemoglobin (oxy- and deoxy-Hb) with high temporal resolution (-10 ms). The phase relationship between oxy- and deoxy-Hb changes is referred to as hemoglobin phase of oxygenation and deoxygenation (hPod) and provides rich information about cerebral hemodynamics and metabolism [14]. The phase of the time series has been computed by the Hilbert transform to bandpass-filtered data [6, 14, 15]. However, the selection of the appropriate frequency of the bandpass filter is not clear. It is also not clear whether multivariate time series such as oxy- and deoxy-Hb signals originate from a common oscillator or independent oscillators.

Recently, a statistical method for extracting oscillators from time series data was developed [16, 17], called the oscillator decompisition method (OSC-DECOMP). Here, we briefly demonstrate OSC-DECOMP using figures. Its mathematical detail will be explained in Section 2.3. Fig 1 shows an example of application to univariate time series data. In this case, OSC-DECOMP determines that the given time series (Fig 1A) is a superposition of four oscillation components (Fig 1C) plus noise, based on statistical model fitting. More precisely, each oscillation component in Fig 1C is the projection onto the horizontal axis (first coordinate) of an oscillator that rotates on a plane with random fluctuation (Fig 1B), and OSC-DECOMP extracts an appropriate number of oscillators from the time series. Thus, OSC-DECOMP also provides the phase of each oscillator (Fig 1D). In this way, OSC-DECOMP enables data-driven investigation of the oscillatory dynamics that underlie the time series data. Notably, OSC-DECOMP determines the number and the frequencies of the underlying oscillators in a data-driven (objective) manner and does not require the arbitrary selection of band-pass filters.

thumbnail
Fig 1. OSC-DECOMP for univariate time series data.

(A) Input data. (B) Schematic of OSC-DECOMP. Each oscillation component is the projection onto the horizontal axis (first coordinate) of an oscillator that rotates on a plane with random fluctuation. (C) Decomposition into four oscillation components. (D) Phase of four oscillators (in degrees).

https://doi.org/10.1371/journal.pcbi.1009985.g001

OSC-DECOMP is also applicable to multivariate time series data [17]. Fig 2 shows an example of application to bivariate time series data. In this case, OSC-DECOMP establishes that the given bivariate time series (Fig 2A) originates from three oscillators. As a result, each time series is decomposed into three oscillation components (Fig 2B and 2C), and the phase of the extracted oscillators is also obtained (Fig 2D). The oscillation components in Fig 2B and 2C are the projection of the extracted oscillators. However, unlike the case of the univariate time series (Fig 1B), when applied to multivariate time series, OSC-DECOMP estimates the projection pattern between each time series and each oscillator based on data, and the projection pattern describes the amplitude and phase modulation. Fig 3 shows the estimated projection pattern. In Fig 3A, the thickness of each line represents the amplitude modulation, whereas the color and text of each line indicate the phase modulation. For example, Fig 3B shows the projection pattern for the first oscillator, in which the black and red vectors correspond to the first and second time series, respectively. Given that the red vector is longer than the black vector, the first oscillator superposes on the second time series more than on the first time series. On the other hand, since the angle between the two vectors is 31°, the second time series has a 31° phase delay compared to the first time series. Similarly, Fig 3C and 3D presents the projection pattern for the second and third oscillators. In this way, the projection pattern is estimated in the form of vectors. Thus, OSC-DECOMP enables the investigation of the common oscillators underlying multivariate time series in a data-driven manner.

thumbnail
Fig 2. OSC-DECOMP for bivariate time series data.

(A) Input data. (B) Decomposition into three oscillation components of the first time series. (C) Decomposition into three oscillation components of the second time series. Although the waveform may look almost the same, the scale is different from (B). (D) Phase of extracted oscilators (in degrees).

https://doi.org/10.1371/journal.pcbi.1009985.g002

thumbnail
Fig 3. Estimated projection pattern between each time series and each oscillator.

(A) Bipartite graph representation. The thickness of each line indicates the amplitude modulation. The color and text of each line represent the phase modulation. (B) (C) (D) The projection pattern for three oscillators.

https://doi.org/10.1371/journal.pcbi.1009985.g003

In this study, we investigated oscillatory changes in infant fNIRS data by using OSC-DECOMP. Specifically, we analyzed fNIRS data (3 minutes) acquired from sleeping 3-month-old infants [10] and extracted oscillators on each channel of each infant. The phase difference between oxy- and deoxy-Hb was then calculated for each oscillator. By comparing the Akaike Information Criterion (AIC) of two state space models, we determined that the time series of oxy- and deoxy-Hb on each channel originate from common oscillatory activity. We also utilized the result of OSC-DECOMP to investigate frequency-specific functional connectivity.

Material and methods

Ethics statement

In the present study, we analyze data that was previously reported in [10, 18]. The Office for Life Science Research Ethics and Safety, the University of Tokyo, approved this study (No. 20–225). The parent(s) of all the infants provided written informed consent prior to initiation of experiments.

Participants

We analyzed the data obtained from 21 full-term infants (11 girls and 10 boys; mean postnatal age: 111.6 days, range: 102–123 days) as they slept naturally.

Data acquisition

We used an fNIRS instrument (ETG-7000; Hitachi Medical Corporation) with 94 measurement channels (47 channels in each hemisphere) to detect the relative concentration changes in oxy-Hb and deoxy-Hb [millimolar-millimeter (mM ⋅ mm)] for 3 min without presenting external stimuli. The sampling rate was set at 10 Hz, resulting in 1,800 points per channel of spontaneous changes in oxy- and deoxy-Hb signals. Two sets of 3 × 10 array probes composed of 15 sources and 15 detectors of NIR light (wavelengths: 785 nm and 830 nm, intensity: 0.6 mW), which defined 47 channels, were mounted on a flexible cap over the frontal, temporal, parietal, and occipital regions of the left and right hemispheres. The distance between the source and detector was set to approximately 2.0 cm. The international 10/20 system of electrode placement was referenced to set the positions of the measurement channels. We previously estimated the locations of channels using a magnetic resonance imaging brain atlas [15, 19]. To evaluate the tendency in broad cortical regions and reduce computational complexity, we selected data from nine channels located in the middle line of the array on each hemisphere for the following analyses (Fig 4).

thumbnail
Fig 4. Arrangement of fNIRS measurement channels.

Circles indicate 94 measurement channels. Eighteen circles filled in gray (L1 to L9 and R1 to R9) show the channels used in the present analysis. For reference, some landmarks (black dots) of the international 10/20 system are shown.

https://doi.org/10.1371/journal.pcbi.1009985.g004

Oscillator decompisition method (OSC-DECOMP)

To extract the underlying oscillators from the bivariate time series of oxy- and deoxy-Hb on each channel, we applied the oscillator decompisition method (OSC-DECOMP) [16, 17]. The MATLAB code is available online at http://www.stat.t.u-tokyo.ac.jp/~t-matsuda/software.html. Here, we briefly review this method. See S1 Appendix for background on state space models and time series decomposition.

Let y = (yt,j) (t = 1, …, N; j = 1, …, J) be a J-dimensional time series data of length N with sampling period Δt. The current fNIRS data corresponds to J = 2, N = 1800 and Δt = 0.1 s. OSC-DECOMP assumes that there are K oscillators underlying y, each of which is represented by a point on a two-dimensional plane that rotates around the origin with random fluctuation (Fig 1B). Let be the coordinate of the k-th oscillator at time t. Then, the stochastic dynamics of this oscillator is modeled as (1) where 0 < ak < 1, 0 ≤ 2πfkΔtπ and is an isotropic Gaussian noise with mean zero and variance : Since the 2 × 2 matrix in the right hand side of (1) represents the rotation through an angle 2πfkΔt about the origin, the parameter fk is interpreted as the frequency of the k-th oscillator. The parameter ak is viewed as the regularity of the k-th oscillator in the sense that the waveform of (or ) is closer to the sinusoidal curve for ak closer to one. Also, the parameter specifies the power of the k-th oscillator because the stationary variance of (or ) is . The phase of the k-th oscillator is defined by , where arg(z) is the argument of the complex number z (Fig 1B).

By using the K oscillators, the observed time series y = (yt,j) is modeled as (2) where are independent Gaussian noise. Namely, the observation yt,j is assumed to be the sum of the K inner products of the vectors (cjk,1, cjk,2) with the oscillator coordinates and observation noise wt,j. Recall that the inner product AB of two vectors A and B is given by where ‖A‖ is the length (norm) of A and ϕ is the angle between A and B. Therefore, the observation model (2) means that the k-th oscillator is superposed on the j-th time series with amplitude multiplied by and phase delayed by arg(cjk,1 + icjk,2) (Fig 2E and 2F). In order to ensure parameter identifiability, we fix (c1k,1, c1k,2) = (1, 0) for k = 1, …, K.

The pair (1) and (2) forms a Gaussian linear state space model [20, 21]. Therefore, the posterior distributions of the oscillator coordinates given y1, …, yN are Gaussian and their mean and covariance can be computed by the Kalman smoother algorithm (see S1 Appendix for details). Thus, we apply the Kalman smoother to compute the time series of K oscillators with credible intervals. To extract oscillators in a data-driven manner, the model parameter is estimated by the maximum marginal likelihood: The confidence intervals can be also constructed. See S1 Appendix for details. In addition, the number of oscillators K is determined by minimizing the Akaike Information Criterion (AIC) [22]. Specifically, we fit the state space model (1) and (2) with K = 1, …, Kmax and then select K by where

Testing of common oscillator hypothesis

We fitted the state space model (1) and (2) with J = 2 to the bivariate time series of oxy- and deoxy-Hb on each channel. This model is based on the assumption that the two time series originate from common oscillators. Whereas this assumption seems reasonable for the fNIRS data, there is another possibility that each time series originates from its own oscillators and thus develops independently with each other. Namely, we can consider the state space model (1) and (2) with J = 1 [16] for each time series, which provides another model for the fNIRS data. Thus, we tested the null hypothesis that these two time series models have the same goodness-of-fit to the data. We used the following Linhart-type test [23] to determine the significance of the AIC difference.

A statistical model for time series data y1, …, yN is given by a probability distribution where θ is an unknown parameter to be estimated from data. Consider two candidate models p(1)(y1, …, yNθ(1)) and p(2)(y1, …, yNθ(2)). Let and be the maximum likelihood estimates of θ(1) and θ(2), respectively. Then, the AIC of these models are Let be the log-likelihood of one-step ahead prediction, where and are defined from the stationary distributions of the models. Motivated by Linhart’s test [23], we construct the test statistic as where Under the null hypothesis that the two models have the same goodness-of-fit to the data, the distribution of z is approximately the standard normal N(0, 1). Therefore, the p-value is calculated as where Φ is the cumulative distribution function of the standard normal distribution. Thus, the null hypothesis is rejected if p is smaller than the significance level α. For verification, we also applied the Wilcoxon signed-rank test to and , which does not require the normality assumption on z.

Frequency-specific functional connectivity

Finally, we investigated the frequency-specific functional connectivity based on the result of OSC-DECOMP.

As a result of OSC-DECOMP, we obtained a bivariate time series for each oscillator. To evaluate the correlation between two oscillators, we employed canonical correlation analysis, which is a statistical method for investigating the relationship between two sets of random variables. The (first) canonical correlation coefficient between two random vectors and is defined as where Corr denotes the Pearson correlation coefficient and (a, b) ≠ (0, 0). This problem is reduced to the generalized eigenvalue problem: where ΣXX, ΣXY, ΣYX, ΣYY are submatrices of the covariance matrix of (X, Y), given by The canonical correlation coefficient r is regarded as a measure of the linear dependence between X and Y.

To determine the functional connectivity for target frequency bands, we first selected an oscillator with a frequency fk in the target frequency band on each channel. If there was more than one oscillator in the target frequency band, we selected the oscillator with the maximum power. Next, for each pair of channels, we calculated the (first) canonical correlation coefficient between the selected oscillators. Namely, we substituted the empirical covariance matrix of the oscillators into the preceding formula. We then considered the channel pairs with a canonical correlation coefficient that exceeded the threshold as functionally connected. Thus, we obtained a network of frequency-specific functional connectivity.

Results

Oscillator decomposition of synthetic data

Before presenting the results on real data, we first validate the performance of OSC-DECOMP on synthetic data with comparison to the conventional method with bandpass filtering. For more experiments, see [16, 17].

We generated bivariate time series from the state space model (1) and (2) with J = 2, K = 3, N = 1800, Δt = 0.1s and This parameter setting is motivated from the property of infant fNIRS data found in the next subsection. For this data, OSC-DECOMP correctly detected oscillators by minimizing AIC and their frequencies were Hz, Hz and Hz. Fig 5 summarizes the result.

thumbnail
Fig 5.

(A) Synthetic data. (B) Estimated oscillation components.

https://doi.org/10.1371/journal.pcbi.1009985.g005

Conventionally, oscillation components in time series data are extracted by using the bandpass filters. Here, we compare the accuracy of OSC-DECOMP and bandpass filtering in reconstructing the first oscillation component (time series of the first coordinate of the first oscillator). Note that the second and third oscillators are difficult to separate by the bandpass filter because their spectra significantly overlap. For bandpass filtering, we used the MATLAB function bandpass with passband 0.01–0.1 Hz. The mean squared error in reconstructing the first oscillation component was 5.1 × 10−4 for OSC-DECOMP and 5.6 × 10−3 for bandpass filtering. Thus, OSC-DECOMP extracts the oscillation component more accurately than bandpass filtering.

We also compare the accuracy of OSC-DECOMP and bandpass filtering in estimating the phase difference between two time series (e.g. oxy- and deoxy-Hb). For the same reason as the previous paragraph, we focus on the phase difference for the first oscillator, which is 230°. For OSC-DECOMP, the estimate of the projection vector was and thus the estimate of the phase difference is . In addition, by using the Hessian of the negative log-likelihood at the maximum likelihood estimate (observed Fisher information [24]) and the delta method [25], the 95% confidence interval of the phase difference was obtained as [230.03°, 243.05°]. We used an extension of the Kalman filter algorithm to compute the Hessian [26]. See S1 Appendix for details. For bandpass filtering, the phase difference in 0.01–0.1 Hz at every time point was computed by applying the Hilbert transform to the filtered signals. From these samples of angular variables, the point estimate and 95% confidence interval of the mean direction were computed by using the MATLAB toolbox CircStat [27], which implements the usual procedures for the von Mises distribution [28]. The results were 209.86° and [205.70°, 214.02°]. Thus, OSC-DECOMP provides more accurate point estimates and confidence intervals of the phase difference than bandpass filtering.

Fig 6 shows the noise sensitivity of OSC-DECOMP. It plots the mean squared error of OSC-DECOMP in estimating the first oscillation component with respect to the variance τ2 of observation noise. It demonstrates that the estimation accuracy of OSC-DECOMP is almost constant as long as τ2 is smaller than 10−3.

Oscillator decomposition of infant fNIRS data

To investigate the oscillatory changes in the 18-ch fNIRS data (3 minutes) taken from 21 sleeping 3-month-old infants, we applied OSC-DECOMP to the bivariate time series of oxy- and deoxy-Hb on each channel.

Fig 7 presents one example, in which four oscillators with frequencies 0.028, 1.84, 1.87, 3.74 Hz are extracted. The results are shown for the first 20 seconds. The first oscillator is considered to reflect hemodynamic changes in response to brain activity, whereas the second and third oscillators are interpreted as mirroring noise and cardiac pulse waves, respectively, as will be explained below. The second oscillator superposes on oxy- and deoxy-Hb with almost the same power, whereas the third and fourth oscillators superpose less on deoxy-Hb compared to oxy-Hb. Given that the frequency of the fourth oscillator is twice that of the third oscillator, the fourth oscillator is considered to be the harmonic of the third oscillator. Fig 8 shows the results for the entire 3 minutes and Fig 9 shows the estimated projection pattern. Fig 10 shows the spectrogram and specta. In Fig 10B, in addition to the periodogram of the raw data, the spectral density functions of the fitted oscillators (ARMA(2,1) model) and their sum are also plotted. See [16] for mathematical details. For oxy-Hb, both spectrogram and periodogram have two peaks around 0 Hz and 2 Hz, which correspond well to the first and second/third oscillators, respectively. For deoxy-Hb, only the peak around 0 Hz is clear. Thus, the spectrogram and spectra do not provide enough information to extract common oscillators underlying oxy- and deoxy-Hb time series. In particular, it is difficult to figure out that there are two types of oscillators with different projection patterns around 2 Hz by visual inspection. OSC-DECOMP can separate these oscillators in a data-driven manner by statistical model fitting and model selection with AIC.

thumbnail
Fig 7.

(A) Infant fNIRS (20 seconds). Red: oxy-Hb; blue: deoxy-Hb. (B) Four oscillation components in oxy-Hb. (C) Four oscillation components in deoxy-Hb. (D) Phase of the extracted oscillators.

https://doi.org/10.1371/journal.pcbi.1009985.g007

thumbnail
Fig 8.

(A) Infant fNIRS (3 minutes). Red: oxy-Hb; blue: deoxy-Hb. (B) First oscillation component in oxy- (upper) and deoxy-Hb (lower). (C) Phase of the first oscillator.

https://doi.org/10.1371/journal.pcbi.1009985.g008

thumbnail
Fig 9. Estimated projection pattern for infant fNIRS data in Figs 7 and 8.

https://doi.org/10.1371/journal.pcbi.1009985.g009

thumbnail
Fig 10.

(A) Spectrogram of infant fNIRS data in Fig 8. Upper: oxy-Hb; lower: deoxy-Hb. (B) Periodogram (black), spectral density functions of the fitted oscillators (red) and their sum (blue) for infant fNIRS data in Fig 8. Upper: oxy-Hb; lower: deoxy-Hb.

https://doi.org/10.1371/journal.pcbi.1009985.g010

Fig 11 summarizes the properties of the extracted oscillators for the 21 infants and 18 channels. Fig 11A presents a histogram of the number of extracted oscillators. Five to seven oscillators were extracted in most cases. Fig 11B and 11C shows a histogram of the frequency of the extracted oscillators. There are three peaks near 0.01–0.1 Hz, 1.6–2.4 Hz, and 3.6–4.4 Hz. The first peak is considered to reflect brain activity, whereas the second peak is interpreted as cardiac pulse waves or mirroring noise, as will be explained below. The third peak is regarded as the harmonic of the second peak. Fig 11D and 11E shows scatter plots of the frequency and power of the extracted oscillators. The power tends to be smaller for higher frequencies.

thumbnail
Fig 11.

(A) Histogram of the number of extracted oscillators. (B) (C) Histogram of the frequency of extracted oscillators. (D) (E) Scatter plot of the frequency and power of the extracted oscillators.

https://doi.org/10.1371/journal.pcbi.1009985.g011

Fig 12 presents the projection pattern of the oscillators with frequencies of 0.01–0.1 Hz or 1.6–2.4 Hz. For oscillators with a frequency of 0.01–0.1 Hz, the phase difference between oxy- and deoxy-Hb is approximately 230 degrees (Fig 12A), which is consistent with the findings by [14]. The norm of the projection vector on deoxy-Hb is distributed in the range of 0.5 to 0.7 (Fig 12B), which indicates that these oscillators superpose less on deoxy-Hb compared to oxy-Hb. Since the slow oscillations between 0.01 Hz and 0.1 Hz are likely to reflect hemodynamic changes in response to the spontaneous activity of the brain, they are referred to as “brain oscillators.” Regarding the oscillators with a frequency range of 1.6–2.4 Hz, which coincide with the frequency of cardiac pulses, the phase difference has two peaks at 0 degrees and 180 degrees (Fig 12C). Fig 12D shows that the norm of the projection vector on deoxy-Hb is approximately 0.2 for oscillators with a phase difference of approximately 0 degrees, whereas it is approximately one for oscillators with a phase difference of approximately 180 degrees. Therefore, there are two types of oscillators in the frequency band of 1.6–2.4 Hz. The first type is superposed on oxy- and deoxy-Hb with almost the same phase, and its effect on oxy-Hb is approximately five times larger than that on deoxy-Hb. The second type is superposed on oxy- and deoxy-Hb with almost the same power and opposite phase. The in-phase changes in oxy- and deoxy-Hb can be generated by cardiac pulses because the associated oscillatory changes in the blood volume produce simultaneous changes in both oxy- and deoxy-Hb. In contrast, the anti-phase changes in oxy- and deoxy-Hb can be generated by measurement noise. According to the Lambert–Beer law [29], oxy- and deoxy-Hb changes are calculated as where ΔOD(λ) is the change in optical density measured at a given wavelength λ, and are estimated changes in concentration of oxy- and deoxy-Hb, respectively, εo(λ) and εd(λ) are the extinction coefficients of oxy- and deoxy-Hb at a given wave length λ, respectively, and L is the optical pathlength. However, the detected intensities of the two wavelengths of light include measurement noise in addition to changes due to absorption by chromophores as where and are real changes in concentration of oxy- and deoxy-Hb, respectively, and N(λ) is applied noise at a given wave length λ. Then, we obtain the following equation as For measurements of ΔOD(λ) at two wavelengths the above equation are given by A few lines of algebra indicates that where E = 1/(εo2)εd1) − εo1)εd2))L. This shows that the measurement noise for each wavelength is transformed to noisy signals of oxy- and deoxy-Hb with opposite sign. Thus, the first type is referred to as a “pulse wave oscillator,” and the second type as a “mirroring noise oscillator.” Fig 12E and 12F presents the histogram of the power of these two types of oscillators. Fig 13 summarizes these findings.

thumbnail
Fig 12.

(A) Histogram of phase difference (0.01–0.1 Hz). (B) Scatter plot of phase difference and norm of the projection vector (0.01–0.1 Hz). (C) Histogram of phase difference (1.6–2.4 Hz). (D) Scatter plot of phase difference and norm of projection vector (1.6–2.4 Hz). (E) Histogram of power of pulse wave. (F) Histogram of the power of the mirroring noise.

https://doi.org/10.1371/journal.pcbi.1009985.g012

thumbnail
Fig 13. Projection pattern of three types of oscillators in infant fNIRS.

https://doi.org/10.1371/journal.pcbi.1009985.g013

In the oscillator model (1), the parameter ak ∈ (0, 1) specifies the width of the spectral peak at the frequency fk. Namely, the width of the spectral peak is smaller for ak closer to one. Fig 14 shows the distribution of ak for each type of oscillators. It indicates that the spectral peak is very sharp for the brain oscillators and pulse wave oscillators. In other words, the waveform is close to sinusoidal for these oscillators. It implies that the brain oscillators cannot be interpreted as non-oscillatory drift, although its frequency is close to zero. On the other hand, for the mirroring noise oscillators, the parameter ak is concentrated around 0.4. Thus, the mirroring noise is considered to be non-sinusoidal.

thumbnail
Fig 14.

(A) Scatter plot of the frequency fk and parameter ak of the extracted oscillators. (B) Histogram of the parameter ak for the brain oscillators. (C) Histogram of the parameter ak for the pulse wave oscillators. (D) Histogram of the parameter ak for the mirroring noise oscillators.

https://doi.org/10.1371/journal.pcbi.1009985.g014

Testing of common oscillator hypothesis

To verify that the time series of oxy- and deoxy-Hb on each channel originate from common oscillatory activity, we compared the Akaike Information Criterion (AIC) of two state space models. Fig 15A shows the histogram of the AIC difference, wherein we have 378 data points corresponding to the 18 channels of 21 infants. This shows that the AIC is smaller for the model with common oscillators in all cases. By using the extension of Linhart’s test (see Material and methods), these AIC differences were found to be significant with p < 10−16 in all cases. The Wilcoxon signed-rank test, which does not require the normality assumption on the test statistic, also demonstrated that these AIC differences were significant with p < 10−12 in more than 90% cases. Therefore, the assumption that common oscillators underlie the oxy- and deoxy-Hb time series is indeed supported by the experimental data. For comparison, we applied the same analysis to pairs of oxy- and deoxy-Hb data that were randomly shuffled across the subjects. The result is shown in Fig 15B. In this case, the AIC is larger for the model with common oscillators in all cases except for three.

thumbnail
Fig 15. Histogram of AIC difference for (A) same channel and (B) different channels.

https://doi.org/10.1371/journal.pcbi.1009985.g015

Frequency-specific functional connectivity

Finally, we examined the frequency-specific functional connectivity using the canonical correlation coefficients as explained in Material and Methods. We considered that this method revealed spatial synchrony when oscillators of pulse waves and mirroring noise were analysed. For each of the three oscillators shown in Fig 13, we focused on the oscillator with the highest power. The average number of infants with oscillators was as follows: brain oscillators = 20.2 per channel, pulse wave = 20.7 per channel, and mirroring noise = 14.6 per channel. The number of infants with mirroring noise was lower compared to the other two oscillators because this noise was not always mixed with the fNIRS signal. Fig 16 shows the frequency distribution of the oscillators. For brain oscillators (0.01–0.1 Hz), oscillators with frequencies between 0.02 Hz and 0.03 Hz were most often observed, and the number decreased as the frequency increased (mean: 0.029 Hz, mean frequency range between 18 channels: 0.024–0.037 Hz). Although the oscillators of the mirroring noise (1.6–2.4 Hz, antiphase) exhibited frequencies between 2.1 Hz and 2.2 Hz (mean: 2.075 Hz, mean frequency range: 2.068–2.084 Hz), the distribution of the pulse waves did not have such a peak (mean: 1.993 Hz, mean frequency range: 1.853–2.074 Hz). Friedman’s tests for the brain oscillators did not show significant differences in the frequencies between the channels after Bonferroni correction was performed for multiple comparisons. However, there were significant differences between the two-channel pairs in the tests for pulse waves (L5 and L9, p < 0.05; L5 and R8, p = 0.005) and mirroring noise (L3 and R3: p < 0.05; L9 and R3: p < 0.05). Fig 17 shows the average of the canonical correlation coefficients for the 21 infants. The results for the brain oscillators are presented in Fig 17 (upper). When the threshold is set to 0.7 or 0.6, there are primarily two types of connections: homologous connectivity (inter-hemispheric connectivity between homologous regions) and short-distance ipsilateral connectivity. In addition, when the threshold was set to 0.5, long-distance ipsilateral connectivity between the frontal and temporal/occipital regions was also found. These results are similar to previously reported findings for adult participants [30]. Fig 17 (middle) shows the results for pulse wave oscillators. The correlation is much larger regardless of the channel locations, and most channels are connected even when the threshold is set to 0.7. Fig 17 (lower) shows the results for the mirroring noise oscillators. The correlation between the channels is extremely small, and no connection appears even when the threshold is set to 0.4.

thumbnail
Fig 16. Histogram of the frequency of the oscillators used in the functional connectivity analysis.

https://doi.org/10.1371/journal.pcbi.1009985.g016

thumbnail
Fig 17. Frequency-specific functional connectivity or spatial synchrony.

The numbers at the top represent the correlation coefficients of the thresholds. Upper: 0.01–0.1 Hz, middle: 1.6–2.4 Hz (in-phase), lower: 1.6–2.4 Hz (antiphase).

https://doi.org/10.1371/journal.pcbi.1009985.g017

To evaluate the aforementioned visual inspections, we categorized channel pairs into the following four groups (Fig 18A) based on a previous study [31]: (i) short-range connectivity (16 pairs), (ii) contralateral-transverse connectivity (25 pairs), (iii) ipsilateral-longitudinal connectivity (20 pairs), and (iv) control (92 pairs). Control connectivity consisted of pairs other than (i), (ii), and (iii). Fig 18B shows averaged values of coefficients within each group. The averaged values of the pulse waves and those of the mirroring noise were the highest and the lowest, respectively, and the values of the brain oscillators were intermediate. They were mostly higher than 0.4. This result indicates that the pulse wave is superposed on most channels with similar frequencies and phases, but the presence of mirroring noise in the signal depends on the measurement channel. Regarding brain oscillators, Friedman’s test revealed that short-range and contralateral-transverse connectivities were significantly larger than ipsilateral-longitudinal and control connectivities, which support the aforementioned assertions. The threshold of 0.7 in Fig 17, represents short-range connectivity; the threshold of 0.6 exhibits short-range and contralateral-transverse connectivities; the threshold of 0.5, which shows short-range, contralateral-transverse, and a part of the ipsilateral-longitudinal and control connectivities, and the threshold of 0.4 exhibits most connectivities.

thumbnail
Fig 18. Mean correlation coefficients of each infant for the four types of functional connectivities.

(A) Short range (red), contralateral transverse (green), ipsilateral longitudinal (blue), and control (black) connectivities. (B) Red, green, blue, and black lines indicate individual data for short-range, contralateral transverse, ipsilateral longitudinal, and control connectivity, respectively. **** p < 0.001, *** p < 0.005, ** p < 0.01, and * p < 0.05 after Bonferroni correction for multiple comparisons.

https://doi.org/10.1371/journal.pcbi.1009985.g018

For comparison, Fig 19 shows the functional connectivity computed by the correlation coefficients for oxy-Hb and deoxy-Hb in the frequency range 0.01–0.1 Hz. For both oxy-Hb and deoxy-Hb, short-range and contralateral connectivities are stronger than other types of connectivities. Thus, the network structure is similar to the functional connectivity of brain oscillators in Fig 17. The value of the correlation coefficients is larger in oxy-Hb than deoxy-Hb, which is compatible to previous studies. They tend to take smaller values than the canonical correlation coefficients shown in Fig 17. It may be due to the denoising effect of OSC-DECOMP. It is an interesting future work to explore this more.

thumbnail
Fig 19. Functional connectivity in the frequency range 0.01–0.1 Hz computed by correlation coefficients.

The numbers at the top represent the correlation coefficients of the thresholds. Upper: oxy-Hb, lower: deoxy-Hb.

https://doi.org/10.1371/journal.pcbi.1009985.g019

Comparison with empirical mode decomposition

Here, we compare OSC-DECOMP with the empirical mode decomposition (EMD) [32]. EMD is a method for analyzing nonlinear and nonstationary time series data by decomposing them into several oscillatory components called the “intrinsic mode functions.” By applying the Hilbert transform, the instantaneous frequency of each intrinsic mode function is computed. Thus, the time-frequency distribution of signal amplitude is obtained and it is called the Hilbert spectrum or Hilbert–Huang transform.

Fig 20 shows the intrinsic mode functions (IMF) obtained by EMD (MATLAB function “emd”) for the infant fNIRS data in Fig 7A. In this case, EMD detected seven and eight components in oxy-Hb and deoxy-Hb, respectively. The first IMF seems to correspond to the pulse wave oscillator. The other IMFs have lower frequencies.

thumbnail
Fig 20.

(A) Intrinsic mode functions in the oxy-Hb. (B) Intrinsic mode functions in the deoxy-Hb.

https://doi.org/10.1371/journal.pcbi.1009985.g020

EMD is an algorithm for decomposing univariate nonstationary time series data into several mode functions by focusing on the local properties of the signal. It does not aim at extracting common oscillators in multivariate time series data. On the other hand, OSC-DECOMP extracts oscillators from both univariate and multivariate stationary time series data based on statistical models. It is an interesting future work to extend OSC-DECOMP to nonstationary data.

Discussion

In this study, we investigated the oscillatory activity in infant fNIRS data by using the oscillator decompisition method (OSC-DECOMP) [16, 17]. Five to seven oscillators were extracted on most channels and their frequency distribution had three peaks at 0.01–0.1 Hz, 1.6–2.4 Hz and 3.6–4.4 Hz. The first peak was considered to reflect brain activity and the phase difference between oxy- and deoxy-Hb for the oscillators were at approximately 230 degrees. The second peak involved two types of in-phase and anti-phase oscillations of oxy- and deoxy-Hb. The former and latter are attributed to cardiac pulse waves and mirroring noise. The third peak was regarded as the harmonic of the second peak. By comparing the Akaike Information Criterion (AIC) of the two state space models, we verified that the time series of oxy- and deoxy-Hb on each channel originate from common oscillatory activity. We also applied the results of OSC-DECOMP to investigate the frequency-specific functional connectivity.

Advantages of OSC-DECOMP

Neural oscillations are relevant in many brain functions [33]. For example, the electroencephalogram (EEG) time series is composed of several oscillation components (e.g., alpha, beta, and gamma), and each oscillation component has its physiological role. Conventionally, neural oscillation analysis is conducted by applying band-pass filters to the time series data [34]. The outputs of the band-pass filters are assumed to represent the target oscillation components. For example, the output of a filter with a pass-band of 8–13 Hz is considered as the alpha component. Although such conventional methods are widely used and easy to implement, there are several associated limitations [35]. In particular, the band-pass filter is often selected arbitrarily, and this affects the final result. For example, even the definition of the alpha band seems to vary among studies, such as 8–13 Hz and 9–12 Hz.

In contrast, OSC-DECOMP extracts the oscillation components in a data-driven (objective) manner using a statistical modeling approach. Both the number and frequencies of underlying oscillators are determined based on data. When applied to multivariate time series data, this method also estimates the projection pattern between each time series and each oscillator, which describes the amplitude and phase modulation. Thus, even if there are oscillators with very close frequencies, OSC-DECOMP can separate them correctly if their projection patterns are sufficiently different. The current result on pulse wave and mirroring noise can be interpreted as a typical example of this. Therefore, OSC-DECOMP enables the investigation of the common oscillators that underly multivariate time series in a data-driven manner. Several studies have focused on the phase difference between the band-pass filtered signals of oxy- and deoxy-Hb. For example, [14] investigated the phase difference at a frequency of 0.05–0.1 Hz in infant fNIRS data and discussed its relationship with preterm birth. Our findings may provide a foundation for such investigations of phase differences. It should also be noted that the framework of state-space models naturally leads to real-time phase estimation [36].

Oscillators underlying infant fNIRS data

The present study based on OSC-DECOMP demonstrated that the brain, pulse wave, and mirroring noise oscillators can be extracted from the fNIRS signals obtained from the cranium of sleeping infants. In the time domain, fNIRS signals with slow components between 0.01 Hz and 0.1 Hz generally include hemodynamic and oxygenation changes in response to brain activity, similarly to fMRI [3739]. The present study revealed that oscillators with frequencies between 0.02 Hz and 0.03 Hz were most common, as shown in Fig 11B and 11D. Furthermore, the phase difference between oxy- and deoxy-Hb chambers was approximately 230 degrees, as shown in Fig 12A and 12B, which reflects complex mechanisms for hemodynamics and oxygenation [14]. This method can serve as a new approach for studying brain activation and connectivity using the decomposed oscillators of fNIRS signals. In addition to the changes induced by brain activity, the slow components of fNIRS signals are subject to systemic influences such as heart rate and blood pressure [40]. Although it was not proven that the brain oscillator was purely neurogenic in the present study, data-driven decomposition of the oscillator will simplify the procedure for determining if some oscillation components are associated with systemic oscillations.

For frequencies higher than 0.1 Hz, fNIRS signals include components associated with respiration (-0.3 Hz) and heartbeat (-1 Hz), in the case of adults [41]. In the present study, although distinct oscillators corresponding to respiration were not observed, the conspicuous oscillators caused by cardiac pulse waves were extracted at approximately 2 Hz, as shown in Fig 11C and 11E. The oscillators included in-phase and anti-phase oscillations of oxy- and deoxy-Hb, as shown in Fig 12C and 12D. The in-phase changes in oxy- and deoxy-Hb data can be attributed to blood volume changes induced by pulse waves [11, 42]. The anti-phase changes in oxy- and deoxy-Hb can be attributed to mirroring noise that is artificially produced during the calculation of the concentration of two chromophores using two-wavelength absorption based on the Lambert-Beer equation [29]. The measurement noise in the intensities of the detected light at each wavelength is transformed to noisy signals of oxy- and deoxy-Hb with opposite signs. The mirroring changes in oxy- and deoxy-Hb can also be caused by insufficient separation of chromophore changes due to the wavelength dependence of the light path length [43, 44]. As we assumed the same optical path length for the two wavelengths, the estimated signals associated with oxy- and deoy-Hb may have cross-talk. This cross-talk is not confined to the higher frequency range, but can also occur in the lower frequency range. Moreover, a small number of oscillators with a phase difference of 180° was observed in the slower frequency range, as shown in Fig 12B.

fNIRS signals are assumed to reflect changes in hemodynamics and oxygenation in cerebral tissue. However, those that originate from extracerebral tissue can have an impact on the fNIRS signals [45]. In the present study, we used data acquired from sleeping 3-month-old infants using a source-detector distance of 2 cm. A previous study involving infants of the same age showed that 2 cm is the optimal inter-optode distance for the detection of hemodynamic changes in response to auditory stimulation [46]. Based on the deep-shallow signal separation method using multi-distant probes and independent component analysis, it was further revealed that the deep tissue contribution during sleeping with a 2 cm separation was 66–79% for oxy-Hb changes [47]. Thus, the brain oscillator in the present study primarily reflects the hemodynamics and oxygenation of cerebral tissue. To gain more detailed information exclusive to brain activation by combining decomposition methods of fNIRS signals in the time domain of oscillations and in the spatial domain of tissue depth, further studies are required.

Frequency-specific functional connectivity in fNIRS and future studies

In this study, we examined frequency-specific functional connectivity by calculating the canonical correlation between two oscillators in different channels (Figs 17 and 18). We found that each of the three types of oscillators (i.e., brain oscillators, pulse waves, and mirroring noise) showed individual patterns of functional connectivity or spatial synchrony. Conventionally, oscillator analysis and functional connectivity of fNIRS signals utilize the coherence between band-pass filtered signals [3, 10, 13, 30, 48]. For example, [30] revealed that the connectivity pattern depends on the frequency band, and short-range and long-range connectivities are manifested in oxy-Hb signals at frequencies below 0.1 Hz. It would be interesting to investigate the relationship between the proposed oscillator-based method and the conventional coherence-based method. The most obvious difference between the two methods is the processing of the two types of Hb signals. The conventional method treats the two types of Hb signals separately, but the proposed method can treat them individually or together. The deoxy-Hb signals were smaller than oxy-Hb signals, and the measurement noise had a greater effect on the coherence between the channels in the deoxy-Hb signals compared to the oxy-Hb signals. As a result, there was a tendency to rely on oxy-Hb signals when summarizing the results obtained using the conventional method. Given that the present method uses common oscillators for the oxy-Hb and deoxy-Hb signals, it has the advantage that the relationship between the channels can be investigated without limiting the oxy-Hb signals. Moreover, Welch’s average periodogram method is often used for calculating coherence, which involves the subjective selection of the window type, window length, and overlap length. In contrast, the proposed method does not require such selections. This reduces the risk that the analysis method will eventually bias the results.

The short-range and contralateral-transverse connectivity of the brain oscillators were higher than the ipsilateral-longitudinal connectivity and the control (Fig 18). These results are consistent with those of previous studies on infants [31] and adults [30]. The short-range and contralateral-transverse connectivity reflect functional relationships based on structural connectivity, such as U-fibers beneath the cortical surface and callosal fibers. Although there are longitudinal fasciculi that connect the anterior and posterior regions, even in infancy [49, 50], the presented and previously reported results consistently showed relatively low synchronization between distant ipsilateral regions. As almost all infants exhibited brain oscillators and the mean frequency of the oscillators was not significantly different between channels, it is possible that the phases of the brain oscillators were not the same in the anteroposterior direction, and the differences caused lower ipsilateral-longitudinal connectivity. This possibility can be examined if we used a longer time series of fNIRS data in a future study.

In the present study, we first determined the brain oscillators for each measurement channel. The functional connectivity was then evaluated by calculating the correlation between the time series of the oscillators. Another way of defining functional connectivity between channels is to determine whether the change in the oxy-Hb and deoxy-Hb signals in two or more channels is caused by common oscillatory activity. For example, each of the two channels measures a set of time series of oxy-Hb and deoxy-Hb signals, and a total of four time series can be used to investigate the presence of a common oscillator. If there is a low-frequency oscillator with power above a certain threshold, these two channels can be considered to exhibit significant functional connectivity. A functional network can also be defined by establishing whether there is a common oscillator among multiple channels. This method can be applied to time series of functional magnetic resonance imaging (fMRI), magnetoencephalography (MEG), electrocorticography (ECoG), and multi-electrode EEG to clarify functional relationships across broad brain regions.

Supporting information

S1 Appendix. Technical appendix.

The supplementary material includes the technical details on (1) state space models and time series decomposition, (2) Kalman filter/smoother algorithm, and (3) Hessian computation and confidence interval.

https://doi.org/10.1371/journal.pcbi.1009985.s001

(PDF)

References

  1. 1. Elwell C. E., Springett R., Hillman E. & Delpy D. T. (1999). Oscillations in cerebral haemodynamics. Implications for functional activation studies. Advances in Experimental Medicine and Biology 471, 57–65. pmid:10659132
  2. 2. Hoshi Y., Kosaka S., Xie Y., Kohri S. & Tamura M. (1998). Relationship between fluctuations in the cerebral hemoglobin oxygenation state and neuronal activity under resting conditions in man. Neuroscience Letters 245, 147–150. pmid:9605477
  3. 3. Obrig H., Neufang M., Wenzel R., Kohl M., Steinbrink J., Einhaupl K. et al. (2000). Spontaneous low frequency oscillations of cerebral hemodynamics and metabolism in human adults. NeuroImage 12, 623–639. pmid:11112395
  4. 4. Toronov V., Franceschni M.A., Filiaci M., Fantini S., Wolf M., Michalos A. et al. (2000). Near-infrared study of fluctuations in cerebral hemodynamics during rest and motor stimulation: temporal analysis and spatial mapping. Medical Physics 27, 801–815. pmid:10798703
  5. 5. Schroeter M.L., Schmiedel O., Watanabe H., Sasaki A. T. & von Cramon D.Y. (2004). Spontaneous low-frequency oscillations decline in the aging brain. Journal of Cerebral Blood Flow & Metabolism 24, 1183–1191. pmid:15529019
  6. 6. Taga G., Konishi Y., Maki A., Tachibana T., Fujiwara M. & Koizumi H. (2000). Spontaneous oscillation of oxy- and deoxy- hemoglobin changes with a phase difference throughout the occipital cortex of newborn infants observed using non-invasive optical topography. Neuroscience Letters 282, 101–104. pmid:10713406
  7. 7. Lu C. M., Zhang Y. J., Biswal B. B., Zang Y. F., Peng D. L. & Zhu C. Z. (2010). Use of fNIRS to assess resting state functional connectivity. Journal of Neuroscience Methods 186, 242–249. pmid:19931310
  8. 8. White B.R., Snyder A.Z., Cohen A.L., Petersen S.E., Raichle M.E., Schlaggar B.L. et al. (2009). Resting-state functional connectivity in the human brain revealed with diffuse optical tomography. NeuroImage 47, 148–156. pmid:19344773
  9. 9. Zhang H., Zhang Y., Lu C., Ma S., Zang Y. & Zhu C. (2010). Functional connectivity as revealed by independent component analysis of resting-state fNIRS measurements. NeuroImage 51, 1150–1161. pmid:20211741
  10. 10. Homae F., Watanabe H., Otobe T., Nakano T., Go T., Konishi Y. et al. (2010). Development of global cortical networks in early infancy. Journal of Neuroscience 30, 4877–4882. pmid:20371807
  11. 11. Tong Y., Lindsey K. P. & deB Frederick B. (2011). Partitioning of physiological noise signals in the brain with concurrent near-infrared spectroscopy and fMRI. Journal of Cerebral Blood Flow & Metabolism 31, 2352–2362. pmid:21811288
  12. 12. Yücel M. A., Selb J., Aasted C. M., Lin P. Y., Borsook D., Becerra L. et al. (2016). Mayer waves reduce the accuracy of estimated hemodynamic response functions in functional near-infrared spectroscopy. Biomedical Optics Express 7, 3078–3088. pmid:27570699
  13. 13. Sasai S., Homae F., Watanabe H., Sasaki A. T., Tanabe H. C., Sadato N. et al. (2011). A NIRS-fMRI study of resting state network. NeuroImage 63, 179–193.
  14. 14. Watanabe H., Shitara Y., Aoki Y., Inoue T., Tsuchida N., Takahashi N. et al. (2017). Hemoglobin phase of oxygenation and deoxygenation in early brain development measured using fNIRS. Proceedings of the National Academy of Sciences 114, E1737–E1744. pmid:28196885
  15. 15. Taga G., Watanabe H. & Homae F. (2018). Spatial variation in the hemoglobin phase of oxygenation and deoxygenation in the developing cortex of infants. Neurophotonics 5, 011017. pmid:29021987
  16. 16. Matsuda T. & Komaki F. (2017). Time series decomposition into oscillation components and phase estimation. Neural Computation 29, 332–367. pmid:27870609
  17. 17. Matsuda T. & Komaki F. (2017). Multivariate time series decomposition into oscillation components. Neural Computation 29, 2055–2075. pmid:28562213
  18. 18. Homae F., Watanabe H., Nakano T. & Taga G. (2011). Large-scale brain networks underlying language acquisition in early infancy. Frontiers in Psychology 2, 93. pmid:21687461
  19. 19. Matsui M., Homae F., Tsuzuki D., Watanabe H., Katagiri M., Uda S. et al. (2014). Referential framework for transcranial anatomical correspondence for fNIRS based on manually traced sulci and gyri of an infant brain. Neuroscience Research 80, 55–68. pmid:24445146
  20. 20. Durbin J. & Koopman S. J. (2012) Time Series Analysis by State Space Methods. Oxford: Oxford University Press.
  21. 21. Kitagawa G. (2010). Introduction to Time Series Modeling. London: Chapman & Hall.
  22. 22. Akaike H. (1980). Likelihood and the Bayes procedure. In Bayesian Statistics, eds. Bernardo J. M. DeGroot M. H. Lindley D. V. and Smith A. F. M. Valencia, Spain, 1–13.
  23. 23. Linhart H. (1988). A test whether two AIC’s differ significantly. South African Statistical Journal 22, 153–161.
  24. 24. Efron B. & Hinkley D. V. (1978). Assessing the accuracy of the maximum likelihood estimator: Observed versus expected Fisher information. Biometrika 65, 457–483.
  25. 25. van der Vaart A. W. (1998). Asymptotic Statistics. Cambridge University Press.
  26. 26. Kitagawa, G. (2020). Computation of the gradient and Hessian of the log-likelihood of the state-space model by the Kalman filter. arXiv:2011.09638.
  27. 27. Berens P. (2009). CircStat: A MATLAB toolboc for circular statistics. Journal of Statistical Software 31, 1–21.
  28. 28. Mardia K. V. & Jupp P. E. (2000). Directional Statistics. New York: John Wiley & Sons.
  29. 29. Villringer A. & Chance B. (1997). Non-invasive optical spectroscopy and imaging of human brain function. Trends in Neurosciences 20, 435–442. pmid:9347608
  30. 30. Sasai S., Homae F., Watanabe H. & Taga G. (2011). Frequency-specific functional connectivity in the brain during resting state revealed by NIRS. NeuroImage 56, 252–257. pmid:21211570
  31. 31. Imai M., Watanabe H., Yasui K., Kimura Y., Shitara Y., Tsuchida S. et al. (2014). Functional connectivity of the cortex of term and preterm infants and infants with Down’s syndrome. NeuroImage 85, 272–278. pmid:23631984
  32. 32. Huang N. E., Shen Z., Long S. R., Wu M. C., Shih H. H., Zheng Q. et al. (2014). The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis. Proceedings of the Royal Society of London. Series A 454, 903–995.
  33. 33. Buzsáki G. (2011). Rhythms of the Brain. Oxford: Oxford University Press.
  34. 34. Cohen M. X. (2014). Analyzing neural time series data. Cambridge: MIT Press.
  35. 35. Siapas A. G., Lubenov E. V. & Wilson M. A. (2005). Prefrontal phase locking to hippocampal theta oscillations. Neuron 46, 141–151. pmid:15820700
  36. 36. Wodeyar A., Schatza M., Widge A. S., Eden U. T., & Kramer M. A. (2021). A State Space Modeling Approach to Real-Time Phase Estimation. bioRxiv. pmid:34569936
  37. 37. Hathout G. M., Gopi R. K., Bandettini P. & Gambhir S. S. (1999). The lag of cerebral hemodynamics with rapidly alternating periodic stimulation: modeling for functional MRI. Magnetic resonance imaging 17, 9–20. pmid:9888394
  38. 38. Snyder A. Z. & Raichle M. E. (2012). A brief history of the resting state: the Washington University perspective. NeuroImage 62, 902–910. pmid:22266172
  39. 39. Mitra A & Raichle M. E. (2016). How networks communicate: propagation patterns in spontaneous brain activity. Philosophical Transactions of the Royal Society B: Biological Sciences 371, 20150546. pmid:27574315
  40. 40. Katura T., Tanaka N., Obata A., Sato H. & Maki A. (2006). Quantitative evaluation of interrelations between spontaneous low-frequency oscillations in cerebral hemodynamics and systemic cardiovascular dynamics. NeuroImage 31, 1592–1600. pmid:16549367
  41. 41. Pierro M. L., Hallacoglu B., Sassaroli A., Kainerstorfer J. M. & Fantini S. (2014). Validation of a novel hemodynamic model for coherent hemodynamics spectroscopy (CHS) and functional brain studies with fNIRS and fMRI. NeuroImage 85, 222–233. pmid:23562703
  42. 42. Themelis G., D’Arceuil H. E., Diamond S. G., Thaker S., Huppert T. J., Boas D. A. et al. (2007). Near-infrared spectroscopy measurement of the pulsatile component of cerebral blood flow and volume from arterial oscillations. Journal of biomedical optics 12, 014033. pmid:17343508
  43. 43. Strangman G., Franceschini M. A. & Boas D. A. (2003). Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for focal changes in oxygenation parameters. NeuroImage 18, 865–879. pmid:12725763
  44. 44. Uludag K., Steinbrink J., Villringer A. & Obrig H. (2004). Separability and cross talk: optimizing dual wavelength combinations for near-infrared spectroscopy of the adult head. NeuroImage 22, 583–589. pmid:15193586
  45. 45. Kirilina E., Jelzow A., Heine A., Niessing M., Wabnitz H., Brühl R. et al. (2012). The physiological origin of task-evoked systemic artefacts in functional near infrared spectroscopy. NeuroImage 61, 70–81. pmid:22426347
  46. 46. Taga G., Homae F. & Watanabe H. (2007). Effects of source-detector distance of near infrared spectroscopy on the measurement of the cortical hemodynamic response in infants. NeuroImage 38, 452–460. pmid:17884584
  47. 47. Funane T., Homae F., Watanabe H., Kiguchi M. & Taga G. (2014). Greater contribution of cerebral than extracerebral hemodynamics to near-infrared spectroscopy signals for functional activation and resting-state connectivity in infants. Neurophotonics 1, 025003. pmid:26157977
  48. 48. Sakakibara E., Homae F., Kawasaki S., Nishimura Y., Takizawa R., Koike S. et al. (2016). Detection of resting state functional connectivity using partial correlation analysis: A study using multi-distance and whole-head probe near-infrared spectroscopy. NeuroImage 142, 590–601. pmid:27521742
  49. 49. Dubois J., Dehaene-Lambertz G., Kulikova S., Poupon C., Hüppi P. S., & Hertz-Pannier L. (2014). The early development of brain white matter: a review of imaging studies in fetuses, newborns and infants. Neuroscience 276, 48–71. pmid:24378955
  50. 50. Feng L., Li H., Oishi K., Mishra V., Song L., Peng Q. et al. (2019). Age-specific gray and white matter DTI atlas for human brain at 33, 36 and 39 postmenstrual weeks. NeuroImage 185, 685–698. pmid:29959046