Functional coupling networks inferred from prefrontal cortex activity show experience-related effective plasticity

Functional coupling networks are widely used to characterize collective patterns of activity in neural populations. Here, we ask whether functional couplings reflect the subtle changes, such as in physiological interactions, believed to take place during learning. We infer functional network models reproducing the spiking activity of simultaneously recorded neurons in prefrontal cortex (PFC) of rats, during the performance of a cross-modal rule shift task (task epoch), and during preceding and following sleep epochs. A large-scale study of the 96 recorded sessions allows us to detect, in about 20% of sessions, effective plasticity between the sleep epochs. These coupling modifications are correlated with the coupling values in the task epoch, and are supported by a small subset of the recorded neurons, which we identify by means of an automatized procedure. These potentiated groups increase their coativation frequency in the spiking data between the two sleep epochs, and, hence, participate to putative experience-related cell assemblies. Study of the reactivation dynamics of the potentiated groups suggests a possible connection with behavioral learning. Reactivation is largely driven by hippocampal ripple events when the rule is not yet learned, and may be much more autonomous, and presumably sustained by the potentiated PFC network, when learning is consolidated.

phases by Peyrache et al. (Peyrache et al., 2009). During the awake epoch the animal faced a learning task, where it had to find the rewarded arm in a Y-shaped maze; this arm was chosen according to a rule (left or right arm, or where the light is on or off) set by the operator. As soon as the rule was consistently learned it was changed. PCA-based analysis of the recorded activity showed that the activity of the learning phase was replayed during the subsequent sleep in some experimental sessions; this replay is at the basis of memory consolidation. Here, we reanalyze the same recordings, with a more sophisticated statistical approach than PCA, based on the inference of functional connectivity between the recorded cells. Our motivation is twofold. First, the use of functional couplings allows us to characterize task-related changes in the activity of the sleep epochs in a more quantitative way than with PCA, and, in addition, to identify more sessions showing replay. Second, we expect that our more precise and extended characterization of replay could make more precise the possible connection with behavioral learning sketched in the study by Peyrache et al. Our statistical approach relies on the inference of graphical models expressing the con-Graphical models: Probabilistic models of interacting variables, in which the distribution of each variable generally depends on a restricted number of other variables. ditional dependencies between the spiking events of the recorded cells through functional couplings. We make use of the maximum-entropy Ising model from statistical mechanics, Ising model: Mathematical model to describe systems of binary interacting components in statistical physics. whose parameters are tuned to reproduce the recorded firing frequencies and pairwise crosscorrelations (Schneidman, Berry, Segev, & Bialek, 2006). This inference approach has been tested on several multielectrode recordings of both in vitro (Barton & Cocco, 2013;Cocco, Leibler, & Monasson, 2009;Ferrari, Obuchi, & Mora, 2017;Schneidman et al., 2006) and in vivo (Barton & Cocco, 2013;Posani, Cocco, Jezek, & Monasson, 2017) neural activity. Working with the inferred couplings rather than considering directly the correlations in the data allows us to refine the analysis of the recording and to unveil modifications in the functional couplings between the two sleep epochs (effective positive or negative potentiation), which are consistent with the functional network derived from the learning epoch. Effective potentiation is supported by a subset of the recorded cells, which we identify by means of an automatized Effective potentiation: A quantity that sums up, over all the pairs of neurons, the increases of the coupling, between the sleep after the task and the sleep before the task, only if there is a correlated coupling increase between the task and the sleep before the task. procedure. Our findings are supported by a large-scale study of about 100 experimental sessions. Despite the variations from session to session, presumably because of the partial and random sampling of cells, we are able to identify in about 20% of sessions a potentiated group. We Potentiated group: A group of cells that sustain most of the potentiated couplings. This group is identified from the largest components of the top eigenvector of the potentiation matrix. then investigate in the data the change in the collective firing properties of the identified potentiated groups and find a strong increase in coactivation for such groups between the two sleep epochs. Hence, these identified potentiated groups are likely to belong to task-related cell assemblies. We then analyze how much the reactivation over time of the potentiated group is related to hippocampal inputs (ripple events), known to be important for memory consolidation. In sessions where the rule has not been learned yet, reactivation can be essentially explained as a fast response to hippocampal ripples. In some of the sessions where the rule was learned, reactivation shows a strong slow dynamical component, often unrelated to ripples, which presumably reflects the existence of a potentiated prefrontal cortex (PFC) synaptic network.

RESULTS
We have reanalyzed recordings of the activity of tens of neurons in the prefrontal cortex of five behaving rats (Peyrache et al., 2009); see Materials and Methods for more details. Each one of the 96 recording sessions is divided into three 30-minute epochs: a Task epoch in which the rat had to learn a cross-modal rule (go left, right, where the light is on, or off, in a Y-shaped Cross-modal rule shift: An introduction of new rules that can be based on spatial (go to the left or to the right) or visual (go where the light is on or off) cues. maze), which was changed as soon as the rat had learned it, and two Sleep epochs, one before (Sleep Pre) and one after (Sleep Post) the Task epoch. Through spike sorting one can identify the same neurons recorded in the different epochs (Sleep Pre, Task, Sleep Post) of a session; the number N of neurons reliably mapped in all three epochs varies from 3 to 56 depending on the session. No mapping could be established between different sessions.

