1. Introduction

The human brain is seen as a complex system due to its functional capabilities and structural organization. Indeed, complex systems’ tools have been successfully employed to neuroimaging data recently, which has allowed researchers to explore intrinsic features of the brain. Of major interest in the literature is the brain’s functional interaction between different regions, which can be investigated with graph theory approaches [1–4]. Recent studies from Blood Oxygen Level Dependent (BOLD) functional Magnetic Resonance Imaging (fMRI) has presented remarkable features of brain networks, such as highly connected hubs [5], hierarchy [6], assortativity [6,7], and small-worldness [8].

In particular, the study of the human brain at rest has gained attention mainly because of its potential for diagnosing the diseased brain [9–13]. In the resting state, the brain lacks external input, so the subjects do not have to perform any task, avoiding the subjective need of computing effort and task demand among subjects and groups. Cerebral hemodynamic fluctuations of the human brain at rest, also known as spontaneous brain fluctuations, are characterized by low-frequency (< 0.1 Hz) oscillations [14–18]. Although the origin of these fluctuations is still subject of debate, most of the resting-state patterns reported so far are observed among brain areas that share the same functionality and/or that are anatomically connected [14,19–21]. This fact supports the hypothesis that resting state networks (RSNs) reflect at least part of neuronal activity [19].

The first report regarding resting state functional connectivity (rs-FC) dates from 1995 in which Biswal et al. correlated the BOLD signal from the right and left somatosensory cortices [14]. Since then, fMRI has been employed as a standard technique to explore rs-FC mainly due to its high spatial resolution and capability of covering the entire brain, including deep areas. A variety of procedures were developed to assess rs-FC from fMRI data, such as seed-based correlation and Independent Component Analysis (ICA) [13,22–25]. In the seed approach, an element or a region of interest (ROI) is arbitrarily chosen as a seed, and a similarity measurement (typically, the Pearson correlation coefficient) between the seed’s time series and the time series of all other regions is performed. This simple yet powerful approach for investigating rs-FC has the advantage of providing a direct interpretation of the results and a good sensitivity to identify RSNs [26]. Its main limitation is the need of a priori information (by choosing one ROI as a seed), which can bias the analysis to the choice. On the other hand, the ICA method does not require a priori information [25]. The ICA approach consists of decomposing the temporal signal onto independent components, from which the interaction between the regions can be inferred [24,27]. In the last years, ICA has become standard in fMRI experiments to investigate RSNs.

Functional connectivity analysis methods have been previously translated and adapted to near-infrared spectroscopy (NIRS) data [28–35], showing that NIRS can also probe spontaneous brain fluctuations in human subjects during the resting state. Due to its high temporal resolution, high portability and low-cost, NIRS can additionally provide bedside and/or longitudinal monitoring of rs-FC in both the healthy and the diseased brain [36–39]. The optical technique can also assess populations that would be hardly possible and significantly expensive to measure with fMRI, such as neonates [38,40]. In addition, NIRS high temporal resolution allows for physiological signal removal quite efficiently [28,29]. On the other hand, NIRS data can only provide information from the cortical regions due to the principles of photon diffusion through tissue. Even in the cortex, the relatively low number of probes in most NIRS systems further limits analysis of the whole brain and makes it difficult to understand the interaction among different cortical regions. Very few studies so far have presented a detailed analysis about resting state functional connectivity of the whole brain using NIRS [29,41–43], and at least some of these studies suggest a high variability, both across different subjects and even for different data sets of the same subject. Even fewer studies attempted to employ complex networks and graph theory algorithms to analyze NIRS data [41,43]. Unlike fMRI, the translation of graph theory to NIRS is not properly established in the literature and lacks further investigation.

In this work, we hypothesized that the use of network algorithms could be employed to further investigate the human brain functional connectivity during the resting state with NIRS data. By properly taking account of NIRS specific features, we aimed to characterize the structure of the interaction among functional NIRS measurements of the whole head with global network parameters. We also analyzed how these parameters varied as function of the strength of the interaction. We considered that the behavior of networks built from NIRS data could provide a more global picture about cortical connections than the standard seed-based approaches. Last, we wanted to investigate whether there were common connectivity patterns among different subjects, and we developed a novel methodology based on the frequency of appearance of links in a network to find out similar features across different NIRS-based networks. We found that global network parameters extracted from NIRS are reproducible and therefore reliable to characterize the healthy brain at rest. Similar to previous findings with seed-based approach in resting state NIRS, networks extracted from total-hemoglobin concentration changes presented the highest reproducibility across different runs and subjects.

