Faster than thought: Detecting sub-second activation sequences with sequential fMRI pattern analysis

Neural computations are often anatomically localized and executed on sub-second time scales. Understanding the brain therefore requires methods that offer sufficient spatial and temporal resolution. This poses a particular challenge for the study of the human brain because non-invasive methods have either high temporal or spatial resolution, but not both. Here, we introduce a novel multivariate analysis method for conventional blood-oxygen-level dependent functional magnetic resonance imaging (BOLD fMRI) that allows to study sequentially activated neural patterns separated by less than 100 ms with anatomical precision. Human participants underwent fMRI and were presented with sequences of visual stimuli separated by 32 to 2048 ms. Probabilistic pattern classifiers were trained on fMRI data to detect the presence of image-specific activation patterns in early visual and ventral temporal cortex. The classifiers were then applied to data recorded during sequences of the same images presented at increasing speeds. Our results show that probabilistic classifier time courses allowed to detect neural representations and their order, even when images were separated by only 32 ms. Moreover, the frequency spectrum of the statistical sequentiality metric distinguished between sequence speeds on sub-second versus supra-second time scales. These results survived when data with high levels of noise and rare sequence events at unknown times were analyzed. Our method promises to lay the groundwork for novel investigations of fast neural computations in the human brain, such as hippocampal replay.


Introduction
Training fMRI pattern classifiers on slow events. In slow trials, participants repeatedly viewed in earlier TRs and activation strengths will be proportional to the ordering of events during the (c) Repetition trials were always fast (32 ms ISI) and contained two visual images of which either the first or second was repeated eight times (causing backward and forward interference, respectively). In both task conditions, participants were asked to detect the serial position of a cued target stimulus in a sequence and select the correct answer after a delay period without visual input. One sequence or repetition trial came after five slow trials. (d) Mean behavioral accuracy (in %; y-axis) in upside-down slow trials. (e) Mean behavioral accuracy in sequence trials (in %; y-axis) as a function of sequence speed (ISI, in ms; x-axis). (f) Mean behavioral accuracy in repetition trials (in %; y-axis) as a function of which sequence item was repeated (fwd = forward, bwd = backward condition). All error bars represent ±1 standard error of the mean (SEM). The horizontal dashed lines in (e) and (f) indicate the 50% chance level. Detecting sequentiality in fMRI patterns following fast and slow neural event sequences.

151
Our first major aim was to test detection of sequential order of fast neural events with fMRI. We 152 therefore investigated above-mentioned sequence trials in which participants viewed a series of five 153 unique images at different speeds (Fig. 1b)  17.60%, p < .001, corrected). Thus, item detection is impacted more by succeeding than preceding 269 activation patterns, leading to the increased dominance of the last item in sequence trials particularly 270 in the fast conditions (Fig. 3a). Importantly, however, both sequence elements still differed from 271 non-sequence items even under conditions of interference (forward: 7.76% and backward: 7.59%, 272 respectively, all ps < .001, corrected), indicating that sequence element detection remains possible 273 under such circumstances. Using data from all TRs revealed qualitatively similar significant effects 274 (p < .05 for all but one test after correction, see SI). Repeating all analyses using proportions of 275 decoded classes (the class with the maximum probability was considered decoded at every TR), 276 or considering all repetition trial conditions, also revealed qualitatively similar results). Thus, brief 277 events can be detected despite significant interference. 278 We next asked which implications these findings have for the observed pattern transitions [cf. Detecting sparse sequence events with lower signal-to-noise ratio (SNR). The results 292 above indicate that detection of fast sequences is possible if they are under experimental control.

293
In most applications of our method, however, this will not be the case. When detecting replay, for 294 instance, sequential events will occur spontaneously during a period of noise. We therefore next 295 assessed the usefulness of our method under such circumstances. 296 We first characterized the behavior of sequence detection metrics during periods of noise.  Fig. 5a).

307
As before, we next fitted regression coefficients through the classifier probabilities of the rest 308 data and, for comparison, to concatenated data from the 32 ms and 2048 ms sequence trials ( Fig.   309 5b-c). As predicted by our modelling approach (Fig. 2e), and shown in the previous section ( Fig.   310 3b), the time courses of regression coefficients in the sequence conditions were characterized by 311 rhythmic fluctuations whose frequency and amplitude differed between speed conditions (Fig. 5c).

