We invited 30 healthy adults from the NYUAD community to participate in this study, recruited from a larger pool who had previously participated in a separate behavioral study with the same task. In order to ensure that we had sufficient error trials for our analyses, we only invited subjects who had symmetric swaps (see below) on at least 5% of trials in the behavioral study. The final sample included in our analyses was 26 subjects [age: 24.33 ± 4.57, range (19–37); four females, 22 males], after discarding three subjects due to poor performance (fewer than 15 correct responses in one or more conditions of interest), and another due to technical challenges transforming MEG data from sensor space to source space. All subjects were right-handed with normal or corrected-to-normal vision, and provided informed written consent in accordance with procedures approved by NYUAD's IRB.

Stimuli displays, timing, and responses were controlled and recorded using the Psychophysics Toolbox (Brainard, 1997) in MATLAB (The MathWorks, Inc., Natick MA, USA). We assessed WM behavior using a delayed continuous report task (Figure 1A). A fixation cross of 0.2 degrees of visual angle (DoV) was presented in the center of the screen during baseline, stimulus presentation, and memory delays. Each trial began with a 1 s fixation period, after which subjects saw three colored discs (each 0.55 DoV, presented at an eccentricity of 4.5 DoV) for 0.2 s. Subjects were asked to hold the colors and locations of all three discs over a blank 2 s memory delay. At the end of the memory delay, instead of asking subjects to report one of the items, we cued subjects to report the location of each of the three items sequentially (in a random order) by presenting a disc at fixation whose color indicated the color of the item subjects were to report. Subjects adjusted the position of a white circle of the same size and of the same eccentricity as the disc using an MEG-compatible response dial and locked in their report using a button box (Current Designs, Philadelphia PA, USA). Each location report was self-paced. Following all three reports, subjects were shown a 1 s feedback screen that indicated the presented and reported locations on that trial. The trial concluded with a self-paced intertrial interval (ITI)—that is, the next trial began when subjects pressed a button. Subjects completed 500 ± 75 (400–720) trials, divided into blocks of 80 trials (seven subjects) or 100 trials (19 subjects). Two subjects completed the experiment in two sessions 1 and 4 days apart; the remaining subjects completed the experiment in a single session.

The memory items were always presented in one hemifield in a given block of trials, and left and right hemifield blocks were presented in a random order. The color of the memory items on each trial were chosen from 180 color equally spaced segments drawn from a circle in CIELAB color space (radius 59°, centered at L = 54, a = 18, and b = −8) (CIE, 2004), with a minimum gap of 15 color segments between memory items on a given trial to avoid errors due to the color similarity of the memory items. Memory locations were randomly sampled from an isoeccentric (semi)circular space (depending on the hemifield block) with a minimum gap of 15° of polar angle. We additionally excluded locations that were 10° of polar angle away from the vertical meridian to avoid well-documented oblique effects (Furmanski and Engel, 2000).

Stimulus displays were projected onto a screen 85 cm away from the subject's head. Subjects were instructed to remain still during the MEG recordings and to avoid blinking or moving their eyes during the cue and delay screens. Subjects were free to move their eyes during the response screens and were encouraged to use the ITI as well as breaks between blocks to blink.

To label trials based on subjects' behavioral performance, we fit a probabilistic mixture model with three components (target, non-target, and guess responses) to the distribution of errors using the C016_fit function (currently called mixtureFit, bayslab.com) (Bays et al., 2009; Schneegans and Bays, 2016):

where θ is the target location (in radians), θ^ is the reported location, γ is the probability of random guess, β is the probability of swap, θi* are the locations of non-targets with i = 1, …, m and m = 2, and ϕ_σ is the circular normal distribution (Von Mises) with zero mean and σ standard deviation. Maximum likelihood estimates of parameters σ, γ, and β were obtained at subject level. Then, trial-wise posterior probabilities that responses were drawn from each of the three mixture components were used to define trial types (see Figure 1C in Bays et al., 2009 and Figure 2 in Schneegans and Bays, 2016). High accuracy (correct) responses were trials with a target probability above 90% in all 3 reports, swaps errors were trials with a probability of non-target response above 70% in at least one report, and low accuracy trials comprised the remaining trials. The more lenient cut-off criterion for swap trials was motivated by the tendency of the model to underestimate swap frequency, while the conservative cut-off in correct trials in all three reports minimized the likelihood that correct trials contain lucky random guesses. Note that the label “low accuracy” does not mean that all three reports were inaccurate or guess responses; it simply indicates that subjects did not report all three items with a high degree of accuracy.

