2.1.

Simulation Analysis

Studies in photon migration in tissues have shown that photons follow a somewhat banana shape trajectory between a source and detector pair. Placing the source–detector with a certain distance apart, we can probe the tissues with a particular depth sensitivity profile similar to a banana shape, as shown as in Fig. 1.

Fig. 1

Geometry of the simulation set up. S, source; D, detector; L1 stands for the skin layer of the head; L2 for skull and cerebrospinal fluid (CSF); and L3 for cortex. Vasculature in L2 has been ignored since it contains negligible amounts of [Hb] and [HbO] variations. The tubes represent the veins running parallel to the source–detector axis. Optical path is illustrated as the semitransparent banana shaped figure.

A signal received by a photodiode includes the changes in the absorption within the underlying tissues. In the near-infrared light spectrum, the highest absorbers are the [Hb] and [HbO] and a spectroscopic measurement of the absorption will yield the concentration changes of these chromophores, which can be calculated from Beer–Lambert’s law, explained in detail in other studies.^1,35–37 Since absorption change within a layer of tissue $l$ at a specific wavelength $\lambda$ is formulated as $\mathrm{\Delta }{\mu }_l^{\lambda }={ϵ}_C^{\lambda }\mathrm{\Delta }{\mathbf{C}}_l$, where ${ϵ}_C^{\lambda }=[{ϵ}_{\mathrm{Hb}}^{\lambda }{ϵ}_{\mathrm{HbO}}^{\lambda }]$ and $\mathrm{\Delta }{\mathbf{C}}_l={(\mathrm{\Delta }{[\mathrm{Hb}]}_l\mathrm{\Delta }{[\mathrm{HbO}]}_l)}^T$ for fNIRS studies, we can define a signal model for the [Hb] and [HbO] concentrations for any $l$ layer in the form of ${\mathrm{Hb}}_l$ and ${\mathrm{HbO}}_l$.

2.1.1.

Signal model

Let us assume that the hemodynamic changes [$h_l^C(t)$] at each layer ($l$) for given chromophore ($C$) can be modeled as a sum of independent signal activity ($\mathbf{S}$) weighted with layer specific weights ($\mathbf{W}$), as shown below:

Eq. (5)

$$h_l^C(t)=\mathbf{WS},$$

where the entries of the matrix $\mathbf{W}$, $w_{l,k}$ are the weights of a specific hemodynamic activity $s_k(t)$ at that specific layer $l$ that are the entries of the $\mathbf{S}$ matrix. For the sake of simplicity, we will assign three signal activities to be present at each layer: $s^b(t)$, brain hemodynamic response function (BRHF); $s^s(t)$, task irrelevant systemic physiological hemodynamic fluctuations; and $s_n(t)$, uncorrelated instrumentation noise, modeled as $\mathcal{N}(0,{\sigma }_n)$. Hence, $\mathbf{S}=[s^b(t)s^s(t)s_n(t)]$. Naturally, the weight of $s^b(t)$ for the first and second layers will be zero (i.e., $w_1^b=w_2^b=0$, $w_3^b=1$) while $0.1\le w_1^s=w_2^s=w_3^s\le 1$.

Traditionally, $s^b(t)$ and $s^s(t)$ are defined as follows:

Eq. (6)

$$s^b(t-{\theta }_i^b)={(t-{\theta }_i^b)}^2{e}^{-\left[\frac{(t-{\theta }_i^b)}{\tau }\right]},$$

Eq. (7)

$$s^s(t-{\theta }_i^s)=\left\{\sin \right[2\pi \frac{f_1}{f_s}(t-{\theta }_i^s)]+\sin [2\pi \frac{f_2}{f_s}(t-{\theta }_i^s)\left]\right\},$$

