Abstract
The brainstem dorsal periaqueductal gray (dPAG) has been widely recognized as being a vital node orchestrating the responses to innate threats. Intriguingly, recent evidence also shows that the dPAG mediates defensive responses to fear conditioned contexts. However, it is unknown whether the dPAG displays independent or shared patterns of activation during exposure to innate and conditioned threats. It is also unclear how dPAG ensembles encode and predict diverse defensive behaviors. To address this question, we used miniaturized microscopes to obtain recordings of the same dPAG ensembles during exposure to a live predator and a fear conditioned context in male mice. dPAG ensembles encoded not only distance to threat, but also relevant features, such as predator speed and angular offset between mouse and threat. Furthermore, dPAG cells accurately encoded numerous defensive behaviors, including freezing, stretch-attend postures, and escape. Encoding of behaviors and of distance to threat occurred independently in dPAG cells. dPAG cells also displayed a shared representation to encode these behaviors and distance to threat across innate and conditioned threats. Last, we also show that escape could be predicted by dPAG activity several seconds in advance. Thus, dPAG activity dynamically tracks key kinematic and behavioral variables during exposure to threats, and exhibits similar patterns of activation during defensive behaviors elicited by innate or conditioned threats. These data indicate that a common pathway may be recruited by the dPAG during exposure to a wide variety of threat modalities.
SIGNIFICANCE STATEMENT The dorsal periaqueductal gray (dPAG) is critical to generate defensive behaviors during encounters with threats of multiple modalities. Here we use longitudinal calcium transient recordings of dPAG ensembles in freely moving mice to show that this region uses shared patterns of activity to represent distance to an innate threat (a live predator) and a conditioned threat (a shock grid). We also show that dPAG neural activity can predict diverse defensive behaviors. These data indicate the dPAG uses conserved population-level activity patterns to encode and coordinate defensive behaviors during exposure to both innate and conditioned threats.
Introduction
In order to survive, animals must rapidly display adaptive defensive behaviors during exposure to threats. The rodent defensive behavioral ethogram includes diverse behaviors, including stretch-attend postures to investigate threats as well as freezing and escape (Perusini and Fanselow, 2015). The brainstem dorsal periaqueductal gray region (dPAG) has been long recognized as a key downstream effector region coordinating innate defensive behaviors (Gross and Canteras, 2012; Perusini and Fanselow, 2015; Tovote et al., 2016). The dPAG's role in escape has been particularly well studied. dPAG optogenetic stimulation causes escape even in the absence of threats, and single dPAG cells are active during escape from threats, such as looming stimuli, which simulate aerial predators (Evans et al., 2018). Recent evidence also shows the dPAG affects other defensive behaviors. Optogenetic activation of the dPAG has been reported to cause freezing (Deng et al., 2016), and dPAG NMDA lesions decrease predator-induced freezing and risk-assessment behaviors, such as stretch-attend and crouch-sniffing postures (Sukikara et al., 2010). Furthermore, intra-dPAG pharmacological manipulations have implicated transmission of cholecystokinin (Vazquez-Leon et al., 2018), nitric oxide (Morato et al., 2004), and glutamate (Gomes et al., 2014) in avoidance of anxiety-inducing open spaces in the elevated plus maze. These reports demonstrate that the dPAG is a key structure mediating a wide array of innate defensive responses, including not only escape, but also freezing, risk assessment, and avoidance of open spaces.
Although the dPAG has been more extensively studied for its role in escape from innate threats, it also has been implicated in contextual fear conditioning. For example, an increase in the activity marker c-fos has been reported in the dPAG during retrieval of contextual fear (Carrive et al., 1997), and a decrease in contextual fear-induced freezing was observed following intra-dPAG injections of cannabinoid receptor agonist (Resstel et al., 2008). Similarly, inhibition of NMDA or nitric oxide transmission in the dPAG before fear retrieval also decreased freezing to fear conditioned context (Aguiar et al., 2014). Conversely, intra-dPAG injections of NMDA (Mochny et al., 2013) and CRF (Borelli et al., 2013) increased freezing to a fear conditioned context. These data show the dPAG can affect conditioned responses during contextual fear retrieval.
Together, the reports discussed above indicate that the dPAG controls diverse defensive behaviors during exposure to innate and conditioned threats. Despite this large body of work, there are relatively few single-unit recordings of dPAG neuronal activity. A prior report shows that dPAG cells are activated during escape or stretch-attend postures induced by a live predator (Deng et al., 2016). Similarly, increased dPAG activity was also observed during escape from looming stimuli (Evans et al., 2018). Despite these pioneering reports, several key questions remain unanswered. Do dPAG cells respond similarly to innate and conditioned threats? Does the dPAG encode any variables related to threat imminence, such as distance to threat or onset of predator movement? How does the dPAG population code, represent, and predict defensive behaviors?
To address these questions, we obtained recordings from large dPAG ensembles during exposure to a live predator and to a fear conditioned shock grid using miniaturized microscopes to record calcium transients. We show that the dPAG has a shared neural activity pattern to represent kinematic and defensive behavioral variables during exposure to innate and conditioned threats. We also show that dPAG cell responses to ongoing defensive behaviors can be used to separate them into discrete functional clusters, and cluster membership is maintained across threat modalities. These results show that innate and conditioned threats can engage similar activity patterns in the brainstem, which may result in the coordination of vital defensive behaviors to ensure survival to different types of threat.
Materials and Methods
Mice
Male mice between 2 and 5 months of age of the C57BL/6J strain (The Jackson Laboratory stock #000664) were used for all experiments. Mice were maintained on a 12 h reverse light-dark cycle with food and water ad libitum. Sample sizes were chosen based on previous behavioral studies with miniaturized microscope recordings on defensive behaviors, which typically use 6-10 mice per group. All mice were handled for a minimum of 5 d before any behavioral task. All procedures conformed to guidelines established by the National Institutes of Health and have been approved by the University of California, Los Angeles Institutional Animal Care and Use Committee.
Rats
Male Long-Evans rats (250-400 g) were obtained from Charles River Laboratories and were individually housed on a standard 12 h light-dark cycle and given food and water ad libitum. Rats were only used as a predatory stimulus. Rats were handled for several weeks before being used and were screened for low aggression to avoid attacks on mice. No attacks on mice were observed in this experiment.
Surgeries
Eight-week-old mice were anaesthetized with 1.5%-3.0% isoflurane and placed in a stereotaxic apparatus (Kopf Instruments). AAV9.Syn.GCaMP6s.WPRE.SV40 were packaged and supplied by UPenn Vector Core. After performing a craniotomy, 100 nl of virus at a titer of 5 × 1012 was injected into the dPAG (coordinates in mm, from skull surface): −4.20 anteromedial, −0.85 lateral, −2.3 depth, 15 degree angle. Five days after virus injection, the animals underwent a second surgery in which two skull screws were inserted and a microendoscope was implanted above the injection site. A 0.5-mm-diameter, ∼4-mm-long gradient refractive index (GRIN) lens (Inscopix) was implanted above the dPAG (−2.0 mm ventral to the skull surface). The lens was fixed to the skull with cyanoacrylate glue and adhesive cement (Metabond; Parkell). The exposed end of the GRIN lens was protected with transparent Kwik-seal glue, and animals were returned to a clean cage. Two weeks later, a small aluminum base plate was cemented onto the animal's head atop the previously formed dental cement. Animals were provided with analgesic and anti-inflammatory (carprofen).
Rat exposure assay
On day 1, mice were habituated to a white rectangular box (70 cm length, 26 cm width, 44 cm height) for 20 min. Twenty-four hours later, mice were exposed to the same environment but in the presence of a Toy Rat for 20 min. Mice were then exposed to an adult rat in this environment on the 2 following days. The rat was secured by a harness tied to one of the walls and could freely ambulate only within a short radius of ∼20 cm. The mouse was placed near the wall opposite to the rat and freely explored the context for 20 min. No separating barrier was placed between the mouse and the rat, allowing for close naturalistic encounters that can induce a variety of robust defensive behaviors.
Contextual fear conditioning test
To better evaluate a broader species-specific defense repertoire in face of a conditioned stimulus, we used a modified version of the standard contextual fear conditioning method (Schuette et al., 2020). Preshock, fear conditioning, and retrieval sessions were performed in a context (70 cm × 17 cm × 40 cm) with an evenly distributed light intensity of 40 lux and a Coulbourn shock grid (19.5 cm × 17 cm) set at the extreme end of the enclosure. Forty-eight hours after rat exposure, mice were habituated to this context and could freely explore the whole environment for 20 min. On the following day, the grid was activated, such that 0.4 mA foot shocks were delivered for 1 s whenever the mouse fully entered the grid zone. Escapable shocks more closely resemble naturalistic, localized threats, which are typically escapable and are not occurring in the entire environment simultaneously. Twenty-four hours later, retrieval sessions were performed in the same enclosure but without shock. Mice could freely explore the context for 20 min during Preshock habituation, fear conditioning, and retrieval sessions.
Behavior and miniscope video capture
All videos were recorded at 30 frames/s using a Logitech HD C310 webcam and custom-built head-mounted UCLA miniscope (Aharoni and Hoogland, 2019). Open-source UCLA Miniscope software and hardware (http://miniscope.org/) were used to capture and synchronize neural and behavioral video (Cai et al., 2016).
Perfusion and histologic verification
Mice were anesthetized with Fatal-Plus and transcardially perfused with PBS followed by a solution of 4% PFA. Extracted brains were stored for 12 h at 4°C in 4% PFA. Brains were then placed in sucrose solution for a minimum of 24 h. Brains were sectioned in the coronal plane in a cryostat, washed in PBS, and mounted on glass slides using PVA-DABCO. Images were acquired using a Keyence BZ-X fluorescence microscope with a 10× or 20× air objective.
Data analysis
All analysis was performed using custom-written code in MATLAB and Python.
Miniscope postprocessing and coregistration
Miniscope videos were motion-corrected using the open-source UCLA miniscope analysis package (https://github.com/daharoni/Miniscope_Analysis) (Aharoni and Hoogland, 2019). They were spatially downsampled by a factor of 2 and temporally downsampled by a factor of 4, and the cell footprints and activity were extracted using the open-source package Constrained Nonnegative Matrix Factorization for microEndoscopic data (CNMF-E; https://github.com/zhoupc/CNMF_E) (Zhou et al., 2018). Neurons were coregistered across sessions using the open-source probabilistic modeling package CellReg (https://github.com/zivlab/CellReg) (Sheintuch et al., 2017).
Behavior detection
To extract the pose of freely behaving mice in the described assays, we implemented DeepLabCut (Mathis et al., 2018), an open-source convolutional neural network-based toolbox, to identify mouse nose, ear, and tail base xy coordinates in each recorded video frame. These coordinates were then used to calculate velocity and position at each time point, as well as classify defensive behaviors in an automated manner using custom MATLAB scripts. Freezing was defined as epochs when head and tail base velocities fell to <0.25 cm/s for a period of 0.33 s. Stretch-attend posture was defined as epochs when the distance from mouse nose to tail base exceeded a minimum threshold for a period of 0.5 s. Approach and escape were defined as epochs when the mouse moved with a velocity >3 cm/s, respectively, away from or toward the rat. Grooming, rearing, or other behaviors were not observed during frames categorized as approach, escape, freezing, and approach.
Factor analysis
In these freely behaving assays, there are a diversity of behaviors and quantifiable behavioral metrics. To assess whether behaviors across assays exhibited similar motifs, we quantified how the 18 behavioral metrics covaried consistently across assays. We achieved this by performing dimensionality reduction to identify five latent factors that reflect behavioral synergies for each assay. Specifically, we performed factor analysis using the behavioral statistics as input variables from a separate large cohort of mice (n = 44) that did not undergo surgical procedures (see Fig. 1G). This large cohort was necessary to do factor analysis, which was performed using MATLAB's factoran function individually for each environment using default parameters (i.e., with the varimax projection) to find the 5 factors with highest variance. The data were recorded from 44 mice which did not have a Miniscope attached. The input to factoran was therefore a 44 × 18 (mice × variables) array of data for each environment.
z-score normalization
We z-scored the neural activity dF/F for all analyses in the study (see Figs. 3–8), as well as the General Linear Model (GLM coefficients in Gaussian Mixture Model (GMM) clustering (see Figs. 7, 8). The z score is calculated by subtracting the population mean from the raw score and dividing this difference by population SD. The z score dF/F of an individual neuron (
GLM
A GLM was fit for each cell. Each GLM mapped kinematic or behavioral variables (inputs) to the cell's z-scored calcium activity (output). The GLM therefore quantifies how kinematic and behavioral variables are represented by each recorded dPAG cell. We used the GLM to quantify how dPAG activity reflects both mouse and rat kinematics, as well as defensive behaviors. In total, there were four kinematic input variables for rat (distance to the rat, mouse velocity, rat velocity, angle from mouse's head to the rat) and four behavioral input variables (the occurrence of approaching, stretching, escaping, and freezing). There were three kinematic input variables for shock (distance to the shock grid, mouse velocity, angle from mouse's head to the shock grid) and four behavioral input variables (the occurrence of the four behaviors).
We processed these input variables to generate additional features before performing linear regression to predict z-scored calcium activity (see Fig. 3F). For kinematic input variables, we computed nonlinear polynomial features so that behaviors could be nonlinearly related to the z-scored calcium activity. In particular, if xi is the ith kinematic input variable, we also used quadratic and/or cubic features xi2 and xi3 as input features to the GLM. The nonlinear features were determined via cross-validation and were second degree polynomials for Rat 1 and Toy Rat Assays, and third degree polynomials for all other assays. Further, to enable historical kinematics to affect the neural data, the kinematics were convolved with a causal kernel. This kernel had the following form:
These kernels are non-zero when x is between
Discrete behaviors were binary, labeled as 1 at all times they occurred and 0 otherwise. To enable the model to detect relations between behaviors and neural activity before and following the behavior, each binarized behavior was convolved with the same kernel functions as for the kinematics. We used 14 total bases. In addition to the 7 bases used for kinematics, we also had bases activating in ([−0.25, 0], [−0.42, −0.04], [−0.7, −0.075], [−1.15, −0.14], [−1.9, −0.24], [−3.15, −0.41], [−5, −0.69]) seconds, with peaks at (−0.075, −0.14, −0.24, −0.41, −0.69, −1.15, −1.9) seconds. These GLM coefficients were only computed at times when the mouse was exhibiting a defensive behavior.
In all instances, GLMs were trained separately for kinematic variables and for discrete behavioral variables. This is because discrete behavioral variables are only active during particular times in the experiment, and it also decreases the high dependency between kinematics and behaviors from the model. In addition to this, we also regressed out the variable distance to threat from dPAG activity before fitting the behavioral GLM. We did so by computing the expected dPAG activity based on the mouse's distance to threat, computed via the kinematic GLM, and subtracted this activity from the recorded dPAG activity. This adjusted activity was used to fit the behavioral GLM used in Figures 3, 4, 6–8. As such, any significant weights in the behavioral GLM cannot be explained as a byproduct of encoding of distance to threat.
These input variables, processed via convolution with kernels and denoted here as x1, x2, etc., were modeled to produce the cell's calcium fluorescence
Throughout the paper, we trained GLMs using all 7 kernels, or only 1 kernel. When assessing the overall model fit and the relative contribution of kinematic and discrete variables to calcium activity (see Figs. 3G,H, 4D,E), we trained the GLM with all 7 kernels because it resulted in the best model fit with the highest cross-validated r2. Because there are many parameters, we used ridge regression to prevent overfitting, with the ridge hyperparameter penalty weight being found via cross-validation.
When comparing GLM weights for each kinematic and behavior input within or across assays, it is not straightforward to use a GLM with 7 kernels, since each input variable is associated with multiple weights that may differently weigh particular kernels in a dataset dependent way. In these instances, we trained a GLM with only 1 kernel, so each input behavior was characterized by a single weight. The kernel was activated in the window [0, 0.25] seconds, with the peak at 0.075 s. We found that ridge regression had a negligible effect in this GLM with relatively few parameters, and therefore did not incorporate ridge regression in 1 kernel GLMs. This enabled direct comparison of GLM coefficients across assays. This GLM was used in Figures 3I-K, 4F–H, 6–8.
For the analyses using GLM in Figures 3H-K and 4E–H, all mice and all cells were used in each assay. In total, there were 7 mice and 715 cells in Rat 1 assay, 6 mice and 625 cells in Rat 2 assay, 8 mice and 737 cells in Fear Acquisition assay, and 6 mice and 621 cells in Fear Retrieval assay.
Relative contribution in GLM
To better understand how well calcium activity could be predicted by behavioral variables in GLM, we used the Pearson correlation coefficient squared, r2 (Engelhard et al., 2019). This coefficient was measured between GLM reconstructions and recorded calcium activity. A higher r2 indicates that the dPAG activity is better predicted by the input variables (kinematics or behaviors). While r2 is equal to the coefficient of determination R2 for training data, this is not the case for testing data.
We computed the
We also computed the reduction in cross-validated r2 of the GLM when removing a variable in the GLM. The reduction in r2 when removing input variable i was calculated as follows:
For Figures 3H and 4E, we computed the relative contribution of each kinematic and behavior input variable to the GLM. The relative contribution for each input variable by variables was calculated following the protocol of Engelhard et al. (2019). This is a normalized version of reduction in
Classification of significant cells for variables in GLM
Using the trained GLMs, we also identified whether cells significantly represented particular kinematic or behavioral variables. In particular, we classified cells that were significantly positively or negatively modulated for the following kinematic and behavioral variables: distance to threat (rat or shock grid), mouse speed, rat speed (Rat Assays only), angle, approaching, stretching, escaping, and freezing (see Figs. 3J, 4G).
To do so, we first bootstrapped the GLM coefficients of each cell by training a 1 kernel GLM with randomly resampled input variables. These bootstrapped GLM coefficients were used to build a null distribution of GLM coefficients for each input variable. A cell was classified as a significantly positive cell for a variable if it had a positive GLM coefficient and its value was >95% of the bootstrapped coefficients (p < 0.05). Similarly, a cell was classified as a significantly negative cell for a variable if it had a negative GLM coefficient and its value was <95% of the bootstrapped coefficients (p < 0.05). These classes of cells were also used in the analyses of Figures 5 and 6.
Distinct representations of distance to threat and defensive behaviors in dPAG
We performed additional analyses to show that dPAG's representation of defensive behaviors is not explained by an encoding of distance to the threat. These analyses enable us to show that defensive behaviors, such as approach, stretch, freeze, and escape, are not explained solely by where they occur relative to the threat.
We compared the activation of cells that encode behaviors during time points in which that behavior was displayed compared with time points without the behavior. Importantly, this comparison was done in the same range of distance to threat (see Fig. 5A). For each behavior, we analyzed events that occurred within a position window during which the behavior occurred most frequently. For approach and freeze, we chose 0-20 cm close to the left “safe” wall, furthest from the threat. For stretch and escape, we chose the position window 20 cm closest to the threat. This allowed us to compare, for example, whether freeze-encoding cells were more activated during freeze rather than nonfreeze time points even when restricting the comparison only to data points that occurred within the same 20 cm segment of length within the environment. This ensured that the differences would not be driven by distance to threat, as the same distance to threat range was used when comparing behavioral and nonbehavioral time points. If this analysis shows, for example, that freeze cells are more activated during freeze time points compared with nonfreeze time points, it would indicate that these cells are indeed encoding freezing, rather than simply representing distance from threat.
For each defensive behavior, we quantified the mean absolute value of the z-scored dF/F whenever the behavior occurred within the position window (see Fig. 5A). The absolute value of the mean z-scored dF/F was computed across a window [−2.5, 5] seconds centered around behavior onset. The mean was then taken across all cells significantly modulated for the defensive behavior. We also defined a nonbehavior event as any period where the mouse stayed within the position window for at least 4 s but did not exhibit the defensive behavior. We computed the mean absolute value of the z-scored dF/F over these 4 s windows. We then compared the mean during behavior and nonbehavior events (see Fig. 5B). If the dPAG modulation for defensive behaviors was only explained by the mouse's distance to the threat and not a unique representation of the defensive behavior, we would expect these mean activities to be statistically similar for both behavioral and nonbehavioral event times. Otherwise, the dPAG represents defensive behaviors beyond the distance to the threat.
We also analyzed whether the cell activation during defensive behaviors is not dependent on the location in which the behavior is displayed. To do so, we quantified the correlation between (1) the distances at which pairwise behaviors occur relative to the threat and (2) the similarity of the pairwise behavior triggered neural responses (see Fig. 5C,D). For a defensive behavior, we analyzed all pairwise occurrences of the behavior in the following way: for each pairwise occurrence, for example, two escapes occurring at distances d1 and d2 from the threat with dF/F activity traces f1 and f2, we computed the difference in distance to threat, a1 = |d1-d2|, and the Pearson correlation between the dF/F traces, r1 = corr(f1, f2). This pair, (a1, r1), comprised one data point in a scatter plot (see Fig. 5C) with the distance between pairwise threats on the x axis, and the correlation in dF/F on the y axis. We performed these calculations on all pairwise instances of a defensive behavior. We then fit a linear model to fit the scatter plot data. If the dPAG representation of the defensive behavior was modulated by the distance to the threat, we would expect high correlations ri when the distances ai are small, and small correlations when the distances are large, resulting in a line with negative slope. In contrast, if the slope of the line is close to zero, then the defensive behaviors are similarly correlated regardless of where they occur relative to the threat (see Fig. 5D).
For the analyses showing the distinct representations of distance to threat and defensive behaviors, only cells that were significantly modulated for each defensive behavior were used (see Fig. 5B,D).
Cosine similarity of GLM coefficients between neurons
We used cosine similarity to quantify the similarity of GLM coefficients across assays. This analysis allows us to use a single number to quantify how similar the GLM weight coefficients are across assays for a given cell. In summary, the analysis involved generating a vector of GLM weights for a cell in one assay (the number of dimensions in this vector was equal to the number of input variables in the GLM). Then, a second vector was created, corresponding to the GLM weights for that cell in a different assay. The angle θ between those two vectors was calculated. If the distribution of GLM weights for a cell is very similar across assays, the angle between these two vectors will be small. Conversely, if the distribution of weights is not similar, the angle will be large. As the cosine of an angle of zero produces the maximum result (cos(0) = 1), a similar distribution of GLM weights across assays will produce a small angle between the vectors described above, which in turn will generate a large cosine similarity. Thus, larger values of cosine similarity correspond to increased similarity of GLM weights across assays.
To compute the cosine similarity of GLM coefficients across assays, we formed a vector of GLM weights across all cells. We then compared these GLM weight vectors across assays (see Fig. 6B–D). We computed the cosine of the angle between GLM weight vectors across assays in the following way. If r and s denote the GLM weight vector in Assay 1 and Assay 2, respectively, then their cosine similarity is as follows:
As cosine similarity compares GLM coefficients, we used a GLM with one kernel for each defensive behavior and threat-related kinematic variable, leading to a single weight. For the analyses showing cosine similarity between two groups of cells in this study, we only compared the cells that were coregistered between two assays (see Fig. 6C,D).
Gaussian mixture model clustering
To identify clusters of cells related to the four behaviors, we used a Gaussian mixture model (GMM). This analysis allowed us to investigate whether dPAG cells belong to functional clusters based on the distribution of their GLM weights, and to study whether cluster assignment was conserved across assays. The underlying goal of this analysis is to study whether cells respond to defensive behaviors similarly across threatening assays. To do so, we applied a GMM to the matrix of 1 kernel GLM weights of behavioral variables (approaching, stretching, escaping, and freezing). This was an N × 4 matrix, where N is the number of cells. We used the GMM function in sklearn.mixture with diagonal covariance matrices and 1000 max iterations. This function implements the expectation-maximization algorithm to learn the clusters. This made a GMM where the major axes of the Gaussians are parallel to the axes of the feature space, which uses less parameters but achieves more flexibility than that of the k-means algorithm. A vector of weights of four behaviors can be modeled by the following:
To quantify model fit, we used the Bayesian Information Criterion (BIC), which is defined by the following:
To assess the significance of the clustering model, we shuffled the GLM weights 10,000 times across behavioral variables and across neurons. After each shuffling iteration, we repeated the clustering and calculated the log-likelihood of the entire clustering model. The distributions of log-likelihood for shuffled data were compared with the log-likelihood of the clustering model using the real data. Clustering on the real data achieved significantly higher likelihood values (p < 0.001, Wilcoxon rank-sum test, all assays).
In GMM clustering, a mouse was excluded if it never exhibited one of the four defensive behaviors, since this would heavily bias clustering. We therefore used 6 mice in Rat 1 and 5 mice in Fear Retrieval, since there was a mouse that did not exhibit freezing in these two assays (see Fig. 7).
Behavioral clusters decoding using GLM coefficients
We investigated whether clusters derived in each assay are preserved across different assays, quantifying whether a coregistered cell in two assays would be grouped into the same type of cluster across assays. To do so, we fit GMM cluster parameters from one assay and applied them to coregistered neurons recorded in a second assay. That is, we learned K clusters, each with GMM parameters
In this analysis, only mice that displayed all four behaviors in both assays were used, and only coregistered cells between two assays were counted. For example, there were 2 mice used in the decoding between Rats 1 and 2 assays since 3 of 5 mice did not exhibit significant freezing in one of these assays, and these 2 mice had 97 coregistered cells in total (see Fig. 8).
Behavior decoding using dPAG neural data
To investigate whether dPAG cells encoded a wide range of behaviors, we performed discrete classification of behaviors performed using multinomial logistic regression. This decoding analysis differs from the GLM because the behavior is decoded from the entire population of recorded cells. This therefore quantifies the extent to which a behavior is represented in the population of dPAG cells. We used scikit-learn's LogisticRegression function (Pedregosa et al., 2011). We additionally regressed out distance to the rat or shock grid before training and classification. Only times where one of the four behaviors (approach, stretch, escape, freeze) occurred were included, with each time point treated as an individual data point. Sample weights of behaviors were balanced equally. For single-assay and across-assay prediction, training and validation were performed using fivefold cross-validation, with a minimum of 10 s between training and validation sets. Behavioral transitions were found between all epochs separated by this window. True and predicted labels for validation sets were concatenated and used to create confusion matrices. Significance testing was performed using permutation tests, where the labels of training sets were shuffled before model-fitting (100 trials), with the null hypothesis that the means of confusion matrix entries were not higher than expected by chance (∼25%).
Prediction of future escape behavior using a neural network model of dPAG activity
We also studied whether dPAG cells were able to predict the occurrence of escapes in the future, that is, if dPAG activity can be used to predict escapes several seconds ahead of actual behavioral onset. We used a three-layer neural network (trained with Keras) to predict escape behavior in the future using cell activity in the rat and shock (shock and extinction days) assays. The first layer used a rectified linear activation function (ReLU), and the remaining layers used a sigmoid activation function. The loss function was the binary cross-entropy loss. Testing data were interleaved with training data using a two-split time series cross-validation. Segments of the two data types were 10 frames apart from each other, with each time point treated as an individual data point. Behavioral transitions were found between all epochs separated by this window. For future event predictions, only escape behavior that did not coincide with the present escape occurrence was considered (any overlapping points were considered as “not escape”). Because there are many more time points when the mouse was not escaping compared with escaping, we computed weighted accuracy to equally weight the classes. Weighted accuracy was calculated using equal weights (of 0.5) for escape and not escape accuracies, keeping a chance level of 0.5.
Cross-assay constrained correlation analysis (CoCA)
We developed CoCA to identify a shared projection of dPAG activity across assays that correlates strongly with behavior within each assay. This analysis was done to identify which features of activity (e.g., encoding of escape or distance to threat) are most strongly conserved across assays in the dPAG population. Because the behaviors and kinematics are different within each assay, this necessitated an optimization approach that jointly considered which behavioral features are important in each assay (across mice) and how the neural data relate to these behavioral features (across assays). We denote calcium fluorescence neural data as X ∈ RkxT and externally observed behavioral data as Y ∈ RpxT, where k is the number of recorded cells for the corresponding mouse, shared across assays, p is the chosen number of behavioral variables for the corresponding assay, shared across mice, and T is the length of a recording session, unique for each session, but shared between neural and behavioral data for the same session. Behavioral variables contained both continuous kinematic variables (e.g., speed and distance from rat) as well as binary defensive behavior variables (e.g., the occurrence of freezing, approach, stretch, and escape). All variables were normalized to zero mean and unit variance.
In order to find a common linear projection of threat across mice and assays, we performed the following optimization with mouse IDs i = 1, 2, …7 and assay IDs j = Shock Grid, RAT. Calcium fluorescence traces of dPAG cells for mouse i were linearly combined after multiplying each cell with weights n1i to nki, where k is the number of cells that were coregistered in both assays. Taking the dot product of the calcium activity for mouse i in assay j, given by Xi,j, and the weights ni = [ni1…k] defined a neural projection for mouse i and assay j, given by Ni,j = (ni)TXi,j. For each mouse, the weights, ni, were the same across assays, so that each cell had the same weight in both assays. The behavioral variables for the shock grid (e.g., x and y position, freezing, stretch-attend posture, speed, etc.) were linearly combined with a set of weights b1 to b7 (as 7 behavioral variables were used for the shock grid). These weights, bShock Grid = [b1Shock Grid, b2Shock Grid, …, b7Shock Grid] were conserved across all mice. Linearly combining the Shock Grid assay behavioral variables resulted in a behavioral projection for mouse i and assay Shock Grid, given by Bi,Shock Grid = (bShock Grid)TYi,RAT. Similarly, 8 behavioral variables from the Rat Assay were linearly combined to produce a behavioral projection Bi,RAT = (bRAT)TYi,RAT using weights bRAT = [b1RAT, b2RAT, …, b8RAT]. We chose the neural weights, ni, and the behavioral weights, bShock Grid and bRAT, to optimize the correlations across all mice and assays as follows:
In order to test whether correlations of testing data were better than expected by chance, correlations were computed between projected behavioral data (using projections fit by training data) and random projections of neural data (1000 trials). We emphasize that these correlations were applied to the testing data; therefore, it was possible for a random projection to have higher correlation than the CoCA projection. Here a one-tailed test was used, as we expected the random projections to have lower correlations between behavioral and neural projections than the ones produced by the actual analysis.
Experimental design and statistical analysis
We have not preregistered this study. Male mice were used. Forty-four mice without miniaturized microscope implants were used for Figure 1. Eight mice with miniaturized microscope implants were used for all other experiments. There were no “control” groups because there were no treatments or manipulations done to any mice. Full statistical details for each comparison, including exact p values, can be found in Results. Significance values are included in the figure legends and Results. Unless otherwise noted, all statistical comparisons were performed by either nonparametric Wilcoxon rank-sum or signed-rank tests. SEM was plotted in each figure as an estimate of variation of the mean. Correlations were calculated using Pearson's method. Multiple comparisons were corrected with the false discovery rate method. All statistical analyses were performed using custom MATLAB scripts.
Results
Exposure to a live predator or a shocking grid induce a conserved topology of defensive behaviors
In order to study the effect of innate and conditioned threats on dPAG activity, we used two assays in which these two threat categories could be compared in similar conditions. In the rat exposure test, a mouse was placed in a long rectangular box containing a live predator (an adult rat). The rat was tethered with a harness to one side of the environment, and could only freely ambulate near the edge of the right side of the environment (Fig. 1A). Mice were exposed to the rat in two 20 min sessions separated by 24 h. In the shock grid assay, the mouse was placed in a similarly long rectangular box with a shock grid at the edge of the right side (Fig. 1B). During the Fear Acquisition day, the mouse received one 0.4 mA foot shock every time it went onto the shock grid. Twenty-four hours later during the contextual Fear Retrieval session, the mouse was placed in the same box in the absence of shocks (Fig. 1C). Mice spent most of the time avoiding both threats by predominantly exploring the corners of the environment located furthest away from the rat or the shock grid (Fig. 1D). In both assays, mice displayed similar defensive behaviors, consisting of escape, risk-assessment approaches, stretch-attend postures, and freezing. In both assays, these behaviors had a similar spatial distribution in the environment, as mice froze far away from the threat, and exhibited fast escape responses mostly from locations near the threat (Fig. 1D). This spatial distribution of behaviors matches the description of influential models, such as the predatory imminence theory (Perusini and Fanselow, 2015), indicating that these assays induced defensive behaviors in a naturalistic manner. Mice displayed increased freezing and risk-assessment related stretch-attend postures during exposure to threat assays compared with controls assays (Fig. 1E). Furthermore, relative to Toy Rat and Preshock controls, mice showed decreased approach velocity and increased escape velocity to threatening stimuli (Fig. 1E; n = 44 mice, Wilcoxon signed-rank test, # freezes toy/rat p = 4.18 × 10−9, Preshock/Fear Retrieval p = 3.37 × 10−10; # stretches toy/rat p = 8.33 × 10−6, Preshock/Fear Retrieval p = 0.0016; approach speed toy/rat p = 2.74 × 10−12, Preshock/Fear Retrieval p = 1.27 × 10−8; escape speed toy/rat p = 0.0002, Preshock/Fear Retrieval p = 1.51 × 10−7). Thus, both the rat and the shocking grid induce defensive behaviors robustly relative to control assays. The Toy Rat controlled for a few features of the real rat, such as similar size, shape, and novelty. However, by necessity, the Toy Rat cannot serve as a control for predator odors or rat movement. The grid during the Preshock exposure controlled for all visual, tactile, and olfactory stimuli of the grid during Fear Acquisition and Fear Retrieval. Importantly, the Preshock grid was used as the control for the Fear Acquisition and Retrieval sessions, but not for the Rat assay.
Characterization of the Rat and grid assays. A, B, Scheme represents the Rat and grid assays. C, Summary of behavioral timeline. D, Distribution of behavioral measures along the length of each assay. Mice spent the majority of the time away from the threat (first row). Freeze bouts were concentrated in the locations furthest away from threats (second row), while risk-evaluation stretch-attend postures tended to occur more near the rat or distributed along the context for the shock grid assays (third row). Approaches started away from the threats (fourth row), while escapes tended to be initiated in the half of the environment near the threat (fifth row). The zone containing the threat is shown in shaded light red box (n = 44 mice, Wilcoxon signed-rank test, # freezes toy/rat p = 4.18 × 10−9, Preshock/Fear Retrieval p = 3.37 × 10−10; # stretches Toy Rat/Rat p = 8.33 × 10−6, Preshock/Fear Retrieval p = 0.0016; approach speed toy/rat p = 2.74 × 10−12, Preshock/Fear Retrieval p = 1.27 × 10−8; escape speed toy/rat p = 0.0002, Preshock/Fear Retrieval p = 1.51 × 10−7) E, Bars represent the mean number of freezes, stretches, approach, and escape speed for Toy Rat/Rat and Preshock/Fear Retrieval assays (n = 44 mice, Wilcoxon signed-rank test, ***p < 0.001). F, Matrix represents transition probabilities between different behaviors for each assay. Warmer colors represent higher probability. From approach (app), the most likely transition is to stretch (str). The highest transition probability from stretch is to escape (esc), and escape generally leads to freeze (frz). This transition probability matrix shows that mice started to approach the threat, did stretch-attend postures, then escaped, and then froze. This temporal ordering was conserved across assays (n = 44 mice). G, Behavioral metrics are segregated into similar factors across assays. A factor analysis was performed using several relevant metrics for rat and shock grid assays. Behavioral measures were generally assigned to the same factor across assays, showing a shared behavioral topography across assays. The top 5 factors summed together account for 81.5%, 70.5%, and 75.2% of total variance in Rat assay, Fear Acquisition, and Fear Retrieval, respectively (n = 44 mice).
Mice also showed a similar temporal structure in the order of these behaviors across assays. In both assays, mice would approach the threat, often followed by stretch-attend postures (Fig. 1F). After approaching the threat, mice escaped with high speed to the opposite corner, where they froze for several seconds before reinitiating this behavioral cycle. To identify whether behavior was similar across assays, we performed factor analysis on relevant behavioral measures in both assays. Factor analysis identified that behavioral measures segregated into factors that were largely conserved across assays (Fig. 1G). In both conditions, the highest variance factor had large weights for behavioral measures related to approach and escape. We also found separate factors that captured freezing, stretch-related behavioral measures, and escape vigor. The similarity of behavior factors across assays indicates that mice exhibit consistent behavioral motifs and similar covariation between behavioral measures in both rat exposure and shock grid. These data indicate that the rat and shock grid assays can be used to compare how the dPAG represents innate and conditioned threat in similar conditions. A large cohort of mice without microscope implants (n = 44) was used in Figure 1 to obtain enough statistical power to perform the factor analysis.
dPAG ensembles represent multiple threat-related features and defensive behaviors
We next studied which threat-related features were represented in the dPAG during rat exposure and shock grid assays. The same mice were exposed to both assays. An AAV vector was used to target expression of the calcium indicator GCaMP6s in the dPAG. The Syn promoter was used to ensure GCaMP6s expression pan-neuronally. We obtained high-quality recordings of dPAG cells and tracked the same cells across several days (Fig. 2F; n = 527 cells, for both, p < 0.001). Motion correction was applied to GCaMP6s fluorescence videos as described previously (Aharoni and Hoogland, 2019) (see Materials and Methods).
Characterization of recording quality in the dPAG. A, A viral vector was injected in the dPAG encoding GCaMP6s under the pan-neuronal syn promoter. A GRIN lens was implanted over the injection site and coupled to a miniaturized microscope to obtain recordings of calcium transients in the dPAG. B, Example traces of simultaneously recorded cells. C, Maximum projection of the dPAG FOV in a representative example mouse. D, Example imaging FOV with dPAG cells coregistered between Rat exposure 1 and 2. E, Bars represent the displacement of coregistered cells between Rat and Fear Acquisition assays as a fraction of the mean neuron diameter (n = 527). F, The peak-to-noise ratio (PNR) and mean peak amplitude correlation values were calculated for coregistered cells between Rat and Fear Acquisition assays. Cell identities were then shuffled within the 10 nearest neighbors 1000 times, and the same correlation measures were calculated for each iteration. The resulting bootstrap distribution was compared with the actual peak-to-noise and mean peak amplitude values, indicated with a blue arrow (n = 527; comparison of actual values to bootstrap distribution; p < 0.001).
We found several individual examples of cells that encoded threat-associated information in the rat exposure assay. These examples included cells that were more active near the rat (Fig. 3A), and during higher periods of mouse or rat speed (Fig. 3B,C). We also found examples of cells with the opposite pattern of responses, that were more active far from the rat and during periods of low mouse or rat speed. Strikingly, we found cells that encoded the angular offset between the mouse's head direction and the rat's body (Fig. 3D). We also found cells that were activated during various defensive behaviors, such as escape, freezing, and risk-assessment stretch-attend postures (Fig. 3E).
dPAG cells encode multiple threat-related features and defensive behaviors in the Rat exposure assay. A, Heat maps represent the activity of representative dPAG cells that were more active either near (top) or far (bottom) from the rat. B, Scatter plots represent example cells in which neural activity was positively (left) or negatively (right) correlated to the mouse's speed. C, Example cells displaying positive (left) or negative (right) correlation with rat speed. D, Polar plots for cells that were active at different mouse head direction angle offsets relative to rat body. Red represents z scores corresponding to different-sized radius in the polar plot. Warmer colors also represent higher z scores. E, Color map represents activity of all recorded cells aligned to time of highest activity relative to behavioral onset. Each row corresponds to an individual dPAG cell. F, Scheme represents how GLMs were built to characterize dPAG activity. Two separate GLMs were performed. The first used continuous kinematic variables (mouse and rat speed, angular offset between threat and mouse, etc.), while the second used binary behavioral variables (approach, stretch, escape, and freeze). Both categories of variables were convolved with kernels producing predictors. These predictors were multiplied by different weights and linearly combined to optimize the reconstruction of the activity of each dPAG neural activity trace using either only kinematic variables or only behavioral variables. G, Example traces represent actual neural activity (blue) and GLM-reconstructed predicted activity for two example traces using GLMs with only kinematic (pink) or behavioral (green) variables. H, Bars represent the mean relative contribution of kinematic and behavioral variables across neurons in the Rats 1 and 2 assays (Rat 1: n = 715; Rat 2: n = 625). I, Color map represents correlation of GLM weights for kinematic and behavioral variables across neurons in the Rat 1 assay. Bold values represent significantly strong correlation of weights between two related variables (n = 715, bold: bootstrapping with random GLM weights, p < 0.001). J, Bars represent the percentage of neurons in the Rat assay that were significantly positively modulated (white), negatively modulated (black), or not modulated (gray) for the analyzed variable based on their GLM weights (n = 715, bold: p < 0.05, bootstrapping with random GLM weights). K, GLM weights for the first and second halves of Rat 1 show significant correlation (n = 7 mice, **p < 0.01, bootstrapping with random GLM weights).
To quantify how strongly these variables were represented in neural activity, we constructed a GLM to quantify how dPAG cells encoded these features (Fig. 3F) (Engelhard et al., 2019). A GLM can quantify how different variables linearly relate to the calcium traces from each cell. We trained two separate GLMs. One predicted dPAG cell activity from continuous kinematic variables, including rat and mouse velocity, distance and angular offset between mouse head direction and rat body position. Another GLM predicted dPAG cell activity during defensive behaviors (i.e., threat approach, stretch-attend posture, freezing, and escape). Importantly, these two GLMs were done in nonoverlapping time points. The kinematic GLM was done in time points that do not have any scored behaviors. The behavioral GLM was done only in time points with scored behaviors. This approach ensures that the kinematic GLM weights are not being influenced by ongoing defensive behaviors because time points with behaviors were not used for this analysis. After completing the kinematic GLM, we used the identified weight for distance to threat and regressed out this variable from the calcium traces. Then, only the time points containing scored behaviors were used for the behavior GLM. We regressed out distance to threat before performing the behavioral GLM, ensuring that the GLM weights for behaviors are not simply reflecting distance from the rat. For example, freezing tends to happen far from the rat (Fig. 1D); thus, cells that simply are active far from the rat may be erroneously identified as cells that are positively modulated by freezing. However, as the influence of distance to threat was regressed out from the traces before performing the behavioral GLM, this problem was avoided. The approach used here ensures that behavioral weights are not reflecting encoding of distance to threat and vice-versa.
Using these GLMs, it was possible to reconstruct the neural activity of single-cell dPAG calcium traces with high accuracy (Fig. 3G). We found that each of the GLM variables was predictive of dPAG activity and had comparable relative contributions to dPAG activity (Fig. 3H; see Materials and Methods; Rat 1: n = 715 cells; Rat 2: n = 625 cells). A positive GLM weight for a particular variable indicates that a cell displays higher GCaMP6s than its own average Ca2+-dependent fluorescence when that variable increases, whereas negative weights indicate lower fluorescence. GLM weights were generally not strongly correlated across variables (Fig. 3I; n = 715 cells, p < 0.001, bootstrapping with random GLM weights). We observed a negative correlation between speed and freezing (−.043), which is an expected result as freezing occurs during the lowest mouse speed epochs. Importantly, distance to rat weights was not strongly correlated with behavioral weights, indicating that cells encoded these features independently. We determined whether cells had significant positive or negative weights for each GLM input variable through a bootstrap of GLM coefficients (at p < 0.05, see Materials and Methods). There were large populations of cells with either positive or negative weights for each kinematic variable and defensive behavior (Fig. 3J; n = 715 cells). Finally, GLM weights were relatively stable over the course of an experiment, and GLMs fit separately on the first and second halves of an experiment displayed significant positive correlation (Fig. 3K; n = 7 mice, p < 0.01, bootstrapping with random GLM weights).
Similarly to the Rat assay, in the shock grid assay, we also observed bidirectional modulation of activity in single dPAG cells that responded to distance from the threat (Fig. 4A) and to mouse speed (Fig. 4B). We also identified cells that encoded angular offset between mouse head direction and the shock grid (Fig. 4C). Analogously to the results observed in the Rat assay, a GLM similar to the one used for the Rat assay was able to reconstruct dPAG neural activity in the shock grid test as well (Fig. 4D), and each GLM variable had comparable relative contributions to dPAG activity (Fig. 4E; Fear Acquisition: n = 737 cells; Fear Retrieval: n = 621 cells). GLM weights were generally not strongly correlated across variables as well (Fig. 4F; n = 737 cells, p < 0.001, bootstrapping with random GLM weights). A large fraction of cells had significant positive and negative weights for each of the variables (Fig. 4G; n = 737 cells), similarly to the GLM in the Rat assay. GLM weights were also stable across the first and second halves of a recording session (Fig. 4H; n = 6 mice, p < 0.01, bootstrapping with random GLM weights). Overall, in both assays, GLM variables were not highly correlated in a consistent manner, and generally weights were only weakly correlated with each other (Figs. 3I, 4F). Nevertheless, in both assays, approach weights were positively correlated with mouse speed. This finding indicates that approach cells may be encoding features related to higher exploratory drive. These data show that dPAG cells encoded a wide variety of defensive behaviors and threat-related measures during exposure to both a live predator and a shock grid, highlighting the richness and complexity of dPAG neural activity.
dPAG cells encode multiple threat-related features and defensive behaviors in the shock grid assay. A, Heat maps represent the activity of representative dPAG cells that were more active either near (top) or far (bottom) from the shock grid. B, Scatter plots represent example cells in which neural activity was positively (left) or negatively (right) correlated to the mouse's speed. C, Polar plots for cells that were active at different mouse head direction angle offsets relative to shock grid. D, Example traces represent actual neural activity (blue) and GLM-reconstructed predicted activity for two example traces using GLMs with only kinematic (pink) or behavioral (green) variables. E, Bars represent the mean of relative contribution of variables across neurons in the Fear Acquisition and the Fear Retrieval assays (Fear Acquisition: n = 737; Fear Retrieval: n = 621). F, Color map represents correlation of GLM weights for kinematic and behavioral variables across neurons in the Fear Acquisition assay. Bold values represent significantly strong correlation of weights between two related variables (n = 737, bold: p < 0.001, bootstrapping with random GLM weights). G, Bars represent the percentage of neurons in the Fear Acquisition assay that were significantly positively modulated (white), negatively modulated (black), or not modulated for the analyzed variable based on their GLM weights (n = 737, p < 0.05, bootstrapping with random GLM weights). H, GLM weights for the first and second halves of Fear Retrieval show significant correlation (n = 6 mice, **p < 0.01, bootstrapping with random GLM weights).
dPAG cells encode distance to threat and defensive behaviors independently
The results above indicate that dPAG cells prominently represent both distance to threat and defensive behaviors. Indeed, the behavioral GLM identified large fractions of cells that were modulated by each of the defensive behaviors. These cells will be referred to as approach cells, escape cells, etc, hereafter.
We next evaluated whether behaviors and distance to threat were encoded independently. To do so, we analyzed GCaMP6s traces that were binned in a narrow range of distance to threat. We then compared whether, within that same bin, behaviorally modulated cells showed higher activity during behavioral time points compared with time points when the behavior did not occur (Fig. 5A,B; nonbehavioral events). For example, for escapes, we used a bin corresponding to 20 cm immediately adjacent to the rat zone, as that is the region that contained most of the escapes. We then measured the average activity of escape cells in that area during escape time points and during nonescape time points (behavioral events and nonbehavioral events, respectively; Fig. 5A,B). Even when restricting the analysis to a narrow bin of distance to threat, escape cells (identified by significant GLM weights; Figs. 3J, 4G) showed higher activation compared with nonescape time points. This result was observed for cells modulated by any of the four behaviors (Fig. 5B; cells that are significantly modulated for approach, stretch, escape, and freeze: Rat 1 assay, n = 506, 476, 448, 557 cells, p = 4 × 10−19, 2 × 10−9, 3 × 10−11, 10 × −14; Fear Acquisition assay, n = 525, 556, 423, 524 cells, p = 4 × 10−69, 5 × 10−69, 5 × 10−52, 10−64, Wilcoxon signed-rank test). Importantly, the behavior GLM was only used in this analysis to identify which cells were behavior modulated. The calcium traces used in Figure 5 were the raw fluorescence traces, and we did not regress out distance to threat from these traces. This result shows that distance to threat cannot explain the behavioral modulation seen in dPAG cells.
dPAG cells encode defensive behaviors and distance to threat independently. A, Scheme represents method to compare activation of a GLM-classified escape cell during escape and nonescape events occurring at a narrow window of 20 cm from the rat. Nonescape events corresponded to time points at which no scored behaviors were observed. This escape cell shows higher activation during escape than nonescape, although both events occur at a similar distance to the rat. The same analysis was done for GLM cells with significant weights for other behaviors, comparing their activation during behavior and nonbehavior events at a fixed distance to the rat. B, Bars represent absolute values of z score dF/F for neurons that were significantly modulated for defensive behaviors during behavioral events (black) or nonbehavioral events (white). These bars were computed over a restricted position window, so that all behaviors occurred at a similar distance to the threat. For approach and freeze, this position window was the 20 cm window closest to the left “safe” wall, while for stretch and escape events, it was the 20 cm window nearest the threat. Nonbehavioral events refer to a 4-s-long period within the position window during which the mouse exhibited no defensive behavior (see Materials and Methods). For each type of event, absolute values of mean z score dF/F across time were taken during these events for neurons that were significantly modulated for defensive behaviors. Even at the same distance to threat, cells encoding a particular behavior were more active during that behavior compared with other time points, showing that encoding of behaviors is not a direct consequence of encoding of distance to threat (Rat 1 assay, cells that are significantly modulated for approach, stretch, escape, and freeze, n = 506, 476, 448, 557; Fear Acquisition assay, cells that are significantly modulated for approach, stretch, escape, and freeze, n = 525, 556, 423, 524; Rat 1 assay, ***p = 4 × 10−19, 2 × 10−9, 3 × 10−11, 10−14; Fear Acquisition assay, ***p = 4 × 10−69, 5 × 10−69, 5 × 10−52, 10−64, Wilcoxon signed-rank test). C, Top, Scheme represents that a representative escape-modulated cell displays a similar magnitude of activation by escape regardless of the distance to the rat at which the escape is initiated. r1 is the correlation between the two cell activity traces f1 and f2, and a1 is the absolute value of difference in the distance in cm between the location of onset and the rat of these two escapes. These two metrics were calculated for all pairs of escapes in the session. Bottom, Scatterplot represents the pairwise correlations between the representative cell's activity during all escapes plotted against the difference in the distance between the location of onset and the rat of these two escapes. The correlation among these points is zero, indicating that activations are not more correlated between escapes that occur near each other compared with escapes that occur far apart. The slope for all escape-modulated cells was calculated, and the average of all scatter plot slopes for all escape-modulated cells in the Rat assay is shown in D. The same analysis was repeated for cells significantly modulated by the other three behaviors. D, Violin plots represent that the slopes of regressions for absolute values of difference in distance to threat and correlations between dF/F are closed to zero in both Rat 1 and Fear Acquisition assays, showing that the activation during behaviors is not modulated by the distance to threat at which the behavior onset occurs (n is the same as in B). E, Cells were classified as being positively, negatively, or not modulated by distance to threat. Bar graphs represent the percentage of cells for each of these groups that is modulated by each of the behaviors. Modulation by distance to threat did not affect the % of cells that encodes each behavior. A cell could be modulated for multiple behaviors, which is why the percentage bars do not add to 100%.
We next investigated whether dPAG behavioral responses were modulated by distance to threat. If this modulation exists, then, for example, an escape cell would show similar responses to escapes that occurred near each other, and would display more different escape responses to escapes that occur far from each other (Fig. 5C). To measure similarity of calcium traces for a given escape cell, we calculated pairwise correlations of all the pairs of escapes. Plotting these correlations against the distances between these escape onsets shows a slope of zero. Thus, escape pairs that occurred near each other did not have more strongly correlated neural activity than escape pairs that occur far for the representative cell shown (Fig. 5C). On average, for all escape cells, the observed slope was also centered at zero, and similar results were also observed for all other behaviors (Fig. 5D; n is the same as in Fig. 5B). These data indicate that behavioral responses were not consistently modulated by the location in which the behavior onset occurs.
Last, cells with positive or negative GLM weights for distance to threat showed similar encoding of various behaviors (Fig. 5E; Rat 1, n = 7 mice; Fear Acquisition, n = 8 mice). For example, in the Rat assay, ∼75% of the cells were significantly modulated by approach, regardless of whether the cells were positively, negatively, or not modulated by distance to threat (Fig. 5E, left). These data also show that the same cells encode behaviors and distance to threat.
Together, the results in Figure 5 show that defensive behaviors and distance to threat are being independently encoded by dPAG cells.
Rat and shock grid assays induce shared dPAG activity patterns during defensive behaviors
We then investigated whether the dPAG represents behavioral and kinematic variables similarly across innate (rat) and conditioned (shock grid) threat assays. We first calculated the overlap of positively and negatively modulated cells for all the variables across both assays. The resulting Venn diagrams (Fig. 6A) show that several behavior-modulated cells maintain their GLM classification across assays (e.g., of 138 cells showing significant positive modulation to approach to threat in the Rat assay, 81 maintain that classification in the shock grid assay).
GLM weights are more similar across threatening assays than across control assays. A, Venn diagrams represent the overlap of coregistered neurons that are positively (top row) or negatively (bottom row) modulated for variables shown in columns in Rat 1 and Fear Acquisition assays. For example, of 103 cells showing positive modulation by distance to threat in Rat 1, 51 also displayed positive modulation by this variable in the Fear Acquisition assay. B, Illustration of the cosine similarity metric. Top, Three example points, relating GLM weights for distance to threat and escape across Shock (Fear Aquisition), Rat, and Preshock assays. The angle between the Shock and Rat vector is smaller, producing a larger cosine, while the angle between the Rat and Preshock is larger, resulting in a smaller cosine similarity measure. This method was used for each cell coregistered across two assays to calculate the similarity between the GLM weights for all variables in the two assays. Only two variables were used in the scheme for easier visualization but was done with all variables in the data plotted in C, D. C, Bars represent cosine similarities between GLM coefficients in Fear Acquisition assay and those in all the other assays. dPAG ensembles in Fear Acquisition experiments have stronger similarities with other threat assays (Rat1, Rat2, Fear Retrieval) than control assays (Toy Rat, Preshock) (cells coregistered between Fear Acquisition and Rat 1, Rat 2, Fear Retrieval, Toy Rat, Preshock, n = 300, 262, 262, 275, 388; U = 5.8, ***p = 6 × 10−9, Wilcoxon rank-sum test). D, Bars represent cosine similarities between GLM coefficients in threat-related assays and control assays. dPAG ensembles have significantly stronger similarities between different threat assays than comparisons between threat and control assays (U = 10.88, ***p = 10−27, Wilcoxon rank-sum test), and between different control assays (U = 5.79, ***p = 7 × 10−19, Wilcoxon rank-sum test; threat vs threat, n = 1404; threat vs control, n = 2463; control vs control, n = 357).
To accurately quantify the similarity across assays of GLM weights for all 7 variables, we used the cosine similarity metric. This measure provides a quantification of the similarity of all GLM weights across two assays for each cell. A representative example in Figure 6B illustrates the use of this measure. In this example cell, the GLM weights for distance to threat and escape are plotted for three assays. These values are more similar between the Rat and shock grid assay than between the Rat and the Preshock assay. Consequently, the cosine of the angle between Rat and shock is larger than that for the angle between Rat and Preshock. In this example, only 2 variables were used for simplicity, but in the actual analysis we used all the variables, as this metric can calculate the cosine of the angle between two vectors in multidimensional space, not just in two dimensions. As an example, we show cosine similarity between Fear Acquisition and all other assays. Cosine similarity is higher between Fear Acquisition and threat assays than control assays (Fig. 6C; cells coregistered between Fear Acquisition and Rat 1, Rat 2, Fear Retrieval, Toy Rat, Preshock, n = 300, 262, 262, 275, 388 cells; U = 5.8, p = 6 × 10−9, Wilcoxon rank-sum test). Similarly, averaging across all assay pairs, cosine similarity is highest between threat-threat assay pairs than threat-control (U = 10.88, p = 10−27, Wilcoxon rank-sum test) or control-control (U = 5.79, p = 7 × 10−19, Wilcoxon rank-sum test) pairs (Fig. 6D; threat-threat, n = 1404 cells; threat-control, n = 2463 cells; control-control, n = 357 cells). This result shows that across the population of dPAG cells, GLM values were highly conserved between threat assays, but not other conditions.
dPAG cell activity patterns segregate into shared clusters across Rat and shock grid assays
We next studied whether the dPAG cells formed functional clusters based on their response patterns to defensive behaviors. We clustered cells based on their behavioral GLM weights using a GMM (see Materials and Methods) (Engelhard et al., 2019). The GMM enables us to compute the probability that a cell belongs to a given cluster. We subsequently classify each cell to the cluster with the highest probability. Cells that have similar GLM weights across the four behaviors are more likely to be classified in the same cluster (see Materials and Methods). This is illustrated in Figure 7A where we plot the approach and escape GLM weights for many cells (each cell is one point). The GMM models 2 clusters (purple and green) via a Gaussian distribution, enabling us to compute the probability a cell belongs to the purple or green cluster. This is equivalent to defining a decision boundary (Fig. 7A, right), where a cell above (below) the decision boundary is in the purple (green) cluster. We used the BIC scores to determine the optimal number of clusters. Significant drop in BIC scores were found when changing the number of clusters from 1-2, 2-3 and from 3-4 (Rat assay, n = 6 mice, p = 0.002, 0.013, 0.018; Fear Acquisition, n = 8 mice, p = 0.001, 0.009, 0.012; Fear Retrieval, n = 5 mice, p = 0.001, 0.020, 0.028, Wilcoxon signed-rank test). However, no significant changes were found when increasing the number of clusters from 4-5 (Rat assay, p = 0.128; Fear Acquisition, p = 0.123; Fear Retrieval, p = 0.249, Wilcoxon signed-rank test). Thus, we used four clusters (Fig. 7A,B; see Materials and Methods). On all assays, cells could be clustered based on their responses to the scored behaviors (approach, stretch, escape, and freeze) (Fig. 7C; Rat, Fear Acquisition, Fear Retrieval n = 640, 737, 561 cells). This clustering achieved significantly higher log-likelihood than on random samples (p < 0.001, bootstrapping with random input variables or cells, all assays). The analysis also shows that cells were most likely to strongly be activated during only one of the four scored defensive behaviors (Fig. 7D,F,H; approach, stretch, escape, freeze clusters: Rat assay, n = 148, 165, 159, 168 cells; Fear Acquisition assay, n = 210, 204, 158, 165 cells; Fear Retrieval assay, n = 141, 161, 81, 178 cells). For example, cells from the freeze cluster in the Rat assay only had large positive GLM weights for freeze (Fig. 7D), and these cells were only strongly activated during freeze, but not other behaviors (Fig. 7E). Across all clusters, we found that cells were also activated for distance to threat. Distance to threat GLM weights was significantly higher than behavioral GLM weights in only the approach and escape clusters across all assays (p = 0.013, Wilcoxon signed-rank test, not shown).
Clustering of dPAG cells according to their activity during defensive behaviors. A, Illustration of GMM. Left, The GLM weights for approach and escape behaviors are plotted for a fabricated set of example cells. The mean and covariance of the Gaussian-distributed cluster participation probabilities (purple and green gradients) are determined by the iterative Expectation Maximization algorithm. Right, These probabilities are then used to define a decision boundary (red dotted line) and classify each cell to its most probable cluster. B, Traces represent BIC as a function of numbers of clusters in a GMM based on GLM coefficients in Rat 1 (left), Fear Acquisition (middle), and Fear Retrieval (right) assays. Increasing the number of clusters until 4 results in a significant decrease in BIC and a significant improvement in clustering performance, across all three assays (Rat assay, *p = 0.002, 0.013, 0.018; Fear Acquisition, *p = 0.001, 0.009, 0.012; Fear Retrieval, *p = 0.001, 0.020, 0.028, Wilcoxon signed-rank test) C, dPAG cells were clustered using a GMM based on their GLM weights for approach, stretch, escape, and freeze. Cells are grouped by the probability of belonging to a particular cluster. Right, Colored vertical lines indicate cluster identity. For example, the first cluster for the Rat assay corresponds to green, which is a cluster activated during approach (see color legend, top right corner). All assays showed one cluster for each defensive behavior (n = 640, 737, and 561 cells for Rat, Fear Acquisition, and Fear Retrieval, respectively). D, Bar plot represents the average z-scored GLM weights across behaviors for each of the four clusters found in the Rat assay. E, Behavior triggered average z-scored dF/F traces for each of the clusters found in the Rat assay. Data are mean ± SEM (Rat assay, approach, stretch, escape, freeze clusters: n = 148, 165, 159, 168). F-I, Characterization of clusters from shock grid assay. F, Bar plot represents the average z-scored GLM weights across behaviors for each of the four clusters found in the shock grid Fear Acquisition assay. G, Behavior triggered average z-scored dF/F traces for each of the clusters found in the shock grid Fear Acquisition assay. Data are mean ± SEM (Fear Acquisition assay, approach, stretch, escape, freeze clusters: n = 210, 204, 158, 165). H, I, Same as in D, E, but for shock grid Fear Retrieval day (Fear Retrieval assay, approach, stretch, escape, freeze clusters: n = 141, 161, 81, 178). ns, not significant.
We next used the cluster definitions from one assay and applied them to coregistered cells in other assays. The numbers represent the fraction of cells in a cluster that also were assigned the same cluster in other assays (Fig. 8). For example, in the matrix corresponding to a clustering model trained in Rat 1, but applied to Rat 2 (Fig. 8, top left matrix), 0.37 in the lower right corner indicates that, among coregistered cells in those two assays, 51% of the cells that belonged to cluster 4 in Rat 1 were also assigned to cluster 4 in the Rat 2 assay (using the cluster four definitions from the Rat 1 assay). Only numbers that are significantly above chance are shown. All numbers in the diagonals in all assay pairs were significant, indicating that dPAG cells can be consistently functionally clustered across conditioned and unconditioned threat modalities based on their responses to ongoing defensive behaviors (first column, trained in Rat 1, n = 97, 285, 160 cells; second column, trained in Rat 2, n = 97, 170, 124 cells; third column, trained in Fear Acquisition, n = 285, 170, 237 cells; fourth column, trained in Fear Retrieval, n = 160, 124, 237 cells, p < 0.05, bootstrapping with random GLM weights).
dPAG ensembles share functional GMM cluster assignments across assays. Cluster definitions for each assay were created using a GMM approach. These definitions were then applied to coregistered cells recorded in the other assays. For the first plot, cluster definitions were created based on GMM clustering in the Rat 1 assay. These same cells were then assigned to clusters based on Rat 1 cluster definitions but using only behavioral variables and neural activity from the Rat 2 assay (this is the second exposure to the rat). Numbers indicate the percentage of cells from the Rat 1 clusters that are assigned to the same cluster using Rat 2 recordings. The GMM was thus trained on Rat 1, and the same cluster definitions were applied to data from the Rat 2 recording. For example, the 0.40 in the first confusion matrix indicates that 40% of the cells assigned to Cluster 1 in Rat 1 are also assigned to Cluster 1 in Rat 2. Warmer colors represent higher percentages of cells. Values are highest (and thus have the warmest colors) in diagonal elements, indicating that, in general, cells had a high probability of being assigned to the same functional GMM cluster in both Rats 1 and 2. Similarly, applying the same Rat 1 cluster boundaries on Fear Acquisition and Fear Retrieval data produced the confusion matrices in the second and third rows of the first column. The second, third, and fourth columns represent GMM cluster boundaries across assays for clusters defined, respectively, using training data from Rat 2, Fear Acquisition, and Fear Retrieval: first column, trained in Rat 1, predict in Rat 2 (n = 97), Fear Acquisition (n = 285), Fear Retrieval (n = 160); second column, trained in Rat 2, predict in Rat 1 (n = 97), Fear Acquisition (n = 170), Fear Retrieval (n = 124); third column, trained in Fear Acquisition, predict in Rat 1 (n = 285), Rat 2 (n = 170), Fear Retrieval (n = 237); fourth column, trained in Fear Retrieval, predict in Rat 1 (n = 160), Rat 2 (n = 124), Fear Acquisition (n = 237) (p < 0.05, bold: bootstrapping with random GLM weights).
dPAG activity can be used to decode ongoing defensive behaviors across assays
The aforementioned results indicate that dPAG cells have conserved patterns of neural activity to represent defensive behaviors during exposure to a live predator and to a shock grid. These data suggest that ongoing defensive behaviors can be predicted from dPAG neural data. To test this hypothesis, we used multinomial logistic regression to classify behaviors into approach, escape, freezing, and stretch and used fivefold cross-validation for both training and testing dataset. To remove the effect of distance to threat on neural activity, before training and testing, we subtracted the linear prediction of cell activity from distance to threat. Using this approach, we show that all ongoing monitored behaviors could be robustly predicted within each assay (Fig. 9, top row). We then tested whether the dPAG has a shared representation of defensive behaviors across assays. Because each assay may have an independent and shared representation of the defensive behavior, data from both assays should be used to isolate a shared representation (Elsayed et al., 2016; Jiang et al., 2020). We thus trained our model on combined data from pairs of Rat 1, shock grid Fear Acquisition, and Fear Retrieval to identify if defensive behaviors could be decoded consistently across assays, indicating a shared representation. Behavioral prediction accuracy was then tested on held out nonoverlapping data from these same two assays averaged over fivefold cross-validation. Across all pairs of assays, prediction accuracy was significantly above chance (0.25, or 25%, p < 0.05), showing that these behaviors are consistently encoded across assays (Fig. 9, bottom row). Significant values are shown in yellow. Consistent above chance decoding both within and across assays indicates conserved neural population representations of behavior during differing sessions. That freezes are mistaken as stretches during Fear Acquisition and Fear Retrieval sessions, but not when combining training with another session suggests that logistic regression trained on pairs of assays finds a shared axis which better generalizes neural representations of behavior, even if diagonal true positive rates are lower. This may also be evidence of nonparallel shared and context-specific axes representing threat-related behaviors. Together, these data show that dPAG cells use conserved activity patterns to represent defensive behaviors elicited by innate and learned threats.
Defensive behaviors can be discriminated from dPAG neural activity. Behaviors can be robustly classified within a session (top row, Rat: n = 713 cells total, Fear Acquisition: n = 747 cells, Fear Retrieval: n = 514 cells) and across multiple assays (bottom row, Rat and Fear Acquisition: 299 cells total, Rat and Fear Retrieval: 198 cells, Fear Acquisition and Fear Retrieval: 260 cells) from cell activity with distance to threat regressed out. Classification was performed using fivefold cross-validated logistic regression. Confusion matrices show the averaged testing accuracy of behavior discrimination across all folds (yellow numbers: p < 0.05, randomized training labels).
dPAG cell activity is sufficient to predict future escape behavior in Rat assay
We next evaluated the predictive power of dPAG cell activity with respect to escape from threats (Fig. 10A). To determine whether dPAG activity patterns could predict escape in threat exposure contexts, we trained a neural network to identify escape occurrences in the present and up to 9 s in the future (example in Fig. 10B). The model was able to predict escape behavior occurrence significantly better than chance level at most 5 s in the future in the Rat assay, 3 s in the Fear Acquisition session, and 4 s in the Fear Retrieval session (Fig. 10C,D; p = 0.025, p = 0.021, and p 0.043, respectively, Wilcoxon signed-rank test; more statistical details can be found in the Fig. 10 legend). Our results thus show that the dPAG activity can predict escape occurrence significantly better than chance level for both conditioned and unconditioned threat stimuli. Importantly, this prediction was not possible in the Toy Rat assay or in the Preshock session for any time lags ≥1 s (Fig. 10C,D, left panels, p = 0.063 and p = 0.337, respectively, Wilcoxon signed-rank test), demonstrating that dPAG activity only predicts future escape from threats, rather than movement away from objects in general. Although escape could be predicted, using this same approach, dPAG cells did not predict any other behavior in the future (data not shown), in agreement with prior work, which most strongly implicates the dPAG with escape (Evans et al., 2018; Tovote et al., 2016). These data suggest a potential role for the dPAG in the moments preceding the escape response.
Escape prediction through a neural network model using dPAG activity. A, The activity for each neuron was calculated using z-scored activity between −5 and 5 s with respect to onset of escape in the Rat assay. B, Example segment represents predicted escapes (pink) and true observed escape (green triangles) for example prediction (pink) and true behavior (green triangles), testing data. C, Line plot indicates weighted accuracy for predicting escape at various lags of testing data (n = 7 mice, from left to right n = 826, n = 878). Escape could be predicted using dPAG activity significantly higher than chance up to 5 s in the future for the Rat assay (from 0 to 9 s: p = 0.018, p = 0.028, p = 0.018, p = 0.018, p = 0.091, p = 0.025, p = 0.180, p = 0.180, p = 0.180, p = 0.180, Wilcoxon signed-rank test), and only at escape onset for the Toy Rat assay (p = 0.018 for 0 s, p = 0.063 for 1 s, and p > 0.05 for all other points, Wilcoxon signed-rank test). The weighted accuracy was calculated using equal weights for samples that show escape and samples that do not show escape (50% each); thus, chance levels of this measurement are 50%. Data were separated in timewise contiguous pieces, which were used as either training or testing. D, Same as in C, but for Preshock (n = 6 mice, 633 cells, p = 0.004 for 0 seconds, p = 0.337 for 1 s, and p > 0.05 for all other points, Wilcoxon signed-rank test), Fear Acquisition (n = 6 mice, 627 cells, from 0 to 9 seconds, p = 0.043, p = 0.080, p = 0.043, p = 0.043, p = 0.043, p = 0.062, p = 0.625, p = 0.812, p = 1, p = 0.812, Wilcoxon signed-rank test) and Retrieval (n = 5 mice, 569 cells, from 0 to 9 seconds p = 0.021, p = 0.021, p = 0.021, p = 0.021, p = 0.248, p = 0.125, p = 0.625, p = 0.125, p = 0.625, p = 0.875, Wilcoxon signed-rank test).
Shared neural representation of distance to threat across assays
The results discussed above show that dPAG cells exhibit consistent activation patterns for multiple behaviors during exposure to a live predator or a shock grid (Fig. 9). We next investigated whether dPAG cells also have a conserved representation of distance to threat across these assays. To do so, we calculated the correlation between neural activity and distance to threat in the Rat 1 assay for all cells coregistered between Rat 1 and shock grid Fear Acquisition assays. Lower distance to threat values indicate locations nearer the rat. Cells with negative and positive correlations to distance to threat were, respectively, classified as rat-activated (red) and rat-inhibited cells (blue) (Fig. 11A). Other cells were classified as neither (i.e., rat-insensitive) (gray). We then plotted the neural activity of these three cell types while mice left the safe side of the environment and started approaching either the rat or the shock grid on Fear Acquisition day. As rat-activated cells were defined as cells that show a positive correlation between activity and proximity to the rat, it was expected that these cells would show increased activity when mice left the safe zone of the environment and started approaching the rat. Interestingly, these same cells also showed increased activity when mice approached the shock grid during Fear Acquisition. Conversely, rat-inhibited cells displayed decreases in activity when mice approached the rat (Fig. 11B,C; n rat-activated = 110 cells, n rat-inhibited = 52 cells, n neither = 137 cells; Rat 1: U = 7.30, p = 3 × 10−13, Fear Acquisition: U = 5.12, p = 3 × 10−7, Wilcoxon rank sum test). We then repeated this same analysis but defined threat sensitivity using data from cells coregistered between Rat 1 and Fear Retrieval day. Cells that were activated or inhibited by proximity to the rat also showed the same pattern of activation relative to distance to the shock grid on Fear Retrieval day (Fig. 11D–F; n rat-activated = 43 cells, n rat-inhibited = 42 cells, n neither = 113 cells; Rat 1: U = 6.00, p = 2 × 10−9, Fear Retrieval: U = 3.82, p = 1.3 × 10−4, Wilcoxon rank sum test). Thus, dPAG cells displayed a shared encoding of distance to threat which is generalizable across rat and shock grid exposures.
dPAG ensembles encode distance to innate and conditioned threats using shared activation patterns. A, Histogram represents correlation between neural activity and distance to threat in the Rat assay on the first day (higher x values correspond to locations more near the rat). Blue represents cells that display positive correlation (r > 0.10) with distance to rat during rat exposure (n = 52). Red represents cells that are more active closer to the rat (r < −0.15) (n = 110). Gray represent remaining cells (n = 137). Cells coregistered between Rat 1 and Fear Acquisition assay are shown. B, Traces represent average activity for dPAG cells relative to when animals leave the safe side of the environment and start approaching the threat. Cells defined as rat-activated (red) cells in A show increases in activity when mice approach both the rat or the shock grid during Fear Acquisition. Conversely, rat-inhibited cells defined in the Rat Assay (A) are also inhibited during approach to the rat and shock grid. C, Quantification of neural activity changes from the traces in B, difference in mean activity of 0-2.5 s after leaving the safe side and 0-2.5 s before leaving the safe side. Rat-activated cells displayed higher activity while approaching both rat or shock grid, and rat-inhibited cells were inhibited during approach to threat in both assays (n same as in A, Rat 1: U = 7.30, p = 3 × 10−13, Fear Acquisition: U = 5.12, p = 3 × 10−7, Wilcoxon rank-sum test). D–F, Same as in A–C, but using shock grid Fear Retrieval session instead of Fear Acquisition session. Cells coregistered between Rat 1 and Fear Retrieval are shown (rat-activated n = 43, rat-inhibited n = 42, neither n = 113, Rat 1: U = 6.00, p = 2 × 10−9, Fear Retrieval: U = 3.82, p = 1.3 × 10−4, Wilcoxon rank-sum test). G, Correlation of example cell (cell 153) activity with distance to threat during Rat 1 (top) and Rat 2 assays (bottom). Lower distance to threat corresponds to locations more near the rat. Each point is an individual time point showing activity of cell 153 in these two assays. H, Scatterplot represents correlation between individual cell correlations of cell activity with distance to threat during rat exposure on days 1 and 2. Blue represents cell 153 from G. Each point corresponds to one cell. I, Correlation of correlations between cell activity with distance to threat for each pair of recording sessions. For example, the 0.58 correlation between Rats 1 and 2 shown in H is depicted at the top left corner of the matrix in I. J, Bar graph represents significant difference between correlation of correlations between cell activity with x position for threat-threat versus threat-control pairings (mean ± SEM; U = 2.71, p = 0.007, Wilcoxon rank-sum test). Color-coded purple and pink boxes, shown in I, represent which values were used to calculate threat versus threat and threat versus control averages.
The correlations between neural activity and distance to rat (Fig. 11A,D) indicate that the activity of dPAG cells may encode distance to threat. Figure 11G shows a representative example cell with a significant negative correlation between distance to threat and neural activity in both Rat 1 (top) and Rat 2 (bottom) assays. This representative example suggests that cells may display similar amounts of correlation between neural activity and distance to threat across Rats 1 and 2 assays. Indeed, we show that correlations of activity and distance to threat in these 2 assays tended to be similar (see plot showing this “correlation of correlations” in Fig. 11H for all cells recorded in both Rats 1 and 2 assays). Rats 1 and 2 were the most similar assays, and thus displayed the highest value of correlation of correlations (r = 0.56, Fig. 11I, top left corner). Similar calculations across all assay pairs suggest that pairs of assays consisting of two threatening assays (e.g., Rat 1 and Fear Acquisition) have higher correlations of activity-distance correlations than pairs with one control and one threatening assay (e.g., Rat 2 and Toy Rat). Values corresponding to threat-threat pairs and threat-control assay pairs are, respectively, outlined in purple and pink in Figure 11I. Threat-threat assay pairs (purple) showed significantly higher correlation of correlations than threat-control (pink) values (Fig. 11J; U = 2.71, p = 0.007, Wilcoxon rank sum test). This result shows that cells tended to conserve their relationship of neural activity and distance to threat across threatening assays more than between threatening and control assays. For example, a cell that was more active near the rat will also tend to be more active near the shock grid post-shock, but not near the Toy Rat or the shock grid during the Preshock session.
Importantly, repeating the analysis from Figure 11 using control assays yielded negative results (Fig. 12). We first defined rat-excited and rat-inhibited cells for all cells that were tracked across Rat 1 and Preshock sessions. Rat-excited and -inhibited cells did not display any consistent changes in activity during proximity to the shock grid in the Preshock session (Fig. 12A–C; n rat-activated = 112 cells, n rat-inhibited = 55 cells, n neither = 173 cells; Rat 1: U = 7.23, p = 5 × 10−13, Preshock: U = 1.30, p = 0.19, Wilcoxon rank sum test). Of course, these cells showed changes in activity during approach to the rat (Fig. 12B) because they were originally defined as ensembles with positive or negative correlation between distance to rat and neural activity (Fig. 12A). Similarly, cells that were excited or inhibited with proximity to the control Toy Rat did not change their activity during approach to the shock grid in the Preshock session (Fig. 12D–F; n toy rat-activated-activated = 64 cells, n toy rat-inhibited = 81 cells, n neither = 253 cells; Toy Rat: U = 6.35, p = 2 × 10−10, Preshock: U = −0.29, p = 0.77, Wilcoxon rank sum test). Together, these data indicate that the dPAG uses a shared representation to encode different threats (Fig. 11), but not safe control objects (Fig. 12).
dPAG cells do not have a shared encoding of distance across threatening and neutral stimuli. A, Histogram represents correlation between neural activity and distance to rat in Rat 1 for all cells coregistered between Rat 1 and shock grid Preshock habituation session. Higher x values correspond to locations nearer the rat. Blue represents cells that display positive correlation (r > 0.10) with x position during rat exposure (n = 55). Red represents cells that are more active closer to the rat (r < −0.15) (n = 112). Gray represents remaining cells (n = 173). B, Traces represent average activity for dPAG cells relative to when animals leave the safe side of the environment and start approaching the rat (top) or the shock grid (bottom). Cells defined as rat-activated (red) cells in A show increases in activity when mice approach the rat, but not when approaching the shock grid during the Preshock assay. Rat-inhibited cells (blue) do not show decreases in activity during approach to the shock grid in the Preshock habituation session. C, Quantification of neural activity changes from the traces in B, difference in mean activity of 0-2.5 s after leaving the safe side and 0-2.5 s before leaving the safe side (n same as in A, Rat 1: U = 7.24, p = 5 × 10−13, Preshock: U = 1.30, p = 0.19, Wilcoxon rank-sum test). D, Same as in A, but using correlation of cell activity with distance to toy rat instead of rat. All cells coregistered between Toy Rat and Preshock habituation session are shown (toy rat-activated n = 64, toy rat-inhibited n = 81, neither n = 253). E, Same as in B, using toy rat-activated and toy rat-inhibited cells as defined in D. F, Same as in C, but for Toy Rat session in place of Rat assay (n same as D, Toy Rat: U = 6.35, p = 2 × 10−10, Preshock: U = −0.29, p = 0.77, Wilcoxon rank-sum test). ns, not significant, ***p < 0.005.
A limitation of this analysis (Figs. 11, 12) is that it required arbitrarily separating a continuous distribution of neural activity correlation with distance to threat into three groups of cells (Fig. 11A). It is also unclear whether other behavioral variables in addition to distance to threat contribute to a shared pattern of neural activation across assays.
To address these concerns, we developed a new method called CoCA. In this method, behavioral variables (e.g., distance to threat, speed, escape, freeze, etc.) are linearly combined to produce a behavioral projection. Similarly, the neural activity variables, which consist of the fluorescence traces from each of the dPAG cells, are linearly combined to produce a neural projection. The weights chosen for these linear combinations were optimized to maximize the correlation between the neural and the behavioral projection. However, the weight selection is constrained in two important ways. First, the neural variable weights were conserved across assays for the same cell. Second, the behavioral variable weights were fixed across mice for the same assay. This method was applied to dissect shared patterns, if they exist, of neural activation across the Rat and the shock grid Fear Acquisition assays (Fig. 13A). The weights for behavioral variables are shown in Figure 13B. The weight distribution shows that distance to threat was the most relevant behavioral variable, but many other variables also had an influence on the behavioral projection. A representative scatter plot of withheld data between the neural and behavioral projections for the Rat and Fear Acquisition assays (Fig. 13C, left) shows that these two projections were highly correlated in this example mouse. A colored matrix shows the correlations between neural and behavioral projections for all mice in each of the assays (Fig. 13C, right). Furthermore, heatmaps of neural projections averaged across mice show that the neural projection had higher values (which correspond to warmer colors) in locations more near the threats (Fig. 13D). The value of the neural projection could also be used to predict position, as it was correlated with distance to threat (Fig. 13E; n = 284 cells; Rat: t = −10.37, p = 1.4 × 10−4, Fear Acquisition: t = −5.57, p = 0.0026, one-sample t test). Similar results were also found applying this analysis to the Rat and Fear Retrieval datasets (Fig. 13F,G; n = 182 cells; Rat: t = −7.70, p = 0.0046, Fear Retrieval: t = −9.87, p = 0.0022, one-sample t test). These results show that the dPAG uses a conserved pattern of neural activation to encode position and other variables during exposure to both a live predator and a fear-conditioned shock grid. Importantly, performing this analysis with control assays revealed that dPAG activity did not use a conserved neural representation for distance to threat and control stimuli (Fig. 14A; 6 mice, n = 309 cells, Rat: t = −12.0, p = 0.0012, Toy Rat: t = 6.27, p = 0.008, one-sample t test; Fig. 14B; 8 mice, n = 381 cells, Fear Acquisition: t = −1.78, p = 0.17, Preshock: t = 3.60, p = 0.037, one-sample t test; Fig. 14C; 7 mice, n = 398 cells, Toy Rat: t = −3.75, p = 0.033, Preshock: t = 7.51, p = 0.005, one-sample t test; Fig. 14D; Rat+Fear Acquisition and Rat+Fear Retrieval: n = 20 sessions, with Rat+Toy Rat: n = 12 sessions, U = 2.65, p = 0.008, with Fear Acquisition+Preshock: n = 16 sessions, U = 2.42, p = 0.016, with toy Rat+Preshock: n = 10 sessions, U = 3.39, p = 0.0007, Wilcoxon rank-sum test). Together, these data showed that dPAG cells have cells that encode distance to control stimuli (Fig. 12D); however, distance to stimuli only uses a shared neural representation when encoding two threatening stimuli (Figs. 11, 13), but not two control stimuli (Figs. 12F, 14).
dPAG displays a shared neural representation of threat proximity across the shock grid Fear Acquisition and Rat assays. A, CoCA shows conserved encoding of behaviors and neural activity, consistent across shock grid (Fear Acquisition day) and Rat exposure assays. Linear projections of behaviors (top) correlate with linear projections of single cell neural activity (bottom). B, Behavioral variables. Each cell coregistered across these two assays corresponded to an individual neural variable used to generate the neural projection. Weights used to produce the linear combinations were selected to maximize the correlation between neural and behavioral projections. Weight selection was subjected to two constraints. First, behavioral variable weights were conserved across mice (see B). Second, weights used to generate the neural projection were conserved across both assays for each cell. B, Weights of CoCA behavioral projector variables for the Rat (top) and shock grid (bottom) assays showing the influence of each behavioral variable on the behavioral projection. C, Left, Representative example correlation of CoCA projection for behavioral and neural projections for Rat (top) and shock grid (bottom) for testing data (from mouse ID#1). Each point is one time point of data. Right, Matrix displaying correlation values of CoCA projection of behavioral data with projection of neural data for testing data for individual mice (p < 0.05 calculated by bootstrapping with random weights). m/r dist., m/r angle, distance between mouse and rat and angular offset between mouse head direction and rat. D, Neural data projection overlaid as heat map for the Rat and shock grid assays, using testing data from all mice for the Rat assay (top) and Fear Acquisition (bottom). E, Bar graph represents correlations between neural projection and distance to threat for cells coregistered between Rat and Fear Acquisition assays. Data are mean ± SEM. n = 6 mice; Rat assay: t = −10.37, ***p = 1.4 × 10−4, Fear Acquisition: t = −5.57, **p = 0.0026 (one-sample t test). F, Similar to D, but for cells coregistered across Rat assay (top) and Fear Retrieval (bottom). G, Bar graph represents correlations between neural projection and distance to threat for cells coregistered in Rat and Fear Retrieval assays. Data are mean ± SEM. n = 4 mice; Rat assay: t = −7.70, ***p = 0.0046, Fear Retrieval: t = −9.87, **p = 0.0022 (one-sample t test).
Shared neural representations found using CoCA are weaker across threat + control and control + control assays than across threat + threat assays. A, CoCA shows opposite encodings of proximity to rat and toy rat. Left, Average heatmaps of neural projections are shown. Red arrows indicate location of the threats in each assay. Right, Average correlations between the shared neural projection identified by CoCA across assays (Rat assay: t = −12.0, p = 0.0012, n = 6 mice, Toy Rat: t = 6.27, p = 0.008, n = 6 mice, one-sample t test). B, Same as in A, using CoCA with Shock (Fear Aquisition) and Preshock (Shock: t = −1.78, p = 0.17, n = 8 mice, Preshock: t = 3.60, p = 0.037, n = 8 mice, one-sample t test). C, Same as in A, using CoCA with Toy Rat and Preshock (Toy Rat: t = −3.75, p = 0.033, n = 7 mice, Preshock: t = 7.51, p = 0.005, n = 7 mice, one-sample t test). D, Average correlations between shared neural and behavioral projections are lower for CoCA with control + threat and control + control (medium gray) versus threat + threat (dark gray) (Rat + Fear Acquisition and Rat + Fear Retrieval: n = 20 sessions, Rat + Toy Rat: n = 12 sessions, U = 2.65, p = 0.008, Shock + Preshock: n = 16 sessions, U = 2.42, p = 0.016, Toy Rat + Preshock: n = 10 sessions, U = 3.39, p = 0.0007, Wilcoxon rank-sum test). ns, not significant, *p < 0.05, **p < 0.01,**p < 0.005.
Discussion
The dPAG has been known for many decades to be a critical node coordinating defensive behaviors to innate threats (Lovick, 2000; Perusini and Fanselow, 2015). More recent evidence also indicates that this region plays an important role mediating defensive behaviors in retrieval of contextual fear conditioning induced either by shock (Resstel et al., 2008; Borelli et al., 2013; Mochny et al., 2013; Aguiar et al., 2014) or acquisition of contextual fear induced by predator odor (Souza and Carobrez, 2016). Here we perform the first comprehensive population-level analysis of the dPAG during exposure to threats. We show that dPAG activity encodes a wide variety of defensive behaviors in a consistent manner during exposure to both a live predator and a fear conditioned shock grid. Cells also encoded angular offset between mouse head direction and threat as well as rat speed. We also show that distance to threat in both assays is encoded using shared patterns of neural activity. Last, we show that dPAG activity can predict escape from a live predator several seconds before escape is observed, which indicates that the dPAG may possibly have a role in computing the decision or timing of the escape. Together, these data show that the dPAG uses largely conserved activation patterns to encode defensive behaviors and distance to threats, both for innate and conditioned threats, indicating that this nucleus may use a common pathway to control survival strategies during exposure to numerous threat modalities.
dPAG representation of distance to threat and defensive behaviors
dPAG cells strongly encoded both defensive behaviors and kinematic parameters, such as distance to threat. These two features are represented independently in the same dPAG ensembles (Fig. 5E). The behavioral GLM was trained after regressing out the distance-to-threat variable; thus, significant GLM weights for behaviors cannot be explained by representation of distance to threat (Figs. 3, 4). Furthermore, GLM weights for distance and behaviors were only weakly correlated (Figs. 3I, 4F). Additionally, responses to behavioral events were not modulated by the distance to threat (Fig. 5A–D). Last, we divided cells into three groups: cells displaying positive, negative, or no modulation by distance to threat (as determined by their GLM weights for this variable). In all three groups, a similar proportion of cells showed significant modulation by approach or other behaviors (Fig. 5E). Similar results were obtained for all behaviors in both assays, showing that distance to threat and behavioral encoding occurred in the same cells independently. Together, these data show that the representation of kinematic and behavioral variables occurs in the same cells, but as independent codes.
Encoding of threat-related features in the dPAG is conserved across assays
Evidence from multiple data streams, such as in vivo electrophysiology and immediate early gene expression studies, have shown dPAG cells are activated by numerous threats, including predator cues (Deng et al., 2016), carbon dioxide (Johnson et al., 2011), fear conditioned contexts (Carrive et al., 1997), and fear conditioned auditory tones (Watson et al., 2016). These data strongly indicate that the dPAG encodes features that may be relevant to influence ongoing defensive behaviors. However, the identity of which threat-related features are represented in dPAG remained unknown. Here we show that dPAG cells represent proximity to threat, predator movement onset, and angular offset between the head direction of the mouse and the location of the threat. One previous report had identified cells that encoded distance to a predator in the dPAG (Deng et al., 2016), but the consistency of this representation across threat modalities has not been reported previously. Our observation that representation of these threat-related variables is conserved across assays is in agreement with the hypothesis that encoding these features may be important to guide defensive behaviors. Importantly, shared patterns of neural activity were used to represent variables across threat assays, but not control assays for both behavioral (Figs. 6D, 8, 9) and kinematic variables (Figs. 11–14).
Furthermore, coregistered cells could be grouped into conserved clusters across assays based on their neural response properties during defensive behaviors. Intriguingly, the four clusters identified in the GMM analysis (Fig. 7C) show that the clusters for each behavior are of a similar size, and thus have a comparable number of cells. These data indicate that a wide variety of behaviors (not just escape) are encoded in the dPAG with comparable functional specialization. This result agrees with data showing, for example, that optogenetic activation of the dPAG produces diverse defensive behaviors (Deng et al., 2016). While dPAG cells encoded various behaviors and threat-related features, that does not imply that this region controls all these behaviors or that it uses this kinematic information to influence behavior.
It is noteworthy that the dPAG is a heterogeneous and genetically diverse structure even within the same column, with various cell types releasing different neurotransmitters, such as nitric oxide, substance P, and cholecystokinin (Silva and McNaughton, 2019). It is possible that this genetic heterogeneity may underlie the diversity of dPAG activity patterns during threat exposure (Fig. 7). Future studies are needed to investigate this tantalizing hypothesis.
dPAG as part of the fear memory circuit
Our data show that the dPAG encodes a wide variety of defensive behaviors, during exposure to innate and conditioned threats. In line with this view, exposure to predatory threats (Canteras and Goto, 1999; Dielenberg and McGregor, 2001; Cezario et al., 2008), predator-conditioned context (Cezario et al., 2008), and shock-conditioned context and auditory tones (Carrive et al., 1997; Watson et al., 2016) activate the dPAG. The literature has focused on the central amygdalar nucleus–ventrolateral PAG pathway as the central player in shock-based fear conditioning (Tovote et al., 2016), and the dPAG has received little attention in this matter. Nevertheless, several studies showed that the dPAG influences fear responses in shock-based fear conditioning (Resstel et al., 2008; Borelli et al., 2013; Mochny et al., 2013; Aguiar et al., 2014). In this regard, the prelimbic cortex is a likely candidate to provide information concerning shock-based fear conditioning to the dPAG. The prelimbic cortex supplies considerable glutamatergic inputs to the dPAG, and prelimbic cortex inactivation reduced freezing to both a tone and a context that had been previously paired with footshock but had no effect on freezing to a predator (Corcoran and Quirk, 2007). Because of the existence of reciprocal connections with forebrain circuits (Motta et al., 2017), the dPAG occupies a particularly privileged position to play a critical role in both conditioned and unconditioned threat imminence continuum processing.
Reexamining the role of the dPAG in escape
Prior work has strongly implicated the dPAG to escape to innate threats (Tovote et al., 2016; Evans et al., 2018). Optogenetic dPAG activation causes escape in the absence of threats (Deng et al., 2016), and dPAG cells have been shown to fire during escape from a live predator (Deng et al., 2016) or a looming stimulus (Evans et al., 2018). dPAG activation has also been observed in humans exposed to imminent threats, such as proximity to a spider (Mobbs et al., 2010). Our results agree with these data, as we also found robust dPAG activation during escape from an innate threat consisting of a predatory rat. Our data also expand beyond published results because we found that dPAG cells fire during escape from a shock grid, which is a conditioned aversive stimulus. These data suggest that dPAG cells may also participate in controlling escape from certain modalities of conditioned threats. Shared patterns of dPAG activation for defensive behaviors across assays were found through correlated GLM weights, functional clustering, and across-assay decoding, indicating that this result is fairly robust.
Current models largely view the dPAG as a downstream effector region, such that the decision to escape is computed in other upstream regions, while dPAG activation merely executes the escape and controls escape vigor (Tovote et al., 2016; Evans et al., 2018). Our data challenge this model, as we provide the first evidence that dPAG ensembles can predict escape from a live predator >3 s before the escape is initiated. This result indicates that the dPAG may be involved in computations that affect the timing and the decision to escape, rather than only serving as a downstream effector region to initiate flight. Two important features of the current work may explain the discrepancies between our results and prior reports. First, we performed population-level analysis with large dPAG ensembles, without which it would not have been possible to predict escape several seconds in the future using dPAG activity. Second, in our work, the mouse voluntarily controlled its proximity to the threat by choosing to approach the rat or the shock grid. In this situation, the mouse has more time to evaluate its action selection, as it also voluntarily decides when to start running away from the rat. Perhaps, in these situations, the dPAG may affect the decision to escape. In contrast, in prior work, escape-inducing looming stimuli were presented at predetermined times, independent of the mouse's choice to approach the threat (Evans et al., 2018). This assay induced rapid and vigorous escape in which the dPAG showed increased activity only during escape execution, but did not affect the decision to escape. Together, these results suggest that the dPAG may influence the decision or propensity to escape in slower, self-guided assays in which the animal has time to freely navigate and decide when to approach, flee, or freeze, while the dPAG may be limited to only executing escape in situations where threats are unpredictably presented to the animal, evoking rapid flight.
Footnotes
This work was supported by the National Institute for Mental Health R00 MH106649 and R01 MH119089 to A.A.; Brain and Behavior Research Foundation Grants 22663, 27654, 27780, and 29204, respectively, to A.A, F.M.C.V.R, W.W, and J.C.K.; National Science Foundation NSF-GRFP DGE-1650604 to P.J.S.; UCLA Affiliates Fellowship to P.J.S.; Hellman Foundation to A.A.; and Fundação de Amparo à Pesquisa do Estado de São Paulo Research Grant 2014/05432-9 to N.S.C. and Grants 2015/23092-3 and 2017/08668-1 to F.M.C.V.R.
The authors declare no competing financial interests.
- Correspondence should be addressed to Avishek Adhikari at avi{at}psych.ucla.edu