Before MEG acquisition, each subject's head shape was digitized using a Polhemus dual source handheld FastSCAN-II. MEG data was recorded continuously using a 208-channel axial gradiometer Yokogawa system (Kanazawa Institute of Technology, Kanazawa, Japan) with a sampling rate of 1,000 Hz and an online low-pass filter of 200 Hz. The continuous MEG data was first noise-reduced using eight magnetometer reference channels located away from the participant's head and using the Time Shifted Principle Component Analysis (TSPCA, block width of 5,000 and 30 shifts) as implemented in MEG160 software (Yokogawa Electric Corporation and Eagle Technology Corporation, Tokyo, Japan). To remove eye blink and heartbeat artifacts, an independent component analysis was performed with fast ICA method in MNE-Python (Ablin et al., 2018). All subsequent preprocessing steps were performed using FieldTrip (Oostenveld et al., 2011). Data was epoched from 1 s pre-stimulus onset to 2.7 s post-stimulus onset (including baseline, stimulus presentation, memory delay, and the first 0.5 s of response), demeaned, and low pass filtered at 140 Hz. In order to clean the data, we automatically rejected trials and channels with high variance using ft_artifact_zvalue.m function with the following settings: jumps artifacts threshold z = 20, muscle artifacts threshold z = 5 and blink artifacts threshold z = 10. Additionally, we manually inspected individual epochs to remove trials that were contaminated with jump, muscle, and eye blink artifacts. Rejected channels were reconstructed using the average of all neighbors. Preprocessing resulted in an average of 54 ± 37 (13–200) removed trials [accounting for the 11 ± 6% of total trials (2–29)], and 2 ± 1 (1–10) rejected channels (out of 208) per participant. The average number of trials available for analysis was 446 ± 69 (357–707): 106 ± 94 (19–439) (23 ± 17%) high accuracy trials, 213 ± 53 (80–328) (48 ± 11%) low accuracy trials, and 126 ± 46 (46–250) (29 ± 10%) swap trials. The experimental setup included trials with stimuli presented in either the left or right hemifield. To analyze all trials without canceling out lateralized neural signals, we took the following approach. For trials with left hemifield presentation, we mirrored (reflected) the MEG sensors across the midline, ensuring that MEG sensors on the left (right) side aligned with the hemisphere contralateral (ipsilateral) to the stimulus presentation (as depicted in Figure 2).

In a subset of eight subjects [age: 23.16 ± 3.54, range (20–31); two females, six males], we used a high-resolution anatomical MRI volume collected during a separate session to constrain the transformation from sensor to source space. The anatomical scan was acquired on a Siemens Magnetom Prisma 3T MRI scanner using a 6.5 min MPRAGE-3D T1-weighted, gradient-echo acquisition sequence (TR: 2,400 ms, TE = 2.22 ms, flip angle: 8 degrees, voxel size: 0.8 mm³, 208 slices, FOV: 256 mm).

In order to study effective connectivity, we performed a time-frequency analysis in source space to identify the network of neural sources actively involved in WM, comprising frontal, parietal, and occipital regions. First, we localized oscillatory sources using the beamformer technique (Jaiswal et al., 2020; Westner et al., 2022). The head model template in MNI space (“standard_bem”) (Oostenveld et al., 2003) was transformed into a subject-specific head shape to compute the volume conduction model using a realistically shaped single shell approximation from FieldTrip's ft_prepare_headmodel.m function (Nolte, 2003). Next, we used ft_prepare_sourcemodel.m function to create a source model with a 1.25 cm isotropic 3D grid and with 1,499 sources within the head model. Finally, source and head models were used to create a forward model (ft_prepare_leadfield.m). For the subset of eight subjects with MRI data, the head model was computed using the structural T1 image, and the source model in standard space was spatially normalized into subject-specific structural image with a non-linear transformation (warp). As a result of this transformation, forward model positions were equivalent across subjects in normalized space. To compute the inverse solution, the cross-spectral density (CSD) matrix was calculated using discrete prolate spheroidal sequences (“dpss” taper in “mtmfft” method in ft_freqanalysis.m function) with all valid trials using baseline (−0.75 to −0.25 s) and delay period (0.2–2.2 s) epochs. The resulting spatial filter was used to localize sources in the frequency domain (computing power) using a beamformer dipole analysis with dynamic imaging of coherent sources (“dics” method in ft_sourceanalysis.m function) for frequencies from 2 to 100 Hz in steps of 2 Hz, with a smoothing of ± 2Hz, and a 5% regularization. For each trial type, source activity (power) during the delay period was computed in 31 time windows of 0.5 s, in steps of 0.05 s, and baseline-corrected (relative change method) with the trial-averaged baseline epoch of 0.5 s (−0.75 to −0.25 s, equally spaced from the beginning and end of the 1 s baseline). Finally, to normalize spatial maps across subjects, we transformed subject-specific maps of power from the head coordinate system to the MNI coordinate system using the inverse of the transformation matrix that was computed when generating the forward model.

