Computational Modeling of Electroencephalography and Functional Magnetic Resonance Imaging Paradigms Indicates a Consistent Loss of Pyramidal Cell Synaptic Gain in Schizophrenia

Background Diminished synaptic gain—the sensitivity of postsynaptic responses to neural inputs—may be a fundamental synaptic pathology in schizophrenia. Evidence for this is indirect, however. Furthermore, it is unclear whether pyramidal cells or interneurons (or both) are affected, or how these deficits relate to symptoms. Methods People with schizophrenia diagnoses (PScz) (n = 108), their relatives (n = 57), and control subjects (n = 107) underwent 3 electroencephalography (EEG) paradigms—resting, mismatch negativity, and 40-Hz auditory steady-state response—and resting functional magnetic resonance imaging. Dynamic causal modeling was used to quantify synaptic connectivity in cortical microcircuits. Results Classic group differences in EEG features between PScz and control subjects were replicated, including increased theta and other spectral changes (resting EEG), reduced mismatch negativity, and reduced 40-Hz power. Across all 4 paradigms, characteristic PScz data features were all best explained by models with greater self-inhibition (decreased synaptic gain) in pyramidal cells. Furthermore, disinhibition in auditory areas predicted abnormal auditory perception (and positive symptoms) in PScz in 3 paradigms. Conclusions First, characteristic EEG changes in PScz in 3 classic paradigms are all attributable to the same underlying parameter change: greater self-inhibition in pyramidal cells. Second, psychotic symptoms in PScz relate to disinhibition in neural circuits. These findings are more commensurate with the hypothesis that in PScz, a primary loss of synaptic gain on pyramidal cells is then compensated by interneuron downregulation (rather than the converse). They further suggest that psychotic symptoms relate to this secondary downregulation.

paradigms, but these data were not analysed in this study. We focused on the measures of auditory perceptions and processing speed as these have the largest effect sizes in PScz (4)).
Use of psychotropic medication was recorded. 90 PScz were taking one or more antipsychotics (including Clozapine; see Table S1), for 15 PScz no medication was recorded and 3 were unmedicated. Some participants were also taking antidepressants, hypnotics (as prescribed) or mood stabilisers -the PScz group had higher proportions in each case. The 20 item Brief Psychiatric Rating Scale (BPRS) was used to rate overall symptoms in the PScz group only.

Data acquisition and paradigms
For EEG recording, participants sat in a semi-reclining chair inside a sound-attenuated chamber, wearing an electrode cap. 64 scalp electrode sites were recorded, according to the 10/20 International System, using Neuroscan Stim acquisition software and a Synamps2 or Synamps2 RT amplifier. Recordings were grounded midway between FPZ and FZ using silver/silver-chloride electrodes and referenced to the nose. Eye movements were monitored by vertical and horizontal electrooculograms (EOGs). Data were continuously sampled at 1000 Hz with a DC/100 Hz band-pass filter (24 dB/octave roll-off). Impedances were kept below 5 kΩ.
For the resting state EEG (rsEEG) paradigm, subjects were asked to remain awake whilst EEG data were recorded. Two recordings were obtained on the same day, one with eyes open, and one with eyes closed. Each recording lasted 5 minutes.
In the mismatch negativity (MMN) paradigm, participants were presented with 1000 auditory stimuli through earphones, of which 800 (80%) were standard tones presented at 75 dB, for 60 msec, at 1000 Hz; 200 (20%) were duration-deviant tones at 75 dB, for 150 msec, at 1000 Hz. All tones had a 5 msec rise/fall time, with a stimulus onset asynchrony of 300 msec. 800 tones were presented in a single block. Participants were asked to ignore the tones while viewing a silent movie. Some of the MMN data have been analysed elsewhere (5), using classical statistical tests and structural equation modeling to find correlations between MMN responses and neurotransmitter levels. Here we obtained effective connectivity estimates using DCM.
In the auditory steady state response (ASSR) paradigm, participants listened to click trains delivered by headphones, at 2.5, 5, 10, 20, 40, and 80 Hz. 75 stimulus trains (trials) -each consisting of 15 clicks, with each click at 72 dB and of 1 ms duration -were delivered at each stimulus frequency. The duration ranged from 6 s per train for 2.5 Hz, to 0.1875 s per train for 80 Hz. The inter-train interval was 0.7 s. Therefore, the durations for 2.5, 5, 10, 20, 40, and 80 Hz were 8.38, 4.69, 2.82, 1.89, 1.42, and 1.19 min, respectively, presented in 6 separate blocks separated by 2 min. The order of the blocks was randomized. In this study, only the data for the 40 Hz ASSR were analysed: data from this and the other frequency bands have been analysed elsewhere (3).
For the resting state fMRI (rsfMRI) paradigm, participants were asked to keep their eyes closed and relax. fMRI data were collected at the University of Maryland Center for Brain Imaging Research using a Siemens 3T TRIO MRI (Erlangen, Germany) system equipped with a 32-channel phase array head coil. Resting state functional Blood Oxygen Level Dependent (BOLD) T2*-weighted images were obtained with axial slices parallel to the anterior-posterior commissure (AC-PC) using a T2*-weighted gradient-echo, echo-planar sequence (TR=2.2 s, TE=2.8 ms, flip angle=13°, FOV=256 mm, 128×128 matrix, 1.94x1.94 in-plane resolution, 4mm slice thickness, 37 axial slices, 444 volumes per run; 1 run).