where $s^b(t-{\theta }_i^b)$ is the brain BHRF modeled as a gamma function with a delay of ${\theta }_i^b$, a linearly increasing value for each channel ranging from 3 to 10 s to assure a variance in the correlation values between each channel;^38,39 $s^s(t-{\theta }_i^s)$ is the systemic fluctuations, typically, the Mayer’s wave (resonance frequency of $f_1$) and the breathing-related hemodynamic fluctuations (resonance frequency of $f_2$), with a randomized delay of ${\theta }_i^s$ for each channel and $f_s$ is the sampling rate. Simulated signals for several channels and the correlation between them can be seen in Fig. 2(a). Figure 2(b) shows the original BHRF and the mixture signal from Eq. (5).

Fig. 2

(a) Simulated HRF signals for different channels with a certain delay. (b) A mixed signal simulation for channel 1 based on Eq. (5) with $s^b(t)$ in green, the mixed signal $h_l(t)$ in blue, and subtracted filtered signal $\widehat{x}(t)$ in red. The parameters of the simulation are $w_3^b=1$, $\tau =1$, $w_3^s=0.1$, $f_1=0.1\mathrm{Hz}$, $f_2=0.25\mathrm{Hz}$, $f_s=2\mathrm{Hz}$, and $s_n(t)=\mathcal{N}(0,0.01)$.

In all the simulations, $w_3^b=1$ while $w_1^s=w_2^s=w_3^s$ and $w_3^s=\{0.1,0.2,\dots ,1\}$.

2.1.2.

Functional connectivity analysis

We have decided to use a rectangular probe geometry to simulate the fNIRS signals, as shown in Fig. 3.

Fig. 3

Rectangular probe geometry used in the simulations. The sources are the red circles while the detectors are the black squares. The source–detector separation is fixed at 2.5 cm. Probe is placed so that channel 1 is on the left, midline of the probe is aligned with the midline of the forehead, and the bottom row is right above the eyebrows. A detailed placement and anatomic localization of the probe can be seen in Ref. 4.

We decided to compute the FC using a PC-based analysis. PC provides a relationship between two variables after removing the overlap from both variables. The diagram in Fig. 4 depicts PC between time series 1 and 2 in the presence of a third time series 3. The PC coefficient between 1 and 2 after removing the influence of 3 ($r_{\mathrm{1,2}|3}$) is as follows:²⁰

Eq. (8)

$$r_{\mathrm{1,2}|3}=\frac{r_{\mathrm{1,2}}-r_{\mathrm{1,3}}{r}_{\mathrm{2,3}}}{\sqrt{(1-r_{\mathrm{1,3}}^2)}\sqrt{(1-r_{\mathrm{2,3}}^2)}}.$$

Fig. 4

The shaded region is a graphical representation of the PC between 1 and 2 in the presence of 3.

This equation can be generalized to compute PC between any two channels $(i,j)$ as $r_{i,j|k}$ in the presence of a common influencer ($k$).

Since the fNIRS signal model assumes a linear addition of the BHRF with systemic fluctuations and noise, we assumed a certain frequency band for the brain-related signals and systemic fluctuations.

2.2.

fNIRS Data

The data were collected in an earlier study.^24,40 Twelve healthy controlled subjects performed the computerized version of the color-word matching Stroop task.⁴¹ The subjects were asked to respond to 90 stimuli presented on a screen every 4 s, in groups of six stimuli within each block. Fifteen blocks were divided into five NS, five CS, and five IcS type of stimuli and presented in a random fashion. There were 20 s of rest within each block. The subjects were asked to respond to stimuli by pressing either the left or the right button of the mouse based on a match or unmatch condition. fNIRS data were collected with a 16 channel dual wavelength continuous wave system with a sampling rate of 1.7 Hz (ARGES Cerebro, Hemosoft Inc., Ankara, Turkey). The source–detector configuration is shown in Fig. 3 with a separation of 2.5 cm. Absorption data collected at each detector are converted to [Hb] and [HbO] via the modified Beer–Lambert’s Law. The protocol was approved by the Ethical Review Committee of Bogazici University.

3.1.

