2.1.

Participants, Protocols, Data Acquisition, and Preprocessing

2.1.2.

Data acquisition

During the protocol, a continuous wave, multichannel functional NIRS (fNIRS) brain imaging system (Cephalogics, Massachusetts) was applied to each subject’s forehead to record cortical hemodynamic alterations during resting state [see Fig. 2(a)]. The wavelengths used to calculate changes of oxygenated hemoglobin ($\mathrm{\Delta }\mathrm{HbO}$) and deoxygenated hemoglobin ($\mathrm{\Delta }\mathrm{HbR}$) were 750 and 850 nm. The sensor array consisted of 18 pairs of two light sources each and 18 detectors with the nearest interoptode distance of 2.5 cm and the second-nearest interoptode distance of 3.5 cm, forming 75 measurement channels and covering the forehead entirely [see Fig. 2(b)]. To find anatomical locations of the optodes, the international 10–20 system of electrodes placement was used to coregister the locations of sources and detectors to a standardized brain atlas. In addition, the reference points of nasion (Na), inion (In), top (Cz), left ear (A1), and right ear (A2) in the 10–20 system were measured using a 3-D digitizer to perform the brain size normalization in subsequent steps. The region of detection covered approximately $10\times 20{\mathrm{cm}}^2$ surface area in the forehead, which mainly covered the Brodmann areas 9, 10, and 11. The detected fNIRS signals were acquired at a sampling frequency of 10.8 Hz.

Fig. 2

Experimental setup and optode/probe geometry: (a) optical probe placement on a partcipant’s forehead and (b) optical optode geometry. Open circles represent source locations and solid circles represent detector locations. The source–detector (SD) separation is 2.5 cm, which indicates that the nearest SD distance is 2.5 cm and the second-nearest SD is 3.5 cm. The total number of SD pairs is 75.

2.2.

Brain Imaging with Atlas-Guided Diffuse Optical Tomography

Our previous study⁵ and many others from the literature^2,4,6 extensively introduced the atlas-DOT algorithm. We briefly describe the method as follows.

First, the predesigned ICBM152 brain template including 176,192 voxels was generated using the MRI structure scan, an average brain image taken from 152 subjects.³¹ The size of the brain template was $229\times 192\times 192{\mathrm{mm}}^3$. A careful segmentation was performed to separate the scalp, skull, gray, and white matters. A finite element surface mesh was generated for each segment and later grown inside each of the four segments (i.e., scalp, skull, gray, and white matters) to form a complete volume, constituted with tetrahedrons. A total of 179,162 vertexes were near-uniformly distributed in Montreal Neurological Institute (MNI) space. Considering the differences of optical properties in the real human brain, each segment was assigned a specific set of optical properties suggested by a previous study.² Table 1 provides the optical properties of different brain layers.

Table 1

Optical properties of head tissues for modeling.

Second, spatial coregistration and head size normalization were performed in order to minimize the confounding by the brain size variation generated as a result of normal aging. As indicated in Sec. 2.1.2, we performed location reference measurements before acquiring data for both young and older adults. Specifically, we performed an affine transformation³⁴ between the standard 10–20 system landmarks in ICBM152 MNI space and the 3-D digitizer measurements of subject-specific locations.^5,31 The coordinates of projected optodes for each individual were then normalized in the standard MNI space in the ICBM 152 brain atlas. Figure 3(a) shows the normalized optical source and detector locations plotted with the ICBM 152 cortical image template. The red spots represent the locations of light sources and the blue spots represent the locations of light detectors.

Fig. 3

