Neuromodulation of brain power topography and network topology by prefrontal transcranial photobiomodulation

A total of 49 healthy human subjects (30 men and 19 women; 26 ± 8.8 years of age) were enrolled from the local community of the University of Texas at Arlington. The experimental protocol was approved by the Institutional Review Board of the University of Texas at Arlington and complied with all applicable federal and National Institute of Health (NIH) guidelines. Written informed consent was obtained from each participant before starting the first measurement. Four participants were removed from the dataset because of self-reported or observed sleepiness during the measurement, resulting in 45 participants being included for data analysis. The participants were instructed to refrain from consuming caffeinated drinks for at least three hours before each experiment.

In this study, tPBM and sham experiments were performed using a continuous-wave laser at 1064 nm (Model CG-5000 Laser, Cell Gen. Therapeutics LLC, Dallas, TX, USA), which was cleared by the Food and Drug Administration. The laser was delivered to each participant's right forehead, with an aperture of 4.2 cm in diameter and a period of 8 min. A sham experiment was performed with the laser power set to be 0.1 W and the laser aperture covered with a black cap. Table 2 lists the key parameters of light delivery used in the study. Note that a penetration rate of 1%–2% was used to estimate light deposition on the cortex [55, 56].

Our tPBM protocol was designed based on previous studies [3, 13–16, 57]. The reason for using the 8 min tPBM period was because it was effective in significantly improving human cognition. The estimated dose density delivered to the human cortex was within the range of positive or photo-stimulatory responses (0.001–10 J cm⁻²) [26, 58].

Participants wore protective goggles throughout the experiment. EEG data were collected using a 64-channel EEG instrument (Biosemi, Netherlands). Each subject wore an EEG cap according to the standard 10–10 EEG electrode placement [59]. Electrode gel was used for keeping the electrode impedances below a pre-set threshold (<10 kΩ). The recorded EEG time series were directed to a computer. We followed the EEG measurement procedures described in [38, 60].

The stimulation protocol (figure 1) consisted of a 2 min baseline (pre), an 8 min stimulation (tPBM or sham), and a 3 min recovery (post) period. The EEG data were acquired at either 256 Hz or 512 Hz; all 512 Hz data were down-sampled to 256 Hz during data pre-processing. tPBM was delivered near electrodes FP2 and AF8, under either sham or active conditions. The study was conducted in a single-blind crossover design, with each subject completing both sham and active tPBM experiments in a random order, with a minimum five day interval between the two experiments.

Zoom In Zoom Out Reset image size

Each EEG dataset represented a 13 min time series of 64 channels during both active and sham tPBM experiments from 45 participants. Because the data processing and analysis included multiple steps in this study, we outline a flow chart in figure 2 to guide the reader through them easily. The five subgroups in data processing are: (a) data pre-processing, (b) PSD-based analysis to obtain frequency-specific power topography, (c) graphical edge formation based on the 'imaginary part of coherence' analysis, (d) GTA-based quantification of global connectivity altered by tPBM, and (e) GTA-based analysis to identify local and nodal graphical metrics changed by tPBM.

Zoom In Zoom Out Reset image size

EEGLAB (an open-source package) was used for data pre-processing. First, EEGLAB's 'filtfilt' function was used to band-pass filter (0.5–70 Hz) the raw EEG data with zero phase distortion, followed by a 60 Hz notch filter to remove line noise. Next, each EEG series was re-referenced by the voltage averaged over all the 64 channels. Robust PCA was then applied to identify and remove significant signal artifacts and outliers from EEG signals [61, 62], followed by independent component analysis [63, 64] to further remove motion artifacts [65, 66], such as eye movements, saccades, and jaw clenching.

To quantify the dose-dependent responses of EEG to tPBM, each artifact-free time series was divided into four temporal sections: (a) the last minute of the 2 min baseline before the onset of stimulation (pre); (b) the first 4 min stimulation period (Stim1); (c) the second 4 min stimulation period (Stim2); and (d) the first 2 min recovery (post).

