Transient brain-wide coactivations and structured transitions revealed in hemodynamic imaging data.

Brain-wide patterns in resting human brains, as either structured functional connectivity (FC) or recurring brain states, have been widely studied in the neuroimaging literature. In particular, resting-state FCs estimated over windowed timeframe neuroimaging data from sub-minutes to minutes using correlation or blind source separation techniques have reported many brain-wide patterns of significant behavioral and disease correlates. The present pilot study utilized a novel whole-head cap-based high-density diffuse optical tomography (DOT) technology, together with data-driven analysis methods, to investigate recurring transient brain-wide patterns in spontaneous fluctuations of hemodynamic signals at the resolution of single timeframes from thirteen healthy adults in resting conditions. Our results report that a small number, i.e., six, of brain-wide coactivation patterns (CAPs) describe major spatiotemporal dynamics of spontaneous hemodynamic signals recorded by DOT. These CAPs represent recurring brain states, showing spatial topographies of hemispheric symmetry, and exhibit highly anticorrelated pairs. Moreover, a structured transition pattern among the six brain states is identified, where two CAPs with anterior-posterior spatial patterns are significantly involved in transitions among all brain states. Our results further elucidate two brain states of global positive and negative patterns, indicating transient neuronal coactivations and co-deactivations, respectively, over the entire cortex. We demonstrate that these two brain states are responsible for the generation of a subset of peaks and troughs in global signals (GS), supporting the recent reports on neuronal relevance of hemodynamic GS. Collectively, our results suggest that transient neuronal events (i.e., CAPs), global brain activity, and brain-wide structured transitions co-exist in humans and these phenomena are closely related, which extend the observations of similar neuronal events recently reported in animal hemodynamic data. Future studies on the quantitative relationship among these transient events and their relationships to windowed FCs along with larger sample size are needed to understand their changes with behaviors and diseased conditions.

Introduction

Spontaneous fluctuations in neural signals of both electrical and hemodynamic forms are characterized with brain-wide functional connectivity (FC) in mammalian brains at rest, known as resting-state networks, RSNs (Fox and Raichle, 2007). These networks have been reported in data acquired by various imaging modalities, including functional magnetic resonance imaging (fMRI) (Beckmann et al., 2005, Biswal et al., 1995), electroencephalography (EEG) (Yuan et al., 2016), magnetoencephalography (MEG) (Brookes et al., 2011), and optical imaging (Eggebrecht et al., 2014). It has been extensively reported that these resting-state FCs exhibit consistencies across participants (Damoiseaux et al., 2006, Yeo et al., 2011), and show significant behavioral correlates (Allen et al., 2014, Hutchison et al., 2013b, Rosenberg et al., 2016, Kelly et al., 2008) as well as clinical abnormalities in major neuropsychiatric disorders (Fox and Greicius, 2010).

Resting-state FCs are observed primarily as patterns of covariations between brain regions calculated with either Pearson’s correlation (Biswal et al., 1995) or blind source separation techniques, e.g., spatial independent component analysis (ICA) (Beckmann et al., 2005). However, the exact neurophysiological mechanisms underlying these correlations remain elusive. In fMRI, increasing evidence have been documented that such resting-state FCs are nonstationary over a scan (Preti et al., 2017, Allen et al., 2014, Chang and Glover, 2010). In contrast to static FCs over the entire scanning window, estimations of FCs over sequential sliding-time windows (either partially overlapped or non-overlapped) within a scan have led to the so-called dynamic FCs (dFC). Similarly, dFCs have been observed in other brain signals, e.g., EEG (Shou et al., 2020). It has been further suggested that dFCs potentially have higher sensitivities than static FCs for correlating with behavioral and disease conditions (Damaraju et al., 2014, Hutchison et al., 2013a). The fact that FCs are time-varying has motivated numerous efforts of animal and human studies to identify potential neurophysiological events responsible for correlation-based FCs observed in various brain signals, leading to several candidate hypotheses.

Brain-wide propagating activity, characterized by wave-like recurring patterns (at mesoscopic and macroscopic levels), is one of these hypotheses behind FCs (Muller et al., 2018). Propagating wave dynamics of large-scale spontaneous neural networks have been reported in animals (Huang et al., 2010, Stroh et al., 2013, Steriade et al., 1993, Luczak et al., 2007) and humans (Halgren et al., 2019, Gu et al., 2021, Raut et al., 2021, Massimini et al., 2004). For example, globally-distributed cortical waves have been recorded using optical imaging in mice at a timescale of seconds (Matsui et al., 2016). Similarly, lateral-to-medial propagating waves in rats (Majeed et al., 2011, Majeed et al., 2009) and global wave patterns in humans (Raut et al., 2021) have been reported in blood-oxygen-level-dependent (BOLD) signals. In humans, more types of propagating activities have been observed, including alternating activation waves between the default mode network (DMN) and task-positive networks (TPNs) (Majeed et al., 2011, Yousefi et al., 2018) and directionally constrained waves along the spatial axis representing the cortical hierarchical organization (Gu et al., 2021). These waves may be the underpinnings of correlational structures in resting-state FCs (Matsui et al., 2019), or at least critical in shaping the dynamics in resting-state FCs (Gu et al., 2021), and potentially with novel information about brain functional organizations in parallel to correlational FCs (Majeed et al., 2011).

The second hypothesis proposes that resting-state FCs are driven by short-lived coactivation patterns (CAPs) (Liu and Duyn, 2013). These CAPs are obtained via a point process analysis focusing on selected fMRI frame data surrounding regional peaks (Tagliazucchi et al., 2011, Cifre et al., 2020). Strong correspondences between these transient events and resting-state FCs are reflected in their spatial patterns (Liu et al., 2013), and CAPs account for significant variances describing resting-state FCs (Tagliazucchi et al., 2016). Recent clustering studies on the entire fMRI frame data further indicate that similar coactivation patterns exist beyond local peaks of fMRI activity in predefined regions-of-interests, leading to the detection of the so-called manifold recurring brain-wide CAPs (Gutierrez-Barragan et al., 2019), while their relationship to resting-state FCs is yet to be established. Moreover, CAPs have been recently reported in neuronal electrical signals in both animals (Matsui et al., 2019) and humans (Ding et al., 2021). An obvious advantage of CAP-based analysis is to allow interrogation of dFCs at fine temporal scales down to single timeframes (Liu et al., 2018b), eliminating the need for a priori selected sliding windows as typically required in correlational methods. Moreover, CAP studies have shown clinical relevance (Liu et al., 2018b, Yang et al., 2021, Marshall et al., 2020, Rey et al., 2021, Kaiser et al., 2019, Zhuang et al., 2018).

The third hypothesis is global brain activity, which is linked to a well-known phenomenon in fMRI data analysis that is currently still under debate (Liu et al., 2017), i.e., global signal (GS). The anticorrelated DMN and TPNs, for example, have shown strongly positive correlation under the influence of GS, resulting in global brain activity patterns (Yousefi et al., 2018). Furthermore, global brain activity has been observed as transient brain states during brain-wide propagating waves (Matsui et al., 2019). GS is defined as the time series of intensity averaged across imaging voxels in fMRI, which is usually regressed out in the preprocessing of fMRI data (Desjardins et al., 2001, Aguirre et al., 1998), as it is believed that physiological noises substantially contribute to it (Liu et al., 2017). However, recent studies have demonstrated that GS contains neuronally relevant information (Murphy and Fox, 2017, Schölvinck et al., 2010) and is primarily modulated by the fluctuating state of arousal (Wong et al., 2013, Liu et al., 2018a). While GS is known for affecting FC analysis (Murphy et al., 2009, Murphy and Fox, 2017, Fox et al., 2009), it also drives dFC patterns (Chang et al., 2009). All these facts support the potential role of global brain activity in generating resting-state FCs.

