Automated voxel classification used with atlas-guided diffuse optical tomography for assessment of functional brain networks in young and older adults

Abstract. Atlas-guided diffuse optical tomography (atlas-DOT) is a computational means to image changes in cortical hemodynamic signals during human brain activities. Graph theory analysis (GTA) is a network analysis tool commonly used in functional neuroimaging to study brain networks. Atlas-DOT has not been analyzed with GTA to derive large-scale brain connectivity/networks based on near-infrared spectroscopy (NIRS) measurements. We introduced an automated voxel classification (AVC) method that facilitated the use of GTA with atlas-DOT images by grouping unequal-sized finite element voxels into anatomically meaningful regions of interest within the human brain. The overall approach included volume segmentation, AVC, and cross-correlation. To demonstrate the usefulness of AVC, we applied reproducibility analysis to resting-state functional connectivity measurements conducted from 15 young adults in a two-week period. We also quantified and compared changes in several brain network metrics between young and older adults, which were in agreement with those reported by a previous positron emission tomography study. Overall, this study demonstrated that AVC is a useful means for facilitating integration or combination of atlas-DOT with GTA and thus for quantifying NIRS-based, voxel-wise resting-state functional brain networks.


Introduction
Atlas-guided diffuse optical tomography (atlas-DOT) has undergone rapid development in recent years. Custo et al. 1 initially introduced the concept of atlas-DOT, after which it was further developed with high-density optodes and successfully validated using computational simulations and functional magnetic resonance imaging (fMRI) performed on human subjects. [2][3][4] Atlas-DOT utilized a finite element technique with either subject-specific models 3,4 or a standard magnetic resonance imaging (MRI) brain template (such as ICBM 256) in the forward calculation. 5,6 A significant improvement of accuracy in image reconstruction and source localization has been achieved with atlas-DOT. 7 More recently, atlas-DOT has been utilized together with the general linear model analysis, which further extended its feasibility of measuring hemodynamic changes in complex brain tasks. Several groups have reported consistent and accurate results of brain hemodynamic changes under visual stimulation, 7 speech, 8 and risk decision-making. 5 It is well known that the human brain is naturally organized into groups of networks, with anatomical brain regions involved in either individualized processing or integration with other brain regions, to accomplish different functions. The brain networks experience changes in their properties under a variety of conditions, such as taking on different physical or mental tasks, facing strong emotions, 9 or undergoing normal aging. [10][11][12] Thus, noninvasive mapping of the human brain's structural and functional connectivity enables researchers to better understand the architecture of the brain and to clearly reveal connectivity changes in entire brain networks. Increasingly, many neuroimaging techniques, such as fMRI, diffusion MRI, electroencephalographic (EEG), and DOT, have been recently used to investigate large-scale brain networks. [13][14][15][16][17] Accordingly, researchers have developed corresponding mathematical frameworks that can better model and analyze brain networks for a variety of imaging modalities. Graph theory analysis (GTA) is one of the analytical methods developed to examine largescale complex brain networks. [13][14][15]18,19 It can provide an uncomplicated and yet powerful mathematical means to characterize the brain networks' topological properties.
Specifically, GTA is a mathematical method for the analysis of complex network systems; it finds applications in the modeling of interactions in information, social, and biological systems. Recently, GTA has been reported in the neuroimaging literature to analyze resting-state functional connectivity (RSFC). 18,20,21 It depicts a brain network as a graph with different anatomical brain regions portrayed as nodes and interactions between brain regions represented as links. Based on GTA, a number of studies have found alterations in network parameters in aged patients. Liu et al. 11 conducted a resting-state positron emission tomography (PET) study in young and aged patients and reported a decrease in network efficiency for older patients, who also maintained a higher clustering coefficient than younger subjects. This result indicates that the brain network becomes more random in nature as humans grow older. The reduction in gray matter as a result of aging has also been associated with a decrease in global efficiency and increase in path length for older subjects. 22,23 The activities of the default mode network have been also reported to decrease in older adults, 24 indicating changes in the network's properties along with aging. Furthermore, the aging process has been implicated in alterations in modularity and number of hubs of the brain network. 25,26 In the field of near-infrared spectroscopy (NIRS), a few recent studies have combined GTA with channel-wise NIRS and revealed topological organization and architecture of largescale, resting-state human brain cortical networks. [27][28][29] However, no report combines atlas-DOT with GTA in brain network research. In addition, little validation exists for the atlas-DOT technique in the assessment of RSFC. The difficulties in using atlas-DOT to form brain networks result from the difficulties of graph formation, the initial necessary step in GTA. Instead of using a channel-wise approach, a data set of resting-state atlas-DOT images is four-dimensional (4-D), consisting of multiple three-dimensional (3-D) cortical volumes in a time sequence. Because a reconstructed atlas-DOT image often consists of a large number of voxels, it is computationally expensive and lacks an appropriate means to generate the connectivity matrices needed in GTA. Routinely used techniques to extract regions of interest (ROIs) in fMRI analysis (such as threshold-, anatomy-, and function-based ROI 30 ) cannot be directly applied to atlas-DOT.
In this paper, we introduced an automated voxel classification (AVC) approach to facilitate graph formation for atlas-DOT images by identifying or grouping unregulated voxel distribution. By using a subject-averaged brain template (ICBM 152), 31 we demonstrated that our method could fit for any brain size. We then provided a solution of voxel classification based on the automated anatomical labeling (AAL) of activations in 116 segmented structures (AAL 116). 32 In this way, the graph formation was guided by the anatomical structure of the brain. To test our approach, we conducted a test-retest assessment by measuring RSFC in 15 young adults. Moreover, we quantified age-related changes in the chosen brain network and compared our results with previously published PET-derived reports. Both reproducibility and cross-modality comparison indicated that our AVC with atlas-DOT is an efficient and feasible tool for assessing functional brain networks.

