Mapping and Deciphering Neural Codes of NMDA Receptor-Dependent Fear Memory Engrams in the Hippocampus

Mapping and decoding brain activity patterns underlying learning and memory represents both great interest and immense challenge. At present, very little is known regarding many of the very basic questions regarding the neural codes of memory: are fear memories retrieved during the freezing state or non-freezing state of the animals? How do individual memory traces give arise to a holistic, real-time associative memory engram? How are memory codes regulated by synaptic plasticity? Here, by applying high-density electrode arrays and dimensionality-reduction decoding algorithms, we investigate hippocampal CA1 activity patterns of trace fear conditioning memory code in inducible NMDA receptor knockout mice and their control littermates. Our analyses showed that the conditioned tone (CS) and unconditioned foot-shock (US) can evoke hippocampal ensemble responses in control and mutant mice. Yet, temporal formats and contents of CA1 fear memory engrams differ significantly between the genotypes. The mutant mice with disabled NMDA receptor plasticity failed to generate CS-to-US or US-to-CS associative memory traces. Moreover, the mutant CA1 region lacked memory traces for “what at when” information that predicts the timing relationship between the conditioned tone and the foot shock. The degraded associative fear memory engram is further manifested in its lack of intertwined and alternating temporal association between CS and US memory traces that are characteristic to the holistic memory recall in the wild-type animals. Therefore, our study has decoded real-time memory contents, timing relationship between CS and US, and temporal organizing patterns of fear memory engrams and demonstrated how hippocampal memory codes are regulated by NMDA receptor synaptic plasticity.


Introduction
The major obstacle in understanding of how memory works in the brain is the lack of description of the real-time brain activity patterns and its organizing principles underlying real-time memory process. Most electrophysiological studies examined single unit activity using peri-event based data averaging methods over multiple trials to characterize response or tuning properties of the recording neurons. This data-averaging practice has, unfortunately, led to the loss of crucial information regarding transient activity patterns and fundamental dynamics that underlies realtime memory code. From a structural perspective, memory is believed to be the storage of acquired information in a form of synaptic connectivity patterns via by NMDA receptor-dependent synaptic plasticity [1][2][3][4][5][6][7]. Yet, from a temporal dynamic perspective, memory traces are real-time transient activity patterns at the moment when the stored information is re-emerged by recall [8,9]. Information encoding different memory traces evolves dynamically from a moment to moment as the brain retrieves them. Currently, little is known regarding the real-time brain activity patterns of associative memories and what real-time memory engram looks like and how they are organized. There is an emerging interest in applying large-scale neural recording and powerful mathematical decoding approaches to seek out brain activity patterns or brain activity maps hidden inside the large datasets [10][11][12][13].
Extensive studies have demonstrated that neurons in the hippocampus encode spatial information about animals' location, categorical information about people and items, and timing relationships [14][15][16][17][18][19][20][21], yet the hippocampus is also known to play a central role in trace and contextual fear memories. Fear memory is associative in nature, and fear memory-associated behaviors are well defined and widely studied across a variety of animal species. A wide range of lesion studies on mice, rats, and primates show that trace fear conditioning, especially with longer trace interval, requires the structural integrity of the hippocampus and related circuits [22][23][24][25][26][27][28][29][30]. In rodents, fear memory can be conveniently assessed by measuring the amount of freezing [23][24][25], startle responses [31,32], or changes in heart rates and heart rate variability [33]. Although behavioral changes indicate the formation of new memories, many mnemonic thoughts and precise contents of memory in animals remains internal and largely inaccessible to outside observers. It has been suggested that fear memory may contain information about CS, US, environmental context, self-location, and the time relationship (or time trace interval) between hearing the tone and expected arrival of foot-shock, etc. [23][24][25][26][27][28][29][34][35][36]. However, by assessing freezing behavior from outside, such detailed information generated in the memory circuits can only be inferred. At present, very little is known regarding many very basic questions such as: are fear memory traces retrieved during the freezing state or non-freezing state of the animals? How can different memory contents or traces come together to generate a holistic memory engram in real-time? What are the temporal organizing principles underlying associative memory engram? How fast can memory traces be retrieved in the mouse hippocampus? In seeking the answers to these fundamental questions, in the present study we combine largescale decoding technologies with inducible NMDA receptor knockout approach to map out hippocampal activity patterns of trace and contextual fear memory traces in the CA1 region. We provide the first detailed description about how the content and temporal dynamics of real-time memory engrams are organized by the NMDA receptor-mediated synaptic plasticity.

