Predictive Feedback, Early Sensory Representations, and Fast Responses to Predicted Stimuli Depend on NMDA Receptors

Stimuli

For our psychophysics experiments, we used biomorphic visual stimuli from Michael Tarr's laboratory (https://sites.google.com/andrew.cmu.edu/tarrlab/resources/tarrlab-stimuli). These are known as greebles. Figure 1A shows examples presented to participants. We used three grayscale greebles for each psychophysics session, and each greeble was personified with a name. We used novel sounds (trisyllabic nonsense words) for the greeble names (e.g., “Tilado,” “Paluti,” and “Kagotu”) from Saffran et al. (1996). The sounds were generated using the Damayanti voice in the “text to speech” platform of an Apple MacBook. To avoid differences in the salience of stimuli, greeble images have similar size (13 degrees of visual angle in height and 8 degrees in width), number of extensions and mean contrast, and greeble names have the same number of syllables and sound level (80 dB SPL).

For the pharmacology experiments, we generated three new triplets of greebles. To control for saliency, each participant rated greeble salience for each of four (three plus one from the psychophysics experiment) triplets; that is, the participant identified whether any of the greebles in a triplet stand out compared with the other two greebles from the same triplet. We proceeded to use triplets which the participant rated all three greebles as being equally salient. We then named each of these greebles with a new trisyllabic nonsense word.

Audiovisual delayed match-to-sample task

Each trial of the task involves the sequential presentation of a sound (trisyllabic nonsense word) followed by a greeble image. We refer to stimuli using the following notation: A1, A2, and A3 correspond to each of the three sounds used (A for auditory); and V1, V2, and V3 correspond to each of the greebles used (V for visual). Using this notation, audiovisual stimulus sequences containing the matching name and greeble are A1-V1, A2-V2, and A3-V3. Audiovisual stimulus sequences containing a nonmatching name and greeble are A1-V2, A1-V3, A2-V1, A2-V3, A3-V1, and A3-V2. We pseudo-randomized names for greebles (i.e., matching sounds and images) across subjects.

Causal strength and transitional probability

We quantified the relationship between a sound (name) and its paired image (greeble) using the strength of causal induction: given a candidate cause C (sound) how likely is the effect E (i.e., how likely is it followed by its paired image). We will represent variables C and E with uppercase letters, and their instantiations with lowercase letters. Hence, C = c⁺/E = e⁺ indicates that the cause/effect is present, and C = c^–/E = e^– indicates that the cause/effect is absent (for brevity, we will shorten variables equal to outcomes, such as C = c⁺ or C = c^– as simply c⁺ or c^–, respectively). The evidence for a relationship can be encoded as a 2 × 2 contingency table for each sound, as in Table 1, where N(c⁺,e⁺) represents the number of trials in which the effect occurs in the presence of the cause, N(c^–,e⁺) represents the number of trials in which the effect occurs in the absence of the cause and so on. Applied to our study, for example, C could be hearing sound A1, and E viewing the paired greeble V1. For this case, N(c⁺,e⁺) would be the number of trials V1 follows A1, whereas N(c^–,e⁺) would be the number of trials V1 follows A2 or A3. The full contingency table for the HP auditory cue A1 and its paired greeble V1 is shown in Table 1. There are analogous contingency tables for the other two auditory cues and their paired greebles.

Based on these contingency tables, we calculated three different measures of causal relationship (ΔP, causal power and causal support) for each trial. ΔP and causal power assume that C causes E. ΔP reflects how the probability of E changes as a consequence of the occurrence of the cause C. Causal power corresponds to the probability that an effect E happened because of cause C in the absence of all other causes. Whereas causal support evaluates whether or not a causal relationship actually exists and calculates the strength of that relationship. To do this, causal support estimates the evidence for a graphical model with a link between C and E against one without a link. For example, let us consider the graphs denoted by Graph 0 and Graph 1 in Figure 2D (adapted from Griffiths and Tenenbaum, 2005). There are three variables in each graph: cause C, effect E, and background cause B. In Graph 0, B causes E, but C has no relationship to either B or E. In Graph 1, both B and C cause E. While calculating ΔP and causal power Graph 1 is assumed, whereas causal support compares the structure of Graph 1 to that of Graph 0. Causal support is defined as the evidence provided from data D in favor of Graph 1, P(D | Graph 1), over Graph 0, P(D | Graph 0), which can be calculated by the following equation: Causalsupport=logP(D|Graph 1)P(D|Graph 0)(1)