Illustration of the AVC algorithm. (a) Locations of optical optodes (sources: red dots; detectors: blue dots) projected on the ICBM 152 brain cotrical template. (b) Spatial distribution of DOT measurement sensitivity after FEM forward modeling, based on the normalized optical optode geometries. (c) Averaged locations of AAL 116 regions by colored spots. They have been separated into six networks: default (yellow), frontal–parietal (light blue), occipital (blue), sensorimotor (navy), cingulo-opercular or limic (orange), and cerebellum (dark red). (d) A 3-D view of the AAL 116 ROIs and the surfaces based on Ref. 32. (e) AVC categorizes all atlas-DOT voxels interrogated by the fNIRS probes in this study into six different brain regions as represented by different colors. (f) The locations or coordinates of the 34 ROIs within the atlas-DOT measured in this study. (g) Time-sequence plots of 34 nodal $\mathrm{\Delta }\mathrm{HbO}$ values from one subject (subject 01) during a 10-min resting-state period. The $x$-axis represents time, and the $y$-axis represents the oxygenated hemoglobin changes of 34 nodes.

Third, the forward modeling of light propagation in the ICBM152 brain template was performed based on the Rytov approximation.^35,36 The light sensitivity matrix (or Jacobin matrix) was created by performing the finite element estimation of light propagation inside the ICBM152 brain atlas using a finite element method (FEM)-based MATLAB^® package called NIRFAST.³⁷ The optical probe geometry shown in Fig. 3(a) generated 75 measurement channels. A depth compensation algorithm with a tested compensation factor of 1.7 was further applied to the sensitivity matrix to project the reconstructed DOT images at more accurate brain layers.³⁸ By visually checking the light sensitivity in the forehead [Fig. 3(b)], we observed that the measurements primarily covered Brodmann 9, 10, 11, and 46, including the rostrolateral prefrontal cortex, dorsolateral prefrontal cortex, and parts of the parietal regions.

Last, to recover spatial–temporal changes of light absorptions resulting from brain resting-state hemodynamic activities, we implemented a 3-D image reconstruction approach repeatedly at different time points. We performed inverse image reconstruction using the Moore–Penrose generalized approach with Tikhonov regularization.³⁹ The reconstructed images were 4-D matrices with 176,192 spatial voxels and 6480 temporal points (10 min with a sampling rate of 10.8 Hz). Consequently, determining absorption changes at two wavelengths led to reconstructed images of relative changes in ΔHbO concentrations, based on spectral decomposition of the extinction coefficients for both wavelengths.⁴⁰

2.3.

Graph Formation with Automated Voxel Classification

To perform GTA, we first needed to generate a proper functional connectivity matrix called graph formation, in which the ROIs are classified as nodes, and interactions between nodes are termed edges.¹⁸ Edges in brain networks could be either binary or weighted. Binary edges indicate the presence or absence of interactions (“0” for absence or “1” for presence), while weighted edges represent interaction strength. Using a threshold can remove weak edges. Edges in functional brain studies are usually quantified by calculating the correlation coefficient between two nodes. In functional brain network analysis, graph formation is usually performed directly from the measured data. Niu et al.²⁹ introduced a method using the channel-wise time sequence data as the graph formation input. However, this method does not work for atlas-DOT, because time-dependent data in atlas-DOT are not given in channels but in voxels. In graph formation for fMRI measurements, voxel-wise time-sequence data are often classified or grouped into certain clusters based on either anatomical locations or predefined functional ROIs.⁴¹ The method we applied here shares a similar idea. Following a widely accepted anatomical template, namely, the AAL atlas,³² we presented or projected the 116 anatomical ROIs on the ICBM 152 brain atlas (see more details in Sec. 2.3.1). We then performed AVC to identify respective anatomical clusters on the ICBM 152-based atlas-DOT. Next, we calculated Pearson correlation coefficients by performing cross-correlations of the time series between each pair of identified clusters to create the adjacency matrices. Figure 4 shows a flowchart describing the overall selection process of AVC.

Fig. 4

Flowchart of AVC.

2.3.3.

Cross-correlation to generate adjacency matrix