It is important to note that these three phenomena (i.e., propagating waves, CAPs, and global brain activity) have shown linkages among each other. As discussed above, propagating waves are of global nature (Stroh et al., 2013). Manifold recurring brain-wide CAPs have been suggested to be phase-locked to GS in rodents (Gutierrez-Barragan et al., 2019, Ma et al., 2020). Furthermore, a recent mice study (Matsui et al., 2016) using a wide-field optical imaging concurrently monitoring calcium and hemodynamic signals indicates all three phenomena in both neuronal spiking and hemodynamic data in which transient CAPs are embedded in globally propagating waves. While these three phenomena have been studied separately in both animals and humans, studies indicating their potential connections have been largely performed in animal models (Gutierrez-Barragan et al., 2019, Matsui et al., 2016, Ma et al., 2020) and the understanding of their interactions in humans remains sparse. To further investigate these structured neurophysiological phenomena that are possibly responsible for the formation of correlation-based FCs in the human brain, the present study utilized a recently proposed brain-wide diffuse optical tomography (BW-DOT) technique (Khan et al., 2021b) to map spatiotemporal structures in hemodynamic signals recorded with a cap-based whole-head continuous-wave functional near-infrared spectroscopy (fNIRS) system (Zhang et al., 2021).

DOT has a higher sampling rate than fMRI, which exempts DOT hemodynamic measurements from aliasing noises caused by breathing and systematic fluctuations in circulation (Lowe et al., 1998, Tong et al., 2011). Moreover, DOT technology is relatively robust to head motions (Aslin and Mehler, 2005), which have shown considerable impacts on FCs estimated from BOLD fMRI (Laumann et al., 2017). Advanced from the previous patch-based DOT systems in reconstructing regional RSNs and/or regional nodes of distributed RSNs (Wheelock et al., 2019), our recent study has demonstrated that this new brain-wide DOT, with its whole-head optode placements providing a field-of-view (FOV) covering most of the cortex (i.e., 109 long-separation channels), can reconstruct a collective set of brain-wide RSNs (Khan et al., 2021b). The whole-head coverage provides an ability to study brain-wide transient events, e.g., CAPs and global propagating waves, and a means for estimating GS that is not available in patch-based DOT systems. We reasoned that these capacities of BW-DOT, together with other denoising approaches (e.g., 8 short-separation channels), would also improve our ability in characterizing temporal dynamics of these neuronal events and better distinguishing those of neuronal or non-neuronal origins.

The present pilot study investigated the spatiotemporal structures in DOT data from thirteen healthy participants. We firstly estimated cortical RSNs based on deoxygenated hemoglobin (HbR), using participant-specific MRI head models and state-of-the-art algorithms for data preprocessing and reconstructions. We then adopted a framewise clustering approach on participant-level spatiotemporal RSN data to obtain CAPs with and without the GS regression. We further analyzed the spatial, temporal, and transitional properties of detected CAPs. Lastly, we investigated the differences of global CAPs (indicating global activity) obtained with and without the GS regression regarding their plausible neural relevance and their relationship to GS time series.

Materials and methods

Participants and experimental protocols

The study was approved by IRB at the University of Oklahoma, and written informed consent was obtained from every participant before the study. A total of 20 participants were recruited for the study without any neurological or neuropsychiatric disorders. Three participants were not able to complete all data recordings. In addition, resting-state data from four participants were excluded from the present study due to poor data quality (primarily excessive motions). Thereafter, resting-state data from thirteen participants (31.7 ± 9.3 years, five females, eight males) were included in the present study. Other analyses have been conducted on the same resting-state data and reported (Khan et al., 2021b), as well as on motor task data performing hand-clenching that were collected with resting data in the same protocol (see descriptions below) and reported (Zhang et al., 2021), which excluded data from additional 3 participants due to excessive motions only during tasks.

Structural head MRI of each participant was acquired at OU Health Sciences Center, using a GE Discovery MR750 whole-body 3-Tesla MRI scanner (GE Healthcare, Milwaukee, WI, USA). The scanning parameters were: FOV = 240 mm, axial slices per slab = 180, slice thickness = 1 mm, image matrix = 256 × 256, TR/TE = 8.45/3.24 ms. Data were acquired from two sessions on the same day, each consisting of a block of resting state (6 min) and a block of a motor task of hand-clenching (5.5 min, these data were not analyzed in the present study). In each resting-state block, the participants sat upright on a comfortable chair in a dark and sound-damped room and were instructed to rest with their eyes open and fixated on a cross mark shown on a computer screen and to stay as still as they could to avoid motion-induced artifacts. E-Prime (Psychology Software Tools, Sharpsburg, PA) software was used to program visual cues and instructions displayed on an LCD monitor located about 90 cm from the participant. Data from fNIRS optodes, auxiliary sensors, and EEG electrodes were captured simultaneously in all sessions, but EEG recordings were not analyzed in the present study. fNIRS data were gathered at a sampling rate of 6.25 Hz using a dual-wavelength (760/850 nm) continuous-wave fNIRS device (NIRx GmbH, Germany). The fNIRS system included 34 LED sources, 31 avalanche photodetectors, and 8 short-separation (SS) detectors (inter-optode distance, IOD: 8 mm), resulting in 109 long-separation (LS) and 8 SS channels. The optodes were attached to an elastic cap (EASYCAP GmbH, Herrsching, Germany) designed with the 10–5 system for electrode/optode placements (Oostenveld and Praamstra, 2001). To maximize FOV while retaining the high-density coverage, the optodes were uniformly placed over the entire scalp (IOD mean±std: 3.65 ± 0.44 cm, min: 2.98 cm, max: 4.56 cm). The pneumatic respiratory belt, a triaxial accelerometer, and a photoplethysmography sensor were worn on the belly, forehead, and left index finger, respectively, to capture chest movements, head motions, and cardiac beats at a sampling rate of 500 Hz using an EEG amplifier (Brain Vision, NC, USA). The Polhemus Patriot handheld electromagnetic digitizer (polhemus.com) was used to record the 3D locations of optodes and head landmarks (i.e., nasion, left, and right preauricular points) before recording sessions.

Preprocessing of fNIRS data

The preprocessing of fNIRS data followed the protocol described in our previous study (Zhang et al., 2021) using the HOMER2 toolbox (Huppert et al., 2009). Data from each session were preprocessed separately. Briefly, raw channel-wise light intensity data were transformed to optical densities (OD) by measuring their log-ratios against the channel-wise means over the entire recording window as baselines (see equation in Section 2.3). Power spectral densities (PSDs) of OD data were calculated using the Welch’s method (window length: 60 s, 50% overlapping). The LS-channel OD data showing no visible heartbeat power peaks in the frequency range of 0.8–1.6 Hz were discarded (indicating possible inadequate optode-scalp contact). As the sampling frequency of DOT (6.25 Hz) was suffiently high, resting OD data were then bandpass filtered at 0.008–0.2 Hz to retain useful information in the higher frequency band up to 0.2 Hz (Smith et al., 2012). Bad segments corresponding to excessive head motions were identified based on the metric of global variance in temporal derivative (GVTD) (Sherafati et al., 2020). Specifically, the time points with GVTD > 3 × the mean GVTD over all time points were labeled, and the recordings of 10 s centered at these time points were excluded from further analyses. After this step, 73.4% ± 7.5 of data in participants were retained. To remove the physiolgocial noise related to superficial tissue absorption, motion and pulse/respiration related variations (Zhang et al., 2019, Zhang et al., 2021, Khan et al., 2021b, Cheong et al., 2020), we computed nuissances regressors based on time courses of 8 SS channels that were evenly distributed across the scalp and removed their contributions to LS-channel OD data using general linear model regression (referred to as the SS regression).

Diffuse optical tomography

DOT forward modeling was conducted using our previously developed framework on participant-specific head models with finite element (FE) volume meshes (tetrahedrons: 866,929 ± 223,964, nodes: 140,774 ± 35,883, Fig. 1C) (Khan et al., 2021b). A linear forward relationship between measured light changes at the scalp nodes due to small perturbations in chromophore concentrations is given by (Srinivasan et al., 2005, White, 2012): where are measurement vectors of OD at two light wavelengths λ1 (760 nm) and λ2 (850 nm), respectively. The OD was defined as the log-ratio of measured light intensity I(t) with reference to the baseline light intensity I, e.g., . The unknown variable denotes vectors of time-varying HbR and oxygenated hemoglobin (HbO) concentration changes at the source nodes within the source space. Because the sensitivity of fNIRS drops considerably with the depth, the source space is made up of FE mesh nodes inside the brain (contained by the pial surface) within 45 mm of the scalp (Dehghani et al., 2009b). The forward operator A (Eq. (1)), known as the sensitivity matrix or the Jacobian matrix, was calculated using the NIRFAST software (Dehghani et al., 2009a) with realistic baseline optical properties of brain tissues (Eggebrecht et al., 2012) and smoothed to remove point-like artifacts from numerical computations (Khan et al., 2021b).

