Longitudinal monitoring of prefrontal cortical ensemble dynamics reveals new insights into stress habituation.

The prefrontal cortex is highly susceptible to the detrimental effects of stress and has been implicated in the pathogenesis of stress-related psychiatric disorders. It is not well understood, however, how stress is represented at the neuronal level in the prefrontal cortical neuronal ensembles. Even less understood is how the representation of stress changes over time with repeated exposure. Here we show that the prelimbic prefrontal neuronal ensemble representation of foot shock stress exhibits rapid spatial drift within and between sessions. Despite this rapid spatial drift of the ensemble, the representation of the stressor itself stabilizes over days. Our results suggest that stress is represented by rapidly drifting ensembles and despite this rapid drift, important features of the neuronal representation are stabilized, suggesting a neural correlate of stress habituation is present within prefrontal cortical neuron populations.

Introduction

Up to 90% of people will experience at least one traumatic experience during their lifetime (Kessler et al., 1995; Kilpatrick et al., 2013), and traumatic experiences are an independent risk-factor for most psychiatric disorders (Rundell et al., 1989; Felitti et al., 1998). Severe and repeated stress exposure can lead to disruptions in emotional processing, executive function, and even general health outcomes (Felitti et al., 1998, Mcewen, 2012). Thus, understanding the effects of stress on brain structure and function and how stressful experiences affect these fundamental processes could reveal new insight into the pathophysiology of stress-related psychiatric disorders.

Stimuli, actions, and memory are thought to be encoded by groups of neurons that fire in a coordinated manner called ensembles. Some groups have suggested that once memories are formed, the ensembles are stabilized and the information from the same group of neurons is extracted during memory recall (Tonegawa et al., 2015; Katlowitz et al., 2018; Hebb, 1949). However, the stability of these ensembles has been called into question with representations at both the sensory and mnemonic level drifting even at the relatively short timescale of days (Lutcke et al., 2013; Rule et al., 2019; Gonzalez et al., 2019; Driscoll et al., 2017). While efforts have been made to understand how sensory stimuli, motor actions, and memories are encoded at the neural level, much less is known about how stress exposure is represented. Even less is known about how stable these ensembles are and how habituation to stress is represented within these ensembles.

In rodents, the prelimbic prefrontal cortex (PL) is interconnected with stress reactive regions such as the basolateral amygdala, hippocampus, and infralimbic prefrontal cortex and is highly stress reactive (Arnsten, 2015; Jacobs and Moghaddam, 2021). Moreover, the PL receives inputs that carry information about the internal and external environment such as the amygdala, insular cortex, and midline thalamic nuclei (Hoover and Vertes, 2007). Thus, the PL is ideally located to coordinate distinct behavioral, physiological and emotional responses to stress exposure (Tovote et al., 2015; Apps and Strata, 2015) by balancing aversive and reward responses with populations of neurons that represent both stressful and rewarding stimuli (Del Arco et al., 2020). Single unit recordings of PL neurons during restraint or tail pinch stress have shown that neurons show a mix of activity changes during the stressor which include subsets of neurons exhibiting phasic increases in activity, sustained activity during the stressor, increased activity at the beginning and end of the stressor, and increase in local field potential theta frequency (Del Arco et al., 2020; Jackson and Moghaddam, 2006).

Foot-shock stress is brief and can be delivered in a repeated fashion which allows for the study of real-time stress effects and their evolution through the duration of the stressor. Historically, one thing that has hindered efforts to record neural activity during foot-shock stress has been generation of electrical artifacts which obscure acquisition of neural activity during the foot-shock itself. In the present study we use in vivo single neuron calcium imaging to track PL neural activity patterns and temporal dynamics within and between multiple unpredictable foot-shock exposures. The key results from our work indicate that while large numbers of PL neurons exhibit dynamic changes in activity during a session of unpredictable foot-shock exposures, cells active during any specific shock is comprised of spatially distinct rapidly drifting neuronal ensembles. We also show, despite this persistent representational drift over multiple days, the PL neural representation of shock experience stabilizes, which could reveal new insights into the mechanisms subserving stress habituation.

Materials and methods

Subjects

All experiments were approved by the Vanderbilt University Institutional Animal Care and Use Committees and were conducted in accordance with the National Institute of Health guidelines for the Care and Use of Laboratory Animals. A total of 19 male C57BL/6J mice older than 8 weeks and obtained from Jackson Labs or bred in our animal facility were used for used in our experiments. Mice were housed in a temperature and humidity-controlled housing facility under a 12h light/dark cycle with ad libitum access to food.

Surgery

Mice were initially anesthetized with 5% isoflurane and then transferred to the stereotax (Kopf Instruments, Tujunga, CA) and kept under 2% isoflurane anesthesia. The hair over the incision cite was trimmed and the skin was prepped with alcohol and iodine scrub. The skull was exposed via a midline incision and treated with the local anesthetic, benzocaine (Medline Industries, Brentwood, TN). For all surgeries, we used a motorized digital software (NeuroStar; Stoelting Co., Wood Dale, IL) to guide a 10 mL microinjection syringe (Hamilton Co., Reno, NV) driven by a Micropump Controller (World Precision Instruments, Sarasota, FL). AAVs encoding GCaMP7f under control of the synapsin promoter (AAVrg-syn-jGCaMP7f-WPRE with titer of 1.8 × 10^13 vg/mL; Addgene, Watertown, MA) were injected into the prelimbic prefrontal cortex at 2 dorsoventral levels (AP:+2.42, ML: + 0.35, DV: 1.6 and 1.9), 300 nL per level. After allowing for diffusion of the virus, the needle was withdrawn, and a 0.5 mm diameter needle attached to a vacuum line was lowered stereotaxically to a depth of 1.55 over PL to create a tract. Thereafter, a 0.5 mm diameter GRIN lens (Inscopix, Palo Alto, CA) was slowly lowered stereotaxically into PL at DV 1.8 mm. The GRIN lens was affixed to the exposed skull using Metabond (Parkel, Edgewood, NY). Following completion of the surgery, 10 mg/kg ketoprofen (AlliVet, St. Hialeah, FL) was administered as an analgesic, and post-operative treatment with ketoprofen being administered 24 and 48 h after the surgery. Animals were allowed to recover for at least 2 weeks, after which a baseplate was installed over the lens in order to dock the miniaturized microscope at an empirically optimized working distance.