Former studies using fNIRS-based GTA to investigate brain networks suggested that correlations of brain hemodynamic changes or fluctuations among different brain regions could represent brain functional connectivity.^27–29 The same strategy or rationale was applied in this study to nodal $\mathrm{\Delta }\mathrm{HbO}$ time series, using the 10-min resting-state fNIRS measurements to establish the adjacency matrix for each participant. The cross-correlation between each pair of nodes was performed for the given time series, as shown in Fig. 3(g), and Pearson correlation coefficients ($R$) were computed to form a $34\times 34$ adjacency correlation matrix [see Fig. 5(a)]. Note that the color in Fig. 5(a) denotes the values of the correlation coefficients: blue indicates a low-correlation coefficient and red represents a high-correlation coefficient. The adjacency matrix was further converted into a binarized matrix by setting a threshold. The correlation coefficient between node $i$ and $j$ $(i,j=1,2,\dots ,34)$ was set to 1, if the correlation value was larger than the given threshold and 0 otherwise [see Fig. 5(b)]. The edge was defined as functional connections ($\text{either}=1$ to represent a significant correlation between two nodes or = 0 to represent no significant correlation).

Fig. 5

(a) The adjacency matrix generated by cross-correlation of $\mathrm{\Delta }\mathrm{HbO}$ values at 34 nodes from one human subject. The $x$- and $y$-axes represent the node numbers; red color represents a value with high correlation between the temporal profiles of two nodes, while blue color represents a low-correlation coefficient. (b) A binary matrix thresholded by a middle sparsity of 0.25. Note that the diagonal elements were set to be zero since they were not the actual two-channel correlations. (c) A 3-D view of spatial representation of nodes and edges generated to show the RSFC. The connected gray lines represent the fNIRS-derived RSFC measured in this study, whereas the unconnected dots show the ROIs classified by our AVC algorithm.

In principle, different settings or selections on the cross-correlation threshold will result in a different binarized matrix. In this study, we applied different thresholds to the adjacency matrix to obtain a sequence of binary matrices. Specifically, the sparsity-based approach was utilized as suggested by Niu et al.²⁹ Sparsity ($S$) for a fixed graph was defined as the number of current existing edges in this graph divided by the maximum possible number of edges in the current graph. In this study, the range from 0.1 to 0.5 (i.e., $0.1<S<0.5$, $\text{interval}=0.01$) was chosen to be the standard threshold sequence, as a previous study reported.¹¹ Then, the threshold sequence was applied to the adjacency matrix to generate 41 binarized network matrices for each subject. Figure 5(c) shows an example from one subject’s data, revealing a spatial representation of nodes and edges from one binary matrix with a middle threshold ($S=0.25$). Colors of nodes represent different brain anatomical regions, whereas gray lines between two nodes represent respective connections. Actually, the connected lines represent the fNIRS-derived RSFC measured in this study, whereas the unconnected dots show the ROIs that our AVC algorithm classified.

2.4.

Generating Brain Network Matrices Using Graph Theory Analysis

Based on the adjacency matrices, we further quantified the resting-state brain network parameters using GTA. Graph theory metrics for functional brain networks were calculated under global and local network characteristics. In this study, we did not include all of the global and local parameters given by GTA because of the processing load required and the purpose of this study. We focused on demonstrating our AVC by using a few major parameters affected by age effects in a previous study.¹¹ Specifically, a set of global parameters investigated were global network metrics, including (i) clustering coefficient ($C_p$), (ii) shortest path length ($L_p$), (iii) global efficiency ($E_g$), (iv) normalized clustering coefficient ($\gamma$), (v) normalized characteristic path length ($\lambda$), and (vi) small-world ($\sigma$). For the local parameters, we focused on the hub information, including (i) nodal degree, (ii) nodal efficiency, and (iii) betweenness centrality. Appendix B includes detailed descriptions of these graph metrics.

2.5.

Statistical Analysis