We calculated causal support using freely available MATLAB code from Griffiths and Tenenbaum (2005). ΔP and causal power were calculated using the following formulas: ΔP=N(c+,e+)N(c+,e+) + N(c+,e−)−N(c−,e+)N(c−,e+) + N(c−,e−)=P(e+|c+)−P(e+|c−)(2) Causalpower=ΔP1−P(e+|c−)(3)

We compared these three measures of causal relationship with the transitional probability (i.e., a comparison between causation and correlation). We calculated the transitional probability of each greeble (V) given prior presentation of its paired sound (A), using the following equation: Transitional probability=N(c+,e+)N(c+)(4)

For each subject, we used the causal relationship value of each condition (HP, MP, and NP) at the end of a testing phase as the starting values of the next testing phase. For example, the starting values of causal relationship for the “under drug” testing phase were equal to the causal relationship values at the end of the “pre/baseline” testing phase. Similarly, the starting values for the “after recovery” testing phase were equal to the causal relationship values at the end of the “under drug” phase (see Fig. 3A,B). We tested whether subjects (1) regained access to already learned and stored predictive information, after they recovered from ketamine dosing, or (2) relearned the predictive information. We mimicked Hypothesis 2 by forcing the starting values of causal relationship (causal power here) “after recovery” to be zero (see Fig. 3C) instead of starting values equal to the causal relationship at the end of the “under drug” phase (Hypothesis 1; see Fig. 3B).

Pharmacology experiments

To manipulate participants' predictions, we administered two drugs, ketamine and DEX, each on a separate day, with at least 1 month intervening. The fixed order of DEX in the first session, and ketamine in the second session, for the first 13 ketamine subjects was IRB-imposed in their consideration of safety profiles of the different medications (registered on NCT03284307); whereas the final 2 ketamine subjects were not administered DEX in a prior session. All subjects were healthy and between 18 and 40 years of age without contraindication to study drugs. We also acquired EEG data throughout pharmacology experiments (see EEG recording, below), to measure electrophysiological activity during PC. A typical pharmacology experiment consisted of three segments: (1) predrug baseline, (2) under drug influence (DEX or ketamine), and (3) after recovery. During the predrug baseline, participants performed the learning phase (200 trials), then the first testing phase (400 trials). Under drug influence, participants performed the second testing phase (200 trials). Under ketamine, participants also performed a third testing phase (200 trials; for details, see Ketamine dosing). After recovery, they performed the last testing phase (400 trials). As a drug control, we tested 3 participants during saline administration. The saline results were similar to the predrug baseline. Because of a protocol-limited maximum time under drug influence, and the need to acquire sufficient M/NM trials for EEG analyses, the pharmacology experiments did not include inversion trials. All other aspects of the task in pharmacology experiments were the same as that in psychophysics experiments. For each of the two pharmacology experiments involving a particular participant, we used three new greeble name and image pairs to rule out any possible contribution of long-term memory.

DEX dosing