We used Welch's method to quantify EEG PSD [67]. Specifically, with the use of the 'Pwelch' function (with a 4 sec window and 75% overlap [68]) in EEGLAB, a PSD curve of artifact-free time series for each EEG channel in each time section was calculated. Frequency-specific PSD bandwidths were then selected to cover the delta (1–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), and gamma (30–70 Hz) bands. Next, the mean power change at each of the five frequency bands (f), ΔmPower^f, during each of the three temporal segments (Stim1, Stim2, and post) was normalized to the last minute of its baseline (pre), as expressed [38]:

$$\begin{align} { }\Delta m{\text{Powe}}{{\text{r}}_i}^{\text{f}}& = \frac{{{\text{PS}}{{\text{D}}^{\text{f}}}_i - {\text{PS}}{{\text{D}}^{\text{f}}}_{{\text{pre}}}}}{{{\text{PS}}{{\text{D}}^{\text{f}}}_{{\text{pre}}}}} \times 100\% \nonumber\\ & = \frac{{\left( {{\text{PS}}{{\text{D}}^{\text{f}}}_i - {\text{PS}}{{\text{D}}^{\text{f}}}_{{\text{pre}}}} \right) \times {f_{{\text{band}}}}}}{{{\text{PS}}{{\text{D}}^{\text{f}}}_{{\text{pre}}} \times {f_{{\text{band}}}}}} \times 100\% \nonumber\\ & = \frac{{{\text{Powe}}{{\text{r}}^{\text{f}}}_i - {\text{Powe}}{{\text{r}}^{\text{f}}}_{{\text{pre}}}}}{{{\text{Powe}}{{\text{r}}^{\text{f}}}_{{\text{pre}}}}} \times 100\% { } \end{align} \tag{ 1 }$$

where superscript 'f' denotes the five frequency bands, subscript 'i' represents the three temporal segments (Stim1, Stim2, or post), subscript 'pre' represents the baseline segment, f_band denotes the bandwidth of a chosen frequency band for PSD calculations, and PSD_i and PSD_pre indicate bandwidth-averaged PSD values. Note that ΔmPower is a relative value or percentage change in the bandwidth-averaged power caused by tPBM or sham treatment (see the first two orange boxes in figure 2). To illustrate the difference in ΔmPower between the two conditions, we calculated the sham-subtracted (ss) and tPBM-induced change in power (ΔmPower_ss) at each electrode for each of the five frequency bands within each of the three temporal periods:

$$\begin{align}{ }\Delta m{\text{Power}}_{{\text{ss}},i}^{\text{f}} = \Delta m{\text{Power}}_{{\text{tPBM}},i}^{\text{f}} - \Delta m{\text{Power}}_{{\text{sham}},i}^{\text{f}}.\end{align} \tag{ 2 }$$

Because EEG data at neighbouring time points and spatial channels are highly correlated, it is necessary to perform advanced statistical analysis to remove such correlations and for multi-variable comparisons. For this purpose, we utilized several functions (including 'ft_freqstatistics') available in the FieldTrip toolbox [69, 70] to perform CBPT for statistical comparisons of changes in EEG power (i.e. ΔmPower_ss ^f as shown in equation (2)) among the 64 electrodes in each of the five bands within each of the three temporal periods. In principle, CBPT has two components. The first is the cluster‐forming algorithm, which converts one high-dimensional observation into a quantifiable summary of the cluster structure. The second one creates a surrogate null distribution, against which the observed data are compared to obtain p-values [71].

Accordingly, we first grouped electrodes as clusters within a given scalp distance (e.g. 4–5 cm), followed by the identification of the EEG channels whose ΔmPower_ss ^f values were significantly different from zero for each electrode at a significance level of 0.05. Second, statistical evaluation was performed by taking the sum of the t-values for each cluster. Third, the summed t-value was compared with a null distribution. The null distribution for both permutation tests and cluster-based correction was obtained by randomly permuting the ΔmPower_ss ^f values 1000 times (see the last orange box in figure 2). The corresponding brain regions with modulated powers were identified in this way.

For GTA-based connectivity quantification, we determined the edges or links of a graphical network between all pairs of EEG electrodes. Since correlations between the phases or amplitudes of these EEG channels are interpreted as functional connectivity between these points [41, 51], we performed amplitude and phase decompositions of the time series for all the 64 channels. The amplitude and phase of an EEG time-point can be represented as a complex number [41, 52]. Moreover, we utilized multiple tapers, namely, Slepian sequences, to taper the EEG signal in the time domain before performing the Fourier transform [72, 73]. This part of the calculation was conducted using the 'ft_freqanalysis' function within the FieldTrip toolbox [69].