We confirmed our classification algorithm by performing two statistical analyses. First, we performed a test–retest analysis to evaluate the reproducibility of quantification of the brain networks between two measurement visits in 15 young adults. We also performed a second analysis to compare the differences of global and local features between young and older adults. Specifically, we computed global efficiency ($E_{\mathrm{glob}}$), shortest path length ($L_p$), clustering coefficient ($C_p$), normalized clustering coefficient ($\gamma$), normalized characteristic path length ($\lambda$), and small-world coefficient ($\sigma$) for all the groups. Normality was performed on all measurements, and data without normality were converted to $z$-values ($Z$) transformation before further analysis.

3.1.

Resting-State Functional Connectivity Matrices

To check the reproducibility of our approach, 10-min rest-state fNIRS measurements were taken from young adults ($n=15$) twice in two separate experiments. A total number of 75 channels allowed us to cover and interrogate most of the prefrontal regions [see Figs. 3(a) and 3(b)]. Our newly developed AVC algorithm was used resulting in 34 nodes within the cortical regions reconstructed using atlas-DOT [see Figs. 3(f) and 5(c)]. All the 10-min temporal profiles from the 34 nodes were extracted [see Fig. 3(g)] and used to create RSFC adjacency matrices, whose elements or pixels were filled with respective Pearson correlation coefficients ($R_{\mathrm{RSFC}}$) between each pair of the 34 nodes. Figures 6(a) and 6(b) show such RSFC matrices group-averaged over the 15 young adults for visits 1 and 2, respectively. Visual inspection showed that the two RSFC matrices were consistent with each other in overall distribution of correlation patterns and correlation strengths. More quantitatively, we compared these two matrices by plotting the two separate $R_{\mathrm{RSFC}}$ sets against each other, as shown in Fig. 6(c). The comparison demonstrated a linear correlation between the two visits $(R=0.58,p<0.0001)$. These results demonstrated that our AVC algorithm was strongly reproducible when used with atlas-DOT to obtain RSFC matrices.

Fig. 6

RSFC matrices averaged over young adults at (a) visit 1 and (b) visit 2. Each data pixel marks the Pearson correlation coefficient ($R_{\mathrm{RSFC}}$) calculated between each pair of 34 nodes. Note that the diagonal elements were set to be zero since they were not actual two-channel correlations. (c) A linear relationship is illustrated between two $R_{\mathrm{RSFC}}$ data sets of young adults at visit 1 ($x$-axis) and visit 2 ($y$-axis). The red line is the fitted line of the two groups $(R=0.58,p<0.001)$, indicating a significant correlation between the two visits.

It is known that GTA allows for the analysis and categorization of topological measures of brain networks into two groups: (a) global network metrics: small-world properties (e.g., clustering coefficient $C_p$, characteristic path length $L_p$, normalized clustering coefficient $c$, normalized characteristic path length $l$, and small-world $s$), efficiency parameters (local efficiency $E_l$ and global efficiency $E_g$), modularity $Q$, hierarchy $b$, and assortativity $r$ and (b) nodal characteristics: nodal degree, nodal efficiency, and nodal betweenness.²⁷ We followed the same grouping strategy to present our analysis in Sec. 3.2.

3.2.

Global Network Characteristics and Small-World Features

To further evaluate the reproducibility of our algorithm, we quantified and compared the global network and small-world features of the brain networks between the two visits of young adults and between young and older adults. Figures 7(a) to 7(c) show the respective values of clustering coefficient ($C_p$), shortest path length ($L_p$), and global efficiency ($E_g$) of young adults obtained from two separate measurements along with the corresponding readings from the older adults. Two-sample $t$-tests for each of the global brain network parameters were performed at each sparsity value from 0.1 to 0.5, as a previous study suggested.¹¹

Fig. 7