Simulation Analysis

The simulated signals for various weights as in Eq. (5) were used to compute the FC matrices. The plots in Fig. 5(a) show a sample of such signals while the errors in estimating the FC matrices are given in Fig. 5(b).

Fig. 5

(a) Simulated HRF signals for different channels with a certain delay. A mixed signal simulation for channel 1 based on Eq. (5) with $s^b(t)$ in red and the mixed signal $h_l(t)$ in blue. The parameters of the simulation are $w_s=0.3$, $\tau =1$, $f_1=0.1\mathrm{Hz}$, $f_2=0.25\mathrm{Hz}$, $f_s=2\mathrm{Hz}$, and $s_n(t)=\mathcal{N}(0,0.03)$. The first plot in (a) shows the raw HbO signal, second row Butterworth low pass filtered signals and last plot shows the averaged regressor signal ${\mathrm{HbO}}_R$ after being passed through a Butterworth type high pass filter, (b) errors in estimating the FC matrices. $E_{\mathrm{CO}}$ is the error with respect to Pearson’s correlation coefficient, $E_{\mathrm{BW}}$ is the error after low pass filtering with a Butterworth filter, and $E_{\mathrm{PC}}$ is the error with respect to PC analysis. The numbers above the bars indicate the percent improvement in the accuracy (decrease in the errors) between the low pass filtered and PC-based connectivity matrices.

The top plot in Fig. 5(a) depicts how difficult it is to observe the presence of the hemodynamic response from the raw [HbO] signal. Figure 5(b) shows how well the PC-based computation of the correlations is closer to real correlation values. As the weight of the systemic fluctuations increases and starts to dominate the whole signal, the accuracy of extracting the real correlation value decreases.

3.2.

Real Data Analysis

A sample of fNIRS data from subject 1 and the corresponding regressor signal (${\overline{\mathrm{HBO}}}_R$) can be seen in Figs. 6(a) and 6(b). The high pass filter setting was set at 0.09 Hz.^49–51 Once the regressor is obtained, the data are segmented into NS, CS, and IcS parts. The choice of these cutoff frequencies was based on the fact that Mayer’s wave is centered around 0.1 Hz with a slight variation from 0.09 to 0.11 Hz.^52–55 So the choice for the cutoff frequency for the high pass filter was based on the lower end of the Mayer’s wave band. Note that no further filtering was applied to the raw data when computing the PCs. Once the regressor is formed from the raw data, then this regressor was used in computing the PC pairs of channels of unfiltered fNIRS data. This way any smearing effect of the spikey artifacts due to low pass filtering was minimized.

Fig. 6

(a) A representative ${\mathrm{HBO}}_2$ data from channel 4 of subject 1 is low pass filtered with a fourth order Butterworth filter with a cutoff frequency at 0.08 Hz.⁴⁹ NS, CS, and IcS episodes are marked with varying colors. (b) The regressor for subject one is obtained by averaging the high pass filtered signal at 0.09 Hz of each channel. (c) Average of $\mathcal{FC}$ matrices for IcS calculated with Pearson’s correlation to low pass filtered data ($n=12$ subjects, $\overline{\mathcal{F}{\mathcal{C}}_{\mathcal{I}}}$) and (d) $\overline{\mathcal{F}{\mathcal{C}}_{\mathcal{I}}}$ matrix calculated by PC.

3.2.2.

GE values for real data

The FC matrices are usually thresholded at a certain cutoff level to leave only a percent of the strongest connections before the computation of GE values. A scan of threshold values leaving the strongest 15 to 25 values yielded the highest significance among three different stimuli at GE values for [HbO] FC matrices ($\mathrm{TH}=21$) and $\mathrm{TH}=24$ for [Hb]. Average of the thresholded [HbO] FC matrices computed after thresholding at these values is shown in Figs. 7(a)–7(d), while their average correlation values are shown in the last column of Table 1. The average of the interference matrices ($\overline{\mathcal{F}{\mathcal{C}}_{\mathcal{I}}-\mathcal{F}{\mathcal{C}}_{\mathcal{N}}}$) in Fig. 7(d) shows an increase of the correlation values in the right dorsolateral prefrontal cortex (dlPFC), which is in line with several of the fNIRS findings.