Behavior and recordings

Mice were habituated to the miniaturized epifluorescence microscope (Inscopix, Palo Alto, CA) for at least 2 days. On test days, mice were attached to the microscope and placed in fear conditioning chambers (Coulbourn Instruments, Allentown, PA) contained in sound attenuating cabinets (Coulbourn, Allentown, PA). Mice were exposed to 20 foot-shocks at 0.5 mA intensity and 2 s duration delivered at pseudorandom order using FreezeFrame software (Actimetrics, Wilmette, IL). Freezing data was obtained from FreezeFrame. Motion data was generated using motion index data obtained from FreezeFrame. The motion index is an internal measure in FreezeFrame used by the software to generate freezing data. The calcium imaging data acquisition system was interfaced with the experimental control software using custom made telephone connector to BNC connector cables to deliver TTL pulses signaling the start of the behavioral session and shock delivery onset to the data acquisition system. Data was acquired at a frame rate of 10 Hz. LED power, gain, and lens focus were empirically adjusted to maximize the quality of the recordings while minimizing the LED power used. Data was acquired continuously throughout the duration of the session. A subset of 3 mice were used for longitudinal registration experiments in which the activity of neurons was tracked across 3 behavioral sessions occurring across 3 days.

For cFOS immunohistochemistry experiments, 6 mice were exposed to 20 foot-shocks as described above and another 6 mice were placed in the behavioral apparatus for the same amount of time but were not exposed to shocks. Mice were perfused for immunohistochemistry 60 min following the end of the session.

Histology

To confirm location of GCaMP7f expression and lens location, mice were anesthetized using isoflurane and transcardially perfused with ice-cold phosphate buffered saline (PBS) followed by 4% paraformaldehyde (PFA) solution in PBS. Brains were extracted and stored overnight in 4% PFA and transferred to a 30% sucrose solution for 2 days for cryoprotection. 40 μm brain sections were taken using a Leica CM3050 S cryostat (Leica Microsystem, Weitzlar, Germany). Slices were then directly mounted on glass slides and VectaShield H-1200 DAPI mounting medium (Vector Laboratories, Burlingame, CA) was applied before coverslipping. Images were taken using an Axio Imager M2 epifluorescent microscope coupled to a fluorescent lamp.

Immunohistochemistry

Mice were processed as described above. Brain sections were washed in Tris-Buffered Saline (TBS) 3X for 10 min. Slices were then incubated in a blocking buffer TBS containing 0.2% Triton X-100 (Fisher Scientific) and 4% horse serum (TBS+) for 90 min. Slices were then incubated overnight in a solution containing the blocking buffer and 1:500 rabbit-anti-cFOS antibody (Abcam, Cambridge, UK) plus anti-NeuN (1: 200, conjugated with rabbit polyclonal Alexa Fluor 488; EMD Millipore). The following day, slices were washed 3X in TBS+, and then incubated for 2 h in a solution containing the blocking buffer and 1:1000 Alexa Fluor 546 Donkey-anti-Rabbit (Invitrogen, Waltham, MA). Slices were subsequently washed and mounted for image analysis.

Images of PL were obtained as described above. Raw 10x images for the green and red channels of the images were then opened with ImageJ. A threshold was set for the images, PL boundaries were drawn, and area of the region was measured, and ImageJ particle counting software was then used to analyze the number of cFOS+ nuclei. Raw counts were then normalized by dividing by the measured region area. Data were analyzed using an unpaired Student's t-test.

Data analysis

Data analysis of Calcium imaging recordings was performed as previously described (Marcus et al., 2020). Briefly, recordings were spatially downsampled by a factor of 2, bandpass filtered, and motion corrected using Inscopix Data Processing software V1.3. Individual Ca2+ traces were extracted using Constrained Nonnegative Matrix Factorization for microendoscopic data (CNMFe) (Zhou et al., 2018). For CNMFe, we used the following parameters: min corr = 0.9, min pnr = 20, gSiz = 20, gSig = 10. Individual extracted traces and corresponding identified neurons were visually confirmed, and false positive traces were excluded from further analysis. Analysis was conducted using custom made MATLAB codes. Individual calcium traces were aligned to the shock onset periods and averaged across 2 consecutive shock periods or 20 shock periods for within session analyses or for the whole session respectively. The resulting average trace was binned into 1 s intervals and the data was Z normalized using the 10 s pre-shock period as baseline. Traces exceeding a Z-score value of 3 or -3 (p = 0.01, two-tailed) for any 2 consecutive bins following shock onset were considered excitatory shock-responsive or inhibitory shock-responsive, respectively. Longitudinal registration across 3 days was conducted using Inscopix Data Processing software using 0.5 minimum correlation as threshold. The resulting aligned traces were individually reviewed for continuity across the 3 recording days and analyzed as described above. For statistical analyses we used Repeated Measures ANOVA followed by Šídák's multiple comparisons post hoc test, two-tailed Student's t-test, and Fisher's exact test as appropriate. Results for post-hoc comparisons are shown in Supplementary Table 1.