Inference of Functional Coupling Networks
We briefly present the approach to model the distribution of activity of the N recorded neurons; see Methods section. The spiking times are binned within small time bins of duration Δt = 10 ms, as illustrated in Figure 1A; see Supporting Information (SI), Section II and Figure S2, for a discussion of the time-bin choice (Tavoni, Ferrari, Battaglia, Cocco, & Monasson, 2017). The activity configuration (σ 1 , σ 2 , ..., σ N ) is a snapshot of the neural activity, where σ i takes values one or zero depending on whether the i-th neuron is, respectively, active or inactive in the time bin. We define f i and f ij as the average values over time bins of, respectively, σ i and σ i σ j : f i represents the probability that neuron i is active in a time bin, and f ij denotes the joint probability that both cells i and j are active in a time bin.
We model the probability distribution of activity configurations as where Z ensures normalization of the distribution. P in Equation 1, called Ising model in statistical physics, is the least constrained (with maximum entropy), default probability distribution Δt; each neuron i is assigned the variable σ i = 1 or 0, depending on whether it is active in the time bin (left). Pairwise correlation indices CI ij are computed from this binned data, and define a network of correlations (middle). The network of statistical couplings J ij defining the Ising model distribution P, Equation 1, is generally sparser (right). Red and blue links correspond, respectively, to CI > 1, J > 0 and to CI < 1, J < 0; the widths are proportional to the absolute values. Links corresponding to CI or J smaller than one tenth of the maximal correlation index or coupling are not shown. (B) Average values and standard deviations of log CI ij over intervals of couplings 0.5 n − 0.25 ≤ J ij < 0.5 n + 0.25, with integer n, for all epochs and sessions. Note the large error bar in J = 0, corresponding to the very large number of pairs i, j carrying vanishing couplings; see Figure 1A.
reproducing these low-order spiking statistics (Schneidman et al., 2006). We look for the Ising model parameters {h i , J ij } such that f i and f ij match, respectively, the average values of σ i (for all neurons i) and σ i σ j (for all pairs of neurons i, j) over P. To do so we use the adaptive cluster expansion (ACE) inference algorithm (Barton and Cocco, 2013; Adaptive cluster expansion: Technique of inference of the Ising model parameters based on an iterative construction of clusters of strongly interacting variables. 2011,2012), which also gives access to the statistical uncertainties {δh i , δJ ij } over the inferred parameters (Methods). Parameters h i define effective local inputs that tune neuronal frequencies. Parameters J ij define the effective pairwise couplings between the cells ( Figure 1A): J ij different from zero expresses the presence of a conditional dependence between neurons i and j, not mediated by other neurons in the recorded population. The conditional average activity of neuron i given the other neuron activities {σ j }, with j = i, reads It is a logistic function of its total input, V i , equal to the sum of the other neuron activities σ j weighted by the couplings J ij , and of the local input h i .
Though effective couplings J ij are abstract quantities defined through Equations 1 and 2 they can be approximated by the logarithms of the correlation indices, CI ij = f ij /( f i f j ). As Correlation index (CI): Ratio of the joint probability that a pair of neurons are active in a time bin, and of the product of their individual spiking probabilities. For independent neurons, CI = 1.
shown in Methods, in the simple case of N = 2 recorded neurons only, J 12 and log CI 12 are equal. For N ≥ 3, log CI ij is only an approximation to J ij , and their difference quantifies the indirect contributions to pairwise correlations, mediated by other cells and not due to direct interactions; see Methods. Figure 1B shows that this approximation is good for most couplings in the recorded sessions, but deviations can be observed in particular for large and positive correlation indices; see Barton and Cocco (2013) for a discussion of the differences between log CI ij and J ij across various neural datasets. Note that the functional networks are sparse: A large fraction (about 75% over all epochs and sessions) of the couplings are regularized to zero by the inference procedure.
Once the coupling and local input parameters are inferred, we may sample the model distribution P through Monte Carlo simulations to check how the statistics of the data are reproduced by the model. The quality of the reproduction of the single-neuron and pairwise spiking probabilities in a time bin is shown in Figure 2A for the Task epoch of one particular session, which we call A. We can then use P to make predictions for higher-order moments, such as triplet firing probabilities and the probability of multiple-neuron firing in a time bin. Results are compared with the same quantities computed from the spiking data in Figure 2B. The quality of the inferred distribution P is then assessed through a cross-validation procedure: We divide the dataset into a train set (three fourths of the time bins) and a test set (one fourth of the bins). The good agreement between the values of observables in the train and test sets in Figures 2A and B confirms the absence of overfitting in our inference (see Figure S3 in Tavoni et al., 2017, for results on the other three out of four possible ways to define training and testing sets). We show in Figure 2C the probabilities of the 2 10 configurations of firing of one subset of 10 cells. The Ising model predictions are in much better agreement with the data than the independent-cell model, which reproduces the single-neuron spiking probabilities f i only. Taking into account pairwise correlations through the effective couplings J ij is therefore crucial to better fit the neural activity distribution.

Comparison of Functional Couplings Across Epochs Shows Learning-Related Potentiation
The distributions (over all sessions) of inferred coupling parameters are similar across epochs; see Figure 3A. In addition, little variation over the magnitudes of couplings is observed from The agreement between the spiking probabilities computed from the data and from the inferred Ising distribution shows that the inference procedure is accurate. Model distribution P was inferred from three fourths of the recorded data and tested on the same data (blue cross) and on the remaining one fourth of the recording (cross-validation, red squares). (B) Probabilities of firing for triplets ( f ijk , left panel) of neurons, and of k neurons to be simultaneously active in a time bin of duration Δt = 10 ms (right panel). The agreement between the data and model multiple-neuron firing probabilities (p(k)) is very good as long as p(k) times the number of time bins in the recording is > 1, that is, provided the recording time is sufficient to sample rare configurations of multiple neuron firing. Same crossvalidation procedure and symbols as in Figure 2A. (C) Probabilities of the 2 10 = 1, 024 activity configurations over a subset of 10 cells in the Task epoch of session A. Blue symbols show the scatter plot for the Ising distribution P (inferred from all recorded data), while cyan symbols correspond to the independent-cell model (with all couplings J ij = 0, and local inputs h i fitted to reproduce the single-neuron probabilities f i ). Similar plots are found for other subsets of 10 cells among the N = 37 recorded cells. session to session. As an illustration we report in Figure 3 the histograms of coupling parameters for session A. Because of the smaller number of data, the histograms are less smooth than the average distribution over all sessions, but span the same ranges of values for J. Similar results hold for local inputs; see SI, Figure S4 (Tavoni et al., 2017).
Despite the overall similarities between the coupling distributions across epochs, subtle patterns can be observed when tracking the changes in the couplings corresponding to the same pairs of cells across the different epochs of the same session. We partition the set of couplings in each epoch into three classes, according to their values J and statistical uncertainty δJ: Couplings reliably inferred as positive, that is, such that J/δJ > 3, define the  Figure 3B).
An important feature emerging from Figure 3B is the presence of task-related effective (positive) potentiation in the functional couplings. This effect is visible from the relative enrichment of [0 + +] (marked with a circle symbol) with respect to the null model (two-tail binomial test, p 10 −5 ; see Methods), while no such enrichment is found for classes [0 − +] and [00+]. In other words, we find that the fraction of pairs of neurons with close-to-zero couplings in Sleep Pre and positive couplings in both Task and Sleep Post is larger than what would be expected from the knowledge of the coupling classes in the Sleep epochs only. Task-related effective negative potentiation, corresponding to the enrichment of [0 − −] (star symbol), is also found, but with a weaker magnitude (p < 10 −5 ).
While the results above were obtained through averaging over all sessions, there are substantial variations in the fractions of pairs in the classes from session to session. We show in Figure 3C three examples, referred to as sessions A, Y, and B. For sessions A and B effective potentiation, represented in particular by class [0 + +] (circle symbols), is clearly visible (with, respectively, p < 10 −5 and p = 0.002). Session Y shows a strong effective negative potentiation, represented in particular by class [0 − −] (star symbol, p < 10 −5 ).
To characterize quantitatively experience-related changes in the functional couplings in each session, we introduce the following session-wide effective potentiation, measuring the amount of potentiation in the couplings from Sleep Pre to Sleep Post, coherently with their values in Task: Summation is restricted to pairs i, j of neurons, whose couplings are significantly different from zero in both Task and Sleep Post (same criterion |J ij |/δJ ij > 3 as for the classes above). The presence of the θ function, θ(u) = 1 if the argument u > 0 and 0 if u ≤ 0, restricts contributions to pairs, whose effective couplings increase from Sleep Pre to Task. In practice positive, respectively, negative contributions to Pot come mostly from classes [0 + +], respectively, The effective potentiations Pot of all 96 recorded experimental sessions are shown in Figure 3D. In comparison we show in Figure 3E the same quantity Pot computed after swapping, in each session, the Sleep Pre and Sleep Post couplings in Equation 3. No session is found to have a large effective potentiation after the swap, see Figure 3E. This simple control provides clear evidence for the fact that large values of Pot capture experience-related changes in the Sleep Post couplings. This empirical observation can be made more precise through the introduction of a null model, in which the correspondence between pairs of neurons across the epochs is removed by reshuffling the neuron indices, and the values of the couplings are randomly drawn from the distributions of Figure 3A. Red curves in Figure 3D show the average value of the effective potentiation within this null model, together with ± one standard deviation (Methods). In the control analysis of Figure 3E, where the Sleep Pre and Post data have been swapped, the effective potentiation is compatible with the expectations of the null model. Conversely, in Figure 3D, where the causal ordering Sleep Pre-Task-Sleep Post has been maintained, some sessions have large and positive Pot more than one standard deviation above the null model average.
A major source of variability in the data is the limited number of randomly sampled neurons. To assess the influence of sampling on Pot, we focus on one particular session (A) with large effective potentiation, and remove cells, one at a time, from the recording. Results are shown in Figure 4A. In most cases removal of one cell does not significantly affect the value of Pot. A substantial decrease is, however, observed for a small number of cells, indicated by the labels in Figure 4A. This result clearly shows that most contributions to Pot come from a restricted subset of the recorded neurons. How many of those relevant cells are or are not well sampled may explain, at least in part, the variability in potentiation values observed across sessions.

