ORIGINAL RESEARCH article
A New Analysis of Resting State Connectivity and Graph Theory Reveals Distinctive Short-Term Modulations due to Whisker Stimulation in Rats
- 1Institute of Experimental and Clinical Pharmacology and Toxicology, FAU Erlangen-Nuremberg, Erlangen, Germany
- 2Department of Radiology, University Hospital Erlangen, FAU Erlangen-Nuremberg, Erlangen, Germany
- 3Faculty of Mathematical & Physical Sciences, Benemerita Universidad Autonoma de Puebla, Puebla, Mexico
Resting state (RS) connectivity has been increasingly studied in healthy and diseased brains in humans and animals. This paper presents a new method to analyze RS data from fMRI that combines multiple seed correlation analysis with graph-theory (MSRA). We characterize and evaluate this new method in relation to two other graph-theoretical methods and ICA. The graph-theoretical methods calculate cross-correlations of regional average time-courses, one using seed regions of the same size (SRCC) and the other using whole brain structure regions (RCCA). We evaluated the reproducibility, power, and capacity of these methods to characterize short-term RS modulation to unilateral physiological whisker stimulation in rats. Graph-theoretical networks found with the MSRA approach were highly reproducible, and their communities showed large overlaps with ICA components. Additionally, MSRA was the only one of all tested methods that had the power to detect significant RS modulations induced by whisker stimulation that are controlled by family-wise error rate (FWE). Compared to the reduced resting state network connectivity during task performance, these modulations implied decreased connectivity strength in the bilateral sensorimotor and entorhinal cortex. Additionally, the contralateral ventromedial thalamus (part of the barrel field related lemniscal pathway) and the hypothalamus showed reduced connectivity. Enhanced connectivity was observed in the amygdala, especially the contralateral basolateral amygdala (involved in emotional learning processes). In conclusion, MSRA is a powerful analytical approach that can reliably detect tiny modulations of RS connectivity. It shows a great promise as a method for studying RS dynamics in healthy and pathological conditions.
Since Biswal et al. (1995) first described intrinsic functional connectivity in the human brain during rest by functional MRI (fMRI), the so called resting state (RS) connectivity has been studied intensively and several large scale neural networks have been detected (see Raichle, 2011 for review). Despite the influence of individual and instantaneous factors such as mood, physiological and cognitive states, RS networks are remarkably robust and stable over time (Braun et al., 2012; Zuo and Xing, 2014; Du et al., 2015) and species (Lu et al., 2012; Sierakowiak et al., 2015; Gozzi and Schwarz, 2016). Thus, they seem to exhibit an evolutionary conserved and fundamental phenomenon of mammalian brain function. However, its biological relevance is still not fully understood. One hypothesis interprets human RS connectivity as a correlate of introspective mental processes (including such processes as mind wandering) that influences behavioral responses to future events (Rosazza and Minati, 2011). Other researchers emphasize the relation of RS topography and strength to the history of network activation and thus support the Hebbian-like Fire-Wire-hypothesis: regions that are co-activated during task performance develop stronger coherence at rest (Corbetta, 2012; Harmelech et al., 2013). Regardless of other interpretations, it has been widely accepted that resting state networks are dynamic in nature (Deco and Corbetta, 2011). They are modulated by prior task activation, which supports the hypothesis that RS functional connectivity plays a role in learning processes and memory consolidation (Albert et al., 2009; Tambini et al., 2010).
Short-term modulations of RS connectivity have also been detected in several human studies after performance of a motor task (Tung et al., 2013; Mary et al., 2016; Muraskin et al., 2016), but to the best of our knowledge, studies using sensory stimuli as RS modulators or animal models have not been published. One reason might be that we lack adequate analytical tools to detect these tiny modulations. Therefore, our goal is to develop a method that is capable of examining RS connectivity during sensorimotor stimulation of the barrel field in the rat. For that purpose we evaluated three graph-theoretical resting state analysis methods in comparison to ICA to determine which analytical methods can best detect significant and reliable short-term resting state modulations. The barrel field model and its underlying functional and anatomical pathways have been intensively investigated with various neurobiological techniques including fMRI (Grinvald et al., 1986; Yang et al., 1997; Hess et al., 2000; Sachdev et al., 2003; Diamond et al., 2008). We used unilateral stimulation of a small set of whiskers after trimming the remaining ones. Unilateral whisker activation usually indicates food or foe; thus the animal reacts with an immediate behavioral and emotional response (Marshall et al., 1971; Prescott, 2011). When compared to the more natural stimulation of all whiskers, trimming and stimulation of the remaining whiskers induces altered functional activation patterns that are shown to be related to plasticity and learning processes (Mirabella et al., 2001; Sellien and Ebner, 2007; Alonso Bde et al., 2008; de Celis Alonso et al., 2012; Albieri et al., 2015). Consequently this sensory stimulation paradigm is expected to modulate resting state connectivity, yet it should still be subtle enough to pose a challenge for the analytical methods under evaluation.
One major weakness in investigating resting state fMRI is the heterogeneity of analytical approaches, especially given that the method used to analyze RS data strongly influences the results (Ma et al., 2007; Long et al., 2008; Rosazza et al., 2012; Cao et al., 2014). Therefore, reproducing and comparing findings across different laboratories and studies is difficult. Classical methods to analyze RS data include: regional homogeneity analysis (Zang et al., 2004), cluster analysis (Van Den Heuvel et al., 2008a; Lee et al., 2012), seed based correlation analysis (SCA) (Greicius et al., 2003) and independent component analysis (ICA) (Calhoun et al., 2001). More recently, graph-theoretical approaches are used, which translate RS data into networks consisting of nodes and edges by cross correlating regional time-courses (Smith et al., 2013). Definition of regions used in the time-course analysis ranges from voxels (Van Den Heuvel et al., 2008b; Fransson et al., 2011) to only a few anatomical regions (13 nodes) (Fair et al., 2008). The large variance of parcellation resolutions further complicates the ability to compare graph-theoretical findings (Fornito et al., 2010; Wig et al., 2011). Differences in parcellation size influence not only the signal to noise ratio of the mean time-course, but larger parcellation regions are also more likely to integrate more functionally distinct parts.
To access the advantages of sophisticated graph-theoretical analysis recently proposed approaches estimate within-component network graphs for human RS networks (Park et al., 2014; Ribeiro de Paula et al., 2017). Ribeiro de Paula et al. (2017) directly translate spatial ICA components into within-component graphs. Here, edge weights do not represent correlation of time courses (as measure for functional connectivity) but statistical t-values resulting from a linear regression with the ICA component time course as predictor.
As illustrated in Figure 1, RS analysis methods can be classified by their empirical regime (data or hypothesis driven procedures) and their integrated coherence (Long et al., 2008), which is the spatial similarity of all correlating time-courses. Hypothesis driven methods usually involve the definitions of regions of interest (ROIs). These ROIs can be defined either anatomically (e.g., atlas regions) or functionally (e.g., the activated regions of event or stimulus related fMRI). Several reviews provide a comprehensive overview of all methods including the discussion of their advantages and disadvantages (Margulies et al., 2010; Lee et al., 2013).
Figure 1. Prominent categories of resting-state fMRI data analysis classified according to their empirical regime and integrated coherence. Data driven methods ranging from intra-regional to long range inter-regional correlations are ReHo, regional homogeneity; Clus, Cluster analysis; ICA, Independent component analysis; GTA vx, voxelwise graph-theoretical analysis; GTA fuP, Graph theoretical analysis based on functional parcellation (Eickhoff et al., 2015; Fan et al., 2016). Hypothesis driven methods are SCA, seed correlation analysis; GTA usrP, graph-theoretical analysis with user defined parcellation. For the last two methods, regions can be defined relying either on anatomical (anROIs) or functional (fuROIs) properties.
In this study we introduce a new method that combines classical seed correlation analysis (SCA) with graph-theory. We chose SCA as fundamental part of the new approach because seed regions can be small and equal in size providing comparable signal to noise ratios and minimizing the effects of averaging functionally different voxels. Additionally, the correlating voxels in target brain regions are not predefined but determined data driven. That means the method finds for each subject the optimal correlation within the target region. This effect should diminish the variability and enhance the reliability of functional connectivity weights between two regions in brain networks across subjects. We compared this new method to two other graph-theoretical methods that rely on the same parcellation scheme (i.e. the same digital brain region atlas) as well as to ICA. Since ICA was the first method used to characterize distinct resting state networks in a data driven fashion, we used it to evaluate the communities obtained with the graph-theoretical methods. We evaluated the efficacy of these methods to characterize short-term modulation by comparing RS before and after unilateral whisker stimulation in rats. For this purpose a version of the network based statistics (NBS) (Zalesky et al., 2010) was adapted to cope with pairwise statistics.
Materials and Methods
Permission for animal experiments was obtained from the ethics committee of the government of Mittelfranken (Ansbach, Germany, 621-2531.31-30/00). fMRI experiments were performed on rats (male Sprague Dawley, 350–450 g, Janvier, France, n = 25). The animals were initially anesthetized with 5% isoflurane for 7 min in a 1:1 mixture of oxygen and pressured air. Immediately after, rat whiskers were trimmed on both sides of the snout with exception of the ones in the C row. Rats were then mounted on a Plexiglas cradle with incorporated mouth mask and with a tooth biting bar where the rat's head was fixed without the need of ear screws (see de Celis Alonso et al., 2012 for details on experimental setup). Lateral openings in the mask allowed the whiskers to move freely.
Afterwards the animals were placed in the scanner and the isoflurane concentration was reduced to ≈1.2%. In order to control the depth of anesthesia the isoflurane concentration was adjusted during the experiment to maintain a respiration rate of approximately 65/min (~38 mmHg ± 10% pCO2). A stable spontaneous respiration rate leads to stable transcutaneous pCO2 during the fMRI measurement (Ramos-Cabrer et al., 2005). Body temperature was maintained at 37°C by a circulating water bath, and physiological parameters (respiration, temperature) were monitored (for more details see Hess et al., 2007, 2011).
Each animal underwent one fMRI-session starting with a resting state (RS) measurement followed by a stimulation experiment (either sham or whisker stimulation) and a second RS measurement and a final anatomical scan. Animals were separated into two groups, one with stimulation of the remaining whiskers after trimming between the two RS measurement (experimental group, n = 13) and one control group prepared and mounted in the same way as the experimental group but without whisker trimming or stimulation (control group, n = 12) (Figure 2).
Figure 2. Experimental design including resting state (RS) analysis methods. fMRI images were taken of the animals during three experimental periods per session: RS pre, the stimulation period and RS post. Thirteen animals were whisker stimulated and 12 were in the control group. The pooled animals of both RS pre periods (n = 25) were used to evaluate the RS analysis methods: MSRA (Multiple Seed Region Analysis): time-courses extracted from small seed regions in the center of mass of anatomical brain structures, voxelwise correlated. SRCC (Seed Region Cross Correlation): time-courses extracted from seed regions (see above), cross-correlated. RCCA (Regional Cross Correlation Analysis): time-courses extracted from anatomical brain structures, cross-correlated. ICA: Independent Component Analysis.
Whiskers were stimulated with an air driven device integrated into the holding cradle. An inverted comb situated 2 cm apart from the left side of the snout was used for monolateral stimulations at a frequency of 6 Hz with an amplitude of 10 mm. All the remaining whiskers from the left side of the snout were stimulated. Space between comb teeth was enough to leave some flexibility for the whiskers to slide in and out, but not to get free, thereby avoiding painful pulling stimuli. Combs were driven from an external console connected to the scanner running a custom programmed user interface developed in LabView (Labview, National Instruments, Austin, TX, USA).
MRI experiments were performed on a 4.7 T/40 cm horizontal bore magnet BioSpec (BRUKER, Ettlingen, Germany). A whole-body birdcage resonator enabled homogenous excitation operating with an actively shielded high-power gradient system (200 mT/m) and a low-noise, actively RF-decoupled 2 × 2 phased array head coil (Bruker Biospin, Ettlingen, Germany) was used to acquire brain images. This configuration enables image acquisition with a high temporal signal-to-noise ratio (tSNR) of approximately 60 (Kalthoff et al., 2011).
RS data (collected during both pre and post stimulation period) consisted of 300 brain volumes in each 10 min (cf. Figure 2) acquired with a T2*-weighted single-shot gradient echo-based Echo Planar Imaging sequence (GE-EPI) covering 22 axial slices of the brain in 2 s (TEef = 25.03 ms, TR = 2,000 ms, in-plane resolution 0.391 × 0.391 mm, matrix 64 × 64, FOV 25 × 25 mm, slice thickness 1 mm). Slice 14 was positioned at bregma −3.48 mm according to Paxinos rat brain atlas (Paxinos and Watson, 2007). As a positioning reference we used the smallest distance between the posterior tip of the corpus striatum and the anterior tip of the hippocampus on the horizontal anatomical reference image.
Both sham and whisker stimulation driven fMRI data (between the two RS measurements) was acquired using the same imaging sequence (GE-EPI). A total of 1602 brain volumes were scanned during this experimental period (cf. Figure 2) with the following scanning parameters: TEef = 25.03 ms, TR = 2,000 ms, in-plane resolution 0.391 × 0.391 mm, matrix 64 × 64, FOV of 25 × 25 mm, 1 mm slice thickness. The stimulation protocol included 100 stimulations (6 Hz) of the vibrissae with duration of 8 sec (4 volumes) and intermediate rests of 24 s (12 volumes). First stimulation sequence started after 8 s. Finally, 22 corresponding anatomical T2 reference images (RARE, RF = 8, TEef = 11.7 ms, TR = 3,000 ms, NEX = 5, in-plane resolution 0.097 × 0.097 mm, matrix 256 × 256, field of view 25 × 25 mm, slice thickness 1 mm) were acquired at identical positions.
Analysis of BOLD Activation Due to Whisker Stimulation
In this study we focused on alterations of RS networks as a result of whisker stimulation. Additionally, the brain structures activated by the whisker stimulation were compared to those involved in resting state connectivity (cf. Figure 2).
BOLD activation induced by whisker stimulation was analyzed using standard procedures described by de Celis Alonso et al. (2012). After appropriate preprocessing including inter-slice-time correction, motion correction, and spatial and temporal smoothing, a general linear model (GLM) analysis was performed (for detailed information see Supplementary Methods). The significantly (FDR, q = 0.05) activated voxels in each defined brain structure were determined for each animal. Finally, the number of activated voxels per brain structure was averaged across animals and expressed as percent of total area size of that brain structure.
Preprocessing of RS Data
We used Brainvoyager QX 2.8 (Brain Innovation B.V. Maastricht, The Netherlands) for initial inter-slice time and motion correction of RS data. Inter-slice time correction was calculated in ascending interleaved scan order with cubic spline interpolation, and motion correction was computed using rigid registration and resampling with trilinear/sinc interpolation. In general motion was low for all animals [mean RMS movement 0.054 (±0.029) mm, which equates 4.7 (±2.6) % of voxel space diagonal]. Maximum RMS movement was 0.145 mm (12.7 % of voxel space diagonal). If not stated otherwise, all subsequent analyses were performed with MagnAn (BioCom, Uttenreuth), an IDL application (Exelis Visual Information Solutions Inc., a subsidiary of Harris Corporation, Melbourne, FL, USA) designed for complex image processing and analysis with emphasis on MR imaging.
After 3D Gaussian smoothing (FWHM 1.17 mm), time series were low pass filtered with a frequency of 0.1 Hz and the global signal mean removed by linear regression. The high resolution anatomical reference images of all 25 animals were co-registered to one selected reference subject using an affine linear transformation algorithm with 6 degrees of freedom and averaged. This mean reference was interpolated to the isotropic voxel size of 0.1 × 0.1 × 0.1 mm3 (to fit requirements of ICA software) and served as the anatomical template. Afterwards, the motion corrected functional RS data of each animal were manually skull-stripped, registered to the 3D anatomical template and linearly interpolated to the same isotropic voxel size.
Group ICA was performed by using GIFT software (GIFT v1.3g; icatb.sourceforge.net). Briefly, the whole dataset consisting of the concatenated preprocessed time series of all animals was reduced by means of PCA and subsequently decomposed in a predefined set of group-independent components using the Infomax algorithm. Finally, the corresponding components for each animal were calculated by back reconstruction (Calhoun et al., 2001). We calculated group ICAs with 15, 20, and 30 components and evaluated the resulting components by visual inspection. Decomposition into 20 independent components revealed the best correspondence to previous published rodent ICA components (Becerra et al., 2011; Jonckers et al., 2011). Thus, we focused on ICA with 20 components.
Functional connectivity correlation matrices were analyzed separately on each animal in its native space. Regions of interests corresponding to distinct brain structures must be defined prior to running graph-theoretical analyses. For this purpose the first volume of each RS time series was used as individual anatomical reference to semi-automatically register (affine, 6 degrees of freedom) an in-house digital 3D rat brain atlas with 179 brain structures according to Paxinos and Watson (2007). The definition of the voxels used to extract the time-courses strongly influences the resulting correlation matrices. Therefore, we pooled the animals of the first RS measurement (RS pre) of both groups (n = 25) to compare three different methods to extract the time-courses for the creation of the correlation matrices (cf. Figure 2).
The method introduced in this study, is a multi-seed region approach (MSRA). It relies on multiple seed correlation maps: the mean time-course of a seed region per brain region was correlated with the time-course of every voxel in the brain resulting in one correlation volume per brain region. Seed regions were defined automatically as the 4 voxels nearest to the center of mass of each atlas brain region (5 voxel in total). The following restrictions were taken into account: (i) the voxels must be located within the border of the brain structure; (ii) seed voxels must be located within the central plane of the brain structure as the native datasets are highly anisotropic in the z-direction (slice thickness). Subsequently, significant correlations within that correlation volume were determined using false discovery rate (FDR, q = 0.05). For each seed region the significant correlation values were averaged per atlas brain structure, resulting in a 179 × 179 asymmetric correlation matrix (Figure 3). This asymmetry simulates directionality because for each pair of brain structures, the significant correlation between the seed region and data driven target voxels in the other brain structure is used to generate the correlation matrix. In contrast to the user defined seed regions, target voxels vary within regions and across animals. Note that the correlation itself between seed region and each target voxel is not directed. The resulting asymmetric correlation matrix does not reflect direction in causality, therefore it is called “pseudo directed.”
Figure 3. Schematic demonstration of the MSRA-correlation matrix computation using the example of the seed region in the center of mask of the hind limb field of the primary somatosensory cortex (S1HL). Significant correlations of the seed‘s correlation map were determined using false discovery rate (FDR, q = 0.05, red outline). For each brain area covered by the red outline the significant correlation were averaged and stored in a matrix line. The matrix line values for all other brain areas are set to 0. This procedure was repeated for every seed region resulting in an asymmetric correlation matrix where every line represents a seed region and its average significant correlation with every other brain area.
The second method, the seed region cross-correlation (SRCC) analysis, simply cross correlates the mean time-courses of all seed regions (identical to those used in MSRA) to each other, resulting in one symmetric correlation matrix per subject.
The third method, a regional cross-correlation approach (RCCA), correlates the mean time-courses of the whole brain structure to each of the other brain structures. However, it is biased by the highly varying numbers of voxels used to create the average time-course per brain structure. This can lead to false correlations, especially in inhomogeneous regions. To weaken this effect, we used the first principal component of each region instead of the simple mean to calculate a cross correlation matrix for each subject (Pawela et al., 2008; Zhong et al., 2009). Briefly, principle component analysis (PCA) is a multivariate technique that replaces the measured variables (here the time-courses of the voxels in the region) by a new set of uncorrelated variables (principle components), arranged in order of decreasing variance. The mean variance that could be explained by the first principal component for all animals and structures was 60 ± 16%.
In accordance with the MSRA approach only significant correlations in the SRCC and RCCA correlation matrices (FDR, q = 0.05) were counted as connections. In case of the MSRA method the FDR correction was applied to the whole brain, in case of SRCC and RCCA to the correlation matrix. The resulting FDR thresholds are comparable (mean r of all animals ± standard deviation: MSRA 0.1754 ± 0.004, SRCC 0.1742 ± 0.003, RPCC 0.1723 ± 0.002). Non-significant connections were set to 0. Finally, the resulting Pearson's r correlation values were transformed into Fisher's z-values to provide normal distributions for subsequent statistical analysis.
For each graph-theoretical approach the individual correlation matrices of all subjects were averaged in order to obtain representative group correlation matrices. These average correlation matrices were transformed into network graphs consisting of vertices (i.e., brain regions) and edges (i.e., functional connectivity between pairs of brain regions) (for detailed information see Supplementary Methods). The networks were visualized in Amira 5.4 (Visage Imaging) using a force-based algorithm (Kamada and Kawai, 1989).
For comparison with the ICA approach networks were subdivided into the same number of communities as ICA components (20) using a heuristic method that is based on modularity optimization proposed by Blondel et al. (2008) and implemented in NWB [NWB Team (2006). Network Workbench Tool. Indiana University, Northeastern University, and University of Michigan, http://nwb.slis.indiana.edu]. The definition of a fixed number of communities resulted in different numbers of connections per approach (MSRA: 600 directed connections, RCCA: 300 undirected connections and SRCC: 385 undirected connections).
Additionally, ICA components were transformed into an “ICA co-activation index” (Rosazza et al., 2012) which is the 179 × 179 product matrix of the summarized ICA component z-scores per brain region. Adding a power factor k to the procedure described above emphasizes either the relative weight of intense (k > 1) or less intense (k < 1) co-activation (for detailed information see Supplementary Methods and Rosazza et al., 2012).
Similarity of the three graph-theoretical methods and their correspondence to ICA was measured using two different approaches: (i) the overall linear Pearson's correlation coefficient r for the comparison of the matrices and (ii) the Jaccard index for the comparison of brain structure lists of the graph-theoretical communities and the list of brain structures comprised in the ICA components. The Jaccard index is defined as the size of the intersection divided by the size of the union of finite sample sets. Each similarity measure was calculated for the graph-theoretical methods to each other (the overall r coefficient on single animals, the Jaccard index on group level) and for each graph-theoretical method with ICA (both similarity measures on group level). For the latter comparison, the average ICA co-activation index matrices for power levels 0.5, 1, and 2 were used to calculate the overall r coefficient and the lists of brain structures comprising the 5 most stable group-ICA components (binarized to z > 0.3) were used to calculate the Jaccard index. For statistical analysis of overall r coefficients paired t-test (among graph-theoretical methods) and ANOVA (group level graph-theoretical methods with ICA) were performed with α < 0.05.
Reproducibility of graph-theoretical approaches and ICA co-activation matrix was evaluated using the normalized variance matrices of all analysis methods. The median of all normalized correlation variances demonstrates the reproducibility of each analysis approach. The confidence-interval of each median was calculated by bootstrapping. Significant effect of methods was tested using the Kruskal-Wallis test. Additionally, mean-subtracted variance matrices were calculated, which show the variability of variance within all pairs of brain areas.
Short-Term Modulation and Paired Network Based Statistics (pNBS)
The sensitivity of the approaches to robustly detect changes in RS networks was evaluated by investigating the short-term modulation of RS networks during whisker stimulation. For this purpose the correlation matrices, calculated from the RS pre and RS post measurements for the two groups (experimental group with and control group without whisker stimulation, n = 13 and n = 12, respectively), were averaged separately resulting in 4 averaged correlation matrices per graph-theoretical approach. From these matrices network graphs were created with an average of 20 directed or 10 undirected connections per node, i.e., the 1,790 strongest directed connections for MSRA and consequently 895 undirected connections for RCCA and SRCC. These conditions resulted in graphs that are fully connected, yet sparse enough to show a topological network organization clearly distinct from that of a random network.
To statistically evaluate the short-term modulation of RS network graphs during whisker stimulation, we implemented a modified version of the network-based statistics (NBS) introduced by Zalesky et al. (2010). The NBS relies on the assumption that group differences in single connections are more likely to be false positive than differences in larger connected components. To each connected component of group differences, a p-value controlled for the family-wise error can be ascribed using permutation testing. For details see Zalesky et al. (2010).
In experimental designs with paired groups, permutations cannot be used. To overcome this limitation we included an additional control experiment with identical experimental parameters with the exception of stimulation between the two RS measurements (cf. Figure 2). Our pNBS approach uses the fisher's z-transformed correlation matrices of all four RS scans per animal (2 for experimental and 2 for control group). For each animal, the differences of correlation values per connection between the two RS correlation matrices of each session (control and experimental) were calculated. These difference matrices were used to calculate the pairwise t-test statistics and for the permutation testing.
The following steps were used to determine modulated connected components (1) Paired t-statistics was computed for both control and experimental group. (2) The 99% quantile of the p-values of the paired t-statistics of the control group was determined to identify a set of supra-threshold links corresponding to 1% hypothetical false positive connections (we hypothesize there are no stimulus dependent connectivity modulations in the control experiment). (3) The same threshold was applied to the paired t-statistics p-values of the experimental group, all connected components above this threshold equal or smaller to the largest component of the control group were eliminated. The size k of the remaining set of connections was stored. In contrast to the traditional NBS these connections might comprise more than one connected component. Finally, a p-value controlling for the family-wise-error was ascribed to the remaining component based on its size using permutation testing.
In each of M permutations, the difference matrices were randomly assigned to either control or experimental group, and the statistics of interest described in steps (1–3) was recalculated. The size k' of the set of supra-threshold links derived from each permutation was determined and stored. The p-value of an observed set of connections is estimated by finding the total number of permutations with k' > k normalized by M.
Positioning of Seed Regions
One critical step during seed region based analysis of RS data is the reproducibility of positioning the seed regions across subjects. The MSRA approach includes an automatic positioning step for all seed regions. Figure 4 demonstrates the reproducibility of this procedure. Although all animals were individually analyzed (including atlas registration and automatic determination of seed regions), the positions of the seed regions after affine registration of all animals remained consistent. The Euclidian distance to the mean centroid per structure was on average 1.4 pixel (0.55 mm) in x-y plane and 0.2 slices (0.2 mm) in z direction.
Figure 4. Consistency of automatically placed seed regions of the MSRA approach. (A) Seed regions in the center of mass of individually matched digital atlas regions of one animal. (B) Overlap of exemplary binarized seed region voxels (n = 25). (C–E) Up-sampled high resolution correlation maps of exemplary seed regions. S1BF, barrel field of primary somatosensory cortex; S2, secondary somatosensory cortex; thMD, mediodorsal thalamus.
Comparison of Methods: Correlation Matrices
The correlation matrix of the MSRA approach differs from both SRCC and RCCA, which are comparable to each other (upper triangles in Figure 5A). This effect is associated with the different distribution of z-values within the correlation matrix. SRCC and RCCA matrices show left shifted distributions with their maxima below z = 0.01, which is below the significance level of the FDR (average over all animals: z = 0.176 and z = 0.174, respectively). However, the histogram of the MSRA approach has a maximum around z = 0.25 (average FDR z = 0.177; Figure 5B). The number of connections is determined by thresholding the correlation matrix with a given z-value. For lower z thresholds, the numbers of connections are similar between SRCC and RCCA, but the MSRA approach has considerably more connections. For threshold values higher than 0.4 this effect is reversed resulting in slightly lower number of connections in MSRA compared to SRCC and RCCA (insert in Figure 5C).
Figure 5. Comparison of graph-theoretical network approaches. (A) Correlation matrices including top 1790 binarized directed connections (lower right triangle). (B) Histograms of Fisher‘s z-correlation values. (C) Number of directed connections as a function of the threshold used for binarization (mean ± SEM, n = 25). Pairwise similarity of graph-theoretical approaches: (D) overall correlation r of correlation matrices (mean ± SEM, n = 25, *p < 0.05, paired t-test) (E) and Jaccard index as a function of the number of directed connections of binarized average matrices. bas, basal ganglia; cxA, association cortex; cxS, sensory cortex; hyp, hypothalamus; lim, limbic system including hippocampus and amygdala; SIn, sensory input; thal, thalamus.
The similarity of the correlation matrices of the different approaches was pairwise determined by the overall r-coefficient of all weighted connections (Figure 5D) and by the Jaccard index of binarized matrices (Figure 5E). The overall r reveals highest similarity between RCCA and SRCC. As expected, the overall similarity between MSRA and SRCC is higher than between MSRA and RCCA networks because the MSRA and the SRCC approach rely on the same seed regions (Figure 5D). The Jaccard index of RCCA and SRCC and RCCA and MSRA networks increases with the number of connections. In contrast, the similarity between MSRA and SRCC networks decreases with the number of connections (Figure 5E).
Of the 20 communities, the number containing at least four nodes is equal for MSRA and SRCC networks (12) but differs for RCCA (15). Five of these RCCA-communities are very small (consisting of only four nodes) indicating slightly higher network segregation. In general, the communities contain nodes that represent anatomically and/or functionally associated brain structures (Figure 6, Supplementary Figure 1).
Figure 6. Communities of the average RS network (n = 25) from the MSRA approach. Node positions are determined using a forced based algorithm. For better visualization the underlying network comprises 1790 strongest connections and only communities that contain at least 4 nodes are shown.
The MSRA network consists of the following communities (Figure 6): Brainstem including medulla with solitary tract and parvicellular and gigantocellular nuclei (CoBS); sensory input including pons and tegmentum (CoSIn); structures of the sensory input (ventral tegmental area, red nucleus, substantia nigra) connected to the corpora mamillaria (CoSIn-com); thalamic structures associated with anterior hippocampus and septum (Cothal) or with mesencephalon and pretectal area (Cothal_l); structures of the sensory cortex (S1, S2, visual and auditory cortex) and the parietal association cortex (CocxS); frontal association areas, cingulum, motor cortex and claustrum (CoFrA); ecto- and entorhinal cortex (link to the limbic system) in community with the temporal association area (Corhin); posterior hippocampus and dentate gyrus (Cohip); limbic structures (mainly amygdala), pallidum (basal ganglia) and piriform cortex as link to limbic system (Coam); hypothalamus and zona incerta associated with the ventromedial thalamus and stria terminalis (Cohyp); nucleus accumbens and olfactory nucleus (Cobas).
The most striking deviation of the SRCC-network communities with respect to the MSRA is that the hippocampal areas are segregated into left and right hemisphere and are more strongly connected to thalamic structures of the same hemisphere (Supplementary Figure 1A). The communities of the RCCA network show even more deviations from the MSRA network communities described above (Supplementary Figure 1B). The most prominent difference is the higher degree of segregation. This mainly includes the cortical sensory structures, (corresponding to the MSRA-CocxS-community), but limbic structures (diagonal band and pallidum, septum and stria terminalis) and the medulla are also affected. On the other hand, the frontal association areas and basal ganglia (corresponding to the communities CoFrA and Cobas in the MSRA-network) are integrated into one large community. Here, the motor cortex is associated with the sensory instead of the frontal association cortex.
In summary, the communities tend to be consistent for all three network analysis methods. However, they each produce some different connections resulting in community variations. The MSRA is more similar to the SRCC at the community level as already shown in Figures 5D,E for the correlation matrices.
Similarity of ICA Components and Network Communities
We used Independent Component Analysis (ICA) to compare these three methods to an established, independent method. The following six most stable non-noise ICA components widely reflect the communities of the graph networks (Figure 7):
ICFrA: This component includes the prelimbic association cortex, cingulum and motor cortex. It corresponds to the community CoFrA of the MSRA network.
ICbas: This component corresponds to the Cobas-Community of the graph network. It includes the olfactory nucleus and the nucleus accumbens as well as the piriform cortex.
IChyp: Comparable to the graph-theoretical Cohyp-Community this IC involves predominantly autonomic regions including hypothalamus and corpora mammillaria.
ICcxS: The ICcxS involves sensory cortical structures (S1, S2, auditory and frontal part of visual cortex) and the parietal and retrosplenial cortices. It reflects the CocxS-community.
ICthal: Mainly thalamic structures associated with hippocampal areas and septum are involved in this component.
ICceb: This component is the most caudal one and covers large parts of the cerebellum, which is connected to structures of the brainstem (solitary tract, reticular nucleus). Apart from the cerebellum, which is omitted in graph-theoretical network analysis, this IC is most likely represented by the brainstem community (CoBS). Of note, the ICA components in the rodent brain found in this study align with the ones already published (Becerra et al., 2011; Liang et al., 2011).
Figure 7. Average (n = 25) RS ICA components used for comparison with the graph-theoretical methods. The color-coded ICA z-score maps of these components are overlaid on the anatomical image. Z-scores represent the correlation between each voxel time-course and the mean time-course of the associated component. For each component, three subsequent coronal slices are shown with the corresponding atlas regions in the row above. Acb, nucleus accumbens; Am, amygdala; Au, auditory cortex; Cb, cerebellum; Cg, cingulum; CoM, corpora mammilaria; CPu, caudate putamen; DG, dentate gyrus; Ent, entorhinale cortex; Fr3, frontal cortex area 3; Gi, gigantocellular reticular nucleus; hcD, dorsal hippocampus; hcV, ventral hippocampus; hyp, hypothalamus; Ins, insula; M1, primary motor cortex; M2, secondary motor cortex; Orb, orbitofrontal cortex; OT, olfactory tubercle; PAG, periaqueductal gray; PCRt, parvicellular reticular nucleus; Pir, piriform cortex; PrL, prelimbic cortex; PtA, parietal association cortex; RS, retrosplenial cortex; S1, primary somatosensory cortex; Sep, septum; SN, substantia nigra; Sol, solitary tract; thG, geniculate thalamus; thL, lateral thalamus; thM, medial thalamus; thRe, reunions thalamic nucleus; Vis, visual cortex; ZI, zona incerta.
The similarity of ICA and network-connectivity was determined by the average r of all weighted connections per graph-theoretical approach and the ICA co-activation index matrix at three different power factors (Figure 8A). Two factor ANOVA reveals a significant effect for the factor “method” (p = 0.017) and for the factor “power factor” (p = 0.048). The MSRA approach shows the highest overall-r values compared to the other two methods indicating more similarity to ICA components. The maximum similarity occurs at power factor 1. Focusing on either weak (power factor < 1) or strong (power factor > 1) connections causes a decrease in similarity (insert, Figure 8A) for the MSRA approach but a further increase for power factor >1.0 for the other two methods.
Figure 8. Similarity of RS ICA components (IC) and graph-theoretical network communities (Co). (A) ICA co-activation index matrix (power factor 1) and its overall correlation with the average correlation matrices of the three graph-theoretical approaches as a function of power factor. (B) Jaccard-Index of brain areas belonging to associated ICA components AND/OR network communities (*p < 0.05, paired t-test, n = 5 communities). (C) Visual overlap of thresholded ICA components (>z = 0.03) and centroids of the brain areas composing the associated MSRA network communities.
Additionally, we compared the binarized ICA components with the network communities obtained by the three different graph-theoretical approaches. Because of the missing cerebellum node in the graph-theoretical analysis ICceb and CoBS were not included in the comparison. The overlap was quantified using the Jaccard index (Figure 8B). For all community/component pairs, the MSRA approach reveals the highest overlap. This effect is manifested in the mean Jaccard index over all five compared community/component pairs. The mean Jaccard index of the MSRA approach is significantly higher than that of the other two methods (p < 0.05, paired t-test, cf. Figure 8B). As illustrated in Figure 8C, ICA components and the MSRA network communities share a great deal of overlap.
In order to quantify reproducibility of the three graph-theoretical approaches and ICA, we used the variance matrix of all animal's correlation (graph-theory) and coactivation (ICA) matrices representing the variance for each connection. The median of connection variances differed significantly for all evaluated methods. It was lowest for the MSRA approach, indicating higher reproducibility compared to the SRCC and RCCA methods (Figure 9A).
Figure 9. Reproducibility of ICA co-activation index (CoAI) and the three graph-theoretical approaches. (A) Medians ± 95% confidence intervals of the variances of all connections demonstrate significantly different reproducibility of analysis methods (p < 0.001, Kruskal-Wallis test). (B) Mean-subtracted variance matrices demonstrate the variance within all pairs of brain areas.
Interestingly, brain structures at different organizational levels show unequal variances, which are especially obvious for both SRCC and RCCA (Figure 9B). In contrast, the variance distribution of the MSRA approach and the co-activation indices is much more homogenous. For all three graph-theoretical approaches, the highest variances occur in connections between cortex and subcortex, especially in thalamo-cortical connections. The MSRA approach also showed that seed regions placed in cortical regions had more reproducible connections to subcortical structures than vice versa (Figure 9B).
Short-Term Modulations of RS Networks Due to Whisker Stimulation
To determine if the graph approaches are sensitive enough to detect RS modulations, we used them to analyze the RS measurements separated by sensory whisker stimulation in the experimental group compared to the controls. The statistical significance of the induced modulations was determined using pNBS.
Compared to the control condition (no stimulation between the two RS measurements), whisker stimulation in the experimental group led to significant (p = 0.032, corrected using pNBS) alterations in network connectivity (Figure 10A). Basically, these alterations occurred between brain structures that had been active during whisker stimulation (Supplementary Figure 2). One exception was the ipsilateral somatosensory barrel field, which did not significantly change connectivity (Supplementary Figure 2, Figure 10A). However, the connections of neighboring somatosensory structures in the CocxS-community and the ipsilateral (left) motor cortex were altered, some with increasing (including ipsilateral upper lip field) but most with decreasing connectivity strength. Additionally, subcortical regions were affected, predominantly hypothalamus and ventromedial thalamus (Cohyp) with decreasing and amygdala (Coam) with increasing connectivity strength (esp. basolateral amygdala). Both hypothalamus and amygdala had weaker connections to the ento- and ecto-rhinale cortex (Corhin), which represents a link to the limbic system (Figure 10A). All four structures with the most altered connections (i.e., the biggest nodes in Figure 10A: the right ventromedial thalamus, right lateral and medial hypothalamus and right basolateral amygdala) were located contralateral to the stimulation side.
Figure 10. Alterations in RS connectivity due to whisker stimulation. (A) Components of connectivity alterations identified by pNBS using the MSRA approach (p = 0.032, FWE corrected) overlaid on the RS network before whisker stimulation (gray edges and nodes). The visualization scheme is the same as in Figure 6. (B) Empirically computed null distribution of component size for all three graph-theoretical approaches. Only the observed component of the MSRA was significant (p < 0.05) after permutation correction. amBL: basolateral amygdala, hyDM: dorsomedial hypothalamus, hyM: medial hypothalamus, S1BF, primary somatosensory cortex barrel field; S1Ulp, primary somatosensory cortex upper lip field; thVM, ventromedial thalamus; l, left hemisphere; r, right hemisphere.
Applying pNBS to graphs obtained with the SRCC and the RCCA method reveals smaller components of altered connectivity compared to the MSRA approach (Supplementary Figure 3). These components overlap those of the MSRA and thus qualitatively confirm the whisker induced RS modulation effect. However, in contrast to the MSRA approach these two methods did not survive permutation testing. This was demonstrated by the calculated empirical null distribution of each method shown in Figure 10B.
Additionally, the initial 1% quantile threshold of the paired t-test between the first and second RS of the control experiment was higher for the MSRA approach (p = 0.0084) than for the other graph-theoretical methods (0.0032 for both SRCC and RCCA). This effect was in concordance with the lower variance across animals shown in Figure 9B and confirms the results of the MSRA approach are more reproducible than SRCC and RCCA.
In summary, the MSRA approach, which uses both hypothesis and data driven properties, is superior for graph-theoretical analysis of RS data. Particularly in combination with the proposed pairwise network based statistical evaluation of changes in connectivity strength this approach facilitates the ability to detect and define the slight differences in RS networks induced by a physiological sensory stimulus.
The aim of this study was to introduce a new method for RS analysis that integrates the advantages of traditional ICA and graph-theoretical analysis. The advantages of this new method were illustrated in an investigation of short-term modulation of rodent RS connectivity using physiological whisker stimulation.
ICA and Seed Region Analysis
Since it is known that modulations of ICA components are hard to quantify and to interpret, recently affords have been made to develop methods that use ICA component time courses to create networks for graph-theoretical analysis. Here the resulting graphs either detail specific human ICA resting state networks such as the Default mode (Park et al., 2014; Ribeiro de Paula et al., 2017) or ICA was used to parcellate the brain high-dimensionally into more or less brain region specific components (Beckmann, 2012). As discussed in Ribeiro de Paula et al. (2017), general drawbacks are the usage of group ICA time courses and an (automatic) template matching which both might be inappropriate for patients with brain abnormalities. The MSRA does not rely on generalized ICA time courses but nevertheless highly overlaps with spatial ICA components confirming some fundamental commonalities. Ribeiro de Paula et al. (2017) also emphasized the importance of being able to perform analyses at single subject level. The novel MSRA approach starts at the single subject level with a pseudo directed multiple seed region correlation analysis steps.
Both ICA and classical seed region analysis (SCA) reflect inter-regional features of functional connectivity, whereupon ICA reflects integrated synchronization among networks, SCA measures the inter-regional (long-range) correlation of a specific area with others. Seed based and ICA based RS networks are highly concordant (Van Dijk et al., 2010; Barkhof et al., 2014). For example, the posterior cingulate cortex seed recapitulates the ICA default mode network in human RS analysis (Fox et al., 2005; Long et al., 2008), and seeds within the anterior insula can be used to render the ICA salience network (Seeley et al., 2009).
Differences between ICA and seed based analysis are also reported. (Ma et al., 2007) compared both analysis methods on human fMRI data and on simulated data with artificially added gaussian distributed random noise and a (human) cardiac signal to simulate structured noise. They found that ICA differs from seed based analysis in terms of extent and location of the detected networks. Additionally, ICA was superior in dealing with structured noise such as aliased cardiac cycles and in reproducibility (due to variations in seed positioning across experiments and different laboratories). Structured noise such as aliased cardiac cycles may be relevant in human studies, but because the heart rate of rats (4.5–7.5 Hz) is so much higher than the sample frequency (0.1 Hz) this was not an issue in our experiment. Reproducibility of the seed based MRSA method was enhanced by automatically positioning the seed in the center of mass of the matched atlas brain regions. Other positions may lead to different results.
Motion, Global Signal Regression and Anti-Correlation
As demonstrated on human fMRI data by several researchers (see (Power et al., 2012) or Goto et al., 2016 for review), motion might produce spurious changes in time-courses of RS in humans. Due to anesthesia and fixation of animals within the scanner motion usually is a minor issue for rodent fMRI. However, we cannot exclude such an effect. Therefore, we calculated the motion quality measures FD (frame wise displacement) and DVARS (global RMS of the differential of all time-courses within the brain) for each animal (Power et al., 2012). Correlating these quality measures showed only a negligible relation between motion and BOLD signal change (Supplementary Figure 4). Nevertheless, sophisticated correction for motion artifacts such as scrubbing (Power et al., 2012) or wavelet despike (Patel et al., 2014) might further improve analysis on animal RS.
Although distinguished anti-correlated (i.e. negative correlations) RS networks are reported for humans (Fox et al., 2005), for rats (Liang et al., 2012a; Schwarz et al., 2013), and mice (Sforazzini et al., 2014), we focused on positive correlations only. Negative correlations are difficult to interpret, especially if global regression is used as a preprocessing step. It is known that global regression introduces artificial negative correlations (Murphy et al., 2009), but on the other hand global signal regression leads to more robust RS networks in rodents with higher spatial specificity (Liang et al., 2012a; Liska et al., 2015).
As seen in the histograms of correlation values for the three graph-theoretical methods, both undirected methods are far more influenced by negative correlations. Omitting negative correlations does not shift the histogram, it produces only a cutoff. The histogram of our MSRA approach instead is shifted completely in the range of positive correlation values. Therefore, we conclude that our method is less dependent on false negative correlations introduced by the global regression.
Nevertheless, the MSRA approach is capable of investigating these anti-correlated networks by simply inverting the sign of the initial correlation matrices. The evaluation of these anti-correlated networks is beyond the scope of this study.
Paired Network Based Statistics (pNBS)
Comparing graphs of functional brain connectivity at different experimental states usually involves mass-univariate test statistics after which the family-wise error rate (FWE) has to be controlled. Choosing a threshold for the p-values of the test statistics is a balance between sensitivity (i.e., true positive rate) and specificity (true negative rate). Since the generic procedures, such as the false discovery rate (FDR) (Genovese et al., 2002), are highly conservative, they result in low sensitivities (i.e., high false negative rate) and thus may not offer sufficient power.
Zalesky et al. (2010) presented a powerful method to control the FWE called network based statistics (NBS). An essential part of this method is permutation testing: the group affiliation of each subject is randomly permutated and the statistic of interest is recalculated. This procedure presupposes that the subjects of both groups are independent and exchangeable. This assumption cannot be made in our experimental paradigm, so we adapted the NBS to match the conditions of paired mass-univariate statistics. For this purpose we introduced a control session without stimulation between the two RS measurements.
Repeated measurements of RS connectivity reveal fair to excellent reliability (Zuo and Xing, 2014; Du et al., 2015). Though we cannot postulate the null hypothesis of no difference between the two RS measurements is true for all connections, intra-individual variations should be low. Thus, we use the control session to define a first level threshold p-value that rejects the null hypothesis only for few connections—in this study 1% of all connections of both RS networks. Though usually higher than the conservative FDR (q = 0.05) this first level p-value can be used as a measure of variation that should be exceeded by our experimental intervention (i.e., whisker stimulation). From this point of view, the reliability of the repeated intra subject RS measurements is crucial. Lower variability in the control study leads to higher first level thresholds with more potential to detect significant differences due to experimental intervention.
Comparing the three graph-theoretical methods, the MSRA approach gave the highest first level threshold and the most significantly modulated connections. Thus, we conclude that the observed low inter-individual variability also reflects a high intra-individual reliability over repeated measurements.
In any animal fMRI studies, anesthesia is an important issue because of its potential side effects on the cardiovascular system and the characteristics of spontaneous neural activity (Nallasamy and Tsao, 2011; Grandjean et al., 2014). The effect of anesthesia on functional connectivity has been investigated by several studies which provide evidence that the connectional architecture of brain networks is preserved at low anesthetic doses (e.g., 1% isoflurane) (Vincent et al., 2007; Greicius et al., 2008; Liang et al., 2012b; Gozzi and Schwarz, 2016) and anesthesia depth should be as low as possible to obtain the major topological features of networks mapped in conscious states. In this study, we controlled anesthesia depth by adjusting to the lowest isoflurane dose while maintaining a constant breathing rate during the fMRI session (see section Materials and Methods). As a result, we detected stable ICA components with remarkable similarity to those described by Becerra et al. (2011) in awake rats.
The Influence of Pseudo Directionality
Comparing the three graph-theoretical methods, we found the most similarity between the undirected methods (SRCC and RCCA). This suggests that the pseudo direction of the MSRA approach was an important factor. The effect of directionality was highlighted by the comparison of MSRA and SRCC. Both methods relied on the very same seed regions and consequently on the same time-courses. The main difference was the correlation procedure: cross-correlation of two predefined regional time-courses for SRCC (leading to undirected networks) and cross-correlating multiple predefined regional time-courses with each voxel time-course for MSRA (leading to pseudo directed networks). Consequently, the MSRA results were more similar to SRCC than to the RCCA method. The higher resemblance of ICA components to MSRA than to SRCC networks was predominantly induced by the pseudo directionality. Additionally, pseudo directionality enhanced reproducibility, which is evident in the lower variability of the MSRA approach compared to both directed graph-theoretical methods.
The pseudo directionality of the MSRA approach is characterized by fixed seed regions but variable target regions. The location of the fixed seed regions are chosen (hypothesis driven) based on anatomical properties provided by the digital atlas of brain structures; whereas locations of the voxels in the target region are determined by the strength of the connection (data driven) reflecting given functional relationships.
This asymmetry may reveal broken reciprocity in brain connectivity. Connectional reciprocity between areas (i.e., source and target regions) is common but clearly not ubiquitous; especially cortico-subcortical connections are non-reciprocal (Rockland, 2015). Additionally, the brain as a processing system for highly complex data does rely on the principle of convergence of the information particular in bottom-up processing. This convergence implicitly introduces asymmetry within the connectivity matrix. Moreover, the processed top-down information frequently modulates the bottom-up information by negative feed-back loops e.g. within the cortex by inhibitory connections and thalamo-cortical loops (Robinson, 2017) or even to sub-cerebral structures like the spinal cord in pain conditions (Fields, 2004).
Methodically, RS analysis is constrained by the different brain area volumes and the MR spatial resolution. In this context e.g. the connections of a small thalamic region into much larger cortical areas will not be identical to those found by placing the seed in the corresponding cortical area (see chapter anesthesia for more details in that particular regard).
The MSRA approach reflects such asymmetries by combining hypothesis and data driven as well as anatomical and functional features. This is an important factor that distinguishes our MSRA method from all other (symmetry forcing) techniques used to analyze RS connectivity.
Pseudo directionality and high reproducibility are MSRA properties that also allow for deeper insides into anesthesia effects. Low anesthesia leads to reduced reproducibility and impaired thalamo-cortical connectivity (Liang et al., 2013; Kirsch et al., 2017). While the MSRA approach is characterized by an improved reproducibility across animals (cf. Figure 9), impaired thalamo-cortical connectivity (represented by enhanced variability of connectivity strength) is evident in all three graph-theoretical methods. Interestingly, the pseudo directionality of the MSRA method revealed increased variability of thalamo-cortical connections for seed regions in the thalamic regions, but, not in the cortical regions.
The thalamus is a highly heterogeneous brain region consisting of small nuclei with functional and anatomical diverse connectivity. It is likely that, placement of seed regions within the thalamus reduces the sensitivity of functional connectivity mapping due to this heterogeneity. The data driven correlations of thalamic voxels to seed regions within the much bigger cortical areas were highly reliable. Consequently, the MSRA approach might be a more valid approach for the investigation of thalamo-cortical connectivity, e.g. at different states of consciousness.
Resting State Connectivity Modulation Due to Whisker Stimulation
Using the new MSRA approach, we were able to detect distinct short-term RS connectivity modulations due to unilateral whisker stimulation. The modulated network consisted of two subnetworks. One network was comprised of the somatosensory and motor cortex, and the other larger one, was comprised of the thalamus, hypothalamus and amygdala. Most of these structures were activated during whisker stimulation, predominantly those belonging to the sensorimotor network. Most connections, especially those between the thalamic and hypothalamic nuclei and the cortical areas, were diminished during RS. This effect is also present in resting state connectivity after the stimulation period. Such a deactivation of resting state connectivity during task performance has been observed in the default mode network (DMN) (Shulman et al., 1997; Greicius et al., 2003). Li et al. (2012) investigated this specific property of the DMN and observed an enhanced extrinsic connectivity between constituent regions together with decreased intrinsic self-inhibition within these very regions. Paradoxically, the combination of both phenomena leads to the observed deactivation patterns. The authors suggest that this dynamic results in an increase of the DMN's sensitivity to sensory inputs and may optimize distributed processing during task performance (Li et al., 2012). This principle might be part of RS modulation in general and is still present in RS connectivity shortly after the stimulation period.
Whisking is an important tool rodents use to seek for food or react to a possible thread. To initiate an appropriate behavior, neural processing of information received from whisking involves not only sensorimotor areas of the cortex but also subcortical structures such as hypothalamus (Mogenson et al., 1980). The strongest modulated structures of the subcortical network were the right ventromedial thalamus, the left dorsomedial and the right medial hypothalamus and the right basolateral amygdala.
The contralateral ventromedial thalamus is part of the lemniscal pathway that mediates afferent excitatory projections from whisker to somatosensory cortex (Yu et al., 2006; Diamond et al., 2008). The hypothalamus is generally involved in the regulation of metabolic processes and the autonomic nervous system, the dorsomedial hypothalamus takes part in the regulation of blood pressure and heart rate (Stotz-Potter et al., 1996). Thus, its decreased connectivity might be a response to stress induced by the stimulation. The medial hypothalamus is part of circuitry involved in motivated, i.e., defensive, behaviors, Swanson (2000) and Canteras (2002). Lesions in the lateral hypothalamus profoundly impair the ability to orient to stimuli on the contralateral side (Marshall et al., 1971). This deficit is not from motor impairments, but rather a lack of responsiveness to the stimulus, which indicates a direct connection between sensory input and the lateral hypothalamus (Marshall et al., 1971; Northrop et al., 2010). In contrast to the decreased connectivity between hypothalamus and cortex, the RS connection between hypothalamus and amygdala (mainly basolateral amygdala) was strongly enhanced by whisker stimulation. The lateral amygdala nucleus plays a dominant role in emotional learning and fear conditioning (Pape and Pare, 2010). It receives sensory inputs from the cortex and the thalamus (Ledoux et al., 1990), controls their strength and interferes with the acquisition of fear memory (Ehrlich et al., 2009).
Our results indicate that resting state modulation due to sensory stimulation reflects the impression of a prior sensation and related motor output, but it also involves neuronal circuits known to serve basic processes like fear conditioning and emotional learning initiated by the stimulus.
The biological function of resting state networks is not completely understood, but one hypothesis interprets resting state connectivity as a functional gate that allows the retention of prior information and may influence prospective task-dependent network recruitment and related behavioral output (Deco and Corbetta, 2011). The role of resting state networks relating to human sensorimotor learning and memory consolidation has been described in this context (Albert et al., 2009; Mazoyer et al., 2009; Tambini et al., 2010; Gregory et al., 2014). However, its relevance for emotional learning and conditioning related to sensory stimuli such as touch or pain should be subject for further investigation.
In this study we introduced a powerful new method to analyze resting state functional connectivity. The MSRA approach integrates classical seed based correlation and modern graph-theory, as well as hypothesis and data driven analysis (anatomically chosen seed and functionally correlating target regions). In comparison to two undirected graph-theoretical approaches, it resembles ICA components best and is characterized by its high specificity and reproducibility and less influence from negative correlations. In combination with an adaptation of the network based statistics to paired samples, it promises to be a powerful tool to investigate short-term modulations of sensory stimuli related resting state connectivity and ultimately impact our understanding of basic brain functions like fear to higher functions such as learning and memory and consciousness.
SK established and performed data analysis and wrote the manuscript, BdCA performed the fMRI measurements, MU initiated the project, and AH conducted the project, contributed to the experimental design and rational and edited the manuscript. All authors discussed the results and commented on the manuscript.
Conflict of Interest Statement
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.
The authors thank Sandra Strobelt and Johannes Kaesser for their excellent technical assistance. This work was supported by the BMBF (grant numbers NeuroImpa 01EC1403C, NeuroRad 02NUK034D).
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2018.00334/full#supplementary-material
Albieri, G., Barnes, S. J., De Celis Alonso, B., Cheetham, C. E., Edwards, C. E., Lowe, A. S., et al. (2015). Rapid bidirectional reorganization of cortical microcircuits. Cereb. Cortex 25, 3025–3035. doi: 10.1093/cercor/bhu098
Alonso Bde, C., Lowe, A. S., Dear, J. P., Lee, K. C., Williams, S. C., and Finnerty, G. T. (2008). Sensory inputs from whisking movements modify cortical whisker maps visualized with functional magnetic resonance imaging. Cereb. Cortex 18, 1314–1325. doi: 10.1093/cercor/bhm163
Biswal, B., Yetkin, F. Z., Haughton, V. M., and Hyde, J. S. (1995). Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn. Reson. Med. 34, 537–541. doi: 10.1002/mrm.1910340409
Braun, U., Plichta, M. M., Esslinger, C., Sauer, C., Haddad, L., Grimm, O., et al. (2012). Test-retest reliability of resting-state connectivity network characteristics using fMRI and graph theoretical measures. Neuroimage 59, 1404–1412. doi: 10.1016/j.neuroimage.2011.08.044
Calhoun, V. D., Adali, T., Pearlson, G. D., and Pekar, J. J. (2001). A method for making group inferences from functional MRI data using independent component analysis. Hum. Brain Mapp. 14, 140–151. doi: 10.1002/hbm.1048
Cao, H., Plichta, M. M., Schafer, A., Haddad, L., Grimm, O., Schneider, M., et al. (2014). Test-retest reliability of fMRI-based graph theoretical properties during working memory, emotion processing, and resting state. Neuroimage 84, 888–900. doi: 10.1016/j.neuroimage.2013.09.013
de Celis Alonso, B., Sergeyeva, M., Brune, K., and Hess, A. (2012). Lateralization of responses to vibrissal stimulation: connectivity and information integration in the rat sensory-motor cortex assessed with fMRI. Neuroimage 62, 2101–2109. doi: 10.1016/j.neuroimage.2012.05.045
Du, H. X., Liao, X. H., Lin, Q. X., Li, G. S., Chi, Y. Z., Liu, X., et al. (2015). Test-retest reliability of graph metrics in high-resolution functional connectomics: a resting-state functional MRI study. CNS Neurosci. Ther. 21, 802–816. doi: 10.1111/cns.12431
Fair, D. A., Cohen, A. L., Dosenbach, N. U., Church, J. A., Miezin, F. M., Barch, D. M., et al. (2008). The maturing architecture of the brain's default network. Proc. Natl. Acad. Sci. U.S.A. 105, 4028–4032. doi: 10.1073/pnas.0800376105
Fan, L., Li, H., Zhuo, J., Zhang, Y., Wang, J., Chen, L., et al. (2016). The human brainnetome atlas: a new brain atlas based on connectional architecture. Cereb. Cortex 26, 3508–3526. doi: 10.1093/cercor/bhw157
Fox, M. D., Snyder, A. Z., Vincent, J. L., Corbetta, M., Van Essen, D. C., and Raichle, M. E. (2005). The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proc. Natl. Acad. Sci. U.S.A. 102, 9673–9678. doi: 10.1073/pnas.0504136102
Fransson, P., Aden, U., Blennow, M., and Lagercrantz, H. (2011). The functional architecture of the infant brain as revealed by resting-state FMRI. Cereb. Cortex 21, 145–154. doi: 10.1093/cercor/bhq071
Goto, M., Abe, O., Miyati, T., Yamasue, H., Gomi, T., and Takeda, T. (2016). Head motion and correction methods in resting-state functional MRI. Magn. Reson. Med. Sci. 15, 178–186. doi: 10.2463/mrms.rev.2015-0060
Grandjean, J., Schroeter, A., Batata, I., and Rudin, M. (2014). Optimization of anesthesia protocol for resting-state fMRI in mice based on differential effects of anesthetics on functional connectivity patterns. Neuroimage 102 (Pt 2), 838–847. doi: 10.1016/j.neuroimage.2014.08.043
Gregory, M. D., Agam, Y., Selvadurai, C., Nagy, A., Vangel, M., Tucker, M., et al. (2014). Resting state connectivity immediately following learning correlates with subsequent sleep-dependent enhancement of motor task performance. Neuroimage 102 (Pt 2), 666–673. doi: 10.1016/j.neuroimage.2014.08.044
Greicius, M. D., Kiviniemi, V., Tervonen, O., Vainionpaa, V., Alahuhta, S., Reiss, A. L., et al. (2008). Persistent default-mode network connectivity during light sedation. Hum. Brain Mapp. 29, 839–847. doi: 10.1002/hbm.20537
Greicius, M. D., Krasnow, B., Reiss, A. L., and Menon, V. (2003). Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proc. Natl. Acad. Sci. U.S.A. 100, 253–258. doi: 10.1073/pnas.0135058100
Grinvald, A., Lieke, E., Frostig, R. D., Gilbert, C. D., and Wiesel, T. N. (1986). Functional architecture of cortex revealed by optical imaging of intrinsic signals. Nature 324, 361–364. doi: 10.1038/324361a0
Harmelech, T., Preminger, S., Wertman, E., and Malach, R. (2013). The day-after effect: long term, Hebbian-like restructuring of resting-state fMRI patterns induced by a single epoch of cortical activation. J. Neurosci. 33, 9488–9497. doi: 10.1523/JNEUROSCI.5911-12.2013
Hess, A., Axmann, R., Rech, J., Finzel, S., Heindl, C., Kreitz, S., et al. (2011). Blockade of TNF-alpha rapidly inhibits pain responses in the central nervous system. Proc. Natl. Acad. Sci. U.S.A. 108, 3731–3736. doi: 10.1073/pnas.1011774108
Hess, A., Stiller, D., Kaulisch, T., Heil, P., and Scheich, H. (2000). New insights into the hemodynamic blood oxygenation level-dependent response through combination of functional magnetic resonance imaging and optical recording in gerbil barrel cortex. J. Neurosci. 20, 3328–3338. doi: 10.1523/JNEUROSCI.20-09-03328.2000
Jonckers, E., Van Audekerke, J., De Visscher, G., Van Der Linden, A., and Verhoye, M. (2011). Functional connectivity fMRI of the rodent brain: comparison of functional connectivity networks in rat and mouse. PLoS ONE 6:e18876. doi: 10.1371/journal.pone.0018876
Kalthoff, D., Seehafer, J. U., Po, C., Wiedermann, D., and Hoehn, M. (2011). Functional connectivity in the rat at 11.7T: impact of physiological noise in resting state fMRI. Neuroimage 54, 2828–2839. doi: 10.1016/j.neuroimage.2010.10.053
Kirsch, M., Guldenmund, P., Ali Bahri, M., Demertzi, A., Baquero, K., Heine, L., et al. (2017). Sedation of patients with disorders of consciousness during neuroimaging: effects on resting state functional brain connectivity. Anesth. Analg. 124, 588–598. doi: 10.1213/ANE.0000000000001721
Ledoux, J. E., Cicchetti, P., Xagoraris, A., and Romanski, L. M. (1990). The lateral amygdaloid nucleus: sensory interface of the amygdala in fear conditioning. J. Neurosci. 10, 1062–1069. doi: 10.1523/JNEUROSCI.10-04-01062.1990
Liang, Z., Li, T., King, J., and Zhang, N. (2013). Mapping thalamocortical networks in rat brain using resting-state functional connectivity. Neuroimage 83, 237–244. doi: 10.1016/j.neuroimage.2013.06.029
Long, X. Y., Zuo, X. N., Kiviniemi, V., Yang, Y., Zou, Q. H., Zhu, C. Z., et al. (2008). Default mode network as revealed with multiple methods for resting-state functional MRI analysis. J. Neurosci. Methods 171, 349–355. doi: 10.1016/j.jneumeth.2008.03.021
Ma, L., Wang, B., Chen, X., and Xiong, J. (2007). Detecting functional connectivity in the resting brain: a comparison between ICA and CCA. Magn. Reson. Imaging 25, 47–56. doi: 10.1016/j.mri.2006.09.032
Margulies, D. S., Bottger, J., Long, X., Lv, Y., Kelly, C., Schafer, A., et al. (2010). Resting developments: a review of fMRI post-processing methodologies for spontaneous brain activity. MAGMA 23, 289–307. doi: 10.1007/s10334-010-0228-5
Mary, A., Wens, V., Op De Beeck, M., Leproult, R., De Tiege, X., and Peigneux, P. (2016). Age-related differences in practice-dependent resting-state functional connectivity related to motor sequence learning. Hum. Brain. Mapp. 38, 923–937. doi: 10.1002/hbm.23428
Mazoyer, B., Houde, O., Joliot, M., Mellet, E., and Tzourio-Mazoyer, N. (2009). Regional cerebral blood flow increases during wakeful rest following cognitive training. Brain Res. Bull. 80, 133–138. doi: 10.1016/j.brainresbull.2009.06.021
Mogenson, G. J., Jones, D. L., and Yim, C. Y. (1980). From motivation to action: functional interface between the limbic system and the motor system. Prog. Neurobiol. 14, 69–97. doi: 10.1016/0301-0082(80)90018-0
Muraskin, J., Dodhia, S., Lieberman, G., Garcia, J. O., Verstynen, T., Vettel, J. M., et al. (2016). Brain dynamics of post-task resting state are influenced by expertise: insights from baseball players. Hum. Brain Mapp. 37, 4454–4471. doi: 10.1002/hbm.23321
Murphy, K., Birn, R. M., Handwerker, D. A., Jones, T. B., and Bandettini, P. A. (2009). The impact of global signal regression on resting state correlations: are anti-correlated networks introduced? Neuroimage 44, 893–905. doi: 10.1016/j.neuroimage.2008.09.036
Northrop, L. E., Polston, E. K., and Erskine, M. S. (2010). Noradrenergic nuclei that receive sensory input during mating and project to the ventromedial hypothalamus play a role in mating-induced pseudopregnancy in the female rat. J. Neuroendocrinol. 22, 1061–1071. doi: 10.1111/j.1365-2826.2010.02049.x
Pape, H. C., and Pare, D. (2010). Plastic synaptic networks of the amygdala for the acquisition, expression, and extinction of conditioned fear. Physiol. Rev. 90, 419–463. doi: 10.1152/physrev.00037.2009
Park, B., Kim, D. S., and Park, H. J. (2014). Graph independent component analysis reveals repertoires of intrinsic network components in the human brain. PLoS ONE 9:e82873. doi: 10.1371/journal.pone.0082873
Patel, A. X., Kundu, P., Rubinov, M., Jones, P. S., Vertes, P. E., Ersche, K. D., et al. (2014). A wavelet method for modeling and despiking motion artifacts from resting-state fMRI time series. Neuroimage 95, 287–304. doi: 10.1016/j.neuroimage.2014.03.012
Pawela, C. P., Biswal, B. B., Cho, Y. R., Kao, D. S., Li, R., Jones, S. R., et al. (2008). Resting-state functional connectivity of the rat brain. Magn. Reson. Med. 59, 1021–1029. doi: 10.1002/mrm.21524
Power, J. D., Barnes, K. A., Snyder, A. Z., Schlaggar, B. L., and Petersen, S. E. (2012). Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion. Neuroimage 59, 2142–2154. doi: 10.1016/j.neuroimage.2011.10.018
Ramos-Cabrer, P., Weber, R., Wiedermann, D., and Hoehn, M. (2005). Continuous noninvasive monitoring of transcutaneous blood gases for a stable and persistent BOLD contrast in fMRI studies in the rat. NMR Biomed. 18, 440–446. doi: 10.1002/nbm.978
Ribeiro de Paula, D., Ziegler, E., Abeyasinghe, P. M., Das, T. K., Cavaliere, C., Aiello, M., et al. (2017). A method for independent component graph analysis of resting-state fMRI. Brain. Behav. 7:e00626. doi: 10.1002/brb3.626
Rosazza, C., Minati, L., Ghielmetti, F., Mandelli, M. L., and Bruzzone, M. G. (2012). Functional connectivity during resting-state functional MR imaging: study of the correspondence between independent component analysis and region-of-interest-based methods. AJNR Am. J. Neuroradiol. 33, 180–187. doi: 10.3174/ajnr.A2733
Sachdev, R. N., Champney, G. C., Lee, H., Price, R. R., Pickens, D. R. III., Morgan, V. L., et al. (2003). Experimental model for functional magnetic resonance imaging of somatic sensory cortex in the unanesthetized rat. Neuroimage 19, 742–750. doi: 10.1016/S1053-8119(03)00093-4
Schwarz, A. J., Gass, N., Sartorius, A., Risterucci, C., Spedding, M., Schenker, E., et al. (2013). Anti-correlated cortical networks of intrinsic connectivity in the rat brain. Brain Connect 3, 503–511. doi: 10.1089/brain.2013.0168
Seeley, W. W., Crawford, R. K., Zhou, J., Miller, B. L., and Greicius, M. D. (2009). Neurodegenerative diseases target large-scale human brain networks. Neuron 62, 42–52. doi: 10.1016/j.neuron.2009.03.024
Sforazzini, F., Schwarz, A. J., Galbusera, A., Bifone, A., and Gozzi, A. (2014). Distributed BOLD and CBV-weighted resting-state networks in the mouse brain. Neuroimage 87, 403–415. doi: 10.1016/j.neuroimage.2013.09.050
Shulman, G. L., Fiez, J. A., Corbetta, M., Buckner, R. L., Miezin, F. M., Raichle, M. E., et al. (1997). Common blood flow changes across visual tasks: II. decreases in cerebral cortex. J. Cogn. Neurosci. 9, 648–663. doi: 10.1162/jocn.1922.214.171.1248
Sierakowiak, A., Monnot, C., Aski, S. N., Uppman, M., Li, T. Q., Damberg, P., et al. (2015). Default mode network, motor network, dorsal and ventral basal ganglia networks in the rat brain: comparison to human networks using resting state-fMRI. PLoS ONE 10:e0120345. doi: 10.1371/journal.pone.0120345
Smith, S. M., Vidaurre, D., Beckmann, C. F., Glasser, M. F., Jenkinson, M., Miller, K. L., et al. (2013). Functional connectomics from resting-state fMRI. Trends Cogn. Sci. 0.17, 666–682. doi: 10.1016/j.tics.2013.09.016
Stotz-Potter, E. H., Willis, L. R., and Dimicco, J. A. (1996). Muscimol acts in dorsomedial but not paraventricular hypothalamic nucleus to suppress cardiovascular effects of stress. J. Neurosci. 16, 1173–1179. doi: 10.1523/JNEUROSCI.16-03-01173.1996
Tung, K. C., Uh, J., Mao, D., Xu, F., Xiao, G., and Lu, H. (2013). Alterations in resting functional connectivity due to recent motor task. Neuroimage 78, 316–324. doi: 10.1016/j.neuroimage.2013.04.006
Van Den Heuvel, M. P., Stam, C. J., Boersma, M., and Hulshoff Pol, H. E. (2008b). Small-world and scale-free organization of voxel-based resting-state functional connectivity in the human brain. Neuroimage 43, 528–539. doi: 10.1016/j.neuroimage.2008.08.010
Van Dijk, K. R., Hedden, T., Venkataraman, A., Evans, K. C., Lazar, S. W., and Buckner, R. L. (2010). Intrinsic functional connectivity as a tool for human connectomics: theory, properties, and optimization. J. Neurophysiol. 103, 297–321. doi: 10.1152/jn.00783.2009
Vincent, J. L., Patel, G. H., Fox, M. D., Snyder, A. Z., Baker, J. T., Van Essen, D. C., et al. (2007). Intrinsic functional architecture in the anaesthetized monkey brain. Nature 447, 83–86. doi: 10.1038/nature05758
Zhong, Y., Wang, H., Lu, G., Zhang, Z., Jiao, Q., and Liu, Y. (2009). Detecting functional connectivity in fMRI using PCA and regression analysis. Brain Topogr. 22, 134–144. doi: 10.1007/s10548-009-0095-4
Keywords: resting state connectivity, whisker stimulation, multi seed correlation, ICA, fMRI, short-term modulation, rats
Citation: Kreitz S, de Celis Alonso B, Uder M and Hess A (2018) A New Analysis of Resting State Connectivity and Graph Theory Reveals Distinctive Short-Term Modulations due to Whisker Stimulation in Rats. Front. Neurosci. 12:334. doi: 10.3389/fnins.2018.00334
Received: 01 December 2017; Accepted: 30 April 2018;
Published: 23 May 2018.
Edited by:Oscar Arias-Carrión, Hospital General Dr. Manuel Gea Gonzalez, Mexico
Reviewed by:Andrea Soddu, University of Western Ontario, Canada
Hyunjin Park, Sungkyunkwan University, South Korea
Copyright © 2018 Kreitz, de Celis Alonso, Uder and Hess. 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 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: Andreas Hess, firstname.lastname@example.org