CA1 ensemble representations for conditioned and unconditioned stimuli
We implanted 128-channel electrodes into mouse hippocampal CA1 regions and recorded large numbers of single units from both control and knockout mice. The location of electrodes in the CA1 was initially assessed by the presence of characteristic oscillations sharp-wave associated ripples and theta (Figure 1), and then confirmed by post-recording histology. We recorded a total of 1147 CA1 units from five control mice and 1153 CA1 units from five inducible-and forebrain-specific NMDA receptor knockout mice. We used the trace fear conditioning paradigm to train both the inducible knockout mice and their control littermates (5-6 month old males). The trace conditioning contained a neutral tone (CS, 2 sec of 2.8 kHz pure tone at 85 dB) which was followed by a delayed mild foot shock (US, 0.75 mA, 0.3 sec) with a fixed 20 seconds ''trace'' time interval between CS and US ( Figure 2A). Conditioning of this long-duration time trace memory is highly sensitive to hippocampal-lesion [23][24][25][26][27][28][29]34,35], and the paradigm permits closer examination of various memories such as the memories for the conditioned tone, shock event, their causal relationship, and time information for the tone and shock intervals, etc. [36]. To facilitate the physiological identification of CS and/or US response units in the CA1 of the hippocampus, we used seven repetitions of CS/US pairing protocol which produced strong contextual and trace fear memories in the control mice as assessed in 1-hr retention tests. We also examined the genetic effects of forebrain-specific NMDA receptor deletion on trace fear conditioning by feeding the inducible knockout mice with doxycycline 5 days prior to fear conditioning. As expected, dox-induced NMDA receptor knockout caused significant behavioral impairment in learning trace fear conditioning as revealed by reduced contextual freezing ( Figure 2B, Student's t test p ,0.05) and trace retention freezing ( Figure 2C, Student's t test p ,0.01).
Our 128-channel electrode arrays [36][37][38][39] permitted us to record simultaneously, on average, the activity of 229610 units per control mouse or 23169 units per mutant mice from the CA1 region of the hippocampus bilaterally as these animals underwent the acquisition and 1-hr retention tests. CS and US stimulation evoked firing changes in many of the recorded CA1 units ( Figure   3A). To characterize responsiveness of the recorded CA1 units, we conducted peri-event spike raster and histogram analyses by averaging neural responses over seven CS or US during learning, or the conditioned tone during recall trials using the stimulus onset as time zero. Peri-event histograms showed that the conditioned tone stimulus (CS) and foot-shock (US) evoked significant changes in firing of some CA1 units during learning. Some units fired selectively to the CS only during learning ( Figure 3B), whereas other units responded only to US ( Figure 3C). Moreover, some units changed their firings to both CS and US during learning, but not during tone recall ( Figure 3D). In addition, we found a certain percentage of the units would respond to both CS and US during learning and the recall tone during retention tests ( Figure 3E). By using the hierarchical clustering analysis [38,40], we examined and compared that the entire CA1 population activity responsiveness between the genotypes (Figure 4). The overview of the all datasets suggest that the CA1 populations in the control mice exhibited categorical and combinatorial coding patterns in responding to fear conditioning, ranging from broadly associative response units (listed on the top of the plot) to a set of specific response units (listed at the lower portion of the plot) ( Figure 4A). Similarly, the mutant mice also showed categorical and combinatorial patterns, but with a lower percentage of responsive cells, especially those combinatorially associative units ( Figure 4B). As a whole, approximately 22.9% (263 out of 1147 units) of the recorded units from the control mice and 12.3% (142 out of 1153 units) of the recorded units from knockout mice responded to CS, whereas 45.9% (523 out of 1147) of the recorded units from controls and 18.2% (210 out of 1153) of units from mutants reacted to CS-paired foot-shock. Among CS-responding units, 5.1% (59 out of 1147 units) was specific to CS (not responding to US or Tone at recall) in the controls, whereas 5.5% (63 out of 1153 units) was specific to CS in the mutants. The total percentages of units belonging to the specific groups (i.e. combining CS-specific, CS-paired US-specific, and recall tonespecific units) from the wild-type and mutant CA1 populations were 33.4% and 24.1%, respectively.
Because associative binding of distinct information is a hallmark feature of Pavlovian fear conditioning, we examined more closely the associativeness of the CA1 responding units to both tone and shock. We found that controls had a total 192 units (out of 1147) responded to both CS and US during learning (16.7%) ( Figure  4A), whereas the mutants had only 60 units exhibited such associative responses (5.2%) ( Figure 4B), reflecting a reduced associative binding between the tone and foot-shock in mutant CA1. Furthermore, the control group had a total of 46 units out of 1147 (4.0%), exhibiting significant responses to both CS and US during learning as well as to the conditioned tone at trace retention tests ( Figure 4A). By comparison, the mutant mice had only 18 such units responding to all three stimuli (1.6%) ( Figure 4B). Interestingly, the tone at trace retention test resulted in 11.4% (131 out of 1147) of the control units changing their firings ( Figure 4A), whereas the mutant mice had 12.8% (148 out of 1153) of the units responded to the recall tone ( Figure 4B). This suggests that although the similar percentage of cells in mutant mice responded to tone during the retention tests, many of these units did not participate in CS and/or US learning.
NMDA receptor knockout prevented the formation of real-time associative memory traces trial, we employed Multiple Discriminant Analysis (MDA) to obtain statistical pattern classifications of neural ensemble activities representing distinct stimulus categories (see the schematic outline of MDA method and steps for obtaining dimensionality reductionbased pattern classification as described in Figure 5). The MDA method provides not only quantitative classification of neural ensemble patterns but also intuitive visualization of these patterns in encoding subspaces [36][37][38]40,41]. CA1 ensemble activity prior to CS stimulus presentation formed the Rest ellipsoid, whereas the activity during CS or CS-paired foot shock formed the CS tone ellipsoid and US shock ellipsoid, respectively ( Figure 6A). Interestingly, CA1 ensemble responses from both control and mutant mice produced reliable classifications for the rest, CS, and US ensemble patterns ( Figure 6A and 6B). However, the discriminant distances between the Rest and US ellipsoids or between the CS (Tone) and US ellipsoids in the mutant group are significantly smaller than those of the control littermates ( Figure  6C, Student's t test p,0.01).
By further coupling MDA with the sliding-window method (using 20 millisecond steps) ( Figure 7A), we investigated dynamic patterns of CA1 cell population ensembles during the learning phase. We observed that in response to the first tone (before the first paired foot shock arrived), CA1 cell population from the control animal produced a CS ensemble trajectory from the Rest ellipsoid ( Figure 7B, left subplot), whereas the CS-paired US (20 seconds after the CS) elicited US ensemble trajectory ( Figure 7B, right subplot). Interestingly, as CS/US pairings were repeated over learning trials, the subsequent conditioned tone produced either robust CS simple traces or CS-to-US associative traces ( Figure 7C left subplot). Moreover, the CS-paired foot-shock frequently evoked US-to-CS associative trajectories instead of simple US trajectories ( Figure 7C right subplot).
By comparison, fear conditioning in knockout mice still triggered simple CS or simple US trajectories at the first CS/US representation ( Figure 7D). However, the repeated CS/US pairings over trials rarely led to US-to-CS associative trajectories or CS-to-US associative trajectories ( Figure 7E). Of five recorded control mice over seven learning trials, foot shock-induced US-to-CS associative traces occurred predominantly, for 25 times out of a total of 35 CS-paired US presentations ( Figure 7F, see red rectangles in the left panel). In contrast, only four such associative traces were observed in the knockout mice ( Figure 7F, right panel). Similarly, while CS-induced CS-to-US associative traces (red rectangles) were frequently observed in control mice (15 times out  Figure 7G, right panel). Therefore, these results demonstrate that disabling of NMDA receptor function prior to learning impaired the formation of realtime associative ensemble traces in the hippocampus.
On-line reverberations during the learning phase are diminished by NMDA receptor knockout To further investigate CA1 activity dynamics during the learning phase, we used the MDA/sliding-window techniques and scanned through the spike rasters during each of seven CS/ US pairing trials (60 second) ( Figure 8A). This sliding-window technique revealed that the CA1 cell population of control mice produced immediate on-line reverberations (we used this term as supposed to off-line replay after the training or during sleep). Examples of the time points at which CA1 ensemble trace reverberated are marked by triangles and diamonds shown in Figure 8A. The reverberated traces included the CS simple traces (blue triangles), US simple traces (red triangles), and US-to-CS associative traces (red diamonds) as well as CS-to-US associative traces (blue diamonds) ( Figure 8B). Moreover, numbers of trace reverberations in control mice became more prevalent as CS/US pairings were repeated ( Figure 8C). By comparison, inducible knockout of the NMDA receptors had drastically reduced reverberation of CS or US traces during the learning phase ( Figure 8D). Statistical analyses show that the control mice exhibited much higher amount of total reverberations in comparison to the mutant mice ( Figure 8E, left plot, Student's t test p,0.001). The differences in reverberation were observed when we analyzed subclasses of various trace types, such as CSpaired simple US traces, US-to-CS associative traces, and CS-to-US associative traces ( Figure 8E, right plot, Student's t test p,0.05, p,0.001, p,0.001 respectively). More importantly, the compositions of reverberated trace types are also different between the genotypes. Trace reverberations in the control mice consisted of 67.4% simple traces (28.1% simple CS, 39.3% simple US traces) and 32.6% associative traces (15.2% CS-to-US traces, 17.4% USto-CS traces) ( Figure 8F, left pie chart). For the knockout mice, trace reverberation was not only lower in numbers but also overwhelmingly belonged to simple US or CS traces (93%) ( Figure  8F, right pie chart).
The group data analyses further revealed that while the control group had more trace reverberations especially over initial multiple trials, the mutant group did not exhibit any increase ( Figure 8G). Furthermore, CS/US pairing repetition over trials seemed to preferentially increase reverberations of associative traces in the control mice ( Figure 8H, blue line). By contrast, the mutant group had little reverberation ( Figure 8H, red line).
In fear conditioning literature, immediate freezing is often used as an indication for learning. Therefore, we examined how ensemble trace reverberation reflects this form of behavioral learning. We found that the reverberation in the control mice is correlated nicely with the amount of immediate freezing during this learning phase (r = 0.84, p,0.05) ( Figure 8I), whereas the mutant mice had low reverberations and low freezing ( Figure 8J). This correlation provides additional support for the validity of the MDA-based decoding approach.
To identify the neurons that underlie these simple and associative trajectories, we used the time points of reverberated traces revealed by MDA methods as time zero for a new round of peri-event raster and histogram analyses. Indeed, we identified many CA1 units which responded at these time points (Figure 9). These distinct units, through the temporal co-spiking dynamics, collectively produced robust real-time ensemble activity peaks for CS ( Figure 9A), US ( Figure 9B), US-to-CS ( Figure 9C), or CS-to-US ( Figure 9D), respectively (raster plots on the left column of Figure 9A-D). To further evaluate these units' contribution to realtime ensemble trajectories, we employed spike shuffle techniques to examine their trajectories in the MDA subspaces. Indeed, when these responsive units were shuffled, each corresponding trajectory was greatly diminished in comparison to that of before shuffle ( Figure 9, prior to shuffle MDA plot on the left, after shuffle MDA plot in the middle). To rule out the overall non-specific effects of random shuffling of spikes, we shuffled an equal number of nonresponsive units and found that it did not produced any significant effect on these trajectories ( Figure 9, right MDA plot column). Collectively, these analyses strongly suggest that these dynamic trajectories were the results of collective co-firing changes of these perspective units. There was a significant difference in contextual freezing between the control (n = 4, 52%66%) and mutant mice (n = 5, 21%65%). Error bars represent SEM; *p,0.05. (C) Impaired trace fear retention in the mutant group as compared to the control group. Freezing prior to recall tone and after the tone at 1-hr trace recall in the control and mutant mice. The tone was presented for seven times (trials) with a 1-3 min random time interval. There was a significant difference in tone-induced freezing between control (51%67%) and knockout mice (31%64%), **p,0.01. doi: 10 [34][35][36]42]. To qualify any ensemble traces observed during learning as memory traces, these traces must re-appear upon cueinduced memory recall. Therefore, we subjected the mice to trace retention test by placing the animals to a non-conditioned chamber different from the original conditioning chamber. 2-sec tone was played in the absence of foot shock, and it was then repeated for seven times with the 1-3 minutes randomized interval between each tone representation. Prior to the exposure to the recall tone, the mice had a low amount of freezing ( Figure 10A, see the top black bar). Shortly after hearing the conditioned tone, the control animals exhibited significantly more freezing ( Figure 10A, see the orange bar). Our MDA/sliding window decoding analysis revealed that the conditioned tone led to a string of CA1 memory traces re-appeared in the trace recall session ( Figure 10A, see triangles and diamonds at the bottom of the raster plot). In the control mice these transient dynamic trajectories included simple CS traces and US traces, as well as US-to-CS associative memory trace and CS-to-US associative traces ( Figure 10B). On average, various memory traces were retrieved with a rate up to 17 memory traces per min in the controls (1-min recall session) ( Figure 10C).
By contrast, the numbers of memory traces retrieved from the knockout mice were much lower, and few occurred in form of associative trajectories ( Figure 10D). On average, the control mice retrieved 10.4 memory traces per recall, verses 5.5 memory traces per trial in the mutant mice ( Figure 10E). Interestingly, the memory trace interval analysis revealed that memory retrieval followed the exponential decay-like process ( Figure 10F). This exponential decay distribution reflected that retrievals of memory traces were formatted closely together in time, with two or more traces temporally clustered together (see Figure 10C). The shortest time interval between two individual memory traces was 0.76 sec. In comparison, the retrieved memory traces from mutant mice were significantly degraded in term of its bursting temporal format ( Figure 10G). A two-sample Kolmogorov-Smirnov test indicated that the inter-trace intervals of control and mutant mice had different distribution (p,0.001).
One unique property of Pavlovian trace conditioning is that it offers the opportunity to examine whether and how the brain generates real-time memory traces of ''what at when'' time information regarding the temporal relationship between hearing the tone and predicting the timing at which US event should follow [36,42]. By plotting the time distribution of various realtime memory traces retrieved over the trace retention test period, we found that the US memory traces exhibited a distinct recall peak in the control mice, centered around 20 seconds after the offset of the 2-sec tone ( Figure 10H, blue line). This showed that the hippocampus in the control mice consistently predicted the arrival time of foot-shock at that particular moment. The mutant mice failed to show the memory trace of ''what at when'' time relationship ( Figure 10H, red line).
On the other hand, the time distribution plot for CS traces over the 1 min recall period in the control mice did not reveal any significant peak, suggesting that this is a specific phenomenon in term of causal relationship between CS hearing the tone and US memory trace retrieved with 20 seconds of time interval ( Figure  10I).
As a group, all five control mice showed consistent retrievals of the US ensemble traces around the 20-sec time interval ( Figure  10J, red rectangles). In contrast, all of the knockout mice failed in predicting US events at this moment in time ( Figure 10K). Therefore, our results demonstrated real-time memory traces encoding ''what event at when'' are critically dependent on the NMDA receptors.

Real-time hippocampal memory traces during contextual recall test
Trace fear conditioning also offered us an opportunity to examine contextual fear memory in the CA1 [8,[43][44][45][46]. However, the animals can recall contextual fear memory at various time points, the traditional peri-event spike histogram analysis could not be readily performed due to apparent lack of the time zero for the averaging spike raster. We took advantage of MDA/sliding window decoding technique and mapped the retrievals of fear memory traces during the entire contextual retention test period. We found that immediately after the control mouse entered the conditioning chamber (in the absence of the tone), fear memory traces were retrieved within a few seconds (see examples in Figure  11A). The re-emergence of fear memory traces preceded the freezing behavior (see the orange bar above the raster and triangle and diamond below the raster in Figure 11A). Similar to those observed in trace retention tests, the retrieved memory traces during contextual recall included various simple and associative memory traces ( Figure 11B). In all five control mice, we have observed the recalling of various memory traces (See the first 60 sec period after the mice entered contextual retention chamber in Figure 11C). In contrast, the retrievals of CS or US-related traces in the knockout mice were significantly reduced during the contextual retention tests ( Figure 11D).
Interestingly, the number of memory traces retrieved during contextual recall is correlated with the amount of contextual freezing within individual control animals ( Figure 11E). The mutant mice had lower amount of freezing and reduced numbers of memory trace retrieved ( Figure 11F). Linear regression analysis of the group data further showed the correlation between the numbers of memory traces retrieved and freezing responses ( Figure 11G, r 2 = 0.73, p,0.01), with control mice at the higher end of the freezing and memory trace recalled and mutants at the lower end of the plot. Calculation of total memory traces retrieved during the entire 5-min contextual recall showed that control mice recalled at 48.663.7 traces, whereas the mutant mice had only 16.065.4 traces ( Figure 11H). We noted that about 20 percent of memory traces retrieved in the control mice were in forms of associative memory traces during contextual retention tests ( Figure  11H, solid portion of the bar), but very few was in the mutant mice. This is consistent with impaired formation of associative memories observed in learning phase in these mutant mice.
To further quantify the temporal effects of NMDA receptors on retrieval dynamics of contextual memory recall, we plotted the inter-memory trace-intervals among all memory traces retrieved during this 5-min retention test period. The control mice recalled fear memories, often with two or more memory traces tightly clustered together ( Figure 11I). The shortest inter-memory time interval was at 1.02 sec. However, the mutant mice lacked exponential decay-like distribution and were significantly impaired

Temporal map of CA1 fear memory codes' contents and associational dynamics
To seek out additional understanding of the temporal organizing patterns in relation with dynamics of associative memories, we analyzed the various relationships among and between complex associative memory traces and simple memory traces during both contextual and traced recall. We first examined how many associative memory traces (i.e. CS-to-US or US-to-CS associative memory traces) were retrieved during the entire 5-minute contextual recall session. We found that the control group retrieved a total of 47 such associative memory traces recalled (shown as red or blue diamonds) ( Figure 12A). Most interestingly, these associative memory traces often retrieved in alternating doublets or more. We referred the alternating retrieval pattern when a given associative memory trace was followed by another associative memory trace. On the other hand, when a given associative memory trace was followed by other simple traces, we termed it as intertwined retrieval pattern ( Figure 12A). By comparison, the mutant group had a much lower number, with only 5 associative memory trace-like patterns detected (all came from one knockout mouse) during this 5-min contextual retention test ( Figure 12B).
Next we examined the temporal associational relationships among simple CS traces and US memory traces retrieved during contextual recall. We found that many CS or US traces were retrieved in a temporally associational manner as evident from the three following analyses: 1) the intertwined retrieval pattern, again, it referred to the temporal format when a given CS or US simple memory trace was followed by a different memory trace (either a simple or associative memory trace). 2) US trace-based alternating retrieval pattern (when a simple US memory trace was followed by a simple CS memory trace); 3) CS trace-based alternating retrieval pattern (when a given CS trace was followed by a US memory trace).
In the control mice, the vast majority of memory traces retrieved occurred not only closely in temporal domain (in clusters), but also using the intertwined format ( Figure 12A). These intertwined pairs (underlined) also exhibited exponential decay ( Figure 12C). In comparison, memory retrieval in the mutant mice had a loose temporal format, with few doublets or triplets ( Figure 12B,E). For those simple traces retrieved in the mutant CA1, they showed little intertwined retrieval patterns and lacked exponential decay distribution ( Figure 12C, see the red histogram). There was a significant difference in the total pairs of intertwined memory recall between the genotypes ( Figure 12D, Student's t test p,0.05). In the control, about half of the all memory traces retrieved occurred as clusters in doublet or triplet formats, whereas the mutant mice generally lacked doublet or triplets ( Figure 12E).
Similarly, in control mice, the vast majority of the simple US traces followed by simple CS memory traces or vice versa (alternating retrieval pattern) within a 6.6 sec time window ( Figure  12F, blue histogram). Again, this temporally alternating recall patterns were impaired in the mutant mice ( Figure 12F, red histogram. Two-sample Kolmogorov-Smirnov test, p = 0.0042). The average occurrence in CS/US alternating recall formats within the interval less than 7 sec in the control mice was 2.7 times per minute, whereas the knockout mice had only about 0.3 time per minute ( Figure 12G, Student's t test p,0.01).
To investigate whether the intertwined and alternating recall dynamics represent a general temporal mechanism for ensuring associative memory recall associational in its time domain, we further analyzed the traced fear recall in both control and knockout mice ( Figure 13A,B). Again, we found that CS-to-US and US-to-CS associative traces were mostly observed in the control mice ( Figure 13A), but rarely in the mutant mice ( Figure  13B). When using 7 sec as the time threshold, about 36% of various memory traces were dynamically retrieved in doublet, triplet or more in the control mice, whereas the mutant mice had little such bursting manner ( Figure 13E). More importantly, these bursting memory retrieval patterns occurred in intertwined formats (underline) in the control mice as confirmed by the exponential distribution of intertwined recalling pairs ( Figure 13C, blue histogram). The mutant mice failed to exhibit any significant intertwined dynamics ( Figure 13C, red histogram). The total numbers of the intertwined recalling occurrences in the control was at 4.1 times per min, vs. 0.9 times per min in knockout mice ( Figure 13D, Student's t test p,0.01).
Finally, our analysis also shows that alternating retrieval patterns (upper line) between CS/US traces were prominently present in the normal hippocampus, as evident from the exponential distribution in the control mice ( Figure 13F). About 75% of alternating traces in which US memory traces were followed closely by a CS memory trace or vice versa, within a time window of 6.7 sec in the control mice, whereas the mutant mice had very limited such recalling patterns ( Figure 13F). The total number of alternating recalls in the control mice was at 3.3 times per min, whereas the number for the mutant mice was at 0.7 times per min ( Figure 13G, Student's t test p,0.01). Taken together, these results showed that retrievals of Pavlovian fear conditioning memories in the control mice were highly associational in its time domain, with many memory traces organized in the intertwined or alternating recalling formats. Inducible NMDA receptor knockout disrupted such temporal organization of memory traces in the CA1.

Discussion
Mapping brain activity patterns and then understanding their underlying meaning are an emerging interest among many brain researchers [10][11][12][13]. By employing large-scale neural recording techniques and pattern classification mathematical tools, we systematically mapped and decoded real-time CA1 fear memory traces during learning as well as trace and contextual recall. The initial successful decoding of real-time fear memory's neural code was first reported in the wild-type mouse hippocampus [13,36]. Here, we have further examined the hippocampal fear memory engrams in inducible and forebrain-specific NMDA receptor knockout mice. From our present analyses of large datasets using dimensionality-reduction pattern classification techniques, we have made follows sets of novel observations that have not been reported previously: First, CS or US stimulation can still generate significant CA1 ensemble representation in the knockout mice, suggesting that a good level of perceptual real-time patterns was produced with discriminative separation even in the absence of the NMDA receptor in the forebrain. Yet, while the hippocampus of control mice readily generated US-to-CS or CS-to-US associative memory traces upon repeated pairing of CS and US stimuli, the mutant mice failed to do so. To our knowledge, this represents the first network-level evidence for the NMDA receptor's role in the formation of real-time associative fear memory traces in the hippocampus. It is interesting to note that CS-paired foot-shock produced a smaller set of units in the mutant mice. This may suggest several possibilities including the reduced perceptual representation of foot shock in the mutant CA1. However, given the highly convergent inputs from the multi-modal cortex to the hippocampus, these reduced CS-paired foot shock responsive cells in the mutant CA1 region may reflect the loss of associative binding or integration with contextual or environmental information with US stimuli. It might be of interest to examine the  Some of the foot-shock memories most likely contain the context or location information that integrates ''what and where'' memory traces. Due to the smaller size of the mouse fear conditioning chamber and high proportion of freezing during learning and recall, place cell property, which usually requires running behavior for characterizations, has prevented us to study in the current protocol. Additional experimental designs [47] will help address this issue in future.
Second, various CA1 ensemble traces underwent on-line reverberation during the training in the control mice. Interestingly, the mutant mice had greatly reduced reverberation. Recent studies have reported that trace fear conditioning caused changes in membrane potentials [48] or synaptic potentiation in the hippocampus and prefrontal cortex [49,50]. The potentiated excitability in neurons or their synapses may contribute to the initiation of pattern reverberations in the network. The strong correlation between on-line memory trace reverberation and immediate freezing supports the idea that both phenomena reflect learning process, with one at neural population level, other at behavioral level. Diminished memory trace reverberation explains well why blockade of NMDA receptor pathway would impair fear memory formation or consolidation [22,31,[51][52][53]. By restricting inducible knockout of the NMDA receptor to the learning or postlearning consolidation window, the distinct temporal roles of the NMDA receptor can be further studied [54][55][56].
Third, our study has provided the activity map on real-time CA1 fear memory contents and their temporal dynamics. We showed that memory retrievals conformed to a set of dynamic patterns or templates. Overall, fear memory traces were retrieved in clustered or burst templates in the CA1 of control mice, at a frequency range of 7.2 to 13.9 times/min during contextual or tone recall. In contrast, the inducible knockout mice had much lower memory content, limiting to the 0.8 to 7.0 traces/min. Its poor memory content is further evident from the near absence of US-to-CS or CS-to-US associative memory traces during contextual or tone recall. This suggests the qualitative differences in holistic memory representation between the mutant mice and their control littermates.
Fourth, we have uncovered two important temporal dynamics in organizing associative memory patterns in the hippocampus: 1) Many distinct memory traces are retrieved in a highly intertwined manner, and invariantly in a memory trace bursting format (as shown by an exponential decay process); 2) Most of simple memory traces in the control mice were also retrieved in an alternating fashion (i.e. retrieving a US trace and then followed by a CS trace, or vice versa within the 7 sec time-domain, also an exponential process). Such intertwined and alternating retrieval dynamics are consistently observed during both the contextual and tone recall in the control mice. This recalling dynamics may serve as a key temporal organizing mechanism in generating the holistic memory engrams associational in its time domain. It will be of great interest to examine how various manipulations of experimental conditions, such as pre-context exposure or reconsolidation protocols, change fear memory engram [57][58][59].
Finally, the intertwined and the alternating retrieving dynamics are almost completely missing in the inducible NMDA receptor knockout mice. Despite CS or US traces could still be detected the mutant CA1 region, their retrieval did not conform to the temporal association templates. The lack of memory of time for predicting the 20-sec CS-US trace interval in the mutant mice provides the first direct evidence that the NMDA receptormediated plasticity is essential for producing ''what at when'' time memory traces that have been described in the hippocampus [21,36]. The working memory for tracking trace fear memory has been reported to involve cholinergic regulation in the hippocampus and other cortical areas [60,61]. Another neuromodulatory candidate mechanism might be dopamine [62][63][64]. This memory of time is likely an emergent result of the hippocampus interactions with the VTA, anterior cingulate cortex, amygdala, perirhinal cortex, and other higher cortical regions [21,[64][65][66][67][68][69].
In summary, by decoding CA1 activity patterns underlying Pavlovian fear conditioning, we have uncovered the real-time neural codes and temporal organizing patterns of trace and contextual fear memory engrams in the presence or absence of the NMDA receptors. Our study has, for the first time, described precisely the multiple actions of the NMDA receptors on regulating dynamic information contents and temporal patterns of real-time fear memory engrams in the brain.

Ethics statement
All procedures were conducted in accordance with the National Institutes of Health guidelines and with the approval of the Committee on Animal Care at the Georgia Regents University (GRU) and Banna Biomedical Research Institute (BBRI).

Production of inducible and region-specific NMDA receptor knockout mice
The inducible NMDA receptor knockout (KO) mice were generated as described previously and maintained on BCF hybrid background (B6 x CBAF1) [7] _ENREF_1. The KO mice are homozygous for the floxed-NR1 gene and heterozygous for the CaMKII-Cre transgene, the NR1-GFP transgene under control of the tet-O promoter, and the tetracycline transactivator (tTA) transgene, which is driven by the b-actin promoter and contains a floxed stop sequence (fNR1/fNR1, Cre/+, tTA/+, and NR1-GFP/+). The littermates lacking the Cre gene (fNR1/fNR1, tTA/+, NR1-GFP/+; or fNR1/+, tTA/+, NR1-GFP/+) were used as control mice. For genotyping, southern blot method was used to detect the floxed NR1 (fNR1) gene and the protocol is the same as described [70]. About 10 mg-purified tail DNA were digested by EcoR I, fractionated by electrophoresis on 0.7% agarose gels and transferred onto Zeta-probe GT membranes (BioRad). A 1.2 kb DNA fragment of 39 NR1 gene probe was labeled by a-32P-dCTP and hybridized to the GT membranes. PCR detection of the Cre, tTA, and NR1-GFP transgenes, approximately 0.5 to1 mg of mouse tail DNA was amplified in PT100 thermal cycler using the programs as follows: 1 minute, 94uC; 45 sec, 55uC; and 1 min, 55uC for 35 cycles. The primers for Cre detection is 59-AGA TGT TCG CGA TTA TC and 59-AGC TAC ACC AGA GAC GG; for tTA detection is 59-CAA TTA CGG GTC TAC CAT and 59-GGT TCC TTC ACA AAG ATC CTC; and for NR1-GFP associative traces triggered by the CS-paired foot shock in all seven trials from the five control mice (left panel) and five knockout mice (right panel). The blue rectangles indicate that the simple US traces produced by paired foot-shock, whereas the red rectangles indicate the US-to-CS associative traces were elicited by the US presentation. (G) The occurrences of CS-to-US associative traces triggered by the tone in all seven trials from the five control mice (left matrix) and five knockout mice (right matrix). The blue rectangles indicate that the simple CS traces produced by paired CS, whereas the red rectangles indicate the CS-to-US associative traces were elicited by the CS presentation during learning. doi:10.1371/journal.pone.0079454.g007 The PCR between fNR1-1 and fNR1-2 is ,160bp band for detecting the wild type NR1, and ,280bp band between fNR1-3 and fNR1-4 for detecting the fNR1 allele.
In our experiments, the inducible switch-off of the NMDA receptor function occurred 5 days before our recording experi-US and US-to-CS traces combined) in the control, but not mutant mice (n = 5 for each group). Student t-test, *p,0.05, error bars represent SEM. (I) The increase in the numbers of learning pattern reverberations was correlated with the increase in the amounts of immediate freezing during this learning phase in control mice (r = 0.84, p,0.05). (J) Lower level of freezing responses and learning pattern reverberations in knockout mice (r = 0.28, p . 0.05). doi:10.1371/journal.pone.0079454.g008  ments when we fed the inducible knockout mice with food pellets containing dox at 6 mg/g. This feeding protocol has been shown to disable NMDA receptor function within 3,5 days in the hippocampus and cortex in freely behaving mice [7]. Our previous study showed that the cre/loxP-mediated deletion occurred in about 97% of CA1 pyramidal cells and the averaged 60% of cortical principal neurons (CaMKII promoter active neurons). We maintained these mice on dox food throughout the experiments until the mice were sacrificed for histological verification of electrodes position.
In vivo recording and spike sorting 128-channel recording arrays were constructed as previously described in details [39] and used to record neural activity from hippocampal CA1 region in freely behaving mice (5-6 month old, males) [36][37][38][39][40]. Each mouse was implanted with two independently movable bundles of 32 steretrodes (64 channels on each side of the hippocampi) to bilateral hippocampi under deep anesthesia using 60 mg/kg ketamine (Bedford Laboratories, OH) and 4 mg/ kg Dormitor (Pfizer Animal Health, NY). The electrode bundles were positioned above the dorsal hippocampi (2.0 mm lateral to bregma and 2.3 posterior to bregma on both right and left sides). After the mice recovered from surgery, the electrodes were advanced slowly, over next five to ten days, in daily increments of about 0.07 mm until the tips of the electrodes reached the pyramidal layers of the hippocampal CA1 region.
The spike activity was recorded using Plexon Systems (Dallas) and then sorted using the MClust 3.3 program (http://redishlab. neuroscience.umn.edu/MClust/MClust.html). First, the recorded data were filed as Plexon system format (*.plx). Before spike sorting, the artifact waveforms were removed and the spike waveform minima were aligned using the Offline Sorter 2.8 software (http://www.plexon.com, Dallas, TX). The aligned data were then saved as files in Neuralynx system format (*.nst). After that, the MClust 3.3 program was used to isolate different spiking units. Only units with less than 0.5% of spike intervals within a 1 msec refractory period and with clear boundaries, judged by Lratio and Isolation distance calculation [71], were included in the present analysis (S. Fig. 1). Identification of different responsive and nonresponsive units from the same electrodes to tone and foot shocks are consistent with the notion that artifacts from external stimuli were removed and did not contaminate the recording results. Moreover, the observation that responsive units responding to foot-shock usually changed firing exceeding the stimulus duration (i.e. foot shock lasted only for 285-msec) further ruled out the contamination of electrical noise or artifacts.
Ten sets of recording data were obtained from five control and five mutant mice for the current analysis. Isolated units were identified. Further analysis is based on the sorted data; therefore, no electric artifacts were included during the shock events. To confirm the recording sites of the electrodes, 10 mA current was applied to each recording electrode for 5 seconds in order to mark the positions of the stereotrode bundles. Nissl staining confirmed the electrode positions. The stability of the ensemble recordings were judged by comparing waveforms and interspike intervals at the beginning, during, and after the experiments.
Only stably recorded units throughout the experiments are included in the present analysis. Recorded units were classified into putative pyramidal cells and fast-spiking putative interneurons. Putative pyramidal cells were defined as units with relatively wide waveforms (trough-to-peak length $300 ms) and lower average firing rates (,5 Hz), whereas putative interneurons were identified as units had relatively narrow waveforms (,250 ms) and higher firing rates (.5 Hz) [37]. Pyramidal cells are known to fire complex-spike bursts. To measure the characteristic burst activity, the proportion of inter-spike intervals that are shorter than 10 msec among all the inter-spike intervals was examined as the burst index.
By calculating waveforms widths and firing rates putative pyramidal cells and fast-spiking putative interneurons, we found that units' overall basic properties were largely similar between the control and knockout animals. Overall firing rate of putative pyramidal cells were about 1.4660.05 Hz during sleep state and 1.7360.06 Hz during awake in the control group, whereas the firing rates for putative pyramidal cells from the KO mice were 1.5160.05 Hz during sleep and 2.0360.08 Hz during awake. The spike widths of putative pyramidal units were 446.4763.58 in CT mice and were 441.4963.76 in KO mice. The burst index for the putative pyramidal cells was 13.4860.44% for the control and 10.1460.41% in KO mice during sleep (Student's t test p,0.001). During wakeful states, burst index was 13.1260.41% for the control and 10.7060.42% for the mutant mice (Student's t test p,0.001), suggesting a slight decrease in bursting upon the forebrain excitatory neuron-specific NMDA receptor knockout. On the other hand, there is no difference in the overall firing rate or waveform width of fast-spiking putative interneurons between the genotypes. The fast-spiking units from the control mice were about 13.4460.85 Hz during sleep and were 14.9260.92 Hz during awake, whereas the fast spiking units from the KO mice were at 13.6560.82 Hz during sleep and 16.0661.03 Hz during awake. The spike widths of putative interneurons were 194.3963.98 in CT mice and were 185.0063.97 in KO mice. This is consistent with the fact that the deletion of the NMDA receptor was restricted to excitatory neurons under the control of the CaMKII promoter [6,7].

Fear conditioning task
The fear conditioning chamber was a square chamber (10" 6 10" 6 15") with a 24-bar shock grid floor, and the trace recall chamber was a distinct semicircular shape chamber with a smooth and opaque floor. Freezing responses of animals in the chamber could be observed by experimenters and were videotaped [5]. Before training, the mouse was habituated in both chambers for five minutes per day, and three days in total. When the mouse was habituated in contextual chamber a tone was played 10 times. trace recall in the control group has the characteristics of exponential decay distribution, indicating it occurred in a 'bursting' manner. (G) Memory trace retrieved in the mutant mice did not show tight an exponential decay process. The inter-trace intervals of control and mutant mice formed different distribution (Two-sample Kolmogorov-Smirnov test, p = 5.3E-7). (H) Time distribution of US traces (include US simple trace and US-to-CS trace) in the control mice showed a distinct peak at about 22 second time-window at which a foot shock would be anticipated. This peak was absent in the mutant group. One way ANOVA test between two curve, p = 1.9E-6. (I) Time distribution of CS traces did not reveal any significant peaks in either control mice or knockout mice. (J) The occurrence of real-time shock memory traces around the 20-sec traced interval time window (18.5-24.5 sec after the offset of the recall tone) for all seven trials in five control mice. Red rectangles represent the occurrences of the correct foot-shock memory traces, whereas blue squares indicate the absence of foot-shock memory traces at the time point. (K) Lack of memory traces for time relationship in predicting foot shock event and its timing in five knockout mice. There is a significant difference in recalling memory of time relationship between CS and US between genotypes. Wilcoxon rank sum test, p,0.01. doi:10.1371/journal.pone.0079454.g010 This protocol facilitates the formation contextual memory as well as trace fear memory in rodents [36,57].
On the training day, the recording began with a 1-h pretraining sleep period in the home cage (a plastic tub where the mouse lived) and then followed by a 3-min pre-training exploration period in the shock chamber and then a 30-min pretraining period during which a 2-sec tone (85 dB continuous tonic sound at 2.8 kHz) was played ten times at random intervals. This in a mutant mouse in the 5-min contextual test. (G) Linear regression analysis shows that at group level, averaged freezing responses also correlated with their averaged numbers of total pattern retrievals (r 2 = 0.73, p,0.01). Each blue dot represents the data from a single control mouse and each red triangle represents the data from a single knockout mouse. (H) The total numbers of memory traces retrieved in the control and knockout mice (Wilcoxon rank sum test, **p,0.01; error bars represent SEM). The filled bar portion represents the associative memory traces retrieved during the contextual retention test. The knockout mice had few associative memory traces retrieved (Wilcoxon rank sum test, p,0.05). (I) Inter-memory tracetime interval analysis revealed that contextual recall in the control mice has the characteristics of exponential decay distribution. (J) Memory trace retrieved in the mutant mice did not show obvious exponential decay process temporal associativeness. There is significant difference between memory trace time distribution from control and mutant mice (two-sample Kolmogorov-Smirnov test, p = 1.2E-7). doi:10.1371/journal.pone.0079454.g011 Figure 12. Temporal map of CA1 fear memory neural codes underlying contextual memory recall. (A) The intertwined memory trace retrievals (defined as retrievals between individually distinct memory traces, see underlined brackets) and alternating retrieval dynamics (defined as retrievals between distinct simple traces, upper brackets) were prevalent in all five control mice during the 5-min contextual recall. Symbols: CS trace, blue triangle; US trace, red triangle; US-to-CS associative trace, red diamond; CS-to-US associative trace, blue diamond. The intertwined retrievals between different traces within 7 sec time window are underlined, whereas alternating retrievals between US and CS or CS and US are marked by upper brackets. (B) Degraded memory codes in five mutant mice during 5-min contextual retention tests. Only KO #1 had five associative traces during the recall, all other mutants failed to produce such traces. Intertwined and alternating retrieval patterns were greatly diminished. (C) The intertwined memory retrieval in the control group (blue plots) during contextual recall follows an exponential decay process, indicating that intertwined memory traces were recalled in temporal clusters. The mutant mice were significantly impaired in the temporal association between recalled memory traces (Red plots) (two-sample Kolmogorov-Smirnov test, p = 0.014). (D) Significant differences in the average occurrences of intertwined memory trace pairs between the control and knockout mice (Wilcoxon rank sum test, **p,0.01; error bars represent SD). (E) Numbers of various temporal structures of memory traces in single, doublet, triplet across the control and mutant mice during contextual recall. (F) The alternating retrieval for US simple trace followed with a simple CS memory trace or CS simple trace followed by a US memory trace in the control mice also exhibited exponential decay distribution, but not in the mutant mice (two-sample Kolmogorov-Smirnov test, p = 0.0042). Trace interval for 50% temporally associational recall in the control is 3.6 sec for control mice and 8.2 sec for mutant mice. (G) Differences in the alternating retrieval rates between the control and knockout mice (Wilcoxon rank sum test, **p,0.01; error bars represent SD). To appreciate temporal associativeness of memory recalled, memory patterns from a control (Mouse #4) and mutant mouse (Mouse #3) were converted into audio clips of ''Pavlovian memory symphony'' (Sound S1 and S2 correspondingly). doi:10.1371/journal.pone.0079454.g012 allowed us to profile the CA1 responses to naïve tone. The animal was then brought back to the home cage for a 30-min break before trace fear conditioning began. During trace conditioning, the conditioned stimulus (tone, 2-sec, 85 dB continuous tonic sound at 2.8 kHz) was used, but now as conditioned stimulus (CS) paired the unconditioned stimulus (a continuous 285-msec foot shock at 0.75 mA) with a 20-sec interval (after the off-set of the tone). This CS-US pairing was repeated for seven times, with 1-3 min random time intervals between each pairing. The mouse was then brought back to the home cage for 1-h post-training rest. The immediate freezing was calculated from the first 30 seconds after each shock.
After a one-hour rest in the home cages, memory retention/ recall tests in the contextual and traced fear paradigms were conducted. In the contextual memory recall test, the mouse was placed back to the shock chamber again for five minutes. After a short break in the home cages, the mouse was put into a novel chamber for the traced (tone-induced) fear memory recall test. The animal was allowed to explore freely for three minutes before the onset of the conditioned tone was played for 2 seconds (repeated seven times with random 1-3 min intervals between trace recall sessions). Freezing behaviors were scored in both retention tests. Throughout the procedures, animal behaviors were recorded simultaneously by videotaping which were synchronized with spike data collection. Freezing response, defined as absence of body movement except for respiration, was counted as the measurement of fear memory. We counted 'freezing' by observing the animal behaviors frame by frame based on the recorded video, the temporal resolution for the video is 30 frames per second and the absolute amount of time that the mice stayed in the continuous frozen state was summed and compared to the total amount of time.

Characterization of CA1 unit responses to fear conditioning stimuli
To determine whether a recorded unit is responsive to a stimulus, we used the time of stimulus delivery as time zeros to calculate a peri-event histogram using a 100 msec window. The neural activities during two seconds before stimulus were used as the baseline to determine a set of confidence intervals of 80%, 95% to determine neural responsiveness. First, the 80% confidence intervals were used to determine the duration of response while there was any significant peak (positive response) or trough (negative response) happened within one second after stimulus. Second, if the peak or trough was out of the 95% confident intervals, and the duration was longer than 0.5 sec, then it is considered a significant neuronal response to the stimulus. To facilitate comparison between units that exhibit different increases/decreases over baseline activities, we examined normalized neural response: Here, f resp,i represents the average firing rate during detected duration, f baseline,i is the average firing rate during baseline before The intertwined memory retrieval in the control group during tone-traced recall exhibited clear exponential decay distribution, but not in the knockout mice (two-sample Kolmogorov-Smirnov test, p = 0.0024). Time window for 50% temporally associational recall is 4.1 sec for control mice and 7.1 sec for mutant mice. (D) Significant reduction in intertwined memory retrievals during tone-traced recall in knockout mice as compared to the control group (Wilcoxon rank sum test, *p,0.05; error bars represent SD). (E) Average number of temporal structures of memory traces as single traces, doublet, triplet the control and mutant mice during trace recall. (F) The alternating retrieval between simple US trace and simple CS memory trace in the control mice exhibited an exponential decay distribution, suggesting that the retrieval occurred as clusters in time domain. Such temporal association is nearly absent in the mutant mice (two-sample Kolmogorov-Smirnov test, p = 0.022). Trace interval for 50% temporally associational recall in the control is 3.8 sec for control mice and 6.6 sec for mutant mice. (G) Differences in alternating retrieval rates during tone-traced recall between the control and knockout mice (Wilcoxon rank sum test, *p,0.05; error bars represent SD). doi:10.1371/journal.pone.0079454.g013 the i th stimulus, and f 0 is the average population activity during rest states. This transformation allows for uniform quantification of the significant changes in firing patterns for units with both lowand high baseline firing rates.

Hierarchical clustering
Hierarchical clustering methods were used to investigate the stimulus responses of the overall population of the simultaneously recorded CA1 units from both control and mutant mice. The procedure was described in our previously research [38,40,41]. This analysis was performed on the logarithm transformed neuronal response H~log 10 (1z R j j). R is an n|k matrix represents the neuronal responses of n units during k sampling points, k~X N i~1 n i corresponding to the all repetitions of N different stimulus while i th stimulus has n i repetitions respectively, and R j j denotes the absolute value of neural responses. To build hierarchy, first, a Euclidean distance matrix was constructed to measure pair-wise distances between different units. An agglomerative hierarchical cluster tree was created from the individual unit, which means in the beginning of the process, each unit to different stimulus is in a cluster of itself. Then the clusters with the shortest distance were sequential combined into a larger cluster and the pair-wise distances between clusters were updated. After the cluster merging completed, a categorical sorting was applied to facilitate the visualization, in which units were sorted by the number of stimuli they responded. After sorting, the units responded to the most stimuli were put on the top, and the nonresponsive units located at the bottom of the matrix.

Dimensionality-reduction based statistical pattern projection methods
We then used Multiple Discriminant Analysis (MDA) methods to project the neural responses of hundred units corresponding to different stimuli into different classes in subspaces [38,40,41]. To account for transient changes that may occur immediately after CS or US stimuli, we computed responses by using firing frequencies in two 250 msec temporal windows around the delivery of the stimuli. Responses during baseline period were characterized by computing the average firing rates during time intervals preceding the CS or US stimuli. Thus, population neural responses to the stimulus CS or US were normalized and projected by MDA to form the CS or US cluster in subspace, in the meantime responses during baseline were projected to form the Rest cluster.
The matrix of the responses from hundreds of simultaneously recorded neurons corresponding to each event was used to compute the between-class scatter matrix and within-class matrix: Here N is the number of categories (different types of stimuli), n i is the number of elements in each category (repetitions for each stimulus), m i is the mean vector for each category, m is the global mean vector from all categories and the symbol t indicates the transpose operator, r is the response of single neuron triggered by the i th stimulus, R i represents the set of population responses. Using these two matrices, the covariance matrix C can be obtained by C~S {1 W : S B . A set of at most N{1 discriminant projection vectors can be determined by computing the eigenvalue decomposition of the covariance matrix C.
For our data sets, the class covariance matrices S W were noninvertible, which is a direct consequence of data under-sampling, since the number of recorded units is much more than the number of repeated trials. In practice, the matrix S W can be rendered invertible using a regularization technique which changes each class covariance matrices based on the following formula: V 0 i~( 1{l)V i zlI, where V i is the covariance matrix for the i th category, l is a regularization parameter between 0 and 1, and I is the identity matrix. We determined the parameter lfor each data set based on the optimization procedure we developed previously [41]; each particular choice was determined by the particular distributions within each data set [36,41].
After computing the N{1 discriminant dimension, we can project the neural responses to low-dimensional encoding subspaces by transfer matrix, whose columns are the corresponding discriminant vectors and were sorted descend according to the eigenvalues. We can use the multivariate Gaussian distribution probability functions P(x)~1 (2p) N=2 V j j 1=2 exp ({(r{m) t V {1 (r{m)=2) to fit the projections for each category. In addition, we used a 20-msec sliding temporal window to monitor the evolution of the population state throughout the duration of the experiment and to identify the occurrences of patterns similar to the ones produced by CS or US stimuli [36,41].
For the tone-induced trace recall, it occurred in a novel environment which was different from the original training chamber. The global MDA analysis showed that tone-triggered responses during both learning and the recall test could still form a single CS cluster that was well separated from other events, indicating that major information carried by the tones during training and the recall was similar despite the contextual difference. However, the contextual difference can be further analyzed by using a second-step MDA analysis under which the CS ensemble trajectories elicited during recall can be separated from the CS ensemble trajectories during training [36].

Statistical criteria for distinguishing trajectory types in MDA subspaces
To determine the individual identity of ensemble trajectories during both learning and retention tests, two measures were developed: 1) 4s around cluster center was defined as the area close to the cluster in MDA subspaces, 2) since the distributions of the distance between the whole trajectories over the given trial and the cluster followed the Gaussian distribution, 2s of the mean distance was used as a boundary for the trace candidate. The first measure provided global control across multiple trials while the second measure took the trial variability into consideration. A trace candidate should be close to the CS or US cluster, and at the same time far away from the Rest cluster. Specially, to quantify that, first, the trajectory could be determined close to the CS or US cluster, if the trace reached within 4s from the cluster center, or if the distance between the trace candidate and the CS or US cluster reached less than 2s below the mean distance of the Gaussian distribution (Bonferroni correction). Second, to be far away from the Rest, the farthest turning points on the trace candidate should be at least 4s from the Rest cluster center, or the distance between the trace candidate and the Rest cluster should reach 2s above the mean distance to the Rest (Bonferroni correction). Finally, movement directions of ensemble trajectories to the CS/US cluster and to the Rest cluster were visually confirmed by rotating the clusters in 3-D plots. By applying the above criteria, trace candidate were determined. If a trace reached CS cluster and then moved to US cluster, it was determined as a CS-to-US associative trace. If a trace reached US cluster and then to CS cluster, it was determined as a US-to-CS associative trace.
To examine the contributions of neurons to transient ensemble patterns, we shuffled spike trains of the perspective units and then performed MDA projections for comparison with that of trajectories obtained from original spike trains prior to shuffle [40]. Perspective units could be the units responsive to memory traces or non-responsive units as a comparison. Surrogate spike trains of perspective units were generated by randomizing spike time, while maintaining their original overall firing rates. It is expected that the information contained in original spike train like ISIs would be abolished by the shuffling process. The same 250 msec-moving window method was then applied on surrogate spike trains for MDA projections while using Identity matrix obtained from original MDA training on the global datasets. To avoid the effects of numbers of shuffled units on MDA patterns, we matched the number of perspective units, specifically the numbers of nonresponsive units with the numbers of responsive units for spike train shuffling of each type of ensemble traces. Statistical criteria for defining trajectory types in MDA subspaces as described above were then applied.

Supporting Information
Sound S1 ''Pavlovian memory symphony'' of recall patterns of fear memory traces from the CA1 of the hippocampus of a wild-type mouse. The audio file reflects the memory patterns generated from the CA1 of the hippocampus of a wild-type mouse (mouse #4 in Figure 12A). To convert the recall patterns of fear memory traces into ''Pavlovian memory symphony'', four distinct types of fear memory traces retrieved during the five-minutes contextual recall session are represented by four different frequency tones (notes). Each tenor C, E, G and a soprano C note correspond to a simple CS trace, a simple US trace, a CS-to-US associative trace, and a US-to-CS associative trace, respectively. The audio clip of this 5-minute ''Pavlovian memory symphony'' is time compressed by a factor of four (into a 1.25-minute clip). (WAV) Sound S2 ''Pavlovian memory symphony'' of recall patterns of fear memory traces from the CA1 of the hippocampus of a mutant mouse. The audio file reflects the memory patterns generated from the CA1 of the hippocampus of a mutant mouse (knockout mouse #3 in Figure 12B). The conversion method is the same as Sound S1. (WAV)