Groups of Neurons Supporting Effective Potentiation Are Replayed in the Sleep Epoch After Learning
The network of couplings supported by the group of neurons identified in Figure 4A are shown for the three epochs of session A in Figure 4B. The effective potentiation from Sleep Pre to Sleep Post and the strong similarity between the densely interconnected networks in Task and Sleep Post are clearly visible. For this session the potentiated couplings are not supported by independent, nonoverlapping pairs of neurons, but are densely interconnecting a restricted group of neurons (see Figure SI in Tavoni et al., 2017, for statistical validation). We emphasize that experience-related change in the correlational structure of Sleep Post is better seen with effective couplings than with pairwise correlations. For session A again, we show in Figure 4C Figure 4A (red circles in Figure 4C). This shows again that the changes experienced by the couplings J ij between the Sleep epochs are positively correlated to their values in Task. Conversely, the same comparison with the CI instead of the couplings J shows a much blurrier picture; see Figure 4C: The changes in CI between the Sleep epochs do not seem correlated with their values in Task.
While removing one cell at a time is an effective procedure to determine which neurons contribute most to Pot, it is computationally demanding. We have therefore developed a fast and fully automatized spectral procedure to directly identify in each session the group of neurons supporting the densest core of strongly potentiated couplings. Our procedure is based on taking the neurons with largest entries in the top eigenvector of the Pot matrix, whose elements are the contributions to Pot (Equation 3) of the pairs (i, j); see Methods for details. The top eigenvectors are shown in SI, Figure S9 (Tavoni et al., 2017), for a few sessions. We generally observe a few large entries, and many small ones. We have set a conservative value for the threshold to retain only large entries, see Methods and SI (Tavoni et al., 2017). For session A, our automatized procedure finds the five cells (the seven neurons identified in Figure 4A, but neurons 11 and 35) that support the couplings contributing most to Pot. This result extends to other sessions: We find a high correlation (0.76 ± 0.26 across all sessions) between the top eigenvectors of the Pot matrices and the "leave-one-out" potentiation vectors (shown in Figure 4A for session A), which require much more computational efforts.
We now show that the groups of neurons identifed with our automatized procedure across the 96 sessions really coactivate in the spiking data. To this aim we consider an extension of the pairwise correlation index CI ij to groups of more than two neurons. We define the assembly coactivation ratio (CoA) of a group G of neurons over the time scale τ through Coactivation ratio (CoA) of a group of neurons: Ratio of the probability that all the neurons in a group are active within a time scale, and of the product of the individual spiking probabilities. CoA is a multineuron generalization of CI.
where f (G) is the probability that all the neurons in the group are active within the time scale τ, and the denominator is the product of the individual spiking probabilities. For a group of independent cells the CoA is on average equal to unity. CoA is a very stringent measure of coactivation, as it takes into account only events in which all the neurons in the potentiated group are active. A less restrictive measure of the activity of the potentiated group will be studied in the next section.
We first show in Figure 5A the CoA of the five-cell potentiated group of session A above. This five-cell group is found to strongly coactivate in Task on a τ 20 − 40 ms time scale, and in Sleep Post on a similar time scale, τ 30 − 50 ms. The five-cell group does not coactivate more than expected by chance in Sleep Pre, which is compatible with the independent-cell hypothesis due to the low firing frequencies (Methods). This result shows that the potentiated group is replayed in Sleep Post. Interestingly, the coactivation of the potentiated group in Sleep Post is much stronger during non-REM-Sleep periods (non-REM), in which hippocampal sharp Non-REM sleep: Stages 1-3 of sleep characterized by absence of dreams and absence of rapid eye movements.
waves are known to be important for memory consolidation ( Figure 5A , right). In addition, the large CoAs of the potentiated group found in Task and Sleep Post are significantly higher than CoAs for random groups of five neurons (SI, Figure S15; Tavoni et al., 2017). Those findings suggest that the five-cell group is (part of) a cell assembly that is reinforced by experience.
For each recorded session we then measure the maximal values (over the time scale τ) reached by the CoA of the group supporting the effectively potentiated couplings in the Sleep Pre and Sleep Post epochs. The ratio of the maximal CoA in Sleep Post over the maximal CoA in Sleep Pre is a measure of the reinforcement of the coactivation between the neurons in the group across the two sleep epochs in a session. Figure 5B shows the scatter plot of the logarithms of the ratios of those two maximal CoA versus the values of the effective potentiations Pot of the groups (defined as the sums of contributions to Pot over the pairs of neurons in the groups) across the recorded sessions. A clear monotonic trend is observed, showing that our estimate of coupling potentiation is a good estimator of the existence of neural groups in the spiking data, which reinforce their coactivation in the sleep epoch following task-learning. We retain all the sessions in which the increase in the logarithms of maximal CoA across the Sleep epochs is larger than 1.5, and label them with letters A to Q; see Figure S7 (Tavoni et al., 2017) for locating sessions A-Q in the potentiation results of Figure 3D . The sizes of the potentiated groups in sessions A-Q range from two to seven cells. It is important to notice that variants of the identified potentiated groups with, say, one more or less cells, can also have large CoAs. Varying the threshold used in the spectral procedure allows us to explore these alternative groups in each session. Examples are provided in SI, Figures S13-S14 (Tavoni et al., 2017). While we have focused above on effective potentiation corresponding to an increase of the couplings across the Sleep epochs, effective negative potentiation, in which couplings get more negative in Sleep Post than in Sleep Pre, may be found in some sessions, such as Y ( Figure 3C), despite being weak on average ( Figure 3B). An analysis of the [0 − −] coupling class in session Y has permitted us to identify a cluster of three cells. We show in Figure 5A that this three-cell group is associated with a decrease of the CoA from Sleep Pre to Sleep Post in non-REM and at short time scales τ ∼ 20 − 40 ms. Equation 3 for Pot can be straightforwardly modified to define the effective negative potentiation; see SI, Section VIII and Figure S16 (Tavoni et al., 2017). In addition to Y we have identified another strongly negatively potentiated session, Z, and two sessions, C and I, showing both positive and negative potentiations across their coupling networks. As it is statistically hard to reliably estimate low CoA values, a systematic study of negative potentiation across all sessions is difficult, and would require longer recordings.