Data preprocessing
All EEG preprocessing and analysis was performed within Matlab (The MathWorks Inc.).
MMN and 40 Hz ASSR data were imported into EEGLAB (6) (https://sccn.ucsd.edu/eeglab/) and band-pass filtered between 1 and 70 Hz, with a notch filter at 59.5-60.5 Hz. rsEEG data were band-pass filtered between 1 and 50 Hz as higher γ frequencies are more vulnerable to artefacts in non-timelocked data. The data were then epoched: rsEEG into 6 s epochs, MMN into epochs from -50 ms to 300 ms around the stimulus onset, and 40 Hz ASSR data into epochs from -100 ms to 600 ms around the stimulus onset. The data then underwent an automated artefact rejection pipeline within EEGLAB, also employing the Multiple Artefact Rejection Algorithm (MARA) (7). First, epochs with amplitudes of >±5 std were rejected, as were epochs containing linear trends (using pop_rejtrend, max slope=5, min R2=0.7) and low (0-2 Hz) or high (20-40 Hz) frequency bands (using pop_rejspec) within power thresholds of -50 to 50 dB and -100 to 25 dB respectively (40 Hz ASSR data used a frequency band of 20-35 Hz given the click train was presented at 40 Hz). Channels were then rejected (using pop_rejchan) on the basis of extreme values in their spectrum (-4 to 6 std), kurtosis (-7 to 15 std) or joint probability (-9 to 7 std). If >50% epochs were rejected, channel rejection was performed first, and epoch rejection second. If >50% epochs were still rejected, or if >20% channels were rejected, the data was discarded. For the MMN, 17/246 subjects ultimately had >50% epochs rejected but no subjects had >20% channels rejected. For the 40 Hz ASSR, 12/259 subjects ultimately had >50% epochs rejected but no subjects had >20% channels rejected.
Independent component analysis (ICA) was then performed, and components with p>0.5 of being artefacts -determined using stored artefact templates within MARA -were then rejected. MARA contains a supervised machine learning algorithm that learns from expert ratings of 1290 components by extracting features from the spectral, spatial and temporal domains, in order to identify artefacts of all kinds (e.g. eye, muscle, loose electrodes, etc).
After MARA, channel and epoch rejection were run a second time using lower thresholds, as sometimes artefacts became apparent following component rejection. Channels were rejected on the basis of the following values in their spectrum (-6 to 5 std), kurtosis (-6 to 9 std) or joint probability (-7 to 7 std). Missing channels were then interpolated and the data were rereferenced to the common average. The data were then converted to SPM format and imported into SPM12 (v7219) (8), available at https://www.fil.ion.ucl.ac.uk/spm/software/spm12/. The MMN and 40 Hz ASSR data were baseline corrected to the first 50 ms of the epoch, and the MMN data were downsampled to 250 Hz. The MMN (standard and deviant conditions) and 40 Hz ASSR were then averaged within each subject using robust averaging (using default options for spm_robust_average, which preserves roughly 95% of data points).
The structural MRI and rsfMRI data were preprocessed in accordance with the Human Connectome Project (HCP) Minimal Preprocessing Pipelines (MPP), adapted to be compatible with "legacy" data (single-band BOLD, no field map or T2-weighted image).
These modifications to the MPP are described in detail in prior publications (9) and are offered as a standard HCP MPP option at https://github.com/Washington-University/HCPpipelines/pull/156. Briefly, the T1-weighted images were first aligned in the volume space to the standard Montreal Neurological Institute-152 (MNI-152) brain template via a single-step combined transformation using the FMRIB Software Library (FSL) linear image registration tool (FLIRT) and non-linear image registration tool (FNIRT) (10).
FreeSurfer (http://surfer.nmr.mgh.harvard.edu/) was then used to segment gray and white matter in the whole brain to produce individual cortical-subcortical anatomical segmentations (11) -specifically the pial and white matter cortical boundaries which were used to define a "cortical ribbon", and a subcortical grey matter voxel mask. These were then used to define the hybrid surface/volume neural file for each subject in the Connectivity Informatics Technology Initiative (CIFTI) "grayordinate" space (12). All images for all subjects were visually inspected to assess the T1-weighted image quality and FreeSurfer grey/white matter segmentation performance. The T1-weighted images were checked for ringing, blurring due to motion, signal inhomogeneity and other artefacts or anomalies. If severe enough to impair segmentation, the subject was omitted. Otherwise, if segmentation failed in a cortical area in >2 consecutive slices, the segmentation was repeated with adjusted parameters and reassessed. If segmentation failure persisted, the subject was omitted. This rigorous procedure excluded 67/269 subjects (Table S1).
BOLD images were first slice time corrected using FSL slicetimer. Next, motion correction was performed by aligning all BOLD data to the first frame of every run via McFLIRT.
Motion corrected BOLD images were then registered to T1w image using FreeSurfer BBR.
Finally, all the linear transformation matrices and nonlinear warp images were combined and the original slice time corrected BOLD images were registered to MNI-152 brain template in a single step. To exclude signal from non-brain tissue a brainmask was applied and cortical BOLD data were converted to the CIFTI gray matter matrix by sampling from the anatomically-defined cortical ribbon. These cortical data were then aligned to the HCP atlas using surface-based nonlinear deformation (12). Subcortical voxels in the BOLD data were extracted using the Freesurfer-defined segmentation to isolate the subcortical volume portion of the CIFTI space.
In addition to the HCP MPP, all BOLD images had to pass stringent quality assurance criteria as previously reported (13) to ensure that all functional data were of comparable and high quality: i) signal-to-noise ratios (SNR) >90, computed by obtaining the mean signal and standard deviation (sd) for a given slice across the BOLD run, while excluding all non-brain voxels across all frames (14); ii) movement scrubbing as recommended by Power et al (15,16). 'Movement scrubbing' refers to the practice of removing BOLD volumes that have been flagged for high motion, in order to minimize movement artefacts, and is a widely used fMRI preprocessing technique. Specifically, to further remove head motion artefacts, as accomplished previously (13), all image frames with possible movement-induced artefactual fluctuations in intensity were identified via two criteria: First, frames in which sum of the displacement across all 6 rigid body movement correction parameters exceeded 0.5 mm (assuming 50 mm cortical sphere radius) were identified; second, root mean square (RMS) of differences in intensity between the current and preceding frame was computed across all voxels divided by mean intensity and normalized to time series median. Frames in which normalized RMS exceeded 1.6 times the median across scans were identified. The frames flagged by either criterion were marked for exclusion (logical or), as well as the one preceding and two frames following the flagged frame. Movement scrubbing was performed for all reported analyses across all subjects. Subjects with more than 50% of frames flagged were removed from subsequent analyses. Next, to remove spurious signal in resting state data we completed additional preprocessing steps, as is standard practice (17): all BOLD timeseries underwent high pass temporal filtering (>0.008 Hz), removal of nuisance signal extracted from anatomically-defined ventricles, white matter, and the remaining brain voxels (i.e. global signal) (all identified via participant-specific FreeSurfer segmentations (18)), as well as 6 rigid-body motion correction parameters, and their first derivatives using previously validated in-house Matlab tools (19). Note that while the removal of global signal is a debated topic in fMRI (20,21), it remains the gold standard for removing spatially pervasive artefact in the brain (although other techniques are emerging (22,23)). BOLD signal within the subject-specific cortical mask was spatially smoothed with a 6 mm full-width-at-halfmaximum (FWHM) Gaussian kernel and dilated by two voxels (6 mm) to account for individual differences in anatomy, as done previously (14,24). The cortical data were then parcellated using the 360-area multi-modal parcellation from Glasser et al. (2016) (25).

