Introduction

High-throughput mapping of IEGs, such as c-Fos or Arc, has been immensely useful in relating complex behaviors in rodents to their underlying brain activation (i.e., ensemble) patterns. The advent of effective intact tissue clearing methods (Chung et al., 2013; Renier et al., 2014; Susaki et al., 2015; Ueda et al., 2020), microscopy advances (Ragan et al., 2012; Reynaud et al., 2015), and automated registration and cell segmentation pipelines (Renier et al., 2016), have made brain-wide activity mapping projects both scalable and feasible. Recent studies have applied these methods to study brain activation during a wide variety of experiences, including after psilocybin, ketamine, and haloperidol administration (Davoudian et al., 2023; Renier et al., 2016), the incubation of palatable food craving (Madangopal et al., 2022), and fear memory (Franceschini et al., 2023). However, with a few exceptions, almost all of these approaches have focused on mapping single ensembles across the brains.

Over the last decade, activity-dependent tagging approaches, such as the tamoxifen-inducible ArcCreER^T2 (Denny et al., 2014) and TRAP1/TRAP2 systems (DeNardo et al., 2019; Guenthner et al., 2013), or the doxycycline-controlled TetTag system (Reijmers et al., 2007), have allowed for the investigation of multiple ensembles—and their shared overlap—during two distinct experiences within the same subjects. Such tagging methods have generated unprecedented insight into the neural mechanisms of memory formation, storage, and retrieval by allowing for the identification of individual memory traces. A memory trace or “engram” (Semon, 1921) refers to the enduring physical substrate of memory, and is defined as a neural population initially active during memory encoding whose reactivation at a later timepoint results in memory retrieval. However, the majority of engram studies have focused on a few focal regions, such as the dentate gyrus (DG; Denny et al., 2014), retrosplenial area (RSP; Cowansage et al., 2014), medial prefrontal cortex (mPFC; Kitamura et al., 2017), and amygdalar subdivisions, such as the lateral amygdalar nucleus (LA), and basolateral amygdalar nucleus (BLA; Reijmers et al., 2007), rather than the investigation of numerous regions across the brain. In part this is because high throughput mapping and analysis of multiple ensembles poses many challenges: 1) Both activity markers need to be clearly immunolabelled across all brain structures under investigation; 2) both labelled markers need to be accurately segmented or automatically counted; 3) a validated strategy needs to be employed for identification of overlapping cell populations; 4) all cell counts and co-labeled cells need to be accurately registered to a standardized atlas space; 5) an infrastructure facilitating easy data management, aggregation, and transformation (e.g., normalizing cell counts by volume) needs to exist; and 6) analytical and visualization methods to intuit functional connectivity between brain regions and how they differ between experimental groupings should be easily implemented.

Here, we developed a workflow incorporating these features to make large-scale analysis of traditionally coronally-sectioned datasets produced by activity-dependent tagging systems more accessible. Intrinsic to this pipeline is SMARTR, a versatile R package that wraps mapping capabilities to the common coordinate framework (CCF) of the Allen Brain Institute (Wang et al., 2020) with functions for statistical analysis and network visualization for multiple ensembles using graph theory. Moreover, this package is versatile, such that the analysis modules can be easily applied to imported mapped datasets from other alternative imaging, segmentation, and registration workflows, including those reliant on alternative atlas ontologies, such as the Unified Anatomical Atlas (Chon et al., 2019). While we have designed and validated this pipeline with the ArcCreER^T2 tagging system (Denny et al., 2014), the steps and analysis software used in this pipeline are generalizable to virtually all widely used dual-ensemble tagging systems.

Beyond its direct application to engram research, this pipeline has broader use in allowing for the comparison of brain-wide activity patterns during any two distinct experiences. To demonstrate the versatility of our workflow, we applied this approach to gain novel insight into brain-wide functional activation patterns underlying the acquisition and expression of LH, a robust behavioral paradigm which models traumatic stress. To our knowledge, wide-scale mapping of ensembles underlying both of these timepoints has not been studied in the context of LH and would allow for greater functional insight into changes occurring between these stages. In this effort, we have identified a number of altered functional connectivity signatures of the ensembles active during exposure to IS, escapable shock (ES), as well as the shared ensemble population between these experiences. Finally, we have presented this work as a reference dataset freely available for download and walkthrough alongside an accompanying tutorial.

Results

Pipeline overview

Our pipeline was designed for mapping and analysis of dual-ensemble datasets generated by the ArcCreER^T2 tagging system (Denny et al., 2014), although the steps are, in principle, generalizable to all dual-ensemble tagging approaches. The ArcCreER^T2 mouse line leverages the activity-dependent expression of the IEG Arc to enable permanent fluorescent tagging of neural ensembles with enhanced yellow fluorescent protein (eYFP) (Fig 1A) (Denny et al., 2014). In this system, injection of a selective estrogen receptor agonist, such as 4-hydroxytamoxifen (4-OHT), opens a transient window for labelling recently active neurons with eYFP. Staining for IEG protein expression (i.e., c-Fos) at a second timepoint allows for the identification of active ensembles driven by two different behavioral experiences (Fig 1B). Following tissue perfusion and brain extraction, multiple ensemble mapping through the pipeline workflow was organized into three different components: 1) immunolabelling and imaging acquisition, 2) image preprocessing and automatic segmentation of dual ensembles and their overlap in ImageJ/FIJI, and 3) atlas registration and downstream analyses in our software package in R (Fig 1C).

Fig 1.

For the immunolabelling and imaging acquisition, we employed an optimized approach for staining Arc-expressing eYFP⁺ cells and c-Fos⁺ cells in thick serial coronal sections based on our previous work (Leal Santos et al., 2021; Pavlova et al., 2018). A traditional serial sectioning approach was used, as it is the most accessible and reliable way to stain two markers at a high signal-to-noise ratio across deep brain structures.