Dynamics of Reactivation: Effects of Hippocampal Ripples and Connection With Behavioral Learning
The previous analysis allowed us to identify effective potentiated groups that are strongly coactivated in Sleep Post in sessions A-Q. A fundamental issue is whether the reactivation of those potentiated groups is mostly triggered by hippocampal inputs (sharp-wave ripples, monitored Sharp-wave ripple complexes: Large activity bursts in the hippocampus composed of large amplitude local field potential deflections (sharp waves) associated with high-frequency field oscillations, called ripples, generated particularly during immobility and non-REM sleep.
in the experiments; see Methods) or reflects the internal dynamics of the PFC network, modified upon learning. To address this question, for each of the above sessions, we define the Reactivation of the potentiated group in time bin t, of the reactivation following a ripple event by a delay τ, normalized by the average reactivation, where the t m s are the times of the N r ripple events. Figure 6A shows the ripple-conditioned reactivations RR(τ) for the Sleep Post epochs of sessions A, B, C, which are representative of the variety of RR patterns found across all sessions. In sessions A and B, a marked reactivation peak is found at short time scales of tens of milliseconds. In addition, this "fast" peak is followed in session B by a long-lasting reactivation, decaying over a few seconds. We stress that fast versus slow reponses to ripples were not studied session by session in Peyrache et al. (2009), which reports only the average response over all sessions. No clear response of the reactivation to ripple events is found in session C on any time scale.
A complementary characterization of the reactivation dynamics, not directly related to ripples, is provided by the following normalized autocorrelation of the reactivation: Autocorrelation of the reactivation (AR): Autocorrelation of the reactivation of the potentiated group.
where the brackets · denote the average over all time bins t. AR captures dynamical scales, irrespective of their origins (inputs from the hippocampus or internal dynamics of the PFC network). The behaviors of AR(τ) are reported for sessions A, B, C again in Figure 6A. For all three sessions we observe a large peak in the autocorrelation at τ = 0, expressing the tendency of neurons in the potentiated group to fire together and decaying over a few tens of in Methods, of the amplitude of the fast (τ = 0) and the slow (averaged over 0.5 s < τ < 1.5 s) components to the ripple-conditioned reactivation for sessions A to Q. Each session is represented by one, two, or three circles of increasing diameters, depending on the number of rules. For each rule the color of the corresponding circle indicates whether the learning-point and rule-changing criteria were reached. Session K is not shown because of the very small number of detected ripples. In session D, both learning-point and rule-changing criteria were met but the rule was not changed in the experiment as it should have been; see SI, Section IX (Tavoni et al., 2017). Average and standard deviation of the fast decay times over the eight sessions with largest fast responses to ripples are τ RR− f ast = 230 ± 180 ms. Average value and standard deviation of the slow decay time over the three sessions (I, B, Q) with a significant slow response to ripples are τ RR−slow = 4 ± 3 s. (C) Amplitudes of the fast (in τ = 0) and the slow (average over 0.5 s < τ < 1.5 s) components in AR(τ) for sessions A to Q. Same color code for learning behavior as in panel B. Fast and slow decay time constants, fitted over all sessions are, respectively, τ AR− f ast = 56 ± 24 ms and τ AR−slow = 3 ± 1.5 s.
milliseconds. This "fast" peak is followed by a slow component decaying over few seconds.
Remarkably, in session C, for which no reactivation following ripples was detected, AR is stronger than for sessions A and B.
The results above were extended through a systematic analysis of sessions A-Q showing strong effective potentiation. The magnitudes of the fast (0 < τ < 50 ms) and slow (0.5 s < τ < 1.5 s) components of the ripple-conditioned reactivation are reported in Figure 6B; results are expressed in terms of Z scores with respect to a null model defined from the behavior of RR(τ) at negative delays τ < 0; see Methods. In Figure 6C we plot the amplitudes of the fast and slow components to the autocorrelation AR(τ).
We now attempt to relate the characterization of sessions in terms of RR and AR to the learning behavior of the rat during the Task epoch. The experimental protocol is described in Methods and in SI, Section IX (Tavoni et al., 2017); see also Peyrache et al. (2009). The rule was changed during the session if the rat had done 10 consecutive correct trials, or made only one error out of 12 successive trials. In most of the selected sessions A-Q, this rule-changing criterion was never reached. In some of the selected sessions, the rule-changing criterion was reached once, and a second rule was introduced; in one session (L), the criterion was met again after the second rule was set, and a third rule was introduced; see SI, Section IX (Tavoni et al., 2017).
In addition, during offline data analysis (Peyrache et al., 2009), a learning point was defined, based again on the success rate of the rat but according to a less stringent criterion: the rule was said to be learned if the rat had done three consecutive correct trials, and afterwards had a success rate larger than 80% over the remaining trials (up to the end of the session end or up to the change of the rule).
Whether the learning-point and the rule-changing criteria were reached defined different learning scenarios, which are illustrated by sessions A, B, C studied above ( Figure 6A). In session A, a new rule was introduced at the beginning of the session, and was neither changed nor learned (learning point not reached) throughout the session. In session B, the rule was the same as in the previous session, and the learning point was reached by the end of the session but the rule was not changed. In session C, the rule was also the same as in the previous session and was changed in the middle of the session since the rat had fulfilled the rule-changing criterion; the learning point was not reached for the second rule by the end of session C. These three sessions can be informally seen as three successive levels of behavioral learning: rule not yet learned in session A (none of the two criteria is reached), intermediary learning stage (success rate has reached intermediary values between the learning and rule-changing points) in session B, and first rule definitively learned (both criteria met) for session C.
Despite the limited number (16) of selected sessions, these three learning scenarios seem to be in correspondence with general features of RR and AR presented in Figures 6B and 6C (the complete list of sessions with detailed results can be found in SI, Section IX, Tavoni Figure 6C (J and L show neither large RR nor large AR). In four of those sessions (C, P, D, K) the rule-changing criterion was reached. The exception is session O, which, however, presents some atypical features of long-lasting learning; see SI, Section IX (Tavoni et al., 2017).
In summary, these observations point to the following possible connection between behavioral learning and the features of the ripple-conditioned reactivation and the auto-reactivation of the neural group identified from the spiking data: (a) When neither the learning-point nor the rule-changing criteria are met, RR shows a fast component, and no slow component, while AR is weak, suggesting that the reactivation of the potentiated group is only due to the ripple inputs; (b) When the learning point is reached but the rule-changing criterion is not met, RR shows both fast and slow components while AR is weak, suggesting the presence of weak, underlying synaptic potentiations able to sustain the activity of the group after its initiation by ripples; (c) When both learning-point and rule-changing criteria are reached, there is no significant RR, but AR is strong; this is compatible with the fact that those sessions have been recorded at a stage of advanced learning, when reactivation of consolidated task-related cell assemblies might have become independent of hippocampal inputs. We emphasize that these observations are highly empirical and speculative. More statistics would definitively be needed to firmly establish and confirm those rules.

DISCUSSION
In the present work we have focused on experience-related modifications to the functional connectivity matrix (between tens of recorded neurons) in the prefrontal cortex (PFC) of rats (Peyrache et al., 2009. Functional connectivity was defined through the introduction of a graphical (Ising) model, accounting for the statistical dependencies betweens spiking events of the neurons in the recorded population. Comparing the functional networks in the two Sleep epochs before and after learning we found, in a substantial fraction of the sessions under investigation, some changes correlated to the functional connectivity during the learning epoch itself. In most of these sessions, we found that a fraction of the couplings became effectively potentiated, and that those couplings were supported by a limited subset of the neuronal cells (the so-called potentiated group). In other words, a group of cells became much more strongly interconnected in the Sleep epoch after learning than before. We have directly verified on the spiking data that neurons in the identified potentiated groups coactivated much more in the Sleep epoch after than before the learning epoch, which is reminiscent of the notion of cell assembly introduced by Hebb (1949) as the basic unit of neural computation and memory. Study of the reactivation dynamics of the potentiated groups allowed us to separate effects due to hippocampal inputs (ripples) or to a putative PFC network, in connection with learning.

Patterns of Changes in Functional Couplings Between Epochs and Potentiated Groups
As a general result we have found that functional couplings define sparse interaction network in each single epoch, the class [000] concentrating most of the pairs. In addition there is an overall correlation between the amplitude of couplings across the different epochs, including Task, which can be seen from the relatively large fractions of pairs in classes While most of our analysis was focused on experience-related modifications to the functional connectivity, other mechanisms may take place. Tonni and Cirelli (2006) have suggested that, during specific phases of sleep (Genzel, Kroes, Dresler, & Battaglia, 2014), small synaptic interactions are erased, a phenomenon called homeostasis. The overall similarity between the distribution of inferred couplings in the two Sleep epochs, with many zero couplings, is somewhat in agreement with this hypothesis. However, it is difficult to distinguish between small couplings and couplings strictly equal to zero; homeostatic changes, if any, would likely fall in the most populated [000] class.
A potential bias in our analysis is that neither the non-REM nor the REM periods have equal REM sleep: Phase of sleep characterized by rapid eye movements and the propensity of the sleeper to dream. durations in the Sleep Pre and Post epochs of the same session ( Figure S1, Tavoni et al., 2017). We have checked the robustness of our estimates for the session-wide effective potentiation, Figure 3D, and for the group potentiation, Figure 5B, under random uniform subsampling of the recorded data in which the duration of the non-REM and REM periods were matched between the two Sleep epochs; see SI, Figure S11 (Tavoni et al., 2017).
The changes in the inferred networks of functional couplings between the Sleep epochs, correlated with the coupling network in the Task epoch, are supported by a subset of the recorded neurons. The identification of these potentiated groups of neurons was done through an automatized spectral analysis of the Pot matrix. We have shown that the groups of potentiated neurons strongly coactivate in the Sleep epoch posterior to learning, and are therefore part of a replayed experience-related cell assembly. It is clear, however, that the notion of potentiated group should be intended in a statistical sense. Slight variations in the composition of the group, such as adding or removing one specific neuron, are associated with large coactivation, as shown in SI, Figures S12 and S14 (Tavoni et al., 2017).
It is a remarkable and somewhat counterintuitive fact that the network of couplings inferred from pairwise coactivation on short time scales, Δt = 10 ms, suffices to predict coactivation patterns between n neurons on longer time scales, τ n × Δt ms. However, even in the case of coactivation events, the repeated spiking of neurons in short bursts generates a sequence of pairwise coactivation events ( Figure 1A), and the coactivated groups appear as strongly interconnected. Robustness of predictions against the global temporal scale of the cell assembly and the activation ordering is an important advantage of the Ising model, because cell assembly can be played and replayed at different time scales (and in direct and reverse orderings).

Functional Couplings: Consequence of Common Inputs or Real Interactions?
As first discussed in the works of Gerstein and collaborators (Aertsen et al., 1989;Gerstein & Perkel, 1969), functional couplings can reflect either synaptic connections or the presence of a transient common input coactivating two or more neurons. Within the Hebbian paradigm, coactivation is a prerequisite to learning, favoring synaptic potentiation, such as through LTP (Buzsáki, 2015). In our data, common inputs could be identified in the transient sharp waves from hippocampus to the prefrontal cortex, occurring preferentially during non-REM sleep. Sharp-wave ripples have been experimentally demonstrated to be essential for memory consolidation (Buzsáki, 2015;Genzel et al., 2014;Girardeau, Benchenane, Wiener, Buzsáki, & Zugaro, 2009). Synaptic potentiation in the cortex has been suggested to take place in the immediately following stage, thanks to spindle oscillations contributing to the shutdown of the transmission from the hippocampus to the prefrontal cortex. The calcium influx taking place during spindle oscillation could facilitate synaptic potentiation between cells in the replayed assemblies (Genzel et al., 2014;Siapas & Wilson, 1998).
An important issue is whether changes in the functional couplings between the Sleep epochs reflect such common inputs, necessary for learning, or "real" plasticity in the synaptic interactions. In their earlier work Peyrache and collaborators estimated, based on principal component analysis, the average reactivation over all sessions, and showed that it occurred within a 2-second time window centered around the sharp-wave event (see Figure 5B in Peyrache et al., 2009). In our study we computed the reactivation of the more precise potentiated groups defined by the Ising model, for each one of the 16 selected sessions. In sessions with a clear response to ripple this response starts within 250 ms from the ripple event, and thus shows a finer temporal resolution. Moreover, despite the difficulty in identifying experience-related cell assemblies, the complexity of rule-changing scenarios, and the experimental limitations in the recordings, different scenarios seem to emerge, depending on the learning stage.
For the four sessions in which the rat has not learned the rule, the ripple-conditioned reactivation RR(τ) of the identified cell assembly decays after a short delay of the order of τ ∼ 200 ms, comparable to the typical duration of sharp-wave ripples. Hence, the strongly interconnected effective network we identified in Sleep Post ( Figure 4D) seemingly mostly accounts for correlations produced by neural coactivation under common hippocampal inputs.
Furthermore, in several sessions with strong effective potentiation and in which the rule (or two rules in session Q) is (are) learned, towards the session end hippocampal ripples evoke a persistent reactivation, lasting several seconds after the ripple event. This effect may signal the existence of an established synaptic potentiation, able to reverberate the activity seconds after the input is over. In other sessions in which the rule was definitively learned and changed, no significant reactivation of the potentiated group following the ripple events was found; however, reactivation in those sessions showed a large autocorrelation, decaying over seconds.
A tentative interpretation of these findings is the following. For sessions in which the rule has not yet been learned, the large coactivation of the experience-related group, as evidenced in Figure 5A, seems to be largely supported by the inputs coming from the hippocampus during the sharp wave ripples, known to be crucial for memory consolidation. Conversely, our finding suggests that when the rule has been learned, reactivation occurs over long time scales with two possible mechanisms: slowly decaying persistence of ripple-induced activity, which is found in sessions where the rule has "just" been learned, or second-long replay periods, unrelated to ripple events, which takes mostly place in sessions where the rule has been learned and changed (more stringent criterion). The presence of long-lasting reactivation suggests the existence of a potentiated synaptic network connecting the PFC neurons. This putative network could be either evoked by ripples or subject to spontaneous excitations, depending on its maturity. A possible interpretation of the absence of ripple-induced activity for sessions in which the rule has been definitively learned is that memory has been consolidated, so ripples are not needed for the passage from short to long-term memory any longer. As stressed in the Results section, these interpretations are highly speculative, as the limited number of selected sessions and the variety of behaviors of animals during those sessions did not allow us to draw any solid statistics.
The above results are also consistent with the finding that when the rule has not been learned, the CoA of the potentiated group is often substantially larger in non-REM than in REM; see for example CoA of session A in Figure 5A. In sessions in which the rule has been learned, on the contrary, there is a still smaller but significative coactivation of the potentiated group also in REM periods, as happens in session B on the time scale of τ ∼ 50 ms (see also Figure S13 in Tavoni et al., 2017, for other sessions). Ripple events are indeed more frequent in non-REM sleep. Finally the presence of a CoA larger than one in Sleep Pre, which further increases in Sleep Post, is also often present when the rule is not new (see Figure S13, Tavoni et al., 2017). New experiments and more data for the different learning and rule shift protocols would be important to confirm these findings.
Note that in a related work (Tavoni, Cocco, & Monasson, 2016), we have shown how the simulation of the inferred Ising model in the presence of an external stimulation allowed us to reveal groups of coactivating cells and to characterize their statistical variations. This external input, introduced as a mathematical tool to scan rare coactivation events in the Ising model of the population activity, mimicks, in a very crude way, hippocampal inputs to the PFC and allows for the investigation of collective activation on time scales larger than the elementary time bin used for inference. The groups of neurons in the experience-related cell assemblies identified by our method and the one in Tavoni et al. (2016) coincide for the two sessions (A and D) common to our set of 16 selected sessions and to our previous work.

Comparison Between Inference Procedures for Ising Model
The couplings defining the graphical model considered in this paper are an extension of the correlation indices first used to quantify functional connectivity (Aertsen et al., 1989;Fujisawa et al., 2008;Schwindel et al., 2014). Informally speaking, couplings can be viewed as a sparse set of "direct" correlations among the population of recorded cells. Even for sessions with few recorded cells, the network of couplings is much sparser than its correlation counterpart (see Figure 1A). We have resorted here to the graphical Ising model, which is the maximum entropy model reproducing the one-and two-cell firing statistics. In this probabilistic framework, different methods exist to infer the couplings parameters.
A competitive inference technique is the standard Boltzmann machine learning algorithm (Hinton & Sejnowski, 1986), broadly used in the analysis of retinal (Schneidman et al., 2006) and cortical (Marre, El Boustani, Frégnac, & Destexhe, 2009) multielectrode recordings. This inference procedure is slow in its naive version but can become efficient with good initial guess of the coupling parameters (Barton and Cocco, 2013), upon replacement of the gradient descent for the minimization of the cross-entropy with an approximate version of the Newton method (Ferrari, 2016), or thanks to improvement specific to the sparse activity of neural population (Broderick, Dudik, Tkacik, Schapire, & Bialek, 2007). Two other promising methods to fit Ising models from data are the pseudolikelihood approach (Aurell & Ekeberg, 2012) and minimal probability flow (Sohl-Dickstein, Battaglino, & DeWeese, 2011). Both approaches use all the data, and not only the first and second moments of the neural activity, to avoid computing the normalization constant Z in the distribution P; see Equation 1. In particular, the minimal probability flow method has been recently applied to multilayer restricted Boltzmann machine to model explicitly the different columns in cortical data (Köster, Sohl-Dickstein, Gray, & Olshausen, 2014).
Here we have used the adaptive cluster expansion (ACE; see Methods and Barton & Cocco, 2013;, which has been shown to accurately reproduce interaction parameters for synthetic data. The Ising distribution inferred with ACE also reproduces the statistics of retinal or hippocampal recordings of the activity of tens to hundreds of recorded neurons, including high-order moments (Barton & Cocco, 2013;Posani et al., 2017). More recently this technique was generalized to the case of nonbinary but multiple-categorial variable, or Potts model, to model coevolution in protein sequences (Barton, De Leonardis, Coucke, & Cocco, 2016). Our inference approach is very fast on these neural data, taking some seconds on a personal computer to infer the input and functional connectivity parameters. It would be possible to use it to identify cell assemblies online. Combined with optogenetics techniques (Reutsky-Gefen et al., 2013) this would open exciting perspectives on the manipulation of cortical cell assemblies in a controlled way.
The Ising model captures the statistics of snapshots of the activities, and, as such, defines symmetric functional couplings J ij = J ji . It can therefore not be used to study the ordering in the dynamical activation of the neurons. Other functional-connectivity-based inference approaches, such as the generalized linear (Pillow et al., 2008;Truccolo, Eden, Fellows, Donoghue, & Brown, 2005), kinetic Ising (Roudi & Hertz, 2011), and integrate and fire (Koyama & Paninski, 2010; models are designed to infer nonsymmetric connectivity matrices from the temporal sequence of spiking events in the neuronal population. In Tavoni et al. (2016), we have inferred the couplings on the same cortical dataset with the generalized linear model and found that they are essentially symmetric, and strongly correlated with their Ising counterparts. One possible explanation is that cell assemblies in the prefrontal cortex may also code for instrinsically nontemporal aspects of the task to be learned, in agreement with the findings of Peyrache et al. (2010).

Comparison With Existing Procedures to Identify Cell Assemblies
Many of the currently available methods to detect and characterize the replay of neural groups or cell assemblies rely on the knowledge of how the neural activity correlates with sensory or internal inputs. For example, place cells in the hippocampus are known to encode location in space, and replay of place-cell assemblies representing behaviorally meaningful trajectories can be determined with template-matching techniques. More precisely, the ordered activation sequences of place cells observed during salient moments, such as sharp-wave ripple events, during sleep or wakefulness, is matched with the sequences of place cells determined by the templates observed during locomotion (Brown, Frank, Tang, Quirk, & Wilson, 1998;Carr et al., 2011;Diba & Buzsáki, 2007;Foster & Wilson, 2006;Lee & Wilson, 2002;Pfeiffer & Foster, 2013). Similarly, in sensory systems cell assemblies can be detected and characterized by studying the neuronal population response to specific stimuli, easily reproducible in experimental settings. An example is provided by the analysis of neural activity patterns following specific sounds in the auditory cortex (Bathellier, Ushakova, & Rumpel, 2012). However, those approaches are not easily applicable to all the regions of the brain. In the prefrontal cortex, for instance, neurons may not be activated in a well-defined temporal order, predictable from the knowledge of external stimuli. Cell assemblies might respond to internal cognitive states, or to a combination of extrinsic covariates and internal states, which are very difficult to determine and control experimentally.
In this context, principal component analysis (PCA) has been used as a way to build approximate templates from the correlational structure of data (top principal components), and to detect reactivation, or replay, of those templates. Though PCA was applied successfully to detect replay Peyrache et al., 2009Peyrache et al., , 2010, it lacks any probabilistic framework and the interpretation of the large entries of the top components is difficult, even with the use of clustering procedures, such as the assembly vector estimation Lopes-dos-Santos et al., 2013. Our analysis significantly extends the PCA of Peyrache et al. (2009Peyrache et al. ( , 2010, as it identifies the neurons participating to replay-related assemblies in a detailed way. Let us stress that the whole approach for computing functional connectivity and identifying cell assemblies is fully automatized, and requires spiking data only. While our approach identifies a single cell assembly that contributes the most to functional coupling potentiation, it could be easily extended to the case of more assemblies; see Methods and Fortunato (2010). Such an extension could be useful for analyzing bigger recordings in the future.
A community detection technique for cell assembly identification, exploiting the Markov stability method, was recently introduced in (Billeh et al., 2014). The method consists of finding a stable partition on a correlation graph and was tested on hippocampal and retinal data. This graph, unlike the coupling networks we infer here, is defined heuristically and does not disentangle direct (giving rise to coupling) from indirect (mediated through other neurons) correlations.

Description of Experiments
Experimental methods were described in detail in Benchenane et al. (2010) and Peyrache et al. (2009) and are summarized in the following. Five Long-Evans male rats weighing 250-300 g at arrival were implanted with tetrodes in the prelimbic (PL) subdivision of the medial prefrontal cortex, and in the intermediate-ventral hippocampus. PL tetrodes were used for recording of single units: signals were band-pass filtered between 600 and 6,000 Hz, and spikes were detected whenever the filtered signal exceeded a manually set threshold. The resulting waveform (1.3 ms long) was fed into an automated spike sorting algorithm (KlustaKwik; Kadir, Goodman, & Harris, 2014). Hippocampal tetrodes were only used for local field potentials, for the detection of theta rhythms and sharp waves. Non-REM was automatically detected, based on power in the cortical delta band (1-4 Hz), hippocampal theta (5-10 Hz), cortical spindles (10-20 Hz), and speed of head motion, by means of a clustering algorithm. A quality check on sleep epochs ensures the absence of systematic biases in Sleep Pre with respect to Sleep Post; see Figure S1 (Tavoni et al., 2017).
The rats performed an attentional set shift task on a Y-maze, which is known to require the function of the medial prefrontal cortex (mPFC) in rats (Birrell & Brown, 2000). Each recording session consisted of a 20-to 30-minute sleep or rest epoch (Sleep Pre epoch) in which the rat remained undisturbed in a padded flowerpot placed on the central platform of the maze, a Task epoch, in which the rat performed the behavioral task described below for 20-40 min, and by a second sleep or rest epoch (Sleep Post epoch; same situation as in Sleep Pre) of 20-30 min. The whole recordings in the Task epoch and in the Sleep phases (both during REM and non-REM periods) were used for our inference.
Rats started each trial in the same arm (the departure arm). One of the two other (choice) arms was illuminated at random (pseudorandom schedule: runs of more than four consecutive trials with the same illuminated arm were avoided, as were repeated bouts of imposed alternation between the two arms). After that, the central platform was lowered, allowing the rat to access the choice arms. Only one of the choice arms was rewarded, according to one of four contingency rules. Two contingency rules were spatially guided (always go to the right arm, or to the left arm); the other two were cue guided (go to the illuminated arm, or to the dark arm). The rule that was employed at any given moment in time was not signaled to the rat in any way, so that the animal had to learn the rule by trial and error. Once the rat reached a criterion of 10 consecutive correct trials, or only one error out of 12 trials, the rule was changed with no further cue warning to the rat. Rule changes were extradimensional, that is, from a spatially guided rule to a cue-guided rule, and vice versa. All five rats learned in a consolidated way the right and light rules (at least 10 consecutive correct trials), whereas only two learned in a consolidated way the left task and one to go where the light is off.

Inference of Ising Model Parameters
Inference procedure.
We have inferred the Ising model parameters with the adaptive cluster expansion (ACE) algorithm (Barton & Cocco, 2013;Barton et al., 2016;, available from https://github.com/johnbarton/ACE. ACE computes an approximation for the (cross-)entropy of the Ising model reproducing the following data: where Z[{h i , J ij }] is a function of the Ising parameters normalizing the distribution P in Equation 1. The ACE procedure recursively builds clusters Γ = (i 1 , i 2 , . . . , i K ) of neurons with increasing sizes K, and selects those whose contributions ΔS Ising (Γ) to the cross-entropy exceed a threshold Θ (in absolute value) (Barton & Cocco, 2013;Cocco & Monasson, 2012). The total cross-entropy is approximated by the sum of the cluster entropies over the set of selected clusters, S Ising ∑ selected Γ ΔS Ising (Γ). The optimal value of the threshold, Θ, is chosen so that the Ising model (with parameters realizing the minimum in Equation 7) reproduces the experimental low-order statistics within the expected sampling accuracy. Choosing smaller values for Θ would require more computational efforts and would overfit the data.

Regularization and statistical error bars.
To regularize the minimization problem in Equation 7 we (a) remove from the datasets neurons spiking less than 10 times in each epoch, that is, such that f i < 10/B, where B is the number of time bins, that is, the duration of the recording divided by the time-bin duration Δt; (b) add to the right-hand side of Equation 7 a term γ ∑ k<l J 2 kl penalizing large couplings. In practice we choose γ = 0.2/B. This choice ensures that the coupling between pairs of neurons that are never active together in any time bin and that would, in principle, be assigned an infinitely large negative value become comparable to the couplings between pairs of neurons with a few common spiking event across the recording.
To quantify the statistical errors on the inferred parameters, we evaluate the matrix of the second derivatives of the entropy S Ising with respect to the Ising parameters {h i , J kl }, also called Fisher information matrix. The squared statistical error bars (δh i ) 2 and (δJ kl ) 2 on, respectively, the inferred local inputs h i and the couplings J kl are given by the diagonal elements of the inverse matrix of the Fisher information matrix, divided by the number of time bins, B (Cocco & Monasson, 2012). Remark that, with the addition of the regularization term, the Fisher information matrix is definite positive, and, therefore, can be inverted.
(i, j) of neurons, Pot ij ≡ p ij,ij,ij , cf. Equation 9. This contribution can be seen as an entry of an N × N-dimensional matrix. This matrix is sparse, square symmetric, and has large entries for neurons (i, j) supporting strongly potentiated couplings. The top eigenvector, v = {v i }, of the matrix is localized over few neurons i, which strongly contribute to Pot (SI, Figure S9, Tavoni et al., 2017). We define the potentiated group as the set of neurons i with components v i larger than a threshold value c (ranging between zero and one for normalized v). In all the experimental sessions considered here, this simple spectral analysis gives at most one large and connected neural group. Spectral graph theory offers efficient methods for dealing with more complex data structures on larger datasets, including more than one largely interconnected group (Fortunato, 2010).
The correlation between the changes in log CoA from Sleep Pre to Post, Δ logCoA, and the potentiation of the group across all sessions (shown in Figure 5B for c = 0.22) varies with c. The best correlations are found in practice in the range 0.15 < c < 0.35, with p values ranging between 10 −30 and 10 −35 (SI, Figure S10, Tavoni et al., 2017). We have arbitrarily set c = 0.22 in between these two limits. This value is also adequate if one imposes in addition that potentiated groups should include at least three cells (SI, Figure S10, Tavoni et al., 2017).

Statistical Significance of the Coactivation Ratio (CoA)
To assess the statistical validity of the CoA defined in Equation 4 for a group G of neurons, we compute the error bar on CoA, shown in Figure 5A. Assuming a Poisson distribution for the coactivation events, the standard deviation of the CoA is estimated to be CoA(τ)/ N G (τ), where N G (τ) is the number of coactivation events for the cells in G over the time scale τ.
Note that simultaneous-firing events (contributing to f (G)) are unlikely to be found, and the CoA is likely to be zero, if the duration of the recording is small, for example, compared with T min = τ/ ∏ i∈G f i (τ). This happens for the five-cell potentiated group of session A for time scales τ ≤ 40 ms in the Task epoch, and for all the values of τ considered in Sleep Pre and Post in Figure 5A.

Analysis of Ripple-Conditioned Reactivation
We define a null model for the average response to ripples after a delay τ, RR(τ) in Equation 5, as follows. For each session, we compute the average value, RR 0 , and the standard deviation, δRR 0 , of RR(τ) over the −3 s < τ < −0.5 s range. The range of delays is sufficiently negative to exclude any inaccuracy in the determination of the ripple times t m . The values of RR 0 and RR 0 ± δRR 0 are shown in Figure 6A (left panels) for sessions A, B, and C. The Z score of ripple-conditioned reactivation for positive delay τ is defined through The value of the Z score in τ = 0 and its average over the 0.5 s < τ < 1.5 s interval are used to estimate the amplitudes of the, respectively, fast and slow components in RR; see Figure 6B.

ACKNOWLEDGMENTS
We thank Georges Debregeas for useful comments about the manuscript. This work was funded by the [EU-]FP7 FET OPEN project Enlightenment 284801.