Methods and Materials
We developed a three-section data processing strategy to perform atlas-DOT with AVC. As represented in Fig. 1, the red, yellow, and blue dashed lines enclose the three major steps: data acquisition and preprocessing, atlas-DOT, and graph formation with AVC, respectively. After generating the brain network matrices (or adjacency matrices), we performed the group-level GTAs and statistical analyses to test the reproducibility of brain connectivity in the young and older adults. All of our analyses were processed using MATLAB ® 2015b (The MathWorks, Inc., Natick, Massachusetts).

Participants and protocols
The study participants comprised 15 young adults (ages 25 to 43) recruited from the University of Texas at Arlington and 21 older adults (ages 65 to 92) recruited from the city of Arlington, Texas. All subjects were right-handed with normal visual ability. No subjects reported any known diseases such as musculoskeletal, neurological, visual, or cardiorespiratory dysfunctions. Written Fig. 1 The work flow for AVC used with atlas-DOT. The first step (marked by the red dashed line box) includes the ICBM 152 brain template generation, data acquisition, and data preprocessing. The second step (marked by the yellow dashed line box), also referred to as atlas-DOT, includes finite element computation for brain atlas, forward modeling, and image reconstruction. The third step (marked by the blue dashed line box) illustrates graph formation with AVC, which includes volume segmentation using predefined AAL with 116 brain regions (AAL 116), voxel classification and ROIs generation, and formation of an adjacency matrix by cross-correlation. To validate our approach, a routine GTA with statistic analysis was followed up to access the local and global network features. Refer to details in the text. consent forms were obtained for all participants before the protocol started. The Institutional Review Board of the University of Texas at Arlington approved this study.
Before the measurement, all participants were instructed to sit comfortably in a quiet room with eyes closed and with minimal body movements. Our protocol included two repeated measurements on the young adults with a time separation of two weeks. We also included one-time recordings for 21 older adults; however, four older adults' data were excluded from data analysis because of head movements, requested recesses, or failure to complete the protocols. All the measurements of 10-min RSFC were performed during daytime in the same experiment room, maintained at a constant temperature, humidity, and light intensity.

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 (ΔHbO) and deoxygenated hemoglobin (Δ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 × 20 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.

Data preprocessing
For data preprocessing, bandpass filtering and a global autocorrelation process were sequentially applied to the channel-wise raw fNIRS data. Specifically, the bandpass filter chosen in this study was 0.02 to 0.3 Hz, as suggested by a previous study, to eliminate the physiological noise generated by heartbeat and respirations. 5 A global signal referencing approach was applied to the channel-wise fNIRS signals to further eliminate the unmodeled global physiology or structure noise. 33 Specifically, for each participant, averaging over 75-channel time sequences generated a global signal. Then, we applied an autocorrelation approach between the globally averaged signal and signals from each individual channel to perform global regression. Any channel-wise signal strongly correlated with the global signal was considered to contain strong artifacts. Global referencing was performed in these channels by deducting the global signal from the individual channel signals. 33

Brain Imaging with Atlas-Guided Diffuse Optical Tomography
Our previous study 5 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. 31 The size of the brain template was 229 × 192 × 192 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. 2 Table 1 provides the optical properties of different brain layers.
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 34 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. 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. 37 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. 38 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. 39 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. 40

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. 18 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. 29 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 timesequence data are often classified or grouped into certain clusters based on either anatomical locations or predefined functional ROIs. 41 The method we applied here shares a similar idea. Following a widely accepted anatomical template, namely, the AAL atlas, 32 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.

Anatomical automated labeling with volume segmentation
Because the ICBM 152 brain atlas provided only an average anatomical image but lacked detailed anatomical features, 32 we had to use the AAL 116-based ROIs to facilitate or define node definitions in our study (as shown in Fig. 4). While the AAL 116 ROIs were originally derived from a single-subject T1-weighted MRI, those ROIs could be still projected or normalized on the ICBM 152 brain template. 32 Relatively, there was little difference in 3-D views of corresponding AAL ROIs on the ICBM 152 brain template and the single-subject T1-weighted brain template. More explanation is given in Appendix A. Specifically, we separated the human brain into 116 anatomical ROIs based on the landmarks of cortical and subcortical structures, as marked by the AAL 116 ROIs. Further, these 116 regions were categorized into six brain networks based on Ref. 22: default mode (yellow), sensorimotor (navy blue), frontal-parietal lobe (light blue), occipital lobe (blue), cingulo-opercular/limbic lobe (orange), and cerebellum (dark red), as identified and presented in Fig. 3(c) using the respective colors. BrainNet Viewer software 42 was used to view or image the three planes.

Voxel classification to identify the region of interest
To perform graph formation within voxel-wise atlas-DOT images, we developed an AVC algorithm to assign each DOT cortical voxel to an appropriate ROI, namely one of the AAL 116 neuroanatomical regions. The selected ROIs then served as graph nodes for graph formation, followed by cross-correlations of different time sequences between each pair of ROIbased nodes. First, both the ICBM152 and the AAL116 templates were transformed to the same coordinate space. Here, we chose the MNI coordinates as our standard space coordinate reference.  As indicated in Sec. 2.2, an affine transformation was performed on the ICBM 152 atlas after the FEM forward calculation, resulting in standard MNI coordinates for each voxel p i (x i , y i , z i , i ¼ 1; 2; : : : ; n, where n is the total number of voxels). The MNI coordinates of the original AAL 116 template were obtained from a previous study. 32 Second, the surfaces of AAL 116 ðT j ; j ¼ 1; 2; : : : ; 116Þ template were extracted using the software ITK, 43 shown in Fig. 3(d), and exported them as vertices (e.g., X i , Y i , Z i , i ¼ 1; 2; : : : ; n, where n is the total number of voxels) and faces (triangles). To efficiently and correctly identify or classify whether the brain voxels were within a certain AAL region, we followed the algorithm called "point-in-polyhedron problem." 44 Specifically, the first operation, the orientation operation, tested whether a point p i falls to either positive or negative sides of a triangle within T j surface. The second operation, point classification operation, classified whether a point p i is on one of the triangle's edges or in its interior. After the voxel classification, we were able to group or classify all the voxels in an atlas-DOT brain image into 116 anatomical regions. For example, Fig. 3(e) shows all atlas-DOT voxels on the cortical template classified into six different brain regions, represented by different colors. Second, the voxels were further selected based on the forward model of light diffusion theory. Because the optical measurements in our study did not cover the whole brain volume, the sensitivity of light interrogation in brain regions deeper than 2 to 3 cm rapidly decayed due to strong light scattering with increasing depth. Therefore, we applied a threshold to eliminate all the voxels whose optical sensitivity (quantified as optical density, OD) was less than two standard deviations of the total OD within an anatomical region. After this selection, the remaining voxels were combined with the results of voxel classification from the previous step.
Using this two-step operation, we were able to quantify or classify 34 ROIs within the 116 anatomical regions. For graph formation, these 34 ROIs were regarded as nodes. Table 2 lists network names, anatomical regions or descriptions, labels, and brain atlas coordinates of these 34 nodes. In addition, Fig. 3(f) plots the averaged coordinates within the 34 ROIs. These 34 nodes in our study were included partially within the default mode network (yellow), frontal-parietal network (light blue), and sensorimotor network (navy blue). Each time sequence, averaged over all the voxels within an ROI (i.e., each of the 34 nodes), was treated as the input of node information to create the adjacency matrix. As an example, Fig. 3(g) shows a time-sequence plot of 34 nodal ΔHbO values from one subject (subject 01). Specifically, the x-axis represents time within 10 min of a resting-state period and the y-axis represents oxygenated hemoglobin changes with the node index from 1 to 34.

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][28][29] The same strategy or rationale was applied in this study to nodal Δ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 × 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; : : : ; 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 (either ¼ 1 to represent a significant correlation between two nodes or = 0 to represent no significant correlation).
In principle, different settings or selections on the crosscorrelation 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. 29 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, interval ¼ 0.01) was chosen to be the standard threshold sequence, as a previous study reported. 11 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 fNIRSderived RSFC measured in this study, whereas the unconnected dots show the ROIs that our AVC algorithm classified.

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. 11 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 (γ), (v) normalized characteristic path length (λ), and (vi) small-world (σ). 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.

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 glob ), shortest path length (L p ), clustering coefficient (C p ), normalized clustering coefficient (γ), normalized characteristic path length (λ), and small-world coefficient (σ) for all the groups. Normality was performed on all measurements, and data without normality were converted to z-values (Z) transformation before further analysis.