Data analysis and modelling: resting state EEG
The power spectrum for each electrode in each subject was computed using the default settings of Welch's method in Matlab (6 s windows with overlaps of 3 s). Each power spectrum was normalised (sensor by sensor) by subtracting the 1/f slope in log space: this gradient was determined using the default settings (bisquare weight function) of robustfit in Matlab. The normalised spectra were averaged over all sensors and all subjects in each group.
The biophysical model used to simulate power spectra is a convolution-based neural mass model (26) in DCM. It is a neural mass model in the sense that it treats ensemble neural activity as a point process (without spatial extension) and only first order statistics (means) of neural population activities influence each other. Convolution-based models (as opposed to conductance-based models, which model ion channels in cell membranes) convolve presynaptic inputs (firing rates) with synaptic response kernels to produce postsynaptic membrane potentials in a given neural population (27,28). The parameters summarising the kernel are the maximum postsynaptic depolarisation (determined by receptor density, or 'connectivity') and membrane time constants. These parameters are listed as connectivity parameters G and time constants T in Figure S1A, and the equations of the model given in Figure S1B. The model transforms the neural population's membrane potential back into a firing rate using a sigmoid operator whose slope is controlled by parameter S. This firing rate becomes the input for another population, determined by the network's intrinsic and extrinsic connectivity. Delays in transmission within the microcircuit (intrinsic) and from one brain area to another (extrinsic) are parameterised by D (not shown). The microcircuit itself ( Figure   S1A) contains four neural populations: spiny stellate cells (ss), superficial pyramidal cells (sp), deep pyramidal cells (dp) and inhibitory interneurons (ii). Extrinsic connections follow known anatomical patterns (29): 'forward' connections project from sp cells to ss and dp cells in the area above, and 'backward' connections project from dp cells to sp and ii cells in the area below. Given all of the above, the dynamic activity of all populations in the network can be computed for any given input (for evoked responses like the MMN) or as a filtered spectral response to endogenous neuronal fluctuations (a mixture of pink and white neuronal noise) as in the rsEEG (26).
There are several forms of 'gain' within a convolutional neural mass model (illustrated in Figure S1B). One is 'somatic gain', the slope S of the sigmoid function converting membrane potential to firing rate. Another is 'synaptic gain', determined by the height (due to connectivity, G) and spread of the synaptic kernel . To assess the synaptic gain of sp and other cells in this work, we use the self-inhibitory G parameters (e.g. sp-sp) for that population, because these apply irrespective of the afferent inputs (unlike other G parameters, e.g. ss-sp, which pertain to specific connections). Note, however, that self-inhibitory connections may parameterise one or more of several physiological effects: i) they control the responsiveness of a population to its inputs, as any mechanism controlling synaptic gain does, e.g. NMDAR function, but also classical neuromodulatory receptors such as cholinergic and dopaminergic receptors, or the connectivity between neurons within a pyramidal population; ii) they may reflect the action of an inhibitory interneuron population in a circuit, e.g. from sp cells projecting to parvalbumin+ (PV + ) interneurons and back to sp cells (30); iii) in the case of interneurons, they may reflect autapses (self-synapses: common on PV + cells (31)). A model cannot distinguish between these mechanisms unless they are explicitly encoded in the model itself; however, some interpretations are more or less likely given what we know of pathology in schizophrenia. For example, a finding of increased sp self-inhibition in PScz is more likely to be due to loss of synaptic gain in those cells than a strengthening of the sp-PV + -sp circuit, which most evidence implies is weakened in PScz (32). The opposite is true for a finding of decreased sp self-inhibition in PScz: its most likely explanation would be a loss of ii inhibition of sp cells.
Similarly, the microcircuit model does not explicitly distinguish between PV + and somatostatin+ (SST + ) interneurons. SST + interneurons are of interest as their cortical markers are just as (if not more) reduced in PScz as PV + (33). That said, the dynamics of the sp selfinhibition connection in the model are faster than those of the sp-ii-sp circuit (see the parameters in Figure S1A), hence the model can potentially distinguish between faster (likely PV + ) and slower (e.g. SST + ) interactions between pyramidal cells and interneurons (30).
(Note however that the empirical priors used for the 40 Hz ASSR analysis accelerated the dynamics of the sp-ii-sp circuit, in order to model the 40 Hz peak.) The microcircuit model used was the same as that used by Shaw et al (34): we termed this spm_fx_cmc_2017. It makes a small adjustment to spm_fx_cmc_2014, in that it replaces a G connection from sp to ss with a connection from sp to dp: see Figure 2C and Figure S1A. We used this model because it is closer to known anatomy (29) than previous versions -e.g.
having distinct but connected superficial and deep pyramidal populations that send ascending and descending projections respectively -but not too complicated to fit to EEG data (e.g. more sophisticated models (35) contain separate interneuron pools for sp and dp cells, and are conductance-based, i.e. contain parameters for specific receptors, but they are harder to fit to EEG (36)). Given our primary interest was in the function of (and connectivity between) pyramidal cells and interneurons in PScz, this model was most suitable (although conductance-based models may furnish more useful pharmacological biomarkers). One key aspect of the model is that it models cortical dynamics only, and hence can model the γ and β peaks generated by superficial and deep cortical layers (respectively) (37-39), but it cannot reproduce an α peak without adding an additional (thalamic) input (40,41). Given the absence of group differences in α in the rsEEG, and for simplicity, we did not model this frequency band (as in previous studies (34,36)).
To simulate the power spectra in controls and PScz, we used the SPM12 (v7219) function spm_induced_optimise, which computes transfer functions (i.e. representations of cortical dynamics in the spectral domain) for parameters in the biophysical models in DCM. The simulated spectra were normalised in exactly the same way as the empirical spectra, subtracting the 1/f gradient using robustfit in Matlab. The standard priors in DCM were used ( Figure S1A).
We did not try to fit the microcircuit model to the empirical rsEEG data (although this is possible (36)) because the choice of sources in rsEEG data can be problematic. We instead used the microcircuit model to simulate the effects of various potential circuit pathologies in PScz, and compared the results to the pattern in the rsEEG. The circuit pathologies were based on reasonable hypotheses about microcircuit connectivity abnormalities in PScz   Model 4 -a <30% gain in interneuron self-inhibition, to model loss of synaptic gain (e.g. hypofunction of NMDARs) on interneurons.
Model 5 -a <30% gain in superficial pyramidal self-inhibition, to model loss of synaptic gain (e.g. hypofunction of NMDARs) on sp cells.