We intravenously administered a 0.5 μg/kg bolus over 10 min, followed by 0.5 μg/kg/h infusion (MedFusion4000 pump at μg/kg/hr). Participants performed the testing phase (under drug influence) during this infusion time, corresponding to stable drug levels according to the pharmacokinetic model for DEX by Hannivoort et al. (2015). We targeted a plasma concentration of DEX that is associated with mild sedation (modified observer's assessment of alertness/sedation [OAA/S]) of 4) (Colin et al., 2017) to control for nonspecific sedative effects, including HCN-1 effects (Yang et al., 2014). The actual sedation achieved was on average slightly deeper than anticipated (modified OAA/S median 3, IQR 2) with a mean plasma concentration of 0.8 (SD 0.33) ng/ml (Colin et al., 2017).

Ketamine dosing

In initial experiments, we tested two doses of ketamine to target the lowest plasma concentration that would modulate NMDARs in the relevant concentration range (<1 μm) (Morris et al., 2017). The first dose corresponded to intravenously administered 0.25 mg/kg ketamine over 10 min, followed by 30 mg/h infusion, corresponding to 0.4 μm (Rugloop software using a Harvard 22 pump). Before testing subjects, lack of nystagmus and visual disturbance was confirmed in all participants. Participants performed the testing phase (under drug influence) during this infusion time, corresponding to stable drug levels according to the pharmacokinetic model for ketamine by Domino et al. (1982, 1984). A second dose was tested with a second bolus of 0.25 mg/kg ketamine over 10 min, followed again by 30 mg/h infusion. Testing again was completed once stable plasma concentrations of ∼0.8 μm were achieved. Ketamine blocked predictions at this second level of ketamine dosing, equating to a minimum plasma concentration of 0.2 µg/ml; range tested: 0.2-0.3 µg/ml. As we found the effective ketamine dose to be 0.2 µg/ml in our first 7 subjects, we targeted that plasma concentration for our remaining subjects. We report data for the 0.2 µg/ml dosing in the manuscript. Three subjects were excluded from “after recovery” testing because of vomiting.

Monitoring

Subjects were monitored during drug exposure according to the American Society of Anesthesiologists guidelines, including electrocardiogram, blood pressure, and oxygen saturation. We monitored arousal level using the modified OAA/S scale (Chernik et al., 1990). We also monitored subjects' experience of dissociation, using subjects' ratings on a dissociation scale ranging from 1 to 10. We asked subjects the question: “On a scale of 1-10, with 1 being completely normal and 10 being an out of body experience, how dissociated are you feeling right now?” The dissociation scores were similar for the predrug baseline (mean ± SD, 1.12 ± 0.35) and the low dose of ketamine (plasma concentration ∼0.2 μg/ml) used for the task (mean 2.65 ± 1.26); whereas dissociation scores were much higher (mean 8.58 ± 1.55) for a higher dose of ketamine (plasma concentration ∼0.4 μg/ml). After the end of the drug infusion, symptoms continued to be monitored. The “after recovery” testing occurred when nausea and nystagmus (resulting from the higher dosing) subsided sufficiently to perform the task. Subjects were generally able to perform the task within 2 h after stopping drug infusion.

EEG recording

We performed high-density EEG recordings using a 256 channel system (including NA 300 amplifier; Electrical Geodesics). After applying the EEG cap with conductive gel (ECI Electro-Gel), we adjusted electrodes so that the impedance of each electrode was within 0-50 kiloohms. We checked electrode impedance before the experiment started, and again before drug administration. Using Net Station, we sampled EEG signals at 250 Hz and, offline, bandpass filtered between 0.1 and 45 Hz.

EEG preprocessing

We combined predrug baseline data from both ketamine and DEX experiments (baseline RT data showed similar results for both drugs), but “under drug” analyses were performed separately for each drug. We performed offline preprocessing and analysis using EEGLAB (Delorme and Makeig, 2004). First, we extracted data epochs −1500 to 3000 ms relative to the onset of the sound and −3000 to 800 ms relative to the onset of the image, for each trial. We then visually inspected each epoch and excluded noisy trials (∼5% of the total trials). Next, we performed independent component analysis using built-in functions of EEGLAB (pop_runica.m) and removed noisy components through visual inspection. We excluded 3 DEX subjects from further analysis because of very noisy EEG data, which after cleaning left too few trials for analysis (conditions with <10 trials). Finally, we performed channel interpolation (EEGLAB function, eeg_interp.m, spherical interpolation) and rereferenced to the average reference.

Time-frequency decomposition

To investigate changes in EEG spectral content, we performed time-frequency decompositions of the preprocessed data in sliding windows of 550 ms using Morlet wavelets, whose frequency ranged from 5 to 45 Hz in 40 linearly spaced steps. Power for each time-frequency point is the absolute value of the resulting complex signal. We dB normalized power (dB power = 10*log10[power/baseline]) to the prestimulus baseline, starting 700 ms before sound onset (i.e., baseline calculated using sliding 550-ms-long wavelets, starting with the wavelet positioned from 700 to 150 ms before sound onset (centered on 425 ms before sound onset) and avoiding the sound-evoked response). For drift-diffusion model (DDM) analysis, we calculated power spectral density and performed divisive baseline correction for each trial.

Electrode selection

We used a data-driven approach, orthogonal to the effect of interest, to select the electrodes of interest based on the task EEG data. In the first step, we averaged power across all electrodes (aligned to image onset) and all sounds/greeble names. This revealed increased alpha power (8-14 Hz) during the delay period compared with baseline (see Fig. 5B). In the second step, we selected the electrode clusters that showed significant change in alpha power during the −925 to −275 ms time window pre-image onset, compared with baseline. This time window ensured that our delay period did not overlap with the sound or image. After cluster-based multiple comparisons correction (6), four different clusters showed significant modulation in alpha power during the delay period [and we selected one electrode (and its four nearest neighbors) in each cluster that showed the greatest modulation, for consistency across clusters; see Fig. 6A]: (1) right frontal (RF electrodes 4, 214, 215, 223, 224); (2) left central (LC electrodes 65, 66, 71, 76, 77); (3) right central (RC electrodes 163, 164, 173, 181, 182); and (4) occipital (OC electrodes 117, 118, 119, 127, 129). Additionally, we found that only alpha power at the RF cluster significantly correlated with the predictive value of sounds.

Alpha power calculation

For each condition, HP, MP, and NP, we averaged alpha power over all five electrodes in a cluster, to calculate the average alpha power of each cluster. To best capture the delay period activity just before image onset, we calculated mean alpha power between −625 and −275 ms (as the wavelet window is centered around each time point, the power estimate before −625 ms and after −275 ms may contain auditory and visual stimulus-related responses, respectively) for each trial aligned to image onset. We also calculated mean delay period alpha power between 600 and 1225 ms for each trial aligned to sound onset. To link EEG power spectral density to behavior using our DDM analysis (see HDDM), we calculated single-trial baseline-corrected (divisive normalization) alpha band power, aligned to image onset.

Source space analysis

We used FieldTrip's beamforming technique to localize sources of the sensor level alpha activity (Oostenveld et al., 2011). This technique uses an adaptive spatial filter to estimate activity at a given location in the brain. We used a source model defined in MNI space for all subjects. Across all sound cues, we used the Dynamic Imaging of Coherent Sources (Gross et al., 2001) algorithm to beamform the delay period alpha activity (window 0-900 ms before image onset). We then calculated the average source power for each ROI of the AAL atlas (Tzourio-Mazoyer et al., 2002). We selected ROIs that showed significant change in source power during the delay period (before image onset) compared with baseline (p < 0.025, Bonferroni-corrected across all cortical AAL ROIs).

Granger causality

We performed nonparametric spectral Granger causality analyses at the source level. Because there were no significant frontal ROIs in the left hemisphere, we restricted our analyses to the right hemisphere only. We used a covariance window of 0-900 ms before image onset and the Linearly Constrained Minimum Variance algorithm (Van Veen et al., 1997) to generate virtual time series for significant frontal and temporal ROIs. We averaged across the significant frontal ROIs (superior frontal gyrus, medial; superior frontal gyrus, dorsolateral; middle frontal gyrus); and there was only one significant temporal lobe ROI (inferior temporal gyrus). We calculated nonparametric spectral Granger causality between the (average) frontal and temporal ROIs using the ft_connectivityanalysis function in FieldTrip for the stable window starting 900 ms before image onset (Wang et al., 2008; Seth et al., 2015). We found Granger causal influence from RF to temporal cortex correlated with the predictive value of sounds only in the β frequency (15-30 Hz) band.

Hierarchical DDM (HDDM)

We used a DDM (Ratcliff and McKoon, 2008), where there are two possible choices (correct/incorrect responses of the match trials) in our PC task. According to this model, decision-making involves the accumulation of evidence (drift process) from a starting point to one of two (upper or lower) thresholds, representing the choices. The accumulation rate is known as the drift rate, v; and the starting point can be biased toward one of the choices (in our study, by the predictive value of the sound), reflected in a bias parameter, z. We used HDDM software (http://ski.clps.brown.edu/hddm_docs/) (Wiecki et al., 2013) for hierarchical Bayesian estimation of the parameters of the DDM. Particularly with fewer trials per condition, this method has been shown to provide more reliable estimates of parameters and is less susceptible to outliers (Ratcliff and McKoon, 2008) than more traditional approaches to DDMs (Vandekerckhove and Tuerlinckx, 2008; Voss and Voss, 2008) .

To directly link the causal relationship between the sound and its paired greeble image to behavior and drift-diffusion parameters, we included the estimates of causal relationship and transitional probability as predictor variables of the bias, z, of the model. That is, we estimated posterior probability densities not only for basic model parameters, but also the degree to which these parameters are altered by variations in the psychophysical measures (ΔP, causal power, causal support, and transitional probability). In these regressions, the bias parameter is given by z(t)=β0+β1×CP(t), where CP is either ΔP, causal power, causal support, or transitional probability, β₀ is the intercept, and t is the trial number. Here, the slope (β₁) is weighted by the value of the psychophysical measure on that trial. The regression across trials allows us to infer how the bias changes depending on the psychophysical measure. For example, if these psychophysical measures are positively correlated to bias, then increased causal strength or transitional probability will yield faster RTs. We fit four different versions of the model: (1) ΔP model, where bias was estimated from the ΔP and updated after each trial; (2) causal power model, where bias was estimated from the causal power and updated after each trial; (3) causal support model, where bias was estimated from the causal support and updated after each trial; and (4) transitional probability model, where bias was estimated from the transitional probability and updated after each trial. We also modeled our data where drift rate varied according to causal strength or transitional probability. We used the deviance information criterion (DIC) for model comparison (Gelman et al., 2013). The DIC is a measure of model fit (i.e., lack thereof) with a penalty for complexity (i.e., the number of parameters used to fit the model to the data) (Spiegelhalter et al., 2002). Models with lower DIC are better models. Models where bias was manipulated instead of drift rate had significantly lower DIC. For the rest of our analyses, we used models where bias varied with causal strength or transitional probability and drift rate was kept constant at the default values, as this yielded the lowest DIC.

We used Markov Chain Monte Carlo chains with 20,000 samples and 5000 burn-in samples for estimating the posterior distributions of the model parameters. We assessed chain convergence by visually inspecting the autocorrelation distribution, as well as by using the Gelman-Rubin statistic, which compares between-chain and within-chain variance. This statistic was near 1.0 for the parameters, indicating that our sampling was sufficient for proper convergence. We analyzed parameters of the best model (model with lowest DIC) using Bayesian hypothesis testing, where the percentage of samples drawn from the posterior fall within a certain region (e.g., >0). Posterior probabilities ≥95% were considered significant. This value is not equivalent to p values estimated by frequentist methods, but they can be coarsely interpreted in a similar manner.

All the model comparisons were estimated on the psychophysics data as these had the greatest number of trials per condition. This ensures robust estimation of the best model. The best-fitting model was then used to analyze data from different conditions: predrug baseline, under DEX, under ketamine, and recovery from ketamine. To directly link EEG power spectral density to behavior and drift-diffusion parameters, we used the HDDM, but now included the RF cluster power estimate (aligned to image onset) in the alpha band as the predictor variable of the bias, z, of the model; that is, in the regression equation above, CP was now alpha power.

ERPs

We used ERPLAB (https://erpinfo.org/erplab) (Lopez-Calderon and Luck, 2014) to run ERP analysis. First, we cleaned epoched data aligned to the sound onset for auditory ERPs using the pop_artmwppth.m function of ERPLAB with a moving window of 200 ms (2%-3% of trials for each subject were excluded). We then averaged over trials to generate an average ERP for each subject. We used 200-0 ms before stimulus onset as baseline. Based on previous literature (Davies et al., 2010; Winkler et al., 2013), we chose channel 9 (Cz electrode) for exemplar auditory ERPs in figures.

EEG signal decoding

We used two different machine learning approaches in this study: (1) a support vector machine (SVM) model to classify HP, MP, and NP trials after sound onset, with EEG power spectral density as input; and (2) a recurrent neural network (RNN) to classify V1, V2, and V3 trials relative to image onset, with “raw” EEG time series data as input.

We used the power spectral density of EEG signals across four frequency bands, theta (5-7 Hz), α (8-14 Hz), β (15-30 Hz), and γ (31-40 Hz), for the first decoding analysis. We calculated the sum of squared absolute power in each frequency band for each electrode cluster (RF, LC, RC, and OC), thus generating 16 features for each trial as the input dataset. Using Scikit-learn (Pedregosa et al., 2011) implemented in Python, for each 20 ms time bin, we trained an SVM model to classify the EEG data into three classes: HP, MP, and NP. We denote yi(t)∈{0,1,2} as the identifier of the condition at time bin t, where 0, 1, and 2 denote HP, MP, and NP, respectively. The SVM classifier is implemented by a nonlinear projection of the training data x(t) feature space X into a high dimensional feature space F using a kernel function ϕ. So with ϕ:X→F being the mapping kernel, the weight vector w can be expressed as a linear combination of the training trials and the kernel trick can be used to express the discriminant function as follows: y(x(t);ζ(t))=aT(t)ϕx(t)+b(t)=∑n=1Nan(t)φ(xn(t),x(t))+b(t)(5) where ζ(t)={a(t),b(t)} is the new parameter at time bin t with a(t) and b(t) as weights and biases of the mapped features space F. We used the radial basis function kernel that allows nonlinear decision boundary implementation in the input space. The radial basis function kernel holds the elements as follows: φ(xn(t),x(t))=exp(−η(t)‖xn(t)−x(t)‖2)(6) where η(t) is a tunable parameter at each time bin. Model hypermeters consisting of regularization penalty (C(t)) and (η) were selected by grid search through 10-fold cross validation. F score at each time bin and for each label was calculated as follows: F−Score=2 × Precision × RecallPrecision + Recall(7) where Precision=TruePositiveTruePositive + FalsePositive(8) Recall=TruePositiveTruePositive + FalseNegative(9)

As we penalized the mapped weights of the classifier at each time bin, we used normalized absolute values of the weights as a measure to deduce each feature's contribution to classify the outputs.

Second, we used a RNN with many-to-many architecture to decode visual stimuli using the EEG time domain signals from the four clusters of electrodes (RF, RC, LC, OC) mentioned above (20 features in total) in successive 20 ms time bins. The proposed RNN consists of a bidirectional long-short term memory (BiLSTM) layer followed by an attention layer, a fully connected with dropout layer, and a softmax layer as output. Considering x(t) as the input at time t, the output of the BiLSTM forward path is calculated as follows: i(t)=σ(WxiT.x(t) + WhiT.h(t−1) + bi)(10) f(t)=σ(WxfT.x(t) + WhfT.h⇀(t−1) + bf)(11) o(t)=σ(WxoT.x(t) + WhoT.h⇀(t−1) + bo)(12) g(t)=tanh(WxgT.x(t) + WhgT.h⇀(t−1) + bg)(13) c((t))=f((t))⊗c((t−1)) + i((t))⊗g((t))(14) h←((t))=o((t))⊗tanh(c((t)))(15)

Where i(t),f(t),c(t),o(t) are the input gate, forget gate, cell gate, and output gate, respectively, σ(x)=11+e−x and tanh(x)=ex−e−xex+e−x,and the output of the backward path is h←(t). The attention layer output can be calculated as follows: u(t)=tanh(Ws.[h⇀(t);h←(t)] + bs)(16) α(t)=exp(u(t)Tuw)∑texp(u(t)Tuw)(17) s(t)=∑t′≤tα(t′)[h⇀(t′);h←(t′)](18) the output of the attention layer was fed to a fully connected layer with dropout followed by a softmax layer, which results in class conditional probability as Equation 19 specifies. ŷ(t)=softmax(WFC.s(t) + bFC)(19)

All BiLSTM, attention, and FC layer weights and biases will be updated through backpropagation in time. Hyperparameters of the model, including the number of hidden units for BiLSTM, attention and FC layers, learning rate, dropout probability, and learning rate of stochastic gradient descent, etc., are optimized via a grid search. The decoding accuracy is calculated with 10-fold cross validation with 60% of trials as training set and 40% as test set by shuffling and stratifying to avoid accuracy bias because of imbalance classes.

To find out how well the model can be generalized in time, we tested the trained model until time t with the data at all time points using the stored parameters of the above RNN model. The results of this analysis are used to see whether information about the visual stimulus can be decoded before stimulus presentation.

To gain insight on how well the model decodes each visual stimulus during the delay period, we trained the model on data at 100 ms after visual stimulus onset and calculated the F score for each class (V1, V2, V3) under all conditions (before drug, under ketamine, and DEX).

Using the same model (trained on data 100 ms after visual stimulus onset), we looked at the proportion of each auditory stimulus (HP, MP, and NP) associated with the correctly decoded trials at 100 ms before visual stimulus onset.

Statistical analysis

We performed statistical analysis of trials from testing phases and used the learning phase only to confirm that the participants learned the correct associations. To only include the trials where the causal power of HP, MP, and NP trials has differentiated (see Fig. 1G), we excluded the first 50 trials of the testing phase in psychophysics experiments and the predrug baseline. We used all available trials for the testing phase of “under drug” and “after recovery,” as the causal power of HP, MP, and NP trials was already differentiated from the beginning (see Fig. 3A,B). For RT analysis, we excluded RTs more than/less than the mean ± 3 SDs for each subject, and we used correct match trials. For delay period EEG analyses, we used both correct match and nonmatch trials.

Because we had clear, a priori predictions about the effects of the different conditions in this study (i.e., the effect of prediction and the effects of the drug conditions), we could apply contrast analysis using linear mixed effects (LMEs) models. Our study conformed to the guidelines set out by Abelson and Prentice (1997); that is, we included the contrast of interest along with a paired, orthogonal contrast and infer significance only when our statistical tests showed that effects were significant for the contrasts of interest and not for the orthogonal contrasts.

EEG data were analyzed using contrast analysis with an LME where the predictive value of sound (HP, MP, NP; varying within subject) and drug condition (before ketamine, under ketamine, under DEX; varying within subject) served as independent variables. For prediction, we used a linear (−1, 0, 1) contrast as our contrast of interest, and a quadratic (−1, 2, −1) contrast as the orthogonal contrast in the analysis. For drug condition, we used a quadratic (−1, 2, −1) contrast as our contrast of interest, and a linear (−1, 0, 1) contrast as the orthogonal contrast in the analysis. Mathematically, our LME can be written as Equation 20 below: DV∼β0 + β1*PredictionC1 + β2*PredictionC2 + β3*DrugC1 + β4*DrugC2+∑i,j=1 & k=5i,j=2 & k=8βk*PredictionCi*DrugCj(20)

Here, DV is the dependent variable (RF alpha power or average GC at beta frequency), PredictionC₁ is the contrast of interest modeling sounds' predictive value, PredictionC₂ is the orthogonal contrast for sounds' predictive value, DrugC₁ is the contrast of interest modeling drug condition, and DrugC₂ is the orthogonal contrast for drug condition. Terms after the summation notation in Equation 20 correspond to interactions between sounds' predictive value and drug condition (all possible interactions, i.e., permutations of contrasts, included).

We used a similar LME with prediction (HP, MP, NP; varying within subject) and drug condition (under ketamine, before drug, under DEX, after recovery; varying within subject) as independent variables to analyze RT data. We used a linear (HP, MP, NP: −1, 0, 1) contrast as our contrast of interest to model RT, and a quadratic (HP, MP, NP: −1, 2, −1) contrast as the orthogonal contrast in the analysis. For drug condition, we used a four-level (under ketamine, before drug, under DEX, after recovery: 3, −1, −1, −1) contrast as our contrast of interest, and two orthogonal contrasts (under ketamine, before drug, under DEX, after recovery: 0, 0, −1, −1 and 0, 2, −1, −1) in the analysis. Both main effects of sounds' predictive value and drug as well as their interaction were included. We used contrast analysis to model our hypotheses. Mathematically, our LME can be written as Equation 21 below. DV∼β0 + β1*PredictionC1 + β2*PredictionC2 + β3*DrugC1 + β4*DrugC2 + β5*DrugC3 + ∑i,j=1 & k=6i,j=3 & k=10βk*PredictionCi*DrugCj(21)

Here, DV is the dependent variable (RT or accuracy), PredictionC₁ is the contrast of interest modeling sounds' predictive value, PredictionC₂ is the orthogonal contrast for sounds' predictive value, DrugC₁ is the contrast of interest modeling drug condition, and DrugC₂ and DrugC₃ are the orthogonal contrasts for drug condition. Terms after the summation notation of Equation 21 correspond to interactions between sounds' predictive value and drug condition (all permutations of two-way interactions included).

To investigate auditory ERPs, we ran nonparametric permutation tests using the ft_timelockstatistics.m function from FieldTrip (Oostenveld et al., 2011). For before and during ketamine conditions, we looked for differences in ERP between sound cues across 50 to 300 ms latency from sound onset. We randomly shuffled condition labels 10,000 times. Alpha value was set to 0.05 for all comparisons.

We used repeated-measures ANOVA (Holm-Sidak–corrected p values) to test the significance of F scores for the three differentially predictive conditions (HP, MP, NP) as well as the significance of each feature's contribution in output classification.

To test the significance of the cross-temporal decoding results, we performed a paired t test between the 20-fold cross validated and 100 times resampled accuracy at each pixel and randomly permuted output labels model at the same pixel. The resulting p values then were corrected for multiple comparisons (Holm-Sidak corrections).