Global characteristics of the GTA. (a)–(c) indicate the clustering coefficient ($C_p$), shortest path length ($L_p$), and global efficiency ($E_g$) quantified from young adults in visit 1 (solid circles) and visit 2 (open circles) and from older adults (pluses). The dashed line on the bottom of panel (a) and (c) marks the sparsity ranges where significant differences in respective network parameters exist between two visits (or measurements) of young adults. The solid lines on the bottom of each panel represent the sparsity ranges where significant differences in respective network parameters exist between young and older adults. (d)–(f) indicate the small-world characteristics, including the normalized characteristic path lengths ($\lambda$), normalized clustering coefficient ($\gamma$), and small-worldness ($\sigma$) quantified from young adults in visit 1 (solid circles) and visit 2 (open circles) and from older adults (pluses), respectively. The solid lines on the bottom of (d) to (f) have the same meaning as those in (a) to (c).

3.3.

Reproducibility in Local Graph Parameter

In this study, network hubs were also quantified by nodal degree ($N_i$), nodal efficiency ($E_{\mathrm{nod}}$), and betweenness centrality ($N_{\mathrm{bc}}$). The nodal metrics were constructed at a sparsity threshold of 0.16, as a previous study suggested,¹¹ to ensure that the networks of both young and older adult groups had the same number of nodes and edges. The hubs were then selected for each of the three nodal parameters (i.e., $N_i$, $E_{\mathrm{nod}}$, and $N_{\mathrm{bc}}$) with respective values larger than one standard deviation of the corresponding average values over all nodes, as previous studies suggested.^27,29 Figures 8(a) to 8(c) demonstrate axial or top views of the hubs determined from the young adults in visits 1 and 2 and from the older adults for all three nodal metrics of degree ($N_i$), efficiency ($E_{\mathrm{nod}}$), and betweenness centrality ($N_{\mathrm{bc}}$). As Niu et al.²⁹ indicated, the betweenness centrality was considered as the major reference for the hub measurements.9

Fig. 8

Hubs determined from (a) young adults in visit 1, (b) young adults in visit 2, and (c) older adults. The hubs are presented based on nodal degree ($N_i$), nodal efficiency ($E_{\mathrm{nod}}$), and betweenness centrality ($N_{\mathrm{bc}}$). Yellow dots represent the hubs within the default mode network, light blue dots represent the hubs within the frontal–parietal network, and dark blue dots represent those within the sensorimotor network. [See Figs. 3(f) and 5(c).]

4.1.

Development of Automated Voxel Classification With Atlas-Guided Diffuse Optical Tomography

Atlas-DOT is an imaging method recently established within the NIRS-based neuroimaging field.^1–8,46 GTA is a network analysis tool that has been commonly used in functional neuroimaging to study brain connectivity and networks.^{13–15,18,19} To date, however, atlas-DOT and GTA have not been used jointly to extract or derive large-scale brain connectivity and networks from the NIRS measurements because of computational challenges posed by atlas-DOT. Specifically, graph formation in GTA is a difficult and computationally expensive process for atlas-DOT because of the large amount of voxels in atlas-DOT images.

The major novelty in our current work was the development of AVC with atlas-DOT, which affords image processing speed and localization accuracy in optical brain imaging. Such an automated classification is critical because it helps to extract and identify cerebral hemodynamic signals from appropriately matched anatomical locations, thus improving the accuracy of interpretations of local and global brain activities. Commonly used methods for registering NIRS data on a human brain template are often based on the EEG 10–20 system,^47–49 which is not directly applicable to the 3-D atlas-DOT images. In addition, routinely used voxel classification in fMRI is based on an affine transformation, which converts the voxel space to MNI space, followed by a voxel matching with a standard MNI structure template.⁴¹ However, this approach works only when the voxels in the brain images are homogeneously distributed or have the same equal voxel size. This approach is not applicable to atlas-DOT for several reasons. First, because atlas-DOT is FEM based, different input parameters (such as the head size, geometry, and shape) to the atlas model would result in a different total number and distribution of voxels in the data space because of the nature of FEM computation and methodology. Second, the spatial distribution of optical densities used for atlas-DOT is variable and heterogeneous, and it depends on the optical properties of the human head with multiple layers. Such heterogeneous multilayer structures result in an unpredictable number and pattern of voxels in a certain region. Therefore, the spatial classification method used by fMRI cannot be transferred directly for voxel classification in atlas-DOT. Furthermore, no universal template can be generated to match the FEM voxels versus specific anatomical structures because each FEM voxel model must be determined according to the optical probe geometries or setups. Thus, we were motivated to address this need and develop the current technique, which is a comprehensive yet efficient solution for the atlas-DOT voxel classification.