Fig. 1.

The framework for reconstructing brain-wide CAPs from individual timeframe data of DOT. Diffuse optical tomography (DOT) consists of (A) the cap-based high-density optode placement system; (B) continuous-wave dual-wavelength (760 mm and 850 mm) light intensity fNIRS measurements; (C) realistic finite-element head models segmented from individual participant MRI data; and (D) inversely estimated cortical distribution of HbR concentration changes. (E) group-level Spatial ICA (gSICA) on temporally concatenated DOT data (D) across N participants. Spatiotemporal regression (STR) performs (F) spatial regression to estimate participant-level time courses of ICs (group-level IC spatial maps as regressors and DOT data as responses); and (G) temporal regression to estimate participant-level spatial maps of ICs (participant-level IC time courses as regressors and DOT data as responses). CAPs: (H) participant-level cluster centroids (k × number of ICs) are obtained via the k-means++ clustering analysis on temporally concatenated IC time courses. The clustering indices are used to average DOT timeframe data to obtain (I) participant-level spatial maps of CAPs; and (J) group-level spatial maps of CAPs.

The inverse solution of Eq. (1) was obtained via a minimum-norm estimate (MNE) (Hamalainen, 1984) with the Tikhonov regularization (Tikhonov, 1963): where the first term denotes the residual errors in fitting measurement data, is the L2-norm regularization term (Michel et al., 2004). The λ(t) is the regularization parameter, which was estimated at each time point using the L-curve method (Hansen and O’Leary, 1993) implemented in the Regularization Tools toolbox (Hansen, 2007) and then smoothed to remove potentially unreasonable values as a result of numerical calculation of the L-curve (Khan et al., 2021b). The minimization of Eq. (2) yielded the inverse solution: where I is the identity matrix and the inverse solution is a vector of joint estimates of HbR and HbO concentration changes, and the subscript ν represents the volumetric solution space. These volumetric data were spatially smoothed using a 6-mm spherical kernel to reject high-frequency spatial noise (Friston et al., 1995). Only HbR data, close to BOLD fMRI (Buxton, 2013), were further analyzed in the following steps.

Estimation of brain-wide resting-state networks

Inverse volumetric data of participants obtained via DOT were firstly mapped to their respective cortical surfaces (Fig. 1D). These participant-specific cortical surfaces (20,484 nodes) were defined as the gray and white matter interface segmented from MRI data using FreeSurfer (Fischl, 2012, Ségonne et al., 2004). The volume-to-surface projection was used as hemoglobin concentration changes close to the gray matter are believed to reflect more neuronal activity (Obrig and Villringer, 2003, Scholkmann et al., 2014). This is different from DOT reconstructions (see Section 2.3) where light propagation involves both cerebral and extracerebral tissues and therefore participant-specific volumetric models are used. The projected cortical surface DOT data from participants (two sessions per participant) were concatenated and subjected to a group-level spatial ICA (Calhoun et al., 2009) (Fig. 1E). The number of ICs was empirically set to 20 (Khan et al., 2021b). The participant-level IC spatial maps (Fig. 1G) and time courses (Fig. 1F) were obtained via the spatiotemporal regression, STR (GIFT, 2020, Beckmann et al., 2009) for preparing the clustering analysis (Section 2.5). This STR approach (also named as dual regression or multi-level regression) has been widely used in functional neuroimaging data (Erhardt et al., 2011). Group-level spatial maps of ICs were then obtained by averaging corresponding IC spatial maps from all participants. These spatial maps, termed DOT RSNs, have shown similarities to fMRI RSNs in our recent study (Khan et al., 2021a, Khan et al., 2021b).

Estimations of CAPs and their spatial, temporal, and transitional properties

