Reduced coupling between offline neural replay events and default mode network activation in schizophrenia

Abstract Schizophrenia is characterized by an abnormal resting state and default mode network brain activity. However, despite intense study, the mechanisms linking default mode network dynamics to neural computation remain elusive. During rest, sequential hippocampal reactivations, known as ‘replay’, are played out within default mode network activation windows, highlighting a potential role of replay-default mode network coupling in memory consolidation and model-based mental simulation. Here, we test a hypothesis of reduced replay-default mode network coupling in schizophrenia, using magnetoencephalography and a non-spatial sequence learning task designed to elicit off-task (i.e. resting state) neural replay. Participants with a diagnosis of schizophrenia (n = 28, mean age 28.2 years, range 20–40, 6 females, 13 not taking antipsychotic medication) and non-clinical control participants (n = 29, mean age 28.1 years, range 18–45, 6 females, matched at group level for age, intelligence quotient, gender, years in education and working memory) underwent a magnetoencephalography scan both during task completion and during a post-task resting state session. We used neural decoding to infer the time course of default mode network activation (time-delay embedding hidden Markov model) and spontaneous neural replay (temporally delayed linear modelling) in resting state magnetoencephalography data. Using multiple regression, we then quantified the extent to which default mode network activation was uniquely predicted by replay events that recapitulated the learned task sequences (i.e. ‘task-relevant’ replay-default mode network coupling). In control participants, replay-default mode network coupling was augmented following sequence learning, an augmentation that was specific for replay of task-relevant (i.e. learned) state transitions. This task-relevant replay-default mode network coupling effect was significantly reduced in schizophrenia (t(52) = 3.93, P = 0.018). Task-relevant replay-default mode network coupling predicted memory maintenance of learned sequences (ρ(52) = 0.31, P = 0.02). Importantly, reduced task-relevant replay-default mode network coupling in schizophrenia was not explained by differential replay or altered default mode network dynamics between groups nor by reference to antipsychotic exposure. Finally, task-relevant replay-default mode network coupling during rest correlated with stimulus-evoked default mode network modulation as measured in a separate task session. In the context of a proposed functional role of replay-default mode network coupling, our findings shed light on the functional significance of default mode network abnormalities in schizophrenia and provide for a consilience between task-based and resting state default mode network findings in this disorder.

Participants were classified as antipsychotic-free if they had been free from antipsychotic treatment for at least 6 weeks for oral or 6 months for depot formulations, consistent with previous approaches 2 . Patients were either taking antipsychotic medication (n = 15, see Supplementary Table 1 for details) or not taking any medication at all. Medication was not stopped for this study for those patients taking medication. Healthy volunteers were not taking any medication.
We assessed psychiatric symptoms with the Positive and Negative Syndrome Scale (PANSS) scale 3 , Montgomery Åsberg Depression Rating Scale (MADRS) 4 , and General Assessment of Function (GAF) 5 . We administered brief measures of IQ (the Wechsler Test of Adult Reading, WTAR 6 ) and working memory (mean of forward and backward Digit Span).