Fig. 7

Average of FC matrices over all subjects for (a) NS, (b) CS, (c) IcS, and (d) interference for IcS-NS are shown. Note that the interference values are in the range of −0.2 to 0.2.

Since not every subject had high correlations for the same channel pairs, the average values are lower for the averaged FC matrices [note the highest value in the colorbar is 0.6 for Figs. 7(a)–7(c)]. The average of thresholded FC matrices shows a clustering trend toward the right dlPFC (channels 9 to 16) as the stimulus became cognitively more demanding. Figure 7(d) is the difference between Figs. 7(a) and 7(c), and elucidates that the correlation values on the right dlPFC increases in strength (an average increase of 0.0028) while a consecutive decrease (an average of $-0.0056$) is observed for the left dlPFC (channels 1 to 8).

The GE values computed from the FC matrices generated by regular correlation (Pearson’s correlation coefficient) of low pass filtered data (via Butterworth filter) did not show any significant differences for various types of stimuli ($p>0.05$). When the correlations are computed via the PC analysis, we observed an increase in the GE values as the cognitive task became more demanding. GE values apparently change with respect to the threshold used for the FC matrices. We swept the threshold values (TH) and observed the significance among the ${\mathrm{GE}}_{\mathrm{N}}$, ${\mathrm{GE}}_{\mathrm{C}}$, and ${\mathrm{GE}}_{\mathrm{Ic}}$ both for [Hb] and [HbO] values. The TH value shows that the highest significance was observed for $\mathrm{TH}=21$ for [HbO], which corresponds to 8.75% of the highest correlations when the diagonals are omitted (21/240) and $\mathrm{TH}=24$ for [Hb] data. In a study by Zhang et al.,³⁰ a strong lateralization effect was observed, favoring the flow of information to the right side. We grouped the detectors into four areas, where L (left) corresponds to the GE computed for detectors from the FC matrices 1 through 8 [$i=1\dots 8$, $j=1\dots 8$ in Eq. (9)], R (right) for detectors from 9 through 16 [$i=9\dots 16$, $j=9\dots 16$ in Eq. (9)], IH corresponds to interhemispheric connectivity and the GE was computed from FC matrices of 1 through 8 with 9 through 16th detectors [$i=1\dots 8$, $j=9\dots 16$ in Eq. (9)], as shown with dark squares in Fig. 7(d). Whole (W) corresponds to the GE computed from the full FC matrix. Table 2 shows the change of GE with respect to different stimuli for different areas of the PFC. As the GE in the left both for [Hb] and [HbO] decreases with varying stimuli difficulty in a statistically nonsignificant manner, the GE in the right both for [Hb] and [HbO] shows a statistically significant increase. This finding is consistent with the literature where the Stroop interference affect was found to be localized bilaterally.^{30,40,41,56,57} Zhang et al.³⁰ claim to see an increase in the information flow from left to right dlPFC during interference. Figure 7(d) and consecutively Table 2 lead us to speculate that connectivity strength increases in the left and right dlPFC, but the connectivity pattern is diffused in the left while being more focused for the right.

Table 2

GE values (rounded) with respect to PFC areas for the three stimuli types. Standard deviations are discarded but statistical significance (p value) was computed by paired two-tailed t-test. n is the number of subjects included in this analysis.

ST	Hb (n=6, TH=24)				HbO (n=10, TH=21)
	L+	R	IH*	W*	L+	R*	IH	W*
N	0.57	0.51	0.43	0.12	0.51	0.48	0.34	0.10
C	0.53	0.52	0.42	0.12	0.49	0.46	0.36	0.11
IC	0.52	0.55	0.47	0.14	0.46	0.56	0.36	0.13