2. Materials and methods

2.1 Subjects and protocol

The experimental data used in this study was previously reported [29,44]. Briefly, data were acquired from 11 healthy-male subjects with an average (standard deviation) age of 35 (12) years old. All subjects were instructed to close their eyes and do nothing (resting state) in a recliner chair inside a quiet dark room. For each subject, 300-sec baseline runs were performed from 2 to 4 times. The experiment protocol was approved by the Institutional Review Board at the Massachusetts General Hospital, where the experiments were carried out. All subjects provided written informed consent.

2.2 NIRS methods

A continuous-wave near-infrared spectroscopy system (CW5, TechEn Inc., Milford, MA) was employed to perform all measurements [45]. The system consisted of 32 laser diodes at 2 different wavelengths (690 and 830nm, emitting ~10mW each) and 32 avalanche photodiodes (APD). The optical probe was designed to cover most of the subject’s head in a configuration that allowed measurements of 48 different source-detector pairs (channels). Each channel had a source-detector distance of 3.0 cm. All optodes were secured with Velcro and foam material. The 48 channels covered frontal (16 channels), parietal (12 channels), temporal (13 channels), and occipital (7 channels) lobes. Blood pressure was independently and simultaneously monitored with a homemade pressure sensor for further background physiological noise removal. The temporal resolution of all (optical and blood pressure) data was 10 Hz, which provided temporal time series of 3000 data points per channel in each run.

For each run of each subject, we first removed motion artifacts by employing an automated algorithm from a standard NIRS data analysis package (HomER) [46]. Briefly, any time point was marked as motion artifact if the signal changed 10 times the standard deviation of the time series, or if it changed by more than 20%, over a short period of time (0.5 s). Motion correction was performed by using a spline interpolation technique [47,48]. Then, channels with signal-to-noise ratio (SNR) less than 25% of the mean SNR of all channels were discarded and not included in the analysis. We discarded 2 subjects from analysis due to this procedure, since they had a significant lower number of channels with good SNR.

For rs-FC data analysis, light intensity for each channel and wavelength was first converted to relative changes in optical density. Then, we performed a global signal regression by applying a Principal Component Analysis (PCA) algorithm and filtering out the first principal component (PC), which has been shown to have the highest correlation with the global average signal [44,49]. The regression of the PC-based global effect estimator minimizes the effects of global artifacts, mainly the ones from systemic physiology (mostly due to cardiac pulsation, respiration, blood pressure and heart rate) [44]. In addition, it has been shown that the PCA regression does not introduce spurious anti-correlations unlike the general global average signal [49]. Next, the time courses were band-pass filtered between 0.009 and 0.08 Hz so that low-frequency interference noise could also be removed [30,50,51]. This low-pass frequency cutoff also removes high frequency physiological noise due to heart beat (~0.8 Hz) and cardiac cycles (~1 Hz). Last, we regressed out the measured blood pressure fluctuations from the filtered signal to remove systemic low-frequency blood pressure contributions, as previously published [29,52,53]. Due to the importance of this step, we discarded 2 other subjects from analysis due to bad blood pressure data. After all regressions of spurious sources in the optical signal, the resultant time series was converted to oxy-hemoglobin (HbO) and deoxy-hemoglobin (HbR) concentration changes by using the modified Beer-Lambert law [46]. Total hemoglobin concentration (HbT) was calculated as the sum of HbO and HbR.

2.3 Network analysis

The NIRS networks were defined using the source-detector pairs and the time courses of the hemoglobin concentrations. Figure 1 depicts the whole process of network construction. Briefly, each channel represented a node in the network. For each NIRS contrast (HbO, HbR, and HbT) we computed the Pearson correlation coefficient, r, between two channels for all possible combinations. The correlation matrix from a single run of a subject was used to define the existence of a link between two nodes of the network [41].

 

Fig. 1 Procedure to construct the topological resting state NIRS networks. (A) Optical geometry of the NIRS experiment showing sources (red crosses), detectors (blue circles) and channels or source-detector pairs (full black circles). For each hemoglobin concentration time-course, we computed the Pearson correlation coefficient, r, among all channels to extract (B) the correlation matrix of a single run of a subject. (C) The adjacency matrix was obtained after choosing a correlation coefficient threshold, r_T. (D) The network was backprojected onto a topographic map of the optical geometry for visualization. (HbO: oxy-hemoglobin, HbR: deoxy-hemoglobin; HbT: total hemoglobin).