Applied learning task and MEG sessions (Dataset A)
The following task description is taken from our published report using the same dataset (Nour et al, (2021)), in which we focused on detection of neural replay and representational similarity structure 7 . Participants performed a validated task during MEG, where they inferred how 8 distinct task pictures were embedded within two separate 'structural sequences' ([AàBàCàD] & [A'àB'àC'àD']), without ever being shown either complete sequence order 8 . Instead, participants viewed scrambled 'visual sequences' containing pictures from both 'structural sequences'. During pre-scan training (visit 1) participants were explicitly shown how 8 training pictures embedded within 3 'visual sequences' could be mapped onto 2 correct 'structural sequences'. Specifically, they were informed that in each 'visual sequence' (e.g., [CàDàC'àD']) only the first and last transitions correspond to correct 'structural relationships' (e.g., [CàD & C'àD'] are 'structural' transitions, but [DàC'] is a 'visualonly' transition), and their understanding was checked subsequently by asking them to explain this 'unscrambling rule'. Once this understanding was established, participants completed 3 Structure Learning sessions in which they practiced how to infer the correct structural relationships between the 8 training pictures; learning not only how the specific training pictures were associated, but also encoding a stimulus-independent mapping from 'visual' to 'structural' sequences (i.e., the 'unscrambling rule'). During MEG (visit 2) participants completed 3 similar Applied Learning sessions, with an entirely new picture set of 8 pictures, thus applying the encoded task schema to novel stimuli.
In both the pre-scan training session (visit 1), and MEG Applied Learning sessions (visit 2) participants were presented with a set of 3 unique 'visual sequences' (e.g., Following each Applied Learning session, we assessed knowledge of structural sequences (i.e., the pictureàstate embeddings) with a 12-question quiz. For each question a single target picture (e.g., B) was presented on the center of the screen for 5 s, followed by the appearance of two probe pictures (e.g., A and D) at the bottom left and right of the screen.
Participants indicated which of the two probe pictures comes later than the target picture in the 'structural sequence' with a button press, with no time restrictions or feedback. Laterality of the correct/incorrect pictures was randomly selected on each trial. Patient and control participants were matched in quiz performance at the end of the final Applied Learning session (immediately prior to a 5 minute post-learning rest session) and also immediately after this rest session. However, patients showed impaired knowledge of inferred sequences after the first and second Applied Learning session (a complete analysis of behavioral results from this dataset is presented in our prior published report, Nour et al, Cell (2021) 7 ).
During MEG, prior to Applied Learning, participants completed a 5 minute rest session with eyes open (PRE, at the start of the scanning session). This was followed by a Stimulus Localizer task, which provided training data for the stimulus decoders required for sequenceness analysis. On each trial of the Stimulus Localizer participants were presented a task picture on the center of the screen (320*320 pixels, 1 s), followed immediately by presentation of a single-word text descriptor, and indicated whether the text matched (50% trials) or did not match (50% trials) the preceding picture with a button press (2 s response window, no feedback, laterality of yes/no response counterbalanced between subjects). Each picture was shown many times (mean 51.0, SD 3.48) in a randomized order. The inter-trial interval was drawn uniformly from 700 -1300 ms on each trial. Following Applied Learning, participants completed a second 5 minute rest session with eyes open (POST), which ended with a 4 th (post-rest) knowledge quiz to assess knowledge retention. At the end of the scan session participants completed a Position Probe task. This differed from the Stimulus Localizer task in one respect alone. Here, each picture was followed by a single number (1, 2, 3 or 4), and participants now indicated whether this number matched the position of the preceding picture with a button press (2 s response window, no feedback, chance accuracy 50%, laterality of yes/no response counterbalanced between subjects). Pictures were presented in randomized order.
The task was implemented in MATLAB (MathWorks) using Cogent (Wellcome Trust Centre for Neuroimaging, University College London) v 1.30.

MEG acquisition
For Dataset A, 3 sensors were not recorded due to excessive noise in routine testing, and the task was projected onto a screen suspended in front of participants. Participants responded using two buttons (L/R) of a MEG-compatible button box (Current Designs) held in the right hand. For Dataset B, MRI data (used for MEG coregistration) were acquired using a Philips Achieva 7T scanner.

MEG Data preprocessing
For the replay detection pipeline (Dataset A), we filtered sensor space data to a pass band of 0.5 to 50 Hz (including a notch filter for power line artefact at 48 -52 Hz) and downsampled to 100 Hz. For the RSN analysis pipeline, we filtered sensor space data to a pass band of 1 to 45 Hz and downsampled data to 250 Hz. In both pipelines excessively noisy MEG data segments (rest data), trials (stimulus localizer data), and sensors were removed, followed by automated detection and removal of artefactual MEG components identified using independent component analysis (ICA) As outlined in our previous study 9 , the use of a slightly amended filter pass band for RSN analysis (1 -45 Hz) ensured RSN-state dynamics were not driven by low frequency sensor drift effects or mainline power noise effects, which the HMM modeling approach is more sensitive to compared to the replay identification methods. Similarly, the higher sampling rate used for RSN analysis pipeline (250 Hz) ensured sufficient resolution of the time embedding to enable good estimation of spectral content for each RSN-state definition.
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 (e.g., eye blink and mains interference) were classified by automated inspection of the spatial topography, time course, kurtosis of the time course and frequency spectrum for all components 8 . Artefacts were rejected by subtracting them out of the data.

Replay analysis pipeline
In our original report we used a Temporally Delayed Linear Modelling (TDLM) framework 10 to quantify extent to which correct (inferred) transitions between task states were expressed in state reactivation time courses (e.g., transition [AàB] is expected to manifest in spontaneous neural reactivation of state A, ( ! ), reliably preceding that of state B, ( " )).
We considered such sequential (time-lagged) reactivation patterns at different time lags, (corresponding to 'neural replay' of inferred structure at different replay speeds), and in a manner that controls for non-specific lagged reactivation patterns and the effect of a background alpha oscillation in MEG data 7 . We found evidence for significant replay of inferred structure in a predominantly forward direction, with maximal evidence (across the sample of patients and controls) at = 40 ms lag 7 , a lag corresponding also to findings from earlier reports in healthy volunteers 8,11 .
In the present study we used this empirically identified replay lag ( = 40 ms) to estimate a separate 'replay probability' time course for each unique task state transition ( This probabilistic output was thresholded at the transition-specific 99 th percentile 9 to provide the time course of 'replay onsets' used throughout this paper.

Source space reconstruction
The RSN analysis pipeline was conducted in source space.  , )).
The transition probability between discrete latent states is governed by a [ , ] transition probability matrix, , which respects the first-order Markov conditional independence property:  The form of the observation model, is a key analysis choice, and may be tailored to the data features best captured by the imaging modality and of most relevance to the experimental question 15 . In the present work we use a HMM with Time Delay Embedding (HMM-TDE) approach, in which the observation model associated with each RSN state reflects the distribution of spectrally-resolved power at each ROI, and coherence between ROIs. This reflects an understanding that spectrally-resolved phase coupling between brain regions is a key mechanism of cross-region communication ('functional connectivity') in the brain 15,16 . The HMM-TDE observation model augments ('temporally embeds') the observation at each time point t with MEG data from 2 adjacent time points before and after , as follows: ,

Further discussion of epoching RSN activation dynamics by replay onsets
As described in main text, for each participant, rest-session, and replay transition, we identified timepoints where 'replay evidence' exceeded the 99 th percentile (where replay evidence is defined as reactivation probability of 1 st task picture multiplied by the lagged reactivation probability of 2 nd task picture). We then used the identified timepoints (replay onsets) to epoch the RSN activation time course and calculated the mean RSN activation epoch for each replay transition in turn (n = 64). Prior to averaging the individual epochs for each transition, we first excluded epochs which contained (transition-specific) suprathreshold replay onsets prior to 0 ms. This ensured that we only included RSN activation epochs that contained an 'event-free baseline' (for the transition in question). Note that we had no similar restriction for time points after 0 ms (such that 0 ms could represent the first of multiple suprathreshold events). We did not include epochs where this event-free baseline overlapped with the previous (transition specific) epoch.

Further discussion of transition-specific epoching procedure
Our transition-specific epoching approach has two key advantages. First Figure 3, wherein we additionally calculate the change in such an effect from pre-to post-learning rest sessions, ∆ -678998: . For this regression analysis we first downsampled the epoched RSN activation time course to 100 Hz for computational expediency). A second advantage of our approach is that by using transition-specific thresholding we ensure that our summary estimates of replay-RSN coupling over transitions are not unduly influenced by differences in output magnitude between different neural decoders (within participants).

Extracting primary axis of variation in evoked RSN dynamics
The primary approach we adopt to analyzing evoked RSN dynamics (illustrated in Figure 3B) is well suited to analyzing activation effects arising at individual time*RSN

Dynamical properties of replay onsets
We present an analysis of replay interval times, conditional on DMN activation, in

Statistical analysis and software
For all analysis we defined outlier participants as those exhibiting effect sizes >3.5 SD ± group median. For replay-evoked MEG analysis, this corresponded to 2 control participants (2 for structure learning performance) and 1 patient (sequenceness effect size meeting outlier criteria), identically to our original replay paper 7 . For analyses restricted to stimulus localizer sessions, where knowledge of task structure was irrelevant, we excluded 2 participants (1 patient, 1 control) for outlier behavioral performance during the concurrent attention task.
We used non-parametric permutation tests to assess the statistical significance of effects where we wished to control for multiple comparisons across time points (i.e., evoked RSN dynamics) or frequencies (i.e., RSN spectral properties). This procedure involved repeatedly computing the effect of interest at each time or frequency sample for each of 500 participantlevel permutations where the nature of the permutation instantiated the relevant null hypothesis.
To test for significant effects greater than 0 (i.e., one-sample hypotheses, e.g., replay-evoked DMN activation > 0 shown in Figure 3B & 3C) in each permutation we flipped the sign of each participant's effect time course with 50% probability (sign-flip permutation). To test for significant differences in these effects between patients and controls (i.e., two-sample hypotheses, as in Figure 3D) in each permutation we randomly assigned 27 participants as 'controls', and the remaining participants as 'patients' (group-membership permutation). For both one-and two-sample tests, for each permutation we then identified the maximal grouplevel effect of interest (for peak-level effects this corresponded to the maximal effect over all time or frequency samples, while for cluster-level this corresponded to the cluster exhibiting the greatest magnitude, as described below). The distribution of maximal effects over all permutations thus defined an empirical null distribution of effect sizes that controlled for multiple comparisons over samples. Any effect in the true (unpermuted) data that exceeded the 95 th percentile of this distribution was deemed statistically significant at PFWE < 0.05. We used a peak-level significance threshold for our primary analysis involving replay-evoked DMN activation for inferred transitions (Figure 3), reflecting the temporally-localized nature of these effects. For one-sample hypotheses the effect of interest was the mean DMN activation over participants (right-tailed hypothesis, Figure 3B & C), while for group-difference hypotheses we used the absolute t value derived from a two-sample t-test (two-tailed hypothesis, Figure   3D). We used an analogous cluster-level significance approach when assessing significance of effects that manifested over a broader range of adjacent time points or frequencies (i.e., non-

Number of replay onsets detected per transition and across groups
For each participant and session (e.g., rest), we defined 'replay onsets' as time points in the replay-probability time course that exhibited high (> 99 th percentile 9 ) evidence of (transition-specific) sequential task state reactivation, preceded by a pre-onset baseline of subthreshold replay evidence (see Supplementary Materials and Methods). Using this criterion, we detected (mean) 87.9 ± SEM 1.91 individual 'replay onsets' per replay transition, across participants. Importantly, there was no difference between patients and controls in the number of identified replay onsets for any transition type (n = 64 two sample t-tests, one for each replay transition type, all P > 0.05, two tailed).

Replay-DMN in pre and post-learning rest sessions separately
The

Coupling between RSN activation and non-specific task state replay
Our primary results use a multiple regression approach to estimate the replay-DMN coupling effect that is specific for task-relevant replay transitions (i.e., that which exists over and above any 'replay'-DMN coupling that might exist for replay of state transitions that do not correspond to a learned task structure). For completeness, here we present the 'non-specific'      The finding that numerically adjacent states exhibit elevated transition probabilities is expected as the numerical assignment of each state is based on the empirical transition matrix derived from the model fitting procedure on Dataset B, as previously described 9 . There was no group difference in the entropy of the RSN transition probability distribution (a measure of the uniformity of RSN state 1-step transition patterns, excluding selftransitions. Controls = 22.86 ± 0.12, patients = 22.67 ± 0.12, t(52) = 1.43, P = 0.16, two sample t-test, two tailed).