Evaluating the reproducibility of current approach
The evaluation of reproducibility was performed in both brain connectivity and graph theory metrics in 15 young adults. First, group-level functional connectivity adjacency matrices from the two visits were compared. Specifically, Pearson correlation coefficients were computed for all pixels in the RSFC matrices. Then a linear fitting between two visits were performed for all the pixels. Second, global (i.e., E glob , C p , L p , γ, λ, and σ) and local [i.e., nodal efficiency (E nod ), nodal degree (N i ), and betweenness centrality (N bc )] parameters of the two visits were compared. Since the fNIRS measurements were sequentially performed, a two-sample t-test was taken at each sparsity threshold for all global parameters, resulting in a range of significant tests along sparsity. A criterion of p < 0.05 was selected to define statistical significance in this study.

Analysis of the age effect on the global and local parameters
We further evaluated our method by performing an age-related analysis between 15 young adults and 21 older adults in both global (E glob , C p , L p , γ, λ, and σ) and local (nodal efficiency, nodal degree, and betweenness centrality) brain network characteristics. The two-sample t-tests for the global metrics ðE glob ; C p ; L p Þ were performed between two groups along all sparsity values (0.1 < S < 0.5). In addition, regions of hubs were compared in nodal efficiency (E nod ), nodal degree (N i ), and betweenness centrality (N bc ).