Download Full Size | PDF

The adjacency matrix, A, was obtained by binarization of the correlation matrix. In the adjacency matrix, two nodes i and j were considered to be linked if the correlation coefficient, a_ij, was higher than or equal to an arbitrary threshold, r_T (i.e., a_ij ≥ r_T). Since there is no standard choice for the threshold, we varied r_T from 0.05 to 0.90 (in steps of 0.05) so that we could analyze the behavior of each network as function of the threshold. For each threshold, the network was backprojected onto a topographic map for visualization purposes.

For network analysis, global network parameters were extracted using standard metrics from graph theory [1,3,4,6]. For each network we calculated: 1) the average number of links per node, known as the average degree (K) of a network. In terms of the adjacency matrix, K can be expressed as $\frac{1}{n}{\displaystyle \sum _i{\displaystyle \sum _j{{\displaystyle a}}_{ij}}}$; 2) the dispersion of the links, measured as the standard deviation of the degree for each node (StdK); 3) the clustering coefficient (CC) in order to provide a metric of how grouped the neighboring nodes are. Specifically, the CC measures the fraction of the collection of three nodes that form a loop, and it can be calculated from the elements of the adjacent matrix as [1]

The network parameters were calculated independently for each run of each subject and each NIRS contrast. Statistics within a subject (intra-subject analysis) were summarized using the mean and the standard deviation of the network parameters across all runs of the same subject for each contrast separately. Similarly, results for the whole group (inter-subject analysis) were performed averaging over each subject’s individual network metrics results. All data analysis and statistics were performed with homemade scripts in MatLab (Mathworks, Inc., Natick, MA).

3. Results

Before analyzing the networks, we further replicated results previously published in the literature with the seed approach [29]. Figure 2 shows HbR correlation matrices (Fig. 2(A)) and their topological resting-state functional connectivity maps (Fig. 2(B)) for 3 different runs of the same subject. For the topographic maps we arbitrarily chose a channel in the parietal lobe as a seed.

 

Fig. 2 Resting-state functional connectivity maps performed with the seed-based approach. A) Correlation matrices derived from deoxy-hemoglobin concentration changes in three different runs for one illustrative subject. The black rectangle highlights the seed chosen for analysis. B) Backprojected functional connectivity topographic maps for the seed chosen in all three runs. The black lines represent the source-detector separation (i.e., channels). The thicker line represents the channel picked as the seed.

Download Full Size | PDF

Although RSNs have been shown to be reproducible with NIRS data, Fig. 2(B) suggests that the connectivity maps can also change considerably over time. As an example, for two specific channels the correlation coefficient varied from −0.86 in the second run to 0.97 in the last run. Across all channels, the average change in correlation was 0.52 between the first and the last runs. The mean correlation coefficients and their dispersion for all subjects in this cohort have been previously reported in the literature, with variations of up to 42% [29]. Although motion artifacts [54,55] and sample size [56] can significantly contribute to contamination in the correlation coefficients (and therefore to large variability in the correlation maps), it is unlikely that the variation in our NIRS-based connectivity maps can be attributed to these effects only due to our preprocessing steps and number of data points in each time series (~3000). The heterogeneity reported in Fig. 2 motivated the search for alternative approaches that could yield less variability in order to better understand brain spontaneous fluctuations at rest, as measured by NIRS.

Given the intra-subject variability in the connectivity maps with the seed-based approach, we first attempted to investigate the NIRS-based network during the resting state over time (i.e., across different runs for the same subject). Figure 3 shows the number of links, N, as function of the correlation threshold for each run of a representative subject. As expected, the number of connected links decays exponentially as the strength of the interaction (as measured by the correlation coefficient) between the regions increases. For all NIRS contrasts, we observed a transition in the behavior of the total number of links as function of the correlation strength. This transition occurs at r_T = 0.7 for HbO and HbT, and at r_T = 0.6 for HbR. For correlation thresholds higher than these transition points, small increases in the threshold make the networks abruptly sparser until the point there is no connection between regions. Interestingly, such behavior is constant across runs for the same subject, and the transition point is very similar for all NIRS contrasts and across all subjects, suggesting that the general behavior of the links in the networks (and therefore the network structure) is approximately constant despite local variations.

 