Following imaging for eYFP⁺ (channel 1) and c-Fos⁺ (channel 2) cells using a scanning confocal microscope (Leica TCS SP8, Leica Microsystems Inc., Wetzlar, Germany) (Fig 1D), we implemented a custom segmentation strategy in ImageJ/FIJI which uses a suite of open-source image processing plugins (Ollion et al., 2013) for accurate cell segmentation in 3D based on the staining morphology characteristics of eYFP⁺ and c-Fos⁺ cells (Fig 1E). The characteristics and positions of segmented eYFP⁺ and c-Fos⁺ cells were saved into .txt files as outputs. Next, co-labeled (eYFP⁺ / c-Fos⁺) cells were identified by comparing all possible overlapping objects using the open-source 3D MultiColoc plugin (Ollion et al., 2013), saving the output as a .txt files, and thresholding by percent volume overlap downstream. The segmentation and colocalization algorithms were efficiently applied in-batch using custom written ImageJ/FIJI macros that are freely available for download, along with example images for trial (https://osf.io/ynqp7/). Importantly, these macros have been generalized to allow for easy parameter modification through a GUI, to account for different imaging resolutions and unique naming conventions.

In the final preprocessing step, since registration is performed on single-plane .tiff images, we collapsed raw images across all z-planes and channels to use as an aggregate reference image for registration to atlas templates from the Allen Institute Mouse Brain Atlas (Lein et al., 2007; Oh et al., 2014). This step was also batch processed through a generalized macro in ImageJ/FIJI that is available for download.

The SMARTR package for registration, analysis, and visualization of multiple ensembles

The SMARTR package (simple multi-ensemble atlas registration and statistical testing in R) includes a high-level application programming interface (API) which facilitates interactive user-friendly registration, mapping of cell counts, and easy statistical analysis and visualization of multiple ensemble datasets (Fig 1F). This standard API was designed around an object-oriented infrastructure. In the package, data is managed using slice, mouse, and experiment objects which also store additional metadata in a hierarchical manner. For example, a slice object contains all imaging-related metadata such as the image ID, image path, channel names, and a matching anterior-posterior (AP) atlas coordinate, alongside registration data and imported segmentation data. Slice objects are hierarchically stored inside mouse objects, which contain additional metadata such as mouse ID, sex, strain, and genotype. Package functions for interactive registration, importation of segmentation outputs, mapping cell counts to a standard atlas space, automated network analysis, and visualization operate directly on these objects. Moreover, there are also a suite of additional capabilities, including data quality checking for outlier counts, removal of regions to omit due to tissue damage, and normalizing co-labeled counts by a particular channel. The organization of the SMARTR package also facilitates alternative pipeline entry points. Import functions exist for alternate segmentation approaches. Compatibility with datasets previously registered and mapped with alternative workflows (Eastwood et al., 2019; Jin et al., 2022; Renier et al., 2016; Yates et al., 2019) is also possible (Fig 1F), allowing users to take advantage of the analysis and visualization capabilities of this package without reliance on prior processing steps. Extensive documentation of the package functions and a detailed tutorial using the datasets from this study are available at package website (https://mjin1812.github.io/SMARTR/).

Identifying ensembles underlying the acquisition and expression of learned helplessness

To demonstrate the strengths and capabilities of our dual-ensemble mapping pipeline, we applied our approach to gain novel insights into the neural activity patterns underlying LH, a robust behavioral paradigm widely used to model aspects of trauma and depression preclinically. The neural correlates of LH have been studied across a variety of model organisms (Maier & Seligman, 2016), and have been largely focused on single or few regions. More recently, there have been investigations into the brain wide neural activity patterns following LH expression (Kim et al., 2016). However, to our knowledge, the neural ensembles underlying both the acquisition and expression of LH have not been studied. This approach naturally allows for study of the overlapping active populations between acquisition and expression of LH, conferring greater insight into functional changes occurring between these stages.

To label active ensembles with eYFP during LH acquisition, we injected ArcCreER^T2 mice with 4-OHT 5 h prior to exposure to a series of inescapable footshocks on one side of a shuttle-box (IS, Inescapable Shock group). (Fig 2A). A second group of ArcCreER^T2 mice underwent an identical procedure without experiencing footshocks (CT, Context Trained group). Following 3 days of dark housing, mice underwent 30 trials in which they could escape footshocks by crossing the opposite chamber of the shuttle box. Ninety minutes following the start of testing mice were sacrificed and brains were extracted for investigation of c-Fos expression.

Fig 2.

During a 3-min habituation period prior to IS testing, locomotor activity was significantly decreased in the IS group (Fig 2B), indicating a fearful association was established with the shuttle box chamber. Across all trials during testing, escape latency was consistently higher in the IS group compared to the CT group (Fig 2C), and average latency to cross for trials 11-30 of the testing phase was significantly greater for the IS group compared the CT Group (Fig 2D).

After sectioning and immunostaining, we targeted our imaging of eYFP⁺ and c-Fos⁺ cells (Fig 2E-H) across whole coronal sections containing the following regions: hippocampus (HPC), amygdala (AMY), insular cortex (INS), entorhinal cortex (ENT), and raphe nuclei (RAmb). With this semi-targeted mapping approach, we ensured that critical regions previously implicated in the expression of LH and stress responses (Ineichen et al., 2022; Kim et al., 2016; Maier & Watkins, 2005; Marques et al., 2022) would be represented in our dataset, along with a myriad of under-investigated regions present within the same atlas plates spread across the rostral-caudal axis. To account for the dendritic/axonal staining pattern of eYFP, a custom algorithm was applied to limit segmentation of cellular processes (Fig 2I). Since this staining pattern is characteristic of other tagging approaches, e.g., TRAP1/TRAP2 (DeNardo et al., 2019; Guenthner et al., 2013), this algorithm is likely generalizable to other activity-dependent tagging mechanisms. Since c-Fos localizes to the cell body, resulting in a punctate filled-spherical staining, a separate optimized algorithm for segmenting this expression pattern was applied (Fig 2J). Next, all possible segmented object overlaps were compared and the proportion of volumetric overlap were calculated (relative to the volume of segmented c-Fos⁺ cells). Colocalized cell counts were identified by thresholding by proportion of volume overlap downstream using SMARTR (Fig 2K). The individual image processing steps for both segmentation algorithms and colocalization analysis were all conducted in-batch using ImageJ/Fiji macros, with details extensively outlined in the Methods. These segmentation parameters are adjustable in a graphical user-interface (GUI) menu to flexibly account for a range of resolutions, making then generalizable to other imaging parameters.

Next, we sought to validate the accuracy of our algorithms by comparing the automatically segmented counts of eYFP⁺ and c-Fos⁺ cells to manual cell counts in a subset of images across all mice. Specifically, automated cell counts in the granule cell layer of the DG of the HPC were compared to manual counts from two independent counters. For eYFP⁺, we found that automated segmentation results were highly correlated with manual cell counts across both annotators (Fig S1A-H, r = 0.85 averaged manual counts), and there was high correlation between independent manual counts (r = 0.92). Performance was further evaluated by independent quantification of false positive (FP, detected by algorithm alone), false negative (FN, detected by annotator alone), and true positive (TP, detected by both annotator and algorithm) cells on 3333 manual counts from Counter A and 3687 manual counts from Counter B (Fig S1I). We found that average precision (P = TP/(TP+FP)) was 0.92 and average recall (R = TP/(TP+FN)) was 0.76. The averaged overall F1 score, calculated as the harmonic mean of precision and recall (1/F1 = (1/P + 1/R)/2), was 0.85. For c-Fos⁺ cells, we similarly found high correlation between automated cell counts and manual cell counts across both annotators (Fig S1J-Q, r = 0.93 averaged manual counts), as well as high correlation between annotators (r = 0.93). We further evaluated the number of FP, FN, and TP cells on 2971 manual counts from Counter A and 3477 manual counts from Counter B (Fig S1R). The c-Fos-optimized algorithm yielded an average precision of 0.96, an average recall of 0.78, and an average F1 score of 0.86. Overall, these results indicate that both algorithms faithfully reflected results that would be obtained from typical manual cell counting approaches and even outperforms metrics reported previously for other automatic segmentation approaches for activity-dependent labelled cells (Franceschini et al., 2023).

Registration and mapping cell counts to a standardized atlas space

Using the SMARTR package, we created mouse and slice objects and auto-populated metadata using a custom script that is available on the package website. We then aligned the preprocessed flattened images (Fig S2A) to the Allen Mouse Brain atlas templates using the registration() function in SMARTR. This function interfaces directly with the WholeBrain (Fürth et al., 2018) and SMART (Jin et al., 2022) packages for interactive user-friendly registration improvement (Fig S2B). We used functions in SMARTR to import all raw external segmentation data, including cell counts from both channels as well as co-labeled cell counts (Fig S2C, c-Fos⁺ cells shown). Segmentation data and registration were next integrated to project counts onto a standardized atlas space (Fig S2D, c-Fos⁺ cells shown). Using package functions for data transformation and quality checking, cell counts per region were aggregated across all slice objects per mouse and normalized by volume of the regions mapped. Outliers were cleaned by dropping regional counts greater or less than two standard deviations from their group mean. Only regions represented across both experimental groups per channel were analyzed. We then applied our suite of analysis and visualization functions to examine the brain-wide activity patterns of eYFP⁺, c-Fos⁺, and reactivated (co-labeled / eYFP⁺ cells) cells.

Learned helplessness acquisition is marked by enhanced sensory and affective processing

We first investigated neural activity differences (eYFP⁺ cells) globally and in selected regions during LH acquisition between CT and IS groups. These included targeted isocortical regions (anterior cingulate area, ACA; agranular insula, AI; RSP), dorsal HPC, ventral HPC, and amygdalar subregions (see Table S3 for all regional acronyms used). Previous studies have implicated these isocortical regions in LH (Bauer et al., 2003; Kim et al., 2016), and the ventral HPC and amygdala serve as core limbic centers, especially in fear and aversive processing. Interestingly, during LH acquisition, global brain-wide neural activation levels did not differ between groups (Fig S3A). Additionally, there were no differences in overall activity amongst all targeted brain regions between groups (Fig S3B-E; see Fig S10 for a brain-wide region activation comparison). Correlation of neural activity with behavioral performance in each experimental group was weak across all isocortical, hippocampal, and amygdalar subregions analyzed (Fig S3F-S), with the exception of the RSP, in which the IS shock group showed strong correlation between eYFP⁺ cell counts and mean escape latency (r = 0.88, p = 0.02).

Next, we examined the cross-correlations in IEG expression across brain regions, as strong co-activation or opposing activation can signify functional connectivity between two regions. Correlations were calculated for eYFP⁺ region counts normalized by volume across all pairs of brain regions for the CT and IS groups, respectively (Fig 3A-B). In both networks, there was a higher proportion of positive correlations in cortical regions. However, the distribution of most significant cortical correlations (p < 0.01) was predominately in motor and somatosensory areas in the IS group, whereas in the CT group it was biased towards the auditory areas and integrative or associative regions, such as the ectorhinal cortex (ECT).

Fig 3.

To better visualize the most relevant co-activations and analyze their topological structure, we thresholded to retain only the strongest (r > 0.9) and most significant (p < 0.01) correlations and constructed co-activation networks, in which connections are represented as edges bridging two nodes or brain regions (Fig 3C-D). We then calculated different topology metrics per node, including degree, clustering coefficient, efficiency, and betweenness centrality(Bullmore & Sporns, 2009). Degree is the simplest measure of connectivity, and represents how many direct connections a single node has. The clustering coefficient represents the tendency of finding a third mutual connection amongst already connected neighboring nodes (i.e., “cliqueishness”), giving insight into how functionally segregated or integrative a given region is amongst its neighbors. Efficiency is inversely related to the average path length between a given node and all other nodes in the network, and betweenness indicates how often a particular node lies on the shortest path between two other nodes, which can indicate a region’s influence over information flow in functional networks.

We first averaged each of these metrics across all nodes to gain insight into each network’s global connectivity patterns (Fig 3E-I). While there was no difference in average degree between groups (Fig 3E), both networks showed characteristic tailed degree distributions that tend to be found in complex biological systems that follow a power-law distribution (Fig 3F) (Barabási & Albert, 1999; Bullmore & Sporns, 2009, 2012; Watts & Strogatz, 1998). There was also no difference in average clustering coefficient, global efficiency, and mean betweenness centrality (Fig 3G-I). As currently there is no consensus on the ideal approach to choose a threshold to optimize retention of important connections and rejection of spurious ones (Garrison et al., 2015), we next examined whether the calculated global metrics might have changed across a wide range of significance thresholds (Fig 3J-M). Trajectories of degree, clustering coefficient, efficiency, and mean betweenness centrality were qualitatively similar between CT and IS groups across a range of alpha thresholds. Thus, global topology metrics reveal no striking differences in propensity for information transfer between IS and CT networks during LH acquisition.

To better investigate which individual regions are most differentially influential between networks, we looked specifically at the distribution of individual nodal degree centrality and betweenness centrality values (Fig 3N-O). The degree distributions revealed that many primary somatosensory (SSp) regions (mouth, SSp-m; lower limb, SSp-Il; barrel field, SSp-bfd; upper limb, SSp-ul), motor regions (primary motor area, MOp) and limbic regions, including the nucleus accumbens (ACB) and BLA, were highly connected in the IS. Comparatively, far fewer sensory regions showed such high degree in the CT networks. Despite the high connectivity amongst mostly sensory regions, IS network betweenness distributions showed that additional key structures such as the AI, piriform cortex (PIR), septohippocampal nucleus (SH), and ENT act as influential information bridges in the IS network. The AI is a critical region involved in pain processing (Labrakakis, 2023), and its central network engagement is largely consistent with experiencing a painful sensory event. In contrast, structures such as the subiculum (SUB), ECT, dDG, and medial septal nucleus (MS) were most influential in CT networks, and betweenness distributions also reveal that the SUB was most influential. The subiculum is well-known to play a role in spatial representation during exploration (Matsumoto et al., 2019; O’Mara et al., 2009), and the MS and dDG are respectively involved in motivated locomotion (Mocellin & Mikulovic, 2021) and contextual processing (Tuncdemir et al., 2019). The overall patterns seen in the individual regional degree and betweenness distributions indicate there is a strong engagement of areas involved in contextual and navigational processing in the CT network, whereas the appears to be greater influences of regions involved in sensory and affective processing of shock (e.g., BLA, AI), and motor reactivity (e.g., MOp) in the IS network.

Learned helplessness acquisition alters ventral CA1, ventral tegmental, and perirhinal connections

To better understand how coordinated activity amongst individual connections differ between groups, we performed a permutation analysis. Regional correlations from the CT group were subtracted from corresponding correlations in the IS group, and compared against a permuted null distribution. We further analyzed only correlations differences that were large (r difference >=1) and highly significant (p < 0.01) (Fig 3P). A number of interesting patterns emerged from this analysis. Specifically, connections involving either the ventral CA1 (vCA1), ventral tegmental area (VTA), or perirhinal cortex (PERI) were strongly positively correlated in the CT group, whereas they became negatively correlated in the IS group (Fig 3Q). This pattern was consistent amongst connections between the vCA1 with associative areas (posterior parietal association areas, PTLp; ENT) and hypothalamic areas (zona incerta, ZI; lateral hypothalamic area, LHA), and between the PERI with primarily limbic, especially amygdalar regions, (AI; BLA; piriform amygdalar area, PAA; cortical amygdalar area, COA), and the VTA with associative regions (temporal association areas, TEa; ECT, PERI). The, VTA is centrally involved in reward processing and functional inversion of its connections may be indicative of the aversive experiencing of inescapable footshocks. Additionally, the vCA1 is known to be involved in processing of affective states and fear retrieval, with the vCA1-LHA projection in particular being enriched in representations of anxiety (Jimenez et al., 2018, 2020), suggesting functional alterations of this particular connection may be important in the development of LH. Finally, the PERI region is involved in object perception and sensory feature processing (Murray & Richmond, 2001), and thus, its strongly altered connectivity may indicate differences in how sensory features of the environment are assigned an affective value in IS versus CT mice.

Learned helplessness expression is characterized by a reversal in functional connectivity of the substantia nigra

We applied the same analysis approach to the mapped c-Fos⁺ dataset, representing cells active during the expression of LH, when mice are tested under ES conditions. Global brain-wide neural activation levels during LH expression did not differ between groups (Fig S4A). Additionally, in the ACA, AI, RSP, dorsal HPC, ventral HPC, and amygdalar subregions (Fig S4B-E), there were no differences in overall activity between groups. Correlation of activity with behavioral performance in each experimental group was weak across all isocortical, HPC, and amygdalar subregions analyzed (Fig S4F-S). Thus, examining only absolute activation patterns did not reveal obvious neural signatures of expression of learned helplessness (see Fig S11 for a brain-wide region activation comparison).

We next examined the functional network signatures. As was performed for eYFP⁺ networks, we calculated regional cross-correlations (Fig 4A-B), and thresholded to retain only the strongest (r > 0.9) and most significant (p < 0.005) correlations to construct networks (Fig 4C-D). Qualitatively, we observed that somatosensory and motor areas appear more integrated with amygdalar and hippocampal structures in CT compared to IS networks (Fig 4C-D). We then calculated the same global network topology metrics previously described. We found that there was no mean difference in degree (Fig 4E), although both networks showed the expected tailed degree distribution (Fig 4F). There was also no difference in clustering coefficient, global efficiency, and mean betweenness centrality (Fig 4H-I). We validated the stable nature of this observation by recalculating these global metrics across a wide range of possible significance thresholds, and observed similar trajectories between CT and IS networks across all metrics (Fig S5A-D). Similarly, to our eYFP⁺ analysis, these findings indicate that the expression of LH is not characterized by global differences in efficacy of information transfer in functional networks.

Fig 4.

To investigate which key influential regions differ between in each network, we plotted the individual region degree and betweenness distributions for CT and IS groups (Fig 4J-K). The most highly interconnected region in the CT network was the nucleus reunions (RE) and region with the highest betweenness value was the diagonal band nucleus (NDB). In contrast, for the IS network, the most highly interconnected regions include the dCA1, and ventral medial nucleus of the thalamus (VM), whereas two basal ganglia structures, the subthalamic nucleus (STN), and the substantia nigra, reticular part (SNr) displayed the highest betweenness centrality. To further investigate how regions showing either the highest degree or betweenness centrality in the IS group are functionally changed, we plotted the distribution of Pearson correlation coefficients for all dCA1, VM, STN, and SNr connections (Fig 4L-O). While there were no differences in the correlation distributions of the dCA1, VM, and STN, in the SNr, there was a strong shift from predominantly negative correlations in the CT group to positive correlations in the IS group (Fig 4O). Altogether, these findings indicate that the VM, a central thalamic relay structure which integrates limbic and prefrontal cortical signals to guide motor behavior (Anastasiades et al., 2021; Monconduit et al., 1999; Sieveritz et al., 2019), and the NDB, a basal forebrain structure capable of mediating motivated movement (Mocellin & Mikulovic, 2021), may be involved in the cognitive processes involved in generating adaptive escape responses in control animals. General reversal of functional connectivity in the substantia nigra when comparing the CT and IS suggest that the substantia nigra may be a critical region involved in gating of motor processing involved in an adaptive crossing response.

Finally, we performed a permutation analysis as previously described to examine pairwise regional correlation changes between groups (Fig 4P-Q). In line with our previous findings, we found a large proportion of the significantly altered connections in the positive direction in the IS group relative to the CT included the SNr. These SNr connections involved a number of sensory, motor, and visual processing regions. Conversely, a number of positively correlated connections involving the PERI region in the CT group shifted to negative correlations in the IS group. Among these connections were various midbrain structures (VTA, MRN) and visuomotor processing centers (SCm), and hypothalamic structures (ZI). As strong correlation reversals of PERI connections between CT and IS groups were similarly seen in ensembles active during learned helplessness acquisition, early modulation of PERI function may be an early functional signature of development of learned helplessness.

Learned helplessness is characterized by dampened reactivation of cells previously active during inescapable shock

We next leveraged the capability of the SMARTR package to study co-labeled eYFP⁺ / c-Fos⁺ cells populations active during both acquisition and expression of learned helplessness. To investigate the proportion of co-labeled cells reactivated amongst the original ensemble active during IS acquisition or context training, co-labeled cells were divided by eYFP⁺ counts to generate reactivation proportions (co-labeled / eYFP⁺ cells). We applied the same analysis approach to the proportion of reactivated cells, as was performed for c-Fos⁺ and eYFP⁺ networks.

Targeted regional analysis of reactivation activity (Fig 5A-E) revealed that a number of regions, including the ACA, AI, dDG, dCA3, and MEA, were significantly decreased in their reactivation proportions (Fig 5F-I). Interestingly, a general analysis of global reactivation proportions shows a trending brain-wide decrease in reactivated activity (p = 0.0560, Fig S6A; see Fig S12 for a brain-wide regional comparisons). Additionally, we found that reactivated activity in LA, a central region involved in associative fear learning (Johansen et al., 2010), showed a strong correlation with mean latency to cross in the IS group (Fig 5J). An additional exploratory analysis of reactivation activity as a proportion of the ensembles active during ES (co-labeled / c-Fos⁺ cells), revealed a similar trending global brain-wide decrease in reactivation proportions (Fig S7A, see Fig S13 for a brain-wide regional comparisons). Overall, these findings indicates that a general global dampening of reactivation activity may be a novel marker of the traumatized helpless state.

Fig 5.

Altered functioning of the insula and central amygdala in reactivated cells characterizes the learned helplessness state

We next investigated functional networks of reactivation activity by plotting regional cross-correlations (Fig 5K-L), thresholding to retain the strongest and most significant correlations (r > 0.9, p < 0.01), and constructing networks (Fig 5M-N). Similar to eYFP and c-Fos networks, examination of global topology metrics revealed no significant differences in degree, clustering coefficient, efficiency, and betweenness centrality (Fig S8A-E), and we confirmed the stability of this observation by seeing similar trajectories of these metrics between groups across a wide range of significance thresholds (Fig S8F-I).

We next plotted the individual region degree and betweenness distributions and found that the most highly interconnected regions in the IS network were the dCA1, AI, MOs, and the primary somatosensory area, nose (SSp-n). The regions with the highest betweenness centrality were the ACA and MOs. For this reason, we further plotted the distribution of all Pearson correlation coefficients for the dCA1, AI, ACA and MOs (Fig 5Q-T). While there were no differences found in the correlation distributions for the dCA1, ACA, and MOs, in the AI, there was a significant increase in proportion of positive correlations in the IS group compared the CT group (Fig R), indicating a general strengthening of functional connectivity of the AI amongst reactivated cells. In an exploratory analysis of co-labeled / c-Fos⁺ cell activity, we also found that AI correlations to the vDG region were preferentially strength in the IS group (Fig S9M).

Finally, we used a permutation analysis, as done previously (Fig 5U-V) to assess regional correlations differences that were large (r difference >=1) and highly significant (p < 0.01) between the two groups. Interestingly, a pattern emerged whereby there was a disproportionate number of strongly negative connections in the CT group that shifted to strongly positive connections in the IS group, many of which involved the CEA, a primary fear and pain processing center. These significantly altered CEA connections included regions such as the infralimbic area (ILA), which is involved in flexible response inhibition and decision-making (Barker et al., 2014), the caudoputamen (CP), a basal ganglia structure involved in motor regulation, and the lateral septum (LS), which is involved in emotional regulation and stress responding (Wirtshafter & Wilson, 2021; Yadin et al., 1993). Altogether, these observations suggest that a functional change in the AI and strong functional reversal of key CEA connections in reactivated neural populations may underlie the process of developing of LH.

Discussion

Dual labelling systems have been invaluable for studies investigating the shared overlapping ensembles active across two separate experiences. Such studies have revealed invaluable insights on the ensembles responsible for memory storage and retrieval i.e., engrams. However, most of these studies have relied on targeted investigation of single or few regions (Cowansage et al., 2014; Denny et al., 2014). The few publications reporting mapping of brain-wide tagged ensembles have relied on intact tissue clearing, 3D imaging, and registration workflows that, while advantageous in their ability to keep tissue intact, are specialized and difficult to implement widely, and/or rely on commercial software (e.g. Imaris) (Franceschini et al., 2023; Roy et al., 2022).

Meanwhile, only one paper, to our knowledge, has reported mapping and analysis of colocalized dual-labelled ensembles (Roy et al., 2022), without introducing a generalizable workflow. Thus far, there have been scant accessible approaches that allow for the scalable brain-wide mapping of multiple ensembles as well as their overlap in traditionally immunolabelled coronal sections. Moreover, there are no standardized tools open-source tools to analyze functional connectivity signatures based on network analysis of brain-wide immediate early gene activity patterns.

Here, we present our workflow for automated segmentation optimized for the identification of eYFP⁺ cells (indicative of previously Arc⁺ cells), c-Fos⁺ cells, and their co-labeled overlap. This process is open-source, as it is conducted entirely in ImageJ/FIJI (Schindelin et al., 2012), and relies on existing image processing toolkits that are freely available . Moreover, we introduce the R package SMARTR (simple multi-ensemble atlas registration and statistical testing in R) which offers a standardized API for the importation of segmentation output from our workflow, and for atlas registration and mapping by interfacing with previously published R packages (Fürth et al., 2018; Jin et al., 2022). However, perhaps the most useful and novel feature of SMARTR is the inclusion of numerous data management, quality checking, graph theory analyses, and visualization functions that are user-friendly, requiring only rudimentary programming knowledge.

We demonstrated the utility of our workflow by applying it to map two ensemble populations active during the experience of inescapable shock (or context) or later during exposure to escapable shock. Our goal was to gain insight into how functional activity networks are changed during both LH acquisition and expression, which has never previously been studied. As the LH paradigm models a fundamental symptomology of trauma-related depression and anxiety states, this work has the major clinical implication of improving our understanding of how these states emerge following trauma. Our analysis of global network properties of eYFP⁺, c-Fos⁺, and co-labeled / eYFP⁺ activity expression revealed no differences in topological metrics, suggesting that the learned helpless state is not characterized by gross deficits or alterations in information transfer during acquisition or expression of learned helplessness.

However, our region and connection specific analyses of eYFP⁺ ensemble networks suggests that, when experiencing inescapable shock compared to context exposure, there is a functional shift from emphasis on navigational and contextual processing to sensory and affective processing. Network-wise, this appears to be represented by a decreased influence of regions such as the subiculum, which is centrally involved with spatial navigation (Matsumoto et al., 2019; O’Mara et al., 2009), to generally high interconnectivity of somatosensory regions, and those potentially involved with affective processing of shock, such as the BLA and AI.

Surprisingly, we did not find strong evidence of absolute activity differences between groups across a targeted or brain-wide regional comparisons, and overall amygdalar activity in IS mice were comparable to that of CT mice during LH training. This is in contrast to previous findings showing Arc gene expression in the amygdala is upregulated 2 h after inescapable shock (Machida et al., 2018). However, a key methodological difference in our approach is that the ArcCreER^T2 x eYFP tagging strategy labels previously Arc-expressing cells with eYFP over a longer window of time, as the serum half-life of 4-OHT is approximately 6 h (Cazzulino et al., 2016). Thus, neurons active during consolidation of the experience in the hours following acute shock or context exploration may also be recruited into the tagged ensemble.

During escapable shock testing, network analysis revealed surprising high influence of two basal ganglia structures, the STN and the SNr, and that there was a striking functional reversal of the substantia nigra between the CT and IS groups. Interestingly, recent work has shown that stressor controllability can change intrinsic dynamics and electrophysiological properties of substantia nigra neurons (Yao et al., 2021). These findings, in conjunction with our results, suggest the possibility that the substantia nigra activity is altered by loss of control over a stressor in a manner that may gate an adaptive motor response such as crossing.

Finally, we examined the overlapping ensembles active during both learned helplessness acquisition and expression (eYFP⁺ / c-Fos⁺) as a proportion of the originally active ensemble population (eYFP⁺). Interestingly, we found a robust decrease in reactivation activity in a number of targeted brain regions, alongside a trending brain-wide decrease. To our knowledge, this pattern of widescale decreased reactivity following exposure to a traumatic event has not been described before, and suggests a potential phenomenon of neural plasticity which acts to suppresses future ensemble reactivity as a major part of the neurobiology underlying LH. We also found general correlative association between LA co-labeled activity and escape latency, suggesting reactivation of LA ensembles initially active during IS involved in regulating response to escapable shock.

Performing network analyses, we identified the AI as a highly interconnected region in the IS networks, which showed a significantly strengthened distribution of functional connections compared to the CT group. As an integrative hub for sensory, affective, and cognitive function, in addition to its role in pain processing, the insula is involved in regulating fear and anxiety behaviors (Casanova et al., 2016; Labrakakis, 2023; Nicolas et al., 2023; Rogers-Carter et al., 2018). The influential positioning of the AI in the IS networks, along with its functional alteration, indicate that reactivated functional activity in the AI may be critical in mediating the helpless state. Network analysis also revealed the ACA as another potential key region, as the ACA displayed the highest betweenness centrality in the IS network. This finding is consistent with a previous study of long-term fear memory using network analysis of c-Fos expression which also found the ACA as an important hub within the fear network (Wheeler et al., 2013). Finally, our permutation analysis revealed a number of key CEA connections that were strongly functionally reversed between IS and CT group, suggesting possible functional targets for future investigation. Given these many novel insights, we believe the strength of SMARTR and our multiple-ensemble mapping workflow lies in its use as a hypothesis-generating toolkit, which may lead to novel insights suggesting future mechanistic studies which may draw more definitive causal conclusions.

Methodological considerations

There are many factors to consider regarding the choice of our workflow over other options. While serial coronal sectioning is not as a rapid as intact tissue immunostaining, this approach allows for adequate antibody penetration for two or more ensemble populations in deep brain structures for high signal-to-noise staining. Intact whole-brain immunohistochemistry (IHC), typically in combination with lipid clearing (Chung et al., 2013; Park et al., 2019; Renier et al., 2014), has been demonstrated to label single ensembles well. However, insufficient penetration of multiple antibodies or purging of epitopes, especially during permeabilization steps in hydrophobic clearing approaches, remains a major technical barrier to widespread adoption in mapping engrams. We have found many of these methods result in confounding staining gradients from the surface to deeper depths of the tissue (Pavlova et al., 2018), leading to inaccurate visualization of multiple ensembles, which confounds study of biologically relevant overlapping populations. Applying intact tissue clearing approaches to brain wide engram mapping remains a tremendous technical hurdle for most.

A second limitation of this workflow is that accuracy of registration alignment depends on the quality of physical sectioning, user annotation, and correction of registrations. Furthermore, accurate registration in the z-plane is limited to a resolution of ∼100 μm. This is because the registration capabilities of Wholebrain (Fürth et al., 2018) and SMART (Jin et al., 2022), the underlying packages used for registration and mapping, require coronal sections to align closely to the sectioning plane of the Allen Mouse Brain Atlas. Thus, registration alignment requires both technical skill and familiarity with brain-wide anatomy. However, this challenge can be mitigated with use of brain matrices, which can help mount an even cutting plane during physical sectioning. Additionally, limited registration accuracy along the z-plane is still comparable to other recent volumetric registration approaches reporting accuracies of ∼300 μm (Franceschini et al., 2023), whereas maximal registration accuracy along the x- and y-planes is much higher due to Wholebrain’s scale-invariant spline-based atlas (Fürth et al., 2018).

Future directions

There are a number of alternative mapping approaches that can account for sectioning angle during registration (Puchades et al., 2019; Song et al., 2020). The structure of SMARTR is modular such that additional pipeline entry points are possible if separate workflows are used for segmentation, registration or mapping. For example, with some data formatting, independently mapped data can be easily imported for integration with SMARTR’s network analysis functions. Similarly, with some simple programming, it is possible to import segmented data from other approaches, such as CellPose (Stringer et al., 2021) or QuPath (Bankhead et al., 2017), and still take advantage of the built-in registration, mapping, and analysis capabilities in SMARTR. We anticipate, based on community feedback, that additional analysis features will be continually added to the package. For example, since the default supported atlas ontology is the Allen Mouse Brain Common Coordinate Framework, the merged Franklin-Paxinos (FP)/ common coordinate framework (CCF) atlas nomenclature (Chon et al., 2019) was recently incorporated to analysis and visualization functions based on collaborator feedback. Finally, while this pipeline was designed to be especially user-friendly to novice programmers, future iterations may be even easier to learn through incorporation of a GUI menu for easy creation and modification of object metadata, running analyses, and visualizations.

Methods

Mice

ArcCreER^T2 (Denny et al., 2014) x R26R-STOP-floxed-eYFP homozygous female mice were bred with R26R-STOP-floxed-eYFP homozygous male mice (Srinivas et al., 2001). All experimental mice were ArcCreER^T2 (+) and homozygous for the eYFP reporter. ArcCreER^T2 (+) x eYFP mice are on a 129S6/SvEv background, as they have been backcrossed for more than 10 generations onto a 129S6/SvEv line. The age of the mice utilized for experimentation ranged from 6 to 7 months.

Genotyping

Genotyping of Cre and eYFP was performed as previously described (Denny et al., 2014). APP and PSEN genotyping were performed as described on the Jackson Laboratory website.

4-Hydroxytamoxifen

Recombination was induced with 4-OHT (Sigma, St Louis, MO). 4-OHT was dissolved by sonication in 10% EtOH/90% corn oil at a concentration of 10 mg/ml. One injection of 200 µl (2 mg) was intraperitoneally administered to each mouse.

Learned Helplessness (LH)

LH was administered as previously described (Brachman et al., 2016). Briefly, in IS phase, mice in the IS group were habituated for 3 min to one side of a two-chamber shuttle box (model ENV 010MD; Med Associates, St. Albans, VT) containing a grid floor connected to a scrambled shock generator (model ENV 414S, Med Associates). The shuttle box was equipped with 8 infrared beams (IR, 4 on each side) for detecting position and activity of the animal, and an automated guillotine door separating the chambers. Training consisted of 140 shocks, each with a 3 s average duration, at 0.5 mA, and with an intertrial interval (ITI) of approximately 15 s. Mice in the CT group were exposed to the shuttle box for an identical duration of time, without receiving shocks.

In the ES phase, all mice were tested in the same shuttle box in which they were trained. Each mouse was placed into the left chamber with the door raised and was allowed to freely explore both chambers for 3 min. IR beam breaks were used as a measure of locomotor activity during this habituation phase. The door then closed automatically, and 30 trials were administered. At the beginning of each trial, the door was raised and 5 s later a foot shock (0.5 mA) was delivered. The subject’s exit from the shocked side ended the trial. If the mouse did not exit after 15 s, the shock was turned off and the trial ended. The door was lowered at the end of the trial. A session consisted of 30 trials, each separated by a 30 s ITI. Escape latencies were computed as the time from shock onset to the end of trial. If the subject failed to make a transition the maximum 15 s was used for the escape latency score. Mean latencies per subject we calculated based on performance across Trials 11-30. Statistical analysis results are listed in Table S2.

Memory Trace Tagging

Five hours before IS or CT, mice were injected with 200μl (2 mg) of 4-OHT intraperitoneally (i.p.). Following behavioral training, mice were dark housed for 3 continuous days (days 2-4). Mice were taken out of the dark on the morning of day 4, cages were changed, and they were returned to the normal colony room prior to behavioral testing. All precautions to prevent disturbances to the ArcCreER^T2 x eYFP mice during dark housing were taken to reduce off-target labeling. Following ES testing, mice were euthanized 90 min after the start of the test to allow for visualization of c-Fos expression.

Statistical Analysis

Behavioral data and summary neural data were analyzed with Prism v.9.0 or v.10.0 using an alpha value of 0.05. The behavioral effects of IS were analyzed using a t-test on averaged latencies between trials 11-30 and, using a two-way analysis of variance (ANOVA) to assess the effects of Group, Trial, and the Group x Trial interaction, where appropriate. All statistical tests and p-values for behavioral data are listed in Table S2. For comparison of regional Pearson correlation value distributions between groups, a two-sample Kolmogorov-Smirnov test was used. For comparisons of neural activation levels across targeted brain regions, multiple t-tests were conducted with a Holm-Sidak correction for multiple comparisons. All statistical tests and p-values for summary neural data are reported in Table S4.

Immunohistochemistry

Mice were deeply anesthetized and transcardially perfused as previously described (Denny et al., 2012). Brains were extracted and postfixed overnight in 4% paraformaldehyde (PFA) at 4°C, and cryoprotected in 30% sucrose / 1X phosphate buffered saline (PBS) solution at 4°C. Serial coronal sections (60 μm) were cut using a vibratome (Leica VT1000S, Leica Biosystems, Nussloch, Germany) and stored in a 0.1 M PBS /0.1% NaN3 solution at 4°C. Cells expressing eYFP and c-Fos were stained using a modified iDISCO-based protocol⁵. Briefly, sections were washed in 1X PBS 3 times for 10 min each, then dehydrated in 50% MeOH for 2.5 h. After, sections were washed in 0.2% Tween-20 in PBS (0.2% PBST) 3 times for 10 min each, and then blocked in a solution of 0.2% PBST, 10% DMSO, and 6% normal donkey serum (NDS) for 2 h. Sections were washed in 1X PBS, 0.2% Tween-20, and 10 µg heparin (PTwH), and then incubated in a primary antibody solution of rat polyclonal anti-c-Fos (1:5000, SySy, Goettingen, Germany) and chicken polyclonal anti-GFP (10mg/ml, 1:2000, Abcam, Cambridge, MA) in PTwH / 5% DMSO / 3% NDS for 48 h at 4°C. Sections were then washed 3 times in PTwH for 10 min each, before being incubated in secondary antibody solution consisting of Alexa 647 conjugated Donkey Anti-Rat IgG (1:500, Abcam, Cambridge, MA) and Cy2 conjugated Donkey Anti-Chicken IgG (1:250, Jackson ImmunoResearch, West Grove, PA) in 3% NDS/ PTwH overnight. The next day, sections were washed in 3 increments of 10 min each in PTwH, then washed in 3 increments of 10 min each in 1X PBS. Sections were mounted on slides and allowed to dry for approximately 20 min before adding mounting medium Fluoromount G (Electron Microscopy Sciences, Hatfield, PA) and a coverslip.

Confocal Microscopy

Coronal brain sections were imaged using a confocal scanning microscope (Leica TCS SP8, Leica Microsystems Inc., Wetzlar, Germany) with 2 simultaneous PMT detectors. Fluorescence from Cy2 was excited at 488 nm and detected at 500-550 nm, and Alexa Fluor 647 was excited at 634 nm and detected at 650-700 nm. Sections were imaged with a dry Leica 20x objective (NA 0.60, working distance 0.5 mm), with a pixel size of 1.08 x 1.08 µm², a z step of 3 µm, and z-stack of 9 µm. Fields of view were stitched together to form tiled images by using an automated stage and the tiling function and algorithm of the LAS X software. Raw images were saved as .lif files storing two channels: Cy2 (eYFP, channel 1) and Alexa Fluor 647 (c-Fos, channel 2).

Automated cell segmentation

Using a custom macro in ImageJ/FIJI, z-stack images of brain sections were converted from .lif to .tif files. A previously developed custom 3D segmentation approach was applied to each label due to their differences in morphology⁷. Because c-Fos localizes to the cell body, resulting in punctate filled-spherical staining, an optimized algorithm for this labelling pattern was applied. c-Fos⁺ cells were identified by first passing the image through a bandpass filter in Fourier space, subtracting the background using a rolling ball algorithm, and identifying the cells using the 3D Local Maxima Fast Filter, 3D Spot Segmentation, and 3D Manager plugins in the 3D ImageJ suite⁸. Due to the dendritic expression of eYFP after labelling with the ArcCreER^T2 x eYFP tagging system, an optimized algorithm to minimize segmentation of cellular processes was applied. eYFP⁺ cells were identified by subtracting the background, blurring the image with a Gaussian kernel, thresholding the image for the 0.5% brightest pixels, and using the thresholded regions as a mask for identifying the cells with the Classic Watershed plugin in the MorphoLibJ suite⁹. The segmentation output was then qualitatively inspected per image, and over-segmentation errors due to structural autofluorescence artifacts e.g., around ventricles and tissue edges, were manually deleted to maximize accuracy of mapping in the final dataset. Co-labeled cells were identified using the 3D MultiColoc plugin in the 3D ImageJ suite⁸, which uses the label images created during segmentation of the individual labels to efficiently identify overlapping objects. To maximize precision of counts, co-labeled cells were additionally filtered by the percent volume of overlap between the objects identified in each individual channel, relative to channel 2. A subset of identified co-labeled cells were manually inspected in ImageJ/FIJI to ensure accuracy. Custom ImageJ/FIJI macros were written to apply these segmentation and colocalization algorithms in bulk to the images. All macros are available for download at https://osf.io/ynqp7/, and were generalized to include adjustable parameters to apply to images across different resolutions.

Manual Cell Quantification

The number of eYFP⁺ and c-Fos⁺ were manually counted by two annotators for the dDG bilaterally across two images per mouse. First, the z-stack image was collapsed into a single maximum projection image using Fiji⁶ to facilitate counting. ROIs were drawn around the granule cell layer and saved for comparison with automated cell-segmentation results. eYFP⁺ and c-Fos⁺ cells were counted within each ROI, with accuracy qualitatively verified by overlaying counts over each section of the original z-stack image. To quantify the number of FN and FP counts, the total number of cells manually quantified per ROI were compared to the collapsed maximum projection image of automated cell counts within the same ROI. Despite the auto segmentation algorithms being applied in 3D, counting FN and FP cells in 2D for validation, followed by qualitative comparison with the original z-stack, allowed for a more time-efficient workflow for manual comparison. Estimated total TP counts per counter were then calculated by taking the average of two derivations: (Total no. manual counts – Total FN counts) and (Total no. auto segmented counts -Total FP counts). Individual calculations for precision, recall, and F1 were performed for each annotator and then averaged.

Package development and reference materials

The SMARTR package was developed to bridge the optimized segmentation and colocalization approach in ImageJ/FIJI with the registration and mapping capabilities of the SMART¹⁰ and WholeBrain¹¹ packages in R. Beyond this, the SMARTR package offers an extensive suite of functions to analyze and visualize the mapped datasets. The package design intrinsically incorporates data management facilitated by its object-oriented structure. Detailed documentation and in-depth descriptions of package objects and functions are available at the package website (https://mjin1812.github.io/SMARTR) and access to the source code is available as a GitHub repository (https://github.com/mjin1812/SMARTR). Additionally, the processed and mapped dataset from the LH experiments is available for download at the package website to follow along an analysis tutorial. SMARTR was written primarily in base R (v3.6.3) and developed and documented with devtools and roxygen2. Built-in functions for network analyses and visualizations were created using the igraph, tidygraph, ggplot2, and ggraph packages. The SMARTR website was built using pkgdown.

Registration & Mapping

Image z-stacks were automatically pre-processed using an ImageJ/FIJI macro to combine fluorescence intensities across both channels and z-planes into a flattened maximum project image. Images were manually assigned to a best-matching anterior-posterior atlas coordinate based on the reference atlas materials available at openbrainmap.org and https://osf.io/cpt5w. Registration with flexible user-correction was then performed in R on these flattened images using the registration() function in SMARTR, which interface with the Wholebrain(Fürth et al., 2018) and SMART(Jin et al., 2022) packages. Misaligned or damaged regions were manually annotated, logged, and excluded from analysis using the exclude_anatomy()function.

Cell counts of eYFP⁺, c-Fos⁺, and co-labeled cells were warped to atlas space using the map_cells_to_atlas()function. The area of each region per image was calculated using Gauss’s formula and aggregated across all sections in a mouse. Volume of mapped images were calculated by multiplying the aggregated areas by the thickness of each z-stack section. Mapped cell counts were then normalized by region volume to yield cells per mm³. Regions expected not to contain IEG-labeled cells, such as the fiber tracts, and ventricles, were excluded from counts in addition to manually curated lists of regions to omit due to tearing, imaging artifacts, etc. Normalized cell counts were stored and saved in the mouse object. Outlier values were detected and removed by finding normalized regional cell counts that were either two standard deviations above or below the group mean.

We improved efficiency by writing a user-friendly notebook for the registration process, and using custom scripts to detect misspellings in annotated region and automate mapping of segmented cell counts onto a standardized atlas space. We have documented and provide examples of all custom scripts in our website tutorial to facilitate use. Comprehensive statistical mapping results for eYFP⁺ cell counts, c-Fos⁺ cell counts, co-labeled / eYFP⁺, and co-labeled / c-Fos⁺ reactivation proportions with FDR corrections are listed in Table S5, Table S6, Table S7, and Table S8 respectively.

Permutation and Network Analysis

Pearson correlations between regions were calculated using rcorr() from the Hmisc package and asymptotic p-values for each correlation were determined using a one-sample t-test. Significantly correlated regions were visualized in heatmaps using conservative alpha values of p < 0.01 (eYFP⁺, co-labeled / eYFP⁺) and p < 0.005 (c-Fos⁺). Each of these thresholds were chosen to ensure an optimal level of edge sparsity in their respective networks.

Permutation analysis was conducted to identify which region correlations differed most between experimental groups. Correlations differences were calculated by subtracting correlations from the CT group from respective correlations in the IS group (rIS – rCT). Each correlation difference was used as a test statistic against a null distribution produce by shuffling the labels of CT and IS mice 1000 times, with the correlation difference recomputed for each shuffle. Correlation differences were compared to each individual null distributions to determine the p-value. Significance correlation differences for all analyses were thresholded using an alpha of 0.01.

Networks were constructed using igraph and tidygraph packages. Because the c-Fos and Arc-driven eYFP labels are both activity dependent, the correlations between regions are akin to functional connections. Thus, significant functional relationships can be visualized by representing regions as nodes and correlations as an edge between two nodes, with the Pearson correlation values determining the weight or thickness of the edge. Since pairwise correlations were calculated for all pairs, the initial functional network is complete, with all nodes fully interconnected. In order to be able to interpret the networks and discover their most salient features, correlations with alpha values below their corresponding heatmap significance thresholds were dropped and additionally filtered to have an r > 0.9 to ensure only the strongest connections were retained. Network topology measures and node centralities were then calculated using functions in the tidygraph and igraph packages. A supplementary community detection analysis was additionally applied to the c-Fos⁺ and co-labeled / eYFP⁺ channels using tidygraph functions. Close community modules were identified by calculating the eigenvector of the modularity matrix for the largest positive eigenvalue and then separating vertices into two community based on the sign of the corresponding element in the eigenvector(Newman, 2006).

Data storage and reference materials

All example data and scripts are downloadable courtesy of the OSF repository sponsored by the Center for Open Science (COS) at https://osf.io/eacn7/.

Acknowledgements

We thank Dr. Sam A. Golden for his feedback on an initial iteration of this manuscript. We also thank members of the Denny laboratory for their constructive comments on this project. This work was supported by the National Institute of General Medical Sciences (T32GM007367 to MJ/Columbia MSTP), the National Institute of Child Health and Human Development (R01HD101402 to CAD), the National Institute on Aging (R21AG064774 to CAD; F30AG084312 to MJ), and the National Institute of Neurological Disorders and Stroke (R21NS114870 to CAD). Representative schematics and figures of behavioral timelines were created with BioRender.com.

Additional information

Author Contributions

MJ, ML, SLS, and CAD conceived of the pipeline and contributed to intellectual development of analysis approaches. MJ programmed, packaged, and created SMARTR and its supporting webpage, with contributions from SO. MJ and AW ran behavioral experiments. MJ and SS performed manual cell count verifications and AF assisted with validation analysis. MJ and SO performed immunohistochemistry and imaging. MJ and SO analyzed the dataset and created visualizations. MJ and CAD wrote the final manuscript with input from all authors.

Ethics

All experiments were approved by the Institutional Animal Care and Use Committee (IACUC) at the New York State Psychiatric Institute (NYSPI). All key resources are listed in Table S1.

Supplemental figures

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Figure S7.

Figure S8.

Figure S9.

Figure S10.

Figure S11.

Figure S12.

Figure S13.