Note: ST, stimulus type. + indicates p<0.1, while * indicates p<0.05.

It should be mentioned that the number of subjects included in the analysis in Table 2 had to be trimmed to achieve statistical significance. Only the highest and lowest GE values obtained from [HbO] data are eliminated while a higher number of subjects had to be eliminated from [Hb] data to achieve same significance levels. This result is in line with the literature, where [HbO] data are favored over [Hb] in cognitive studies due to its higher sensitivity to cognition-related hemodynamic changes.

The average interference values computed from GE for [Hb] and [HbO] data are shown in Fig. 8(a). Interference of the behavioral data (as reaction rates) shows a positive (negative) correlation for the interference values of the [HbO] ([Hb]) data, as shown in Fig. 8(b).

Fig. 8

(a) Average interference of GE (${\mathrm{GE}}_{\mathrm{IC}}-{\mathrm{GE}}_{\mathrm{N}}$) from [Hb] ($n=6$) and [HbO] ($n=10$) data for the different areas of PFC (+ for $p<0.1$, * for $p<0.05$). (b) Correlation with the interference of the behavioral data (reaction rates $R{T}_{\mathrm{Ic}}-R{T}_{\mathrm{N}}$) with interference of GE (${\mathrm{GE}}_{\mathrm{Ic}}-{\mathrm{GE}}_{\mathrm{N}}$) for the whole PFC.

4.1.

On the Accuracy of Low Pass Filtering

We have shown in the first two scenarios that filtering with a low pass filter does not improve the correlation estimation since we do not have access to the frequency characteristics of the systemic fluctuations. The shortcoming of a correlation computation after signals is low pass filtered, which can be proven by deductive reasoning as follows.

Assume that the $i$’th detector signal [$x_i(t)$] is modeled as the sum of the brain activity, $s_i^b(t)$, and nuisance term, $s^n(t)$, [$x_i(t)=s_i^b(t)+s^n(t)$] and $j$’th detector signal as [$x_j(t)=s_j^b(t)+s^n(t)$]. Assuming that there is an overlap in the spectrum of these signals, as depicted in Fig. 9, then a low pass filter applied to the sum signal (spectrum with $|X_i(f)|$) will inevitably include a piece of the nuisance spectrum.

Fig. 9

Spectral composition of the signal from the $i$’th detector, $|X_i(f)|$, can be given as the sum of the spectrum of the band-limited brain signal, $|S_i^b(f)|$, with the spectrum of the wide-band nuisance signal, $|S^n(f)|$, which is the sum of systemic fluctuations and other types of noises. Typical filter response is superimposed on the last graph.

Since low pass filters have nonideal characteristics, a leakage of unwanted signals will be present in the filtered signal. This will result in inaccurate estimation of the correlation coefficient. This same problem of a background activity dominating the channels exists in EEG analysis, where due to volume conduction effect, channels are superimposed with this background activity. This must be removed if accurate correlations are to be computed. Our results show a deviation from the original connectivity values even after signals are filtered with various filter types. Error in FC matrices for the Butterworth type of filter produces very similar yet lower error values compared to no preprocessed data, as expected. As the contamination weight increases (i.e., increase in $w_i^s$) so do the errors in computing the correlation values, as shown in Fig. 5(b).

4.2.

On the Accuracy of PC-based Approach

