Transcranial photobiomodulation-induced changes in human brain functional connectivity and network metrics mapped by whole-head functional near-infrared spectroscopy in vivo

: Transcranial photobiomodulation (tPBM) with near-infrared light on the human head has been shown to enhance human cognition. In this study, tPBM-induced eﬀects on resting state brain networks were investigated using 111-channel functional near-infrared spectroscopy over the whole head. Measurements were collected with and without 8-minute tPBM in 19 adults. Functional connectivity (FC) and brain network metrics were quantiﬁed using Pearson’s correlation coeﬃcients and graph theory analysis (GTA), respectively, for the periods of pre-, during, and post-tPBM. Our results revealed that tPBM (1) enhanced information processing speed and eﬃciency of the brain network, and (2) increased FC signiﬁcantly in the frontal-parietal network, shedding light on a better understanding of tPBM eﬀects on brain networks. Our work demonstrates the feasibility of using whole-head fNIRS to study cortical network reorganization induced by tPBM. The analysis methods based on cerebral hemodynamics, Pearson’s correlation coeﬃcients, and graph theory were used to track changes in network connections and graphical metrics between the rPFC, where tPBM stimulation was administered, and other cortical regions. The hemodynamic-, FC-, and GTA-derived results revealed that during the later period of tPBM stimulation, local connections were enhanced, and faster information processing was observed by the shortened pathlength of FC and increases in global (E g ) and local eﬃciency (E loc ) network parameters. During resting state the brain is not idle, so there are neural networks engaged in processing of cognitive information that was previously acquired (e.g., memory consolidation) as well as processing of homeostatic information from internal body functions. After stimulation, distant connections had increased FC strength, particularly for FPN, which supports cognitive ﬁndings showing that tPBM improves executive These observations indicate that there was a time dependence in the change of FC and graphical network patterns during the administration of tPBM. Our ﬁndings demonstrated the excellent and novel applicability of whole-head fNIRS in tandem with the non-invasive tPBM intervention to investigate or characterize cortical network eﬀects. These observed eﬀects may shed light on the underlying large-scale network mechanism and potential application of tPBM for cognitive enhancement in healthy subjects and patients with brain injury and disorders, such as stroke, traumatic brain injury, Alzheimer’s disease, Parkinson’s disease, anxiety, and depression