Data analysis and modelling: mismatch negativity
The MMN data were plotted (using electrode Fz, as is standard practice (42)) as groupaveraged waveforms for the standard and deviant tones separately ( Figure S2B) and as the traditional 'difference' waves, i.e. deviant-standard waves, illustrating the mismatch effect ( Figure 3A). Group differences in these waveforms were assessed using t-tests at each timepoint using an α of p<0.05 (uncorrected for multiple comparisons), and are displayed on the figures.
To more robustly illustrate the mismatch effects and group differences, incorporating other electrodes and cluster-based correction for multiple comparisons, the MMN sensor-level data were analysed in SPM12, after smoothing using an 8x8x8 mm FWHM Gaussian kernel ( Figures 3B, 3C and S2C).
The sources of the MMN response are well characterised, e.g. using conjoint EEG-fMRI studies (43). These are bilateral sources in primary auditory cortex (A1), secondary auditory cortex (superior temporal gyrus, STG) and inferior frontal gyrus (IFG). To ensure these source locations were adequate for our modelling, we source localised the MMN data using were close to the standard sources ( Figure S2D). Note that this source reconstruction was not used for the DCM analysis, however: source reconstruction (using the same priors) was performed as part of DCM model inversion. Joint optimisation of source locations and model parameters performs better than separate optimisations because it suppresses posterior correlations between source location parameters and biophysical parameters such that they explain away different parts of the data.
Standard DCM analyses estimate subject-specific parameters in four microcircuit connections (G): sp-sp, ii-ii, ii-sp and ii-dp. We estimated two extra G parameters given the strong possibility of differences in connectivity from pyramidal cells to interneurons in PScz (32), i.e. sp-ii and dp-ii. To avoid overparameterising the model, however, we fixed parameters elsewhere. First, we estimated whole group posterior means for synaptic delays (D) and time constants (T) and then used these as priors, fixing the D and T parameters to these values ( Figure S1A) in most models. Insodoing, we are not claiming that D and T are unlikely to be altered in PScz. We are merely assuming that there are very likely to be differences in G parameters in PScz, and that we wish to give the model the best chance of detecting them by reducing its degrees of freedom elsewhere. Second, we constrained all estimated G parameters to be the same (within but not between subjects) in every cortical area except for sp-sp self-inhibition. There are two reasons to adopt this approach: i) it regularises the model by allowing all cortical areas to inform the same parameters, reducing the overfitting that would occur when estimating six times as many parameters from one sixth of the amount of data per parameter, ii) post mortem studies of interneurons in PScz (33,44) indicate that pathology is broadly distributed across the cortex rather than focused in a single area, as one might expect from a disorder involving many genes of small effect. Nevertheless, NMDAR (and other receptor) subtypes are differentially distributed along a cortical hierarchical gradient (45), so if for example certain subtypes -preferentially located in higher or lower hierarchical areas -are more affected in PScz (46), the model could detect these spatial differences in sp-sp connectivity. Furthermore, sp-sp connectivity could also differ between standard and deviant conditions, permitting the model to infer condition-specific neuromodulatory effects on pyramidal cells (e.g. from NMDA (47) or cholinergic (48) receptors). Given we were applying a brain-wide constraint to most G parameters, we used the whole group posterior means ( Figure S1A) as priors, to reduce model fitting problems.
In total, six microcircuit models were evaluated (shown in Figure 3D): Models 4Ga, 4Gb, 4Gc, 4Gd -these four models estimated four microcircuit connectivity (G) parameters each, consisting of sp-sp and ii-ii self-inhibitory connections and then different combinations of two of the four connections between sp, ii and dp cells. Delays (D) and time constants (T) were fixed to their group posterior estimates (see Figure S1A).
Model 6G -this model estimated all six G parameters of interest, but fixed D and T to their group posterior estimates.

