Contextual experience modifies functional connectome indices of topological strength and efficiency

Stimuli presented at short temporal delays before functional magnetic resonance imaging (fMRI) can have a robust impact on the organization of synchronous activity in resting state networks. This presents an opportunity to investigate how sensory, affective and cognitive stimuli alter functional connectivity in rodent models. In the present study we assessed the effect on functional connectivity of a familiar contextual stimulus presented 10 min prior to sedation for imaging. A subset of animals were co-presented with an unfamiliar social stimulus in the same environment to further investigate the effect of familiarity on network topology. Rats were imaged at 11.1 T and graph theory analysis was applied to matrices generated from seed-based functional connectivity data sets with 144 brain regions (nodes) and 10,152 pairwise correlations (after excluding 144 diagonal edges). Our results show substantial changes in network topology in response to the familiar (context). Presentation of the familiar context, both in the absence and presence of the social stimulus, strongly reduced network strength, global efficiency, and altered the location of the highest eigenvector centrality nodes from cortex to the hypothalamus. We did not observe changes in modular organization, nodal cartographic assignments, assortative mixing, rich club organization, and network resilience. We propose that experiential factors, perhaps involving associative or episodic memory, can exert a dramatic effect on functional network strength and efficiency when presented at a short temporal delay before imaging.