The principle of AVC was based on the “point-in-polyhedron problem”⁴⁴ approach. The basic idea was to identify whether a specific voxel was (i) inside, (ii) outside, or (iii) just on a particular surface. Specifically, our approach focused on the location of each atlas-DOT voxel with respect to the surface of several predefined brain structures. Thus, no matter what the voxel distribution was, once the voxel coordinates in the MNI space were assigned, the classification of DOT voxels to an appropriate anatomical ROI would be automatically achieved. Because this morphological operation or classification with atlas-DOT is efficient, it can be applied to other atlas-DOT-based brain imaging applications.

4.2.

Test–Retest Reliability of Resting-State Functional Connectivity Using Automated Voxel Classification With Atlas-Guided Diffuse Optical Tomography

Another innovative aspect of this study was that it demonstrated the good reproducibility of AVC with atlas-DOT in assessing RSFC. We found strong consistency between two matrices of RSFC taken from two separate visits in two weeks (Fig. 6). A recent report by Niu et al.²⁷ revealed a better reproducibility of the RSFC correlation, with a higher test–retest correlation coefficient $(r=0.93,p<0.0001)$ than ours $(r=0.58,p<0.0001)$. One possible reason for this difference resulted from the fact that our two test–retest measurements were separated by two weeks, whereas the measurements in Ref. 27 were obtained from only the same set of measurements taken at two time periods.

Furthermore, the global network matrices obtained from the two visits have shown an overall high consistency [see Figs. 7 and 8], as two-sample $t$-tests also proved. Specifically, except for the clustering coefficient in the sparsity range from $0<S<0.24$, global network features demonstrated good reproducibility along the sparsity range of 0.1 to 0.5. This finding confirmed that during a resting state, the prefrontal cortex activities were noticeably repeatable.²⁷ Similar hub locations observed in both visits also validated the reproducibility of local network matrices. For example, we observed bilateral network hubs in middle frontal and superior frontal gyrus within the three nodal metrics of degree, efficiency, and betweenness centrality in both visits. All these brain regions were consistent with those reported and observed previously by fMRI and PET brain network studies.^11,18,21

4.3.

Comparison of Age Effect on Resting-State Functional Connectivity Using Automated Voxel Classification With Atlas-Guided Diffuse Optical Tomography

Another innovative aspect of this study was that it clearly revealed the age effect on the resting-state brain network observed by combining AVC with atlas-DOT. Although both fMRI and PET have observed and demonstrated age effects on the brain networks by Lieu et al.¹¹ and Meunier et al.,¹² no report on this topic exists using optical methods. In this study, we further investigated partial default mode networks and frontal–parietal networks of young adults (ages 25 to 43) and older adults (ages 65 to 92). Our findings were compatible with previous studies,^11,50,51 which showed a higher global efficiency in young adults than in older adults. It is believed that the anatomical alterations of aging are responsible for these changes, especially after the age of 65.⁵² Our present results also revealed a difference in clustering coefficient between the older and young adults. Similar findings from previous reports^11,25 suggest that a possible reorganization of frontal network occurs as people age. It is also suggested that information processes are less economical in older adults, especially in the frontal and temporal cortical and subcortical regions.²⁵ Several other reports have also indicated an increase of shortest path length with an increase in age resulting from the loss of long-distance connections and interconnected hubs.^25,53 However, we found no differences in this network metric, possibly because our optical probe setup covered a limited brain region (i.e., prefrontal), leading to a shortage of long-distance brain network measurements.