Model 6G,D,T -this model estimated all six G parameters of interest and also D and T (using their group posterior estimates as priors).
Each of these models used the same macroscopic structure ( Figure 2C, right), i.e. forward and backward connections linking areas in adjacent hierarchical levels, and lateral connections linking bilateral areas at the same level (not shown). Note that only forward, backward and self-inhibitory connections could show condition-specific differences between groups, i.e. differences between standard and oddball event related potentials.
MMN model fitting was performed using the 'spatial' or 'IMG' forward model in DCM for evoked response potentials (ERPs). DCM fits the first up to eight modes of the prior predicted covariance in sensor space (49) ( Figure S3A). In practice, the first 3-4 modes usually capture most of the variance (as in Figure S3A, right), so R 2 values generated for the MMN data are based only on the first four modes. A boundary elements head model (50) was used to approximate the brain, cerebrospinal fluid, skull and scalp surfaces.
Fitting non-linear neural mass models to EEG data is an ill-posed problem, and can lead to difficulties in optimisation, such as models getting stuck in local optima. Empirical Bayes for DCM (51) can circumvent this issue by performing model fitting iteratively, each time using the group level posteriors over parameters as priors for each subject's parameters in the next model inversion, yielding more robust and efficient parameter estimates (52). In practice, it substantially improved model fit: model inversions performed with and without it had R 2 values of ~0.7 and ~0.6 respectively. Local optima can also be avoided by fitting one fully parameterised model and then pruning unnecessary parameters using Bayesian model reduction (see below), rather than fitting many models with reduced numbers of parameters which are more prone to local optima problems (53). One subject failed model inversion (added to 17 who failed preprocessing: Table S1).
Following model inversion, models were formally compared using Bayesian model selection (54) ( Figure 3E). We used the protected exceedance probability (55) as the metric of success: i.e. the probability that any one model is more frequent than the others, above and beyond chance. The R 2 for the first four modes of the winning model was also computed ( Figure 3E).
Group differences between Con and PScz or Rel in their R 2 values were also assessed using ranksum tests (as the distributions were skewed; Figure S3C).

Group differences in parameters and relationships between parameters and other measures
were analysed using Parametric Empirical Bayes (PEB; spm_dcm_peb) (53,56). PEB is a hierarchical Bayesian version of a general linear model that can assess which DCM parameters differ between groups or covary with a continuous measure. One critical advantage of PEB (over performing classical statistical tests on parameters) is that it takes account of not just each parameter's expectation but also its covariance: therefore parameters that cannot be estimated with high confidence (e.g. in models that fit less well) contribute less to the inference. PEB is superior to standard model inversion at both selecting the correct model and estimating the model parameters (using simulated data) (51).
A second key aspect of the PEB procedure is the use of Bayesian model reduction (57) and Bayesian model averaging (53) (spm_dcm_peb_bmc). In essence, these steps prune away unnecessary parameters and then obtain posterior estimates for the remaining parameters by averaging over the remaining models 1 . The Bayesian model reduction procedure allows hypotheses about model parameters to be formally tested. This is done by first defining a model space of different parameter groupings that may best describe the modelled effect. In the MMN analysis, two model spaces were used: the first asked which combinations of extrinsic or self-inhibitory connections could best account for the mismatch effect ( Figure   S4A, left), and the second asked which combinations of intrinsic (microcircuit) connections could best account for the overall group differences across conditions ( Figure S4B, left).
These model spaces were based on the assumption that microcircuit parameters would not differ between conditions, but messages passed between areas (or pyramidal gain in specific areas) may well do so. In Figure S4A the first model (column) is a 'null' model, the second contains only 'forward' connections, the third contains only 'backward' connections, etc. The 1 Note that a correction for multiple comparisons is only appropriate when using classical statistics to reject the null hypothesis. In Bayesian inference (including PEB) there is no intimation that any hypothesis is rejectedthe inferences are based upon the relative evidence or marginal likelihood for the hypotheses in question. In other words, if an effect has a posterior probability of 95%, this means there is a 5% possibility that the effect is not necessary to explain the data. Applying a threshold to a posterior probability (as in this paper) is a device to simplify the reporting of effects for which there is very strong evidence. The number of effects for which there was strong evidence may be small or large -and does not affect the posterior probability reporting threshold (or require a correction for the number of effects for which evidence was assessed). whole model space consists of different combinations of forward, backward, and selfinhibitory connections in either temporal or frontal areas.
Following Bayesian model reduction, the posterior probabilities of these models accounting for the group difference effect (PScz > Con) are plotted in the middle. This shows that the PScz > Con effect is best described by Model 7, containing only the left and right IFG selfinhibitory connections. Instead of simply using the parameters estimated from the winning model (as there may not always be a clear winner), posterior parameter estimates are obtained by averaging over all models, weighted by their model evidence (Bayesian model averaging).
Note that if there is no evidence for some models, parameters unique to those models will be redundant: these parameters are eliminated following Bayesian model reduction. Conversely, parameters that are present in all probable models are highly likely to explain the group