Coherence, a widely used connectivity measure, is a frequency-domain function equivalent to the time-domain cross-correlation function. The coherence coefficient is a normalized quantity between 0 and 1 and is computed mathematically for the frequency of ω as follows [41]:

$$\begin{equation}co{h_{xy}}\left( \omega \right) = \frac{{\left| {{S_{xy}}\left( \omega \right)} \right|}}{{\sqrt {{S_{xx}}\left( \omega \right){S_{yy}}\left( \omega \right)} }},\end{equation} \tag{ 3 }$$

where S_xx and S_yy denote the power estimates of the signals x and y, respectively, and S_xy represents the averaged cross-spectral density term of the two signals. These terms were calculated using complex values obtained using the multitaper method.

Because of volume conduction in the human brain, signals generated in one region can be detected by several electrodes, resulting in an artificially high coherence value among these channels. To overcome this effect, the magnitude of the operation can be removed from equation (3), the imaginary part of S_xy is considered, and the cross-spectral density of the signals with a phase difference of 0 or 2π is set to zero. This method is called the 'imaginary part of coherence', which explicitly removes instantaneous interactions [41]. Studies have shown that this method is excellent to minimize the volume conduction issue in EEG data analysis [74].

The pairwise connectivity values for all pairs of electrodes (64 in this study) can be represented by an n × n (i.e. 64 × 64) adjacency matrix, where n is the number of nodes (i.e. 64 channels). The FieldTrip toolbox facilitates the computation of the imaginary part of coherence for all pairs of channels using the 'ft_connectivityanalysis' (the first two green boxes in figure 2).

In this study, each temporal segment was divided into 10 sec epochs, and the adjacency matrices generated for all epochs in each frequency band were averaged for each of the three temporal segments and five frequency bands for each subject individually. These averaged matrices were binarized by varying the sparsity level and used for GTA-based global and nodal connectivity analyses, as described below.

The GTA enables the exploration of topological changes in brain networks through pairwise functional connectivity between channels. A network can be characterized based on functional segregation, functional integration, and centrality. Previous studies showed tPBM-induced alterations in graph measures of the brain network with different setups and protocols [36, 40]. However, these studies focused on residual modulation in post stimulation with 19 channels.

In this study, we used GRETNA [75], a widely used GTA toolbox, to quantify the global and nodal graphical metrics of the human brain network for individual subjects under active and sham tPBM in four temporal segments and five frequency bands. This step was repeated 19 times to assess the chosen metrics under a sparsity range of 5%–95% with an increment of 5% (see the last two green boxes in figure 2).

To examine tPBM-induced effects, three global graphical measures were chosen for analysis: synchronization (S), global efficiency (GE), and small-worldness (SW). The group-level values for each global measure at each sparsity level were statistically compared between the active tPBM and sham conditions using a paired t-test (pink boxes in figure 2).

Two previous publications showed alterations only in global network metrics [40] and within each hemisphere or between two hemispheres [36]. In this study, our analysis quantified five nodal graphical metrics, namely, nodal clustering coefficient (nCC), nodal efficiency (nE), nodal local efficiency (nLE), betweenness centrality (BC), and degree centrality (DC). The respective definitions of graph metrics are listed in supplementary material A.

Although GTA was performed on the 64-channel EEG, resulting in tPBM-induced changes in nodal network metrics, the 64 nodal locations were too dispersed to identify the cortical regions on the human scalp. Thus, we focused on ten local sections according to prefrontal, central, temporal, parietal, and occipital regions [76]. Thus, 64 nodes were grouped into ten clusters with 6–10 electrodes in each cluster (figure 3).

Zoom In Zoom Out Reset image size

At the subject level, each nodal graphical metric within each cluster area was obtained by averaging the specific metric over all the electrodes within the respective region (for each of the three temporal segments and five frequency bands). To compare the changes induced by tPBM and sham, nodal measures for each temporal segment were baseline-normalized by subtracting the corresponding baseline (pre) values from those in each of the three time windows (Stim1, Stim2, and post). Next, for each cluster region, group-level (n = 45) and baseline-subtracted nodal metric values were compared between the active and sham conditions using paired t tests. To correct for multiple comparisons, false discovery rate (FDR) correction was performed for ten regions with a corrected significance level of 0.05 (see the two vertical gray boxes in figure 2).