Introduction
Photobiomodulation (PBM) entails the application of low-power, high-fluence light in the red to near-infrared (NIR) range (usually between 630-1100 nm) to modulate mitochondrial respiration in a non-destructive and non-thermal manner. Transcranial PBM (tPBM) refers to PBM directed at the cerebral cortex with the purpose of enhancing cerebral oxygenation and cognitive function [1][2][3][4][5][6]. Approximately 1-2% of NIR light between 660-940 nm penetrates through various layers of the human scalp and skull and reaches the cerebral cortex, several centimeters below the scalp's surface [7,8]. The purported mechanism of PBM relies on photon absorption by cytochrome-c-oxidase (CCO), which is the terminal enzyme in the mitochondrial respiratory chain that is necessary in cerebral oxygen utilization for energy metabolism. As CCO activity increases, more oxygen is consumed while available metabolic energy increases through simultaneous production of ATP via mitochondrial oxidative phosphorylation [3,4,6]. The resulting increase in CCO activity leads to metabolic and hemodynamic alterations in the brain that can facilitate neuroprotection and cognitive enhancements [9]. In several reviews by Hamblin [6,10], use of tPBM by various groups was reported as a possible intervention for acute and chronic stroke, traumatic brain injury, Alzheimer's disease, Parkinson's disease, depression, anxiety, and cognitive enhancement.
Several neuroimaging modalities have been used for functional brain mapping concurrently with tPBM administration such as electroencephalography (EEG), single photon emission computed tomography (SPECT), and functional near-infrared spectroscopy (fNIRS) [11][12][13][14]. Among them, fNIRS measures noninvasively changes of oxyhemoglobin concentration (∆[HbO]) and deoxyhemoglobin concentration (∆[HHb]) resulting from neurovascular coupling secondary to neuronal activation. Of the few fNIRS studies with concurrent tPBM administration reported to date, the number of fNIRS channels used was limited (≤ 20) and only provided mapping of cerebral hemodynamics over the prefrontal cortex [13,14].
Prior fNIRS studies have demonstrated its reliability and popularity by researchers to investigate resting state functional connectivity (FC) based on Pearson's correlation coefficients (PCC) or graph theory analysis (GTA) [15][16][17][18][19]. PCC quantifies temporal cross-correlations between any available pair of spatially remote cortical/brain regions [18,20], while GTA characterizes topological properties of global and local networks among brain regions. Neuroimaging studies under resting and/or task-based states have taken FC as an index to describe the relationship between hemodynamic activation patterns of different brain regions and thus to reflect the level of functional communication between them [20][21][22][23][24]. On the other hand, GTA-derived graphical metrics can be used to compute alternate measures of statistically significant connectivity strength as an alternate FC analysis approach [25,26]. In particular, GTA features small-world networks that quantify the communication efficiency between nodes can be defined as fNIRS channels or groups of channels clustered within functionally distinct brain regions [27][28][29]. GTA has been applied for the assessment of small-world properties in prior resting state and task-based neuroimaging studies. However, no study to date has explored the effect of tPBM on small-world network properties although such analyses have been applied to other paradigms [15][16][17]19,[30][31].
In this work, we utilized a 111-channel fNIRS system with a large field of view for the first-time to quantify changes in 1) FC pair-wise patterns and 2) brain network topographical metrics induced by tPBM on healthy adult subjects. The tPBM was directed to the right prefrontal cortex (rPFC), encompassing the right frontal polar (rFP) and right dorsolateral prefrontal cortex (rDLPFC). It is known from our prior studies that tPBM induces significant increases of ∆ [HbO] in the PFC during and after tPBM delivery on the right forehead of human participants [5,13,14]. Moreover, a recent fNIRS study demonstrated that high-definition transcranial direct current stimulation (HD-tDCS) increases FC strength and small-world metrics during and after the stimulation using both FC and GTA methods [17]. Complementary to the prior work, the focus of the present study was to examine tPBM-evoked changes of FC pairwise and network patterns during and after tPBM based on 111-channel ∆[HbO] time series. We hypothesized that the application of tPBM would increase cortex-wide FC strength relative to the stimulation site and enhance GTA-derived small-world metrics, despite different neuromodulation mechanisms between tPBM and HD-tDCS. If the hypothesis is proven true, this study will strongly support the possibility of using tPBM as an alternate neuromodulation intervention for future biomedical applications.