Fig. 3 Number of all connected links (N) in the network as function of the correlation threshold, r_T, for all NIRS contrasts of a representative subject. In the figure, each color represents a different 5-minute run from the same subject.

Download Full Size | PDF

In order to investigate the reproducibility of the network at different time intervals, we computed the network parameters for each run of each subject. The results for a representative subject can be seen in Fig. 4. As expected, K and CC are inversely proportional to the correlation threshold, since the number of links decreases exponentially with increasing threshold. The distribution of the links is highly nonlinear, as indicated by StdK. We typically observe a high dispersion of links until r_T = 0.3, potentially due to the presence of noise in lower correlation coefficients. After this threshold, the distribution of links decreases monotonically. Regarding the behavior of the network for long-range interactions, depicted by the diameter and the average pathlength, we observed that both parameters increase with the correlation threshold, as also expected since the network gets sparser with increasing correlation threshold.

 

Fig. 4 Average degree (K), standard deviation of the degree (StdK), clustering coefficient (CC), characteristic pathlength (L) and the diameter (D) as function of the correlation threshold for a representative subject over four different runs. All network parameters are highly reproducible over time.

Download Full Size | PDF

Although the intra-subject variability in functional connectivity maps with the seed-based approach can be quite high (Fig. 2), the network parameters were significantly less variable than the correlation coefficients across the different runs. We quantified the variability across runs of the same subject by calculating the standard deviation of the network parameters over all the runs of the subject for all correlation thresholds. The maximum standard deviation of the average degree for the subject in Fig. 4 was 15% (at r_T = 0.3 in the HbO network). Similarly, the CC showed a maximum standard deviation of 18% (at r_T = 0.55 in the HbO network), while the maximum standard deviation of L and D (whose values are at least 1) was 0.32 (at r_T = 0.3 in the HbO network) and 1.2 (at r_T = 0.3 in the HbO network), respectively.

The low variability in the network parameters was found in all of our subjects. The greatest standard deviation across runs of the same subject was 25% for both K (at r_T = 0.6 in the HbO network) and CC (at r_T = 0.45 in the HbR network). The highest standard variation of L and D were 0.34 (at r_T = 0.55 in the HbR network), and 1.5 (at r_T = 0.4 in the HbO network), respectively. Among all NIRS contrasts, we found that HbO has the lowest reproducibility (i.e., larger standard deviations across runs) for K and D in the intra-subject analysis. HbR presented the lowest reproducibility for the CC.

We also hypothesized that the stability of the network parameters still holds across different subjects. Figure 5 shows the mean and standard deviation of the network parameters for all NIRS contrasts and for all subjects. Similarly to what we saw for individual subjects, the behavior of the network parameters as function of the threshold is similar for the whole group. We did not observe significant differences in how the network parameters vary as function of the threshold between the NIRS contrasts. As expected, the variability among different subjects is higher than the intra-subject variability, but still surprisingly low. Table 1 summarizes the largest variability found for each of the global network metrics. One can note that HbT provides the lowest range and standard deviation among different subjects for most of the network metrics we calculated, similar to what have been previously reported with seed-based approaches [29]. Overall, the global network parameters appear to be a reliable metric to characterize functional connectivity during the resting state with NIRS.

 

Fig. 5 Network parameters for HbO (red circles), HbR (blue squares) and HbT (green triangles) data as function of the correlation threshold across all subjects. The error bars represent the standard deviation across all runs of all volunteers.

Download Full Size | PDF

Table 1. Variability of the network parameters across all subjects. Variability was quantified by the largest range and standard deviation of the parameters at a given threshold, r_T. The maximum variation represents the ratio of the standard deviation and the mean value.

View Table

Last, we were interested in analyzing the common links of the networks across different subjects in order to find NIRS-based network patterns for the group. To do so, we arbitrarily increased r until the average degree of the network was 20, which corresponds to approximately 40% of the maximum average degree. We arbitrarily chose this average degree for illustrative purposes, since we wanted to visually analyze the topographic network. We repeated this process for every run of each subject. The resultant network of each subject was created with the links that repeated in at least 67% (or 2/3) of all runs for that subject. This choice guarantees that we would always look at the most frequent links, regardless of the number of runs for each subject.