As described in sections 2.5 and 2.6, the baseline-normalized values of ΔmPower^f (see equation (1)) for each group of tPBM and sham conditions among the three temporal segments (Stim1, Stim2, and post) and in five frequency bands were calculated. To demonstrate clear statistical differences in ΔmPower^f between the two stimulation conditions, baseline-normalized and sham-subtracted topographical maps of ΔmPower_ss ^f (%) values (see equation (2)) over 64 channels were achieved, as shown in figure 4, for all five frequency bands. In addition, after CBPT for 64-channel statistical comparison, the electrode sites/clusters that were significantly affected by tPBM are superimposed on the topographies in figure 4 with '*' denoting p < 0.01 and '×' denoting p < 0.05.

Zoom In Zoom Out Reset image size

These results illustrate a significant, dose-dependent increase in EEG rhythm powers at 8–13 Hz and 13–30 Hz during the last 4 min of tPBM (Stim2). Specifically, the increase in alpha ΔmPower_ss was seen as two major clusters of channels in the bilateral frontal and left parietal-occipital regions, whereas the increase in beta ΔmPower_ss was mainly seen as one cluster of electrodes in the central/parietal region of the scalp. The enhanced alpha ΔmPower_ss remained in the affected locations during the post-tPBM period, whereas the significant increase in beta ΔmPower_ss ceased during the recovery time. Furthermore, delta power was reduced in the frontal, left temporal, and occipital regions during tPBM, and in the right frontal region during recovery.

Following the steps given in sections 2.7–2.9, adjacency matrices for all three temporal segments and five frequency bands were generated. These matrices were further binarized for different sparsity values, resulting in the GTA-derived graphical networks. In this study, we identified three global network metrics (S, GE, and SW) that were significantly altered by tPBM with respect to the sham condition and only in the beta band. As shown in figure 5, the three rows illustrate the respective global metrics during Stim1, Stim2, and the post period under both active and sham stimulation.

Zoom In Zoom Out Reset image size

These results clearly show that tPBM significantly reduces the global synchronization, GE, and SW of the network connectivity of the human brain. Specifically, significant decreases in synchronization and GE occurred during Stim1 and the recovery period with more sparsity values, while a significant reduction in SW appeared in Stim2 with more sparsity units. We also confirmed that there was no significant difference between the pre-stimulation baselines under tPBM and sham conditions for any of the three global network metrics.

After performing the analysis steps given in sections 2.9 and 2.10, cluster-averaged, baseline-subtracted values for each of the five nodal graphical metrics (i.e. nCC, nE, nLE, BC, and DC) were obtained for each of the ten spatial clusters under both tPBM and sham conditions. After performing paired t-tests with FDR correction for ten spatial clusters (i.e. p < 0.05, FDR corrected), we identified the clusters whose nodal metric values were significantly altered by tPBM for each of the five metrics at all three temporal periods and only in the beta rhythm band. Topographical representations of the results for the five nodal metrics are shown in figure 6.

Zoom In Zoom Out Reset image size

Based on figure 6, we made the following observations. (a) During Stim1, significant increases in nE, BC, and DC were observed in the right frontal region near the tPBM stimulation site. (b) During Stim2 and post stimulation, significant changes occurred in the bilateral frontal regions for all five nodal metrics. More specifically, tPBM significantly decreased the clustering coefficient and nLE, whereas stimulation significantly increased the other three nodal metrics. (c) Combined temporal and spatial results revealed that both nE and DC were initially enhanced by tPBM in the right frontal region during Stim1, followed by expansion of this enhancement to the contralateral side during Stim2, which persisted during the post-tPBM period. (d) In the case of BC, unilateral enhancement in the right frontal region remained during the entire stimulation time (Stim1 and Stim2) and then expanded to the contralateral side in the post. (e) On the other hand, during Stim2 and post-stimulation, significant decreases occurred in the bilateral frontal regions for the nCC and nLE. (f) During the same time periods for the same two metrics (nCC and nLE), decreases were observed in the left parietal and occipital regions. (g) Moreover, the left temporal region showed a reduction in nE and BC only for Stim2. (h) The only significant modulation in the right occipital region was a decrease in BC in Stim2.