Data analysis and modelling: 40 Hz auditory steady state response
The timecourses of the 40 Hz ASSR waveforms in electrode Fz are plotted in Figure 4A: ~40 Hz oscillations are superimposed on more sustained baseline changes (known to occur during auditory click trains (58)). PScz and Rel response baselines diverge (t-tests at each timepoint) from Con around 250 ms, but we restricted our modelling to well-replicated group differences in ~40 Hz power. Each subjects' EEG data were transformed into measurements of phase and power in the frequency domain using a Morlet wavelet in spm_eeg_tf. For the Fz sensor-level analyses, the spectra were normalised in order to assess the power at 40 Hz relative to the background 1/f slope. For normalisation, we again subtracted the 1/f gradient in log space (computed using robustfit in Matlab applied to the 10-30 Hz and 50-70 Hz ranges, on a subject-specific basis): example gradients of group averaged data from electrode Fz are shown in Figure S2E. γ peak frequency was the frequency in which the maximum (normalised) power within the 35-45 Hz range occurred (plotted in Figure 4B). The unnormalised group averaged time-frequency plots are shown in Figure S2F, and the normalised plots in Figure 4C. Permutation testing assessed whether the number of observed The neuronal model for DCM had to be adjusted in order to model the high amplitude ~40 Hz peaks in the virtual electrode data from left and right A1 (dotted lines in Figure S3B). To fit these peaks we introduced a Gaussian 'bump' at 40 Hz (of width w≤4 Hz) into the background neuronal noise, to reflect the specific 40 Hz input (similar to a previous model (41) using a 10 Hz bump added to neuronal noise to simulate a thalamic input). As in the MMN analysis, we estimated group mean values for T (time constants) and also S (the slope of the sigmoid activation function) and J (the contribution of neuronal populations to the EEG signal) in an initial model inversion, and then used these as empirical priors ( Figure   S1A). These changes meant that sp and ii time constants were slightly shorter, there was increased gain in the transformation from membrane potential to firing rates, and that ss cells (which receive the 40 Hz thalamic drive) made more contribution to the measured EEG signal (respectively). We fixed T to its priors because without being fixed, they often became biologically implausible. We did not fix D (delays) because these might have important effects on peak frequencies in this paradigm. As in the MMN analysis, we estimated six microcircuit (G) parameters but constrained all except sp-sp self-inhibition to be identical (within subjects) in A1 bilaterally. Last, we increased the precision of the data (hE), which was necessary to force the model to fit the unnatural peaks in the 35-50 Hz range. We attempted to fit a wider range of frequencies , but the model fitted the lower frequencies at the expense of the 40 Hz peak, ignoring the key data feature, so the narrower  Figure S3B. The model fits reasonably well -it often cannot capture the full amplitude of the ~40 Hz peak, but usually has its peak at the right frequency: one exception is subject 0254 (bottom right), whose unusual peaks at 45 Hz are too far from 40 Hz for the model to accommodate. R 2 values for the winning model are shown in Figure 4E (right), and group differences between Con and PScz or Rel in their R 2 values were also assessed ( Figure S3C). R 2 values were more variable in the 40 Hz ASSR than the other paradigms, because in the 40 Hz ASSR the model was trying to fit a small data feature very precisely. 19 subjects failed model inversion (added to 12 who failed preprocessing: Table S1).
The model space used for the analysis of A1 microcircuit (G) parameters within PEB is shown in Figure S5A (first). The first two rows denote sp-sp self-inhibition on the left and right, and the next five rows containing white boxes denote the remaining microcircuit parameters which were identical bilaterally. The (column-wise) model space is constituted by different combinations of sp and ii self-inhibitory parameters and ii input and output parameters. Following Bayesian model reduction, no model of the group differences dominates (middle), however all models with any evidence for them contain sp-sp selfinhibition parameters, so the posterior probability (after Bayesian model averaging) that these parameters contribute to the group difference is very high (right).