Time courses of 20 DOT RSNs (i.e., ICs) from all participants were concatenated and subjected to a k-means clustering analysis to find distinct CAPs with DOT RSNs as building blocks. To remove scale discrepancies between participants and ICs, the time course of each IC was normalized in each participant to obtain z-scores. After that, all participant-level z-score data for corresponding ICs were temporally concatenated (Fig. 1F). The k-means++ clustering using the Euclidean distance metric (Arthur and Vassilvitskii, 2006) was applied to these concatenated z-scored IC timeframe data (20 ICs × # of time points from all participants) on the time dimension to yield k clusters. The IC timeframes (each is a vector of 20 elements) from the same participant that were classified into the same cluster were firstly averaged to generate the participant-level centroids (Fig. 1H), which represented the geometric center of the subspace expanded by all timeframe vectors belonging to the same cluster in a 20-dimensional space. Next, the participant-level spatial map of a CAP (Fig. 1I) was obtained as an average of the DOT timeframe data belonging to the CAP from clustering. The group-level spatial maps of CAPs (Fig. 1J) were obtained via averaging participant-level CAPs across all participants. Note that the participant-level CAPs were obtained in order to perform statistical analyses on spatial, temporal, and transitional metrics discussed in this section to characterize CAPs using participants as samples. To evaluate the spatial link among all CAPs at the group level, the distances between their cluster centroid vectors were calculated and projected onto a 3D space using a multidimensional scaling tool from MATLAB, i.e., cmdscale.m (Seber, 2009). Clustering was repeated with different k values (6≤ k ≤30) (see most results in Supplementary Fig. 2A). Percentage-wise explained variance and explained variance gain were calculated with the method adopted for neural data (Shilling-Scrivo et al., 2021) using the k-means_opt function in MATLAB File Exchange (Landtsheer, 2018). Their curves were plotted to help choose an optimal cluster size, i.e., k (Supplementary Fig. 2B). With an elbow area in the range of k=4–7, the explained variance gain between consecutive k values showed a decreasing pattern and was consistently less than 2% beyond k=7. We, therefore, reported results from k=6 as they showed the highest stability (after performing the clustering several times using different seed points) and explained 87.1% of the variance.

Following clustering, entire recording data in each participant were characterized as a sequence of recurring CAP events. To examine dynamic patterns in these sequences, multiple temporal, dynamic, and spectral metrics of CAPs were computed on each participant and then summarized to generate group-level statistics. An occurrence of a CAP was defined as a set of continuous timeframes belonging to the same CAP. The occurrence rate of a CAP was calculated as the ratio of its total occurrences to the total occurrences of all CAPs. The mean dwell time of one type of CAP was calculated as the mean time duration from all its occurrences. The mean interval time of one type of CAP was calculated as the mean time duration between two neighbored occurrences of the CAP. To probe the relationship between a pair of CAPs, the one-step transition probability of one type of CAP to another type of CAP (directional) was calculated as the number of such transitions divided by the total number of transitions in these CAP sequences. To calculate PSDs of CAPs, we adopted a concept from (Gutierrez-Barragan et al., 2019) to construct the continuous-time course of each CAP in each participant. Briefly, these participant-level CAP continuous-time courses were obtained via regressing participant-level normalized (z-scores) cortical DOT data with the group-level CAP spatial maps as regressors. Then CAP PSDs were calculated using the Lomb-Scargle method (Horne and Baliunas, 1986), which could better handle DOT data of temporal discontinuities after preprocessing (Section 2.2) than the popular Welch’s method (Welch, 1967).

To examine the intra-cluster and inter-cluster spatial similarities among CAPs, the metric of spatial correlation was used as all CAP spatial maps were defined on the standard cortical surface model. To study temporal correspondences of different CAPs, two types of co-occurrence analyses were utilized. The framewise co-occurrence analysis was applied to two different sets of clustering results usually obtained from the same original data but with different processing conditions (e.g., with or without the SS regression). Specifically, we calculated, at the resolution of single timeframes, the percentage ratio of overlapped occurrences of two CAPs each from different conditions to the total occurrences of one of two CAPs (i.e., the reference CAP of investigational interests). The framewise co-occurrence analysis was not possible to investigate temporal correspondences of different CAPs from a single clustering analysis. We, therefore, performed the co-occurrence analysis on predefined time windows, termed as windowed co-occurrence analysis to be distinguished from the one above. Specifically, within a predefined window, the ratio between the numbers of occurrences of two investigated CAPs was calculated, termed as the windowed co-occurrence ratio (WCOR). These two temporal metrics were calculated in data from each participant, and their statistics were summarized at the group level.

The mixed-effect model with multiple recording sessions (i.e., 2 in the present study) per participant (Holmes, 1988, Beckmann et al., 2003) was adopted to perform statistical analyses on CAP metrics discussed above. Briefly, each metric value being tested was separately calculated in each recording session data and then averaged across two sessions for each participant. These participant-level metric values were then tested, unless otherwise mentioned, using two-tail paired t-tests with Bonferroni correction (Bonferroni, 1936).

Relationship between CAPs and GS

To probe the relationship between CAPs and GS, we tested whether the occurrences of CAPs locked toward specific phases of GS as such a phase-locking phenomenon had been reported between CAPs and GS in animal fMRI data (Gutierrez-Barragan et al., 2019). Firstly, GS was calculated as the mean time course of all cortical node time courses within the source space in each participant (Section 2.3), similar to the calculation of GS in BOLD fMRI data (Zarahn et al., 1997). Secondly, to enhance infraslow activities (Gutierrez-Barragan et al., 2019), the obtained GS time course from each participant was bandpass filtered at 0.01–0.03 Hz and converted into an analytical signal using the Hilbert Transform (Huang et al., 1998) to yield estimations of instantaneous phases and amplitudes. The instantaneous phases were tabulated for each cluster and pooled across all participants. Then their circular distributions and means were plotted using the MATLAB CircStat toolbox (Berens, 2009).

Results

Brain-wide DOT CAPs

Fig. 2A illustrates six coactivation patterns (CAP1–6) obtained with DOT RSNs. The DOT RSNs as building blocks to these CAPs are shown in Supplementary Fig. 1. The spatial distributions of these CAPs are structured, brain-wide, and of hemispheric symmetry. Three distinct spatial patterns are identified among them: the dorsal-ventral (DV) pattern (i.e., CAP1, 2), the anterior-posterior (AP) pattern (i.e., CAP4, 5), and the global pattern (i.e., CAP3, 6). The DV CAPs show the opposite HbR changes (positive vs. negative ΔHbR) along the dorsal-ventral axis. In the orthogonal direction of the dorsal-ventral axis, the separation of positive and negative HbR changes approximately runs along the inferior frontal sulcus in the frontal cortex, cuts through the sensorimotor cortex above the deep lateral fissure, and then extends toward the posterior end of the brain roughly along the lateral occipital sulcus. On the middle wall, the separation is on the upper boundary of the inferior temporal cortex. The AP CAPs show the opposite HbR changes along the anterior-posterior axis. The separation of positive and negative HbR changes appears roughly along the postcentral sulcus and cuts through the posterior end of the temporal cortex in the orthogonal direction of the anterior-posterior axis. The global CAPs exhibit either global positive or global negative HbR changes with diminished changes toward the tip of the temporal cortex.

Fig. 2.

Brain-wide DOT CAPs. (A) Group-level spatial maps of HbR CAPs when k = 6, showing three anticorrelated CAP pairs, i.e., anterior-posterior (CAP4 and CAP5), dorsal-ventral (CAP1 and CAP2), and global pairs (CAP3 and CAP6); (B) group-level cluster centroids of 6 CAPs with ICs as building blocks. Large IC weights in each row (i.e., greater than 60% of the largest absolute value in the row) are denoted by solid circles (same colors used for anticorrelated CAP pairs); (C) spatial correlation (SC) between all possible pairs for the CAPs in (A). solid red dots: SC values < −0.90, indicating three anticorrelated pairs; (D) mean SC values of anticorrelated pairs and non-anticorrelated pairs of CAPs show statistically significant differences using a t-test (*****: p < 1e-5, corrected); (E) the 3D distance map of the cluster centroids among 6 CAPs obtained with a multidimensional scaling tool from Matlab.

Visual inspections reveal anticorrelated spatial patterns of two CAPs in the DV, or AP, or global groups. This is further confirmed by the quantitative metric of pair-wise spatial correlation (SC) between all CAPs (Fig. 2C). The SC between CAP1 and CAP2, i.e., the pair with DV patterns, is −0.96, −0.91 for CAP4 and CAP5, i.e., the pair with AP patterns, and −0.97 between CAP3 and CAP6, i.e., the pair with global patterns. The mean SC of these pairs (i.e., −0.95) is statistically significantly (t-test, t(13)=6.6, p < 1e-5, corrected) lower than the mean SC (i.e., −0.01) of all other possible pairs of CAPs (Fig. 2D). The anticorrelated nature of these pairs is also evident from the weight matrix (Fig. 2B), wherein each row is the centroid vector for a cluster (i.e., a CAP) and each element in the row represents the weight coefficient of each building block (i.e., IC or DOT RSN). Specifically, the anticorrelated pairs have their significant weight coefficients (i.e., > 60% of the largest) mostly contributed by the same ICs (denoted by solid circles in Fig. 2B), but with opposite polarity in the paired CAPs.

In addition to characteristic patterns in single CAPs and paired CAPs, the 3D distance map illustrating all cluster centroids further reveals a well-structured spatial relationship among all six CAPs (Fig. 2E). Firstly, the three anticorrelated pairs have an almost orthogonal structure, with each CAP pair occupying one of three axes in the 3D space. Secondly, two CAPs in any anticorrelated pair are polarized separately toward either positive or negative directions of the dimensional axis. Lastly, the three-dimensional axes illustrated as the direct links between the three paired CAPs are then across over the same spatial point in the center. This spatial structure indicates that the dynamics of moment-to-moment brain states could be described in a 3D space by this collection of six CAPs, with each of them representing one spatially orthogonal pattern. The Euclidean distance between the projected centroids of the AP pair (2.50) was only about half of those from DV (5.07) and global (4.94) pairs.

In general, for k ranging from 6 to 20, all three spatial patterns of CAPs described above (i.e., AP, DV, and global) are retained (Supplementary Fig. 2A), albeit with a few changes. Some new emerged patterns may be due to the sub-patterns split from the original CAPs obtained when k equals 6. For example, the CAP9 from the analysis with k between 9 and 20 shows a close-to-global positive pattern with the temporal pole as moderately negative. This might be a sub-pattern in the original DV pattern with lowered positive/negative boundary. Furthermore, the CAP4 with k between 10 and 20 shows a DV pattern with elevated positive/negative boundary as compared with the original DV pattern. Some other CAPs retain such as DV or AP patterns but with reduced magnitudes (e.g., CAP18 from k between 18 and 20). The detection of sub-patterns might suggest that there are more detailed ingredients in the space represented by six CAPs, but potentially less representative at the group level.

Temporal, dynamic, and spectral measures of DOT CAPs

The occurrence rates of the AP pair (23.8% and 23.3% for CAP4 and CAP5, respectively) are significantly higher than all other CAPs (paired t-test, t(12)≥4.9, p < .001, corrected, Fig. 3A). The CAP2 from the DV pair occurs the least at 9.2%, which is significantly lower than all other CAPs (paired t-test, t(12) ≥5.5, p < .01, corrected) except the global CAPs. There is no significant difference in the dwell times among CAPs, which is about 4.02 s on average (Fig. 3B). The interval times of these CAPs are shown in Fig. 3C. The AP pair’s interval times (12.9s and 13.2s for CAP4 and CAP5, respectively) are considerably lower than all other CAPs (paired t-test, t(12) ≥ 3.1, p < .05, corrected). Except for the DV pair (paired t-test, t(12)=3.5, p < .01, corrected), there was no significant difference in interval times between two CAPs within a pair. Fig. 3D illustrates the one-step transition probability matrix among 6 CAPs. Firstly, transitions involving the AP pair of CAPs show higher transition probabilities (enclosed by the solid red lines) than other transitions. Specifically, the top 5% highest CAP transitions all involve either CAP4 or CAP5 (red dots). Furthermore, while the mean transition probability between all possible CAPs is 3.33%, the mean transition rate involving either CAP4 or CAP5 is 4.63%, significantly greater than the overall mean transition rate (t-test, t(46)=2.5, p < .05, corrected). Secondly, except for the AP pair, the mean CAP transitions (1.4%) within the DV pair (between CAP1 and CAP2) and the global pair (between CAP3 and CAP6) were significantly lower than the overall mean transition rate (t-test, t(40)=3.1, p < .01, corrected). The CAP PSDs across all participants revealed a significant amount of power in the infraslow frequency band between 0.005 and 0.03 Hz (Fig. 4) overlapping with the frequency band of GSs showing most powers, i.e., 0.008–0.06 Hz (Fig. 4). It is noted that the cutoff frequency at the low boundary (i.e., 0.005Hz in CAP PSDs) is caused by the bandpass filter used in the preprocessing. Without it, CAP PSDs should show ~1/f power spectra that are typically observed in BOLD signals (He et al., 2010).

Fig. 3.

Temporal and dynamic measures of DOT CAPs: (A) occurrence rates; (B) dwell time (dashed line: mean of all CAPs); (C) interval time; (D) one-step transition probability matrix of all CAPs from time t to t+1. Red dots: transition probabilities greater than 5%. ****: p < .0001, ***: p < .001, **: p < .01 *: p < .05, all corrected. Standard deviation bars were calculated across participants.

Fig. 4.

PSD of GS and DOT CAPs.

Temporal discontinuities of DOT data caused by bad segment rejections in preprocessing (see Section 2.2) might potentially affect the values of temporal and dynamic metrics of CAPs investigated above. An analysis was performed to investigate the potential effects by creating a dummy CAP (CAP7) in places where data segments were removed and then all temporal metrics discussed above were re-calculated (Supplementary Fig. 7). The results are similar to the ones in Fig. 3, which indicates that the effect of missing data on the temporal and dynamic metrics used in the present study does not change the distinct patterns observed in various CAPs.

Global CAPs phase-locked to GS

We hypothesized a potential linkage between GS and CAPs because 1) GS and all CAPs retain most signal power below 0.1Hz; 2) the occurrences of positive and negative global CAPs are likely leading to GS peaks and troughs; 3) no global CAPs are detected if GS is removed from DOT data in the analysis for all k values (see Supplementary Fig. 3A for k = 6). Therefore, we specifically investigated the quantitative correspondences between GS phases and CAP occurrences (Section 2.6). Fig. 5 illustrates the circular occurrence histograms of each CAP across all participants as functions of the instantaneous GS phase. Note that a phase of 0° and 180° represents the peak and trough of GS, respectively. Interestingly, the occurrences of positive and negative global CAP patterns are concentrated most pronouncedly around the GS peaks and troughs, i.e., CAP3 average phase (±std) at −6° ± 55° and CAP6 at 181° ± 55°, respectively, which are significantly different from a uniform circular distribution (Rayleigh test; CAP3: p = 0, CAP6: p <1e-69, corrected). Although the occurrences of other four CAPs also indicate non-uniform circular distributions, their standard deviations are much wider (CAP1: ± 74°, CAP2: ± 78°, CAP4: ± 79°, and CAP5: ± 77°), and their magnitudes of the resultant vectors indicating the average phases are much lower. The mean normalized magnitude of the resultant vectors at the centered average phases of global CAPs is 0.46, which is significantly greater than the mean magnitude of 0.082 from the other four CAPs (t-test, t(4)=10.7, p < .001, corrected).