Results
Context exposure significantly reduced network strength. Rats were assigned one of 3 experimental groups that included a naïve control group (no pre-scan exposure; n = 11, with 5 males), a cage only group (context exposure group; n = 10, with 5 males), and a cage + social stimulus group (context/social exposure; n = 9, with 5 males) (Supplementary Figure 1; Details in "Materials and methods" section).
We observed a significant reduction in node strength at all graph densities in the context and context/social groups relative to controls (group × density interaction F 16,216 = 8.5 p < 0.00001; Tukey p < 0.05 for graph densities 5-45% for both prescan manipulation groups compared to controls) (Fig. 1A,B). No differences between the context only and context/social groups was observed. Differences in node strength between context, context/ social and controls was also observed at the nodal level in left and right cingulate (left: F 2,27 = 11.7, p = 0.009; right: F 2, 27 Fig. 1C indicated that context and context/social stimuli both reduced node strength bilaterally in these nodes (p < 0.05).
Context exposure significantly reduced global network efficiency and clustering. No differences in the small world index was observed between the groups ( Fig. 2A). However, we observed significant differences in the global weighted values for the clustering coefficient (group × density interaction F 16,216 = 4.8 p < 0.00001; Tukey p < 0.05 for 10-45% graph densities for both prescan manipulation groups compared to controls) and average characteristic path length (group × density interaction F 16,216 = 4.4 p < 0.00001; Tukey p < 0.05 for 5-45% graph densities for both prescan manipulation groups compared to controls). As expected from the low global clustering and high average path lengths, global efficiency at 10% graph density was significantly lower in context and context/social groups relative to controls (F 2,27 = 7.6, p = 0.002) (Fig. 2D). Compared to global efficiency, and based on the size of nodes in 3D connectome maps shown in Fig. 2E-G, average values for local efficiency were similarly yet less robustly reduced at the individual node level (Fig. 2E-G). Analysis of local Scientific Reports | (2020) 10:19843 | https://doi.org/10.1038/s41598-020-76935-0 www.nature.com/scientificreports/ efficiency values across all nodes in all rats (as a function of local clustering coefficient values) suggested relatively minor reductions local efficiency in context and context/social groups compared to controls (Fig. 2H). This suggested that while context exposure produced robust reductions in global efficiency associated with reduced node strength, clustering, and increased path lengths, this globally reduced efficiency was not widespread in all nodes. Based on this possibility, we hypothesized that the particular arrangement of edges connecting nodes may be important in preserving local efficiency. If so, there would be relative few differences in edge quantifiers that support reduced local efficiency. To investigate this, we calculated search information, which estimates the  www.nature.com/scientificreports/ amount of information (bits) needed to randomly trace shortest paths between node pairs 25 . We plotted search information as a function of path transitivity, which estimates the density of detours lying between node pairs 25 . This analysis indicated no differences in terms of accessibility or degree of difficultly in randomly finding shortest paths between node pairs in control, context, and context/social exposure rats (Fig. 2I). We next plotted the number of edges (termed 'hops') lying along shortest paths between a pair of nodes as a function of routing efficiency (the average of inverse shortest path length) [26][27][28] . As anticipated, differences in efficiency were subtle and limited to single 'hops' which had slightly lower routing efficiency in context and context/social groups relative to controls (Fig. 2J). Finally, we note that in randomized functional connectivity matrices the estimated values for clustering coefficient and characteristic path lengths were lower than in original matrices, although these preserved node (D) Global efficiency at a threshold of 10%. (E-G) 3D connectome maps with nodes scaled by local efficiency in control (E), context exposed (F), and context/social exposed (G). Sale bar represent edge weights values (Pearson r threshold 0.15). (H) Distribution of local efficiency vs. clustering coefficient values for all nodes and all subjects in control, context and context/social exposed rats. (I) Distribution of search information vs. path transitivity values for all edges and all subjects in control, context and context/social exposed rats. (J) Distribution of number of steps or hops between node pairs vs. routing efficiency values for all edges and all subjects in control, context and context/social exposed rats. Matrix density threshold for (D-J) is 10%. In (B,C), blue, orange, red line plots are for randomized (r) versions of functional connectivity matrices. In (H-J), gray color dots represent data points for nodes or edges of randomized graphs. Data are mean ± standard error. *Context and **context/social group significantly different from control (ANOVA with Tukey-Kramer post hoc test). www.nature.com/scientificreports/ strength differences between the groups (Fig. 2). Also, while global efficiency values were greater in randomized versus original matrices (data not shown), local efficiency and path transitivity values were much lower in rewired networks (Fig. 2H,I).
No significant effects of context or context/social exposure on additional network metrics. No effect of stimuli on global modularity was observed (Supplementary Figure 2A). The number of total affiliations or groups did not vary between the 3 groups (Supplementary Figure 2B). To further investigate the role of nodes within the identified modules, we used an approach published by Guimerà and Amaral 29 and updated for functional MRI data sets by Meunier et al. 30 (Supplementary Figure 2C-F). There was no effect of cage or social exposures on any of these nodal role assignments. Finally, modularity index values were significantly lower in randomly rewired functional networks compared to original versions. Also, there were lower fMRI hubs in randomly rewired versus original networks and the distribution of participation coefficient values indicated that rewired networks only had nodes with roles as non-hub kinless and non-hub connectors.
No effects of the stimuli was observed on assortative mixing or the rich club index ( Supplementary Figure 3A,B). Assortativity index values were lower in randomized networks versus original versions. Similarly, rich club index values at a specific range of k-level values between 10 and 16 were lower in the randomized versus original functional networks.
We carried out the targeted node removal process 31-34 based on descending rank of either node strength, node degree, betweenness centrality, or eigenvector centrality scores (Supplementary Figure 3C Context and context/social exposure resulted in a switch in high eigenvector centrality node locations from cortex to subcortical regions. The group-averaged 3D connectome maps in Supplementary Figure 3E suggested that both context and context/social stimuli prior to scanning reduced node eigenvector centrality values. We sorted average node eigenvector centrality values from highest to lowest and selected the top 10 in each group. Results are summarized in Fig. 3. The top 10 nodes with highest eigenvector centrality values (which we refer to here as 'hubs') in control rats were cortical. These included cortical areas that had reduced node strength in context and context/social groups. Many of these same regions showed a significant reduction in eigenvector centrality and this effect was observed in the trunk (F 2,27 = 4.2, p = 0.03) and shoulder regions of primary somatosensory cortex (left: F 2,27 = 11.9, p = 0.0009), and parietal (left: F 2,27 = 3.2, p = 0.05), cingulate (left: F 2,27 = 6.3, p = 0.006; right: F 2,27 = 3.7, p = 0.04), retrosplenial (left: F 2,27 = 3.6, p = 0.04) Figure 3. Context and context/social exposure resulted in a switch in high eigenvector centrality node locations from cortex to subcortical regions. Eigenvector centrality rankings of top 10 nodes in control rats and in context exposed rats (arrows above bar graphs indicate top 10 nodes in control vs. context and context/ social exposed rats). Network density set at 10%. Data are mean ± standard error. *Context and **context/social group significantly different from control (ANOVA with Tukey-Kramer post hoc test). Abbreviations: L, left; R, right; S1, primary somatosensory; Tr, trunk; Sh, shoulder; Cg, cingulate, Pt, parietal; RSR, rostral retroplenium; V2, secondary visual cortex; M1, primary motor cortex; VMH, ventromedial hypothalamus; LH, lateral hypothalamus; PM, posteromedial thalamus; AHA, anterior hypothalamuc area; Me, medial amygdala; CA3v, ventral CA3 of hippocampus. www.nature.com/scientificreports/ and primary motor cortical regions (left: F 2,27 = 3.8, p = 0.04). The group-averaged 3D connectome maps for the familiar context and context/social exposure groups in Supplementary Figure 3E also suggest that cortical hubs are significantly suppressed. Indeed, these maps suggest that nodes in deep subcortical areas have the highest eigenvector values. Upon sorting nodes from highest to lowest eigenvector values in these groups, we observed that the top 10 nodes were subcortical-largely hypothalamic. Among these regions, the anterior hypothalamic area had eigenvector centrality values that were greater in context and context/social exposed compared to control animals (left: F 2,27 = 5.1, p = 0.01).
An arousing saline injection stimulus prior to imaging produced network topological changes that are diametrically distinct from familiar context exposure. We reanalyzed an MRI data set in which rats were scanned either 1 or 24 h after an intraperitoneal injection of sterile saline solution. The stimulus in this case was unexpected by rats, arousing due to handling during injection, and in one group of rats the injection was administered at the same temporal delay before setup for imaging and image acquisition as the context and context/social exposures of the present study. These groups were compared to a naive control group of age, sex, weight, vendor and imaging protocol matched rats that were scanned directly upon arrival at the imaging facility and did not receive any pre-scan intervention. We included this reanalysis as a way to assess if an arousing stimulus produced topological changes comparable to those produced by context and context/social exposure in the present study. As a first step in validating this comparison we observed that values for node strength, small world index, clustering coefficient, characteristic path length, modularity, assortativity and rich club index in control rats imaged at 4.7 T were almost identical to those calculated for functional connectivity data in control rats imaged at 11.1 T (see control rat data in Fig. 4 versus Figs. 1, 2). 3D connectome maps in Fig. 4A show stronger overall connectivity and average node strength in rats injected 1 h prior to image compared to a naive control group (that did not receive any injection or manipulation prior to imaging). Node strength was much greater in 1 h vs. 24 h and control rats across several network density thresholds (group × density interaction F 16,168 = 2.6 p < 0.0015; Tukey p < 0.05 for 5-45% graph densities for 1 h pre-scan injection vs. controls) (Fig. 4B).
The small world index did not differ between 1 h injected rats versus 24 h and control animals (repeated measures ANOVA F 16,168 = 1.4, p = 0.1; Fig. 4C). Injecting rats 1 h prior to scanning resulted in a significantly greater clustering coefficient (group × density interaction F 16,168 = 21.8 p < 0.00001; Tukey p < 0.05 for 5-45% graph densities for 1 h pre-scan injection vs. controls ), modularity (group × density interaction F 16,168 = 7.2 p < 0.00001; Tukey p < 0.05 for 5-45% graph densities for 1 h pre-scan injection vs. controls) and assortative mixing (group × density interaction F 16,168 = 8.4 p < 0.00001; Tukey p < 0.05 for 5-45% graph densities for 1 h pre-scan injection vs. controls) than control rats. Interestingly, the injection 1 h prior to scanning did not alter characteristic path length. Rich club organization differed between the groups (for k-level 6-8 F 2,21 > 6.9, p < 0.005 and for k level = 14-24 F 2,21 > 6.7, p < 0.01; one way ANOVA) ( Fig. 4D-H). In sum, the response to a general arousing stimulus such as an intraperitoneal injection administered 10 min prior to sedation and 50-60 min prior to fMRI scanning caused an increase in functional connectivity network measures. This is in stark contrast with the network changes caused by exposure to a context and context/social groups. To further determine if the general arousing injection procedure produced network changes that are distinct from context or context/social exposure, we analyzed node strength in regions that had reduced node strength under these experimental conditions (Supplementary Figure 4). This additional analysis revealed differences in local connectivity strength between the two studies. The 1 h injection increased node strength in cingulate and insular cortices and dorsomedial striatum, and reduced node strength in the rostral retrosplenial cortex, but had no effect on any other region that were responsive to presentation of a familiar context (Supplementary Figure 4; Statistical results in Supplementary Table 1). In addition, the injection procedure significantly increased node strength in a different set of structures. These included the anterior, basal, central, lateral, medial amygdala nodes and bed nucleus of stria terminalis, prelimbic cortex, primary and secondary motor cortices, caudal piriform cortex, lateral septal nucleus, medial preoptic area, central grey, nucleus accumbens, ventromedial striatum, ventral pallidum (Statistical results in Supplementary Table 2).

Discussion
The data indicate that exposure to a familiar context either in the absence or presence of an additional social stimulus resulted in a significant reduction in connectivity strength in rat brain. Analyses of 72 bilateral brain regions further indicated that node strength values were reduced in areas of the cortex known for their roles in long-term memory. The regions included left and right cingulate, parietal, insular, perirhinal, retrosplenial, and primary somatosensory cortices, and dorsal striatum. In addition to changes in node strength, we observed a significant reduction in weighted measures of clustering and efficiency. This suggests that because of the overall reduction in node strength there is a significant loss of strongly weighted triad structures (weighted clusters) in rats exposed to a familiar context, and this in turn is associated with a reduced potential for information transfer. Finally, while in control rats the top 10 nodes with highest eigenvector centrality values were cortical, in rats exposed to a familiar context the top 10 hub nodes were subcortical. The top 10 hub regions in controls were the same structures that showed reduced connectivity strength in context-exposed rats. Thus, exposure to a familiar context modified neural interactions in such a way that suppressed synchronous cortical hubs, which resurfaced in hypothalamic regions. In our current study design we investigated the effect of a transient social interaction prior to imaging and hypothesized that network topology would distinguish the familiar context from the unfamiliar, novel social stimulus. However, rats presented with a novel social stimulus rat within the familiar test cage had very similar values to the context alone in all global and local network quantifiers. We interpret this lack of effect of the social stimulus in two possible ways. In the first scenario, the regions involved might prioritize the processing of previously encountered stimuli versus novel stimuli. Synaptic changes in response Scientific Reports | (2020) 10:19843 | https://doi.org/10.1038/s41598-020-76935-0 www.nature.com/scientificreports/ to the familiar stimulus generate more robust or faster neuronal responses than stimuli that have yet to induce plasticity related changes in distributed central synapses 24,36,37 . Second, both social and contextual stimuli alter connectivity patterns in slightly overlapping regions and this results in an inability to distinguish subtle differences between the two experimental conditions. In the latter scenario fMRI is unable to detect neuronal activity patterns in microcircuits considerably below voxel resolution. Among unanswered questions in this study, two stand out as critical to address in future studies. The first relates to neuronal activity mechanisms underlying the observed reductions in cortical node strength, efficiency, clustering, and cortical hubness. A second question relates to the lack of effect of the unfamiliar social stimulus relative to familiar context exposure 38 . The first question could be addressed in future studies using local field potential recordings in multiple cortical regions of rats tested under similar experimental conditions as our present neuroimaging study. It is well established that the frequency bandwidth of spontaneous BOLD signal oscillations falls within previously characterized infraslow frequencies (< 0.5) that cover a large expanse of rat cortex 39,40 . Studies in primate visual cortex have also shown that normalized spectral frequencies in the band limited power for gamma (40-100 Hz) share a high degree of statistical dependency (mutual information) with the fMRI BOLD signal while alpha (8)(9)(10)(11)(12), beta (18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30), and higher frequency (0.9-3 kHz) multi-unit activity much less so 41 . Changes in BOLD signal amplitudes scale in a step wise manner with changes in gamma www.nature.com/scientificreports/ power 41 . Thus, somatodendritic processing at band limited power in gamma and infraslow ranges might correspond with synchronous BOLD oscillations underlying functional connectivity. Alterations in gamma band power for a particular frequency might indicate large scale changes in neural processing, which produces asynchronous cortical BOLD activity. In other words, a switch from gamma to another band limited power spectral frequency, i.e., theta, involved in the processing of information related to a familiar context might underlie reduced connectivity strength in cortical regions. However, it is important to note that this viewpoint is incomplete and research on the large synaptic activity underlying BOLD signal fluctuations is ongoing 42 . Although band limited power for gamma shows strong statistical dependency with BOLD, using a sliding window approach others have shown that the relationship between spontaneous fluctuating BOLD signals and underlying neural signals is more complex 43 . Using this approach, strong correlations with spontaneous BOLD activity are see in gamma, theta and beta bands 43 . Work in human subjects have also reported a closer relationship between alpha and regions of the salience and dorsal attention networks, and theta and frontal parietal and default mode networks 44 . Thus, the relationship between BOLD and underlying large scale synaptic activity may be behavior-, cognition-, and affective state-specific.
Regarding the second open question stated above, there have been prior reports suggesting a strong impact of familiarity to an environmental context in habituation (learning) processes involved in social interaction and recognition in rats 38 . Behavioral changes observed in response to novel social interactions may depend on familiarity with the associated context (59). The observed reductions in strength, clustering, efficiency and cortical hubs following familiar context exposure could thus correspond to reconsolidation involving changes in underlying synaptic activity which can affect inter-regional communication in distributed cortical neuronal populations 44 . What is perplexing is the shift to a less efficient topology. One possible explanation for this change could be that in order to maintain the fidelity of synaptic changes related to the already experienced contextual environment, network topology transiently moves to a less efficiency state as a way of minimizing novel-stimulus interference. Reductions in network efficiency during learning has been previously reported 45 . This would prevent the novel social stimulus or any unfamiliar or unrelated stimuli from modifying overall network topology during a state of reconsolidation 44 . However, its important to emphasize that the applied graph algorithms may have not been sensitive to modality-specific changes in functional networks. Future refinements to stimulus presentations at short time delays prior to scanning may better isolate memory versus non-memory effects on functional networks. The brain regions shown to have greater node strength and eigenvector centrality in control relative to stimulus exposed rats have been studied in part for their roles in learning and memory. Specifically, these regions may play significant roles in spatial learning, and affective and episodic memories. However, their roles in social learning and recognition are not established. Parietal cortex is known to play important roles in long-term spatial memory functions in rodents and primates 46 . Anterior cingulate is involved in establishment of long-term valence associations, such as pain and reward 47,48 . Hippocampally-driven reconsolidation processes may involve communication with retrosplenial cortical circuits for long-term storage 49,50 . Insular cortex plays an important role in memory of aversive experiences 51,52 . It also plays significant roles in reward based learning as does the striatum [53][54][55] . The somatosensory cortex shares input specific synaptic inputs with the hippocampus related to sensory-memory consolidation 56 .
Reduced network efficiency and clustering are inversely related to an increase in the average number of intermediate steps separating any pair of nodes in a network (increased characteristic path length). We note that this should not be taken here to imply anatomical steps or axonal connections. Instead, we infer from these network quantifiers variations in the degree of correlated BOLD activity across regions. Regions with highest statistical correlations in BOLD thus represent shortest routes of communication compared with regions with intermediate, low, or weakest pairwise correlations 21,44 . A short (and efficient) path would be characterized by regions linked by the highest or strongest statistical correlations in BOLD. Although an underlying neurophysiological mechanism remains unclear, synchronous activity plays a central role in the temporal entrainment of spike delays to theta rhythms in distributed neuronal populations communicating cortical regions and hippocampus [57][58][59] . Such spike-timing adjustments to behaviorally-driven oscillations could play a significant role in information transfer and long-term memory and plasticity 44 . As indicated above, a change in rhythmic activity may impact band limited power of specific frequencies in cortex and this has the potential of desynchronizing BOLD activity when shifting from a baseline state of rest to a state where a behavioral task was carried out shortly before imaging. At this time, the mechanisms remain unclear. However, there is support in the literature for a reduced global efficiency associated with high rates of learning 45 .
We further assessed quantifiers that while not having solid grounding in known neurophysiological mechanisms, may help deepen our understanding of the structure of networks in relation to their efficiency. Thus, we did not observe differences in the degree of difficulty finding edge sequences that form random shortest paths between nodes (search information) nor the number of edges that constitute detours in randomly transited paths between nodes (path transitivity). We did observe that edges along single short 'hops' or individual steps between two nodes constituted slightly higher efficient routes for information transfer in control relative to context exposed rats. This is likely associated with a higher number of close local connections (higher clustering coefficient) with short distances (inversely related to high correlation coefficients edge-connected nodes) in controls relative to context-exposed groups. In sum, calculations of shortest path lengths and efficiency may reveal interesting and novel views of the mechanisms of learning and memory, both in terms of organization of edge weights and nodal connectivity patterns involved.
Contrary to previous results in human subjects undergoing specific cognitive and learning tasks 60,61 , we did not observe changes in global modularity in response to familiar or novel stimuli. Our analysis of node roles also indicated no overall differences in nodal cartographic assignments according to Guimerà and Amaral 29 or Meunier and colleagues 30 . Consistent with a lack of differences in modular structures, we also observed no differences in assortativity and rich club indices (suggesting no differences in degree distributions). As a way to Scientific Reports | (2020) 10:19843 | https://doi.org/10.1038/s41598-020-76935-0 www.nature.com/scientificreports/ assess whether changes in node strength and hub reallocations impact network robustness, we used a targeted node elimination strategy. No differences in network robustness were observed either in response to graded removal of high-to-low strength nodes, or betweenness or eigenvector centrality nodes (we also removed nodes according to their degree values and observed no differences). A random failure model also indicated no differences in overall network robustness (data not shown). A close examination of composite 3D connectome maps in which the nodes were weighted by eigenvector centrality values (Supplementary Figure 3E) showed stark differences in nodal sizes between control and context-exposed groups. In the latter, node sizes for deep subcortical structures in the hypothalamus appear enlarged. The cortical areas with highest centrality in control rats were the same regions with reduced node strength in context-exposed groups. This suggested that these regions were specifically involved in the reorganization of BOLD correlations induced by the familiar contextual stimulus. The role of hypothalamic structures in the context groups is somewhat surprising as these regions are usually not included as seed regions in functional connectivity studies. Deeper structures are subject to large blood vessel pulsations that are typically well above the resting state frequency power. Undersampled Nyquist frequencies can alias into the measured fMRI signal and drive hypothalamic 'synchronous' activity in the presence of suppressed cortical connectivity strength. However, despite this possible non-neuronal mechanism there are studies that have observed a role for hypothalamic nuclei in functional connectivity in models of pain, in human studies of migraine headache, and in studies of obesity where increases in hypothalamic functional connectivity have been reported [62][63][64] .
We reanalyzed a database of male rats imaged under experimental conditions similar to the present study. A subset of these rats were given an intraperitoneal sterile injection 10 min before anesthetic induction and preparations for imaging. As a comparison, we also analyzed a group of rats undergoing the same injection procedure 24 h before imaging and a naive group never administered any injection. Results from this analysis indicated that rats that received a saline injection at a short delay before imaging had high connectivity strength, along with other quantifiers that strongly supported a general increase in functional connectivity involving to a large extent cortical and thalamic areas. A possible mechanism for the observed increase in connectivity strength is the release of corticosterone which can act upon memory networks via glucocorticoid receptors in hippocampal, insular, and basolateral amygdala neurons [65][66][67] . Also, the rapid release of catecholamines 68 which may act upon distributed regions of the cortex in a less input selective manner compared to excitatory glutamatergic inputs 21 . The importance of this 'positive control' study is that it offers initial evidence that changes observed in response to the familiar context stimuli and the novel social stimulus are not reproduced by another arousing or stress-inducing stimulus, We should note that these animals did have stable respiratory rates and arterial oxygen saturation 23 . Therefore, systemic physiological instability is unlikely to have caused the increases in overall network strength.
We should note several limitations of the present study. First, sedatives and anesthetics exert varied effects on functional connectivity and BOLD signal activity based on their underlying molecular targets 69 . These include the magnitude of the BOLD signal, the strength of correlations 40 , the specificity of connectivity patterns, and these are related in part with the direct actions of isoflurane in blood vessel smooth muscle cells and small arteriolar vessel relaxation as well as the neuronal actions of this anesthetic 70 . Use of dexmedetomidine is currently more widely accepted than use of isoflurane but it is no less active in altering baseline conditions via actions on cerebrovascular parameters and neuronal activity mechanisms. We have conducted awake rat imaging in the past in the context of cue stimulated reward associations 71 , but the experience related to awake acclimatization involving repeated restraint sessions and repeated isoflurane inductions would alter the background upon which the contextual and social stimuli were presented. We should, however, note that while animals are under low sedation and low level anesthesia during scanning, the stimulus is presented with a temporal delay prior to setup and scanning that would allow the generation of neuronal responses and preservation of some degree of network topology. A second limitation is the non-specific simplified stimulus presentation model and its loose association with memory mechanisms. With regards to the association of the current stimulus paradigm with memory, it is important to note that rats were exposed over 3 days to the same test cage and then placed into the same test cage for 10 min prior to scanning. This pre-scan experience resulted in robust changes in node strength and other connectome quantifiers relative to unexposed animals. The same outcome was observed in two groups of rats undergoing the same context-exposure regime, those without any additional stimuli and those also presented with a novel unfamiliar social stimulus rat. This supports the possibility that the effect of exposing rats to the pre-scan familiar context was reproducible within the same cohort of rats, both male and females. We should also note that the word 'memory' is used here to denote the preservation of a pattern of network topology in stimulus groups that is absent in control unexposed groups. The paradigm used here did not generate a memory 'response' by the rats prior to or after scanning and this weakens the direct linkage between the network differences between groups and their role in memory structure in the brain. Regarding the experimental paradigm itself, we have previously applied other validated system-selective paradigms to assess learning and memory that allowed a priori hypothesis testing of the brain regions involved in memory processing 18,19 . We anticipate future studies using similar paradigms but with the added modification that the memory-retrieving process or the memory-signaling behavior be carried out or measured during a short temporal delay prior to imaging. Finally, approaches to selectively modify the activity of neuronal circuits to measure differences between control and stimulus groups was not conducted here. Optogenetic and chemogenetic approaches along with direct deep brain electrical stimulation have been applied to study the function of hippocampal circuitry in relation to broader network topology has been studied previously 72 . Our objective in future studies will be to integrate variants of these approaches to quantitative assessments of memory networks in control and in aging or neurological disease rodent models.

Materials and methods
Subjects. Thirty adult female and male Long Evans rats (225-350 g on arrival) were obtained from Charles River laboratories (Wilmington, NC). An additional 12 juvenile female and male rats (100-150 g upon arrival; 6 per sex) were used as social stimuli prior to imaging. Rats were housed in sex-matched pairs in a temperature and humidity-controlled vivarium with a 12 h light cycle (lights on at 0700 h) and food and water provided ad libitum. All procedures received prior approval from the Institutional Animal Care and Use Committee of the University of Florida and followed all applicable NIH guidelines.
Experimental design. Rats were assigned to one of 3 experimental groups, which included a naïve control group (no pre-scan exposure; n = 11, with 5 males), a cage only group (context exposure group; n = 10, with 5 males), and a cage + social stimulus group (context/social exposure; n = 9, with 5 males). The context and context/ social exposure groups were acclimated for 30 min a day for 3 consecutive days prior to imaging to a clear plastic test cage (dimensions: 40 cm 3 ) surrounded along the exterior of its walls with a black foam padding to minimize light inside the test cage (Supplementary Figure 1A). Acclimation was carried under low noise and dim light conditions in the MRI facility room near the scanner. On the day of imaging a subset of context exposed rats were presented with an unfamiliar juvenile rat of the same sex inside the familiar test arena for 10 min (Supplementary Figure 1B). Exploratory interactions between rats were visually confirmed (Supplementary Figure 1C).
In a cohort of rats tested separately but in the same test cages we observed that test rats interact with the unfamiliar rat through repeated bouts of exploratory sniffing, grooming, and pursuing inside the test cage (Supplementary Figure 1D). Typically, during a second exposure 1 h later (which corresponds to the time in which rats are scanned) the test rats reduce exploratory behaviors towards the now familiar rat (Supplementary Figure 1D). In the present study, each rat was scanned 50 min after initial social exposure (Supplementary Figure 1E). For the context group, rats were exposed to an empty test arena (familiar context) for 10 min and imaged 50 min later (Supplementary Figure 1E). To control for effects of context and context/social exposure, we included a control group of rats imaged without any experimental manipulation other than the same imaging setup process that rats in the other two groups underwent (Supplementary Figure 1E). All handling of rats, acclimation, pre-scan exposures, and image data collection were carried out by the same investigator to minimize novelty stress effects. . Using this protocol, we previously reported that rats had stable breathing rates (45-65 beats/min), heart rate, and relative blood oxygen saturation for over 1 h post-dexmedetomidine 73 . Animals were placed on a custom plastic bed with a respiratory pad underneath the abdomen. Respiratory rates were monitored continuously, and rectal temperature was maintained at 36.5-37.5 ºC using a warm water recirculation system (SA Instruments, Inc., New York). For each rat we acquired a 10-min high-resolution T2 weighted anatomical scan followed by a 10-min functional magnetic resonance imaging (fMRI) scan without any stimulation (under 'resting state' conditions). Functional MRI scans were acquired at least 50 min after dexmedetomidine administration. A T2-weighted Turbo Rapid Acquisition with Refocused Echoes (TurboRARE) sequence was acquired with the following parameters: effective echo time (TE) = 41 ms, repetition time (TR) = 4 s, RARE factor = 16, number of averages = 12, field of view (FOV) of 24 mm × 18 mm and 0.9 mm thick slice, and a data matrix of 256 × 192 and 25 interleaved ascending coronal (axial) slices covering the entire brain from the rostral-most extent of the prefrontal cortical region and caudally to the upper brainstem/cerebellum. Saturation bands were applied around the brain to suppress non-brain (muscle) signal.

Magnetic resonance imaging.
Functional images were collected using a single-shot spin-echo echo planar imaging (EPI) sequence 74 with the following parameters: TE = 15 ms, TR = 2 s, 300 repetitions, FOV = 24 × 18 mm and 0.9 mm thick slice, and a data matrix of 64 × 48 with 25 interleaved ascending coronal slices in the same space as the T2 anatomical. A fat saturation pulse was used to suppress chemical shift artifact and phase correction parameters were acquired during receiver gain adjustment for Nyquist ghost correction. Respiratory rates, isoflurane concentration, body temperature, lighting, and room conditions were kept constant across subjects 18 .
Image pre-processing. All fMRI and anatomical image pairs were reoriented to LPI orientation prior to processing. Binary masks outlining rat brain boundaries were generated in MATLAB using Three-Dimensional Pulsed Coupled Neural Networks (PCNN3D) on the high resolution T2 Turbo-RARE anatomical scans and on fMRI scans 75 . N4 bias field correction 76 was used to remove B1 RF field inhomogeneities and reduce FOV intensity variations in T2 anatomical scans prior to alignment to a parcellated rat brain template 77 . The binary 'brainonly' masks were multiplied by T2 anatomical and fMRI scans to remove voxels outside the brain. The cropped T2 anatomical brain images were aligned to the rat brain T2 template using the FMRIB Software Library (FSL 6.0.3; https ://fsl.fmrib .ox.ac.uk/fsl/fslwi ki/) linear registration tool, FLIRT 78 . Linear transform matrices were stored and later applied to functional scans to align these with the same template space.
Resting state processing was carried out using in house UNIX bash scripts that sequentially called tools in Analysis of Functional NeuroImages (AFNI) 79 , FSL, and Advanced Normalization Tools (ANTs) 80  www.nature.com/scientificreports/ images. First, in AFNI we used 3dDespike to remove time series spikes, 3dTshift for slice timing correction, and 3dvolreg for motion and linear drift correction. Cropped fMRI scans were split into 300 individual images and the first in the series was linearly registered to the corresponding cropped T2 anatomical using FSL FLIRT and then in ANTs we used antsRegistrationSynQuick.sh to warp the fMR image to the same cropped anatomical. Linear and warping matrices were stored and applied to the remaining fMR images in the series. After alignment between the fMRI and corresponding anatomical scans was improved, linear registration matrices transforming fMRI scans to the T2 rat brain template space was applied. Finally, fMRI scans were merged into a single image time series file. In the template space, we regressed white matter and CSF (ventricle) signals, applied bandpass filtering to time series images (0.009-0.12 Hz) and spatial smoothing (0.5 mm FWHM), and voxel time series L2 normalization 79 .
Extraction of spontaneously fluctuating fMRI signals for functional connectivity analysis. A total of 144 region of interest (ROI) masks, divided into 72 left and 72 right ROI's, were included in our analyses. Individual ROI masks were generated using a previously published rat brain parcellation and its accompanying T2 weighted image template 77 . To generate the ROI seeds we first produced left and right hemispheric versions of the parcellations, split these into individual left and right image files, extracted center of mass local coordinates for each image, and at these spatial locations we created images with a 0. Application of graph theory algorithms to analyze brain wide functional connectivity networks. Functional connectivity graphs were analyzed using Brain Connectivity Toolbox 7 . Pairwise correlations were reorganized into symmetrical graphs containing 144 nodes and 10,152 edges. Graph edge weights were normalized to a maximum value of 1 (while retaining their normal distribution). To avoid any bias in the calculated network measures, due to inter-subject variations in edge densities, we applied a threshold that maintained an equal edge density across subjects. A percentage of top ranked z values were retained and z values under the set threshold were zeroed. Using this procedure, network measures were calculated for 9 different edge density thresholds ranging from 5-45% (in 5% steps). Following normalization and thresholding, graphs (G) were binarized such that Gij = 1 indicates the existence of a 'connection' between any nodes i and j. For each node, the sum of above-threshold connections in half of the binarized graph provides the node degree (total nodal connections). For non-binarized (weighted) graphs, the sum of above-threshold edge weights (z values) provides the node strength (strength of total nodal connections). All network calculations were carried out using weighted undirected graphs. Also, for all brain functional network quantifiers we also estimated these for randomly rewired matrices (details below).

Small world coefficient.
To assess network integration and efficiency, we first analyzed the clustering coefficient (CC; an index of the average number of connected neighbor pairs j k of a node i that are connected to each other) 82 , characteristic path length (L; the lowest average number of edges between any pair of nodes i and j in graph G), and the small world coefficient (σ; which is > 1 for small world networks with CC greater than random networks and L near random) 6 . For L, we first substituted matrix weights by their inverse values such that highest z values correspond to shortest distances. Global efficiency was determined as the average inverse of L. To calculate σ, we randomly swapped (10 times) all edges in each normalized weighted graph (Gorig) to generate a randomized version of the original graph (Grand) that preserved degree and strength distributions 7,83 . CC and L were then determined for original and randomized graphs and σ calculated as CCorig/CCrand/Lorig/Lrand 84 .
Assortativity and rich club index. We next analyzed the tendency of assortative vs. dissortative mixing of nodes 85 . The assortativity index is a Pearson correlation coefficient comparing node strength values between pairs of edge-connected nodes. Positive r values indicate connectivity between pairs of nodes with similar strengths( e.g., high strength nodes pairs with high and low with low), while negative r values indicate cross-pairings between low and high strength nodes. We also generated and analyzed weighted rich club coefficient curves (Φ) at an edge density of 10% 86 . Nodes belonging to a rich club subnetwork have an above-chance tendency to tightly connect with each other and share the bulk of connections within the network 86,87 . The approach creates subgraphs containing nodes with strength values at a predetermined degree value, k. For each k-level subgraph, Φ w is calculated as the ratio of the sum of edge weights that connect nodes of degree k or higher (W > k) to the total possible number of summed edge weights 87 : www.nature.com/scientificreports/ Modularity and node cartography. Functional connectivity has been previously reported to follow a modular organization in which subgroups of nodes show a greater degree of intra-group functional connectivity than inter-group connectivity 31 . To assess the effect of context exposure and social interaction on this segregated pattern of activity we calculated a probabilistic version of the modularity index 88 . The agglomerative procedure starts with a random grouping of nodes and iteratively moves nodes into groups which maximize a modularity index, Q. The final module assignments for each node (e.g., a community affiliation vector) was taken as the median of 1000 iterations of the modularity maximization procedure. The role each node plays within their assigned module can vary according to their within-module degree (relative to other nodes in the same module) and their level of participation in connections with nodes in other modules. To investigate the role of intramodule nodes, we analyzed the within-module degree z-score (z) and the participation coefficient (P) 29 . A node i in module mi may play a central role if it has a within-module node degree (d i ) greater than the average node degree (δ) for that module, The difference between individual intra-module node degree and average intra-module node degree is variance-normalized to produce a z-score which allows the assignment of connector hubs at zi > 2.5 (zi > 1 in functional connectivity networks) 29,30 . The participation coefficient is calculated as: where dim is the within-module node degree and di is the overall network degree for a node I 29 . Within-module nodes with different levels of participation in connections with nodes in other modules were classified as follows: P < 0.05 are ultra-peripheral, 0.05 < P < 0.62 are peripheral nodes, 0.62 < P < 0.80 are non-hub connectors, and P > 0.80 are non-hub kinless 29 .
Centrality and network robustness. Functional connectivity networks are thought to contain high degree 'hub' nodes that either share connections with a large proportion of nodes, connect with other high degree nodes, and/or have high transit of connections between pairs of nodes in a network 89 . To measure the effect of context exposure and social interaction on putative hub node configurations we calculated betweenness and eigenvector centrality. A node i has high betweenness centrality if it lies within many shortest paths between any nodes j k (relative to all shortest paths between j k) 89 . Eigenvector centrality considers not only the number of connections node i possesses (quantity) but places greater weight on how many of these are with high strength nodes (quality) 31 .
In addition, we tested network robustness at an edge density of 10% using a random and preferential node removal strategy, as previously proposed 5,35 . For random node removal ('random failure'), we used a random permutation approach that iteratively removed individual nodes, and at each removal step calculated the size of the largest component. The size of the largest connected component (S) was averaged over 1000 iterations and proportionally scaled between 0 and 1. For preferential node removal ('targeted attack'), nodes are first ranked from highest to lowest values for node degree, strength, betweenness centrality, or eigenvector centrality, and then removed in a step wise manner while calculating S at each step. We compared the results of random and targeted attack in rat fMRI networks with synthetically generated random and ring lattice networks with the same number of nodes and edge densities.

Functional connectivity network visualizations. Functional connectivity networks were visualized in
BrainNet viewer 81 . Center coordinates for each ROI were generated in BrainNet viewer based on parcellations on a previously published rat brain template 77 . A 3D whole brain surface mesh file of the rat brain template (in *.byu format) was generated using an image binarization command in FSL (fslmaths) 78 and mesh construction in ITKSNAP 90 . Several 3D brain maps were generated. A first series of connectivity strength maps were produced with the size of nodes (spheres) weighted by node strength and the edges (lines connecting node pairs) weighted by the Pearson coefficient between pairs of nodes (lower bounds threshold set at z = 0.15). Additional maps were generated with node color and edges representing the module assignment (e.g., community affiliation vector) and node size weighted by either the modularity index, local efficiency, betweenness, or eigenvector centrality scores.
Statistical analyses. Statistical analyses of network measures and data plotting were carried out in MAT-LAB. For local network measures at a single graph density threshold, one way analysis of variance (ANOVA) was used with pre-scan manipulations as independent variables (p values were adjusted using a false discovery rate correction method). For global network measures at 9 different graph density thresholds we used a repeated measures ANOVA. Where appropriate, Tukey-Kramer tests were used for post hoc multiple comparisons (p < 0.05).
Functional MRI reanalysis: comparing results with an alternative pre-scan intervention. We reanalyzed a subset of control scans that were previously collected as part of a different set of studies. This additional analysis was carried out to determine the degree of specificity of the present results by assessing whether a different intervention prior to imaging exerted a similar or different effect on resting state functional con- www.nature.com/scientificreports/ nectivity and connectome quantifers. We determined the effect of a single saline injection administered either 1 h (n = 6) or 24 h (n = 10) prior to fMRI scanning 23 and compared these two groups to a naive group of rats scanned with no pre-scan intervention (n = 8) 91 . Saline injections were applied over a lower abdominal quadrant using sterile packaged tuberculine (1 cc) syringes and 26 gauge needles, at a 45° angle into the intraperitoneal layer. Following injections, rats were returned to their home cages and scanned according to previously reported methods 23,91 . All image processing, network calculations and statistical analyses and plotting were carried out as described in the previous sections.