Agglomerative hierarchical clustering was conducted using available MATLAB functions. The post-shock Z-scored traces corresponding to the analyzed periods were concatenated together with each row corresponding to the traces from a single neuron. Linkage was then calculated using Ward's method (inner squared distance) with Euclidean distance for metric. Then pairs of neurons in close linkage proximity were clustered and then these clusters were clustered with other clusters in close proximity, this process was repeated iteratively until a hierarchical tree (a dendrogram) was built. The clustergram function was used to generate the dendrogram. 3% of the maximum linkage was used as a threshold to determine the different clusters which are identified by distinct colors within the dendrogram and its associated heatmap. For visualization, these clusters and the corresponding portion of the heatmap was manually separated using Adobe Illustrator. The average and standard error of the mean for each cluster's time course was plotted in GraphPad Prism V9.1.

Principal component analysis (PCA) was conducted using GraphPad Prism V9.1. The post-shock Z-scored traces corresponding to the analyzed periods were labeled with their corresponding epoch and analyzed together. The obtained PCA scores for the first 3 components for each of the analyzed epochs were plotted. Centroid distances in 3d space between the different epochs graphed were calculated and graphed. More specifically, centroid distances were calculated between the early, mid, and late-shock epochs or between early and late epochs across the 3 shock sessions.

Our graph theory analysis was based on the methods and fragments of MATLAB codes used by Jimenez et al. (2020). First, we created a binary matrix with the corresponding calcium spike onsets for each individual calcium trace; time resolution was 10 frames per second. The binary traces were aligned around shock onsets averaged into blocks of 2. The resulting time locked binary traces were binned into 1 s bins and the pre-shock period was discarded. Pearson correlation coefficients were calculated between each of the binary traces to create and adjacency matrix. Strongly correlated pairs were selected by selecting pairs with Pearson correlations >0. 5. Graph objects were generated in MATLAB from this adjacency matrix. The number of correlated pairs for each recording was quantified and the number was normalized by dividing by the total number of neurons for that field of view. Neurons that were previously identified as (+) shock-responsive were identified in the adjacency matrix and the number of correlated pairs between the binary traces for these neurons and all other neurons or other shock-responsive neurons were quantified and normalized as described above.

Decoder analysis for each session was conducted using available MATLAB functions. We used Z-scored post-foot-shock neuronal data for the first 2 individual trials and last 2 individual trials. Data sets were split into training and testing data sets by randomly selecting individual traces; the testing data set consisted of 1/3 of the traces from the full data set. The number of trials from early and late epochs were equally represented in both data sets. PCA was conducted across the temporal dimension; we kept 10 out of 11 total dimensions. Afterwards, row vectors were created by concatenating randomly picked post-PCA traces. We created 150 row vectors from each epoch. We applied ICA to the resulting matrices to obtain 100 features from each row vector. Afterwards, we trained a fine gaussian support vector machine (SVM) decoder by feeding it the resultant training data matrix together with trial identities. Afterwards, 150 row vectors were randomly picked from the testing data set and fed into the trained decoder and accuracy was calculated; this was repeated for 100 iterations to obtain an accurate depiction of the decoder accuracy. Accuracy was calculated by obtaining the proportion of correct predictions to total prediction attempts. This whole process was performed 28 times to assess the variability of the procedure. As a control, we trained the decoder with data in which the trial identities had been shuffled and we obtained the empirical chance decoder accuracy.

Results

Distinct foot-shock-induced cellular response patterns in the PL

The PL is a stress reactive region, however the temporal dynamics and patterns of activity exhibited by PL neurons remains poorly understood. To examine PL stress reactivity, we first exposed mice to repeated unpredictable foot-shocks and subsequently quantified immunoreactive nuclei to the activity marker cFos. As shown in Fig. 1A, foot-shock stress increased the number of cFos-positive cells relative to non-shocked control mice (t(10) = 2.444, p = 0.0346), thus, confirming the stress reactivity of the PL. To gain insight into the temporal dynamics of PL neuron activity in response to aversive stimuli we used in vivo single-cell calcium imaging and recorded neuronal responses to repeated unpredictable foot-shock in male mice. Fig. 1B depicts GRIN lens placement and example field of view of neurons detected using the CNMFe approach (Zhou et al., 2018). Fig. 1C summarizes the experimental trajectory from lens implantation to recordings. A total of 752 neurons from 7 mice were used for analysis. We first averaged calcium responses over all 20 foot-shocks (10 s pre-shock baseline and 10-s post-shock) for each neuron and Z-scored data to the baseline activity of each neuron. Sorting the Z-scored activity for all neurons in order of average post-shock magnitude shows shock exposure bidirectionally modulated neural activity (Fig. 1D). By averaging all Z-scored.

Fig. 1