312
To quantify the magnitude of this effect, we calculated frequency spectra of the time courses of 313 the regression coefficients in rest and concatenated sequence data ( Fig. 5d; using the Lomb-Scargle 314 method [e.g., 60] to account for potential artefacts due to data concatenation, see Methods). This 315 analysis revealed that frequency spectra of the sequence data differed from rest frequency spectra 316 in a manner that depended on the speed condition ( Fig. 5d- this led to a step-wise reduction of the inserted sequence signal from 80% to 0%, relative to the 329 SNR obtained in the experimental conditions reported above.

330
As expected, differences in above-mentioned standard deviation of the probability gradually 331 increased with both the SNR level and the number of inserted sequence events when either fast or 332 slow sequences were inserted (Fig. 5f). In our case this led significant differences to emerge with 333 one insert and an SNR reduced to 12.5% in both the fast and slow condition ( Fig. 5g; comparing 334 against 0, the expectation of no difference with a conventional false positive rate α of 5%; all ps 335 false discovery rate (FDR)-adjusted).

336
Importantly, the presence of sequence events was also reflected in the frequency spectrum of 337 the regression coefficients. Inserting fast event sequences into rest led to power increases in the 338 frequency range indicative of 32 ms events (∼ 0.17 Hz, Fig. 5f Participants were rewarded with 3 cents for each oddball (i.e., stimulus presented upside-down) 459 that was correctly identified (i.e., hit) and punished with a deduction of 3 cents for (incorrect) 460 responses (i.e., false alarms) on non-oddball trials (i.e., when stimuli were presented upright). In 461 case participants missed an oddball (i.e., miss), they also missed out on the reward. Auditory were randomly assigned to each participant following this pseudo-randomized procedure.

499
To investigate the influence of sequence presentation speed on the corresponding neural ac-500 tivation patterns, we systematically varied the ISI between consecutive stimuli in the sequence. Poisson distribution with λ = 1.9 and truncated to an interval from 1 to 5. Thus, across all trials, the targets appeared more often at the later compared to earlier positions of the sequence. This was 532 done to reduce the likelihood that participants stopped to process stimuli or diverted their attention 533 after they identified the position of the target object. The serial position of the alternative response 534 option was drawn from the same distribution as the serial position of the target. As for the oddball 535 trials, auditory feedback was presented immediately following a response. The coin sound indicated 536 a reward of 3 cents for correct responses, whereas the buzzer sound signaled incorrect or missed 537 responses (however, there was no deduction of 3 cents for incorrect responses or misses). Together, 538 participants could earn a maximum reward of 2.25 Euro in this task condition.

539
Repetition trials We included so-called repetition trials to investigate how decoding time course or 3) or later (7, 8, or 9). This was done to ensure that the task was reasonably easy to perform. for Human Development in Berlin, Germany. The scanning procedure was exactly the same for both 595 study sessions. For the functional scans, whole-brain images were acquired using a segmented  run, the task began after the acquisition of the first four volumes (i.e., after 5.00 s) to avoid partial 607 saturation effects and allow for scanner equilibrium. We also recorded two functional runs of restingstate fMRI data, one before and one after the task runs. Each resting-state run was about 5 minutes 609 in length, during which 233 functional volumes were acquired. After the functional task runs, two 610 short acquisitions with six volumes each were collected using the same sequence parameters as Conversion of data to the brain imaging data structure (BIDS) standard.  procedure for which data from seven task runs were used for training and data from the left-out run 708 (i.e., the eighth run) was used for testing. This procedure was repeated eight times so that each 709 task run served as the training set once. We trained an ensemble of five independent classifiers,

773
The resulting brain maps of voxel-specific t-values resulting from the estimation of the de-774 scribed t-contrast were then combined with an anatomical mask of occipito-temporal brain regions.

775
All participant-specific anatomical masks were created based on automated anatomical labeling of 776 to the whole-brain smoothing performed during preprocessing, voxel activation from brain regions 785 outside the anatomical mask but within the sphere of the smoothing kernel might have entered the 786 anatomical mask (thus, in principle, also including signal from surrounding non-gray-matter voxels).

787
Finally, we combined the t-maps derived in each cross-validation fold with the anatomical masks.

788
All voxels with t-values above or below a threshold of t = 3 (i.e., voxels with the most negative across all correctly answered, upright slow trials (Fig. 2a). The mean accuracy scores of all partici-798 pants were then compared to the chance baseline of 100%/5 = 20% using a one-sided one-sample 799 t-test, testing the a-priori hypothesis that classification accuracy would be higher than the chance we modelled this response function as a sine wave that was flattened after one cycle, scaled by an 819 amplitude and adjusted to baseline. The model was specified as follows: whereby A is the response amplitude (the peak deviation of the function from baseline), f is the We fitted the four model parameters (A, f , d and b) to the mean probabilistic classifier evidence Considering the case of two sine waves with identical frequency but differing by a temporal shift sin(2πf δ 2 ) can be written as sin(2π δ 2λ ), which illustrates that the term monotonically increases until 837 δ > λ 2 . Second, the frequency term has to be adapted as follows: The flattening of the sine waves 838 to the left implies that the difference becomes positive at 0 rather than δ 2 , thus undoing the phase 839 shift and stretching the wave by 1 2 δ TRs. The flattening on the right also leads to a lengthening of 840 the wave by an additional 1 2 δ TRs, since the difference becomes 0 at 2πf + 2πf δ, instead of only 841 2πf + 2πf δ 2 . Thus, the total wavelength has to be adjusted by a factor of δ TRs, and no phase 842 shift relative to the first response is expected. The difference function therefore has frequency instead of f , and Equation 4 becomes -2A sin(2πf δ 2 ) sin(2π f 1+f δ t). We can now apply Equation

844
3 to the fitted response function as follows wherebyf ,d,b andÂ indicate fitted parameters. 846 We determined the relevant TRs in the forward and backward periods for sequence trials by 847 calculating δ depending on the sequence speed (the ISI the model fitted the data reasonably, we inspected the fits of the sine wave response function for 853 each stimulus class and participant using individual parameters (Fig. S2).

854
Detecting sequentiality in fMRI patterns on sequence trials. In order to analyze the neural 855 activation patterns following the presentation of sequential visual stimuli for evidence of sequentiality,   . and corresponding classification probabilities that were assigned their sequential position within the 859 true sequence that was shown to participants on the corresponding trial.

860
The main question we asked for this analysis was to what extend we can infer the serial order

874
The slopes were then averaged at every TR separately for each participant and sequence speed 875 across data from all fifteen sequence trials (Fig. 3b). Here, if later events have a higher classification 876 probability compared to earlier events, the slope coefficient will be negative. In contrast, if earlier 877 events have a higher classification probability compared to later events, the slope coefficient will be [129]). 895 We hypothesized that sequential order information of fast neural events will translate into order 896 structure in the fMRI signal and successively decoded events in turn. Therefore, we analyzed 897 the fMRI data from sequence trials for evidence of sequentiality across consecutive measurements.

898
The analyses were restricted to the expected forward and backward periods which were adjusted , also adding an interaction term for the two effects. Again, 958 the model included both by-participant random intercepts and slopes (Fig. 4c) We then inserted 1 to 6 sequence events into the pre-task resting state period by blending TRs   Attentive processing of the visual stimuli was a prerequisite to study the evoked activation patterns 3 in visual and ventral temporal cortex. We therefore excluded all participants that performed below 4 chance on either or both the repetition and sequence trials of the task. To this end, we removed all 5 participants with a mean behavioral accuracy below the 50% chance level from all further analyses 6 (Fig. S1a). We compared the relative proportion of misses and false alarms for each of the eight 7 functional task runs in the experiment. To this end, we conducted a LME model with trial type 8 (miss, false alarm), session (first, second) and session run (run 1-4) as fixed effects and included 9 by-participant random intercepts and slopes. As shown in Fig. S1b, misses (M = 0.55%) consis-10 tently occurred more frequently than false alarms (M = 0.30%), F 1,501.00 = 4.1, p = .04, which 11 was consistent across task runs (no effects of session or run, ps ≤ .70). Our classification was per-12 formed using a leave-one-run-out approach. In order to examine whether the accuracy of behavioral 13 performance on slow trials was stable across all task runs of the study, we conducted a LME model 14 that included the eight task runs as the fixed effect of interest as well as random intercepts and 15 slopes for each participant. The results showed no effect of task run indicating that the accuracy of 16 behavioral performance was relatively stable across task runs, F 1,92.72 = 0.13, p = .72 (Fig. S1c). 17 We examined whether behavioral accuracy on sequence trials was influenced by either the sequence   Figure S4: Classification time courses on sequence trials. Time courses (in TRs from sequence onset; x-axis) of (a) mean linear regression coefficients (slope), (b) mean correlation coefficients (Kendall's τ ), (c) mean step size between probability-ordered within-TR events, and (d) mean decoded serial event position with maximum probability for each sequence presentation speed (in ms; panels / colors). Shaded areas represent ±1 SEM. The blue and red rectangles indicate forward and backward period, respectively. Red dots indicate significant differences from baseline (horizontal gray line at zero; all ps ≤ .05, FDR-corrected for 38 comparisons). 1 TR = 1.25 s.
As reported in the main text, we verified that the sequentiality effects observed on sequence 71 trials (Fig. 3b) are not only driven by the event with the maximum probability but that sequentiality 72 is also present if the event with the maximum probability is removed. Examining the mean slope 73 coefficients within the expected forward and backward period (adjusted by considering only four 74 sequence events) after removing the event with the maximum probability showed that we could still 75 find evidence for sequential ordering (Fig. S5a). Significant forward ordering in the forward period 76 was still evident at sequence speeds of 512 and 2048 ms (ts ≥ 3.99; ps ≤ .001, FDR-corrected;  Additional analyses of repetition trials 91 We conducted two additional analyses for the data on repetition trials. First, we analyzed the effect of 92 event duration (number of repetitions) on event probability in more detail by calculating the average 93 event probability for each event type (first, second, and averaged non-sequence) as a function of 94 event duration (number of repetitions). Importantly, while we focused only on the two repetition 95 conditions with the highest degree of interference before, we now also included the data from all 96 intermediate repetition trial types. As before, we averaged the probabilities for each serial event 97 type but this time as a function of how often each item type was repeated in any given trial. Then, 98 in order to test how likely we were in decoding each serial event type (first, second, non-sequence), 99 when each item was only shown briefly once, we conducted three independent pairwise two-sample   Period Regression slope f Figure S5: Effects of sequence item removal on sequentiality metrics. (a, c, e) Time courses (in TRs from sequence onset; x-axis) of mean slope coefficients of a linear regression between serial event position and classifier probability (y-axis) for each speed condition (in ms; colors) on sequence trials after removal of (a) the sequence item with the highest classification probability, (c) the first sequence item, (e) the last sequence item. (b, d, f) Mean slope coefficients (y-axis) as a function of time period (forward versus backward; x-axis) and sequence speed (in ms; colors) after removal of (b) the sequence item with the highest classification probability, (d) the first sequence item,   Figure S6: Effects of event duration (element repetition) Average probability (in %; y-axis) as a function of the number of item repetitions (i.e., total event duration), separately for event types (first, second, and out-of-sequence events; colors) based on data of all TRs.
We asked whether we would be more likely to decode items that were part of the sequence actually 141 shown to participants (within-sequence items) as compared to items not part of the sequence (out-142 of-sequence items). To this end, we assessed if the serial events 1 and 2 were more likely to be 143 decoded in the repetition trials than other events. As before, we identified the item with the highest 144 classifier probability at every TR of each trial and then calculated the relative frequency of each 145 item in the decoded sequence of events. These frequencies were then averaged separately for each 146 repetition condition across all trials and participants. Next, using paired t-tests, we performed two 147 statistical tests: First, we tested how well we were able to decode a single briefly presented item 148 in a 32 ms sequence compared to items that were not presented, when the item is followed by a 149 statistical representation that could mask its activation pattern (short → long trials). Second, we 150 tested how well we were able to decode a single briefly presented item (first serial event) in a 32 ms 151 sequence compared to items that were not part of the sequence, when the item (last serial event) 152 is followed by a random statistical signal, for example, during an ITI (long → short trials).  an interaction between event type and event duration, F 2,915 = 17.90, p < .001 (see Fig. 4d).

171
In order to further characterize the origin of this interaction, we also conceived a reduced model 172 that did not include the data from out-of-sequence events. The results of this reduced model again from the averaged out-of-sequence items (M = 9.75%, SD = 3.05%), t (39) = 9.36, p < .001 while 184 the average probability of the first event was also higher compared to the out-of-sequence items, correction).
We also analyzed the trial-wise proportion of transition types between consecutively decoded 188 events using data from all 13 TRs following stimulus onset. This analysis revealed that in the short 189 → long condition the mean trial-wise proportion of forward transitions (M = 6.50) was higher than  Post-hoc comparisons indicated that sequence items had a higher mean probability than out-of-206 sequence (9.55%) items (both ps < .001, Tukey-correction for three comparisons), while the second 207 (16.77%) and first (16.77%) within-sequence event type also differed (p = .01, Tukey-correction 208 for three comparisons). Repeating the analysis for the forward and backward interference conditions 209 using data from all TRs again revealed smaller but qualitatively similar effects, with a main effect 210 of event type (first, second, out-of-sequence), F 2,43.34 = 55.42, p < .001, an interaction between 211 event type and duration, F 2,105.00 = 37.72, p < .001, and no main effect of duration (number of 212 repetitions), F 1,35.70 = 0.08, p = .78 (see Fig. S8c). Post-hoc comparisons indicated that in the 213 forward interference condition the longer second event had a higher probability (19.20%) compared 214 to both the out-of-sequence (M = 9.74%) and the short, first event (M = 11.42%, ps < .001, 215 Tukey-correction for three comparisons). As reported in the main text, when using data from all 216 TRs, the short first event did not differ from the out-of-sequence events (p = .13, Tukey-correction 217 for three comparisons). In the backward interference condition, in contrast, there was no difference 218 between the long first (16.25%) and short second event (14.34%, p = .22, Tukey-correction for three 219 comparisons) but significant differences between both within-sequence items and the averaged out-220 of-sequence (9.36%) items (ps < .001, Tukey-correction for three comparisons). We also repeated