Figure 6 shows the resultant network with the most frequent links across all subjects for each NIRS contrast. Overall, we found a total of 68, 40 and 67 connections common to all subjects (i.e, connections that appeared in at least 6 of 7, or 85%, of all subjects) for HbO, HbR and HbT networks, respectively. Such number of connections represents an average degree of approximately 3, indicating that the common pattern must not be densely connected and therefore information transfer should be efficient across most of the brain regions. One quantitative measurement of such efficiency in networks can be performed by the giant component [57,58]. Briefly, the giant component of a network is the largest subgraph of the network in which starting from an arbitrary node it is possible to reach any other. HbO, HbR and HbT resultant networks had giant components of 75%, 41% and 81% of all nodes, respectively. Interestingly, HbT networks had the highest reliability in terms of global network parameters (both inter- and intra-subject analysis) as well as the largest giant component among all NIRS contrasts.

 

Fig. 6 Most frequent links across all subjects for all NIRS contrasts. The colorbar indicates the normalized frequency of each link across all subjects. Black dots represent the nodes (i.e., the channels) of the network. Hubs in the networks are represented by the bigger gray dots.

Download Full Size | PDF

Regarding location of the most frequent connections, all NIRS contrasts exhibited links between symmetrical hemisphere regions. The parietal and occipital lobes have approximately 6.5% and 28% of all possible inter-hemisphere connections linked, respectively. The number of regions connected in the frontal lobe is lower, with only 3.1% of all possible inter-hemisphere connections linked in the resultant network.

Last, we also questioned whether there were highly connected nodes (hubs) that could affect or drive NIRS networks at rest. In order to find the hubs in each network, we first computed a topological measure of centrality based on the weighted degree, w, of each node, and then we found the nodes with w ≥ 90% of the maximum w of each network [59]. Here, the value of 90% was arbitrarily chosen considering previous studies of centrality [59,60]. We found 4, 5 and 8 hubs in the HbO, HbR and HbT networks (gray nodes in Fig. 6). Most of the hubs were located in the frontal and parietal lobes, with a slightly predominance in the left hemisphere.

4. Discussion

In this work, we aimed to present and to validate a novel approach based on graph theory to deeply investigate the interactions of spontaneous hemodynamic activity between different regions of the brain during the resting state with NIRS. Several studies – mostly with fMRI data – show evidences that hemodynamics throughout the brain is dynamically connected through some sort of network structure. Graph theory has recently been employed as a useful tool to analyze the topology of functional connectivity brain hemodynamics in fMRI [8,21,61,62]. Similarly, networks from human brain electrophysiology have also been previously analyzed with magnetoencephalography and electroencephalography [63–65].

Distinct from most BOLD-fMRI works, NIRS offers a direct measurement of hemoglobin concentration changes. However, NIRS measurements are limited by the short depth of light penetration – it ranges from ~0.5 to 2 cm in adult humans, – which is just enough to probe the cerebral cortex [66]. Since most of resting state networks reported with BOLD-fMRI has cortical activation [27], NIRS networks could also be employed to investigate interactions between hemodynamic regions during the resting state. Unlike most of NIRS studies in the literature, we took advantage of a large spatial coverage from a previous data set to get hemodynamic information from the whole head. The spatial extension of our probe was crucial to analyze the interaction among different regions of the brain with the network approach. On the other hand, it does not allow measurements of short channels to be used as regressors for removal of extra-cortical contributions. However, we have shown in the past that by independently measuring systemic signals, such as blood pressure oscillations, most of systemic bias present in the NIRS signal could be removed [29].

Few studies so far have used NIRS to study rs-FC considering the whole head [29,41,43], and even fewer combined NIRS data and graph theory [41,43]. A previous study found that NIRS-based networks are as efficient as random networks to globally process information [41], which is in agreement with our high giant components estimates (75%, 41% and 81% for HbO, HbR and HbT networks, respectively). The same authors showed previously that the network parameters were reproducible among subjects, similar to what we have found here (Fig. 5).

Our present study extends on previous findings by showing that the network parameters are also stable within subject variability, as opposed to variations previously reported in rs-FC with NIRS using the seed approach [29]. Some authors have previously discussed that the seed correlation approach is biased to variations in the seed location and size [27]. Since NIRS lacks a high spatial resolution, different runs could be highly affected by variations in the optode displacement during the experiment, and even more across subjects. Figure 2 shows that these variations do exist, and they can be quite large. In addition, compared to most of rs-FC approaches performed with NIRS so far, including ICA, the network approach has the advantage of characterizing global information and structure with well-developed and highly explored global network metrics. By using the network approach we could find a remarkably and consistent dynamics of global network parameters as function of the strength of coherence at every run for each subject (Fig. 4). Despite large local variations over time, global network features seem to remain stable over time. Over all NIRS contrasts, HbT provided the most reproducible networks both within subjects and across subjects. However, given the small sample size in this pilot study, future validation of the technique in a larger group of healthy subjects is needed. It would also be interesting to compare the global network metrics measured in healthy subjects with the parameters of the diseased and injured brain. We are inclined to believe that the powerfulness of network analysis coupled with some intrinsic features of NIRS enriches the potential to investigate the human brain at rest and to future diagnose the diseased brain.

