Abstract
Spontaneous neural activities are endowed with specific patterning characterized by synchronizations within functionally relevant distant regions that are termed as resting-state networks (RSNs). Although the mechanisms that organize the large-scale neural systems are still largely unknown, recent studies have proposed a hypothesis that network-specific coactivations indeed emerge as the result of globally propagating neural activities with specific paths of transmission. However, the extent to which such a centralized global regulation, rather than network-specific control, contributes to the RSN synchronization remains unknown. In the present study, we investigated the contribution from each mechanism by directly identifying the global as well as local component of resting-state functional MRI (fMRI) data provided by human connectome project, using temporal independent component analysis (ICA). Based on the spatial distribution pattern, each ICA component was classified as global or local. Time lag mapping of each IC revealed several paths of global or semi-global propagations that are partially overlapping yet spatially distinct to each other. Consistent with previous studies, the time lag of global oscillation, although being less spatially homogenous than what was assumed to be, contributed to the RSN synchronization. However, an equivalent contribution was also shown on the part of the more locally confined activities that are independent to each other. While allowing the view that network-specific coactivation occurs as part of the sequences of global neural activities, these results further confirm an equally important role of the network-specific regulation for its coactivation, regardless of whether vascular artifacts contaminate the global component in fMRI measures.
Introduction
Once considered to be a noisy, stochastic process, spontaneous activity of the cortical neuron is now understood to be by no means random but is endowed with specific patterning that reflects the functional architecture of the underlying network at the level of micro- or meso-circuits (Tsodyks et al., 1999; Kenet et al., 2003). Over the last two decades, it has become apparent that this is analogously true at the level of large-scale networks that are defined using resting-state functional magnetic resonance imaging (rs-fMRI) (Biswal et al., 1995; Fox and Raichle, 2007). The spatial patterns identified as areas with synchronous oscillation of the blood oxygenation-level dependent (BOLD) signal are termed as resting-state networks (RSNs) (Fox et al., 2005). These networks are closely related to anatomical connectivity among the neural subsystems that have been revealed by a wide variety of visual, sensorimotor, and cognitive task paradigms (Vincent et al., 2007; Zhang et al., 2010). However, the neurophysiology of the phenomenon, or the mechanism that controls and coordinates the intrinsic synchronization across distributed neural systems, largely remains to be established.
While conventional rs-fMRI analyses based on seed-based correlation or independent component analysis (ICA) implicitly assume that the spatial distribution of the synchronous neural activity is temporally constant, animal studies have revealed that spontaneous neural activity is spatiotemporally structured, and propagating waves of activity have been recorded in a variety of species (for review, see Muller and Destexhe, 2012). Neuronal membrane potential in the cortex undergoes a spontaneous transition between up and down states in the absence of sensory inputs (Steriade et al., 1993; Lampl et al., 1999; Petersen et al., 2003; Shu et al., 2003). Population activity of the neurons during the up state manifests as propagating waves not only within a part of the cortex (Petersen et al., 2003; Ferezou et al., 2007; Xu et al., 2007; Civillico and Contreras, 2012), but also throughout the entire brain (Stroh et al., 2013). The spatiotemporal dynamics of the low-frequency oscillation have also been identified by examining the repetitive spatiotemporal patterns (Majeed et al., 2009, 2011; Takeda et al., 2016; Belloy et al., 2018; Abbas et al., 2019) or by analyzing the time lag structures of the rs-fMRI data (Mitra et al., 2014, 2015a; Amemiya et al., 2016; Matsui et al., 2016). An intriguing hypothesis proposed by one of those studies is that the RSN synchronization indeed emerges as the result of several independent global propagations of spontaneous neural activity (Mitra et al., 2015a). Using synthetic time series embedded with the measured time lag structures of the rs-fMRI data, Mitra et al. (2015a) showed that the functional connectivity (FC) matrix representing the RSN synchronization could be reconstructed to a fair approximation. In support of this idea, a recent animal study also showed that a global wave of spontaneous neuronal activity propagating across the networks contributes to within-network coactivations of the neurons that correspond to RSN synchronization (Matsui et al., 2016). Based on these findings it follows that seemingly independent RSN activity can be viewed as being controlled by a single centralized mechanism, through global wave(s) of activity that regulate and constrain the relative relationships of the network activity by determining the order and timing of the activation of each network. However, it remains unclear if this is the sole mechanism that gives rise to the RSN synchronization, as suggested by those studies (Mitra et al., 2015b; Matsui et al., 2016). An alternative, thought non-exclusive, origin of the synchronization would be network-specific coactivations among the neural populations confined within each network (Mohajerani et al., 2013; Ponce-Alvarez et al., 2015). In the neural system, it is generally supposed that diverse physiological mechanisms coexist for rhythm generation and population synchronization for which different levels of integration interact closely with each other (Ivanchenko et al., 2008; Harris-Warrick, 2010; Wang, 2010). For example, in the respiratory central pattern generator of the mammals, rhythm generation is dependent on the endogenously oscillatory neurons that serve as pacemaker, as well as the pattern of synaptic connections within the network that forms a network pacemaker (hybrid pacemaker-network mechanism) (Calabrese, 1998; Rybak et al., 2004, 2007; Sohal et al., 2006; Johnson et al., 2007).
It seems possible, therefore, that multiple mechanisms – namely, global propagation and local synchronization – contribute to the emergence of the coherent RSN activity that characterizes the functional architecture of the brain. In order to address this question, it is imperative to evaluate not only the paths but also the whole picture of the traveling waves. We thus started by identifying the signal time course of the global waves by applying the temporal ICA to rs-fMRI data. In fMRI, virtually all applications of ICA use spatial rather than temporal ICA. Although spatial ICA is suitable for the separation of the spatially distinct activations from each other, temporal ICA would be more appropriate if the aim is to find functionally independent and spatially overlapping activities (Smith et al., 2012), such as what we assume to be multiple global waves.
In contrast to previous studies focusing on estimating the paths of traveling waves by analyzing the signal time lag (Mitra et al., 2014, 2015a; Amemiya et al., 2016; Matsui et al., 2016), direct detection of the traveling waves enables us to infer the likelihood that each mechanism contributes to the emergence of network synchronization, as well as to map the magnitude of each type of activity in each region. Moreover, identification of individual traveling waves allows accurate estimation of the time lag structures, in contrast to previous studies employing a decomposition approach, regardless of the validity of the assumptions that the time lags of multiple waves can be linearly superposed, or that the paths of traveling waves are spatially independent. By comparing the correlation matrix of both the global and the local component to that of the FC matrix, contribution of each type of activity to the RSN synchronization that characterizes the resting state FC was evaluated.
Materials and Methods
Overview
A summary of the analysis is presented as a schematic in Figure 1 to provide an overview of the study. We used data from the WU-Minn Human Connectome Project (HCP) young healthy adults (ages 22–35) S1200 release that provides paired dataset of the same group of subjects (day 1 and day 2). All preprocessing and data analyses were performed for each dataset, respectively, in the same way.
FIGURE 1
HCP Data
Data of 50 subjects who underwent 3 T resting-state fMRI sessions without quality control issues, and whose mean framewise displacement was less than 0.2 mm were included for the analysis. The number of subject was determined by the maximum size of the data that could be processed by a core program for ICA, Multivariate Exploratory Linear Decomposition into Independent Components (MELODIC) (Beckmann and Smith, 2004). HCP imaging and pre-processing protocol have been previously described in detail (Glasser et al., 2013; Smith et al., 2013; Van Essen et al., 2013; Griffanti et al., 2014; Salimi-Khorshidi et al., 2014). In brief, resting-state fMRI data were acquired using a single customized Siemens 3 T scanner housed at Washington University in St. Louis, using a standard 32-channel receive head coil, with 2.0 mm isotropic spatial resolution, 0.72 s repetition time (TR), and 1200 frames, i.e., 14.4 min per run. For each subject, and for each session, two runs with reversed phase encoding directions, RL or LR, with the order counterbalanced across each of two sessions, were acquired (WU-Minn HCP 1200 Subjects Data Release Reference Manual), and the geometric distortions were corrected using spin echo field map EPI scans (Glasser et al., 2013). The data were then subjected to spatial ICA using MELODIC with automatic dimensionality estimation. Using FMRIB’s ICA-based X-noisifier (FIX) (Griffanti et al., 2014; Salimi-Khorshidi et al., 2014) that is a machine learning classifier trained on HCP data, spatially specific noise components were identified and removed for each run. Then 24 movement regressors were further regressed out of the data (Smith et al., 2013; Griffanti et al., 2014; Salimi-Khorshidi et al., 2014).
Data Analysis
Further Pre-processing
Further pre-processing and analysis of the data were performed using tools from SPM12 software1, AFNI libraries2 and in-house scripts written and implemented in Matlab 9.3 (MathWorks, Natick, MA, United States). Linear trends were removed from the HCP data that had been processed with subject-level ICA noise reduction (sICA + FIX), and the data were band-pass filtered at 0.01–0.1 Hz. The pre-processed data were temporally concatenated across runs to create a single 4D dataset of 120,000 timepoints for test and re-test dataset, respectively.
Temporal-ICA
For temporal ICA decomposition of the data, we employed a strategy adapted from Smith et al. (2012) to perform group-wise spatial ICA in advance of the final temporal ICA. The spatial ICA allows further identification of artifact components at the group level, as well as to achieve a high-dimensional functional parcelation of the group data, which reduces the dimensionality of the data to feed into temporal ICA. The overall ICA analysis is described as follows:
where X is the data matrix of size Voxels × Time points, Ss are the spatial maps estimated by spatial ICA, K is the number of spatial ICA components that were subsequently fed into temporal ICA, after removing noise components. L is the number of temporal ICA components. St is the decomposed time series (ICA sources), and At is the central mixing matrix of temporal ICA. E combines noise and artifact aspects of the data (Smith et al., 2012).
Group-wise spatial ICA was performed using MELODIC with automatic dimensionality estimation. Of the 61 and 62 ICs generated from dataset 1 and 2, three ICs were classified as artifacts on the basis of their spatial features, respectively. Specifically, activation patterns clearly outlining the intensity edges of the gray matter were classified as noise (Smith et al., 2012; Salimi-Khorshidi et al., 2014; Pruim et al., 2015). For the remaining components, functional nodes’ time series were computed using dual regression technique (Filippini et al., 2009), and fed into temporal ICA.
For temporal ICA, we used Icasso algorithm (Himberg et al., 2004) to estimate the most appropriate decomposition yielding a set of reproducible IC clusters. For all possible dimensions or number of components, Icasso was ran with both resampling mode that uses random initial condition as well as bootstrapping of the data for 100 times, which pooled all temporal ICA estimates using FastICA (Hyvärinen, 1999) with tanh non-linearity and a symmetric decorrelation approach. We chose the decomposition yielding the maximum number of clusters of reproducible components that gives a stability index Iq larger than 0.5. Iq is computed as the difference between the average intracluster similarities and average intercluster similarities, which reflects the compactness and isolation of a cluster (Himberg et al., 2004). For each dataset, 28 and 30 reproducible clusters were found, respectively.
Identification of the Global Waves
For all temporal ICs, time series of each run, once concatenated to be subjected to a temporal ICA was deconcatenated so that following analyses can be performed for each run separately unless otherwise noted. Firstly, spatiotemporal patterns or paths of traveling waves were estimated for each IC by computing the relative time lag t that gives the best positive fit between each voxel’s time series and the time-shifted (±6.3 s or ±9 TR) IC (time series) using cross-correlation analysis. As in our previous study, 12 s limit of propagation delay was set to include whole-brain vascular time lag that can range up to nine seconds (Amemiya et al., 2016). Parabolic interpolation (Meijering, 2002) was further applied to locate the peak time lag t’ using the extremum t, as well as the two nearest points given by the cross-correlation analysis.
The magnitude map of each IC was then computed as the Pearson’s correlation coefficients between each voxel’s time series and the IC time series that was shifted as much as t’ using sinc interpolation to give the maximum correlation. Classification of ICs was based on the spatial distribution pattern of each component. We performed template matching using 21 RSN templates (Smith et al., 2012), as well as a whole-brain signal template that is the correlation map of the whole-brain mean signal, averaged over 100 runs (Supplementary Figure S1). The whole-brain signal was computed as the average signal within a gray matter mask was created by thresholding MNI template at 10% or larger probability of being gray matter. Pearson’s correlation coefficients were computed between each IC and each of the 22 templates within the mask. Any component that is more similar to the whole-brain signal template than any of the RSN templates (i.e., giving a greater correlation coefficient with a whole-brain signal template) was classified as global. Note that our study focuses particularly on examining the existence of coactivations restricted within each RSN, in addition to the globally propagating activities that were assumed to exist throughout the brain and treated as such in the analysis of previous studies. In this context, it would certainly make sense to identify any component whose distribution is restricted within any functionally distinct areas as being local as opposed to a functionally and spatially less specific global (or more precisely semi-global) pattern.
Using the time-shifted global ICs as regressors, linear regression analysis was performed for each voxel’s time series. The global component was then computed as the integral of all global ICs by summing up the shifted time series multiplied by the corresponding regression coefficients. The rest of the signal was classified as local contribution (Figure 1).
Estimation of the Likelihood That Each Mechanism Gives Rise to the RSN Synchronization
In order to evaluate the contribution of each mechanism to the emergence of RSN synchronization, we compared regions of interest (ROI)-wise correlation matrices given by each component using a set of 132 ROIs provided as part of the CONN toolbox (Whitfield-Gabrieli and Nieto-Castanon, 2012) that were originally defined from FSL Harvard-Oxford Atlas maximum likelihood cortical or subcortical atlas and cerebellar parcelation from AAL atlas.
For each subject’s each run, FC matrix was obtained by computing the Pearson’s correlation coefficient between each possible couple of ROI’s mean time series (i.e., global + local component). Similarly, three types of correlation matrices were computed by correlating the time series of (1) global component, (2) local component, (3) global component reconstructed to reflect only time lag of each global IC, respectively. We correlated each matrix against the FC matrix run-wisely, using the correlations above the diagonal of each matrix, transformed to Fisher’s Z and tested by using a two-tailed t-test over runs against the null hypothesis of no correlation. Next, we examined whether the time lag itself contributes to the RSN synchronization by computing the correlation between the FC matrix and the correlation matrix of the global component that was reconstructed without implementing the magnitude difference. In the presence of time lag, even when the global component is composed of a single IC, the spatial difference of its magnitude can contribute to the characterization of the correlation matrix, let alone the global component composed of multiple ICs. It is therefore important to eliminate the effect to determine if time lag is the source of synchronization. The magnitude of each IC was adjusted to reflect the contribution of each IC that is computed as the root mean square of the mixing matrix of the temporal ICA.
In order to further confirm the relationship between the signal synchronization and the time lag of the global component, for each global component, Pearson’s correlation coefficient was computed between the FC matrix and the matrix of relative time lag that is the difference of the time lags between given ROIs.
Contribution of the local component was also assessed by computing the correlation between the FC matrix and the correlation matrix of the local component in the same way. Correlation between the correlation matrices were computed using the upper triangle of each correlation matrix. The threshold of the statistical significance was set at p = 0.05, and the Bonferroni correction was used to control for the multiple comparisons.
Origin of the Time Lag
To estimate the origin of the time lag that characterizes the global component, all magnitude (correlation coefficient) and time lag maps of the global ICs that were averaged across subjects were compared with those of the whole-brain signal, respectively. Partial correlation analysis was also performed to control the effect of vascular time lag that was measured using dynamic susceptibility contrast enhanced perfusion imaging (Amemiya et al., 2016). We also compared the time lag structure of the local component and that of the whole-brain signal by computing the Pearson’s correlation between the time lag maps with each local IC. A correction for the spatial degrees of freedom was given via Gaussian random field theory and empirical smoothness estimation, which estimated the number of independent resels or resolution elements to be 103.
Results
Identification of the Global Waves
Of the 28 and 30 reproducible ICs for dataset 1 and 2, 7, and 10 were classified as global IC based on the pattern of spatial distribution for each dataset, respectively (Figure 2 and Supplementary Figure S2). The magnitude of the global IC is shown as Pearson’s correlation coefficient between the time-shifted global IC and the time series of each voxel, with the corresponding time lag structures showing the paths of each global component (Figure 2 and Supplementary Figure S2). All magnitude and time lag maps shown were obtained by averaging the resulting maps across all subjects’ all runs. Time lag maps of the global components showed structural paths of each signal, which is consistent with previous studies suggesting the existence of multiple global waves of activity in the resting state; c25, c07 and c05 show early regions in the rostral and lateral part of the frontal lobes and delayed regions in the medial part of the frontal lobes, insular and inferior frontal gyrus and occipital lobes, while the pattern is almost opposite for c10. C12’s path is characterized by early regions in the sensorimotor, auditory, and visual cortex, as well as delayed regions in the association cortex and posterior cingulate cortex, while c07 shows the opposite pattern. C15 resembles c05, c07, c25 pattern, but the delay in the dorsal attention network is more conspicuous (Figure 2).
FIGURE 2
However, there were significant correlations among the paths of global signals (Supplementary Figures S6A,B). Some global ICs also showed apparent similarity to the whole-brain signal not only in its spatial distribution characterized by symmetrical involvement of the dorsal cerebral cortex with predominantly high magnitude in the occipital lobes, but also in its propagation pattern: time lag: |r| = 0.42 ± 0.21 (re-test, 0.37 ± 0.23); magnitude: r = 0.58 ± 0.15 (re-test, 0.55 ± 0.12) (Figure 2 and Supplementary Figure S2). Partial correlation analysis controlled for the perfusion time lag also showed sometimes reduced but still significant correlation between the time lag maps of the global ICs and the whole-brain signal: |r| = 0.40 ± 0.23 (re-test, 0.39 ± 0.22). These results suggest that even if multiple global waves of activity coexist in the resting state, the paths of the traveling waves can be substantially overlapped, as can the spatial distribution, which is not simply explained as the result of common background vascular perfusion.
Time lag maps of the global ICs detected from the test dataset were well replicated by the analysis of the re-test dataset. Supplementary Figure S6C demonstrates that all global ICs of the test data were significantly positively correlated with at least one global IC of the re-test data.
Consistent with previous studies (Mitra et al., 2014; Amemiya et al., 2016), magnitude and latency of the whole-brain signal were not significantly correlated: r = −0.10, p = 0.32 (re-test, r = −0.051, p = 0.62). The majority of the global components showed significant correlation between magnitude and time lag (p < 0.05): c07, r = 0.41; c10, r = −0.30; c12, r = −0.48; c12, r = −0.48; c25, r = −0.27 (re-test, c05, r = −0.29; c07, −0.44; c12, r = −0.31; c18, r = −0.36; c22, r = −0.24; c24, r = −0.25; c26, r = −0.42), which might be caused by the attenuation of the waves of activity during the process of transmission.
Some of the global ICs showed anteroposterior propagation that might correspond to the pattern detected using electroencephalogram in sleeping humans: c22 of the test dataset as well as c22 and c26 of the re-test dataset show early regions in the rostral compared with caudal part of the cerebral cortex.
Contribution of the Global Waves to the RSN Synchronization
The correlation matrix of the global component reconstructed with the detected global ICs showed significant correlation with the FC matrix: r = 0.31 ± 0.13 (re-test, 0.32 ± 0.11), p < 0.001 (Figures 3A,B and Supplementary Figures S3A,B). Significant correlation was found even when the global component was reconstructed without considering the spatial difference of its magnitude: r = 0.22 ± 0.13 (re-test, 0.27 ± 0.11), p < 0.001 (Figure 3C and Supplementary Figure S3C). Furthermore, significant negative correlation between the strength of synchronization (FC) and the relative time lag was also shown for all global waves: r = −0.19 to −0.39, p < 0.001 (re-test, r = −0.14 to −0.40, p < 0.001) (Figures 3E, 4 and Supplementary Figures S3E, S4), which suggests that the time lag of the global component can contribute to the RSN synchronization.
FIGURE 3
FIGURE 4
Characteristics of the Local Component
The spatial distribution of the magnitude of the 28 and 30 reproducible local ICs is shown in Figure 5 and Supplementary Figure S5, respectively. Each magnitude map of the local ICs showed significant synchronization within functionally relevant structures, which would correspond to spatial maps for the temporally independent functional modes (Smith et al., 2012). In other words, as previously well-explored in Smith et al. (2012), the local ICs could be considered as functional “modes” that in some cases could subdivide and/or reorganize the currently standard spatial RSNs.; e.g., c01-03 contain visual cortex (predominantly extrastriate areas) and ventral sensorimotor cortex; c04 contains right orbitofrontal cortex in addition to sensorimotor and visual cortex; c06 contains extrastriate cortex and basal ganglia; c08 contains sensorimotor cortex, c09, c17, c18 involve frontotemporal network nodes; c11 mainly involves association cortex; c13 and contains primary visual cortex and default mode network nodes in the angular gyrus and posterior cingulate cortex; c14 contains visual cortex and dorsolateral prefrontal cortex; c16 involves relatively widespread areas mainly involving the sensorimotor and primary visual cortex in addition to basal ganglia; c19 contains; c20 and c21 contain dorsal attention network and frontal lobes; c23 involves primary visual cortex and frontal lobes; c24 and c28 contains salience network and frontoparietal network nodes; c26 contains auditory and sensorimotor networks; c27 involves sensorimotor and visual cortex. Figure 3D and Supplementary Figure S3D show correlation matrices of the local component that are fairly similar to the FC matrix (r = 0.41, p < 0.05; re-test, r = 0.42, p < 0.05), suggesting an equivalent or even larger contribution of the local activity to the RSN synchronization that characterizes the FC matrix. The magnitude of the activity of the local component relative to the whole signal was 0.70 for both test and re-test dataset.
FIGURE 5
Discussion
By applying temporal ICA to the rs-fMRI data, we have identified several global or semi-global waves of slow oscillation that are temporally independent yet spatially overlapping with each other. Although the correlation matrix of the global component showed substantial correlation with the FC matrix, an equivalent or even greater contribution of the local component was also shown. The results indicate that while global waves of activity, although being less spatially homogenous than what was assumed to be, could contribute to the emergence of the RSN, which is partly consistent with previous studies suggesting that within-network synchronization can arise from the time lag of the global waves (Mitra et al., 2015a; Matsui et al., 2016), this does not exclude the contribution of local activity that are more confined within functionally relevant structures.
Multiple Waves of Activity
Although the number of the global IC was slightly varied depending on the dataset, there was substantial overlap among the global ICs detected within or across the temporal ICAs for test and re-test datasets. Moreover, the majority of the global waves detected across the temporal ICAs were significantly correlated with the path of the whole-brain (global mean) signal, which would be the most robust representation of the global signal. Partial correlation analysis that controlled for the effect of vascular time lag also confirmed that significant correlation was still found with the global waves. While these results might support the finding that the global neural activity has a predominant path of propagation (Matsui et al., 2016), they also suggest the existence of multiple overlapping paths of neural oscillations. Some paths of global ICs that share common features with those obtained in the previous studies: c25, c05, c07, c10, and c12 would correspond to the thread 2 reported in Mitra et al. (2015a) that shows a contrast between the rostral and lateral part of the frontal lobes vs. medial part of the frontal lobes, insular and inferior frontal gyrus and occipital lobes, while c15 to thread 8 (note that the polarity of the threads can be inverted). However, there were some other time lag structures representing the paths of propagating neural activity that were not obtained by merely decomposing the measured total time lag as multiple orthogonal components (i.e., threads in Mitra et al., 2015a) or as independent components (Amemiya et al., 2016). Specifically, some global ICs showed anteroposterior propagation that might correspond to the pattern detected using electroencephalogram in sleeping humans (Massimini et al., 2004) or with calcium imaging as well as BOLD imaging in anesthetized mice (Stroh et al., 2013; Matsui et al., 2016). While physiological basis or significance of such global activity remains to be known, all these data further support the view that spatiotemporal pattern of BOLD signal could reflect large-scale dynamics of underlying neuronal activity.
Origin of the Time Lag and Synchronization
Given that cerebral vascular time lag is quite uniform across subjects (Amemiya et al., 2016), existence of multiple paths of traveling BOLD signal suggest the existence of multiple waves of neural signal (Mitra et al., 2015a; Amemiya et al., 2016). Although BOLD represents hemodynamic response to neural activity that is necessarily influenced by the characteristics of the underlying vasculature (Amemiya et al., 2012; Bandettini, 2014), the path of globally propagating activity has been shown to coincide with that of the neuronal calcium signal in mice (Matsui et al., 2016). Assuming that the same holds true for non-anesthetized awake human data, the propagating pattern of global oscillations that are characterized by structured and smooth gradation can be seen as corresponding to a gradual propagation via short-range corticocortical connections. In addition, small time lag observed between distant regions across RSNs may reflect the presence of a mechanism controlling the initiation of spreading activity, mediated via long-range connections in a rapid manner. Such activities might help integrate the spontaneous oscillation of the cortex across RSNs in the whole brain, which is analogous to the concept of synfire chains in synchronous mode (Abeles, 1982, 1991), in which groups of neurons are organized into chains, and the architecture enables precisely timed sequences of spikes to form a propagating wave of activity (Abeles, 1982, 1991; Diesmann et al., 1999).
Alternatively, it would also be possible to assume that the time lag of the global oscillations mainly reflects vascular dynamics for some ICs. Indeed, whole-brain signal and vascular perfusion are known to share similar spatiotemporal characteristics (Amemiya et al., 2014, 2016; Tong et al., 2017). It might be important to note, however, that the source of the time lag is not necessarily identical to the source of the signal, so even if the time lag were totally non-neural in origin, that does not mean that the origin of the global signal is non-neural. This is because perfusion time lag can also be reflected in the time lag of BOLD signal of neural origin (Roc et al., 2006; Amemiya et al., 2012). Therefore, for the ICs whose time lag maps are similar to that of perfusion, e.g., like whole-brain signal, a measurement that is independent of neurovascular coupling would be preferable and perhaps essential for a more precise prediction of the spatiotemporal profile of the underlying neural activity.
Nevertheless, even if the contribution of some global oscillations to the apparent network synchronization were an artifact (Tong et al., 2015), the results of the present study suggest that network-specific synchronization does exist besides such component, which is consistent with the growing evidence supporting the link between BOLD and electroencephalographic or magnetoencephalographic measures of resting state activity (Goldman et al., 2002; Yan et al., 2009; Brookes et al., 2011; Tagliazucchi et al., 2012). Moreover, the present study indicates that within-network synchronization is dependent on local neural activities that are temporally independent to each other, which necessitates the presence of a mechanism that would be conceived as network-specific pacemaker irrespective of the contribution of the global oscillations.
Technical Issues
In the present study, two large HCP datasets acquired from the same 50 subjects were used to test the reproducibility of the analysis. In temporal ICA, each component’s independence is optimized for the axis of time. Therefore, the temporal dimensionality or timepoints of the data should be large enough. Although there is no good reason to assume that the number of timepoints should be as large as the number of voxels in original data, when the dimension was reduced during the process of group spatial ICA in advance to temporal ICA. Rather the problem lies in that it is generally difficult and practically impossible to know in advance how many timepoints are needed for an ICA, which is particularly dependent on the non-gaussianity of the data. This is why post hoc analysis is generally considered important to validate an ICA.
It is also theoretically apparent that higher spatiotemporal resolution is preferable for a better mapping of the spatiotemporal characteristic of the data. However, identification of the global component is not likely dependent on the spatiotemporal resolution of the BOLD fMRI. This is because ICs are classified according to the pattern of spatial distribution, which is not dependent on the temporal or spatial resolution itself let alone the speed of the traveling waves. Therefore, granting that the time lag maps would become more accurate if the sampling rate or spatial resolution of the data were increased, given that the maps obtained in the present study represent structured and highly similar patterns even across studies for the whole-brain signal (Amemiya et al., 2014; Mitra et al., 2014; Tong et al., 2017), such contribution would be negligible compared with other factors, at least for the range of the neural band of 0.01–0.1 Hz. The same holds true for the slice-time correction. HCP do not recommend us doing slice timing correction for the dataset, because while the effect of the slice timing correction is limited for the short TR (0.72s), slice timing correction interacts with movement correction in ways that have not ever been appropriately addressed in available tools. However, given the fact that it is impossible to align the subjects’ head in exactly the same position for every scan for all the subjects, and that the acquisition was performed by using simultaneous multislice imaging with a slice thickness of 2 mm, we consider that the small slice timing differences were expected to be canceled out during the course of spatial normalization and group averaging of the time lag maps, which was confirmed by high correlation among the whole-brain signal time lag maps.
For the preprocessing, we did not apply global signal regression (GSR). Although GSR is a useful process to remove physiological noise like motion artifact, it eliminates any global signal regardless of the origin and can distort the resulting connectivity or activation measures in a complex way (Saad et al., 2012; Gotts et al., 2013; Glasser et al., 2016, 2018; Tang et al., 2019). Therefore, GSR and related approaches still remain controversial. Given that the study aim is to understand the possible contribution of the global signal, we consider it important not to apply GSR for our analysis.
In the present study, detection of the global IC was based on the spatial distribution pattern of each component that was judged by template matching using RSNs as well as the magnitude map of the whole-brain signal, which enabled us to classify an IC into global or local component without setting a threshold for the spatial coverage of that component. Although intuitively, an IC showing a larger spatial coverage would be considered as a global component, such classification is practically very difficult, because there is no objective definition regarding the coverage of a global component, as was the case in the original definition of the RSNs.
Rather, it is important to note that the present study, like previous studies, focuses on identifying slow waves with fixed patterns of propagation. While such an approach is advantageous in exploring the most robust representation of the phenomenon, an analysis allowing more spatiotemporally complex and dynamically changing patterns of propagation will probably reveal a more precise picture of the inter-network activities that may contribute to the integration of the network-specific activities constituting the functional architecture of the brain.
Statements
Data availability statement
The datasets for this study are openly available at https://www.humanconnectome.org/.
Author contributions
SA contributed to the conception and design of the study, performed the statistical analysis, and wrote the first draft of the manuscript. SA and HT contributed to the interpretation of the data. All authors contributed to the manuscript revision, and read and approved the submitted version.
Funding
This work was supported by JSPS Grant-in-Aid for Scientific Research (C 18K07707), and grants from the Takeda Science Foundation, the Tokyo Society of Medical Sciences, and the Japanese Society of Neuroradiology to SA.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnsys.2019.00065/full#supplementary-material
References
1
AbbasA.BelloyM.KashyapA.BillingsJ.NezafatiM.SchumacherE. H.et al (2019). Quasi-periodic patterns contribute to functional connectivity in the brain.Neuroimage191193–204. 10.1016/j.neuroimage.2019.01.076
2
AbelesM. (1982). Local Cortical Circuits: an Electrophysiological Study.Berlin: Springer.
3
AbelesM. (1991). Corticonics: Neural Circuits of the Cerebral Cortex.New York, NY: Cambridge University Press.
4
AmemiyaS.KunimatsuA.SaitoN.OhtomoK. (2012). Impaired hemodynamic response in the ischemic brain assessed with BOLD fMRI.Neuroimage61579–590. 10.1016/j.neuroimage.2012.04.001
5
AmemiyaS.KunimatsuA.SaitoN.OhtomoK. (2014). Cerebral hemodynamic impairment: assessment with resting-state functional MR imaging.Radiology270548–555. 10.1148/radiol.13130982
6
AmemiyaS.TakaoH.HanaokaS.OhtomoK. (2016). Global and structured waves of rs-fMRI signal identified as putative propagation of spontaneous neural activity.Neuroimage133331–340. 10.1016/j.neuroimage.2016.03.033
7
BandettiniP. A. (2014). Neuronal or hemodynamic? Grappling with the functional MRI signal.Brain Connect.4487–498. 10.1089/brain.2014.0288
8
BeckmannC. F.SmithS. M. (2004). Probabilistic independent component analysis for functional magnetic resonance imaging.IEEE Trans. Med. Imaging23137–152. 10.1109/tmi.2003.822821
9
BelloyM. E.NaeyaertM.AbbasA.ShahD.VanreuselV.van AudekerkeJ.et al (2018). Dynamic resting state fMRI analysis in mice reveals a set of Quasi-periodic patterns and illustrates their relationship with the global signal.Neuroimage180(Pt B), 463–484. 10.1016/j.neuroimage.2018.01.075
10
BiswalB.YetkinF. Z.HaughtonV. M.HydeJ. S. (1995). Functional connectivity in the motor cortex of resting human brain using echo-planar MRI.Magn. Reson. Med.34537–541. 10.1002/mrm.1910340409
11
BrookesM. J.WoolrichM.LuckhooH.PriceD.HaleJ. R.StephensonM. C.et al (2011). Investigating the electrophysiological basis of resting state networks using magnetoencephalography.Proc. Natl. Acad. Sci. U.S.A.10816783–16788. 10.1073/pnas.1112685108
12
CalabreseR. L. (1998). Cellular, synaptic, network, and modulatory mechanisms involved in rhythm generation.Curr. Opin. Neurobiol.8710–717. 10.1016/s0959-4388(98)80112-8
13
CivillicoE. F.ContrerasD. (2012). Spatiotemporal properties of sensory responses in vivo are strongly dependent on network context.Front. Syst. Neurosci.6:25. 10.3389/fnsys.2012.00025
14
DiesmannM.GewaltigM. O.AertsenA. (1999). Stable propagation of synchronous spiking in cortical neural networks.Nature402529–533. 10.1038/990101
15
FerezouI.HaissF.GentetL. J.AronoffR.WeberB.PetersenC. C. (2007). Spatiotemporal dynamics of cortical sensorimotor integration in behaving mice.Neuron56907–923. 10.1016/j.neuron.2007.10.007
16
FilippiniN.MacIntoshB. J.HoughM. G.GoodwinG. M.FrisoniG. B.SmithS. M.et al (2009). Distinct patterns of brain activity in young carriers of the APOE-epsilon4 allele.Proc. Natl. Acad. Sci. U.S.A.1067209–7214. 10.1073/pnas.0811879106
17
FoxM. D.RaichleM. E. (2007). Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging.Nat. Rev. Neurosci.8700–711. 10.1038/nrn2201
18
FoxM. D.SnyderA. Z.VincentJ. L.CorbettaM.Van EssenD. C.RaichleM. E. (2005). The human brain is intrinsically organized into dynamic, anticorrelated functional networks.Proc. Natl. Acad. Sci. U.S.A.1029673–9678. 10.1073/pnas.0504136102
19
GlasserM. F.CoalsonT. S.BijsterboschJ. D.HarrisonS. J.HarmsM. P.AnticevicA.et al (2018). Using temporal ICA to selectively remove global noise while preserving global signal in functional MRI data.Neuroimage181692–717. 10.1016/j.neuroimage.2018.04.076
20
GlasserM. F.SmithS. M.MarcusD. S.AnderssonJ. L.AuerbachE. J.BehrensT. E.et al (2016). The human connectome project’s neuroimaging approach.Nat. Neurosci.191175–1187. 10.1038/nn.4361
21
GlasserM. F.SotiropoulosS. N.WilsonJ. A.CoalsonT. S.FischlB.AnderssonJ. L.et al (2013). The minimal preprocessing pipelines for the human connectome project.Neuroimage80105–124. 10.1016/j.neuroimage.2013.04.127
22
GoldmanR. I.SternJ. M.EngelJ.Jr.CohenM. S. (2002). Simultaneous EEG and fMRI of the alpha rhythm.Neuroreport132487–2492. 10.1097/01.wnr.0000047685.08940.d0
23
GottsS. J.SaadZ. S.JoH. J.WallaceG. L.CoxR. W.MartinA. (2013). The perils of global signal regression for group comparisons: a case study of autism spectrum disorders.Front. Hum. Neurosci.7:356. 10.3389/fnhum.2013.00356
24
GriffantiL.Salimi-KhorshidiG.BeckmannC. F.AuerbachE. J.DouaudG.SextonC. E.et al (2014). ICA-based artefact removal and accelerated fMRI acquisition for improved resting state network imaging.Neuroimage95232–247. 10.1016/j.neuroimage.2014.03.034
25
Harris-WarrickR. M. (2010). General principles of rhythmogenesis in central pattern generator networks.Prog. Brain Res.187213–222. 10.1016/B978-0-444-53613-6.00014-9
26
HimbergJ.HyvarinenA.EspositoF. (2004). Validating the independent components of neuroimaging time series via clustering and visualization.Neuroimage221214–1222. 10.1016/j.neuroimage.2004.03.027
27
HyvärinenA. (1999). Fast and robust fixed-point algorithms for independent component analysis.IEEE Trans. Neural Netw.10626–634. 10.1109/72.761722
28
IvanchenkoM. V.ThomasN.SelverstonA. I.RabinovichM. I. (2008). Pacemaker and network mechanisms of rhythm generation: cooperation and competition.J. Theor. Biol.253452–461. 10.1016/j.jtbi.2008.04.016
29
JohnsonS. M.WiegelL. M.MajewskiD. J. (2007). Are pacemaker properties required for respiratory rhythm generation in adult turtle brain stems in vitro?Am. J. Physiol. Regul. Integr. Comp. Physiol.293R901–R910. 10.1152/ajpregu.00912.2006
30
KenetT.BibitchkovD.TsodyksM.GrinvaldA.ArieliA. (2003). Spontaneously emerging cortical representations of visual attributes.Nature425954–956. 10.1038/nature02078
31
LamplI.ReichovaI.FersterD. (1999). Synchronous membrane potential fluctuations in neurons of the cat visual cortex.Neuron22361–374. 10.1016/s0896-6273(00)81096-x
32
MajeedW.MagnusonM.HasenkampW.SchwarbH.SchumacherE. H.BarsalouL.et al (2011). Spatiotemporal dynamics of low frequency BOLD fluctuations in rats and humans.Neuroimage541140–1150. 10.1016/j.neuroimage.2010.08.030
33
MajeedW.MagnusonM.KeilholzS. D. (2009). Spatiotemporal dynamics of low frequency fluctuations in BOLD fMRI of the rat.J. Magn. Reson. Imaging30384–393. 10.1002/jmri.21848
34
MassiminiM.HuberR.FerrarelliF.HillS.TononiG. (2004). The sleep slow oscillation as a traveling wave.J. Neurosci.246862–6870. 10.1523/jneurosci.1318-04.2004
35
MatsuiT.MurakamiT.OhkiK. (2016). Transient neuronal coactivations embedded in globally propagating waves underlie resting-state functional connectivity.Proc. Natl. Acad. Sci. U.S.A.1136556–6561. 10.1073/pnas.1521299113
36
MeijeringE. (2002). A chronology of interpolation: from ancient astronomy to modern signal and image processing.Proc. IEEE90319–342. 10.1109/5.993400
37
MitraA.SnyderA. Z.BlazeyT.RaichleM. E. (2015a). Lag threads organize the brain’s intrinsic activity.Proc. Natl. Acad. Sci. U.S.A.112E2235–E2244. 10.1073/pnas.1503960112
38
MitraA.SnyderA. Z.TagliazucchiE.LaufsH.RaichleM. E. (2015b). Propagated infra-slow intrinsic brain activity reorganizes across wake and slow wave sleep.eLife4:e10781. 10.7554/eLife.10781
39
MitraA.SnyderA. Z.HackerC. D.RaichleM. E. (2014). Lag structure in resting-state fMRI.J. Neurophysiol.1112374–2391. 10.1152/jn.00804.2013
40
MohajeraniM. H.ChanA. W.MohsenvandM.LeDueJ.LiuR.McVeaD. A.et al (2013). Spontaneous cortical activity alternates between motifs defined by regional axonal projections.Nat. Neurosci.161426–1435. 10.1038/nn.3499
41
MullerL.DestexheA. (2012). Propagating waves in thalamus, cortex and the thalamocortical system: experiments and models.J. Physiol. Paris106222–238. 10.1016/j.jphysparis.2012.06.005
42
PetersenC. C.HahnT. T.MehtaM.GrinvaldA.SakmannB. (2003). Interaction of sensory responses with spontaneous depolarization in layer 2/3 barrel cortex.Proc. Natl. Acad. Sci. U.S.A.10013638–13643. 10.1073/pnas.2235811100
43
Ponce-AlvarezA.DecoG.HagmannP.RomaniG. L.MantiniD.CorbettaM. (2015). Resting-state temporal synchronization networks emerge from connectivity topology and heterogeneity.PLoS Comput. Biol.11:e1004100. 10.1371/journal.pcbi.1004100
44
PruimR. H.MennesM.van RooijD.LleraA.BuitelaarJ. K.BeckmannC. F. (2015). ICA-AROMA: a robust ICA-based strategy for removing motion artifacts from fMRI data.Neuroimage112267–277. 10.1016/j.neuroimage.2015.02.064
45
RocA. C.WangJ.AncesB. M.LiebeskindD. S.KasnerS. E.DetreJ. A. (2006). Altered hemodynamics and regional cerebral blood flow in patients with hemodynamically significant stenoses.Stroke37382–387. 10.1161/01.str.0000198807.31299.43
46
RybakI. A.AbdalaA. P.MarkinS. N.PatonJ. F.SmithJ. C. (2007). Spatial organization and state-dependent mechanisms for respiratory rhythm and pattern generation.Prog. Brain Res.165201–220. 10.1016/s0079-6123(06)65013-9
47
RybakI. A.ShevtsovaN. A.PatonJ. F.PierreficheO.St-JohnW. M.HajiA. (2004). Modelling respiratory rhythmogenesis: focus on phase switching mechanisms.Adv. Exp. Med. Biol.551189–194. 10.1007/0-387-27023-x_29
48
SaadZ. S.GottsS. J.MurphyK.ChenG.JoH. J.MartinA.et al (2012). Trouble at rest: how correlation patterns and group differences become distorted after global signal regression.Brain Connect.225–32. 10.1089/brain.2012.0080
49
Salimi-KhorshidiG.DouaudG.BeckmannC. F.GlasserM. F.GriffantiL.SmithS. M.et al (2014). Automatic denoising of functional MRI data: combining independent component analysis and hierarchical fusion of classifiers.Neuroimage90449–468. 10.1016/j.neuroimage.2013.11.046
50
ShuY.HasenstaubA.McCormickD. A. (2003). Turning on and off recurrent balanced cortical activity.Nature423288–293. 10.1038/nature01616
51
SmithS. M.BeckmannC. F.AnderssonJ.AuerbachE. J.BijsterboschJ.DouaudG.et al (2013). Resting-state fMRI in the human connectome project.Neuroimage80144–168. 10.1016/j.neuroimage.2013.05.039
52
SmithS. M.MillerK. L.MoellerS.XuJ.AuerbachE. J.WoolrichM. W.et al (2012). Temporally-independent functional modes of spontaneous brain activity.Proc. Natl. Acad. Sci. U.S.A.1093131–3136. 10.1073/pnas.1121329109
53
SohalV. S.Pangratz-FuehrerS.RudolphU.HuguenardJ. R. (2006). Intrinsic and synaptic dynamics interact to generate emergent patterns of rhythmic bursting in thalamocortical neurons.J. Neurosci.264247–4255. 10.1523/jneurosci.3812-05.2006
54
SteriadeM.McCormickD. A.SejnowskiT. J. (1993). Thalamocortical oscillations in the sleeping and aroused brain.Science262679–685. 10.1126/science.8235588
55
StrohA.AdelsbergerH.GrohA.RühlmannC.FischerS.SchierlohA.et al (2013). Making waves: initiation and propagation of corticothalamic Ca2+ waves in vivo.Neuron771136–1150. 10.1016/j.neuron.2013.01.031
56
TagliazucchiE.von WegnerF.MorzelewskiA.BrodbeckV.LaufsH. (2012). Dynamic BOLD functional connectivity in humans and its electrophysiological correlates.Front. Hum. Neurosci.6:339. 10.3389/fnhum.2012.00339
57
TakedaY.HiroeN.YamashitaO.SatoM. A. (2016). Estimating repetitive spatiotemporal patterns from resting-state brain activity data.Neuroimage133251–265. 10.1016/j.neuroimage.2016.03.014
58
TangY.ZhouQ.ChangM.ChekroudA.GueorguievaR.JiangX.et al (2019). Altered functional connectivity and low-frequency signal fluctuations in early psychosis and genetic high risk.Schizophr Res.210172–179. 10.1016/j.schres.2018.12.041
59
TongY.HockeL. M.FanX.JanesA. C.FrederickB. (2015). Can apparent resting state connectivity arise from systemic fluctuations?Front. Hum. Neurosci.9:285. 10.3389/fnhum.2015.00285
60
TongY.LindseyK. P.HockeL. M.VitalianoG.MintzopoulosD.FrederickB. D. (2017). Perfusion information extracted from resting state functional magnetic resonance imaging.J. Cereb. Blood Flow Metab.37564–576. 10.1177/0271678x16631755
61
TsodyksM.KenetT.GrinvaldA.ArieliA. (1999). Linking spontaneous activity of single cortical neurons and the underlying functional architecture.Science2861943–1946. 10.1126/science.286.5446.1943
62
Van EssenD. C.SmithS. M.BarchD. M.BehrensT. E.YacoubE.UgurbilK.et al (2013). The WU-minn human connectome project: an overview.Neuroimage8062–79. 10.1016/j.neuroimage.2013.05.041
63
VincentJ. L.PatelG. H.FoxM. D.SnyderA. Z.BakerJ. T.Van EssenD. C.et al (2007). Intrinsic functional architecture in the anaesthetized monkey brain.Nature44783–86. 10.1038/nature05758
64
WangX. J. (2010). Neurophysiological and computational principles of cortical rhythms in cognition.Physiol. Rev.901195–1268. 10.1152/physrev.00035.2008
65
Whitfield-GabrieliS.Nieto-CastanonA. (2012). Conn: a functional connectivity toolbox for correlated and anticorrelated brain networks.Brain Connect.2125–141. 10.1089/brain.2012.0073
66
XuW.HuangX.TakagakiK.WuJ.-Y. (2007). Compression and reflection of visually evoked cortical waves.Neuron55119–129. 10.1016/j.neuron.2007.06.016
67
YanL.ZhuoY.YeY.XieS. X.AnJ.AguirreG. K.et al (2009). Physiological origin of low-frequency drift in blood oxygen level dependent (BOLD) functional magnetic resonance imaging (fMRI).Magn. Reson. Med.61819–827. 10.1002/mrm.21902
68
ZhangD.SnyderA. Z.ShimonyJ. S.FoxM. D.RaichleM. E. (2010). Noninvasive functional and structural connectivity mapping of the human thalamocortical system.Cereb. Cortex201187–1194. 10.1093/cercor/bhp182
Summary
Keywords
fMRI, resting-state network, spatiotemporal dynamics, spontaneous neural activity, neuronal pathway tracing
Citation
Amemiya S, Takao H and Abe O (2019) Global vs. Network-Specific Regulations as the Source of Intrinsic Coactivations in Resting-State Networks. Front. Syst. Neurosci. 13:65. doi: 10.3389/fnsys.2019.00065
Received
15 August 2019
Accepted
14 October 2019
Published
29 October 2019
Volume
13 - 2019
Edited by
Nicola Toschi, University of Rome Tor Vergata, Italy
Reviewed by
Ines R. Violante, University of Surrey, United Kingdom; Janine Dale Mendola, McGill University, Canada
Updates
Copyright
© 2019 Amemiya, Takao and Abe.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Shiori Amemiya, amemiya-tky@umin.ac.jp
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.