Results
We performed careful analysis to test the reproducibility and feasibility of our AVC method. The reproducibility analysis of the method was based on the measurements from young adults, with each of the individuals having participated in the fNIRS recordings twice within a period of two weeks. The method's feasibility analysis tested the age effect on the brain networks by comparing them between the young and older adults.

Resting-State Functional Connectivity Matrices
To check the reproducibility of our approach, 10-min rest-state fNIRS measurements were taken from young adults (n 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 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 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. 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. 27 We followed the same grouping strategy to present our analysis in Sec. 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. 11

Significant differences in global features between two visits of young adults
Figures 7(a) to 7(c) present the values of clustering coefficient (C p ), shortest path length (L p ), and global efficiency (E g ) determined from two visits of young adults. Overall, results illustrated a high reproducibility of the global characteristics of the brain networks. As Figs. 7(a) and 7(c) show, significant differences in C p and E g between visits 1 and 2 existed only within a short range of the sparsity. L p showed no significant difference within the given sparsity range between visits 1 and 2.

Significant differences in global features between young and older adults
For better comparison, we also plotted global features of the brain networks from older adults together with those from the young adults. Figures 7(a) and 7(c) clearly show age-related changes in C p and E g . Quantitatively, the corresponding global features of young adults from the first visit were used to carry out statistical assessment. Specifically, older adults revealed a significantly larger clustering coefficient (C p ) within the sparsity range between 0.1 and 0.47 [see Fig. 7(a)]. On the other hand, young adults had a significantly larger global efficiency (E g ) than older adults within the sparsity range of 0.1 to 0.31 [see Fig. 7(c)]. Meanwhile, we observed no significant difference in shortest path length (L p ) between the two age groups [see Fig. 7(b)]. Our findings were consistent with a previous study 11 that investigated resting-state brain networks in 115 young adults (average age ¼ 35) and 110 older adults (average age ¼ 54) using PET. Their results indicated a global efficiency decline with an increase in clustering coefficient in older adults. 11

Small-world functional network
Functional networks of the human brain have small-world characteristics. 21,45 Compared to random networks, the smallworld characteristics in the human brain include large local connectivity and an approximately identical shortest path length between any two nodes in the network. Figures 7(d) to 7(f) show the normalized characteristic path length (λ), normalized clustering coefficient (γ), and small-worldness (σ) taken from young adults in visit 1, visit 2, and from older adults, respectively. We observed that λ in all three groups approached 1 with sparsity S ¼ 0.1 to 0.5 [see Fig. 7(d)], and so did γ.
Noticeably, no significant difference existed in λ, γ, and σ between two visits of young adults, revealing excellent reproducibility in the small-world features between the two measurements. Regarding age effects on the small-world characteristics of the brain networks, our results, furthermore, demonstrated that the young adults have significantly larger λ than older adults over a range of sparsity (0.12 < S < 0.24).
In addition, the young adults have clearly exhibited larger values of γ (0.1 < S < 0.41) and σ (0.1 < S < 0.3) over a large range of sparsity than older adults.

Reproducibility in Local Graph Parameter
In this study, network hubs were also quantified by nodal degree (N i ), nodal efficiency (E nod ), and betweenness centrality (N bc ). The nodal metrics were constructed at a sparsity threshold of 0.16, as a previous study suggested, 11 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 nod , and N 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 nod ), and betweenness centrality (N bc ). As Niu et al. 29 indicated, the betweenness centrality was considered as the major reference for the hub measurements.  The current study used the betweenness centrality as the marker of the hubs.

