Participants
The participants were 17 patients with intractable epilepsy who were implanted with depth electrodes to delineate a potentially surgically treatable epileptogenetic zone. Demographics information and neuropsychological scores are presented in Supplementary Table 1. Electrode placements were determined solely on the basis of clinical treatment criteria. The follow-up studies (Extended Data Figs. 1 and 10) included 33 healthy controls (26 female participants; mean age: 31 ± 7 years old) and 5 additional patients with epilepsy (2 female participants; mean age: 38 ± 12 years old). All participants volunteered for the study by providing informed consent according to a protocol approved by the UCLA Medical Institutional Review Board (IRB).
Neural recordings
Patients were stereotactically implanted with 7–12 Behnke-Fried electrodes with 40-µm diameter microwire extensions (eight high-impedance recording wires and one low-impedance reference wire per depth electrode) that capture local field potentials and extracellular spike waveforms46. Microwire electrophysiology data were amplified and recorded at 30 kHz on a Blackrock Microsystems recording system or at 32 kHz on a Neuralynx recording system (Cheetah 5.0).
Microelectrode localizations
Prior to data collection, each microelectrode location was confirmed by an expert neurosurgeon (I.F.) based on the patient’s postoperative computed tomography (CT) scan with visible electrode artifacts overlaid on a co-registered preoperative T1 structural MRI (BrainLab software). For descriptive purposes (Fig. 1d), we additionally used the following procedure to transform locations from each participant’s ‘native brain space’ to the standard Montreal Neurological Institute (MNI) space. First, each participant’s MRI and CT images were co-registered using the FSL ‘flirt’ function. Second, the MRI image was: (1) segmented into the grey matter, white matter, and cerebrospinal fluid probability maps; (2) resampled (1 × 1 × 1 mm voxel size); and (3) normalized to the 152 T1-weighted MNI template using the nonlinear transformation algorithm implemented in the Statistical Parametric Mapping toolbox (SPM12, Wellcome Department of Cognitive Neurology, London, UK). Third, the same transformation parameters were applied to the participant’s CT image. MNI coordinates for each microelectrode were extracted manually from the normalized CT overlaid on the normalized MRI from a given participant using the FSLeyes software.
General procedure
Before the main experiment (typically 1–2 days prior), a screening experiment was conducted to find 6 stimuli (images of people) associated with robust and preferential responses of single neurons in the MTL. These six images were then used during the main experimental task (Fig. 1b), which was introduced to the patients as a follow-up of the screening study without mentioning that the stimuli would be presented in a specific order. At the end of the main experiment, we asked the participants to answer the following questions: “Have you noticed any pattern in the sequence of images shown in any of the phases? If yes, what was it?”; “Have you had any special strategy during this study?”. None of the participants reported noticing any pattern that was relevant to the experimental manipulation (Fig. 1b). We used the Psychophysics Toolbox to control the timings of stimuli presentation and register behavioural responses47.
Screening session
During screening, approximately 120 images were repeatedly shown to the patients on a laptop computer (taking around 40 min). These images showed people, animals, objects and landmarks that were partly selected based on the participant’s preferences (for example, favourite actors, musicians, places, etc.). The experiment consisted of eight blocks, each with a different instruction (for example, block 1: “Determine whether each image shows a person or not”; block 2: “Determine whether each image shows a plant or not”; etc.). Each image was presented exactly once during each block, for the duration of 1 s, against a black background. The order of stimuli presentation was random. Participants indicated their responses using two assigned keys on a hand-held game pad.
Experimental task
The main study consisted of three parts: pre-exposure (PRE), exposure (E1–E6), and post-exposure (POST; Fig. 1c). During PRE (121 stimuli presented), all images were displayed in a pseudo random sequence (60 direct and 60 indirect graph-transitions; on average, each direct transition was presented 7 times and each indirect transition 9 times). The task was to determine whether each image showed a male or female (gender task). The participants used the right and left arrow keys on a laptop keyboard to indicate their responses. During the six subsequent exposure phases (121 trials in each phase), the order of stimuli was still randomized but restricted by the topological structure of the pyramid graph (Fig. 1c) so that only images directly linked on the graph were shown immediately after another. The starting location was selected randomly in each experiment phase. The behavioural task during all exposure phases was to determine whether a given image was mirrored or not when compared to PRE (Fig. 1c; the participants used the right and left arrow keys on a laptop keyboard to indicate their responses). During each phase, 61 images were ‘normal’ and 60 were ‘mirrored’. The order of mirrored and normal images was random. The POST phase was the same as PRE (all stimuli presented in a pseudo random sequence, without the ‘pyramid rule’; on average, each direct transition was presented 7 times and each indirect transition 8 times). Behavioural instructions displayed in the beginning of each phase emphasized that the participants should try to respond as quickly and accurately as possible. The first trial in each phase (that is, the beginning of a sequence) was discarded from the analyses, so effectively each phase consisted of 120 trials. The experiment in all phases was self-paced, that is: (1) a given image was displayed for as long as it took the participant to respond; and (2) the participants could have had breaks between phases for as long as they needed. All stimuli were displayed against a grey background. During a randomized inter-trial interval (1-3 s), a black ‘fixation’ circle was displayed in the middle of the screen. After each stimulus presentation, the participants received feedback (“correct!” or “incorrect” in relation to the currently performed task) displayed for 500 ms. All trials (correct and incorrect) were included in the analysis of electrophysiological data, as the behavioural tasks were unrelated to the main research question. Behavioural accuracy of responses during PRE and POST was near-perfect indicating that the gender task was easy for all the participants (Extended Data Fig. 1a). Accuracy in the ‘mirror task’ was lower but improved over the course of the study (Extended Data Fig. 1a; this task was supposed to be more challenging to maintain the participants’ attention).
Spike sorting
Automated spike detection and sorting were performed using the WaveClus3 software package in MATLAB48. We then manually reviewed each unit for inclusion by evaluating the waveform’s shape, amplitude, inter-spike intervals, and firing consistency across study phases. We rejected units that were likely contaminated by artifacts, in keeping with field-standard spike evaluation criteria49. For electrodes with multiple putative units that passed this inclusion check, we merged units whose waveform features could not be well-separated in principal components space, retaining for analysis a combination of single- and multi-units.
Single-neuron analyses
For each neuron and each stimulus presentation, we selected a time window around the stimulus onset (from −1 to +2 s). Then we calculated the number of spikes in 0.1 s time bins, smoothed (moving sum: ± 0.25 s) and baseline-corrected the data (subtracted the mean activity in the −0.5 to 0 s time window). The ‘response window’ was defined from 0.1 to 1.2 s after the stimulus onset. For a given neuron, the ‘preferred stimulus’ was the image associated with the strongest mean response in the response window during PRE. Depending on the position of the preferred stimulus on the pyramid, the remaining images were labelled as ‘direct’ or ‘indirect’ (Fig. 1b). This assignment was used across all study phases. ‘Selective neurons’ were defined as cells that during PRE: (1) responded significantly stronger to the preferred stimulus in the response window versus baseline; and (2) responded significantly stronger to the preferred stimulus than to the remaining stimuli combined. ‘Relational neurons’ were defined as cells that: (1) were selective (see above); (2) responded significantly stronger to the direct than indirect stimuli during E5 and E6; and (3) responded significantly stronger to direct stimuli during E5 and E6 than during E1 and E2. The ‘diminishing selectivity neurons’ were defined as cells that: (1) were selective; and (2) responded significantly weaker to the preferred stimulus in E5 and E6 than in E1 and E2. All the above criteria were tested with the Wilcoxon signed-rank tests (one-sided) with a P value threshold of 0.05. The above procedure was repeated 1,000 times, with random permutations of the stimulus or phase labels, depending on which criterion was tested. These permutations informed how many neurons of a given type are expected in a given brain region by chance. The empirical P value was calculated as the number of permutations with more neurons of a given type than the number of neurons actually detected, divided by the total number of permutations. If this value was less than 0.05, we concluded that a given brain region contained a significant proportion of a given neuron type (Extended Data Fig. 2 and Supplementary Table 3). To analyse combined responses of all relational neurons (Extended Data Fig. 4), we calculated the difference between each neuron’s mean responses to direct minus indirect stimuli and preferred minus non-preferred stimuli. This was done for each study phase separately. Then, we peak-normalized and baseline-corrected (−0.5 to 0 s) those differences and extracted the mean from the 0.1 to 1 s time window. For line plots showing the mean responses of individual neurons (Fig. 2a,b and Extended Data Figs. 3a and 9b), we used 0.01 s bins and the ±0.25 s moving sum. For plots showing multiple neurons (Figs. 1f and 2c), we z-scored and baseline-corrected (−0.5 to 0 s) the data from each neuron (for illustration purposes, we used the ± 0.2 s moving sum and heat maps were additionally smoothed with ± 0.1 s moving average).
Neural population analyses
To decode stimulus identity during each image presentation, we used the Poisson naive Bayes classifier, as implemented in the Neural Decoding Toolbox50. The spiking activity of each neuron was extracted from the −1 to +2 s time window relative to the stimulus onset. Data was binned (0.1 s) and smoothed (moving sum: ±0.25 s). The decoder was run on the summed spiking activity in the 0.1 to 1 s time window (Extended Data Fig. 5). The main analysis focused on posterior probabilities assigned by the decoder to the image actually presented (actual), images directly linked to that stimulus on the graph (direct), and images linked indirectly (indirect). The analysis was performed for each recording session separately (different stimuli), but the resulting probability distributions were combined across all sessions and image presentations. The classifier was trained on the data from PRE and tested on all subsequent phases. For testing in PRE, we used the ‘leave-one-trial-out’ cross-validation. Kolmogorov–Smirnov tests were used to compare cumulative distribution functions (CDFs) of posterior probabilities (one-sided). To reconstruct the entire pyramid graph (Fig. 3f,g and Extended Data Figs. 6b and 7b), we calculated Euclidean distances between mean responses of each neuron to each pair of images across all relevant neurons (neurons that stopped firing during the late study phases were excluded; distances were z-scored; bin-size: 0.1 s; baseline-correction: −0.5 to 0 s; moving sum: ± 0.15 s; time window: 0.1 to 1 s). Then, we compared the resulting neural distance matrixes to three templates (Fig. 3f). In the geodesic template, distances between each pair of nodes corresponded to the number of edges of the shortest path connecting the nodes. In the Euclidean template, distances between nodes 1–5, 4–3 and 6–2 were calculated from the Pythagorean theorem (right triangles: 1–5–6, 4–3–1, 6–2–1). The remaining distances corresponded to the shortest path (see above). In line with the previous literature2,51, the successor template (ST) was calculated as the negative of the matrix exponential of the adjacency matrix A:
$${\rm{ST}}={-{\rm{e}}}^{A}=-\mathop{\sum }\limits_{n=0}^{\infty }\frac{{A}^{n}}{n!}$$
The above metric provided a slightly better fit to the data than a related index that defines the relationships between states (Extended Data Fig. 7b):
$$\mathop{\sum }\limits_{n=0}^{\infty }{\gamma }^{n}{A}^{n}={(I-\gamma A)}^{-1}$$
Here, entries aij for each An correspond to the number of possible paths of length n between objects i and j and a discount factor is 0 < γ < 1 (refs. 2,19). To illustrate most faithful 2D representations of the respective distance matrixes, we used the multidimensional scaling analysis (MDS; ‘mdscale’ function in MATLAB; criterion: ‘sammon’). Because MDS can only be performed on matrices with positive entries, we normalized the matrixes by adding the absolute value of the matrix’s minimum plus a constant of 0.1. The similarity between neural distance matrixes and each template was calculated as the Spearman correlation (Fisher-transformed). Because the aim was to test how this similarity changes over the course of the study (unconfounded by any potential pre-existing similarity), for each phase, we subtracted the degree of similarity in PRE (Fig. 3g and Extended Data Figs. 6b and 7b). To obtain null distributions of correlation coefficients, the above procedure was repeated 10,000 times with random permutations of the nodes’ positions. P values were calculated as the number of permutations with higher correlation coefficients than the one actually detected, divided by the total number of permutations. If in none of the permutations the correlation was above the actual value, the P < 0.0001 range is reported.
Replay analysis
We analysed sessions that contained at least three selective hippocampal–entorhinal neurons from the same hemisphere, whose preferred stimuli from PRE mapped onto three-element pyramid trajectories (Fig. 4c, one seed neuron, one direct neuron and one indirect neuron forming a connected path). For each spike of the seed neuron, we checked whether the direct and indirect neurons fired at least once in the 0 to 30 ms time window. The above situation had to occur at least five times to be included in the analysis (that is, n < 5 was considered insufficient for robust statistical inference). There were 536 such putative replays in B1 (break after PRE) and 1,012 in B2–B7 (breaks after exposure phases). If the direct neuron fired significantly earlier than the indirect neuron, the replay was labelled ‘congruent’ (Fig. 4c; we analysed latencies of the first ‘direct spikes’ versus latencies of the first ‘indirect spikes’ across all repetitions of a given replay52; Wilcoxon signed-rank test, one-sided, with a P value threshold of 0.05). If the opposite was true, a replay was labelled ‘incongruent’. To obtain P values for the comparisons between proportions of congruent and incongruent replays throughout the study, we randomly shuffled spikes from the direct and indirect neurons (1,000 permutations) and used the resulting null distribution as reference.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.