Fig. 5.

GS phase histograms of DOT CAPs. Each plot is a circular histogram of occurrences of a CAP pooled from all participants as a function of the GS instantaneous phase at the resolution of 20°/bin. Solid black line: magnitude of the resultant mean vector; arc: 95% confidence intervals of estimated mean vectors; 0°: GS peaks; and 180°: GS troughs.

Plausible neuronal relevance for global CAPs

Repeating the above ICA and clustering analysis after removing GSs in participant-level data, our results indicate no detection of global CAPs (Supplementary Fig. 3A). Due to the debates on the origins of GS in the past decade (Murphy and Fox, 2017), it is important to understand whether our detected global CAPs are of neuronal relevance. Therefore, we have performed further analyses to search evidence of neuronal relevance for our detected global CAPs. While removal of GS is a well-established procedure in BOLD fMRI in dealing with non-neuronal systematic fluctuations, it is more difficult in fNIRS, particularly for patch-based systems (Wheelock et al., 2019), due to the absence of whole-brain coverage. On the other hand, non-neuronal hemodynamic systemic fluctuations in fNIRS have been often estimated using the SS channels (ICD: <15mm) (Gagnon et al., 2012, Goodwin et al., 2014) as such channels are only sensitive to superficial extracerebral tissues, e.g., scalp (Saager and Berger, 2008). The SS regression used in the present study based on such an estimate of systemic fluctuations has demonstrated to be effective in removing non-neuronal signals recorded with LS channels (Zhang et al., 2021). As a result, the detected global CAPs with the SS regression in preprocessing provide the first set of evidence toward their potential neuronal relevance.

To further support this evidence, one way is to test whether the absence of SS regression in preprocessing would result in additional global CAPs due to non-neuronal systematic fluctuations, and whether these additional global CAPs are separatable from two original global CAPs in Fig. 2. Therefore, we repeated the same analysis on the same fNIRS dataset but without the SS regression (k=6 to 20). An empirical metric, i.e., the percentage ratio of global CAPs (PRGC) out of all CAPs as a function of k (see details about PRGC in Supplementary Fig. 5), was calculated to evaluate the difference in the number of global CAPs with and without the SS regression. The PRGC metric data indicate a large difference in the number of global patterns only with the k value between 8 and 11, which supports our hypothesis that non-neuronal systematic fluctuations lead to additional global CAPs. As our previous results indicate that CAPs primarily exist in anticorrelated pairs (Fig. 2A), we performed further studies on CAPs obtained without the SS regression at k = 8 (Fig. 6A). All eight CAPs were consistently observed across participants (Supplementary Fig. 6). From the spatial patterns, CAP4 and CAP6 form the AP pair homologous to CAP4 and CAP5 in Fig. 2A, respectively, and CAP5 and CAP7 form the DV pair homologous to CAP2 and CAP1 in Fig. 2A, respectively. The results from the framewise co-occurrence analysis indicate that the occurrences of CAP4 and CAP6 in Fig. 6 are dominantly overlapped with the occurrences of CAP4 and CAP5 in Fig. 2A, respectively. Similarly, the occurrences of CAP5 and CAP7 in Fig. 6A are dominantly overlapped with the occurrences of CAP2 and CAP1 in Fig. 2A, respectively. These temporal correspondences establish the reliability of DV and AP brain states and the robustness of results from the clustering analysis in detecting them under different conditions (i.e., with or without the SS regression). Such reliability and robustness provide strong confidence about the observations below on global CAPs where changes are expected between the conditions with or without the SS regression.

Fig. 6.

Investigation of neuronal origins of global CAPs. (A) Group-level CAPs without SS regression at k=8; (B) GS phase histograms of these eight CAPs; (C) framewise co-occurrences between 8 CAPs in (A) and 6 CAPs in Fig. 2A. Percentage value on each bar: the ratio of co-occurrences between a Fig. 6A CAP (each sub-plot) and a Fig. 2A CAP (each bar in a sub-plot) to total occurrences of the Fig. 6A CAP; (D) windowed co-occurrence ratios (WCORs) between the numbers of occurrences of two paired CAPs (CAP2 vs. CAP3 for GS peaks and CAP1 vs. CAP8 for GS troughs) for variable window lengths. The WCOR metric was calculated twice for each pair, each with one of two CAPs as the reference CAP. X/Xs: mean dwell time and mean transit time of each reference CAP. The inset shows an enlarged view of these plots for the window lengths less than 200 s.

