Dynamics of fMRI patterns reflect sub-second activation sequences and reveal replay in human visual cortex

Neural computations are often fast and anatomically localized. Yet, investigating such computations in humans is challenging because non-invasive methods have either high temporal or spatial resolution, but not both. Of particular relevance, fast neural replay is known to occur throughout the brain in a coordinated fashion about which little is known. We develop a multivariate analysis method for functional magnetic resonance imaging that makes it possible to study sequentially activated neural patterns separated by less than 100 ms with precise spatial resolution. Human participants viewed five images individually and sequentially with speeds up to 32 ms between items. Probabilistic pattern classifiers were trained on activation patterns in visual and ventrotemporal cortex during individual image trials. Applied to sequence trials, probabilistic classifier time courses allow the detection of neural representations and their order. Order detection remains possible at speeds up to 32 ms between items (plus 100 ms per item). The frequency spectrum of the sequentiality metric distinguishes between sub- versus supra-second sequences. Importantly, applied to resting-state data our method reveals fast replay of task-related stimuli in visual cortex. This indicates that non-hippocampal replay occurs even after tasks without memory requirements and shows that our method can be used to detect such spontaneously occurring replay.

Reporting for specific materials, systems and methods automated anatomical labeling of brain surface reconstructions from T1-weighted reference images created with Freesurfer's recon-all (Dale et al., 1999, NeuroImage) as part of the fMRIPrep workflow. Main statistical analyses were conducted using LME models employing the lmer function of the lme4 package (version 1.1.21, Bates et al., 2015, Journal of Statistical Software) implemented in custom code in R (version 3.6.1, R Core Team, 2019). Model fitting was performed using NLoptr, an R interface to the NLopt library for nonlinear optimization (Johnson et al., 2019) employing the COBYLA (Constrained Optimization BY Linear Approximation) algorithm (Powell, 1994, Springer;Powell, 1998, Acta Numerica). An overview of all custom code used in the study is available on the project website https://wittkuhn.mpib.berlin/highspeed/ and the corresponding GitLab repository at https://git.mpib-berlin.mpg.de/wittkuhn/highspeed/. Please see the data availability statement below for the details on the custom code used to analyze the individual datasets.
We publicly share all code and data used in this study. An overview of all the resources is publicly available on our project website: https://wittkuhn.mpib.berlin/ highspeed/. The source code of the website is available at https://git.mpib-berlin.mpg.de/wittkuhn/highspeed. All individual datasets can be found at https://gin.gnode.org/lnnrtwttkhn. Please note, that each dataset is associated with a unique URL and Digital Object Identier (DOI). Data and code management was realized using DataLad [version 0.13.0; 142, for details, see https://www.datalad.org/]. Please note, that instead of separate data and code repositories, we share all data in modularized units alongside the code that created the data, usually in a dedicated code directory in each dataset. This approach allows to better establish the provenance of data (i.e., a better understanding which code produced which data), loosely following the DataLad YODA principles (for details, see the chapter "YODA: Best practices for data analyses in a dataset" in the DataLad handbook [version 0.13; 143], available at https://handbook.datalad.org/).
No statistical methods were used to predetermine the sample size but it was chosen to be larger than similar previous neuroimaging studies (e.g., Schuck & Niv, 2019, Science;Momennejad et al., 2018, eLife;Tambini et al., 2013, PNAS).
Four participants were excluded from further analysis because their mean behavioral performance was below the 50% chance level in either or both the sequence and repetition trials suggesting that they did not adequately process the visual stimuli used in the task. This exclusion criterion was pre-established based on previous piloting results of the behavioral task.
No direct replication of the experiment was performed. The current study is itself a conceptual replication of Schuck & Niv, 2019, Science.
For the task condition that showed sequences of five images we ensured that all possible sequences were chosen equally often across all participants. Given 120 possible sequential combinations in total and 15 sequences per participant, the sequences were distributed across eight groups of participants. Sequences were randomly assigned to each participant following this pseudo-randomized procedure.
Blinding is not relevant to our study as we treat data from all participants equally.