Hubs found differently between young and older adults
In the older adult group, as shown in Fig. 8(c), two, one, and three hubs were identified under nodal degree, nodal efficiency, and betweenness centrality, respectively. In contrast, five, four, and six hubs were identified from young adults in visit 1; similarly, five, four, and six hubs were also found in visit 2.
In particular, the two visits identified very similar hubs (four in nodal degree, four in nodal efficiency, and five in betweenness centrality). These results revealed an overall decline of hub numbers in the older adult group. Furthermore, changes in regional nodal characteristics between the two age groups were also noted. Specifically, hubs in default mode networks were observed only on the right side of the middle frontal gyrus (MFG.R) and the left dorsolateral region of the medial superior frontal gyrus (SFGmed. L) in the older adult group, whereas both sides of the middle frontal gyrus (MFG.L and MFG.R) and the dorsolateral region (SFGmed.L and SFGmed.R) were observed in the young adult group. In addition, the bilateral hub pairs (MFG.L, MFG.R, SFGmed.L, and SFGmed.R) were observed in young adults, while unilateral hub (SFGmed.L) was observed in older adults for betweenness centrality. Table 3 summarizes and presents all the hub locations.

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][2][3][4][5][6][7][8]46 GTA is a network analysis tool that has been commonly used in functional neuroimaging to study brain connectivity and networks. [13][14][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. 41 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" 44 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-DOTbased brain imaging applications.

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. 27 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. 27 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

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. 11 and Meunier et al., 12 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. 52 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. 25 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 smallworldness 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.

Consideration of Spatial Leakage
In recent years, much research in the context of EEG-and MEGbased RSFC has focused on possible confounding effects resulting from spatial leakage. [54][55][56] This concern is essential when the measurement data are reconstructed into sourcespace. 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 (∼5 to 6 cm). 55 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. 55 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. 54

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. 57 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. 58 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.

Conclusion
In this study, we introduced an AVC approach that can facilitate graph formation for atlas-DOT images by grouping unregulated or unequally sized FEM voxels into anatomically meaningful ROIs within the human brain. To have an accurate human brain template for use with AVC, we chose a subject-averaged brain template (i.e., ICBM 152) 31 along with the AAL 116 template 32 to mark the brain ROIs. Thus, the graph formation was guided by the anatomical structure of the brain. To demonstrate the usefulness of AVC, we conducted a test-retest assessment by measuring RSFC in 15 young adults and quantified age-related changes in the measured brain network. Our results derived from AVC with atlas-DOT were in good agreement with those reported by a previous PET study. Overall, this study confirmed that AVC with atlas-DOT has great potential to serve as an efficient and feasible tool to quantify voxel-wise, resting-state functional brain networks.
This is the average inverse shortest path length dði; jÞ between two node i and j. n is the number of nodes in the network N. This is the minimum number of edges that link any two nodes of the network. dði; jÞ is the shortest path length between node i and node j.

Normalized clustering coefficient (γ)
E Q -T A R G E T ; t e m p : i n t r a l i n k -; x 2 ; 6 3 ; 4 1 5 γ ¼ C real p ∕C rand p : This is the mean of all clustering coefficients over all nodes in a network. C real p and C rand p are the clustering coefficient in a real network and random networks.

Characteristic path length (λ)
E Q -T A R G E T ; t e m p : i n t r a l i n k -; x 2 ; 6 3 ; 3 2 6 λ ¼ L real p ∕L rand p : This is the average of the shortest path lengths between any nodes of the network. L real p and L rand p are the characteristic path length in a real network and random networks. Not that current study averaged 100 random networks in the calculation.

Small worldness (σ)
E Q -T A R G E T ; t e m p : i n t r a l i n k -; x 2 ; 6 3 ; 2 2 3 σ ¼ γ∕λ: The normalized characteristic clustering coefficient is divided by the normalized characteristic path length. A network is considered small-world if σ > 1. This is the number of edges linked directly to a particular node. a ij is the i'th row and j'th column element in the adjacency matrix. It is defined as the inverse of the harmonic mean of the minimum path length between a particular node and all other nodes in the network, where dði; jÞ is the shortest path length between node i and node j. This is the number of shortest paths between any two nodes that run through a particular node (i). δ mn is the total number of shortest paths from node m to node n and δ mn ðiÞ is the number of shortest paths from node m to node n that pass through node i. A high N bc denotes large impacts of this node on the information flow over the whole network.
within a distance of approximately 16 to 20 mm, revealing that most of voxels affected by the "spatial leakage" are within 1.5 to 2 cm from the seed site.
To specifically estimate the effect of "spatial leakage" in our current study, we investigated individual Euclidean distances across all the 34 AAL ROIs that we used in the network formation and presented the results in Fig. 11. In this figure, x-and y-axes are the identified 34 ROIs, while blue color represents the ROI-to-ROI distance being larger than 20 mm and yellow color indicates the distance shorter than 20 mm. We can see that 22 out of 561 (¼34 × 33∕2) ROI-to-ROI distances are shorter than 20 mm, implying that the "spatial leakage" effect is not too much significant (about 3.9% of total ROIs). However, the "spatial leakage" effect could vary depending on locations of the seed because the light propagation in the human brain is not homogeneous. More rigorous analysis should be performed in future studies to confirm this preliminary conclusion. Fig. 10 (a) Pearson's correlation coefficients, R, between the seed region (marked by the black circle) and all other voxels covered by the optical optodes, were computed and rendered on the 3-D dummy brain template. (b) It plots the R values versus corresponding Euclidean distances from the seed to other voxels. California, Los Angeles. His current research interests include resting-state functional brain imaging, and multimodality neuronal disorders studies and algorithm development.
Mary Cazzell received her BS and PhD in nursing from Marquette University and the University of Texas at Arlington, respectively. She is the director of Nursing Research and Evidence-Based Practice, Cook Children's Medical Center, Fort Worth, Texas. Her research area focused on behavioral measurement of risk decision-making and identification of related neural correlates, using optical imaging, across the lifespan.
Olajide Babawale received his bachelor of engineering degree in electrical and electronics engineering from Covenant University, Nigeria, and his master's degree in biomedical engineering from the University of Texas at Arlington. He is a second-year doctoral student at the University of Texas at Arlington. His current research work includes multimodal neuroimaging using electroencephalographic (EEG) and functional near-infrared spectroscopy (fNIRS) modalities, EEG functional connectivity analysis using source-space imaging approaches, and fNIRS analysis of cognitively impaired children.
Hanli Liu received her MS and PhD degrees in physics from Wake Forest University, followed by postdoctoral training in tissue optics at the University of Pennsylvania. She is a full professor of bioengineering at the University of Texas at Arlington. Her current expertise lies in the field of near-infrared spectroscopy of tissues, optical sensing for cancer detection, and diffuse optical tomography for functional brain imaging, all of which are related to clinical applications.