Four CAPs, i.e., CAP1, CAP2, CAP3, and CAP8, exhibit global patterns in Fig. 6A. Via visual inspections, CAP1 and CAP3 resemble the two original global CAPs in Fig. 2A (CAP6 and CAP3, respectively). Two additional global CAPs, i.e., CAP2 and CAP8, show moderately anticorrelated HbR changes at the anterior part of the temporal lobe, which are slightly different from the two global CAPs in Fig. 2A. The pair of CAP2 and CAP8 in Fig. 6A also preserved hemispherically symmetric anticorrelated patterns, like all other CAPs in both Figs. 2A and 6A. Furthermore, the four global CAPs in Fig. 6A all indicate strong phase-locking to GS, while four non-global CAPs (i.e., AP and DV pairs) are not phase-locked to GS (Fig. 6B), consistent with those observed in Fig. 2B. The co-occurrence analysis (Fig. 6C) indicates that the occurrences of CAP3 in Fig. 6A are dominantly overlapped with the occurrences of CAP3 in Fig. 2A. The occurrences of CAP1 in Fig. 6A are primarily overlapped with the global state CAP6 and a non-global state CAP5 in Fig. 2A. On the contrary, both CAP2 and CAP8 in Fig. 6A are significantly overlapped with four and three CAPs (>10% for each) in Fig. 2A, respectively, including both global and non-global CAPs. These spatial and temporal correspondences indicate two types of global CAPs in Fig. 6A are spatiotemporally distinguishable. Specifically, the pair of CAP1 and CAP3 are potentially of neuronal relevance (like the original global CAPs in Fig. 2A). In contrast, the two additional global CAPs, i.e., CAP2 and CAP8, are introduced possibly due to non-neuronal systematic fluctuations.

As the two types of global CAP pairs (i.e., the original global CAP pair and the additional global CAP pair) are both phase-locked to GS peaks and troughs, we next investigated if the two types of global CAP pairs occupied different subsets of GS peaks (for the positive global ones: CAP2 vs. CAP3) and troughs (for the negative global ones: CAP1 vs. CAP8) using a windowed co-occurrence analysis (see Section 2.5). Specifically, we centered a window, which varied from 0.5 to 100 s in the step of a single timeframe (0.16 s), to a GS peak (or trough) with the original global positive CAP (or the original global negative CAP) as the reference to calculate values of the metric WCOR and then repeated the calculation with the additional global positive CAP (or the additional global negative CAP) as the reference. We focused the resulted group-level WCOR values (Fig. 6D) on two window sizes: <=6 s that is halfway between the mean dwell time (~4 s, Fig. 3B) of four global CAPs and the mean time between a pair of neighboring peak and trough of GS (8.7 s); and >40 s as the mean interval time of the four global CAPs is about 40 s (Fig. 3C). The WCOR values corresponding to the window size <=6 s are below 12% regardless of the type of the reference CAP, which indicates that each GS peak or trough and its surrounding timeframes are mainly occupied by one type of global CAPs (see examples of overlay of GS peaks, troughs, and CAPs in Supplementary Fig. 4). In other words, two types of global CAPs are phase-locked to different subsets of GS peaks and troughs. For window size >40 s, the WCOR values are roughly in the range between 30% and 40% regardless of the type of the reference CAP. Moreover, these values plateau between 40 and 100 s window sizes, which indicate that each type of CAPs largely occur in clusters. In other words, one type of phase-locked global positive or negative CAP largely occurs consecutively. These observations indicate that the phase-locked occurrences of two types of global CAPs are temporally separable in relation to GS peaks/troughs, which further support that the two types of global CAPs have different origins.

Participant-level brain-wide CAPs

We investigated if CAPs could be detected at the participant level. Fig. 7A illustrates the spatial maps of six CAPs obtained from data in one participant. Although spatially noisier, they show a strong resemblance to the group-level CAPs (Fig. 2A). Three pairs of distinguished patterns are identifiable, including AP (CAP4 and CAP5), DV (CAP1 and CAP2), and global patterns (CAP3 and CAP6), which further indicate anticorrelated features in each pair (Fig. 7B and C). Moreover, the participant-level spatial maps of all six CAPs are matched with their corresponding group-level maps with statistically significant spatial correlations greater than zero (paired t-test, p < .00001, corrected, Fig. 7D and Supplementary Fig. 3C). Furthermore, all six CAPs have been identified in all participants with only one exception of one session on one participant (see Supplementary Fig. 4). These findings suggest that FC variations in the brain can be consistently observed at the participant level.

Fig. 7.

Participant-level brain-wide DOT CAPs. (A) Spatial maps of CAPs from a representative participant corresponding to those in Fig. 2A; (B) spatial correlation (SC) between all possible pairs for the CAPs in (A); solid red dots: SC values < −0.79, indicating three anticorrelated pairs; (C) mean SC values of anticorrelated pairs and non-anticorrelated pairs of CAPs show statistically significant differences using a t-test (***: p < .001, corrected); (D) SCs between the group-level and participant-level spatial maps of CAPs show significantly higher values than zero (paired t-test, *****: p < .00001, corrected).

Discussion

Brain-wide coactivations indicating recurring brain states

The present study for the first time reports the detection of brain-wide CAPs in hemodynamic signals recorded with a wearable whole-brain high-density optical imaging system operating in the near-infrared range, known as DOT, on healthy human participants. Our results of DOT CAPs show distinctive spatiotemporal features of CAPs that resemble findings in BOLD fMRI data. Firstly, our data indicate that resting-state fluctuations of whole-brain hemodynamic signals can be accounted with a relatively small number of brain-wide functional states, i.e., 6 CAPs, which is consistent with the number of CAPs reported in recent fMRI literature, e.g., 6 (Yang et al., 2021) and 8 CAPs (Janes et al., 2020). Small numbers of brain-wide functional states have been similarly reported in animal BOLD data, e.g., 3 (Liang et al., 2015) and 6 CAPs (Gutierrez-Barragan et al., 2019). These CAP numbers are at the similar level as the number of dynamic brain states estimated from cross-RSN functional connectivity (FC) in resting-state fMRI, e.g., 7 (Allen et al., 2014). Since the CAPs in the present study are defined as a weighted combination from a set of DOT RSNs (Supplementary Fig. 1), each CAP represents interactions among multiple RSNs (and transient in nature due to their short dwell times), which is, in concept, similar to the cross-RSN FC (Allen et al., 2014). However, it is important to note that brain states defined in cross-RSN FCs are estimated based on the metric of correlation over preselected time windows, which are different from the CAP analysis based on individual timeframes. Interestingly, multiple RSNs have been reported to transiently change together in the concept of CAPs in both human (Yang et al., 2021) and animal fMRI data (Gutierrez-Barragan et al., 2019). On the other hand, large numbers of CAPs have also been reported, e.g., up to 16 (Chen et al., 2015), 20 (Karahanoğlu and Van De Ville, 2015), and 30 CAPs (Liu et al., 2013), from selected fMRI timeframe data showing local intensity peaks in a priori anatomically-seeded regions. Moreover, the CAPs based on pre-determined regions of interest have established close relationships to RSNs (Liu et al., 2013). The CAPs reported in the present study are identified using a differernt approach via clustering all timeframe data without being restricted by polarity and intensity of local hemodynamic signals, which allow regionally unbiased detections of tranisent recurring patterns. Therefore, the CAPs in the present study are believed to be more representive of brain-wide functional states, which has been termed as manifold CAPs (Gutierrez-Barragan et al., 2019) as each CAP is defined by a subspace (the mathemical definition of a brain state) spanned by the collection of timeframes assigned from clustering in the total 20-dimensional space (i.e., 20 ICs as building blocks in the present study).

Secondly, the spatial maps of brain states in CAPs show bilaterally symmetrical patterns. In humans (Yeo et al., 2011, Smith et al., 2009) and anesthetic monkeys (Vincent et al., 2007), hemodynamic resting-state FC often exhibits bilateral symmetry across two cerebral hemispheres, which is even preserved in humans with complete agenesis of the corpus callosum (Tyszka et al., 2011). The majority of coactivation studies with BOLD fMRI have also reported CAPs with symmetrical spatial patterns in humans (Liu et al., 2013, Liu et al., 2018a, Karahanoğlu and Van De Ville, 2015, Chen et al., 2015) and animals (Gutierrez-Barragan et al., 2019, Liang et al., 2015). It is important to note that spatial patterns of ICs, as the building blocks of CAPs, mostly exhibit asymmetrical patterns (Supplementary Fig. 1), albeit in complementary pairs (e.g., the IC2-IC20 pair for the dorsal attention network and the IC13-IC15 for the somatomotor network), which is consistent with previously reported DOT RSNs (Khan et al., 2021b). This observation indicates that the bilateral symmetry of CAPs could result from temporally coordinated activations of multiple RSNs. It is plausible that the bilateral coordinations across hemispheres are transient, which are better detected in the framewise clustering analysis but not in the static ICA analysis. This hemispheric symmetry of CAPs in the present study is preserved despite increasing the k value (Supplementary Fig. 2A), ruling out the possibility of it being forced due to a small number of CAPs, e.g., 6.