Foot-shock stress modulates prelimbic neuronal ensembles. (A) Left. Coronal diagram depicting the location of the prelimbic cortex and images of cFOS + neurons at 10x magnification from foot-shock unexposed and exposed mice. Right. Foot-shock exposure increases the number of cFos + nuclei/per area in PL (n = 6 per group; two-tailed unpaired Student's t-test; *p < 0.05). (B) Left. Coronal diagram of GRIN lens positioned above GCaMP7f expression in PL. Center. Maximal projection image of Ca2+ transients and extracted cell contours from representative mouse. Right. Photograph of GCaMP7f expression in coronal section of PL underneath tract left by GRIN lens. Blow out photograph depicts the approximate focal plane of the recordings. (C) Experimental timeline. (D) Average neuronal calcium activity for 20 foot-shocks Z-scored and sorted in ascending order for all recorded individual neurons and (E) Average population response for all recorded neurons showing increase in average neuronal activity following foot-shocks. The neuronal population shock-response is significantly larger than that seen for that calculated for temporally permuted traces. (n = 752 neurons; N = 7 mice; Two Way ANOVA; P = 0.0006). Post-hoc comparison results can be found in supplementary Table 1 (F) Proportion of neurons with no change and significant increases ((+) Responsive) and decreases ((−) Responsive) in activity following foot-shock. We detected a few neurons that crossed the (+) and (−) significance threshold in response to the foot-shock (biphasic). Inset, proportion of response types for temporally permuted traces. Fisher's exact test for comparisons between shock-response types and the equivalent for temporally permuted neurons (****P < 0.0001). (G) Left. Average population response for (+) Responsive neurons. Right. Example of calcium transients for neurons increasing their activity in response to the foot-shock. (H) Left. Average population response for (−) Responsive neurons. Right. Example of calcium transients for neurons decreasing their activity in response to the foot-shock. For in vivo recordings, n = 752 neurons from 7 mice.

Traces, we found a net increase in overall PL activity following foot-shock exposure relative to pre-shock baseline (Fig. 1E and S1). Importantly, this average increase in activity to foot-shocks is significantly larger than that seen for the permuted-trace control (F(1, 31542) = 11.11, p = 0.0009, Fig. 1E). We next grouped cells based on their activity pattern in response to shock exposure. Neurons with post-shock responses exceeding ±3 Z-scores (p = 0.01, two-tails) in any 2 consecutive 1s bins during the 10s post-shock period were considered to have significant excitatory or inhibitory responses ((+) responsive and (−) responsive), respectively. We chose this statistical threshold as it has been successfully used by other groups studying in vivo neuronal activity in the PL (Siciliano et al., 2019; Burgos-Robles et al., 2009). ∼28% of neurons were (+) responsive (Fig. 1F and G) while ∼22% were (−) responsive (Fig. 1F and H). A small group of neurons (∼2%) responded to the foot-shock in a biphasic manner where they crossed both the (+) and (−) significance threshold during the post-shock period. ∼49% of neurons were non-responsive. These proportions were significantly different than those obtained from permuted traces analyzed in a similar manner (Fisher's exact test, Fig. 1F inset).

We next used agglomerative hierarchical clustering to identify distinct activity patterns within our neuronal population. Analysis using Euclidean distances and Ward linkages revealed 7 clusters. Visual representation of the cluster analysis and extracted time course data for all 7 clusters is shown in Fig. S2A. These data reveal several distinct patterns of activity embedded within the PL response to foot-shock with most clusters changing their activity to shock mostly differing in the magnitude of the response and timing (delayed and sustained vs. transient, for example). A similar analysis of permuted traces using the same parameters yielded 5 clusters with the magnitude of changes within each cluster being substantially smaller than for shock-aligned traces (Fig. S2B). Overall, these data indicate that PL neurons exhibit multifaceted and distinct activity patterns in response to acute aversive stimuli.

Within session analysis of PL neuron responses to foot-shock exposure

While analysis of averaged responses to foot-shock revealed several distinct response patterns within PL neurons, these data do not reveal whether PL neuronal responses show any within session adaptations. To address this question, we first examined cellular activity in response to shocks throughout the session by averaging the activity across 2 sequential foot-shocks for each neuron for 10 trial blocks in total (derived from 20 foot-shocks). This analysis revealed that most neurons respond significantly (±3 Z-scores) to the shock at some point during the session with ∼90% of neurons significantly increasing, and ∼63% significantly decreasing, their activity during a least one of the post-shock periods analyzed (Fig. 2A). The proportions of (+) responsive cells using this analysis (∼90%) is substantially higher than for the averaged session data presented in Fig. 1 (∼28%), suggesting that some neurons are highly responsive to shock presentation at one or a few points during the session, but not active for most others. In this situation, the average responses may not meet the 3 Z-score threshold, resulting a lower percentage of (+) responsive cells in the averaged data.

Fig. 2

Within session changes in neuronal activity patterns in response to foot-shock. (A) Top. Proportion of neurons with (+) responses to any of the 10 shock blocks analyzed. Bottom. Proportion of neurons with (−) responses to any of the 10 shock blocks analyzed. (B) Top. Representation of 20 shocks delivered through the session and analyzed blocks of 2 (colored) for the early, middle, and late epochs of the session. Bottom. The magnitude of the population response to the foot-shock increases by mid-session. (C) Top. The magnitude of the population shock-response increases by the middle of the session (mid-shock) for (+) Responsive neurons and returns to a similar response magnitude to that of early-shock by the end of the session (late-shock). Inset. Proportion of neurons with (+) responses during the 3 different foot-shock epochs. Bottom. Venn diagram showing minimal overlap between (+) shock-responsive neuronal populations at the 3 different epochs of the session. (D) Top. The magnitude of the population shock-response for (−) Responsive neurons revealed a significant time × epoch interaction (P = 0.0002). Post-hoc analysis did not reveal any differences between epoch differences. Inset. Proportion of neurons with (−) responses during the 3 foot-shock epochs. Bottom. Venn diagram showing minimal overlap between (−) shock-responsive neuronal populations at the 3 different epochs of the session. Repeated measures ANOVA where appropriate followed by Šídák's multiple comparisons post hoc test, *P < 0.05. Post-hoc comparison results can be found in supplementary table. (E) Sorting neuronal activity changes in response to foot-shock for the epoch of interest in ascending order reveals lack of stable response pattern during the other 2 sessions. (F) Plots showing the PC score trajectories for the 3 epochs for the first 3 PCs. (G) The distance for 3d PC score centroids for the 3 epochs is similar between each other (left) and from the origin (right). For in vivo recordings, n = 752 neurons from 7 mice.

To explore this hypothesis, we compared the activity patterns of neurons between early, middle, and late shock exposures (Fig. 2B). Analysis of the average population activity to foot-shock stress revealed a significant effect of epoch (F(2,47313) = 7.398, p = 0.0006, Fig. 2B) with middle and late shocks showing higher overall Z-score changes relative to early shocks. This suggests that overall PL responsivity to shock increases through the session. This increase in average activity of PL neurons could be due to an increase in the number of neurons responding to the shock during middle and late exposures, or an overall enhancement in the magnitude of shock responses by the end of the session, or both. To address this question, we analyzed significant changes in activity of single neurons by focusing on cells that exceed the 3 Z-score threshold during the post-shock period during early, middle, and late time epochs of the foot-shock stress session. This analysis revealed a significant difference in the magnitude of neuronal responses between epochs and a significant time × epoch interaction for (+) responsive neurons (F(2, 574) = 3.449, p = 0.0324 and F(40, 11480) = 2.233, p < 0.0001, respectively, Fig. 2C). To assess the contribution of the number of responsive neurons, we analyzed the proportion of neurons responding shocks at the different time epochs. The Venn diagram for this population of neurons demonstrates that more neurons respond positively to late shocks relative to both early and mid-session shocks; a relatively small number of neurons respond positively during at least 2 epochs. Importantly, only 5% of cells were active during all epochs (Fig. 2C).

Similar analysis of neurons that decrease their activity under 3 Z-scores revealed a significant time × epoch interaction (F(40, 5820) = 2.017, p = 0.0002, Fig. 2D). As shown in the Venn diagram, more neurons respond with inhibition to both mid and late shocks relative early shocks; a relatively small number of neurons respond in an inhibitory fashion during at least 2 epochs and, again, only 2% of cells were inhibited during all epochs (Fig. 2D). These data indicate that the ensemble of neurons activate during the early, middle and late shock exposures are comprised of largely distinct groups.

To further confirm the lack of stability of the neural representation of foot-shocks across epochs, we sorted the Z-scored data for early shock presentation in ascending order (left to right) of shock-evoked activity and visualized the activity of these sorted cells across all middle and late epochs. Similarly, sorting by the activity of middle or late shock periods revealed a similar lack of stability across time (Fig. 2E). These data support our hypothesis that PL neurons show robust changes to only a subset of foot-shock exposures, such that by the end of the session ∼90% of neurons had demonstrated some significant increase in activity to shock exposure.

These results suggest that mid to late increases in PL net neuronal activity are mediated by a combination of increases in the proportion of neurons responding to late shocks and to a lesser degree the magnitude of the responses. To further examine whether the resulting shock responses across the 3 distinct epochs are different from each other we applied principal component analysis to the post-shock period of the Z-scored traces (Fig. 2F). Graphical representation of the three epochs in principal component space reveal that the Euclidean distance between the centroids of these groups is equidistant from each other and from the origin (Fig. 2G). Thus, the ensembles representing the 3 different time epochs are mostly distinct from each other in the overall properties of their activity.

Between session analysis reveals neuronal correlates of stress habituation

It has been shown that repeated stressors can lead to habituation (Herman, 2013; Grissom and Bhatnagar, 2009; Patel et al., 2005; Campmany et al., 1996). It is not known, however, how habituation is reflected in PL neuronal dynamic responses to foot-shock stress. To address this question, we recorded the activity of PL neurons longitudinally across 3 days in a subset of mice. As shown in Fig. 3A, foot-shock stress increases the average magnitude of neuronal activity by the end of the Day 1 session as shown by a significant group difference between early and late epoch and a time × epoch interaction (F(1,686) = 5.686, p = 0.0174 and F(20,13720) = 2.95, p < 0.0001, respectively). By Day 3, however, PL neurons fail to exhibit a net increase in average activity by the end of the session, demonstrating a between session habituation of the escalating response over time observed only on Day 1. This is further supported by a lack of a significant difference when directly comparing Day 1 early and Day 3 late epochs with repeated measures ANOVA (F (1, 686) = 1.037, P = 0.3089).As shown in Fig. S3, we saw a similar pattern of habituation behaviorally as suggested by distinct changes in post shock motion index on day 1 between early and late shocks but similar post-shock motion index between early and late shocks by Day 3. We did not observe significant shock-induced freezing responses (Fig. S3).

Fig. 3

Neuronal response during habituation to foot-shocks across 3 days. (A) Average population response for longitudinally registered neurons across 3 days. The magnitude of the neuronal population response to foot-shocks increases by the end of the session on Day 1, but not on Days 2 and 3. Repeated measures ANOVA followed by Šídák's multiple comparisons post hoc test. Post-hoc comparison results can be found in supplementary table. (B) Pie charts. Repeated foot-shocks across days increases the proportion of neurons that respond to the foot-shock early in the session and decreases the proportion of neurons that respond to foot-shock by the end of the session. Fisher's exact test. Right. The ratio of the proportion neurons responding to the foot-shock early in the session to the that of neurons responding to the foot-shock late in the session increases across 3 days. (C) Heatmaps of post shock neuronal activity (Z-score) following early and late epoch foot-shocks across 3 days were generated and sorted according to the activity of each of the epochs. Sorting of the neuronal activity following foot-shocks did not show any patterns in the neuronal responses for any of the other epochs. *p < 0.05 and **p < 0.01. n = 344 neurons from 3 mice.

The lack of increase in the average PL neuronal activity after Day 2, may be the result of an overall decrease in the magnitude of neuronal responses to the shock, an increase in the magnitude of cells decreasing activity to the shock, or a redistribution of the proportion of cells that increase and decrease their activity to the shock. Thus, we analyzed shock-induced changes in activity of single neurons by focusing on neurons that exceed an absolute value of 3 Z-scores in response to the shock during early and late epochs across the 3 days of foot-shock stress. As shown in Fig. S4A, analysis of shock-responsive neurons did not reveal any differences in their net activity between early and late epochs for Days 1 and 2. For Day 3, however, the magnitude of the shock response was higher for the late epoch as indicated by a main effect in epoch (F(1, 154) = 5.085, p = 0.0255).

Given that we did not see differences in the magnitude of the net response by Day 3 and there is only a slight increase in the magnitude of the shock response for (+) responsive neurons between early and late epochs on Day 3, the loss of increase in net activity for the late epochs for Days 2 and 3 should be due to redistribution of the populations of cells responding to the shock. To test this, we classified the proportions of cells responding to shocks in different manners to early and late shocks across the 3 days. More specifically, cells responding in either (+) or (−) manner to early, late, or both shocks. For simplicity, we excluded a small proportion of cells that responded in mutually opposite manner to early and late shocks and cells that responded in both (+) and (−) manner in the same post-shock period. As shown in Fig. 3B, the proportion of cells responding to early shocks increased through the days relative to Day 1, while the proportion of cells responding to late shocks decreased by Day 3 (pairwise Fisher's exact test). Moreover, analysis of the ratio of (+) responsive neurons responding to late shocks relative to early shocks revealed an increase across days which suggests that the number of (+) responsive neurons between early and late epochs becomes more similar across days (pairwise chi square test, Fig. 3B). It is important to note that the overall proportion of neurons responding in any manner to foot-shock stress did not differ significantly across days (Fig. S4B). Thus, habituation to repeated foot-shock stress is associated with a reduction in the differential proportion of activated neurons between early and late shock responses in PL, which could in turn blunt within session increases in net PL activity on Days 2 and 3, relative to Day 1. Furthermore, this reduction in the differential response is due to a redistribution of the responsive population of neurons as there was no difference across days in the overall proportion of neurons with any type of response.

Our data thus far indicate that there is considerable and rapid within session ensemble drift, however the degree to which this spatial drift is present between sessions is unclear. It is possible the increasing similarity between the representation of early and late shock exposures described above could be due to stabilization of the neuronal populations responding to shock exposure, i.e. reduced representational drift. To test this hypothesis, we aligned the activity for early, middle, and late epochs for Days 1–3 and sorted the activity by each of the epochs. As shown in Fig. 3C, sorting by activity of any of the epochs did not reveal any pattern in the activity of the other epochs in any of the 3 Days. Thus, responses in any of the epochs do not predict responses in any other epoch and this lack of pattern holds across the 3 days of foot-shock stress suggesting that the representation of the shock response drifts within and between days. These data do not support the hypothesis that the increasing similarity in PL ensemble response patterns between early and late shock across days is mediated via reduced representational drift.

Stabilization of foot-shock representation across days

To further investigate how increasing neuronal similarities within sessions occurs in the context of persistent representational drift across days, we utilized several approaches to evaluate changes in activity between early and late shock exposures across days. We first applied principal component analysis to the post-shock period of the Z-scored traces of the early and late epochs across the 3 days of stress. The data for the first 3 principal components, which account for 38% of the variance within the data cumulatively, was plotted in 2D component space (Fig. 4A). The Euclidean distance between the centroids for early and late epochs for each foot-shock session was calculated. This analysis revealed that the distance between the centroids for Day 3 dropped relatively to Days 1 and 2 (Fig. 4A). This suggests that the patterns of post-shock neuronal activity between early and late epochs become more similar by Day 3 relative to the other days. These findings mirror our earlier results indicating more similar population neuronal dynamics between early and late shocks on day 3 relative to Day 1 (see Fig. 3).

Fig. 4

Neuronal response habituation to foot-shocks across 3 days in PL. (A) Left and center. Plots showing the PC score trajectories for the early and late epochs across 3 days for the first 3 PCs in response to foot-shocks. The distance between the PC score trajectories for early and late foot-shock epochs is reduced by Day 3. Right. The distance between the centroids for early and late epochs is dramatically reduced by Day 3. (B) Graphs obtained from graph theory analysis showing the correlations between all neurons (top) and between (+) Responsive neurons (bottom) in the field of view in response to shock on Day 1 (left) and Day 3 (right). (C) Quantification of correlations following foot-shock from graph theory analysis between any of the neurons in the FOV (left), between any neurons correlated with (+) Responsive neurons in the FOV, and only between (+) Responsive neurons in the FOV. The increase in correlated pairs in response to foot-shock is blunted after the first day of foot-shock stress. Two way ANOVA followed by Šídák's multiple comparisons post hoc test comparing early and late epochs for each day. *p < 0.05 and **p < 0.01. n = 344 neurons from 3 mice. (D) SVM decoder workflow (left). Depiction of N-dimensional hyperplane used by decoder to classify data into Early and Late shock categories. (E) Decoder accuracy when attempting to predict if PL data corresponds to Early Shock or Late Shock trials. Decoder is able to decode Early Shock from Late Shock trials on Day 1, but not for subsequent days. n = 28 iterations; two-tailed unpaired Student's t-test; ****p < 0.0001.

To examine in more detail how neuronal activity in the PL changes across repeated days of stress exposure, we used a graph theory approach (Jimenez et al., 2020; Bullmore and Sporns, 2009) to determine how neuron-neuron firing dynamics change across repeated days of foot-shock exposure. In this analysis, neurons with significant correlations in their spike activity during the post-shock period were considered to form correlated pairs. Analysis of mutually correlated neurons revealed a significant epoch (early vs. late) x day interaction (F(2, 4) = 15.99, p = 0.0124, Fig. 4B and C). Post hoc comparisons showed that on Day 1 of foot-shock stress, the number of mutually correlated neurons increased between the early and late epochs of the session (p = 0.0367), however, by the 3rd day of shock exposure there was a reduction between the early and late epochs (Fig. 4C). Further analysis focusing on the correlated activity of neurons that have been classified as shock-responsive for the corresponding epoch revealed similar increases in the number of correlated neuronal pairs between early and late epochs in Day 1, but not Days 2 and 3, between shock-responsive neurons with any neuron type, and between shock-responsive neurons with only other shock-responsive neurons (epoch × day interaction: F(2, 4) = 7.769, p = 0.0419 and F(2, 4) = 17.01, p = 0.0111, respectively; Fig. 4B and C). These results suggest that repeated sessions of foot-shock stress are associate with a blunting in the progressive increase in correlated activity that occurs within session on Day 1, which is most striking when examining correlations between only shock responsive neurons. The failure to increase correlated activity within sessions across days may represent a neural correlate of stress habituation within PL neuronal ensembles.

Our results suggest that PL neural dynamics on Day 1 are distinct between early and late epochs of stress exposure, and they become more similar on subsequent days. Thus, we hypothesized that PL neuronal activity could be used to predict the stress epoch on Day 1 but not on subsequent days. Indeed, using an SVM decoder trained for each day independently (Fig. 4D) revealed that the identity of stress epochs is decodable from each other on Day 1 (t(54) = 7.515, p < 0.0001), while on subsequent days, the decoders performed at chance levels (Fig. 4E). This further suggest that the PL representation of shock stabilizes over days and provides additional evidence of habituation at the cellular level.

Discussion

In the present study using unpredictable foot-shock stress we expand on previous findings that show that PL neurons are stress responsive (Marcus et al., 2020; Del Arco et al., 2020; Jackson and Moghaddam, 2006). The use of calcium imaging across time allowed us to track individual neurons within and between foot-shock exposure sessions. Using this approach, we show that the PL representation of foot-shock exposure lacks stability across time; most neurons (∼90%) show positive responses to foot-shock at some point within the first session, however we found minimal overlap in the population of neurons responding during early, middle, and late shock exposures within the same session. These data reveal a rapid and persistent representational drift within the PL ensemble response to shock exposure. We also found evidence for between session adaptation in the neuronal representation of the shock exposure. Specifically, the number of neurons responding to the shock increased across the first session. However, repeated sessions were associated with an attenuated within-session activity increase, as reflected in a smaller proportion of late shock responsive neurons on day 3 relative to day 1. Consistent with this adaptation, our graph theory analysis revealed that the increase in correlated activity across the day 1 session was absent on days 2 and 3. Lastly, we also found that the representation of the shock at the PL neuron population level appears to stabilize across days as evidenced by closer distance between early and late activity in principal component (PC) space on day 3 relative to day 1. This stabilization occurs despite persistent representational drift. Together these finding provide new insight into how PL neurons represent aversive experience across time and could provide new insights into the mechanisms of stress adaptation and habituation.

Our results are in line with previous studies that show that stress increases activity of prelimbic neurons during acute stressors as shown by increased expression of the activity marker cFos and in vivo single neuron activity. While cFos only reveals the subgroup of neurons that increase activity at some point during the stress session, our recording of in vivo neuronal activity shows that neurons increase, decrease, or modulate their activity in a biphasic manner in response to the shock itself. These changes in activity cannot be attributed to chance, as applying similar analyses to permuted traces revealed only a small proportion of responsive neurons using the same statistical thresholds.

Our data shows that most PL neurons exhibit significant changes in activity during at least one epoch contrasts to our averaged data shown, where <50% of neurons showed significant average changes in activity across all 20 trials. These data suggest subgroups of neurons within the PL are responsive to only a sub-set of shock exposures. Indeed, the degree of overlap in neuronal activity between early, middle, and late epochs was relatively low for both shock-excited and shock-inhibited neurons. These data indicate that the ensemble of activated neurons in response to shock-exposure migrates spatially across repeated stimulus presentations. These data might suggest the experience of the stress may not be represented by a static population but rather a migrating population that represents similar information. This hypothesis is partially supported by our PCA which revealed that the centroid distance within PC space between the representations of the shock is similar distance from the origin and between each of the different epochs (early, middle, and late). Previous studies have suggested that ensembles that encode for sensory stimuli and well-trained motor actions are stable (Katlowitz et al., 2018; Lutcke et al., 2013; Perez-Ortega et al., 2021). However, recent studies have demonstrated that even engrams that have been hypothesized to be fairly stable can drift across days or show a combination of drifting populations coexisting with a much smaller proportion of stable ensembles (Marks and Goard, 2021; Driscoll et al., 2017; Gonzalez et al., 2019). Surprisingly, in our current study the observed representational drift occurs on a much faster time scale-minutes rather than days- and remains persistently unstable across days.

Although we have not investigated the mechanisms by which such a migrating ensemble could arise, several hypotheses can be posed based on previous studies. Within the amygdala, sparse neuronal ensembles are allocated to fear engrams based on stochastic neuronal activity and gene transcription levels (Cowansage et al., 2013; Kim et al., 2014; Yiu et al., 2014). If this represents a conserved mechanism for ensemble allocation, it is possible that early shock responsive cells are those with intrinsically high levels of activity, while shock experiences later within the session activate a distinct population based on activity levels at that specific time. One could also envision a mechanism that excludes previously shock responsive neurons to facilitate “migration” of the ensemble within a neuronal population. For example, once neurons have been allocated as shock responsive, activity-dependent plasticity mechanisms such as long-term synaptic depression or long-term depression of intrinsic excitability could reduce the probability previously activated neurons from participating in the shock responsive ensemble at later time points. Together, such mechanisms could explain our data that ∼90% of neurons meet activation criteria at some point during a foot-shock session, but that the proportion of neurons active during any one period (early, middle, late) is relatively low, and that there is little overlap between ensembles at these different time points. The functional consequences of this process could be advantageous in that it enables a maximal number of cells to participate in the shock representation and could therefore maximize the neural substrates for computational processes that would require strong representation of the aversive experiences including associative learning, for example. This apparent instability could also allow ensembles to “try” different configurations while maintaining stable representations and simultaneously accounting for variations in internal states and external environment (Kappel et al., 2015). Moreover, representational drift may provide the required flexibility to integrate new information and update memories in real time (Mau et al., 2020). Furthermore, recent studies have shown that although drifting representations lead to loss of information contained within individual neuron dynamics, information can be decoded from the population level even when many of the responsive neurons differ from those activated days before (Driscoll et al., 2017; Rule et al., 2019).

Repeated stress exposures can lead to habituation of behavioral and neuroendocrine responses to stress (Grissom and Bhatnagar, 2009). It is postulated that this habituation is protective against the toxic effects of stress exposure (Herman, 2013). Our longitudinal recordings of PL neurons across 3 days of stress exposure reveal evidence for habituation at the neuronal population activity level. Our data suggest that habituation to stress is associated with several adaptations within the PL. 1) The within session increase in net neuronal activity observed within the first shock session is blunted on days 2 and 3 of shock exposure, which is associated with an increase in the proportion of cells with inhibitory responses, and a reduction in proportion of cells with excitatory responses, on days 2 and 3. 2) Our PCA revealed that by day 3, the centroid distance between the representations of the stressor across each epoch (early and late) is reduced, which suggests that the representation of the stressor itself becomes more similar across days. Thus, early and late shock exposures are represented by “more similar” neuronal dynamics on day 3 than on day 1. 3) The increase in the degree of similarity between early and late epochs of the shock session on day 3 is also reflected in our graph theory analysis, which shows that repeated stress leads to blunting of the increase in shock-induced synchronized activity observed between early and lake shocks on day 1. Specifically, on day 1, correlated activity increases strongly between shock responsive neurons within the session, but this within session increase is absent on day 3. Moreover, after Day 1, within-session changes in PL neuronal activity become blunted to the point that they cannot be decoded from each other when using a decoder approach that is able to predict the identity of the stress epoch (early vs late) on Day 1. Likely, the decoder relies on features such as neuron-neuron correlations and the proportion of neurons that respond to the shock, among other things that are different early in the session vs late in the session in Day 1 for reliable decoding. As the features the decoder relies on become more similar between epoch on the other days, the decoder fails. Overall, our data support the notion that PL PFC neuronal dynamics undergo a form of habituation to repeated homotypic stress via a progressive blunting of within session increases in overall activity and correlated activity to the point of becoming indistinguishable from each other (based on the failure of PL activity to decode early vs late shock on days 2 and 3). These data are supported by previous studies demonstrating that habituation to repeated stressors requires the dorsal PFC function (Adamec et al., 2012; Katz et al., 2009; Weinberg et al., 2010; Diorio et al., 1993).