April 2020
We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response. The final sample consisted of 36 participants (mean age = 24.61 years, SD = 3.77 years, age range: 20 -35 years, 20 female, 16 male). All participants were screened for MRI eligibility during a telephone screening prior to participation and again at the beginning of each study session according to standard MRI safety guidelines (e.g., asking for metal implants, claustrophobia, etc.). None of the participants reported to have any major physical or mental health problems. All participants were required to be right-handed, to have corrected-to-normal vision, and to speak German fluently. Furthermore, only participants with a head circumference of 58 cm or less could be included in the study. This requirement was necessary as participants' heads had to fit the MRI head coil together with MRI-compatible headphones that were used during the experimental tasks.
Participants were recruited from an internal participant database or through local advertisement. Any potential self-selection bias, if present, cannot be explicitly ruled out since participants freely chose whether they wanted to participate and contact the experimenter based on the public advertisement and announcements sent through the participant database. These biases are, if present, unlikely to affect the results since the experiment was conducted in a within-subjects design (i.e., all participants experienced all conditions). The pseudo-randomized procedure used to assign all 120 possible stimulus sequences to participants (see above section on "Randomization") is also unlikely to interact with any self-selection bias, if present. The main effects investigated in this study (fMRI patterns of fast activation sequences) can be considered general and not specific to a population of young and healthy individuals with high education.
The ethics commission of the German Psychological Society (DGPs) approved the study protocol (Reference number: NS 012018).
The study consisted of two experimental sessions. In each session, we acquired four functional task runs of about 11 minutes during which participants performed the main task in an event-related design. We also recorded two functional runs of resting-state fMRI data, one before and one after the task runs.
In each session, we acquired four functional task runs of about 11 minutes during which participants performed the main task in an event-related task design. During each functional run, participants performed trials from three different task conditions. Slow trials started with a waiting period of 3.85 s during which a blank screen was presented. After a fixation dot for 300 ms, a stimulus was shown for 500 ms followed by a variable ISI during which a blank screen was presented again. The duration of the ISI was drawn from a truncated exponential distribution with a mean of 2.5 s and a lower limit of 1 s. Behavioral responses were collected during a fixed time period of 1.5 s after each stimulus onset. During sequence and repetition trials a target cue was shown for 1000 ms followed by a blank screen for 3850 ms. A short presentation of a gray fixation dot for 300 ms signaled the onset of the upcoming sequence of visual objects. All objects in the sequence were presented briefly for 100 ms. The ISI for each trial was determined based on the current sequence speed (32, 64, 128, 512, or 2048 ms) and was the same for all stimuli within a sequence. The sequence of stimuli was followed by a delay period with a gray fixation dot that was terminated once 16 s since the onset of the first sequence object had elapsed. Subsequently, the name of the target object as well as the response mapping were presented for 1.5 s. Slow trials were interleaved with sequence and repetition trials such that each of the 120 slow trials was followed by either one of the 75 sequence trials or 45 repetition trials. For resting state scans, participants were asked to stay awake and focus on a white fixation cross presented centrally on a black screen. Each resting-state run was about 5 minutes in length, during which 233 functional volumes were acquired.
Behavioral performance was assessed by measuring correct button presses (accuracy) and response times. Mean behavioral accuracy was used as the primary indicator to establish that the subjects were performing the task as expected. For the functional scans, whole-brain images were acquired using a segmented k-space and steady state T2*--weighted multiband (MB) echo-planar imaging (EPI) single-echo gradient sequence that is sensitive to the blood-oxygen-level dependent (BOLD) contrast (sequence specification: 64 slices in interleaved ascending order; anterior-to-posterior (A-P) phase encoding direction; repetition time (TR) = 1250 ms; echo time (TE) = 26 ms; voxel size = 2 x 2 x 2 mm; matrix = 96 x 96; field of view (FOV) = 192 x 192 mm; flip angle (FA) = 71 degrees; distance factor = 0%; MB acceleration factor 4). Slices were tilted for each participant by 15 degrees forwards relative to the rostro-caudal axis to improve the quality of fMRI signal from the hippocampus (cf. Weiskopf et al., 2006, NeuroImage) while preserving good coverage of occipitotemporal brain regions. These sequence parameters were the same for task and resting state acquisitions. After the functional task and resting-state runs, two short acquisitions with six volumes each were collected using the same sequence parameters as for the functional scans but with varying phase encoding polarities, resulting in pairs of images with distortions going in opposite directions between the two acquisitions which were used for distortion correction. In addition, a whole-brain spoiled gradient recalled (GR) field map with dual echo-time images ( Whole-brain images were acquired. , using brain-extracted versions of both T1-weighted volume and template. All classification analyses were performed on functional images that were co-registered to the individual T1-weighted reference. Here, a reference volume and its skull-stripped version were generated using a custom methodology of fMRIPrep. The BOLD reference was then co-registered to the T1-weighted reference using bbregister (FreeSurfer) which implements boundary-based registration (Greve et al., 2009, NeuroImage). Co-registration was configured with nine degrees of freedom to account for distortions remaining in the BOLD reference. We included the following nuisance regressors estimated during preprocessing with fMRIPrep in the first-level GLMs used to determine the functional ROIs for feature selection: the frame-wise displacement for each volume as a quantification of the estimated bulk-head motion, the six rigid-body motion-correction parameters estimated during realignment (three translation and rotation parameters, respectively), and six noise components calculated according to the anatomical variant of CompCorr (Esteban et al., 2018, Nature Methods).
No volume censoring was performed.
The main analyses were performed using multivariate leave-one-run-out cross-validated pattern classification. Feature selection was performed by combining functional ROIs based on thresholded t-maps with anatomical masks of predefined brain regions. Details about the multivariate classification analysis and the creation of functional ROIs can be found in the section "Multivariate modeling and predictive analysis" below. Details about the creation of anatomical ROIs can be found in the section "Anatomical location(s)" below.
The following effects were tested: Behavioral accuracy in the three task conditions compared to chance; cross-validated classification accuracy on slow trials compared to chance in occipito-temporal and hippocampal data; peak in classification probabilities of the true class on slow trials compared to all other classes and chance in occipito-temporal and hippocampal data; mean TR-wise linear regression slope coefficients as a function of sequence speed and time period compared to chance, Pearson's correlation between time courses of regression slopes predicted by modeling approach and based on data between and within participants; mean serial event position as a function of sequence speed and time period compared to