Temporally delayed linear modelling (TDLM) measures replay in both animals and humans
Abstract
There are rich structures in offtask neural activity which are hypothesized to reflect fundamental computations across a broad spectrum of cognitive functions. Here, we develop an analysis toolkit – temporal delayed linear modelling (TDLM) – for analysing such activity. TDLM is a domaingeneral method for finding neural sequences that respect a prespecified transition graph. It combines nonlinear classification and linear temporal modelling to test for statistical regularities in sequences of taskrelated reactivations. TDLM is developed on the noninvasive neuroimaging data and is designed to take care of confounds and maximize sequence detection ability. Notably, as a linear framework, TDLM can be easily extended, without loss of generality, to capture rodent replay in electrophysiology, including in continuous spaces, as well as addressing secondorder inference questions, for example, its temporal and spatial varying pattern. We hope TDLM will advance a deeper understanding of neural computation and promote a richer convergence between animal and human neuroscience.
Introduction
Human neuroscience has made remarkable progress in detailing the relationship between the representations of different stimuli during task performance (Haxby et al., 2014; Kriegeskorte et al., 2008; Barron et al., 2016). At the same time, it is increasingly clear that resting, offtask, brain activities are structurally rich (Smith et al., 2009; Tavor et al., 2016). An ability to study spontaneous activity with respect to taskrelated representation is important for understanding cognitive process beyond current sensation (Higgins et al., 2021). However, unlike the case for taskbased activity, little attention has been given to techniques that can measure representational content of resting brain activity in humans.
Unlike human neuroscience, representational content of resting activity is studied extensively in animal neuroscience. One seminal example is ‘hippocampal replay’ (Wilson and McNaughton, 1994; Skaggs and McNaughton, 1996; Louie and Wilson, 2001; Lee and Wilson, 2002): during sleep, and quiet wakefulness, place cells in the hippocampus (that signal selflocation during periods of activity) spontaneously recapitulate old, and explore new, trajectories through an environment. These internally generated sequences are hypothesized to reflect a fundamental feature of neural computation across tasks (Foster, 2017; Ólafsdóttir et al., 2018; Pfeiffer, 2020; Carr et al., 2011; Lisman et al., 2017). Numerous methods have been proposed to analyse hippocampal replay (Davidson et al., 2009; Grosmark and Buzsáki, 2016; Maboudi et al., 2018). However, they are not domain general in that they are designed to be most suited for specific needs, such as particular task design, data modality, or research question (van der Meer et al., 2020; Tingley and Peyrache, 2020). Most commonly, these methods apply to invasive electrophysiology signals, aiming to detect sequences in a linear track during spatial navigation task (Tingley and Peyrache, 2020). As a result, they cannot be directly adapted for analysing human resting activity collected using noninvasive neuroimaging techniques. Furthermore, in rodent neuroscience, it is nontrivial to adapt these algorithms to even small changes in tasks (such as 2D foraging). This may be a limiting factor in taking replay analyses to more interesting and complex tasks, such as complex mazes (Rosenberg et al., 2021).
Here, we introduce temporal delayed linear modelling (TDLM), a domaingeneral analysis toolkit, for characterizing temporal structure of internally generated neural representations in rodent electrophysiology as well as human neuroimaging data. TDLM is inspired by existing replay detection methods (Skaggs and McNaughton, 1996; Davidson et al., 2009; Grosmark and Buzsáki, 2016), especially those analysis of population of replay events (Grosmark and Buzsáki, 2016). It is developed based on the general linear modelling (GLM) framework and can therefore easily accommodate testing of ‘secondorder’ statistical questions (van der Meer et al., 2020), such as whether there is more forward than reverse replay, or is replay strength changing over time, or differs between behavioural conditions. This type of question is ubiquitous in cognitive studies, but is typically addressed ad hoc in other replay detection methods (van der Meer et al., 2020). In TDLM, such questions are treated naturally as linear contrasts of effects in a GLM.
Here, we show TDLM is suited to measure the average amount of replay across many events (i.e. replay strength) in linear modelling. This makes it applicable to both rodent electrophysiology and human neuroimaging. Applying TDLM on noninvasive neuroimaging data in humans, we, and others, have shown it is possible to measure the average sequenceness (propensity for replay) in spontaneous neural representations (Wimmer et al., 2020; Nour et al., 2021; Liu et al., 2019; Liu et al., 2021a). The results resemble key characteristics found in rodent hippocampal replay and inform key computational principles of human cognition (Liu et al., 2019).
In the following sections, we first introduce the logic and mechanics of TDLM in detail, followed by a careful treatment of its statistical inference procedure. We test TDLM in both simulation (see section ‘Simulating MEG data’) and real human MEG/EEG data (see section ‘Human replay dataset’). We then turn to rodent electrophysiology and compare TDLM to existing rodent replay methods, extending TDLM to work on a continuous state space. Lastly, using our approach we reanalyse rodent electrophysiology data from Ólafsdóttir et al., 2016 (see section ‘Rodent replay dataset’) and show what TDLM can offer uniquely compared to existing methods in rodent replay analysis.
To summarize, TDLM is a general, and flexible, tool for measuring neural sequences. It facilitates crossspecies investigations by linking largescale measurements in humans to singleneuron measurements in nonhuman species. It provides a powerful tool for revealing abstract cognitive processes that extend beyond sensory representation, potentially opening doors for new avenues of research in cognitive science.
Results
Temporal delayed linear modelling
Overview of TDLM
Our primary goal is to test for temporal structure of neural representations in humans. However, to facilitate crossspecies investigation (Barron et al., 2021), we also want to extend this method to enable measurement of sequences in other species (e.g. rodents). Consequently, this sequence detection method has to be domain general. We chose to measure sequences in a decoded state space (e.g. posterior estimated locations in rodents [Grosmark and Buzsáki, 2016] or time course of taskrelated reactivations in humans [Liu et al., 2019]) as this makes results from different data types comparable.
Ideally, a general sequence detection method should (1) uncover structural regularities in the reactivation of neural activity, (2) control for confounds that are not of interest, and (3) test whether this regularity conforms to a hypothesized structure. To achieve these goals, we developed the method under a GLM framework, and henceforth refer to it as temporal delayed linear modelling, that is, TDLM. Although TDLM works on a decoded state space, it still needs to take account of confounds inherent in the data where the state space is decoded from. This is a main focus of TDLM.
The starting point of TDLM is a set of n time series, each corresponding to a decoded neural representation of a task variable of interest. This is what we call the state space, X, with dimension of time by states. These time series could themselves be obtained in several ways, described in detail in a later section (‘Getting the states’). The aim of TDLM is to identify taskrelated regularities in sequences of these representations.
Consider, for example, a task in which participants have been trained such that n = 4 distinct sensory objects (A, B, C, and D) appear in a consistent order $:A\to B\to C\to D$ (Figure 1a, b). If we are interested in replay of this sequence during subsequent resting periods (Figure 1c, d), we might want to ask statistical questions of the following form: 'Does the existence of a neural representation of A, at time T, predict the occurrence of a representation of B at time T+ $\u2206t$?' and similarly for $B\to C$ and $C\to D$.
In TDLM, we ask such questions using a twostep process. First, for each of the n^{2} possible pairs of variables X_{i} and X_{j}, we find the linear relation between the X_{i} time series and the $\u2206t$shifted X_{j} time series. These n^{2} relations comprise an empirical transition matrix, describing how likely each variable is to be succeeded at a lag of $\u2206t$ by each other variable (Figure 1e). Second, we linearly relate this empirical transition matrix with a taskrelated transition matrix of interest (Figure 1f). This produces a single number that characterizes the extent to which the neural data follow the transition matrix of interest, which we call ‘sequenceness’. Finally, we repeat this entire process for all $\u2206t$ of interest, yielding a measure of sequenceness at each possible lag between variables and submit this for statistical inference (Figure 1g).
Note that, for now, this approach decomposes a sequence (such as $A\to B\to C\to D$) into its constituent transitions and sums the evidence for each transition. Therefore, it does not require that the transitions themselves are sequential: $A\to B$ and $B\to C$ could occur at unrelated times, so long as the withinpair time lag was the same. For interested readers, we address how to strengthen the inference by looking explicitly for longer sequences in Appendix 1: Multistep sequences.
Constructing the empirical transition matrix
In order to find evidence for statetostate transitions at some time lag $\u2206t$, we could regress a timelagged copy of one state, ${X}_{j}$, onto another, ${X}_{i}$ (omitting residual term ε in all linear equations):
Instead, TDLM chooses to include all states in the same regression model for important reasons, detailed in section ‘Moving to multiple linear regression’:
In this equation, the values of all states ${X}_{k}$ at time t are used in a single multilinear model to predict the value of the single state ${X}_{j}$ at time $t+\u2206t$.
The regression described in Equation 2 is performed once for each ${X}_{j}$, and these equations can be arranged in matrix form as follows:
Each row of X is a time point, and each of the n columns is a state. $X\left(\u2206t\right)$ is the same matrix as X, but with the rows shifted forwards in time by $\u2206t$. ${\beta}_{ij}$ is an estimate of the influence of ${X}_{i}\left(t\right)$ on ${X}_{j}\left(t+\u2206t\right)$. $\beta $ is an $n\times n$ matrix of weights, which we call the empirical transition matrix.
To obtain $\beta $, we invert Equation 3 by ordinary least squares regression:
This inversion can be repeated for each possible time lag ( $\u2206t=1,2,3,\dots $), resulting in a separate empirical transition matrix β at every time lag. We call this step the firstlevel sequence analysis.
Testing the hypothesized transitions
The firstlevel sequence analysis assesses evidence for all possible statetostate transitions. The next step in TDLM is to test for the strength of a particular hypothesized sequence, specified as a transition matrix,T. Therefore, we construct another GLM which relates T to the empirical transition matrix, β. We call this step the secondlevel sequence analysis:
As noted above, $\beta $ is the empirical transition matrix obtained from firststage GLM. It has dimension of $n$ by $n$, where $n$ is the number of states. Each entry in $\beta$ reflects the unique contribution of state i to state j at given time lag. Effectively, the above equation models this empirical transition matrix $\beta $ as a weighted sum of prespecified template matrices, $T}_{r$. Thus, $r$ is the number of regressors included in the secondstage GLM, and each scalar valued $Z\left(r\right)$ is the weight assigned to the $r$ th template matrix. Put in other words, $T}_{r$ constitutes the regressors in the design matrix, each of which has a prespecified template structure, for example, $T}_{auto$, $T}_{const$, $T}_{F$, and $T}_{B$ (Figure 1h).
$T}_{F$ and $T}_{B$ are the transpose of each other (e.g. red and blue entries in Figure 1b), indicating transitions of interest in forward and backward direction, respectively. In 1D physical space, $T}_{F$ and $T}_{B$ would be the shifted diagonal matrices with ones on the first upper and lower off diagonals. $T}_{const$ is a constant matrix that models away the average of all transitions, ensuring that any weight on $T}_{F$ and $T}_{B$ reflects its unique contribution. $T}_{auto$ is the identity matrix. $T}_{auto$ models selftransitions to control for autocorrelation (equivalently, we could simply omit the diagonal elements from the regression).
Z is the weights of the secondlevel regression, which is a vector with dimension of 1 by $r$. Each entry in Z reflects the strength of the hypothesized transitions in the empirical ones, that is, sequenceness. Repeating the regression of Equation 5 at each time lag ($\mathrm{\Delta}t=1,2,3,\dots$) results in time courses of the sequenceness as a function of time lag (e.g. the solid black line in Figure 1f). $Z}_{F$, $Z}_{B$ are the forward and backward sequenceness, respectively (e.g. red and blue lines in Figure 1g).
In many cases, Z_{F} and Z_{B} will be the final outputs of a TDLM analysis. However, it may sometimes also be useful to consider the quantity:
$D$ contrasts forward and backward sequences to give a measure that is positive if sequences occur mainly in a forward direction and negative if sequences occur mainly in a backward direction. This may be advantageous if, for example, $Z}_{F$ and $Z}_{B$ are correlated across subjects (due to factors such as subject engagement and measurement sensitivity). In this case, $D$ may have lower crosssubject variance than either $Z}_{F$ or $Z}_{B$ as the subtraction removes common variance.
Finally, to test for statistical significance, TDLM relies on a nonparametric permutationbased method. The null distribution is constructed by randomly shuffling the identities of the n states many times and recalculating the secondlevel analysis for each shuffle (Figure 1g). This approach allows us to reject the null hypothesis that there is no relationship between the empirical transition matrix and the taskdefined transition of interest. Note that there are many incorrect ways to perform permutations, which permute factors that are not exchangeable under the null hypothesis and therefore lead to false positives. We examine some of these later with simulations and real data. In some cases, it may be desirable to test slightly different hypotheses by using a different set of permutations; this is discussed later.
If the time lag $\mathrm{\Delta}t$ at which neural sequences exist is not known a priori, then we must correct for multiple comparisons over all tested lags. This can be achieved by using the maximum Z_{F} across all tested lags as the test statistic (see details in section 'Correcting for multiple comparisons'). If we choose this test statistic, then any values of Z_{F} exceeding the 95th percentile of the null distribution can be treated as significant at $\alpha =0.05$ (e.g. the grey dotted line in Figure 1g).
TDLM steps in detail
Getting the states
As described above, the input to TDLM is a set of time series of decoded neural representations or states. Here, we provide different examples of specific state spaces (X, with dimension of time by states) that we have worked with using TDLM.
States as sensory stimuli
The simplest case, perhaps, is to define a state in terms of a neural representation of sensory stimuli, for example, face, house. To obtain their associated neural representation, we present these stimuli in a randomized order at the start of a task and record wholebrain neural activity using a noninvasive neuroimaging method, for example, Magnetoencephalography (MEG) or Electroencephalography (EEG). We then train a model to map the pattern of recorded neural activity to the presented image (Figure 1—figure supplement 1). This could be any of the multitude of available decoding models. For simplicity, we used a logistic regression model throughout.
The states here are defined in terms of stimulievoked neural activity. The classifiers are trained at 200 ms poststimulus onset. For example, the stimuli are faces, buildings, body parts, and objects. Source localizing the evoked neural activity, we found that the activation patterns of stimuli in MEG signal are consistent with those reported in fMRI literature. For faces, activation peaked in a region roughly consistent with the fusiform face area (FFA) as well as the occipital face area (OFA). Activation for building stimuli was located between a parahippocampal place area (PPA) and retrosplenial cortex (RSC), a region also known to respond to scene and building stimuli. Activation for body part stimuli localized to a region consistent with the extrastriate body area (EBA). Activation for objects was in a region consistent with an objectassociated lateral occipital cortex (LOC) as well as an anterior temporal lobe (ATL) cluster that may relate to conceptual processing of objects. Those maps are thresholded to display localized peaks. The full unthresholded maps can be found at https://neurovault.org/collections/6088/. This is adapted from Wimmer et al., 2020.
In MEG/EEG, neural activity is recorded by multiple sensor arrays on the scalp. The sensor arrays record wholebrain neural activity at millisecond temporal resolution. To avoid a potential selection bias (given the sequence is expressed in time), we choose wholebrain sensor activity at a single time point (i.e. spatial feature) as the training data fed into classifier training.
Ideally, we would like to select a time point where the neural activity can be most truthfully read out. This can be indexed as the time point that gives the peak decoding accuracy. If the state is defined by the sensory features of stimuli, we can use a classical leaveoneout crossvalidation scheme to determine the ability of classifiers to generalize to unseen data of the same stimulus type (decoding accuracy) at each time point (see Appendix 2 for its algorithm box). In essence, this crossvalidation scheme is asking whether the classifier trained on this sensory feature can be used to classify the unseen data of the same stimuli (Figure 2a, b).
After we have identified the peak time point based on the crossvalidation, we can train the decoding models based on the multivariate sensor data at this given time.
Specifically, let us denote the training data, $M$, with dimension of number of observations, $b$, by number of sensors, $s$. The labels, Y, have dimension of $b$ by 1. The aim here is to obtain the classifier weights, W, so that $\mathrm{Y}\approx \sigma \left(\mathrm{M}\mathrm{W}\right)$. $\sigma$ is the logistic sigmoid function.
Normally we apply L1 regularization on the inference of weights (we will detail the reasons in section ‘Regularization’):
Next, we translate the data at testing time (e.g. during rest), R, from sensor space to the decoded state space:
where R is the testing data, with dimension of time by sensors, and X is the decoded state space, with dimension of time by states.
States as abstractions
As well as sequences of sensory representations, it is possible to search for replay of more abstract neural representations. Such abstractions might be associated with the presented image (e.g. mammal vs. fish), in which case analysis can proceed as above by swapping categories for images (Wimmer et al., 2020). A more subtle example, however, is where the abstraction pertains to the sequence or graph itself. In space, for example, grid cells encode spatial coordinates in a fashion that abstracts over the sensory particularities of any one environment, and therefore can be reused across environments (Fyhn et al., 2007). In human studies, similar representations have been observed for the location in a sequence (Liu et al., 2019; Dehaene et al., 2015). For example, different sequences have shared representations for their second items (Figure 2). These representations also replay (Liu et al., 2019). However, to measure this replay we need to train decoders for these abstract representations. This poses a conundrum as it is not possible to elicit the abstract representations in the absence of the concrete examples (i.e., the sensory stimuli). Care is required to ensure that the decoders are sensitive to the abstract code rather than the sensory representations (see Appendix 2 for algorithm box of selecting time point for training abstract code). Useful strategies include training classifiers to generalize across stimulus sets and ensuring the classifiers are orthogonal to sensory representations (Figure 2—figure supplement 1; details in Liu et al., 2019). One way that excludes the possibility of sensory contamination is if the structural representations can be shown to sequence before the subjects have ever seen their sensory correlates (Liu et al., 2019).
TDLM can also be used iteratively to ask questions about the ordering of different types of replay events (Figure 2d). This can provide for powerful inferences about the temporal organization of replay, such as the temporal structure between sequences, or the repeating pattern of the same sequence. This more sophisticated use of TDLM merits its own consideration and is discussed in Appendix 3: Sequences of sequences.
Controlling confounds and maximizing sensitivity in sequence detection
Here, we motivate the key features of TDLM.
Temporal correlations
In standard linear methods, unmodelled temporal autocorrelation can inflate statistical scores. Techniques such as autoregressive noise modelling are commonplace to mitigate these effects (Colclough et al., 2015; Deodatis and Shinozuka, 1988). However, autocorrelation is a particular burden for analysis of sequences, where it interacts with correlations between the decoded neural variables.
To see this, consider a situation where we are testing for the sequence $X}_{i}\to {X}_{j$. TDLM is interested in the correlation between $X}_{i$ and lagged $X}_{j$ (see Equation 1). But if the $X}_{i$ and $X}_{j$ time series contain autocorrelations and are also correlated with one another, then ${X}_{i}\left(t\right)$ will necessarily be correlated with ${X}_{j}\left(t+\mathrm{\Delta}t\right)$. Hence, the analysis will spuriously report sequences.
Correlations between states are commonplace. Consider representations of visual stimuli decoded from neuroimaging data. If these states are decoded using an nway classifier (forcing exactly one state to be decoded at each moment), then the n states will be anticorrelated by construction. On the other hand, if states are each classified against a null state corresponding to the absence of stimuli, then the n states will typically be positively correlated with one another.
Notably, in our case, because these autocorrelations are identical between forward and backward sequences, one approach for removing them is to compute the difference measure described above ($D={Z}_{F}{Z}_{B}$). This works well as shown in KurthNelson et al., 2016. However, a downside is it prevents us from measuring forward and backward sequences independently. The remainder of this section considers alternative approaches that allow for independent measurement of forward and backward sequences.
Moving to multiple linear regression
The spurious correlations above are induced because ${X}_{j}\left(t\right)$ mediates a linear relationship between ${X}_{i}\left(t\right)$ and ${X}_{j}\left(t+\mathrm{\Delta}t\right)$. Hence, if we knew ${X}_{j}\left(t\right),$ we can solve the problem by simply controlling for it in a linear regression, as in Granger causality (Eichler, 2007):
Unfortunately, we do not have access to the ground truth of $X$ because these variables have been decoded noisily from brain activity. Any error in ${X}_{j}\left(t\right)$ but not ${X}_{i}\left(t\right)$ will mean that the control for autocorrelation is imperfect, leading to spurious weight on $\beta}_{ij$, and therefore spurious inference of sequences.
This problem cannot be solved without a perfect estimate of X, but it can be systematically reduced until negligible. It turns out that the necessary strategy is simple. We do not know ground truth ${X}_{j}\left(t\right)$, but what if we knew a subspace that included estimated ${X}_{j}\left(t\right)$? If we control for that whole subspace, we would be on safe ground. We can get closer and closer to this by including further coregressors that are themselves correlated with estimated ${X}_{j}\left(t\right)$ with different errors from ground truth ${X}_{j}\left(t\right)$. The most straightforward approach is to include the other states of $X\left(t\right)$, each of which has different errors, leading to the multiple linear regression of Equation 2.
Figure 3a shows this method applied to the same simulated data whose correlation structure induces false positives in the simple linear regression of Equation 1, and by the same logic, so too in crosscorrelation. This is why previous studies based on a crosscorrelation (Eldar et al., 2018; KurthNelson et al., 2016) cannot look for sequenceness in forward and backward directions separately, but have to rely on their asymmetry. The multiple regression accounts for the correlation structure of the data and allows correct inference to be made. Unlike the simple subtraction method proposed above (Figure 3a, left panel), the multiple regression permits separate inference on forward and backward sequences.
Oscillations and long timescale autocorrelations
Equation 2 performs multiple regression, regressing each ${X}_{j}\left(t+\mathrm{\Delta}t\right)$ onto each ${X}_{i}\left(t\right)$ whilst controlling for all other state estimates at time t. This method works well when spurious relationships between ${X}_{i}\left(t\right)$ and ${X}_{j}\left(t+\mathrm{\Delta}t\right)$ are mediated by the subspace spanned by the other estimated states at time t (in particular, ${X}_{j}\left(t\right)$). One situation in which this assumption might be challenged is when replay is superimposed on a large neural oscillation. For example, during rest (which is often the time of interest in replay analysis), MEG and EEG data often express a large alpha rhythm, at around 10 Hz.
If all states experience the same oscillation at the same phase, the approach correctly controls false positives. The oscillation induces a spurious correlation between ${X}_{i}\left(t\right)$ and ${X}_{j}\left(t+\mathrm{\Delta}t\right)$, but, as before, this spurious correlation is mediated by ${X}_{j}\left(t\right)$.
However, this logic fails when states experience oscillations at different phases. This scenario may occur, for example, if we assume there are travelling waves in cortex (Lubenov and Siapas, 2009; Wilson et al., 2001) because different sensors will experience the wave at different times and different states have different contributions from each sensor. MEG sensors can be seen as measures of local field potential on the scalp, which contain background neural oscillations. In humans, this is dominantly alpha during rest.
In this case, ${X}_{i}\left(t\right)$ predicts ${X}_{j}\left(t+\mathrm{\Delta}t\right)$ over and above ${X}_{j}\left(t\right)$. To see this, consider the situation where $\mathrm{\Delta}t$ is $\frac{1}{4}\tau$ (where $\tau$ is the oscillatory period) and the phase shift between ${X}_{i}\left(t\right)$ and ${X}_{j}\left(t\right)$ is pi/2. Now every peak in ${X}_{j}\left(t+\mathrm{\Delta}t\right)$ corresponds to a peak in ${X}_{i}\left(t\right)$ but a zero of ${X}_{j}\left(t\right)$.
To combat this, we can include phaseshifted versions/more time points of $X\left(t\right)$. If dominant background oscillation is at alpha frequency (e.g. 10 Hz), neural activity at time T would be correlated with activity at time T + $\tau$. We can control for that by including $X\left(t+\tau \right)$, as well as $X\left(t\right)$, in the GLM (Figure 3b). Here, $\tau$ = 100 ms if assuming the frequency is 10 Hz. Applying this method to the real MEG data during rest, we see much diminished 10 Hz oscillation in sequence detection during rest (Liu et al., 2019).
Spatial correlations
As mentioned above, correlations between decoded variables commonly occur. The simplest type of decoding model is a binary classifier that maps brain activity to one of two states. These states will, by definition, be perfectly anticorrelated. Conversely, if separate classifiers are trained to distinguish each state’s representation from baseline (‘null’) brain data, then the states will often be positively correlated with each other.
Unfortunately, positive or negative correlations between states reduce the sensitivity of sequence detection because it is difficult to distinguish between states within the sequence: collinearity impairs estimation of β in Equation 2. In Figure 3c, we show in simulation that the ability to detect real sequences goes down as the absolute value of a spatial correlation goes up. We took the absolute value here because the direction of correlation is not important, only the magnitude of the correlation matters.
Ideally, the state decoding models should be as independent as possible. We have suggested the approach of training models to discriminate one state against a mixture of other states and null data (Liu et al., 2019; KurthNelson et al., 2016). This mixture ratio can be adjusted. Adding more null data causes the states to be positively correlated with each other, while less null data leads to negative correlation. We adjust the ratio to bring the correlation between states as close to zero as possible. In Figure 3d, we show in simulation the ensuing benefit for sequence detection. An alternative method is penalizing covariance between states in the classifier’s cost function (Weinberger et al., 1988).
Regularization
A key parameter in training highdimensional decoding models is the degree of regularization. In sequence analysis, we are often interested in spontaneous reactivation of state representations, as in replay. However, our decoding models are typically trained on taskevoked data because this is the only time at which we know the ground truth of what is being represented. This poses a challenge insofar as the models best suited for decoding evoked activity at training may not be well suited for decoding spontaneous activity at subsequent tests. Regularizing the classifier (e.g. with an L1 norm) is a common technique for increasing outofsample generalization (to avoid overfitting). Here, it has the added potential benefit of reducing spatial correlation between classifier weights.
During classifier training, we can impose L1 or L2 constraints over the inference of classifier coefficients, $W.$ This amounts to finding the coefficients, $W$, that maximize the likelihood of the data observations under the constraint imposed by the regularization term. L1 regularization can be phrased as maximizing the likelihood, subject to a regularization penalty on the L1 norm of the coefficient vector:
L2 regression can be viewed as a problem of maximizing the likelihood, subject to a regularization penalty on the L2 norm of the coefficient vector:
where M is the task data, with dimension of number of observations, $b$, by number of sensors, $s$. Y is the label of observations, a vector with dimension of $b$ by 1. $\mathrm{P}\left(\mathrm{Y}\mathrm{M},\mathrm{W}\right)=\sigma \left(\mathrm{M}\mathrm{W}\right)$, and $\sigma$ is the logistic sigmoid function.
We simulate data with varying numbers of true sequences at 40 ms lag and find that the beta estimate of sequence strength at 40 ms positively relates to the number of sequences. We also find that L1 weight regularization is able to detect sequences more robustly than L2 regularization, while L2 performs no better than an unregularized model (Figure 3e). The L1 models also have much lower spatial correlation, consistent with L1 achieving better sequence detection by reducing the covariances between classifiers.
In addition to minimizing spatial correlations, as discussed above, it can also be shown that L1induced sparsity encodes weaker assumptions about background noise distributions into the classifiers as compared to L2 regularization (Higgins, 2019). This might be of special interest to researchers who want to measure replay during sleep. Here, the use of sparse classifiers is helpful as background noise distributions are likely to differ more substantially from the (awake state) training data.
Statistical inference
So far, we have shown how to quantify sequences in representational dynamics. An essential final step is assessing the statistical reliability of these quantities.
All the tests described in this section evaluate the consistency of sequences across subjects. This is important because even in the absence of any real sequences of taskrelated representations spontaneous neural activity is not random but follows repeating dynamical motifs (Vidaurre et al., 2017). Solving this problem requires a randomized mapping between the assignment of physical stimuli to task states. This can be done across subjects, permitting valid inference at the group level.
At the group level, the statistical testing problem can be complicated by the fact that sequence measures do not in general follow a known distribution. Additionally, if a statetostate lag of interest ($\mathrm{\Delta}t$) is not known a priori, it is then necessary to perform tests at multiple lags, creating a multiple comparisons problem over a set of tests with complex interdependencies. In this section, we discuss inference with these issues in mind.
Distribution of sequenceness at a single lag
If a statetostate lag of interest ($\mathrm{\Delta}t$) is known a priori, then the simplest approach is to compare the sequenceness against zero, for example, using either a signedrank test or onesample t test (assuming Gaussian distribution). Such testing assumes the data are centred on zero if there were no real sequences. We show this approach is safe in both simulation (assuming no real sequences) and real MEG data where we know there are no sequences.
In simulation, we assume no real sequences, but state time courses are autocorrelated. At this point, there is no systematic structure in the correlation between the neuronal representations of different states (see later for this consideration). We then simply select the 40 ms time lag and compare its sequenceness to zero using either a signedrank test or onesample t test. We compare falsepositive rates predicted by the statistical tests with falsepositive rates measured in simulation (Figure 4a). We see the empirical false positives are well predicted by theory.
We have tested this also on real MEG data. In Liu et al., 2019, we had one condition where we measured resting activity before the subjects saw any stimuli. Therefore, by definition these sensory stimuli could not be replayed, we can use classifiers from these stimuli (measured later) to test a falsepositive performance of statistical tests on replay. Note, in our case, that each subject saw the same stimuli in a different order. They could not know the correct stimulus order when these resting data were acquired. These data provide a valid null for testing false positives.
To obtain many examples, we randomly permute the eight different stimuli 10,000 times and then compare sequenceness (at 40 ms time lag) to zero using either a signedrank test or onesample t test across subjects. Again, predicted and measured falsepositive rates match well (Figure 4b, left panel). This holds true across all computed time lags (Figure 4b, right panel).
An alternative to making assumptions about the form of the null distribution is to compute an empirical null distribution by permutation. Given that we are interested in the sequence of states over time, one could imagine permuting either state identity or time. However, permuting time uniformly will typically lead to a very high incidence of false positives as time is not exchangeable under the null hypothesis (Figure 4c, blue colour). Permuting time destroys the temporal smoothness of neural data, creating an artificially narrow null distribution (Liu et al., 2019; KurthNelson et al., 2016). This false positive also exists if we circular shift the time dimension of each state. This is because the signal is highly nonstationary. Replays come in bursts, as recently analysed (Higgins et al., 2021), and this will break a circular shift (Harris, 2020). State permutation, on the other hand, only assumes that state identities are exchangeable under the null hypothesis, while preserving the temporal dynamics of the neural data represents a safer statistical test that is well within 5% falsepositive rate (Figure 4c, purple colour).
Correcting for multiple comparisons
If the statetostate lag of interest is not known, we have to search over a range of time lags. As a result, we then have a multiple comparison problem. Unfortunately, we do not as yet have a good parametric method to control for multiple testing over a distribution. It is possible that one could use methods that exploit the properties of Gaussian random fields, as is common in fMRI (Worsley et al., 1996), but we have not evaluated this approach. Alternatively, we could use Bonferroni correction, but the assumption that each computed time lag is independent is likely false and overly conservative.
We recommend relying on state identitybased permutation. To control for the familywise error rate (assuming $\alpha =0.05$), we want to ensure there is a 5% probability of getting the tested sequenceness strength ($S}_{test$) or bigger by chance in *any* of the multiple tests. We therefore need to know what fraction of the permutations gives $S}_{test$ or bigger in any of their multiple tests. If any of the sequenceness scores in each permutation exceed $S}_{test$, then the maximum sequenceness score in the permutation will exceed $S}_{test$, so it is sufficient to test against the maximum sequenceness score in the permutation. The null distribution is therefore formed by first taking the peak of sequenceness across all computed time lags of each permutation. This is the same approach as used for familywise error correction for permutations tests in fMRI data (Nichols, 2012), and in our case it is shown to behave well statistically (Figure 4d).
What to permute
We can choose which permutations to include in the null distribution. For example, consider a task with two sequences, $Seq1:A\to B\to C\to D$ and $Seq2:E\to F\to G\to H$. We can form the null distribution either by permuting all states (e.g. one permutation might be E $\to F\to A\to B$, H $\to C\to E\to D$), as implemented in KurthNelson et al., 2016. Alternatively, we can form a null distribution which only includes transitions between states in different sequences (e.g. one permutation might be D $\to G\to A\to E$, H $\to C\to F\to B$), as implemented in Liu et al., 2019. In each case, permutations are equivalent to the test data under the assumption that states are exchangeable between positions and sequences. The first case has the advantage of many more possible permutations, and therefore may make more precise inferential statements in the tail. The second case may be more sensitive in the presence of a signal as the null distribution is guaranteed not to include permutations which share any transitions with the test data (Figure 4e). For example, in Figure 4e, the blue swaps are the permutations that only exchange state identity across sequences, as in Liu et al., 2019, while the red swaps are the permutations that permit all possible state identity permutations, as in KurthNelson et al., 2016. Note that there are many more different state permutations in red swaps than in blue swaps. We can make different levels of inferences by controlling the range of the null distributions in the state permutation tests.
Cautionary note on exchangeability of states after training
Until now, all nonparametric tests have assumed that state identity is exchangeable under the null hypothesis. Under this assumption, it is safe to perform state identitybased permutation tests on $Z}_{F$ and $Z}_{B$. In this section, we consider a situation where this assumption is broken.
More specifically, take a situation where the neural representation of states $A$ and $B$ is related in a systematic way or, in other words, the classifier on state $A$ is confused with state $B$, and we are testing sequenceness of $A\to B$. Crucially, to break the exchangeability assumption, representations of $A$ and $B$ have to be systematically more related than other states, for example, $A$ and $D$. This cannot be caused by lowlevel factors (e.g. visual similarity) because states are counterbalanced across subjects, so any such bias would cancel at the population level. However, such a bias might be induced by task training.
In this situation, it is, in principle, possible to detect sequenceness of $A\to B$ even in the absence of real sequences. In the autocorrelation section above, we introduced protections against the interaction of state correlation with autocorrelation. These protections may fail in the current case as we cannot use other states as controls (as we do in the multiple linear regression) because $A$ has systematic relationship with $B$, but not other states. State permutation will not protect us from this problem because state identity is no longer exchangeable.
Is this a substantive problem? After extensive training, behavioural pairing of stimuli can indeed result in increased neuronal similarity (Messinger et al., 2001; Sakai and Miyashita, 1991). These early papers involved long training in monkeys. More recent studies have shown induced representational overlap in human imaging within a single day (KurthNelson et al., 2015; Barron et al., 2013; Wimmer and Shohamy, 2012). However, when analysed across the whole brain, such representational changes tend to be localized to discrete brain regions (Schapiro et al., 2013; Garvert et al., 2017), and as a consequence may have limited impact on wholebrain decodeability.
Whilst we have not yet found a simulation regime in which false positives are found (as opposed to false negatives), there exists a danger in cases where, by experimental design, the states are not exchangeable.
Source localization
Uncovering temporal structure of neural representation is important, but it is also of interest to ask where in the brain a sequence is generated. Rodent electrophysiology research focuses mainly on the hippocampus when searching for replay. One advantage of wholebrain noninvasive neuroimaging over electrophysiology (despite many known disadvantages, including poor anatomical precision, low signalnoise ratio) is in its ability to examine neural activity in multiple other brain regions. Ideally, we would like a method that is capable of localizing sequences of more abstract representation in brain regions beyond hippocampus (Liu et al., 2019).
We want to identify the time when a given sequence is very likely to unfold, so we can construct averages of independent data over these times. We achieve this by transforming from the space of original states, $X}_{orig$, to the space of sequence events, $X}_{seq$. First, based on the transition of interest, $T$, we can obtain the projection matrix, $X}_{proj$:
If we know the state lag within sequence, $\mathrm{\Delta}t$ (e.g. the time lag give rise to the strongest sequenceness), or have it a priori, we can obtain the timelagged matrix, $X}_{lag$:
Then, we obtain state space with sequence event as states by elementwise multiply $X}_{proj$ and $X}_{lag$:
Each element in $X}_{seq$ indicates the strength of a (pairwise) sequence at a given moment in time. At this stage, $X}_{seq$ is a matrix with number of time points as rows (same as $X}_{orig$), and with number of pairwise sequences (e.g. A>B; B>C; etc.) as columns. Now on this matrix, $X}_{seq$, we can either look for sequences of sequences (see Appendix 3), or sum over columns (i.e. average over pairwise sequence events), and obtain a score at each time point reflecting how likely it is to be a sequence member (Figure 5a).
We can use this score to construct averages of other variables that might covary with replay. For example, if we choose time points when this score is high (e.g. 95th percentile) after being low for the previous 100 ms and construct an average timefrequency plot of the raw MEG data aligned to these times, we can reconstruct a timefrequency plot that is, on average, associated with replay onset (Figure 5b). Note that although this method assigns a score for individual replay events as an intermediary variable, it results in an average measure across many events.
This approach is similar to spiketriggered averaging (Sirota et al., 2008; Buzsáki et al., 1983). Applying this to real MEG data during rest, we can detect increased hippocampal power at 120–150 Hz, at replay onset (Figure 5b, c). Source reconstruction in the current analysis was performed using linearly constrained minimum variance (LCMV) beamforming, a common method for MEG source localization. This is known to suffer from distal correlated sources (Hincapié et al., 2017). A better method may be Empirical Bayesian Beamfomer for accommodating correlated neural source as a priori (O'Neill, 2021).
TDLM for rodent replay
So far, we have introduced TDLM in the context of analysing human MEG data. Relatedly, its application on human EEG data was also explored (Appendix 4: Apply TDLM to human wholebrain EEG data). Historically, replaylike phenomena have been predominantly studied in rodents with electrophysiology recordings in the hippocampal formation (Davidson et al., 2009; Grosmark and Buzsáki, 2016; Tingley and Peyrache, 2020). This raises interesting questions: how does TDLM compare to the existing rodent replay methods, can TDLM be applied to spiking data for detecting rodent replays, and what are the pros and cons? In this section, we address these questions.
Generality of graph vs. linebased replay methods
Given that TDLM works on the decoded state space, rather than sensor (with analogy to cell) level, we compared TDLM to rodent methods that work on the posterior decoded position (i.e., state) space, normally referred to as Bayesianbased methods (Tingley and Peyrache, 2020). (Note that these methods are typically Bayesian in how position is decoded from spikes [Zhang et al., 1998] but not in how replay is measured from decoded position.) Two commonly used methods are Radon transform (Davidson et al., 2009) and linear weighted correlation (Grosmark and Buzsáki, 2016).
Both methods proceed by forming a 2D matrix, where one dimension is the decoded state (e.g. positions on a linear track), and the other dimension is time (note that the decoded state is embedded in 1D). The methods then try to discover if an ordered line is a good description of the relationship between state and (parametric) time. For this reason, we call this family of approaches ‘line search’.
The radon method uses a discrete Radon transform to find the best line in the 2D matrix (Toft, 1996) and then evaluates the radon integral, which will be high if the data lie on a line (Figure 6a). It compares this to permutations of the same data where the states are reordered (Tingley and Peyrache, 2020). The linear weighted correlation method computes the average correlation between the time and estimated position in the 1D embedding (Figure 6b). The correlation is nonzero provided there is an orderly reactivation along the state dimension.
Both methods are applied to decoded positions, where they are sorted based on the order in a linearized state space. TDLM also works on the decoded position space, but instead of directly measuring the relationship between position and time, it measures the transition strength for each possible state to state transitions (Figure 6c).
This is a key difference between TDLM and these popular existing techniques. To reiterate, the latter rely on a continuous parametric embedding of behavioural states and time. TDLM is fundamentally different as it works on a graph and examines the statistical likelihood of some transitions happening more than others. This is therefore a more general approach that can be used for sequences drawn from any graph (e.g. 2D maze, Figure 6d), not just graphs with simple embeddings (like a linear track). For example, in a nonspatial decisionmaking task (KurthNelson et al., 2016), all states lead to two different states and themselves can be arrived at from two other different states (Figure 6e). Existing ‘line search’ methods will not work because there is no linear relationship between time and states (Figure 6f).
Multiscale TDLM
While continuous spaces can be analysed in TDLM by simply chunking the space into discrete states, TDLM in its original form may potentially be less sensitive for such analyses than techniques with builtin assumptions about the spatial layout of the state space, such as the linear relationship between time and reactivated states (Appendix 5 ‘Less sensitivity of TDLM to skipping sequences’). In essence, because TDLM works on a graph, it has no information about the Euclidean nature of the state space, while techniques that make assumptions about the linear relationship between space and time benefit from these assumptions. For example, detecting state 1 then state 5 then state 10 counts as replay in these techniques, but not in TDLM.
However, TDLM can be extended to address this problem. For continuous state spaces, we first need to decide how to best discretize the space. If we choose a large scale, we will miss replays that occur predominantly within a spatial bin. If we choose a small scale, we will miss transitions that jump spatial bins. A simple solution is to apply TDLM at multiple different scales and take an (varianceweighted) average of the sequenceness measures across different scales. For example, when measuring replay at the same speed, we can average events that travel 5 cm in 10 ms together with events that travel 10 cm in 20 ms.
Specifically, to perform multiscale TDLM, we discretize position bins at multiple widths. This generates rate maps at multiple scales (e.g. 5 cm, 10 cm, 20 cm, 40 cm), and hence a multiscale state space. For each replay speed of interest, we apply TDLM separately at each scale, and then take a varianceweighted average of replay estimates over all scales.
where $\beta}_{i$ is the sequence strength of given speed (i.e. statetostate lag) measured at scale $i$, $V}_{i$ is the variance of its $\beta}_{i$ estimator, and $n$ is the number of scales. In the end, statistical testing is performed on the precision weighted averaged sequence strength, $\beta}_{M$, in the same way as we do in the original TDLM.
It is easy to see why this addresses the potential concerns raised above as some scales will capture the 1 > 2 > 3 transitions, whilst others will capture the 1 > 10 > 20 transitions: because the underlying space is continuous, we can average results of the same replay speed together, and this will reinstate the Euclidean assumptions.
Applying multiscale TDLM to real rodent data (place cells in CA1)
We demonstrate the applicability of multiscale TDLM by analysing CA1 place cell spiking data from Ólafsdóttir et al., 2016. In Ólafsdóttir et al., 2016, rats ran multiple laps on a 600 cm Z maze and were then placed in a rest enclosure for 1.5 hr (Figure 7a). The Z maze consists of three tracks, with its ends and corners baited with sweetened rice to encourage running from one end to the other. The animal’s running trajectory was linearized, dwell time and spikes were binned into 2 cm bins and smoothed with a Gaussian kernel (σ = 5 bins). We generated rate maps separately for inbound (track 1 > track 2 > track 3) and outbound (track 3 > track 2 > track 1) running (see details in section ‘Rodent replay dataset’).
As in Ólafsdóttir et al., 2016, cells recorded in CA1 were classified as place cells if their peak firing field during track running was above 1 Hz with a width of at least 20 cm (see an example in Figure 7b). The candidate replay events were identified based on multiunit (MU) activity from place cells during rest time. Periods exceeding the mean rate by three standard deviations of MU activity were identified as possible replay events. Events less than 40 ms long, or which included activity from less than 15% of the recorded place cell ensemble, were rejected (see an example of putative replay event in Figure 7c), and the remaining events were labelled putative replay events.
We analysed data from one full recording session (track running for generating rate map, postrunning resting for replay detection) from Rat 2192 reported in Ólafsdóttir et al., 2016. Following the procedure described above, we identified 58 place cells and 1183 putative replay events. Replay analysis was then performed on the putative replay events, separately for inbound and outbound rate maps given the same position has a different decoded state depending on whether it was during an outbound or inbound run.
A forward sequence is characterized by states from the outbound map occurring in the outbound order or states from the inbound map occurring in the inbound order. Conversely, a backward sequence is when states from the inbound map occur in the outbound order or states from the outbound map occur in the inbound order. Candidate events were decoded based on a rate map, transforming the ncells * ntime to nstates * ntime. Each entry in this state space represents the posterior probability of being in this position at a given time. Replay analysis was performed solely on this decoded state space.
Note that TDLM is applied directly to the concatenated rather than individual replay events. This is because TDLM is a linear modelling framework. Applying TDLM on each single replay event and then averaging the beta estimates (appropriately weighted by the variances) is equivalent to running TDLM once on the concatenated replay events. It quantifies the average amount of replay across many events, which is different compared to existing replay methods that focus on single replay events. Because TDLM addresses statistical questions in linear modelling, it does not require secondary statistics to ask whether the ‘counts’ of individual events are more likely than chance or more likely in one situation than another.
During the whole sleep period, TDLM identified a significant forward sequence for the outbound map with a wide speed range around from 1 to 10 m/s (Figure 7d, left panel), consistent with recent findings from Denovellis, 2020 on varying replay speed (similar results were obtained for inbound map, not shown here for simplicity). In our analysis, the fastest speed is up to 10 m/s, which is around 20× faster than its free running speed, representing approximately half a trackarm in a typical replay event, consistent with previous work (Lee and Wilson, 2002; Davidson et al., 2009; Karlsson and Frank, 2009; Nádasdy et al., 1999).
Secondorder inferences
As pointed out by van der Meer et al., 2020, there are two types of statistical questions: a ‘firstorder’ sequence question, which concerns whether an observed sequenceness is different from random (i.e. do replays exist?); and a ‘secondorder’ question, which requires a comparison of sequenceness across conditions (i.e. do replays differ?). Because it is embedded in a linear regression framework, TDLM is ideally placed to address such questions. There are two ways of asking such questions in linear modelling: contrasts and interactions. We explain them with examples here.
Linear contrasts
After fitting a regression model, resulting in coefficients for different regressors, we can test hypotheses about these coefficients by constructing linear combinations of the coefficients that would be zero under the null hypothesis. For example, if we want to test whether effect A is greater than effect B, then we can compute the linear contrast A – B (which would be zero under the null hypothesis) and perform statistics on this new measure. If we want to test whether replay increases linearly over five conditions [A, B, C, D, E], we can compute the linear contrast −2*A – B + 0*C + D + 2*E (which would be zero under the null hypothesis) and perform statistics on this new measure. Statistics (within or across animals) can operate with these contrasts in exactly the same way as with the original coefficients from the linear model. Here, we demonstrate this by showing in our example dataset that there was a greater preponderance for forward than backward replay. We construct the contrast (forwards – backwards) and test it against zero using a multiplecomparisoncontrolled permutation test (Figure 7d, right panel, pink line). By constructing a different contrast (forwards + backwards), we can also show that the total replay strength across both types of replays was significant (Figure 7d, right panel, green line).
Interactions
A second method for performing secondorder tests is to introduce them into the linear regression as interaction terms, and then perform inference on the regression weights for these interactions. This means changing Equation 2 to include new regressors. For example, if interested in how reactivations change over time, one could build new regressors ($Xtim{e}_{k}\left(t\right)$), obtained by elementwise multiplying the state regressor, e.g. ${X}_{k}\left(t\right)$ with time indices ($Xtim{e}_{k}\left(t\right)={X}_{k}\left(t\right).\ast \mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}$). Now the firstlevel GLM is constructed as (omitting residual term ε, same as Equation 2):
Example regressors in the design matrix can be seen in Figure 7e. The first regressor, ${X}_{k}\left(t\right)$, is one of the state reactivation regressors used in standard TDLM. The second regressor, $Xtim{e}_{k}\left(t\right)$, is the same as ${X}_{k}\left(t\right)$ multiplied by time. (There are k regressors of each form in regressor matrix.) Here, we chose to demean the time regressor before the interaction, so the early half of the regressor is negative and the late half is positive. This has no effect on the regression coefficients of the interaction term, but, by rendering the interaction approximately orthogonal to ${X}_{k}\left(t\right)$, it makes it possible to estimate the main effect and the interaction in the same regression.
Note that the interaction regressor is orthogonal to the state reactivation regressor, so it will have no effect on the firstorder regression terms. If we include such regressors for all states, then we can get two measures for each replay direction (sequence effect and time effect). The first tells us the average amount of replay throughout the sleep period (first order). The second tells us whether replay increases or decreases as time progresses through the sleep period (second order).
Orthogonal tests in regions of interest
When examining forward–backward replay above, we did separate inference for each replay speed, and then performed multiple comparison testing using the maxpermutation method (see section 'Statistical inference'). We now take the opportunity to introduce another method common in human literature.
To avoid such multiple comparison correction, it is possible to select a ‘region of interest’ (ROI), average the measure in question over that ROI, and perform inference on this average measure. Because we are now only testing one measure, there is no multiple comparison problem. Critical in this endeavour, however, is that we do not use the measure under test or anything that correlates with that measure as a means to define the ROI. This will induce a selection bias (Kriegeskorte et al., 2009). In the example in Figure 7f, we have used the average replay (forwards + backwards) to select the ROI. We are interested in speeds in which there is detectable replay on average across both directions and the whole sleep period (Figure 7d, right panel, green shaded area). If we select our ROI in this way, we cannot perform unbiased inference on firstorder forward or backward replay because forward and backward regressors correlate with their sum (Figure 7f, statistical inference in the red rectangle is biased). However, we can perform unbiased inference on several secondorder effects (Figure 7f, statistical inference in the green rectangle). We can test (forwards – backwards) assuming the difference of terms is orthogonal to their sum (as it is in this case). Further, we can test any interaction with time because the ROI is defined on the average over time and the interaction looks for differences as a function of time. When we perform these tests in our example dataset (Figure 7f, green rectangle), we confirm that there are more forward than backward replay on average. We further show that forward replay is decreasing with time during sleep, and that backward replay is increasing with time. Their difference (forwards – backwards) is also significant.
In addition to the timevarying effect, we can also test the spatial modulation effect, that is, how replay strength (at the same replay speed) changes as a function of its spatial content. For example, is replay stronger for transitions in the start of track compared to the end of the track? As an illustrative example, we have used the same ROI defined above and test the spatial modulation effect on forward replay. Note that this test of spatial modulation effect is also unbiased from the overall strength of forward replay, and thereby no selection bias in this ROI, as well.
For visualization purposes, we have first plotted the estimated strength for each pairwise forward sequence (Figure 8a), separately within each scale (from 1 to 4, with increasing spatial scales). The pairwise sequences are ordered from the start of the track to the end of the track. Alongside the pairwise sequence plot, we have plotted the mean replay strength over all possible pairwise transitions (in red) in comparison to the mean of all control transitions (in grey; as expected, they are all around 0). Note that we cannot perform inference on the difference between the red and grey bars here because they have been selected from a biased ROI. It is simply for illustration purposes. We have therefore put them in red squares to match Figure 7f.
To formally test the spatial modulation effect, we can use the exact same approach as outlined above in section 'Linear contrasts'. Here, we test a linear increase or decrease across different transitions. We take the linear contrast weight vector, $c$ ([2,1,0,1,2] for the largest scale, [3:3] for the next scale, [5:5] for the next scale, and [12:12] for the smallest scale), and multiply these by the beta estimates of the transitions:
If this new measure, $contrast$, is different from zero, then there is a linear increase/decrease from one end of the track to the other. Note that this new contrast is no longer biased by the ROI selection as each transition contributed equally to the ROI selection, but we are now comparing between transitions. Inference on this contrast is therefore valid. We have therefore put them in green boxes to match Figure 7f (Figure 8b, c).
Within the larger two scales, these contrasts are significantly negative (tested against permutations in exactly the same way as the ‘mean’ contrasts). Since we are still in the linear domain, we can now just average these contrasts across the four scales and get a single measure for spatial modulation of replay. This average measure is significantly negative (Figure 8b). Hence, on average, forward replay is stronger at the beginning of the track.
We can do the same thing for backward replay. We found an opposite pattern, that is, strength of backward replay is stronger at the end of the track, and similarly, it is not significant in the smallest scale and becomes significant in the largest scale, and also significant on average across all scales (Figure 8c). Again, since we are in the linear domain, we can further contrast these contrasts, asking if this effect is different for forward and backward replay. We found that the difference is indeed significant (Figure 8d). This set of results is consistent with previous rodent literature (Diba and Buzsáki, 2007). Note that we would like to stress again that this analysis is not about a single replay event but is testing for average differences across all replay events.
Notably, extra care needs to be exercised for secondorder questions (compared to firstorder ones). Problems can emerge due to biases in secondorder inference, such as in behavioural sampling (e.g. track 1 may be experienced more than track 2 during navigation; this creates a bias when evaluating replay in tack 1 vs. track 2 during rest). Such issues are real but can be finessed by experimental design considerations of a sort commonly applied in the human literature. For example:
Ensure that biases that might occur within subjects will not occur consistently in the same direction across subjects (e.g. by randomizing stimuli across participants).
Compare across conditions in each subject.
Perform a random effects inference across the population by comparing against the betweensubject variance.
Such approaches are not yet common in rodent electrophysiology and may not be practical in some instances. In such cases, it remains important to be vigilant to guard against these biases with TDLM as with other techniques. If these approaches are feasible, the machinery for computing secondorder inferences is straightforward in a linear framework like TDLM.
Generality of TDLM
We have now discussed the applicability of TDLM in relation to human MEG, as well as in rodent electrophysiology (with comparisons to standard replay detection methods). A preliminary attempt at detecting replay in human EEG is also shown in Appendix 4. We believe that this establishes TDLM as a domaingeneral sequence analysis method: TDLM works at the level of decoded state space, rather than the sensor/cell level of the data. It can be applied to a wide range of data types and settings in both humans and rodents, stimulating crossfertilization across disciplines. It is based on the GLM framework, and this lends it flexibility for regressing out potential confounds while offering an intuitive understanding of the overall approach.
In this section, we discuss the generality of TDLM.
States
TDLM assesses the statistical likelihood of certain transitions on a graph. In its original form, TDLM works on discrete states (i.e. nodes in the graph). Continuous spaces can be incorporated by chunking them into discrete spaces. Furthermore, by averaging the same replay speeds measured at multiple scales of discretization (see section ‘TDLM for rodent replay’), the statistical benefits of an assumption of a Euclidean geometry can be recovered.
Time length
The longer the time length, the more accurate the estimates in TDLM. This is because TDLM assesses sequence evidence based on a GLM framework, where time length is the sample size. Higher sample size will lead to more accurate estimates. In the case of rodent analysis, we recommend applying TDLM to aggregated replay events rather than to a single event because this results in (1) more time samples for estimation and (2) more activated states in the analysis time framework. Unlike other techniques which search for a single replay in a single event, this aggregation can be implemented without losing generality as TDLM is able to handle multiple sequences in the same data with respect to different directions, contents, or speeds. Furthermore, by aggregating linearly across all replay events of the same condition, it provides a natural measure for comparing replay strength, speed, and direction across different experimental conditions.
TDLM has already proved important in human experiments where complex state spaces have been used (Wimmer et al., 2020; Liu et al., 2019; Liu et al., 2021a; KurthNelson et al., 2016). We expect this generality will also be important as rodent replay experiments move beyond 1D tracks, for example, to foraging in 2D, or in complex mazes.
Discussion
TDLM is a domaingeneral analysis framework for capturing sequence regularity of neural representations. It is developed on human neuroimaging data and can be extended to other data sources, including rodent electrophysiology recordings. It offers hope for crossspecies investigations on replay (or neural sequences in general) and potentially enable studies of complex tasks in both humans and animals.
TDLM adds a new analysis toolkit to the replay field. It is especially suited for summarizing replay strength across many events, comparing replay strength between conditions, and analysing replay strength in complex behavioural paradigms. Its linear modelling nature makes it amenable to standard statistical tests and thereby allows wide use across tasks, modalities, and species. Unlike alternative tools, we have not shown TDLM applied to individual replay events.
The temporal dynamics of neural states have been studied previously with MEG (Vidaurre et al., 2017; Baker et al., 2014). Normally such states are defined by common physiological features (e.g. frequency, functional connectivity) during rest and termed resting state networks (e.g. default mode network [Raichle et al., 2001]). However, these approaches remain agnostic about the content of neural activity. The ability to study the temporal dynamics of representational content permits richer investigations into cognitive processes (Higgins et al., 2021) as neural states can be analysed in the context of their roles with respect to a range of cognitive tasks.
Reactivation of neural representations has also been studied previously (Tambini and Davachi, 2019) using approaches similar to the decoding step of TDLM or multivariate pattern analysis (MVPA) (Norman et al., 2006). This has proven fruitful in revealing mnemonic functions (Wimmer and Shohamy, 2012), understanding sleep (Lewis and Durrant, 2011), and decisionmaking (Schuck et al., 2016). However, classification alone does not reveal the rich temporal structures of reactivation dynamics. We have described the application of TDLM mostly during offtask state in this paper. The very same analysis can be applied to ontask data to test for cued sequential reactivation (Wimmer et al., 2020) or sequential decisionmaking (Eldar et al., 2020). For example, the ability to detect sequences ontask allows us to tease apart clustered from sequential reactivation, where this may be important for dissociating decision strategies (Eldar et al., 2018) and their individual differences (Wimmer et al., 2020; Eldar et al., 2020). TDLM, therefore, may allow testing of neural predictions from process models such as reinforcement learning during task performance (Dayan and Daw, 2008), which have proved hard to probe previously (Wimmer et al., 2020; Nour et al., 2021; Liu et al., 2019; Liu et al., 2021a).
In the human neuroimaging domain, we have mainly discussed the application of TDLM with respect to MEG data. In Appendix 4, we show that TDLM also works well with EEG data. This is not surprising given EEG and MEG are effectively measuring the same neural signature, namely local field potential (or associated magnetic field) on the scalp. We do not have suitable fMRI data to test TDLM. However, related work has suggested that it might be possible to measure sequential reactivation using fMRI (Schuck and Niv, 2019), but particular methodological caveats need to be considered (e.g. a bias from last events due to slow haemodynamic response) (Wittkuhn and Schuck, 2021). We believe that TDLM can deal with this, given it models out nonspecific transitions, although further work is needed. In future, we consider it will be useful to combine the high temporal resolution available in MEG/EEG and the spatial precision available in fMRI to probe regionspecific sequential computations.
In the rodent electrophysiology domain, we show what TDLM (its multiscale version) has to offer uniquely compared to existing rodent replay methods. Most importantly, TDLM works on an arbitrary graph and its generality makes replay studies in complex mazes possible. Its linear framework makes the assessment of timevarying effect on replay (Figure 7) or other secondorder sequence questions straightforward. In future work, a promising direction will be to further separate process noise (e.g. intrinsic variability within sequences) and measurement noise (e.g. noise in MEG recording). This might be achieved by building latent statespace models as have been explored recently in rodent community (Maboudi et al., 2018; Denovellis, 2020).
Together, we believe that TDLM opens doors for novel investigations of human cognition, including language, sequential planning, and inference in nonspatial cognitive tasks (Eldar et al., 2018; KurthNelson et al., 2016), as well as complicated tasks in rodents, for example, forging in 2D mazes. TDLM is particularly suited to test specific neural predictions from process models, such as reinforcement learning. We hope that TDLM can promote an acrossspecies synthesis between experimental and theoretical neuroscience and, in so doing, shed novel light on neural computation.
Materials and methods
Simulating MEG data
We simulate the data so as to be akin to human MEG.
Task data for obtaining state patterns
Request a detailed protocolWe generate ground truth multivariate patterns (over sensors) of states. We then add random Gaussian noise on the ground truth state patterns to form the task data. We train a logistic regression classifier on the task data so as to obtain a decoding model for each of the state patterns. Later we use this decoding model to transform the restingstate data from sensor space (with dimension of time by sensors) to the state space (with dimension of time by states).
Rest data for detecting sequences
Request a detailed protocolFirst, to imitate temporal autocorrelations and spatial correlations commonly seen in human neuroimaging data, we generate the rest data using an autoaggressive model with multivariate (over sensors) Gaussian noise and add a dependence among sensors. In some simulations, we also add a rhythmic oscillation (e.g. 10 Hz).
Second, we inject a sequence of state patterns in the rest data. The sequences follow the ground truth of state transitions of interest. The statetostate time lag is assumed to follow a gamma distribution. We vary the number of sequences to be injected in the rest data to control the strength of sequences.
Lastly, we project the rest data to the decoding model of states obtained from the task data. TDLM will then work on the decoded state space.
An example of the MATLAB implementation is called ‘Simulate_Replay’ from the Github link: https://github.com/yunzheliu/TDLM (copy archived at swh:1:rev:015c0e90a14d3786e071345760b97141700d6c85), Liu, 2021b.
Human replay dataset
Task design
Request a detailed protocolParticipants were required to perform a series of tasks with concurrent MEG scanning (see details in Liu et al., 2019). The functional localizer task was performed before the main task and was used to train a sensory code for eight distinct objects. Note that the participants were provided with no structural information at the time of the localizer. These decoding models, trained on the functional localizer task, capture a sensorylevel neural representation of stimuli (i.e. stimulus code). Following that, participants were presented with the stimuli and were required to unscramble the ‘visual sequence’ into a correct order, that is, the ‘unscrambled sequence’ based on a structural template they had learned the day before. After that, participants were given a rest for 5 min. At the end, stimuli were presented again in random order, and participants were asked to identify the true sequence identity and structural position of the stimuli. Data in this session are used to train a structural code (position and sequence) for the objects.
MEG data acquisition, preprocessing, and source reconstruction
Request a detailed protocolWe follow the same procedure that has been reported in Liu et al., 2019. We have copied it here for references.
'MEG was recorded continuously at 600 samples/s using a wholehead 275channel axial gradiometer system (CTF Omega, VSM MedTech), while participants sat upright inside the scanner. Participants made responses on a button box using four fingers as they found most comfortable. The data were resampled from 600 to 100 Hz to conserve processing time and improve signaltonoise ratio. All data were then highpassfiltered at 0.5 Hz using a firstorder IIR filter to remove slow drift. After that, the raw MEG data were visually inspected, and excessively noisy segments and sensors were removed before independent component analysis (ICA). An ICA (FastICA, http://research.ics.aalto.fi/ica/fastica) was used to decompose the sensor data for each session into 150 temporally independent components and associated sensor topographies. Artefact components were classified by combined inspection of the spatial topography, time course, kurtosis of the time course, and frequency spectrum for all components. Eyeblink artefacts exhibited high kurtosis (>20), a repeated pattern in the time course and consistent spatial topographies. Mains interference had extremely low kurtosis and a frequency spectrum dominated by 50 Hz line noise. Artefacts were then rejected by subtracting them out of the data. All subsequent analyses were performed directly on the filtered, cleaned MEG signal, in units of femtotesla.
All source reconstruction was performed in SPM12 and FieldTrip. Forward models were generated on the basis of a single shell using superposition of basis functions that approximately corresponded to the plane tangential to the MEG sensor array. LCMV beamforming (Van Veen et al., 1997) was used to reconstruct the epoched MEG data to a grid in MNI space, sampled with a grid step of 5 mm. The sensor covariance matrix for beamforming was estimated using data in either broadband power across all frequencies or restricted to ripple frequency (120–150 Hz). The baseline activity was the mean neural activity averaged over −100 ms to −50 ms relative to sequence onset. All nonartefactual trials were baseline corrected at source level. We looked at the main effect of the initialization of sequence. Nonparametric permutation tests were performed on the volume of interest to compute the multiple comparison (wholebrain corrected) pvalues of clusters above 10 voxels, with the null distribution for this cluster size being computed using permutations (n = 5000 permutations)'.
Rodent replay dataset
Data description
Request a detailed protocolThis data is from Ólafsdóttir et al., 2016. We analysed one full recording session (track running for generating rate map, postrunning resting for replay detection) from Rat 2192.
Task description
Request a detailed protocolIn Ólafsdóttir et al., 2016, rats ran multiple laps on a Z maze and were then placed in a rest enclosure. The two parallel sections of the Z (190 cm each) were connected by a diagonal section (220 cm). Animals were pretrained to run on the track. At the recording session, rats were placed at one end of the Ztrack. The ends and corners of the track were baited with sweetened rice to encourage running from one end to the other. In each session, rats completed 20 full laps (30–45 min). Following the track session, rats were placed in the rest enclosure for 1.5 hr.
Preprocessing
Request a detailed protocolFollowing Ólafsdóttir et al., 2016, when generating rate maps we excluded data from both the ends and corners because the animals regularly performed nonperambulatory behaviours there. Periods when running speed was less than 3 cm/s were also excluded. Running trajectories were then linearized, dwell time and spikes were binned into 2 cm bins and smoothed with a Gaussian kernel (σ = 5 bins). We generated rate maps separately for inbound (track 1 > track 2 > track 3) and outbound (track 3 > track 2 > track 1) running.
As in Ólafsdóttir et al., 2016, cells recorded in CA1 were classified as place cells if their peak firing field during track running was above 1 Hz and at least 20 cm wide. The candidate replay events were identified based on MU activity from place cells during rest time. Only periods exceeding the mean rate by three standard deviations of MU activity were identified as putative replay events. Events less than 40 ms long or which included activity from less than 15% of the recorded place cell ensemble were rejected.
We analysed data from one full recording session (track running for generating rate map, postrunning resting for replay detection) of Rat 2192 reported in Ólafsdóttir et al., 2016. Following the procedure described above, we have identified 58 place cells and 1183 putative replay events. Replay analysis was then performed on the putative replay events, separately for inbound and outbound rate maps.
Code availability
Request a detailed protocolSource code of TDLM can be found at https://github.com/yunzheliu/TDLM.
Data availability
Request a detailed protocolData are also available at https://github.com/yunzheliu/TDLM.
Appendix 1
Multistep sequences
TDLM can be used iteratively. One extension of TDLM of particular interest is: multistep sequences. It asks about a consistent regularity among multiple states.
So far, we introduced methods for quantifying the extent to which the statetostate transition structure in neural data matches a hypothesized taskrelated transition matrix. An important limitation of these methods is that they are blind to hysteresis in transitions. In other words, they cannot tell us about multistep sequences. In this section, we describe a methodological extension to measure evidence for sequences comprising more than one transition: for example, $A\to B\to C$.
The key ingredient is controlling for shorter subsequences (e.g. $A\to B$ and $B\to C$) in order to find evidence unique to a multistep sequence of interest.
Assuming constant statetostate time lag, $\mathrm{\Delta}t$, between A and B, and between B and C. We can create a new state space AB by shifting B up $A\to B$, and elementwise multiply it with state A. This new state AB measures the reactivation strength of $A\to B$, with time lag $\mathrm{\Delta}t$. In the same way, we can create a new state space, BC, AC, etc. Then we can construct the same firstlevel GLM on the new state space. For example, if we want to determine the evidence of $A\to B\to C$ at time lag $\mathrm{\Delta}t$, we can regress AB onto state time course C at each $\mathrm{\Delta}t$ (Equation 1). But we want to know the unique contribution of AB to C. More specifically, we want to test if the evidence of $A\to B\to C$ is stronger than $X\to B\to C$, where X is any other state but not A. Therefore, similar to Equation 2, we need to control CB, DB, when looking for evidence of AB of C. Applying this method, we show that TDLM successfully avoids false positives arising out of strong evidence for shorter length (see simulation results in Appendix 1—figure 1a, and see results obtained on human neuroimaging data in Appendix 1—figure 1b). This process can be generalized to any number of steps.
TDLM, in its current form, assumes a constant intrasequence statetostate time lag. If there is a variability between state transitions TDLM can still cope, but not very elegantly. Assume there is a threestate sequence, $A\to B\to C$, with intrasequence variance. TDLM will need to test all possible combinations of statetostate time lags in $A\to B$ and $B\to C$. If there are $n$ number of time lag of interest in either of the two transitions, TDLM will then have to test $n}^{2$ possible time lag combinations. This is a large search space and one that increases exponentially as a function of the length of a sequence.
We note that this analysis is different from a typical rodent replay analysis which assesses the overall evidence for a sequence length (Davidson et al., 2009; Grosmark and Buzsáki, 2016). TDLM asks if there is more evidence for A > B > C, above and beyond evidence for B > C, for example. If the main question of interest is ‘do we have evidence of A > B > C in general’, as normally is the case in the rodent replay analysis (Davidson et al., 2009; Grosmark and Buzsáki, 2016), we should not control for shorter lengths. Instead, we can simply average the evidence together, as implemented in KurthNelson et al., 2016.
Appendix 2
Pseudocode of sensory code and abstract code crossvalidations
Appendix 3
Sequences of sequences
We have detailed the use of either sensory or abstract representations as the states in TDLM. We now take a step further and use sequences themselves as states. Using this kind of hierarchical analysis, we can search for sequences of sequences. This is useful because it can reveal temporal structure not only within sequence, but also between sequences. The organization between sequences is of particular interest for revealing neural computations. For example, the forward and backward search algorithms hypothesized in planning and inference (Penny et al., 2013) can be cast as sequences of sequences problem: the temporal structure of forward and backward sequence. This can be tested by using TDLM iteratively.
To look for sequences between sequences, we need first to define sequences as new states. To do so, the raw state course, for example, state B, needs to be shifted up by the empirical withinsequence time lag $\mathrm{\Delta}t$ (determined by the twolevel GLM), to align with the onset of state A, if assuming sequence $A\to B$ exist (at time lag $\mathrm{\Delta}t$). Then, we can elementwise multiply the raw state time course A with the shifted time course B, resulting in a new state AB. Each entry in this new state time course indicates the reactivation strength of sequence AB at a given time.
The general twolevel GLMs framework still applies, but now with one important caveat. The new sequence state (e.g. AB) is defined based on the original states (A and B), and where we are now interested in a reactivation regularity, that is, sequence, between sequences, rather than the original states. We need therefore to control for the effects of the original states. Effectively, this is like controlling for main effects (e.g. state A and shifted state B) when looking for their interaction (sequence AB). TDLM achieves this by including timelagged original state regressors A, B, in addition to AB, in the firstlevel GLM sequence analysis.
Specifically, let us assume that the sequence state matrix is $X}_{seq$, after transforming the original state space to sequence space based on the empirical withinsequence time lag $\mathrm{\Delta}{t}_{w}$. Each column at $X}_{seq$ is sequence state, denoted by $S}_{ij$, which indicates the strength of sequence i > j reactivation. The raw state i is $X}_{i$, and the shifted raw state j is $X}_{jw$ (by time lag $\mathrm{\Delta}{t}_{w}$).
In the first level GLM, TDLM ask for the strength of a unique contribution of sequence state $S}_{ij$ to $S}_{mn$ while controlling for original states ($X}_{i$ and $X}_{jw$). For each sequence state $ij$, at each possible time lag $\mathrm{\Delta}t$, TDLM estimated a separate linear model:
Repeat this process for each sequence state separately at each time lag, resulting in a sequence matrix $\beta}_{seq$.
At the secondlevel GLM, TDLM asks how strong the evidence for a sequence of interest is compared to sequences that have the same starting state or end state at each time lag. This secondlevel GLM will be the same as Equation 5, but with additional regressors to control for sequences that share the same start or end state. In simulation, we demonstrate, applying this method, that TDLM can uncover hierarchical temporal structure: state A is temporally leading state B with 40 ms lag, and the sequence A > B tends to repeat itself with a 140 ms gap (Appendix 3—figure 1a). One interesting application of this is to look for theta sequence (Mehta et al., 2002; McNaughton et al., 2006; Buzsáki and Moser, 2013). One can think of theta sequence, a welldocumented phenomenon during rodent spatial navigation, as a neural sequence repeating itself in theta frequency (6–12 Hz).
In addition to looking for temporal structure of the same sequence, the method is equally suitable when searching for temporal relationships between different sequences. For example, assuming two different types of sequences, one sequence type has a withinsequence time lag at 40 ms; while the other has a withinsequence time lag at 150 ms (Appendix 3—figure 1b, left and middle panel); and there is a gap of 200 ms between the two types of sequences (Appendix 3—figure 1b, right panel). These time lags are set arbitrarily for illustration purposes. TDLM can accurately capture the dynamics both within and between the sequences, supporting a potential for uncovering temporal relationships between sequences under the same framework.
Appendix 4
Apply TDLM to human wholebrain EEG data
An autocorrelation is commonplace in neuroimaging data, including EEG and fMRI. TDLM is designed to specifically take care of this confound, and, on this basis, we should be able to work with EEG and fMRI data. We do not have suitable fMRI data available to test TDLM but are interested to investigate this in more depth in our future work. We had collected EEG data from one participant to test whether TDLM would *just* work.
The task was designed to examine ontask sequential replay in decisionmaking by Dr. Toby Wise. This is a ‘Tmaze’ like task, where a participant needs to choose a left or right path based on the value received at the end of the path. We could decode seven objects well on the wholebrain EEG data using just raw amplitude feature (same with our MEGbased analysis) and could detect fast backward sequenceness (peaked at 30 ms time lag) during choice/planning time (Appendix 4—figure 1), similar to our previous MEG findings (KurthNelson et al., 2016). As this result is from one subject, we are cautious about making an excessive claim, but nevertheless we believe that the data show the TDLM approach is highly promising for EEG data.
Appendix 5
Less sensitivity of TDLM to skipping sequences
In a linear track where replays only go in a single direction, it is possible that TDLM is less sensitive compared to the linear correlation or the Radon method, given the latter assumes a parametric relationship between space and time. For example, if only the first and last states are activated, but not the intermediate states, the existing methods will report replay, but TDLM will not, because in existing methods space and time are parametric quantities (Appendix 5—figure 1). In contrast, TDLM only knows about transitions on a graph.
Data availability
No new data is used or generated in the current paper. Data relevant for current paper is available at https://github.com/YunzheLiu/TDLM (copy archived at https://archive.softwareheritage.org/swh:1:rev:015c0e90a14d3786e071345760b97141700d6c85). This dataset is from Ólafsdóttir et al., 2016.
References

Online evaluation of novel choices by simultaneous representation of multiple memoriesNature Neuroscience 16:1492–1498.https://doi.org/10.1038/nn.3515

Repetition suppression: a means to index neural representations using BOLD?Philosophical Transactions of the Royal Society B: Biological Sciences 371:20150355.https://doi.org/10.1098/rstb.2015.0355

Crossspecies neuroscience: closing the explanatory gapPhilosophical Transactions of the Royal Society B: Biological Sciences 376:20190633.https://doi.org/10.1098/rstb.2019.0633

Cellular bases of hippocampal EEG in the behaving ratBrain Research Reviews 6:139–171.https://doi.org/10.1016/01650173(83)900371

Memory, navigation and theta rhythm in the hippocampalentorhinal systemNature Neuroscience 16:130–138.https://doi.org/10.1038/nn.3304

Decision theory, reinforcement learning, and the brainCognitive, Affective, & Behavioral Neuroscience 8:429–453.https://doi.org/10.3758/CABN.8.4.429

Auto‐regressive model for nonstationary stochastic processesJournal of Engineering Mechanics 114:1995–2012.https://doi.org/10.1061/(ASCE)07339399(1988)114:11(1995)

Forward and reverse hippocampal placecell sequences during ripplesNature Neuroscience 10:1241–1242.https://doi.org/10.1038/nn1961

Granger causality and path diagrams for multivariate time seriesJournal of Econometrics 137:334–353.https://doi.org/10.1016/j.jeconom.2005.06.032

Replay comes of ageAnnual Review of Neuroscience 40:581–602.https://doi.org/10.1146/annurevneuro072116031538

Decoding neural representational spaces using multivariate pattern analysisAnnual Review of Neuroscience 37:435–456.https://doi.org/10.1146/annurevneuro062012170325

ThesisUncovering temporal structure in neural data with statistical machine learning modelsUniversity of Oxford.

Awake replay of remote experiences in the HippocampusNature Neuroscience 12:913–918.https://doi.org/10.1038/nn.2344

Representational similarity analysis  connecting the branches of systems neuroscienceFrontiers in Systems Neuroscience 2:4.https://doi.org/10.3389/neuro.06.004.2008

Circular analysis in systems neuroscience: the dangers of double dippingNature Neuroscience 12:535–540.https://doi.org/10.1038/nn.2303

Overlapping memory replay during sleep builds cognitive schemataTrends in Cognitive Sciences 15:343–351.https://doi.org/10.1016/j.tics.2011.06.004

Viewpoints: how the Hippocampus contributes to memory, navigation and cognitionNature Neuroscience 20:1434–1447.https://doi.org/10.1038/nn.4661

Path integration and the neural basis of the 'cognitive map'Nature Reviews Neuroscience 7:663–678.https://doi.org/10.1038/nrn1932

Replay and time compression of recurring spike sequences in the HippocampusThe Journal of Neuroscience 19:9497–9507.https://doi.org/10.1523/JNEUROSCI.192109497.1999

Beyond mindreading: multivoxel pattern analysis of fMRI dataTrends in Cognitive Sciences 10:424–430.https://doi.org/10.1016/j.tics.2006.07.005

Coordinated grid and place cell replay during restNature Neuroscience 19:792–794.https://doi.org/10.1038/nn.4291

The role of hippocampal replay in memory and PlanningCurrent Biology 28:R37–R50.https://doi.org/10.1016/j.cub.2017.10.073

Forward and backward inference in spatial cognitionPLOS Computational Biology 9:e1003383.https://doi.org/10.1371/journal.pcbi.1003383

Neural representations of events arise from temporal community structureNature Neuroscience 16:486–492.https://doi.org/10.1038/nn.3331

Awake reactivation of prior experiences consolidates memories and biases cognitionTrends in Cognitive Sciences 23:876–890.https://doi.org/10.1016/j.tics.2019.07.008

On the methods for reactivation and replay analysisPhilosophical Transactions of the Royal Society B: Biological Sciences 375:20190231.https://doi.org/10.1098/rstb.2019.0231

ThesisThe Radon Transform: Theory and ImplementationTechnical University of Denmark.

Progress and issues in secondorder analysis of hippocampal replayPhilosophical Transactions of the Royal Society B: Biological Sciences 375:20190238.https://doi.org/10.1098/rstb.2019.0238

Localization of brain electrical activity via linearly constrained minimum variance spatial filteringIEEE Transactions on Biomedical Engineering 44:867–880.https://doi.org/10.1109/10.623056

ConferenceAdvances in neural information processing systemsConference on Neural Information Processing Systems. pp. 1473–1480.

Episodic memory retrieval success is associated with rapid replay of episode contentNature Neuroscience 23:1025–1033.https://doi.org/10.1038/s415930200649z
Article and author information
Author details
Funding
Wellcome (098362/Z/12/Z)
 Raymond J Dolan
Wellcome (104765/Z/14/Z)
 Timothy E Behrens
Wellcome (219525/Z/19/Z)
 Timothy E Behrens
James S. McDonnell Foundation (JSMF220020372)
 Timothy E Behrens
Wellcome (212281/Z/18/Z)
 Caswell Barry
Max Planck Society
 Yunzhe Liu
Wellcome (203139/Z/16/Z)
 Timothy E Behrens
Wellcome (091593/Z/10/Z)
 Raymond J Dolan
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Matthew A Wilson for help with rodent theta sequence analysis. We thank Elliott Wimmer and Toby Wise for helpful discussion and generous sharing of their data. We thank Matt Nour for helpful comments on a previous version of the manuscript. YL is also grateful for the unique opportunity provided by the Brains, Minds and Machines Summer Course. We acknowledge fundings from the Open Research Fund of the State Key Laboratory of Cognitive Neuroscience and Learning to YL, Wellcome Trust Investigator Award (098362/Z/12/Z) to RJD, Wellcome Trust Senior Research Fellowship (104765/Z/14/Z), and Principal Research Fellowship (219525/Z/19/Z), together with a James S McDonnell Foundation Award (JSMF220020372), to TEJB; and Wellcome Trust Senior Research Fellowship (212281/Z/18/Z) to CB. Both Wellcome Centres are supported by core funding from the Wellcome Trust: Wellcome Centre for Integrative Neuroimaging (203139/Z/16/Z), Wellcome Centre for Human Neuroimaging (091593/Z/10/Z). The Max Planck UCL Centre is a joint initiative supported by UCL and the Max Planck Society.
Ethics
Human subjects: The human dataset used in this study was reported in Liu et al 2019. All participants were recruited from the UCL Institute of Cognitive Neuroscience subject pool, had a normal or correctedtonormal vision, no history of psychiatric or neurological disorders, and had provided written informed consent prior to the start of the experiment, which was approved by the Research Ethics Committee at University College London (UK), under ethics number 9929/002.
Version history
 Received: January 26, 2021
 Accepted: June 6, 2021
 Accepted Manuscript published: June 7, 2021 (version 1)
 Version of Record published: July 28, 2021 (version 2)
Copyright
© 2021, Liu et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 4,473
 views

 836
 downloads

 27
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
 Neuroscience
Normal aging leads to myelin alterations in the rhesus monkey dorsolateral prefrontal cortex (dlPFC), which are positively correlated with degree of cognitive impairment. It is hypothesized that remyelination with shorter and thinner myelin sheaths partially compensates for myelin degradation, but computational modeling has not yet explored these two phenomena together systematically. Here, we used a twopronged modeling approach to determine how agerelated myelin changes affect a core cognitive function: spatial working memory. First, we built a multicompartment pyramidal neuron model fit to monkey dlPFC empirical data, with an axon including myelinated segments having paranodes, juxtaparanodes, internodes, and tight junctions. This model was used to quantify conduction velocity (CV) changes and action potential (AP) failures after demyelination and subsequent remyelination. Next, we incorporated the single neuron results into a spiking neural network model of working memory. While complete remyelination nearly recovered axonal transmission and network function to unperturbed levels, our models predict that biologically plausible levels of myelin dystrophy, if uncompensated by other factors, can account for substantial working memory impairment with aging. The present computational study unites empirical data from ultrastructure up to behavior during normal aging, and has broader implications for many demyelinating conditions, such as multiple sclerosis or schizophrenia.

 Neuroscience
How is new information organized in memory? According to latent state theories, this is determined by the level of surprise, or prediction error, generated by the new information: a small prediction error leads to the updating of existing memory, large prediction error leads to encoding of a new memory. We tested this idea using a protocol in which rats were first conditioned to fear a stimulus paired with shock. The stimulus was then gradually extinguished by progressively reducing the shock intensity until the stimulus was presented alone. Consistent with latent state theories, this gradual extinction protocol (small prediction errors) was better than standard extinction (large prediction errors) in producing longterm suppression of fear responses, and the benefit of gradual extinction was due to updating of the conditioning memory with information about extinction. Thus, prediction error determines how new information is organized in memory, and latent state theories adequately describe the ways in which this occurs.