Despite the novel findings presented herein, there are several limitations to the current study. We only studied the neural activity of male mice during repeated foot-shock stress. Future studies examining potential sex differences in PL neuronal responses to stress will be important to conduct. Examination of the neuromodulatory factors contributing to ensemble drift and the genetic identity or projection targets of neurons involved in the adaptation response are all intriguing questions that could be visited in future studies.

Impairments in adaptation to stress are thought to lead to psychiatric disorders such as depression, anxiety disorders, and PTSD. Furthermore, these disorders are often associated with impaired cognitive flexibility which can further hinder recovery. Stress adaptation is a dynamic process which requires broad circuit-wide adaptations (Russo et al., 2012). We postulate that high levels of rapid and persistent representational drift in the PL may facilitate computational processes such as associative learning by maximizing the neuronal population eligible to participate in such processes by virtue of having, at least briefly, contributed to the sensory representation of the stressor. We also provide new insight into population level habituation within the PL and suggest that stabilizing neuronal representations of shock experience over time may be an important neuronal correlate of behavioral and endocrine habituation to repeated homotypic stress. It is possible that drift-induced maximization of populations participating in the shock representation may be required for optimal habituation. Overall, these data may provide new insight into the mechanisms by which PFC responds to stress and how it could participate in important stress adaptive processes such as habituation.

Funding

These studies were supported by grants K08MH126166 (L.R.V.) and R01MH107435 (S.P.) and a the Brain & Behavior Research Foundation (L.R.V.).

CRediT authorship contribution statement

Sachin Patel: Conceptualization, Methodology, Writing – original draft, Writing – review & editing, Visualization. Keenan Johnson: Investigation, Visualization. Danielle Adank: Investigation. Luis E. Rosas-Vidal: Conceptualization, Methodology, Investigation, Writing – original draft, Writing – review & editing, Visualization.

Declaration of competing interest

S.P. is a scientific consultant for Psy Therapeutics, Janssen Pharmaceuticals, and Jazz Pharmaceuticals unrelated to the present work. All other authors declare no conflicts of interest.