Since there is some sort of dependence between the brain signals and the contaminating signals both in the time and frequency domains, it would be extremely difficult to try to separate them from each other with conventional signal processing tools, such as digital filtering, PCA, or ICA. All these methods require several conditions of independence in time or frequency domain. Kalman filtering would work if there were access to systemic fluctuation data with a detector placed closer to the source (short detector), as has been proposed by several investigators.^6,13,63–66 Hence, the only solution if one wants to compute the correlation under such dependence assumption and a fixed source–detector separation is the use of a PC approach, where one can compute the correlation under a controlled variable. This is exactly what the simulation shows. Its accuracy exceeds by far both the filtered scenarios by almost two to three folds, as shown in Fig. 5(b). Yet even this approach cannot reach a zero level accuracy due to its use of a filter (this time a high pass filter used when computing the nuisance signal that will be used as the controlling variable). Hence, the success of this method depends on how good the nuisance signal can be provided to the algorithm. One solution is the use of a short detector, which will pick up signals only from the skin. Similarly, the use of a pulse oximetry signal even from the fingertip might prove to be useful as long as one can account for the delays in the systemic activity. Similar to the low pass filter’s performance, the PC-based approach starts to fail as the contamination level increases.

One alternative approach could be the use of this common nuisance signal (${\mathrm{HBO}}_R$) in regressing it from the main detector signal using a general linear model (GLM) based approach. This might be possible but it will introduce other computational steps before actually computing the correlation (first the GLM approach, then reconstruct the HBO signal with the $\beta$‘s that are statistically inferred from the GLM approach, then compute the correlation). Hence, it will undoubtedly introduce more errors even due to the use of numerical methods.

4.3.

Real Data Analysis

Having shown that if one has access to only far detector data and no means to measure some sort of systemic activity, then PC-based connectivity analysis is the best choice. Our results show that with conventional means of computing, the connectivity yielded statistically nonsignificant findings in our fNIRS connectivity analysis. A very significant difference is observed among the GE for the three different stimulus types, as shown in Table 2. Several studies have shown a right dominance for the Stroop activity.³⁰ Our results also confirm a shift toward right dlPFC connectivity as the stimulus becomes more challenging [from neutral stimulus (NS) to incongruent stimulus (IcS)], as shown in Table 2. As the GE in the left PFC decreases for increasing task complexity, the right PFC picks up and the interhemispheric GE increases, as shown in Fig. 8(a) although not significantly. There is a clear dominance of the left dlPFC, as indicated by many activation and connectivity studies.^40,67–74 Astolfi et al.⁷⁵ showed a bilateral connectivity in the Stroop task with a “predominance of outflow from right premotor and prefrontal cortical areas.” In both the [Hb] and [HbO] GE results, we see a drop in the incongruent stimuli compared to the neutral stimuli in the left side and a consequent increase in the interhemispheric and whole head interference results [Fig. 8(a)]. This can be interpreted as the aggregation of shorter paths to one longer path (leading to an increase in GE) for the [HbO] signal.

Similar to some fMRI studies, there is a negative correlation between the interference observed in reaction rates and the interference computed by GE values, as shown in Fig. 8(b), for [HbO] at the left and interhemispheric connectivity metrics. A significant positive correlation is observed for the [HbO] for whole head connectivity analysis ($r=0.729$, $p=0.0259$) while a significant negative correlation ($r=-0.745$, $p=0.0213$) is observed for interhemispheric connectivity from [Hb] data, as shown in Fig. 8(b). A positive correlation between the interference computed for the reaction rates and GE values means that as the cognitive challenge increases, sequential short paths are replaced by one long path, a direct connection between distant areas. An exact opposite is observed for the [Hb] interference, where the negative correlation means that a smaller interference in reaction times gives higher interference in GE values. An explanation can be as follows: as the cognitive challenge increases, brain regions working in coherence increase, leading to a simultaneous demand of oxygenated blood to those regions (hence an increase in the GE observed for [HbO]), yet the venous side is not necessarily as coordinated leading to a drop of efficiency. Neurovascular coupling literature has united on the finding that both neurons and astrocytes may lead to a vasodilation of the arteriole smooth muscle cells in response to glutamate releases during neuronal activation.^76,77 We might also infer that [Hb] activity is mostly regulated on a regional basis in a passive way and independent of the unification required for responding to a cognitive challenge, much similar to a balloon effect.^78–80 That type of demand can be explained only by the increase in [HbO] activation and its GE metrics.