Participants
In total, 14 patients with schizophrenia (5 female) and 14 gender and age matched healthy controls (HCs) participated in the current study. Table 1 summarizes the demographic characteristics of our sample. Parts of the HC sample have been previously published32. For this current analysis, HCs were selected to match the schizophrenia sample regarding age and gender at the group level (no pair-wise matching). Participants with schizophrenia were recruited through the Department of Psychiatry, Psychotherapy and Psychosomatics of the RWTH Aachen University Hospital. HCs were recruited via public flyers and the RWTH Aachen University community. All participants met the following inclusion criteria: age between 18 and 55 years, no neurological diseases, normal or corrected-to-normal vision, and all fulfilled MR-scanning criteria. In addition, for the patient group schizophrenia had to be the primary psychiatric diagnosis. Comorbid affective disorders and schizoaffective disorders were excluded. The diagnosis was based on the International Classification of Diseases (ICD-10) by an experienced psychiatrist or psychologist. In addition, as part of the screening procedure of the study, the diagnosis was ascertained by a psychologist using the German version of the Structured Clinical Interview for DSM-IV Axis I Disorders (SCID-I). HCs were excluded if they had any history of psychiatric disorders. The presence of psychiatric diagnoses of first-degree relatives was an exclusion criterion as well.
During participation, all patients were taking psychotropic medication. The most frequently used psychotropic class was atypical antipsychotics (n = 19), followed by antidepressants (n = 5), and anticonvulsants (n = 1), and typical antipsychotics (n = 1). Most patients took multiple kind of medication. The mean antipsychotic dose in olanzapine equivalents was 21.39 mg per individual (converted according to the Defined Daily Doses Method by Leucht et al.33). See Table 1 for further characteristics of the patient group, including duration of illness, age of onset and PANSS scores. All but one participant were right-handed as assessed by the Edinburgh Handedness Inventory.
The experimental protocol was approved by the local Ethics Committee of the RWTH Aachen University Hospital. All participants gave oral and written informed consent, and no adverse events occurred. Participants were reimbursed with 10 Euro per hour for study participation.
Materials
Neuropsychological and psychopathological assessment
The neuropsychological tests administered included the Trail Making Test (TMT) A and B, a crystallized verbal intelligence estimation (Wortschatz-Intelligenztest, WST), and the digit span task of the Wechsler Memory Scale-Revised.
The TMT-A assesses psychomotor speed, requiring participants to connect 25 numbers in ascending order. TMT-B evaluates cognitive flexibility and working memory, with participants connecting numbers and letters alternately. An executive control measure was derived from the ratio of part B to part A34. Verbal intelligence was assessed using the WST, which consists of 42 items, with six words each, where each item contains five nonsense words and one existing word. The participant is asked to identify the existing word (in German language). A total of 31 correct identifications corresponds to a mean verbal IQ of 10135. Memory performance (forward and backward) was assessed with the Wechsler Memory Scale—revised (WMS-R). In the forward task, participants repeat a sequence of numbers in the order verbally presented, starting with three digits and progressing through six trials. The backward task requires participants to repeat sequences in reverse order.
Self-reports assessing psychopathology that were administered in both patients and controls included German versions of the Beck Depression Inventory-II (BDI-II36), as well as the Toronto Alexithymia Scale (TAS-20)37. The BDI-II comprises 21 multiple-choice questions to assess the severity of depressive symptoms. It is one of the most commonly used psychometric measures for assessing depression, demonstrating high validity with an internal consistency of alpha = 0.84 in a German sample of individuals diagnosed with depression38. The TAS-20 consists of 20 items, which are to be rated on a 5-point Likert scale (1 = strongly disagree, 5 = strongly agree), measuring difficulties in identifying and describing emotions. Moreover, the Global Assessment of Functioning (GAF) was used to measure the general level of functioning in participants. Mental, social, and occupational functions are conceptualized on a continuum from mental health to illness. The GAF scale is divided into 10 functional levels, each with 10 points, ranging from 100 (highest level of performance) to 1 (lowest level of performance). The rater assigns a single score that most accurately reflects the patient’s general level of functioning, taking into account both the severity of symptoms and the extent of impairment. If symptom severity and degree of impairment fall into different categories, the lower level is chosen.
In addition, psychopathological ratings have been performed (by BSH), namely for depression and schizophrenia symptoms, the Hamilton Rating Scale for Depression (HAM-D) and PANSS (Positive and negative affect scale for schizophrenia) were administered to patients. The PANSS allows the assessment of positive, negative and general psychopathological symptoms with a total of 30 items, each rated on a 7-point scale. It is administered in the form of an interview.
Stimuli and task
Stimuli and tasks are reported elsewhere in more detail30. The experimental task consisted of targets and flankers (hereafter referred to as crowds) of face stimuli. Specifically, crowd stimuli consisted of a circularly shaped crowd of six colored pictures of faces (3 female, 3 male), all either expressing the same emotion (happy, neutral or fearful) or the faces were made unrecognizable and pixelated (scrambled). The face crowds defined the four experimental conditions of the variable “crowd” (happy, neutral, fearful, scrambled). The face crowd stimulus was presented for 300 ms, followed by the target stimulus, which was presented for 500 ms. The target stimulus consisted of the same circularly shaped crowd as the preceding flanker stimulus. In addition, the target stimulus also depicted another face in the center of the crowd (the target face depicted one of the three emotions: happy, neutral, fearful) and was followed by a blank screen for 450 ms. The participants’ task was to identify the target emotion on a seven-point scale ranging from extremely fearful through neutral (middle icon) to extremely happy or vice versa. The final decision had to be made within 3500 ms. A jittered fixation cross was shown for a mean duration of 875 ms before the next trial. An overview of the entire task sequence is shown in Fig. 1.
Each trial commenced with the display of a fixation cross, indicating the position of the upcoming target and facilitating gaze fixation. After 500 ms, a facial crowd stimulus was introduced (for 300 ms), followed by the replacement of the fixation cross with a target face, presenting the crowd-target combination for 500 ms. A subsequent blank screen (450 ms) was included to prevent any overlap between response-related movement and the signal of interest. Following this, a response screen was presented for 3500 ms.
Each face crowd condition (happy, neutral, fearful, scrambled) comprised 120 trials (each stimulus shown twice), totaling 480 trials. The order of trials was pseudo-randomized. The experiment was divided into four parts of 120 trials each and was presented via MR-compatible goggles using Presentation® (version 14, Neurobehavioral Systems Inc., San Francisco, CA). Responses were given in the form of button presses, recorded via an MR-compatible response system. All stimuli were selected from a standardized stimulus set and were comparable in terms of age, emotion intensity, luminance, and emotion valence32,39.
Procedure
Before participating, all participants provided written informed consent after having received detailed oral and written information on the study. All procedures were in accordance with the Declaration of Helsinki and were approved by the University’s Ethics Committee (Ethics Committee at the RWTH Aachen Faculty of Medicine).
Simultaneous EEG-fMRI data acquisition
In the current study, fMRI data was acquired using a Siemens 3 Tesla MR scanner (Siemens TimTrio®, Siemens Medical Systems, Germany). An echo planar imaging sequence was used with the following parameters: TR/TE = 2000/30 ms, 76˚ flip angle, 3.125 × 3.125 × 3.4 mm3 voxel size, 64 × 64 matrix, 200 × 200 mm2 FOV, 33 3.4 mm-thick axial AC-PC slices with 0.51 mm gap. Data were recorded in ascending sequential slice acquisition with whole brain coverage using a standard 12-channel head coil.
EEG data were simultaneously acquired with the fMRI data using an MR-compatible EEG system (BrainAmp MR, Brain Products GmbH, Gilching, Germany) at a sampling rate of 5000 Hz using MR-compatible EEG-caps (Easycap GmbH, Germany) with 64 Ag-AgCl electrodes (extended 10–20 system). An additional electrocardiogram (ECG) electrode was placed under the left collarbone of the participant. Electrodes FCz and FPz served as the recording reference and ground channel, respectively. Impedance levels were below 10 kΩ at start of the recording (except ECG). From the amplifier, the digitized data were transmitted via fiber optic cables outside the scanner room to an USB interface located in the control room. To facilitate the removal of MR-induced artifacts from the EEG data, the sampling clocks of the EEG and MRI systems were synchronized by a SyncBox (Brain Products GmbH, Gilching, Germany).
A trained specialist positioned the EEG caps on the participants. To ensure valid standard positions, the electrode Cz was placed halfway between the nasion and the inion, and was right-left-centered.
Data processing and analysis
We describe the (pre-)processing steps of each modality (behavioral, EEG, fMRI) separately below. The pipeline of our imaging fusion analyses may be found in Fig. 2. In fMRI-EEG fusion, we leverage the high spatial precision of fMRI (identifying regions related with face emotion processing) and the high temporal resolution of EEG (capturing early neural activity related to face emotion processing) to offer complementary and informative constraints on model parameters. Briefly, multimodal fusion was partly based on a pipeline described by Wei and colleagues40:
-
1.
Whole-brain fMRI analysis: We identified brain regions significantly activated by our face emotion processing paradigm, using whole-brain analysis of the acquired fMRI data (unimodal analysis using a flexible factorial GLM in SPM).
-
2.
Coordinate Extraction: MNI coordinates of the activated regions were extracted. These identified regions form the network structure for the subsequent DCM analysis.
-
3.
Network Node Definition for DCM Analysis: We conducted a DCM analysis of the simultaneous acquired EEG data by setting the prior source locations to the activated regions identified in step 1/step 2, in order to obtain the posterior densities of the neuronal parameters. A neural mass model was allocated to each of the regions. These models were connected via possible forward and backward connections. ERPs over 0–400 ms post-target were fitted to the models. The effects of each task condition (TAF, TAH, TAN) on the ERPs were contrasted against the TAS condition, and modulatory effects on these contrasts were modeled for all extrinsic connections (i.e., between the sources) of the given model structure.
-
4.
Bayesian Model Selection (BMS): We applied BMS to identify the best-fitting model (that most likely generated the ERPs). The results showed that Model 6 provided the best fit among the six competing candidates.
-
5.
Post-hoc analyses on the winning model.
Step 1: Whole-brain fMRI analysis: We conducted an unimodal flexible factorial whole-brain analysis, examining the main effect of the task condition (F-contrast). This revealed regions significantly activated by the experimental paradigm, thresholded at p < 0.001 with cluster-level correction at k = 34, corresponding to a cluster-level FWE of p < 0.05. Step 2: Coordinate Extraction: We extracted MNI coordinates of these activated regions, excluding the cerebellum. Step 3: Network Node Definition for DCM Analysis: Using the extracted coordinates, we defined nodes for the cortical network in Dynamic Causal Modeling (DCM) analysis. Six model structures were created, incorporating combinations of forward, backward, and bilateral connections among the four brain regions identified in the whole-brain fMRI analysis. Step 4: Bayesian Model Selection (BMS): We applied BMS to identify the best-fitting model. The results showed that Model 6 provided the best fit among the six competing candidates.
Behavioral data analysis
Behavioral variables of interest included the number of correctly identified target emotions, duration until the first response, and intensity ratings of target emotions (happy, fearful). Note that trials with responses <150 ms were discarded. The count variable “correctly identified target emotions” was transformed using Tukey’s Ladder of Power transformations to normalize the data. Next, we computed three linear mixed models using the R package “lme4” with “subject” as a random effect to estimate the above stated dependent variables. Independent, fixed effect variables were “group” (HC, ISZ), and the repeated measures variables “crowd” (fear, happy, neutral, scrambled) and “target” (fear, happy, neutral). Note that for intensity ratings, the ‘target’ variable consisted of only the fear and happy level. P-vales of post-hoc comparisons (two-sided) were adjusted using Tukey method, implemented in the “emmeans” R package.
EEG data preprocessing and analysis (unimodal)
Offline EEG data preprocessing was performed using the BrainVision Analyzer software (version 2.0, Brain Products GmbH, Gilching, Germany). The preprocessing pipeline included the following steps: 1) removal of MR-induced artifacts from the raw EEG signal using the sliding average procedure with 21 templates for the calculation of the correction template 4>; 2) removal of cardioballistic artifacts by a semi-automatic algorithm41: R-peak markers were set at highly correlated (r = 0.7) time points with above-threshold amplitudes (0.4–1.7) and inspected before template subtraction; 3) filtering of the EEG data using an infinite impulse response filter (IIR, 70 Hz, 48 dB slope); 4) and resampling to 250 Hz; 5) independent component analysis (ICA), employing the restricted biased Infomax algorithm, was utilized to effectively eliminate artifacts. Components indicative of eye blinks and movements, residual gradient, and muscle artifacts were identified and removed based on their power spectrum and topography. Finally (6) data were re-referenced to a common average reference. The subsequent processing steps were done in EEGLAB42, a MATLAB-based toolbox: A limited number of poor contact electrodes from individual measurements were screened for visual inspection and were subsequently substituted using spherical interpolation techniques.
The continuous EEG recording was then segmented from −1000 to 1750 ms time-locked to the onset of the four different conditions of the crowd stimuli. Baseline correction was applied to the epochs from −200 to 0. Segments from all conditions were visually inspected, and those containing muscular activities or non-physiological artifacts were rejected.
For the subsequent steps, only correct trials were considered. EEG segments were averaged for each subject and for each condition. Specifically, P1 amplitudes were extracted and averaged at electrodes P7/PO7 and P8/PO8 at 90–120 ms post-stimulus (i.e., crowd) onset. N170 amplitudes were extracted at electrodes P7 and P8 at 160–200 ms post-stimulus onset. In addition, of particular interest to the current study was to investigate the processing of the target stimuli after crowds. For this purpose, another baseline correction was applied to the data from 250 to 350 ms (post-crowd onset) before averaging. Subsequently, again P1 (at 390–420 ms) and N170 (at 470–520 ms) ERPs were extracted by averaging (condition- and subject-wise). The selection of electrode positions and time windows was guided by previous literature and topographic plots (see supplementary Figs. 1 and 2). Since the main interest of our study was to study the effect of task-irrelevant faces on the processing of subsequent target faces, we focus on reporting ERP amplitudes in response to targets stratified by preceding crowds. In addition, we report ERP amplitudes of crowd faces in supplementary materials. The P1 and N170 amplitudes of the targets after crowds (and crowd) ERPs were each analyzed with separate repeated measures ANOVAs. The models each included the within-subjects factors “crowd” (happy, neutral, fearful, scrambled) and “electrode site” (left, right) as well as the between-subjects factor “group” (HCs, ISZ). In cases of violations of the sphericity assumption, Greenhouse Geenhouse-Geisser correction was applied. Post-hoc comparisons were Bonferroni-corrected.
fMRI data preprocessing and analysis (unimodal)
The MATLAB-based toolbox SPM12 was used for preprocessing (https://www.fil.ion.ucl.ac.uk/spm/software/spm12/). Six motion parameters were estimated and used to realign the functional images to their mean. Subsequently, the functional images were co-registered to the structural image that was segmented into tissue components (gray matter, white matter, and cerebrospinal fluid) and normalized to the Montreal Neurological Institute (MNI) standard brain template. Functional images were smoothed with an 8 mm Gaussian kernel.
After preprocessing, parameter estimates reflecting the BOLD signal change for correctly answered target onsets within crowd conditions were calculated for each subject. In the following, these are referred to as target after fear (TAF), target after happy (TAH), target after neutral (TAN), and target after scrambled (TAS). In addition, for each subject, a separate regressor for all incorrect trials was included. Four T-contrasts were created, each testing for the main effect of the four conditions, which were further used for the group-level analysis. Subsequently, a flexible factorial design with the task conditions for the different groups and a subject factor was estimated on the second level. Deviations from sphericity were corrected for by variance components for within-subjects measures and heteroscedasticity between subjects and conditions. Finally, an F-test was carried out to test the effect of condition (i.e., main effect across groups).
Identifying regions significantly activated by the experimental paradigm, using whole brain (SPM) analysis of fMRI data was done in order to create a network architecture for subsequent DCM analysis (c.f., ref. 40). Results were thresholded at p < 0.001 and cluster-level corrected at k = 34, corresponding to a cluster-level FWE of p < 0.05.
fMRI-informed EEG DCM analysis
DCM is a hypothesis-driven approach to compare hypotheses about the mechanisms (in terms of neuronal coupling) that underlie the regional responses detected in conventional analyses43. Here, we applied DCM for ERPs based on a neural mass model aimed at elucidating source activity. This model integrates the dynamics of three cortical layers, including an excitatory subpopulation in the granular layer, an inhibitory subpopulation in the supra-granular layer, and a subpopulation of deep pyramidal cells in the infra-granular layer. Specifically, the DCM generates a predicted ERP by simulating the response of a network of interconnected sources to sensory input. Each source is modeled as a point source, consisting of the three subpopulations, each assigned to a specific cortical layer. Based on extrinsic connectivity rules44, a hierarchical cortico-cortical network is organized, comprising several cortical sources, each modeled using the neural mass model. The strengths of extrinsic connections among these cortical sources are estimated using Variational Bayes. All DCM steps were performed using SPM12.
Briefly, informed by the fMRI analysis of task condition effects (see above; F-contrast on the effect of condition as part of a flexible factorial design), we conducted a DCM analysis. Specifically, six brain regions (MNI coordinates) were identified by this F-contrast, including the bilateral occipital fusiform gyrus (OFG), left insula, V1, and bilateral cerebellum. Four of these regions—aligned with previous findings on emotional face processing—were selected for the subsequent DCM analysis. Using these four regions as prior coordinates, we constructed six models of cortical connectivity (see Step 3 in Fig. 2), each representing varying configurations of forward and backward connections among the ROIs. In the simplest model (model 1), an external visual input was entered at V1, which projected forward to the bilateral OFG, and from the left OFG to the left insula. In models 3 and 5, bilateral OFGs had backward connections to V1, with an additional backward connection from the left insula to the left OFG in model 5 only. Models 2, 4, and 6 also introduced lateral connections between the left and right OFG.
ERPs over 0–400 ms post-target were fitted to the models. Activity of each cortical source was modeled with a single equivalent current dipole method based on the individual forward head model45. The effects of each task condition (TAF, TAH, TAN) on the ERPs were contrasted against the TAS, and modulatory effects on these contrasts were modeled for all extrinsic connections (i.e., between the sources) of the given model structure. In other words, the task-specific connectivity effects are relative to the “baseline-connectivity” of the TAS.
Bayesian model selection (BMS) was performed to identify the model structure within a confined model space that most likely generated the ERPs, assuming that—because of the basic nature of the tasks—that this model structure was the same for all subjects (i.e., a fixed-effects analysis was performed; step 4 of Fig. 2)43. Nevertheless, a random-effect analysis was also performed for comparing model evidence, an approach that admits different models for different subjects who may have performed the task with different strategies43. The model evidence that accounts for both model accuracy and model complexity was used to determine posterior probabilities of the models (by normalizing the model evidence to the respective model space). Based on observed data, a large number of free parameters of the neural dynamics were estimated46. The extrinsic connectivity between cortical sources is commonly the most informative parameters that explain the changes in ERP data47. From the winning model, the estimated modulatory connectivity parameters were subjected to a 2 (SZ, HC) × 3 (TAF, TAH, TAN) repeated measures ANOVA.