Thirdly, DOT CAPs show strongly anticorrelated spatial patterns in pairs, i.e., AP, DV, and global pairs (Fig. 2). In fMRI literature, strong brain-wide anticorrelation patterns are widely existing in resting-state FCs and CAPs. In humans, the phenomenon of anticorrelation in RSNs has been observed using the conventional windowed correlational method, e.g., between sensory regions and default mode regions (Allen et al., 2014) and between distributed task-positive and task-negative brain networks (Fox et al., 2005) (Liu et al., 2013). In animals, the phenomenon of anticorrelation has been reported with CAPs in resting mice (Gutierrez-Barragan et al., 2019). The anticorrelations might be expressed because of separating neuronal processes responsible for competing interests (Fox et al., 2005), such as task focus vs. stimulus-independent thought (Giambra, 1995). Although there have been doubts that the phenomenon of anticorrelation might be artificially introduced due to the GS regression (Murphy et al., 2009), the anticorrelation patterns have been reported in the absence of GS regression (Chai et al., 2012). The results in the present study also show strong anticorrelations in the CAP spatial patterns obtained in the absence of GS regression (Supplementary Fig. 3A and B, see more discussions in Section 4.3).

Lastly, our results further indicate one pair of DOT RSNs showing global positive and negative patterns, which has been reported when hemodynamic data are analyzed using techniques other than correlation-based methods (Yousefi et al., 2018, Matsui et al., 2019, Ma et al., 2020). For example, the global coactivation patterns have been identified in human hemodynamic data without removing GS using a quasiperiodic pattern searching algorithm (Yousefi et al., 2018). Transient global coactivations have been observed in searching globally propagating waves in both hemodynamic and neural signals from rodents (Matsui et al., 2019). A recent human EEG study (Ding et al., 2021) further identifies global coactivations in neuronal electrical signals.

Spatiotemporal metrics of CAPs indicate structured transitions among brain states

The philosophy of coactivation analysis in studying dynamics of spontaneous brain fluctuations is principally different from conventional correlational-based static/dynamic FC methods. The coactivation analysis utilizes single timeframe data as the basic unit of analysis rather than entire or windowed recordings (Liu et al., 2018b). Even in dFCs, the selection of windows is limited to the minimal data length mathematically required to warrant the estimation accuracy of correlations and its non-optimal selection may bias reconstructed brain networks (Liu et al., 2013). Coactivation analysis is, however, independent of the choice of window size. All of our detected CAPs show consistent mean dwell times at 4.02s (Fig. 3B), which is significantly lower than window sizes that are usually used in dFC fMRI studies, e.g., 120s (Savva et al., 2019). The dwell time in the present study is in close agreement with the previous coactivation studies using fMRI, e.g., 4.0–7.6s in humans (Karahanoğlu and Van De Ville, 2015), 4s in mice (Gutierrez-Barragan et al., 2019), and ~10s in rats (Liang et al., 2015). These dwell times are on the order of hemodynamic responses (Friston et al., 2000) and several orders greater than the dwell times of CAPs based on electrical signals, e.g., EEG (Ding et al., 2021, Li et al., 2022), and brain states studied by other methods, e.g., microstates (Lehmann et al., 1987, Michel and Koenig, 2018), which range from tens to hundreds of milliseconds. We speculate that the large difference between the dwell times of transient brain state events between hemodynamic and electrical signals could be mainly determined by the “sluggish” nature of the cerebral hemodynamic response following neuronal activity, which also supports the fact that the least variations are observed in the dwell time than all other temporal and dynamic metrics among different CAPs (Fig. 3 and see more discussions below) in the present study. However, more work is needed to verify this speculation in future studies especially if a quantitative relationship between hemodynamic and electrical-based CAPs could be established with concurrently recorded EEG and fMRI signals (Yuan et al., 2016).

Beyond the metric of dwell time, all other temporal and dynamic metrics indicate significant differences among the six CAPs. These differences are consistent across different metrics, as well as with the spatial metric, and collectively suggest the existence of a structured transition pattern among the six brain states represented by CAPs. Interestingly, the anterior-posterior CAP pair is heavily involved with transitions among all brain states. This important role of the anterior-posterior CAPs in dynamic brain state transitions is supported by their significantly higher occurrence rates and significantly lower interval times as compared with other CAPs. This role is further supported by the 3D distance map (i.e., the spatial metric) between all CAPs (Fig. 2E), where the distance between two opposite anterior-posterior CAPs is significantly shorter than other anticorrelated CAP pairs. Furthermore, the direct transitions between the two anticorrelated CAPs within the dorsal-ventral pair, as well as within the global pair, are among the lowest, which further suggest the anterior-posterior CAP pair might serve as the intermediate brain states during the transition of brain states in opposite directions.

We speculate that these observed structured transitions among CAPs could indicate the presence of globally propagating waves that have been well reported in fMRI (Raut et al., 2020, Raut et al., 2021, Majeed et al., 2011), optical data (Matsui et al., 2016, Matsui et al., 2019), MEG data (Takeda et al., 2021), and electrical data (Massimini et al., 2004, Huang et al., 2010, Takeda et al., 2021). Globally propagating waves made up of a sequence of brain-wide CAPs in both hemodynamic and neuronal spiking data predominantly originating from the anterior brain and propagating toward the posterior end of the brain have been reported in mice (Matsui et al., 2016). A cortical traveling wave originating from the prefrontal-orbitofrontal cortical regions and traveling to the posterior cortical regions is observed in sleeping human participants in a high-density (178–186 channels) EEG study (Massimini et al., 2004). Also, a global anterior-to-posterior phase shift in gamma oscillations (40 Hz) is reported with MEG data (Ribary et al., 1991). Our detected anterior-posterior anticorrelated CAPs might be two moment-to-moment snapshots of brain states with the most representative spatial patterns along these anterior-posterior propagating waves. Similarly, we speculate that the detected dorsal-ventral anticorrealated CAPs might reflect representative brain states of cortical dorsoventral propagating waves, e.g., in fMRI of rats and humans (Majeed et al., 2011) and as alpha waves in electrocorticography (Bahramisharif et al., 2013). However, it is cautioned that the present study does not provide direct evidences of propagations and future studies are needed to integrate the coactivation analysis together with other sequence analysis methods (Majeed et al., 2011) to further clarify existing propagation pathways among brain states indicated by CAPs.

Plausible origins of CAPs and GS

We believe that our detected DOT CAPs reflect nonstationary fluctuations in spontaneous neuronal activity in human brains (Foster and Wilson, 2006, Logothetis et al., 2012, Matsui et al., 2019) and, therefore, are of neuronal origins. Firstly, as both CAPs and dFCs study temporal fluctuations in hemodynamic data in resting brains, the neurophysiological evidence of dFC has been extensively reported in both humans (Chang et al., 2013, Tagliazucchi et al., 2012) and animals (Lu et al., 2007, Thompson et al., 2013). More recently, a study on lightly anesthetized mice has reported CAPs in simultaneously recorded calcium (i.e., neuronal) and hemodynamic signals, as well as their temporal correspondences, providing direct neuronal evidence of CAPs (Matsui et al., 2019). Studies using other brain signals (i.e., electrical) have directly reported neuronal CAPs in animals (Matsui et al., 2016, Mohajerani et al., 2013, Vanni and Murphy, 2014) and humans (Ding et al., 2021). As discussed above (see Sections 4.1 and 4.2), the spatial, temporal, and transitional data of CAPs in our present study corresponds well with hemodynamic CAPs with known neuronal origins.