Finally, we have observed a significant decrease in small-worldness in older adults, implying a degeneration process, in good agreement with previous studies.^11,12,18 Our findings have also identified several hubs of the brain networks located within the frontal cortical regions in both age groups, as reported by several previous publications.^11,18 Meanwhile, a decrease in the number of hubs in the older adults was observed as compared to the young adults. Specifically, the nodal betweenness in default mode regions (such as middle frontal gyrus and superior frontal gyrus) was diminished in the older adults.

4.4.

Consideration of Spatial Leakage

In recent years, much research in the context of EEG- and MEG-based RSFC has focused on possible confounding effects resulting from spatial leakage.^54–56 This concern is essential when the measurement data are reconstructed into source-space. Since such connectivity images or maps are formed by solving an ill-posed inverse problem, there are inherent correlations among adjacent vertices. Each reconstructed brain vertex can be expressed as a linear combination of the channel-space data. This confounding issue is more severe for correction in MEG- and EEG-based connectivity studies than fNIRS-derived brain network mappings because the optical sensitivity falls off much faster (within 1 to 2 cm) than MEG ($\sim 5$ to 6 cm).⁵⁵ It is still worthwhile to estimate the spatial leakage that may cause misleading conclusions in our study. For this purpose, we have added Appendix C that describes our estimation approach to quantify the spatial leakage. Appendix C presents the results and draws a basic conclusion. This extra section mainly implies that the spatial leakage is short-distanced, and the nearby region within 1 to 2 cm could be affected. Future investigation is needed in order to understand better how this distance-related artifactual connectivity can largely affect the brain network interpretations. One possible solution for correction of spatial leakage is to obtain a measurement of spatial leakage artifact by an “empty-room” recording and then to normalize the actual human brain data with respect to the artifactual data.⁵⁵ Another approach, named geometric correction scheme, was based on subtraction of the spatial leakage geometry model from the reconstructed image data and approved efficient to minimize the spatial leakage effect.⁵⁴

4.5.

Limitations of this Study and Suggestions for Future Works

Because the light penetration depth through the human brain tissue (not including scalp nor skull) is not more than 1 cm, the atlas-DOT technique can image only the brain regions mainly in gray matter.⁵⁷ Thus, our atlas-DOT-based graph formation can provide brain network nodes and other metrics within only 1 cm below the cortical surface. Imaging a deeper brain layer requires larger separations between optical sources and detectors, which may lead to a significant signal-to-noise ratio decrease.⁵⁸ Suggestions for future studies include simultaneous or parallel measurements of EEG recordings and fNIRS, which may help us better understand the relationship between neuroelectrical and cerebrovascular functions and their networks.

Another limitation of this study was associated with physiological noises mixed within the measured fNIRS signals. Such noises were caused by superficial layers of the human head, such as the scalp and skull.^59,60 To reduce these noises, we applied a global referencing method, which seemed promising because our results showed high consistency and reproducibility. However, we suggest further examining and confirming the noise effect using other approaches. Future studies could use a computational method such as independent component analysis^27,61 or an experimental approach such as short-distance measurement regression^{2,4,7,8,46,57} to minimize global signal or noise effects.

It is also highly suggested to increase the number of optodes in future studies in order to cover the entire cortical region. Because this study focused primarily on the development of AVC, we chose to cover only the frontal regions of the subjects’ heads to achieve efficient data acquisition. The cortical region detected in this study was only partial: about one-quarter of the entire brain cortical area. Node extraction in this study contained 34 nodal regions (out of 116 regions) in the AAL 116 atlas. Thus, further studies intending to rigorously determine resting-state functional brain networks using atlas-DOT should increase the number of recording optodes and cover the entire cortical region. In this way, researchers can optimally utilize AVC with atlas-DOT and obtain the most accurate results of resting-state brain network analysis.