As shown in figure 4, under the eyes-closed resting-state condition, tPBM significantly increased the power of alpha oscillations in clusters over the bilateral frontal, left parietal, and left occipital regions, as well as the beta power over the bilateral central and parietal regions during the second 4 min of stimulation. These observations agree with previous reports on eyes-open tPBM experiments [60, 77].

A significant increase in alpha power over the frontal-parietal regions clearly confirmed the ability of tPBM to neuromodulate the frontoparietal network, which is an executive network facilitating rapid instantiation of new tasks [78]. It is reported that alpha rhythm is associated with awareness [79] and cognitive functions, such as memory encoding and attention [80–82]. The same experimental protocol with a 1064 nm laser was used previously and demonstrated a significant improvement in human cognitive performance [3, 13–15, 57]. Moreover, the presence of stronger beta waves has been linked to better cognitive ability, as reported in several studies [83, 84]. Thus, the beneficial outcome of frontal tPBM in human cognition can be, at least partially, attributed to its significant modulation of electrophysiological alpha and beta powers in the frontoparietal network. In addition, we attributed the reduction in delta power in the channel clusters during tPBM to the sensation of laser-induced warmness in the forehead skin [38].

As shown in figure 5, tPBM significantly changed three global graph measures, namely, synchronization, GE, and SW, in the beta wave only. We discuss each of these changes as follows.

Our observation that brain network synchronization in the beta band was significantly reduced during Stim1 and recovery agrees with a recent study that reported the effects of tPBM on brain network synchronization with 850 nm LEDs [36]. In addition, a behavioral study attributed a decrease in synchronization in the healthy human brain to awareness and cognitive processing [85].

Similarly, the GE was significantly reduced by tPBM compared with sham during Stim1 and post. This reduction indicates a decrease in brain network integration, which may imply less efficient or more complex information paths in the network. In other words, tPBM may increase the energy and wiring costs of information flow owing to the trade-off between network efficiency, energy, and wiring costs [86, 87]. This could be an indicator of increased brain complexity related to higher cognitive function [88]. This observation supports the expected benefit of tPBM, namely, the beneficial effects of tPBM on cognitive improvement.

SW is calculated as the ratio of normalized integration to normalized segregation [50]. Thus, the reduction in this metric could be attributed to a significant reduction in global integration (as reflected by a reduction in GE) and/or a significant increase in global segregation of the brain network caused by tPBM. These observed significant effects of tPBM on SW taking place only in Stim2 could result from the resistance of resting-state networks against changes in network composition as well as the dose-dependent nature of tPBM-induced effects on neural activity [89, 90]. However, a possible explanation for the lack of significant alterations in synchronization and GE in Stim2 could be the high variability in the functional topography of the frontoparietal network [78].

In addition, the observation that tPBM altered only the beta-wave oscillation in the EEG graphical network was in agreement with other studies [36, 49]. Several studies have shown the role of the beta band in different brain networks [91] and cognitive functions [92, 93].

Significant reductions in nCC and nLE in bilateral frontal regions and in the left parietal and left occipital regions during Stim2 and post implied that the clusters of nodes in these regions of the network became less segregated during the 2nd 4 min and post period of tPBM. In other words, tPBM facilitated less separation and more integration of nodal graphical connectivity. Consistently, nE (reflecting nodal integration) increased in the right frontal region in the beginning of tPBM. This increase indicates enhanced integrity of the nodes over this region in the information flow [94] and parallel information transfer. Furthermore, during the last 4 min of tPBM, the bilateral frontal regions showed an increase in integration into the functional network, while the integrity of the left temporal region in the network was significantly reduced. This phenomenon implies that tPBM stimulated more network integration at the beta rhythm in the frontal regions, with the cost of reducing network integration in the left temporal segment/cluster.

BC represents the fraction of all shortest paths in the network that pass through a particular central node [95]. A large value of BC means a large impact of this central node on information flow over the network. Also, DC quantifies the number of links from nodes in a specific region to other nodes in the same or other regions [50]. Increases in both BC and DC in the frontal region during and post tPBM suggested that the frontal regions became more prominent in connecting disparate parts of the network. The enhancement of BC and DC revealed that these cortical regions could be prominently stimulated by tPBM for more information connections at beta oscillations. It has also been reported that a working memory task during the encoding phase triggers similar increases in DC over the frontal regions of the beta band [93].