The similar behavior of the network across runs can be viewed as a conservation of the network properties over time. The seed-based topographic map can be seen as a ‘picture’ of the system at a certain time from interactions of one particle with the others, similar to a microstate in a many-particle system. Different times will show different ‘pictures’ of the system, from different microstates. The global network properties define the macrostate of the system, which can be assessed by different microstates. This statistical perspective is an alternative approach to explain the local variations in seed-based approaches as well as the global stability of the network properties.

Regarding the network methodology, most network approaches published in the neuroimaging literature are built with a fixed parameter (in most cases, it is used the correlation threshold or the average number of links) [4,5,63]. Here, we chose to study the behavior of the networks as function of a free parameter. Therefore, we treated the strength of interaction as a degree of freedom of the system, which is mathematically represented by the Pearson correlation coefficient. The threshold can be seen as the order parameter of our NIRS-based networks. Indeed, Fig. 3 shows an abrupt change in the number of links of the network as function of threshold. After the threshold reaches 0.6-0.7, the slope of the decay changes abruptly, suggesting a different behavior of the network after this point.

Last, we wanted to further understand the topology of the interactions between the different regions in our data. Despite the fact that the connectivity pattern can be unique to every subject (and can even change over long periods of time due to brain plasticity), we hypothesized that there must be common links across our population cohort. In order to find the most repeatable links we introduced a simple way to quantify the frequency of each link. While our approach enabled us to find the most connected regions and the hubs of the network, it was based on arbitrary assumptions for the parameters. Due to the small number of subjects, we chose to be very conservative in our choices, so that our procedure could indeed represent a common network across all volunteers. More importantly, we certified that our choices did not add bias in our results, since changing the values of the parameters did not remove any of the links that we reported in this study. An alternative choice of parameters added more links that were not presented in our choice. Nevertheless, these assumptions are limited due to our sample size, and they could be adjusted accordingly in future studies with a larger number of subjects.

The frequency map of common links across the population presented in Fig. 6 is elucidative because it provides interesting insights from NIRS-based functional connectivity at rest. First, we have found that most of the hubs are located in frontal and parietal regions, and they are not symmetrically distributed between the hemispheres. Similar to previous studies with BOLD-fMRI [67], there are more hubs in the left hemisphere. Brain asymmetry is a subject of debate in current neuroscience, and we need more studies to better understand and validate our results, but it is interesting to find this asymmetry as a consequence of our approach. In addition, we observed that several of the most frequent links connects regions that resemble the most published networks with BOLD-fMRI during the resting state [14,20,21,24]. However, the lack of spatial resolution in our data does not allow us to properly compare our maps of frequency with the BOLD-fMRI networks. A future study that simultaneously combines NIRS and fMRI is currently ongoing in order to further understand the similarities between our NIRS-based derived networks and the standard fMRI-based networks.

5. Conclusion

In summary, we have investigated network properties of the human brain at rest by applying network theory on NIRS data from the whole head, including frontal, parietal, temporal and occipital lobes from both hemispheres. We found that global network parameters extracted from NIRS data are a reliable metric to characterize the human brain organization at rest. Among all NIRS contrasts, our results showed that the HbT provides the most robust networks with the highest reproducibility over time and among different subjects. When we analyzed the most frequent links across all volunteers, we found connections among different brain areas that are in agreement with the fMRI literature. Brain areas that are symmetric located have high density of links and long clusters. Most of the connected regions are located in the frontal and parietal regions. Overall, our present work enhances the feasibility of NIRS and networks to investigate resting state connectivity of the whole head. It also contributes to the understanding of how the brain works at rest and how it is organized as a complex network. Future studies with a larger number of healthy subjects are desirable to test some of the arbitrary choices we made in this work. The powerfulness of network analysis coupled with some intrinsic features of NIRS, such as portability, low cost, and high-temporal resolution, enriches the potential to investigate spontaneous brain activity and to future diagnose diseased brain.