Participants
Nineteen adults (5 females, age = 31.7 ± 9.5 years) were recruited for this study. All subjects were without any neurological or psychiatric disorders (self-reported). All experimental procedures, including a written consent required prior to participation in this study, were approved by the Institutional Review Board of the University of Texas at Arlington (IRB# 2017-0859).

Experimental instruments and procedures
A continuous wave (CW) fNIRS system (LABNIRS; OMM-3000, Shimadzu Corp., Kyoto, Japan) was used in this study to measure cerebral hemodynamic responses with near infrared light diode sources (780, 805, and 830 nm) and photomultiplier detectors at a sampling frequency of 10.1 Hz. The optode geometry consisted of 32 sources and 34 detectors with a separation of 3 cm resulting in 111 source-detector channels ( Fig. 1(A)). Anatomical cranial reference points (nasion, inion, left and right preauricular points and vertex) and optode locations were recorded for each subject using a 3D digitizer (FASTRAK, Polhemus VT, USA). Montreal Neurological Institute (MNI) coordinates for each source and detector location were calculated using the statistical parametric mapping NIRS_SPM software, which provided the Brodmann area corresponding to each fNIRS channel [32]. All subjects' MNI coordinates were subsequently averaged together. Thus, the 111 channels covered cortical areas of the following 12 regions of interest (ROIs): left and right frontopolar prefrontal cortex (lFP; rFP), left and right dorsolateral prefrontal cortex (lDLPFC; rDLPFC), Broca's area, left and right premotor cortex (lPMC; rPMC), left and right primary motor and somatosensory cortical (lM1/S1; rM1/S1) areas, Wernicke's area, and left and right somatosensory association cortex (lSAC; rSAC). These 12 regions served as clusters for FC analyses in the following sections. Fig. 1. Experimental set-up and protocol. A) 111-channel layout with twelve regions of interest covered by the optode geometry: frontopolar (FP) (red), dorsolateral prefrontal cortex (DLPFC) (yellow), Broca's area (green), premotor cortex (PMC) (light blue), primary motor and somatosensory cortical (M1/S1) areas (dark blue), somatosensory association cortex (SAC) (pink), and Wernicke's Area (gold). B) 1064-nm laser. C) The experimental protocol randomized placebo (PBO) and tPBM treatment for subjects. For the protocol, there was a period of at least one week between the two experiments to avoid any carry-over effect. The red-colored spots present Channels 3,4,7, and 8 on the right side and Channels of 1,2,5,6, and 9 are on the left side.
Transcranial PBM was administered using an FDA-cleared 1064-nm, CW laser (Model CG-5000 Laser, Cell Gen Therapeutics, Dallas, Texas) ( Fig. 1(B)). The laser's aperture delivered a well-collimated beam with an area of 13.6 cm 2 at a maximum power of 3.5 W and a laser power density of 0.25 W/cm 2 [3,5]. The tPBM was applied by noncontact delivery to the right forehead of each subject at a frontal site (near the Fp2 location of the international 10-10 EEG system). Subjects and experimental operators wore laser protection goggles throughout the duration of the experiment.
Subjects participated in two sessions that were randomized for tPBM or placebo stimulations. If the tPBM session was first, there was a waiting period of at least one week for the placebo treatment to avoid any carry-over effects. A separation of 2-3 days was common before the active tPBM was given if the placebo session was initiated first. During the stimulation period, the CW laser was administered on the right forehead with an output of 3.5 W, whereas placebo had negligible laser output by turning the laser on and then off within 3 seconds. Subjects sat with their eyes closed for 8 minutes of pre-placebo/tPBM period, 8 minutes during a placebo/tPBM period, and 4 minutes of post-tPBM/placebo ( Fig. 1(C)).

Data preprocessing
FNIRS data were preprocessed using Matlab 2012b (MathWorks, Natick, MA, USA) and the open-source package Homer 2.0 [33]. FNIRS data were detrended based on the baseline by a linear least-squares fit that was subtracted from the data [18]. The raw intensity data were then low-pass filtered using a 3 rd order Butterworth filter at a cut-off frequency of 0.2 Hz to remove large portions of physiological noise, including heartbeat (1-1.5 Hz) and respiration (0.2-0.5 Hz) [34] as done in prior work [23]. Principal component analysis was utilized to remove motion artifacts and global hemodynamic fluctuations, such as those from Mayer waves (∼0.1 Hz) [35,36], that may overlap with hemodynamic response frequencies. The first two principal components were removed from all fNIRS channel data in order to remove these global artifacts [33,34]. Channels located near the branches from the middle cerebral artery or the superficial temporal artery and temporal muscle were removed due to large signal contamination [37,38]. These removed channels are shown in gray in Fig. 1(A). Optical density data were converted into changes in oxyhemoglobin and deoxyhemoglobin concentrations relative to baseline, ∆[HbO] and ∆[HHb], using the Modified Beer-Lambert Law with an estimated differential pathlength factor of 6.0 for each wavelength, an estimate used in Homer 2.0 [39]. Changes in ∆[HHb] values were plotted and found to have smaller amplitudes and lower signal-to-noise ratios while maintaining qualitatively similar inverse patterns to ∆[HbO]. Thus, beyond plotting the corresponding ∆[HHb] time series to verify the qualitative data trends, we excluded ∆[HHb] from further data analysis, as done in prior neuroimaging studies [21,[40][41][42]. values were relatively much smaller than ∆[HbO], they were excluded for further FC analyses. These data processing steps were repeated for only ∆[HbO] temporal changes of all 111 channels, which were further used to derive tPBM-induced alterations in human brain FC and network metrics across the human whole cortex. To quantify channel-wise significant differences in ∆[HbO] changes between tPBM and placebo stimulations during and post-stimulation, a paired t-test was performed at each time point at the significance level of p < 0.05 using a two-tailed t-test.

Correlation analysis to quantify functional connectivity
FC was calculated for each period: pre-, during, and post-tPBM/placebo. Whole-head PCC was performed by calculating the Pearson correlation values, R, between each pair of all 111 channels, thus yielding an 111 × 111 matrix per period [18]. Pearson R values were then converted into Fisher's z-scores. A larger correlation R value and thus a larger z-score implied a stronger FC strength. Paired t-tests were performed on the z-transformed R values across subjects. False discovery rate (FDR) was employed for multiple comparisons correction with q = 0.05 and α = 0.05 with only significant results shown [43]. Only significant connections to the rPFC were displayed in topographic images that were generated using BrainNet Viewer, an open-source Matlab package [44].

Graph theory analysis to determine topographical network metrics
Graph theory analysis (GTA) was applied to investigate tPBM-induced changes in global topological network organization derived from the fNIRS measurements of the whole head. In the analysis, channels were considered as nodes and connections between channels were considered as edges of the network. The FC matrix was then thresholded into a binary matrix in terms of sparsity (S), the number of current existing edges divided by the total possible number of edges in the current matrix [15,30,45]. In this study, we selected and quantified seven topological properties/metrics to study network patterns within a range of S levels (S; 0.05 < S < 0.50; increment = 0.05). Five small-world property metrics were calculated including: clustering coefficient (C p ), characteristic path length (L p ), normalized clustering coefficient (γ), normalized characterized path length (λ), and small-worldness (σ). C p is the average of the clustering coefficients of all nodes. C p is a measure of network segregation by signifying the likelihood that two nodes significantly connected to a third node are also significantly connected to each other, thus forming a connected triangular cluster. C p (i) of a certain node i is defined as the following [29]: where N is the total node number, and d ij is the shortest path length between node i and node j. Longer L p values indicate weaker connections [46].
Normalized clustering coefficient, is the mean of all clustering coefficients over all nodes in a network [27,29]: where N is the total node number, and d ij is the shortest path length between node i and node j. Longer L p values indicate weaker connections [46]. Normalized clustering coefficient, γ, is the mean of all clustering coefficients over all nodes in a network [27,29]: is the average of clustering coefficients over all nodes in a network, quantifying the extent of local group formation within a network. C random p is the mean cluster coefficient of matched random networks that preserve the same number of nodes, edges, and degree distribution as the real network [47].
Normalized characteristic path length, λ, is the average of the shortest path lengths between any nodes of the networks and is defined as [27,29]: where L real p is the average of the shortest path lengths between any pair of nodes in the network, quantifying the capability of parallel information propagation within a network. L random p is the mean characteristic path length of matched random networks that preserve the same number of nodes, edges, and degree distribution as the real network [47].
Two additional parameters were quantified to measure the efficiency of small-world networks, namely global efficiency (E g ) and local efficiency (E loc ). They are defined as [27]: where N is the total node number, d ij is the shortest path length between node i and node j, and E(G i ) is the global efficiency of the subgraph composed of the nearest neighbors of node i. E g describes the efficiency of a parallel information transfer in the network, whereas E loc describes how efficient the communication is between the first neighbors of i when i is removed [27]. All global metrics were calculated for each subject at each period (i.e., pre-, during, and post-tPBM) for each stimulation type (placebo/tPBM) using GRETNA [48]. For the stimulation period, GTA metrics were only quantified for the last half period of tPBM (i.e., 5-8 min) to optimally observe its effects since tPBM has shown gradual alterations in cerebral hemodynamics and electrophysiology [5,11,49]. Paired t-tests were used to test statistical (p < 0.05) differences for each global metric between placebo and tPBM interventions for each period. Each global metric was compared at each S level, 0.05 < S < 0.50 at 0.05 increments. All statistical assumptions were verified, including normality. In particular, the channels close to the stimulation site are 2,3,4, and 7 ( Fig. 1(A)); the corresponding ∆[HbO] trends clearly exhibit significant increases during the last 2-4 minutes of tPBM compared to those with the placebo treatment. On the other hand, most channels present non-significant or little change in ∆[HHb] between tPBM (blue solid lines) and placebo (blue dotted lines), except for only those in channels 3 and 4. Thus, we utilized only ∆[HbO] time series for functional connectivity analysis hereafter. Specifically, following the same data processes used for 8 channels shown in Fig. 2, we acquired 111 sets of time-dependent ∆[HbO] alternations during and after tPBM/placebo stimulations from each of the 19 subjects. All of these temporal series of ∆[HbO] permitted our next step to quantify tPBM-evoked changes in human brain functional connectivity and network metrics.

Functional connectivity to the right prefrontal cortex derived with PCC
FC was quantified between all channels within the rPFC to every other channel on the cortex based on PCC. Channels within the rPFC were chosen due to its proximity to the stimulation site during tPBM exposure. Two types of comparisons were made to show differences in FC (1) between tPBM and placebo during and post stimulation periods (Fig. 3) and (2) between three pairs of periods, namely, pre-versus during tPBM, pre-versus post-tPBM, and during versus post-tPBM, under the true stimulation (Fig. 4). Paired t-tests were performed to determine significant differences in FC R-values between tPBM and placebo stimulation. Positive t-values indicated that FC was greater during tPBM than placebo, whereas negative t-values denoted the opposite.
During the pre-stimulation period, there were no significant (p > 0.05) differences between tPBM and placebo. During the 8-min stimulation period, tPBM reduced FC with respect to the placebo treatment, as indicated by the blue lines in Fig. 3(A). There were 5 significant (p < 0.05) The channels closest to the stimulation site are illustrated by a green box around the corresponding plots. The blue-shaded area marks the 8-min stimulation. Channels 1,2,5, and 6 were located on the left side of the forehead, while channels 3,4, 7, and 8 were placed on the right cortex ( Fig. 1(A)). The black lines indicate significant differences in ∆[HbO] values between tPBM and placebo stimulations at respective times when a paired t-test was performed at each time point at the significance level of p < 0.05.  connections in total to the rPFC: rFP to lFP (1 connection), rFP to rDLPFC (2), rDLPFC to lFP (1), and rDLPFC to l PMC/SMA (1). In contrast, during the 4-min post-stimulation period, there was greater FC strength across multiple cortical regions (or clusters) induced by tPBM treatment with respect to the placebo as marked by the red lines in Fig. 3(B). There were 12 significant (p < 0.05) connections in total to the rPFC: rDLPFC to lDLPFC (2 connections), rDLPFC to lPMC/SMA (1), rDLPFC to rPMC/SMA (1), rDLPFC to lM1/S1 (1), and 7 within the rDLPFC.

Graphical network metrics analyzed by GTA
To show significant differences in global graphical network metrics between tPBM and placebo treatments, statistical paired t-tests were used for each of the three periods, for each sparsity (S) level, to investigate the effect of tPBM on seven global graphical metrics derived from GTA: C p , L p , λ, γ, σ, E g , and E loc. During the pre-stimulation period, the results showed no significant (p > 0.05) difference between tPBM and placebo treatments for any of the metrics at any S level.
During the post-stimulation period, the results showed that only σ had significant differences between the two treatments, occurring at S = 0.1 (p = 0.05), S = 0.15 (p = 0.04), and S = 0.25 (p = 0.05) (Fig. 6). Additionally, statistical paired t-tests were used to determine significant differences between pre-and post-stimulation for each GTA metric for each sparsity, S, level with true tPBM. However, the results showed no significant (p > 0.05) differences in GTA metrics for pre-versus post-tPBM.

Effects of tPBM on hemodynamics
Our hemodynamic results (Fig. 2) clearly demonstrated that 8-min tPBM delivered onto the right forehead of human subjects significantly increased ∆[HbO] during the last 4 min in the prefrontal regions, near the tPBM site, as compared to placebo stimulation. Following the same data processing procedure, all temporal series of ∆[HbO] from 111 channels under different (pre-, during, and post-tPBM/placebo) stimulation periods were obtained and used to quantify tPBM-evoked changes in human brain FC and network metrics.
Most of fNIRS-based functional connectivity studies have utilized Pearson's correlation coefficients to quantify the synchronization level between hemodynamic signals across different brain regions at the resting state [18,30,45,50], where intrinsic brain oscillations are major characteristics. The focus of our study was to examine how 8-min tPBM perturbs the hemodynamic rhythms of the human brain from its resting state during and 4 min after the stimulation. The 8-min tPBM may evoke a significant ∆[HbO] increase in certain cortical regions, but the same tPBM may enhance hemodynamic synchronization between other cortical regions. These two types of altered regions may not necessarily be the same regions: One region may show a significant increase in ∆[HbO] amplitude (univariate type of change) while a different region may show a significant correlation in ∆[HbO] temporal waveforms with another region (covariance type of change). In other words, a change of hemodynamic activity in one brain region does not provide information on how that region interacts with another region in a network. To analyze neural networks interactions requires the use of covariance analyses, such as FC analyses, that determine how the activity of a region correlates with the activity of other regions that are part of a neural network [51]. The relationship between these two types of tPBM-evoked significant brain regions will be investigated in our upcoming studies.

Effects of tPBM on functional connectivity
To the best of our knowledge, this study is the first to evaluate the effects of tPBM on whole-head FC in humans using large field of view fNIRS and tPBM simultaneously. By examining the time dependence on FC for each period, this study shows that tPBM has two different effects on FC during and after stimulation. Seed-based PCC was performed to quantify and map FC between all channels within the rPFC, which encompasses the rFP and rDLPFC, and to all other ROIs due to rPFC's proximity to the location of tPBM administration.
During the stimulation period, the FC map exhibited decreased, localized connectivity between the rPFC and the lPFC and lPMC/SMA (Fig. 3(A); Fig. 4(C)). Decreased FC from tPBM stimulation may be due to a regional increase in ∆HbO, as seen in Fig. 2. Very consistent with this study, several prior studies reported that tPBM increases ∆[HbO] during human brain stimulation [5,13,14,49] and regional cerebral blood flow (rCBF) in animal and human studies [12,52,53]. Increased rCBF can be attributed to tPBM's stimulation of the synthesis of nitric oxide (NO), an endogenous vasodilator [54]. In an animal study, tPBM at 808 nm increased NO concentration and rCBF in mice. However, the NO concentration and rCBF did not increase when NO synthase was blocked [53], indicating that NO synthesis is integral in rCBF changes during tPBM stimulation. The decrease in localized FC during stimulation may therefore be partially attributed to the increase of rCBF in the rPFC region only, which caused differences or desynchronization in hemodynamic variability patterns relative to others in non-stimulated ROIs.
The post-stimulation FC maps exhibited enhanced global connectivity between the rPFC and all other ROIs (Fig. 3(B); Fig. 4(A); Fig. 4(B)). These results suggest that post-tPBM hemodynamic recovery influences coactivation between the rPFC and ROIs that are involved in resting state networks (RSN), as also seen in tDCS studies [17,55]. The rFP and rDLPFC are associated with high cognitive and executive functions; their direct connectivity in RSN and other ROIs is related to a top-down modulation of attention and working memory as seen in increased coactivation of frontal and parietal regions [20,24,[55][56][57][58]. The default mode network (DMN), which comprises the rFP and rDLPFC, is a network that is highly active at rest and linked to human cognition, including the integration of cognitive and emotional processing [20,24,57]. The frontal-parietal network (FPN), also known as the central executive network, encompasses the FP, DLPFC, and SAC and is associated with executive functions, attention, and motor control. The FPN is highly integrated with other brain networks like the DMN, motor network, and language network [20,24,57,[59][60][61]. The FPN's integration with other networks is evident in this study with the extensive connection to the motor network, comprising of PMC/SMA, M1/S1 [20,24], and the language network anchored in Wernicke's and Broca's areas [61,62].
Overall, the most significant changes in FC were found in the FPN and regions that belong to other brain networks and are highly integrated with the FPN. These findings are significant since they unambiguously demonstrate that tPBM modulates FC in distinct and remote functional brain networks of the human brain. The preferential involvement of the FPN is consistent with cognitive findings showing that tPBM improves executive function [63].

Study of cortical network changes during tPBM
To investigate effects induced by a brain stimulation, such as tDCS and tPBM, a common approach is to make a direct comparison of either behavioral or neurophysiological outcomes between preand post-stimulation for each individual, followed by a statistical test at a group level. In order to better understand the mechanism of brain stimulations, however, it is important to monitor neurophysiological progression and alterations during respective stimulations with respect to the pre-stimulation conditions. It is often technically challenging, unfortunately, to record brain signal changes concurrently during the stimulations because of potential electromagnetic interference or limited room on a human subject's head to set up concurrent stimulation and sensing devices. In our study, fNIRS measurements did not interfere with tPBM light delivery, nor were fNIRS measurements interfered by tPBM. We were therefore able to record hemodynamic changes during and after tPBM to help improve our current understanding of the evolution and immediate post-stimulation effects of tPBM over many cortical regions.
The potential contamination from the 1064-nm laser used during tPBM stimulation to concurrent fNIRS recording is very low because of low quantum efficiency of the photo multiplier tube (PMT) used in our fNIRS system, LABNIRS. All detectors equipped in LABNIRS utilize multialkali PMTs [64], which have a quantum efficiency less than 0.1% at any wavelength equal to or longer than 900 nm [65]. Thus, LABNIRS detectors would not have any ability to sense or detect any 1064-nm light during tPBM.

Effects of tPBM on global connectivity network metrics
GTA was used to quantify network connectivity changes induced by tPBM across the entire cortical area mapped by fNIRS. TPBM induced significant changes during the 2nd half of stimulation period (5-8 min) (Fig. 5) and post-stimulation (Fig. 6) across the entire cortical network that can be characterized by the network's small-world features, including five small-world parameters: clustering coefficient (C p ), characteristic path length (L p ), normalized characteristic path length (λ), normalized clustering coefficient (γ), and small-worldness (σ) [26][27][28]46]. During the 5-8 min interval of tPBM stimulation, L p was significantly lower than the placebo treatment at several S levels ( Fig. 4(C)). L p signifies the ability to rapidly combine pieces of specialized information from different brain regions or functional integration [25,46]. Shortened L p indicates that during the 5-8 min tPBM interval, there is enhanced or faster information/integration processing [66,67]. Economic small-worldness, σ, was exemplified (Fig. 5) after tPBM with σ >1 [25]. However, σ is significantly lower after tPBM stimulation compared to placebo. These results are similar to an EEG study that showed brain networks with reduced levels of small-worldness after tDCS stimulation [68].
Two additional parameters, global (E g ) and local efficiency (E loc ), were also calculated to determine how efficiently small-world networks exchange information. TPBM treatment resulted in significantly higher efficiency metrics for both E g and E loc during the stimulation. E g was significantly higher than the placebo due to its susceptibility to shorter L p , as also seen in this study. A high E g assures effective functional integration of information between and across remote regions [27,46,67]. During the 5-8 min tPBM stimulation interval, E loc was also significantly higher than the placebo. E loc measures the ability of fault tolerance of a network, or how efficient communication is between nodes when one node is removed [27,28,46]. Typically, a small-world network is characterized by high E g and E loc [27] but this is not seen after tPBM stimulation. This is in contrast with a prior fNIRS study that noted enhanced E loc after anodal tDCS, but not during stimulation as seen in this study [19]. This difference could be due to different underlying mechanisms affecting neuroplasticity between tPBM and tDCS.
Overall, two different analysis methods (i.e., PCC and GTA) provided different aspects or views of tPBM-induced effects on human brain FC. During tPBM stimulation, tPBM resulted in desynchronization of hemodynamic oscillations between the intervention site and non-stimulated ROIs, while tPBM enabled faster functional integration/communication and improved efficiency of small-world brain networks. During the 4-min post-tPBM period, while the cortex gradually recovered to its RSN, FC in the FPN became significantly enhanced; connections from the prefrontal cortex to other distinct and remote ROIs, such as Broca's and Wernicke's Area, were also augmented in the human brain.

Understanding of the tPBM mechanism of action in the whole head
Our previous studies have shown experimental evidence that tPBM is able to enhance oxidized CCO concentration and ∆[HbO] near the stimulation site during and after 8-min tPBM [5,13,49]. It was concluded that this increased metabolic and hemodynamic activity near the light delivery site subsequently results in more oxygen consumption, simultaneous mitochondrial energy production via oxidative phosphorylation, and increased regional cerebral blood flow (rCBF) [3,4,6]. It was unknown whether such a local metabolic/hemodynamic stimulation would affect global cortical networks. In this work, for the first time, we demonstrated that tPBM with 1064-nm laser can (1) increase hemodynamic FC significantly in the frontal-parietal network and (2) enhance information processing speed and efficiency of the brain network. Based on these results, we propose the mechanism of action for tPBM in the whole head as follows.
Step 1: tPBM increases rCBF during the stimulation, as expected in line with prior animal and human studies [12,52,53]. During this period, the observed decreases in FC (see Figs. 3(A) and 4(C)) can be attributed to an imbalance and desynchronization between rCBFs from the treated rPFC and other non-stimulated ROIs.
Step 2: Toward the end of stimulation, significant increases in ∆[HbO] at the stimulation site (Fig. 2) have progressively reached balance and synchronization with those in other ROIs through the brain networks across the human head. Thus, during the post-stimulation time, we observed significant improvements in global FC (Figs. 3(B), 4(A), and 4(B)) and reduced small-worldness (σ) (Fig. 6). These results suggest hemodynamic coactivation between the rPFC and other ROIs that are involved in the same resting state networks, as also seen in tDCS studies [17,55].

Limitations and future work
A few limitations should be noted for this study. First, fNIRS is limited to measuring only cortical brain areas and has a lower spatial resolution in comparison to fMRI [33,34]. Second, this study would have improved FC and GTA metrics during the post-stimulation period if subjects were measured for at least 5 min after tPBM stimulation. The time was shortened in consideration of the subject's comfort while wearing the fNIRS optodes. However, an additional 1-4 min would likely not have put undue burden onto the subject. Typically resting-state FC scans are most reliable when measured between 5-13 min [50,69]. Third, the neuromodulation by tPBM involves a relatively slow hemodynamic process that is progressive [5,11] and thus could not be rigorously represented or quantified by either FC or GTA, both of which are more appropriate for analyzing resting state FC without much external intervention during a period of 5 min or longer. Lastly, handedness of the subject was not taken into consideration. However, prior studies have shown that there are differences in brain connectivity patterns between left and right handed subjects and should be considered in the future, as in the studies reported in Refs. [23,70].
For future work, more advanced methods to analyze time-varying dynamic FC need to be explored to account for time-evolving features of neuromodulation by tPBM. Next, given the complex effects of tPBM on the brain, future studies may be better informed by multi-modal imaging. For example, since rCBF is expected to change with tPBM, concurrent fNIRS measurements with diffuse correlation spectroscopy or arterial spin labeling fMRI would be helpful to map the spatiotemporal extent of these changes [71]. In addition, concurrent mapping of fNIRS with electrical activity from the neuronal networks with multi-channel EEG would be useful in deconvolving the tPBM-induced relative changes in hemodynamics versus cerebral electrophysiological activity [11].

Conclusion
Our work demonstrates the feasibility of using whole-head fNIRS to study cortical network reorganization induced by tPBM. The analysis methods based on cerebral hemodynamics, Pearson's correlation coefficients, and graph theory were used to track changes in network connections and graphical metrics between the rPFC, where tPBM stimulation was administered, and other cortical regions. The hemodynamic-, FC-, and GTA-derived results revealed that during the later period of tPBM stimulation, local connections were enhanced, and faster information processing was observed by the shortened pathlength of FC and increases in global (E g ) and local efficiency (E loc ) network parameters. During resting state the brain is not idle, so there are neural networks engaged in processing of cognitive information that was previously acquired (e.g., memory consolidation) as well as processing of homeostatic information from internal body functions. After stimulation, distant connections had increased FC strength, particularly for FPN, which supports cognitive findings showing that tPBM improves executive function [63]. These observations indicate that there was a time dependence in the change of FC and graphical network patterns during the administration of tPBM. Our findings demonstrated the excellent and novel applicability of whole-head fNIRS in tandem with the non-invasive tPBM intervention to investigate or characterize cortical network effects. These observed effects may shed light on the underlying large-scale network mechanism and potential application of tPBM for cognitive enhancement in healthy subjects and patients with brain injury and disorders, such as stroke, traumatic brain injury, Alzheimer's disease, Parkinson's disease, anxiety, and depression [10,14,52].

Funding
National Institute of Mental Health (RF1MH114285).