Data analysis and modelling: resting state fMRI
We restricted the rsfMRI analysis to being directly comparable to the MMN and 40 Hz ASSR analyses. We therefore chose to analyse effective connectivity in exactly the same network involved in the MMN (and 40 Hz ASSR, in part), i.e. bilateral A1, STG and IFG. The sources used for the MMN lie within areas A1, A4 and 44 of the Glasser parcellation (25) (Figure 1), so these were chosen to be the nodes of the network. We did not use subject-specific sources from the MMN source reconstruction, because the MMN sources are well-established, and if a subject's EEG source reconstruction lay outside these areas, it is arguably more likely that their EEG source reconstruction was inaccurate than that a different cortical network activated in their MMN. Thus, using the same (individually parcellated) rsfMRI sources for each subject took account of individual variation in anatomy whilst ensuring the rsfMRI analysis was robust and reproducible.
For the DCM analysis we used the movement-scrubbed rsfMRI data with global signal regression. Although global signal likely contains neuronal signal that is relevant to schizophrenia (61), we removed it in order to reduce global physiological confounds as much as possible. No subjects were excluded due to low SNR, although doing so (for the ten subjects with SNR<25) did not change any of the results. Model structure was assumed to be as it was in the MMN analysis (Figure 1 values (and their differences between groups) are shown in Figure S3C.
Again, PEB was used to analyse group differences in parameters. The model space for the PScz > Con analysis is shown in Figure S6A (left). It comprises different combinations of forward, backward and self-inhibitory connections. There is no clear winning model (middle) but all likely models contained bilateral IFG self-inhibitory connections, hence these connections are inferred to be highly likely to explain group differences (right).
For comparative purposes, in addition to the spectral DCM analysis, we performed a standard functional connectivity analysis ( Figure S6B), i.e. computing the Pearson (zero lag) correlation between two parcellations' timeseries. We also computed their (zero lag) covariance ( Figure S6C), because comparing correlations in Con and PScz can be problematic as their variances may differ (61), and the corresponding zero lag covariance in the neuronal states estimated by DCM ( Figure S6D). To see whether spectral DCM selfinhibition parameters might be recapitulated in basic data features, we also assessed group differences in the coefficient of variation (standard deviation normalised by the mean, computed using an unbiased method) of BOLD signal in each node ( Figure S6E if one group is older (66) or taking (antipsychotic) medication, or if one group moves more (respectively); ii) spectral DCM connectivity includes non-zero lags, which indicate direction of connectivity (62,67); iii) and -most importantly -spectral DCM also estimates selfinhibition within cortical regions.
We also performed the same DCM and PEB analyses without global signal regression, to ascertain the effects of this preprocessing step on the results ( Figure S7).

Data analysis and modelling: parameter sensitivity and cross-paradigm analysis
Numerous studies have shown that DCM can accurately recover group-specific parameter changes when the ground truth is known: for example, inferring the known physiological effects of isoflurane anaesthesia from rodent cross-spectral data (68), inferring antiNMDA receptor antibodies have specific effects on NMDA receptor parameters from human rsEEG (36), and reliably recovering parameter values (using DCM and PEB) from simulated MMN EEG data (52) and from simulated rsfMRI data (62).
Nevertheless, we performed parameter sensitivity analyses on the six estimated microcircuit (G) parameters to address two potential weaknesses in the MMN and 40 Hz ASSR analyses.
In the MMN analysis, constraining condition-specific (mismatch) effects to only sp-sp selfinhibition amongst the G parameters -along with extrinsic connections -might miss such effects in other G parameters. In the 40 Hz ASSR analysis, the data were afforded high precision during model fitting in order to force the model to fit the unnatural 40 Hz peak: this may lead to overfitting of other data features, however, and spurious results. We therefore selected two subjects from each paradigm with EEG responses typical of the group average. We simulated virtual electrode ('LFP') data from a single cortical area using spm_induced_optimise and either spm_gen_erp (for the MMN) or spm_csd_mtf (for the 40 Hz ASSR). In each case we used these subjects' posterior G parameter estimates, varying each G parameter by ±30% in turn; in the MMN, we also fixed delays (D) and time constants (T) to the values used in Figure S1A. The simulated data are shown in Figure S8, with the key data features that the parameters showing group differences ought to explain circled in red.
We used the PEB framework to assess whether the abnormal parameters in the PScz group related to each other across modalities, i.e. across EEG and fMRI. If so, this would support the use of multiple paradigms in combination to estimate these 'biomarker' parameters. We selected the left and right IFG self-inhibition parameters from the MMN and left and right A1 self-inhibition parameters from the 40 Hz ASSR, and their counterparts from the rsfMRI. For each of the four rsfMRI parameters separately, we used the PEB framework to reveal which of all MMN or 40 Hz ASSR parameters related to it. Thus the posterior expectation and covariance in the EEG parameters, but only the expected value of the rsfMRI parameter, was used by the PEB model. (We reasoned that the rsfMRI parameters were likely estimated with greater confidence, given their simpler models and more consistently high R 2 values compared with the EEG models ( Figure S3C), hence it would be better to include the EEG parameters' covariance in the PEB model.) Two of the four analyses yielded significant results: R IFG self-inhibition across rsfMRI and MMN ( Figure S9A), and R A1 selfinhibition across rsfMRI and 40 Hz ASSR ( Figure S9B).
For comparative purposes, we also performed classical analyses of the relationships between the same rsfMRI and EEG parameters. The Pearson correlations between the rsfMRI and 40 Hz ASSR or MMN self-inhibition parameters are shown in Figure S9C. Both EEG paradigms -especially the MMN -show some clustering of the parameters around their prior values (0 in the 40 Hz ASSR, and 0.5 in the MMN). Posterior parameter values that are further from the priors are typically estimated with greater confidence. The advantage of the PEB framework is that it will (optimally) upweight the contribution of these more confident parameters to group-level inference.

Excluding benzodiazepine effects on rsEEG
The group differences in rsEEG power were unlikely to be driven by participants taking benzodiazepines, as daily use of these drugs was an exclusion criterion, and removal of these subjects (5 Con, 16 PScz) did not substantially alter the effects, although Bonferronicorrection now shifted the gamma difference to trend-level significance. Compared with Con  Figure S7A) to reduced synaptic gain within IFG. Given this, and given GSR in PScz affords a relatively uniform transform of the data (61), the results with GSR are probably more reliable, especially given their consilience with the EEG results ( Figure 6).

Relationships between rsfMRI and EEG data features and self-inhibition parameters, and the possibility of 'model-based biomarkers'
We assessed relationships of rsfMRI data features to MMN amplitude and 40 Hz ASSR Only bilateral A1 and R IFG would survive correction for 6 comparisons, however. None of these relationships were significant in Con (all P>0. 2) or Rel (all P>0.1). There were no relationships between 40 Hz ASSR γ power and rsfMRI BOLD coefficient of variation in A1 in any group (all P>0.07).
Finally, in post hoc analyses, we asked whether self-inhibition parameters -that were altered in PScz -had predictive validity across paradigms, thus licensing their use as 'model-based biomarkers'. We used PEB to assess relationships in PScz between self-inhibition parameters in IFG in the MMN, or in A1 in the 40 Hz ASSR, and their corresponding self-inhibition parameters in the rsfMRI, in participants whose EEG and rsfMRI were <100 days apart (MMN n=37, 40Hz ASSR n=41), as some participants underwent EEG and rsfMRI many months or even years apart. Within PScz, there was evidence of associations between selfinhibition in R IFG in the MMN and rsfMRI (P>0.95; Figure S9A) and in R A1 in the 40 Hz ASSR and rsfMRI (P>0.95; Figure S9B). However, there were no such relationships in the left hemisphere or in Con.
To explore these relationships further we assessed their Pearson correlations ( Figure S9C). This clustering of EEG estimates around their prior values must be finessed (by using more informative data or paradigm design) to assess cross-paradigm self-inhibition reliability in all participants. If parameters show cross-paradigm validity, hierarchical PEB can be used to estimate them from multiple paradigms simultaneously (51), and thereby provide modelbased biomarkers for psychosis.

Reduced self-inhibition relates to Digit Symbol performance in Con across three tasks
We also assessed relationships of DCM parameters to a cognitive measure: Digit Symbol task performance ( Figure S10). In the MMN, across both conditions, reduced self-inhibition in bilateral IFG in Con was associated with greater information processing speed (i.e. Digit Symbol Substitution Task score; P>0.99 and P>0.95, Figure S10A), but there were no such relationships in PScz (all P<0.95).
In the 40 Hz ASSR, in Con and PScz, Digit Symbol performance was associated with reduced pyramidal self-inhibition in right and left A1 respectively (both P>0.99) and reduced deep pyramidal to inhibitory interneuron (dp-ii) connectivity (P>0.99; Figure S10B).
In the rsfMRI, Digit Symbol performance was associated with decreased R IFG selfinhibition in Con (P>0.99; Figure S10C) but not in PScz, who instead showed decreased L A1-STG and increased L STG-A1 connectivity (both P>0.99, not shown).
Thus, processing speed (i.e. Digit Symbol performance) in Con was associated with disinhibition (or increased synaptic gain) across three paradigms, which speaks to the predictive validity of the synaptic biomarkers furnished by DCM. This possible relationship between cognitive performance and synaptic gain complements a recent rsfMRI DCM study of healthy adults (n=602) that found fluid intelligence was related to decreased self-inhibition in dorsal attention and salience networks (66). R IFG self-inhibition's relationship to Digit Symbol performance (in Con) also mirrors the finding that the global functional connectivity of a nearby region of left lateral prefrontal cortex relates to fluid intelligence (71). In PScz, however, self-inhibition and cognition were only associated in R A1 (in the 40 Hz ASSR). Interestingly, a rsEEG DCM study in Down syndrome (n=36) found intelligence was associated with self-inhibition in V1 (72). Why cognition relates to self-inhibition in primary sensory regions in these groups is unclear: it may be estimated more efficiently there.      conditions. The P values denote ranksum tests (not corrected for both comparisons).

B -Standard (above) and deviant (below) MMN responses for Con (n=94), PScz (n=96) and
Rel (n=42) at electrode Fz. Red bars denote group differences in Con vs PScz, green bars denote group differences in Con vs Rel, both P<0.05 (t-tests at each timepoint, uncorrected).
C -Left: The lower plot shows the overall Con > PScz effect at sensor level, displayed at   Right: Interestingly, the Rel > Con contrast indicated a loss of sp input into interneurons but also a disinhibition of sp cells -the former effect was also seen at P>0.75 in PScz, but the latter effect was reversed in PScz at P>0.75: see Figure 3C for an analysis combining all three groups. All effects shown are also present without the addition of age, sex and smoking covariates (P>0.95). 40 Hz ASSR microcircuit parameters: Rel > Con 40 Hz ASSR microcircuit parameters: Scz > Con     reproduce key data features in the MMN and 40 Hz ASSR paradigms. In each paradigm, two subjects representative of the group average were selected, and their posterior G (intrinsic connectivity) parameters used to simulate data. In each simulation, one G parameter in turn was progressively changed by ±30% -results of increases and decreases are shown in orange and purple, with darker shades representing bigger changes. The left and right panels correspond to the two subjects in each paradigm, and the six plots in each panel to the six G parameters: sp-sp, ii-ii, ii-sp, ii-dp, sp-ii and dp-ii connectivities. Note that the first two are self-inhibitory connections (synaptic gain).
A -Sensitivity analysis in the MMN paradigm. The ERP data showed a reduced MMN amplitude ( Figure 3A) and P100 amplitude ( Figure S2A) but no significant change in MMN latency in PScz; in PScz, sp-sp was increased ( Figure 3F). In both subjects' simulations, increasing sp-sp (i.e. reducing synaptic gain) reduces the amplitude of the negative deflection (circled in red) without changing the latency, as it also does to the P100. Increasing ii-ii has similar but smaller effects on these amplitudes, but shortens the latency. Decreasing ii-sp likewise reduces the amplitudes, but lengthens the latency. The remaining parameters cause more substantial changes in the latter half of the waveform. The waveform itself is compressed in time compared with empirical data, because extrinsic connections -which cause a ~10 ms delay per connection -are not included in the simulation of a single area.
B -Sensitivity analysis in the 40 Hz ASSR paradigm. The power spectral data showed a reduction in ~40 Hz power ( Figure 4C). In both subjects' simulations, increasing sp-sp and decreasing sp-ii (the parameter changes in PScz vs Rel and PScz+Rel vs Con respectively - Figure 4G) both reduce ~40 Hz power (circled in red). Reducing ii-sp has similar but less powerful effects on sp-ii. Changing ii-ii changes the peak frequency, but its effects on power are only prominent at precisely 40 Hz. Dp-ii and ii-dp connections have minimal effects.

Figure S9 -Resting state fMRI & EEG DCM parameter relationships in PScz
All analyses are conducted on PScz whose rsfMRI and EEG paradigms were recorded <100 days apart (for MMN, n=44; for 40 Hz ASSR, n=40). A -In PScz, increased self-inhibition in R IFG (estimated from rsfMRI) relates to both loss of superficial pyramidal cell gain in R IFG specifically and a net increase of inhibition of pyramidal cells (i.e. stronger pyramidal inputs to interneurons and stronger returning connections except for ii-sp, whose decrease is smaller than the sp-ii increase) in the MMN analysis. All effects (except the ii-dp increase) shown are also present without the addition of age, sex and smoking covariates (P>0.95).
B -In PScz, increased self-inhibition in R A1 (estimated from rsfMRI) relates to increased sp self-inhibition in R A1 specifically and stronger sp input to interneurons. It alsocounterintuitively -relates to greater self-inhibition of and weaker dp input to interneurons.
All effects (except the dp-ii decrease) shown are also present without the addition of age, sex and smoking covariates (P>0.95).
NB With the inclusion of chlorpromazine dose equivalents as a covariate the specific R IFG effect (in A) and R A1 effect (in B) are lost, however, and likewise if PScz with rsfMRI and EEG paradigms >100 days apart are included.
C -These plots show the correlations between self-inhibition parameters across EEG and rsfMRI paradigms within PScz. The left two plots show correlations between 40 Hz ASSR and rsfMRI estimates of self-inhibition in bilateral A1, the right two plots show correlations between MMN and rsfMRI estimates of self-inhibition in bilateral IFG. Relationships in R A1 and R IFG approached statistical significance, and are analysed in more detail using PEB in panels A and B above. There were no correlations between the variance of these parameters across these combinations of paradigms (all P>0.3). B -PEB analyses in Con (left) and PScz (right) groups, showing Digit Symbol score relates to reduced sp self-inhibition (i.e. increased synaptic gain), and lower dp-ii connectivity.
C -PEB analysis of the rsfMRI showed that in Con, Digit Symbol score related to reduced self-inhibition in R IFG. PScz results are not shown.
The relationship between Digit Symbol (i.e. processing speed) performance and reduced selfinhibition is remarkably consistent across paradigms in Con, but not in PScz. All effects shown are also present without the addition of age, sex and smoking covariates (P>0.95), and also with inclusion of chlorpromazine dose equivalents as a covariate.