Combining all these observations, we concluded that tPBM facilitated a reduction in local segregation, increases in nodal integration and centrality of frontal regions, and the growth of connection links between nodes in these frontal regions compared to other regions. These results agree with the observed changes in the flow of information reported in a tPBM-evoked causal connectivity study [96].

The alpha and beta powers of the human brain, especially in the frontoparietal network, are believed to be related to cognitive functions, such as memory encoding and attention [80–82]. Our observations clearly demonstrated that tPBM increased the alpha and beta powers in the frontal-central-parietal regions, indicating the underlying association between tPBM and the enhancement of human memory.

However, our results showed that tPBM altered EEG graphical network metrics only in the beta band, which was consistent with the results given in [36]. According to [91], beta oscillations in the prefrontal region appear to serve as short-term memory executors and focus enhancers during executable tasks. In addition, the beta rhythm in the temporal lobe plays an important role in long-term memory retrieval [97]. Memory retrieval starts in the temporal lobe, passes through different parts of the neocortex, and stops in the prefrontal cortex [91]. There is always a balance in the information flow in this order, even during the resting state. We speculated that tPBM enables neuromodulation of beta oscillations and the corresponding network connectivity globally across the scalp and regionally in several nodal regions. The significant modulation of beta-wave connectivity in the human brain may be an underlying association between tPBM and the enhancement of human cognition. This speculation may be validated through PET imaging using fluorodeoxyglucose [98].

As shown in table 1, only two recent studies reported tPBM-induced modulations of global network metrics [40] and alterations in brain connectivity between the two hemispheres [36]. It would be helpful to compare the results from the two studies with those found in this study. Detailed summaries and comparisons are given in supplementary material B.

There is a relatively unresolved debate regarding the basic mechanism of action of tPBM. One side of the debate states that a red or NIR photon can photo-dissociate nitric oxide (NO), a molecule that may inhibit CCO by noncovalent binding [99, 100]. This view emphasizes the greater effects of tPBM on diseased/damaged cells because unhealthy cells are more likely to have inhibitory concentrations of NO. On the other hand, the other opinion highlights tPBM-induced effects on CCO, which stimulates ATP production by increasing mitochondrial membrane potential and oxygen consumption [26, 58]. This view enables to explain the tPBM-induced augmentation of human cognition in healthy humans besides in patients with brain disorders [26, 58, 101]. This side of opinions has been supported by several publications [58], and now by two other independent studies that also reported significant improvement of working memory in healthy young [102] and older adults [103] by repeated tPBM. However, the two opinions on the basic mechanism of tPBM have coexisted for years, without one outweighing the other. Further confirmation of each remains to be explored in future studies.

First, the international 10–10 electrode placement system in this study was not strictly followed on the human head because a clear area with a 4.2 cm diameter was needed for tPBM light delivery. The EEG cap was shifted 1–2 cm backwards. Second, our power spectral and connectivity analyses were performed in the sensor space. Source space analysis can be conducted to observe cortical and subcortical regions in the brain that are affected by tPBM. Third, the current study was based on EEG signals of the tPBM-treated human brain in the resting state. Concurrent assessments of changes in cognitive enhancement and brain connectivity after tPBM would provide quantitative correlation and association between functional connectivity and behavioral effects of tPBM.

As for future work, there is a lack of systematic studies on the optimization parameters of tPBM, such as the wavelengths used, light irradiance, stimulation dose, and sites, as well as the effective period after stimulation. In particular, tPBM has a hormetic dose response, which is characterized by stimulation of a biological process at a low dose and inhibition of that process at a high dose [26, 58]. Namely, photo-stimulatory or photo-inhibitory effects occur with low (0.001–10 J cm⁻²) and high (>10 J cm⁻²) optical energy (dose) density [104–106]. However, the documented dose range for positive responses results from in vitro experiments and is quite broad (0.001–10 J cm⁻²). It is difficult to select an optimal irradiance for non-invasive human treatments with tPBM. Accordingly, more research is needed with both healthy humans and selected populations of patients before tPBM becomes an effective device to treat patients with brain disorders and for healthy aging in the rapidly growing aging population.