Sources of interest for the DCM analysis were selected using a combined theory- and data-driven approach. Based on previous literature (Sreenivasan and D'Esposito, 2019), we predefined regions-of-interest (ROI) using the Brainnetome Atlas (Fan et al., 2016): superior frontal gyrus (SFG, areas 1–14), middle frontal gyrus (MFG, 15–26), intraparietal sulcus (IPS, 127–142), and early visual cortex (EVC, 203–206). Then, we used the time frequency matrices of power (31 time windows by 50 frequency bins) to select the sources of interest within each ROI. For each source, we performed a massive univariate one-way analysis of variance (ANOVA) with trial type as a factor with three levels (high accuracy, low accuracy, and swap trials), getting as output a time frequency matrix of F-statistics. Then, for each source, we summed the resulting uncorrected F-statistics over time and frequency, and selected the source with highest collapsed F-statistic within each ROI. All ROIs were considered separately for contralateral and ipsilateral hemispheres, except in SFG. There, given the proximity of the sources to the midline, we selected the source with the highest F-statistic bilaterally. Sources coordinates were (Figure 1B): SFG [F_collapsed = 3,576; MNI: (8, 18, 68)], MFG [contralateral: F_collapsed = 3,399, (−43, 30, 30); ipsilateral: F_collapsed = 3,760, (45, 30, 30)], IPS [contralateral: F_collapsed = 5,487, (−43, −58, 43); ipsilateral: F_collapsed = 4,129, (58, −33, 55)] and EVC [contralateral: F_collapsed = 3,108, (−18, −95, 5); ipsilateral: F_collapsed = 3,161, (20, −95, −8)]. Frontal and parietal nodes fall inside of a recent meta-analysis of fMRI studies about delay activity in WM (Li et al., 2022), while visual nodes fall in the occipital pole, from where previous studies decoded memory contents during delay period (Serences et al., 2009). Timecourses of selected sources of interest were reconstructed using a linear constrained minimum variance beamformer (“lcmv” method in FieldTrip) (Veen et al., 1997) with a fixed orientation projected along the first PCA axis.

We used a DCM with cross-spectral density as data feature (range 2–100 Hz in steps of 1 Hz), as implemented in SPM12 (Friston et al., 2012; Moran et al., 2013), with the canonical microcircuit model (CMM) based on the 4-population conductance-based neural mass model (Moran et al., 2011; Bastos et al., 2012). The model comprises two components: the neuronal model characterizes the dynamics of hidden neural states,

and the observation model describes how neural activity gives rise to observed data,

where x is the hidden neural state, u is the inputs to the neuronal population, θ are the model parameters, Γ_x are the stochastic noise terms in neuronal model, y is the observed MEG data, h is the function mapping neural states to observed data, involving a lead matrix C, and ϵ are stochastic noise terms in the observation model.

In the conductance-based model used in the current study, neural dynamics are parametrized as the summed active and passive currents across the membrane (Pereira et al., 2021). Each neuronal population was modeled based on the Morris-Lecar model,

where C is the membrane capacitance, V is the membrane potential, g_k is the conductance of channel k, V_K is the reversal potential for channel k, u is the applied input current, Γ_V is the stochastic term modeling Gaussian noise. k comprises passive leak current, and active currents of excitatory (Na⁺) and inhibitory (Cl⁻) ion flow mediated by fast AMPA and GABAA receptors, depending on the type of neuronal population. Hence, membrane potential V and conductance of channels g_k comprise the hidden states x of the model in each neuronal population. Please refer to Pereira et al. (2021) for a full description of the mathematics of conductance based models, and see Moran et al. (2013) for a broad overview of neural mass models. Within each source (brain region), four neuronal populations are modeled: spiny stellate cells, superficial pyramidal cells, inhibitory interneurons, and deep pyramidal cells (Bastos et al., 2012). Neuronal populations within each source were connected as follows: (i) excitatory connections between excitatory spiny stellate and superficial pyramidal cells and from superficial pyramidal to deep pyramidal cells, (ii) recurrent connections between excitatory populations and the interlaminar population of inhibitory interneurons, and (iii) inhibitory self-connections of each population, modeling synaptic gain (Friston et al., 2019).

The connectivity between sources, the seven brain regions described above, was modeled through the matrices A = [A{1}, A{2}], the intrinsic connectivity matrices representing the strength of feedforward A{1} and feedback A{2} connections between different neural sources (see model architecture below). Types of connections were defined based on the specific populations they projected to. Feedforward connections originated from superficial pyramidal in the projecting source to spiny stellate in the target source, feedback connections originated from deep pyramidal in the projecting source to both inhibitory interneurons and superficial pyramidal cells in the target source, and lateral (interhemispheric) connection comprised both feedforward and feedback connections. Since our focus of interest was not the comparison between different model architectures of extrinsic connectivity (inference on model architecture), but rather the performance-dependent modulation of effective connectivity (inference on model parameters), we defined a model architecture with the canonical functional hierarchy. Building upon research on the functional organization of brain networks (Bastos et al., 2015; Michalareas et al., 2016; Nee and D'Esposito, 2016; Gratton C. et al., 2018; Gratton G. et al., 2018; Marek and Dosenbach, 2018) and previous DCM studies about WM (Dima et al., 2014; Heinzel et al., 2017), prefrontal cortex was placed at the top and visual cortex at the bottom of the hierarchy, designating posterior-to-anterior connections as feedforward, anterior-to-posterior connections as feedback, and interhemispheric connections between regions at the same level of the hierarchy as lateral connections (Figure 1C). Finally, the experimental effects of interests, the baseline corrected estimates of effective connectivity during delay period, were modeled with the matrix B, comprising all the connections defined in matrices A: 12 feedforward, 12 feedback, and six lateral connections (Figure 1C). Note that there is no circularity in using the same data to define ROIs and to study their interactions with DCM, since the purpose of DCM is to test hypothesis about the mechanisms underlying the experimental effects observed with conventional analysis (Stephan et al., 2010). Massive univariate statistical analyses (e.g., statistical parametric maps) to define ROIS are standard practice in several DCM studies (Moran et al., 2011; Dima et al., 2014; Auksztulewicz and Friston, 2015; Heinzel et al., 2017; Jung et al., 2018; Friston et al., 2019; Adams et al., 2021).

For each subject, we inverted two DCMs for each trial type (high accuracy trials, low accuracy trials, and swap trials), encoding the baseline corrected estimates of effective connectivity (matrix B) during early and late delay. All time intervals had a length of 800 ms [baseline (−900 to −100 ms), early delay (400–1,200 ms), and late delay (1,200–2,000 ms)], and were separated from stimulus presentation (0 to 200 ms) or report (from 2,200 ms) by at least 100 ms. The separation of delay period in early and late time windows was motivated by previous studies in non-human primates showing that the dynamic properties of WM coding were stronger during the initial portions of the memory delay (Murray et al., 2017), and the modulation of memory load showed distinct spectral profiles during the early and late portions of the memory delay (Buschman et al., 2011) as well as differences in effective connectivity in frontoparietal networks (Pinotsis et al., 2019). Models were inverted to fit the cross-spectral density through the optimization of variational free energy under Laplace approximation (Friston et al., 2007) to optimize log scaling parameters around default priors (spm_fx_cmm.m). Please see Zeidman et al. (2023) for a comprehensive technical review of the Bayesian inference scheme employed to fit DCM models.

We performed statistical inference on the parameters of the optimized DCM model using PEB (Friston et al., 2015; Zeidman et al., 2019). Compared to standard statistical inference, PEB has the advantage of taking into account the estimated uncertainty (variance) about model parameters, downweighing the contribution of less certain parameter estimates on group effects. The PEB approach is also advantageous in the context of model inversion and comparison, since the inversion of the same “full” DCM per subject/condition avoids the issue of fitting multiple models and performing Bayesian model comparison over models, which can lead to different DCMs falling into distinct local optima. PEB is a Bayesian hierarchical model over parameters with individual DCMs as a first level, and group effects on DCM parameters as a second level (see Zeidman et al., 2019 for details).

The first level (Equation 1) represents the observed MEG data Y of the i-th subject modeled with the dynamic causal model Γ with parameters θi(1) and noise term εi(1). The second level (Equation 2) shows that parameters of interest from first level θ⁽¹⁾ are modeled at group level through a GLM with design matrix X and group-level parameters θi(2), with zero-mean noise to account for between-subject variability as random effects εi(2). In our study, parameters of interest in PEB analysis (within-subjects design matrix) were the baseline corrected estimates of extrinsic (between-source) effective connectivity (matrix B). We performed pairwise comparisons between trial types with a group-level design matrix with an intercept (column of 1's) encoding commonalities, and a regressor of interest encoding the pairwise contrast between trial types (e.g., high accuracy vs. low accuracy). We report PEB parameters with strong evidence (posterior probability, Pp > 0.95) (Kass and Raftery, 1995), and a large effect size, assessed through its predictive accuracy in a leave-one-out cross-validation (spm_dcm_loo.m). In this procedure, a PEB model was constructed with all subjects but one, predicting the excluded subject's covariates, and iterated through each participant serving as the left-out subject. The predictive accuracy was quantified as Pearson's correlation between PEB covariates and cross-validation predictions. We applied a statistical threshold of p = 0.05, uncorrected.

All scripts to invert DCMs and perform PEB analysis are available online: https://github.com/asantoangles/dcm_working_memory. Source timecourses used as input for DCM, in addition to inverted DCMs and PEB results, are available on Open Science Framework: https://osf.io/hbgu2/?view_only=8bc743e2068d47f0889f7b556beb2068.

Absolute mean error was 24.4° ± 5.17° (14.2°, 34.5°) for all trials, indicating that subjects were successfully able to complete the task. As expected by the definition of trial types, performance in high accuracy trials [11.2° ± 1.9° (7.4°, 16°)] was better than low accuracy trials [21.5° ± 4.9° [(14.3°, 33.8°); two-sample Wilcoxon test: W = 354, z = −6.12, p < 0.0001, effect size r = 0.85], and swap trials [51.6° ± 5.2° (43.4°, 68°); W = 351, z = −6.17, p < 0.0001, r = 0.85]. Performance in low accuracy trials was significantly better than swap trials (W = 351, z = −6.17, p < 0.0001, r = 0.85). Additionally, we observed no significant differences in absolute mean error between trials with left vs. right hemifield stimulus presentation [W = 734, p = 0.8; left: 25.2° ± 5.7° (13.2°, 36.5°), right: 23.7° ± 4.9° (13.9°, 32.5°)] and no difference in the proportion of left/right trials observed in distinct trial types (high accuracy trials: 0.49/0.51, low accuracy trials: 0.5/0.5, swap trials: 0.51/0.49; all pairwise comparisons between trials with a Chi-square test with Yates continuity yielded p-values > 0.8). The similar performance across hemifields helped justify our decision to flip sensors across the midline for MEG analysis (see Methods).

Correlation between empirical and predicted amplitude (modulus) was excellent (r = 0.9907 ± 0.004). Model predictions of within-ROI power spectral density (PSD) matched empirical data in all frequencies. Between-ROIs cross-spectral density (CSD) predictions were below the empirical values, although empirical CSD differences between trial types were generally reproduced by model predictions (Supplementary Figure S1B). Coefficient of determination also showed a good model fit [R² = 0.97 ± 0.01 (0.46–0.99)]. PEB results summarized below describe the effects with strong evidence (Pp > 0.95) an out-of-sample predictive accuracy with p < 0.05, uncorrected (see Supplementary material for a detailed description of PEB results).

Low accuracy trials, in comparison with high accuracy trials (Figure 2A), during early delay, showed reduced feedback signals in ipsilateral frontooccipital and medial to superior frontal connections, and increased interhemispheric connections between parietal and frontal areas. During late delay, low accuracy trials showed reduced feedforward signals in ipsilateral frontoparietal and frontooccipital connections, and reduced lateral connections between early visual areas. Swap trials, in comparison with high accuracy trials (Figure 2B), during early delay, showed reduced feedback connectivity in ipsilateral frontooccipital and increased superior frontal to ipsilateral parietal. Feedback connectivity was also reduced in superior frontal to contralateral early visual cortex during the whole delay period. Low accuracy trials, in comparison with swap trials (Figure 2C), during early delay, showed reduced feedback ipsilateral frontoparietal connectivity and reduced feedforward connectivity in ipsilateral early visual cortex to superior frontal. During late delay, low accuracy trials showed increased feedback connectivity in several frontooccipital connections: ipsilateral middle frontal to early visual cortex, and superior frontal to ipsilateral and contralateral visual areas.