In addition, our present study has provided indirect evidence about the plausible neuronal relevance of two global CAPs. These evidences include: (1) a pair of original global CAPs (Fig. 2A) detected with an efficient preprocessing technique in place, i.e., the SS regression, to remove non-neuronal systematic hemodynamic fluctuations (Cheong et al., 2020, Zhang et al., 2019, Zhang et al., 2021); (2) a pair of additional global CAPs observed in the absence of SS regression (corresponding to non-neuronal systematic hemodynamic fluctuations); (3) distinguishable patterns in spatial maps, temporal occurrences, and relations to GS between the pair of original global CAPs and the pair of additional non-neuronal global CAPs; and 4) consistent detections of global CAPs in conditions with a wide range of selection of parameters (Figs. 2A, 6A, 7A, and Supplementary Fig. 2A). Further more, the neuronal relevance of this pair of original CAPs is supported by recent literature on transient electrical global patterns in rats (Schwalm et al., 2017) and electrical global CAPs in humans (Ding et al., 2021).

The present study further indicates that GS contains neuronally relevant information due to the close relationship between GS and global CAPs, i.e., global CAPs phase-locked to GS peaks and troughs (Figs. 5 and 6B) and no detection of global CAPs after regressing GS out (Supplementary Fig. 3A). Given that GS is calculated as the spatial mean of DOT data in each time frame, it is not a surprise that the positive and negative global CAPs are corresponding to GS peaks and troughs, respectively, as no cancellations between regional activations occur due to the globally same magnitude polarity. The neuronal relevance of GS has been suggested in recent literature, in which GS of fMRI has been linked to vigilance (Wong et al., 2013, Chen et al., 2020), glucose metabolism (Thompson et al., 2016), and arousal (Liu et al., 2018a). In addition, our recent whole-head EEG-fNIRS study in humans has shown that the fNIRS GS is negatively associated with EEG vigilance in resting states (Chen et al., 2020). Thus, this finding in the present study suggests that care need to be exercised when removing GS of hemodynamic data (Schölvinck et al., 2010) in either fMRI or fNIRS.

It is also important to note that, after removing GS, majority of observations regarding CAPs discussed have been reserved except the loss of two original global CAPs. Specifically, in results at k=6, the CAP pairs of AP patterns (CAP1, 6) and DV patterns (CAP2, 4) (Supplementary Fig. 3A) are still identified. Hemispherical symmetric patterns in all CAPs (Supplementary Fig. 3A) and strong anticorrelations within the AP and DV pairs (Supplementary Fig. 3B) are also observed. Similarly, the CAP pair of AP patterns has more important roles in the transitions of brain states defined by all CAPs (Supplementary Fig. 3D) and the CAP pair of DV patterns has the lowest occurrence rates and the highest interval times (Supplementary Fig. 3F). All CAPs show similar PSDs (Supplementary Fig. 3E) as observed in Fig. 4 and similar participant-level detections (Supplementary Fig. 3C).

Relationship among coactivations, structured transitions, and global activity patterns

Coactivations, brain-wide coordinated transitions, and global activity patterns (or global signals), separately, have all been well documented previously in hemodynamic literature, mainly with fMRI as described above. The co-existences of any two of these three phenomena have been reported as well. (Matsui et al., 2016) have further demonstrated all three phenomena in both neuronal spiking and hemodynamic data in mice. The present study extends the previous research work in animals and reports the co-existence of these phenomena in human brains and their relationship. Specifically, our study focuses initially on the investigation of coactivation patterns via a framewise clustering technique on independent component sources identified in DOT data, which leads to the detections of brain-wide patterns indicating recurring functional brain states. A pair of identified brain states actually exhibit global positive (coactivation) and negative (co-deactivation) patterns, leading to the definition of global activity patterns and their existences as one kind of recurring functional brain state. We further demonstrate that these two global activity patterns are phase-locked to global signals and are responsible for the generation of a significant subset of global signal peaks and troughs, while non-neuronal systematic fluctuations in breathing and circulation are responsible for the generation of other global signal peaks and troughs. Finally, we present evidence about the existence of a structured transitions among identified brain states, in which most transitions among brain states involve the pair of anterior-posterior CAPs (Fig. 3D). Furthermore, brain states with both anterior-posterior and dorsal-ventral patterns might reflect representative spatial patterns during the well-reported globally propagating waves in the anterior-posterior and dorsal-ventral directions, respectively. Global activity pattern has been observed as snapshots of brain states in globally propagating waves (Matsui et al., 2016). It is therefore plausible that all identified CAPs represent recurring brain states resulted from structured propagations. By showing the relationship among coactivations, brain state transitions, and global activity patterns, the present study brings together a broad range of studies in animals and humans to form the basis for understanding these phenomena on a unified foundation.

Limitations

The present study was limited in the following aspects. Firstly, the mapping of coactivation with our BW-DOT technique only covers the neocortex, but not deep brain structures (e.g., hippocampus and amygdala) due to exponential light attenuations as a function of penetration depth and therefore significantly reduced sensitivity in deep parts of the brain (Dehghani et al., 2009b), which is the most significant drawback as compared with fMRI. In terms of density of optode distributions, our BW-DOT system is lower than the reported ultrahigh-density patch-based DOT systems (Eggebrecht et al., 2012), which is a trade-off between the whole-head coverage and ultrahigh spatial resolution (Khan, 2021). Relatively low spatial density of optodes might be a contributing factor to the small number of brain states (i.e., 6) identified in the present study, which should be investigated in future studies after addressing the technical challenge of having an ultrahigh-density whole-head-coverage DOT systems. It is further noted that our whole-head coverage also indicate lower sensitivity in the ventral parts of the temporal cortex (Khan, 2021, Zhang et al., 2021). Secondly, the sampling frequency of BW-DOT at 6.25 Hz was several times higher than typical fMRI scans (~0.5 Hz), which was useful for capturing faster hemodynamics (Polimeni and Lewis, 2021) and reducing aliasing noise. However, it was still fundamentally limited in revealing fast brain dynamics, e.g., in alpha (Bahramisharif et al., 2013, Shou et al., 2020) and gamma bands (Ribary et al., 1991). Thirdly, like fMRI, DOT interrogates human brain activity indirectly on secondary hemodynamic responses of neuronal events and, therefore, neuronal origins of brain-wide CAPs and global brain activity, as well as their propagation characteristics, need to be further investigated using primary response data of neuronal events, i.e., electrical activity (Ding et al., 2021), with technologies like EEG, MEG, and ECoG. In particular, to directly investigate the relationship between brain-wide neuronal events and their hemodynamic correspondences would need simultaneous recording techniques, including concurrent EEG and fMRI (Yuan et al., 2016, Yuan et al., 2012) or concurrent EEG and fNIRS (Khan et al., 2021b, Ahn and Jun, 2017, Chen et al., 2020). Fourthly, it should be noted that our findings only provide indirect evidence to the neuronal relevance of original global CAPs as they are different from brain-wide patterns caused by non-neuronal systematic variations in hemodynamic signals. Their neuronal origins still need to be directly interrogated with simultaneously measured neuronal activities, such as electrophysiology, in human, using concurrent recording techniques discussed above. Moreover, the physiological relevance of these global CAPs needs to be established via identifying their behavioral or cognitive correlates (Preti et al., 2017). Finally, these brain-wide hemodynamic phenomena observed in this pilot study need to be replicated using a larger number of participants.

Conclusion

The present pilot study utilized a wearable brain-wide DOT technology and data-driven methods to investigate spatio-temporal dynamic patterns in hemodynamic signals from spontaneous human brain activity. Our results demonstrate that brain-wide spontaneous hemodynamic signals from DOT can be described by a few co-(de)activation patterns representing recurring and transitioning brain states. Global positive and negative brain states, different from global non-neuronal patterns, are identified as two co-(de)activation patterns indicating transient activations over the entire cortex, which are phased-locked with a subset of global signal peaks and troughs. Collectively, these results indicate the existence and interaction of three phenomena, i.e., CAPs, global brain activity, and brain-wide structured transitions, in human hemodynamic data, which extend existing findings on similar transient neuronal events recently reported in animal fMRI data. Our study has important implications in understanding the spontaneous activities and the fundamental organization of resting state networks in the human brain.