Examining resting-state functional connectivity in first-episode schizophrenia with 7T fMRI and MEG.

Schizophrenia is often characterized by dysconnections in the brain, which can be estimated via functional connectivity analyses. Commonly measured using resting-state functional magnetic resonance imaging (fMRI) in order to characterize the intrinsic or baseline function of the brain, fMRI functional connectivity has significantly contributed to the understanding of schizophrenia. However, these measures may not capture the full extent of functional connectivity abnormalities in schizophrenia as fMRI is temporally limited by the hemodynamic response. In order to extend fMRI functional connectivity findings, the complementary modality of magnetoencephalography (MEG) can be utilized to capture electrophysiological functional connectivity abnormalities in schizophrenia that are not obtainable with fMRI. Therefore, we implemented a multimodal functional connectivity analysis using resting-state 7 Tesla fMRI and MEG data in a sample of first-episode patients with schizophrenia (n = 19) and healthy controls (n = 24). fMRI and MEG data were decomposed into components reflecting resting state networks using a group spatial independent component analysis. Functional connectivity between resting-state networks was computed and group differences were observed. In fMRI, patients demonstrated hyperconnectivity between subcortical and auditory networks, as well as hypoconnectivity between interhemispheric homotopic sensorimotor network components. In MEG, patients demonstrated hypoconnectivity between sensorimotor and task positive networks in the delta frequency band. Results not only support the dysconnectivity hypothesis of schizophrenia, but also suggest the importance of jointly examining multimodal neuroimaging data as critical disorder-related information may not be detectable in a single modality alone.

Introduction

Schizophrenia is frequently described as a heterogeneous disorder of brain connectivity as it is characterized by altered functional and structural cortical network integration (Bassett et al., 2012; Friston and Frith, 1995; Stephan et al., 2009). These aberrant interactions between brain regions may be induced by abnormal physiological and neuronal processes (Bassett et al., 2012; Friston and Frith, 1995; Stephan et al., 2009; Valli et al., 2011). Resting-state functional connectivity measures are often implemented to examine these altered connections between functional networks. In the absence of a task being performed, resting-state functional connectivity measures the temporal coherence between spatially separate regions of the brain (Biswal et al., 1995; Fox and Raichle, 2007; Friston and Frith, 1995). It is believed that examination of the functional interactions between brain networks in schizophrenia can be used to better characterize the disorder as well as further elucidate the underlying mechanisms of the disorder.

Numerous functional magnetic resonance imaging (fMRI) studies examining functional connectivity in schizophrenia have been reported. These studies have indicated that connectivity is altered in schizophrenia and may be manifested via the disorder's clinical symptoms (Fornito et al., 2012; Friston and Frith, 1995; Stephan et al., 2009). Although reported fMRI functional connectivity abnormalities generally implicate alterations involving the task positive network, as well as the default mode network (DMN), connectivity abnormalities are heterogeneous across studies (Birur et al., 2017; Garrity et al., 2007; Keshavan et al., 2008; Kuhn and Gallinat, 2013; Mwansisya et al., 2017; Pettersson-Yeo et al., 2011; Whitfield-Gabrieli et al., 2009; Yu et al., 2014). Therefore, as Houck and colleagues have pointed out, although fMRI provides excellent spatial resolution, sole reliance on fMRI to explore dysconnectivity is limited by the slow temporal resolution of the hemodynamic response (Houck et al., 2017). While fMRI measures hemodynamic changes, magnetoencephalography (MEG) enables the examination of the faster neural oscillations believed to underlie the blood oxygenation level dependent (BOLD) response of fMRI (Hall et al., 2014). However, just as fMRI has limited temporal resolution, the spatial resolution of MEG is limited due to the inferences necessary for signal localization and projection onto the brain (Hall et al., 2014; Proudfoot et al., 2014). Hence, as the complementary modalities of fMRI and MEG respectively provide high spatial and temporal resolution, they can be used together to extend dysconnectivity findings in schizophrenia.

Theoretical and experimental evidence implicate post-synaptic currents as a common origin of both the fMRI and MEG signals, suggesting that while these signals have different properties, they share a common underlying electrophysiological process (Hall et al., 2014). In addition, while low frequency oscillations are thought to support interactions between distant brain regions, high frequency oscillations are thought to be restricted to localized networks (Siegel et al., 2012; von Stein and Sarnthein, 2000). Congruent with this notion, in non-human primates, a coupling between low-frequency oscillations (<20 Hz) and BOLD connectivity in a thalamo-cortical network was demonstrated during a resting state (Wang et al., 2012). Furthermore, several studies have compared the spatial maps of functional connectivity obtained with BOLD and MEG and found a degree of spatial similarity between them (Brookes et al., 2011a; Brookes et al., 2011b; de Pasquale et al., 2010). Interestingly, the best spatial match between fMRI and MEG occurs with neural oscillations in the beta frequency (Brookes et al., 2011b).

Recent multimodal functional connectivity studies have been implemented using complementary fMRI and MEG (Brookes et al., 2011b; Cetin et al., 2016; Houck et al., 2017). Only three studies have evaluated MEG resting state functional connectivity in schizophrenia (Cetin et al., 2016; Cousijn et al., 2015; Houck et al., 2017), and all three combined MEG with fMRI. One study (Cousijn et al., 2015) was limited to the evaluation of fronto-hippocampal functional connectivity. Another study (Houck et al., 2017) evaluated whole brain functional connectivity and reported a pattern of both hypo- and hyperconnectivity in frontal-DMN networks and within the frontal networks with MEG, and hypoconnectivity in temporal-occipital and frontal-occipital networks and within the occipital networks with fMRI in patients compared to controls. Using the same cohorts of chronic patients with schizophrenia and controls, the third study (Cetin et al., 2016) demonstrated improved classification of patients using combined fMRI/MEG methods, proving the importance of examining both the hemodynamic and electrophysiological functional connectivity effects when characterizing schizophrenia (Brookes et al., 2011b; Cetin et al., 2016; Houck et al., 2017).

In addition to characterizing schizophrenia using both hemodynamic and electrophysiological functional connectivity measures, it is important to examine the different stages of the disorder (e.g., premorbid, prodromal, onset/deterioration, residual) to further understand and characterize the pathophysiological mechanisms of schizophrenia. While examining all stages of the disorder are important, examination of first-episode patients with schizophrenia is especially important as studies have demonstrated that the early stages of the disorder are critical in determining the overall course and outcome of the disorder (Larsen et al., 1996). Moreover, studying first-episode patients allows for the examination of intrinsic abnormalities of the disorder rather than the confounding effects (e.g., medication, disorder chronicity on brain structure/function) introduced when investigating later stages of the illness (Hulshoff Pol and Kahn, 2008; Navari and Dazzan, 2009).

The purpose of this study is to examine resting-state functional connectivity in a sample of healthy controls and first-episode patients with schizophrenia using both 7 Tesla (7 T) fMRI and MEG. Using independent component analyses to extract resting-state networks from both fMRI and MEG, we hypothesize that patients will exhibit connectivity alterations in comparison to controls in both resting state fMRI and MEG mapping to task positive and default mode networks. Because low frequency oscillations are thought to support interactions between distant regions, we further hypothesize that functional connectivity alterations will be more similar between fMRI and MEG in the low frequency bands. Likewise, we also hypothesize that connectivity alterations examined in MEG will vary with different frequency bands (Cetin et al., 2016; Houck et al., 2017). The consistent finding that patients tend to exhibit aberrant spontaneous brain activity in the theta and delta EEG/MEG frequency bands (Ranlund et al., 2014), also suggests that MEG functional connectivity differences between patients and controls will be strongest at these low frequencies.

Materials and methods

Participants

Minimally treated first-episode patients with schizophrenia were recruited from the University of Alabama at Birmingham (UAB) emergency room, inpatient units, and various outpatient clinics. Informed consent to participate in this UAB Institutional Review Board (IRB) approved study was obtained following evaluation of competency to provide informed consent (Carpenter et al., 2000). Additionally, healthy controls matched on age, gender, parental socioeconomic status (SES), years of education, and smoking status were enrolled in the study.

Diagnoses were established with a review of medical records and evaluation by two board certified psychiatrists (ACL and NVK) and confirmed using the Diagnostic Interview for Genetic Studies (Nurnberger Jr. et al., 1994). Exclusion criteria for this study included the presence of major medical or neurological conditions, history of head trauma with loss of consciousness, substance abuse within six months of imaging (excluding nicotine), pregnancy, MRI contraindications, and use of medication altering brain function. Exclusion criteria for healthy controls also included a personal history of Axis I disorders or in first-degree relatives.

A total of 23 first-episode patients with schizophrenia were enrolled in the study. Of these, fMRI resting-state data for 19 patients remained for analysis as one subject was too large, one subject was unable to complete the scan due to scanner intolerability, one subject was excluded due to poor scan quality, and one subject did not complete a resting-state fMRI scan. MEG data for 19 patients remained for analysis as one subject was too large, one subject was lost to follow-up, one subject was hospitalized long-term, and data for one subject lacked head position sensor information. A total of 24 matched healthy controls were enrolled in the study. Of the 24 controls enrolled, fMRI resting-state data for 21 subjects remained for analysis as three subjects did not complete a resting-state fMRI scan. MEG data for 24 controls remained for MEG analysis.

Symptom severity was assessed using the Brief Psychiatric Rating Scale (BPRS) (Overall and Gorham, 1962), the Scale for the Assessment of Negative Symptoms (SANS) (Andreasen, 1984a), and the Scale for the Assessment of Positive Symptoms (SAPS) (Andreasen, 1984b). The Repeatable Battery for the Assessment of Neuropsychological Status (RBANS) (Randolph et al., 1998) was used to assess cognitive function for healthy controls and patients.

Functional MRI analysis

fMRI scanning parameters

All scans were performed at the Auburn University MRI Research Center on a whole-body 7 T Siemens MAGNETOM MRI scanner (Siemens Healthineers, Erlangen, Germany), with a 32-channel head coil. High-resolution structural scans were acquired using a 3-dimensional T1-weighted magnetization-prepared rapid acquisition gradient-echo sequence (MPRAGE; 256 slices, repetition time/echo time/inversion time [TR/TE/TI] = 2000/2.89/1050 ms, 7° flip angle, 190 mm field of view, 0.7 mm isotropic voxels, base resolution = 256, sagittal acquisition, GRAPPA acceleration factor = 2). Resting-state fMRI scans were acquired using a gradient recalled echo-planar imaging sequence (TR/TE = 3000/28 ms, 70° flip angle, 200 mm field of view, 37 slices, 0.85 mm × 0.85 mm × 1.8 mm voxels, iPAT GRAPPA acceleration factor = 3, base resolution = 234, interleaved acquisition, A > P phase encode direction, 1 ms echo spacing, gap = 1.08 mm). Six-minute resting-state scans consisted of 120 volumes. Field maps were acquired using a gradient-recalled echo (GRE) sequence (TR/TE1/TE2 = 400/4.92/7.38 ms, 60° flip angle, 200 mm field of view, 36 slices, 3.1 mm × 3.1 mm × 3.0 mm voxels, base resolution = 64, interleaved acquisition, gap = 0.75 mm).

fMRI preprocessing

Data preprocessing of resting-state scans was performed with SPM12 (Wellcome Trust Centre for Neuroimaging, London, UK; http://www.fil.ion.ucl.ac.uk/spm/) and the “CONN” Connectivity Toolbox (Whitfield-Gabrieli and Nieto-Castanon, 2012) using standard preprocessing steps. More specifically, preprocessing steps included realignment and unwarping using phase maps, slice timing correction, coregistration to anatomical space using the T1-weighted MPRAGE structural image and using the first functional image as reference, normalization to Montreal Neurological Institute (MNI) space using the first functional volume as reference, artifact detection, and Gaussian smoothing [5 mm full-width half-maximum (FWHM)]. Resulting preprocessed functional resting-state data maintained a spatial resolution of 0.855 mm × 0.855 mm × 2.88 mm.

fMRI group independent component analysis

The Group ICA of fMRI Toolbox (GIFT; http://mialab.mrn.org/software/gift) was used to perform group-level spatial independent component analysis (ICA) (Calhoun et al., 2001). A subject-specific data reduction principal component analysis (PCA) first reduced the data to 75 principal components followed by a group data reduction retaining 50 principal components (Allen et al., 2014). In order to carry out PCA in a memory efficient-manner, PCA was carried out using the expectation maximization algorithm (Allen et al., 2014; Roweis, 1998). The infomax algorithm (Allen et al., 2014; Bell and Sejnowski, 1995) was then utilized for group-level spatial ICA to generate 50 spatially independent components (ICs) (See Supplementary Fig. 1 for estimation of the number of ICs using minimum description length criteria). In order to measure component stability/quality, the infomax ICA algorithm was repeated 20 times in ICASSO and the most representative run was used in subsequent steps (Allen et al., 2014; Damaraju et al., 2014; Himberg et al., 2004). Subject-specific spatial maps and time courses were estimated using GICA back-reconstruction (Erhardt et al., 2011).

Supplementary Fig. 1

Estimate of the number of independent components. The mean number of independent components for fMRI (left) and MEG (right) was estimated to be 32 and 49 for all subjects, respectively, using minimum description length criteria. For fMRI, the minimum and maximum number of components estimated using the minimum description length criteria was 17 and 55, respectively. For MEG, the minimum and maximum number of components estimated using the minimum description length criteria was 12 and 110, respectively.

fMRI component identification

In accordance with previous studies (Allen et al., 2011; Cetin et al., 2016; Houck et al., 2017; Lottman et al., 2017), ICs were classified as resting-state networks (RSNs), rather than artifacts, based on spectral properties and functional relevance. The power spectra for each component were characterized by dynamic range and the ratio of low frequency to high frequency power. Dynamic range can be defined as the difference between maximum power and minimum power at frequencies to the right of the power spectra peak (Allen et al., 2011; Cetin et al., 2016; Houck et al., 2017). The ratio of low frequency to high frequency power can be defined as the ratio of integral spectral power below 0.1 Hz to integral spectral power between 0.15 and 0.167 Hz (Allen et al., 2011; Cetin et al., 2016; Houck et al., 2017). Components were also evaluated for functional relevance via visual inspection. Components were inspected and classified as non-artifactual or RSNs based on the following criteria: (1) peak activations should occur primarily in grey matter and correspond anatomically to known functional brain networks, (2) low spatial overlap with vascular, ventricular, motion and susceptibility artifacts, (3) time courses predominantly characterized by low frequency fluctuations (Allen et al., 2011; Cetin et al., 2016; Cordes et al., 2000; Houck et al., 2017; Lottman et al., 2017). Of the 50 components extracted, 30 were identified as RSNs (Fig. 1).

Fig. 1

Composite fMRI and MEG spatial maps. Composite maps of the 30 fMRI (left column) and 33 MEG (right column) resting-state independent component spatial maps. Components are categorized into auditory, sensorimotor, visual, task positive, default-mode, and subcortical networks. Colors in each map represent a different component and are indicated with the component number. The number of components in each network is indicated behind the network name. Peak activations of individual components can be found in Table S1 (fMRI) and Table S2 (MEG) in Supplementary Material.

MEG analysis

MEG scanning parameters

Resting-state MEG recordings were collected in a shielded room below ground level using a 148-channel whole head magnetometer (MagnesTM 2500 WH, 4-D Neuroimaging) at the University of Alabama at Birmingham. Recordings were sampled at a rate of 1000 Hz, with a bandpass filter of 0.1 to 200 Hz. Fiducial head-sensor coils were used to ensure the subject's head did not move during a scan. Additional head-sensor data were collected using a 3D digitizer and later used for registration between each subject's MEG and anatomical MRI. During the 5-min resting-state scan, subjects were placed in the supine position with eyes closed. Visual inspection of the real-time display showed no MEG correlates of sleep. Subjects were also clearly responsive immediately before and after the scanning period as they all instantly opened their eyes upon instruction.

MEG preprocessing

MEG data were analyzed using Brainstorm (Tadel et al., 2011). Preprocessing of the raw data included removal of line noise data using a notch filter (60, 120, 180, 240, 300 Hz) and high-pass filtering at 0.3 Hz. Cardiac and motion artifact was removed from the signal using an independent component analysis where the signal was decomposed into 20 independent components using JADE (joint approximate diagonalization of eigenmatrices) (Cardoso and Souloumiac, 1993). Additional artifacts were removed using a signal-space projection method.

Resting-state MEG data was then registered to individual subject high resolution structural MRI scans. Prior to MEG-MRI registration, cortical surfaces and subcortical structures from each structural MRI were extracted using the FreeSurfer image analysis suite (http://surfer.nmr.mgh.harvard.edu/) and subsequently inspected using FreeSurfer's QA Tools. Data identified as unsatisfactory during quality assurance were manually corrected and reprocessed. A pseudo-individual anatomy for each subject missing structural MRI scans (n = 5) was created by warping the default ICBM 152 template anatomy in Brainstorm to the digitized head points. A 15,000-vertex cortical surface source model was used for head modeling.

Forward volume head modeling of the magnetic field was defined using an overlapping-sphere method where voxels were placed on a 5 mm3 grid of the brain (Huang et al., 1999). Inverse source modeling was performed using a linearly-constrained minimum variance (LCMV) beamformer. Beamformer projection was performed using an unconstrained dipole orientation for each subject separately. Data covariance matrices generated for each subject were regularized using a median eigenvalue approach where the tail of the eigenvalues spectrum is flattened to the median value.

MEG Hilbert envelope computation

The Hilbert transform was applied to the resting-state time courses in which the signals were first filtered into delta (1–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), and gamma (30–50 Hz) frequency bands using a band-pass filter followed by application of the Hilbert transform. The absolute value of the transformed signals was subsequently computed in order to provide a measure of the envelope of the signal. The Hilbert envelope values were then standardized using a 1/frequency compensation where values at each frequency bin were multiplied by the frequency value. Following spectrum normalization, the Hilbert envelope was downsampled to 1 Hz in accordance with previous studies (Brookes et al., 2011b; Cetin et al., 2016; Houck et al., 2017). Therefore, scans consisted of 301 volumes for each frequency band. In order to enable comparison between subjects and groups, downsampled source envelope data for each subject were spatially interpolated from the individual subject head model to the MNI ICBM 152 default MRI template in Brainstorm and spatially smoothed (5 mm FWHM).

MEG Group independent component analysis

Group-level spatial ICA was performed on the MEG data in the GIFT Toolbox where each frequency band – delta (1–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), and gamma (30–50 Hz) – was input as a different session for each individual. Similar to fMRI group ICA analysis, MEG data were first reduced to 75 principal components using a subject-specific data reduction PCA and then reduced to 50 principal components via group data reduction (Allen et al., 2014). PCA was carried out using the expectation maximization algorithm (Allen et al., 2014; Roweis, 1998). The infomax algorithm (Allen et al., 2014; Bell and Sejnowski, 1995) was then utilized for group-level spatial ICA to generate 50 spatially independent components (ICs) (See Supplementary Fig. 1 for estimation of the number of ICs using minimum description length criteria). Component stability/quality was measured by repeating the infomax ICA algorithm 20 times in ICASSO (Allen et al., 2014; Damaraju et al., 2014; Himberg et al., 2004). Subject-specific spatial maps and time courses were estimated using GICA back-reconstruction (Erhardt et al., 2011).

MEG component identification

MEG ICs were classified as RSNs, rather than artifacts, following the same method described in the previous section: fMRI Component Identification. Of the 50 components extracted, 33 were identified as RSNs (Fig. 1).

Functional connectivity analysis

Functional connectivity analysis

Functional connectivity analysis for fMRI and MEG was performed using the spatial maps and corresponding time courses extracted during the ICA analysis. For the fMRI data, framewise displacement was regressed from the subject-specific RSN time courses prior to connectivity analyses. Framewise displacement was computed as the absolute frame-to-frame displacement of the brain from six realignment parameters using a radius of 50 mm to convert angle rotations to displacements (Lottman et al., 2017; Power et al., 2012). Therefore, individual displacement values for each volume of the time course (i.e., 120 frame-to-frame displacement values) resulted. Framewise displacement was regressed from the subject-specific RSN time courses prior to functional connectivity analysis. It is important to note that mean framewise displacement did not differ between groups [HC: 0.19±0.09 mm (range: 0.07–0.44 mm); SZ: 0.21±0.14 mm (range: 0.04–0.47 mm); t = −0.78, p = .44]. No participants exceeded mean framewise displacement of 0.5 mm. The number of head motion outliers (defined as FD >0.5 mm) did not differ between groups (HC: 5.67±6.51%, SZ: 6.79+/− 7.49%, t = −0.51 p = .62). Subject-specific fMRI and MEG RSN time courses were detrended and despiked using 3dDespike (http://afni.nimh.nih.gov/afni) and then filtered with a fifth-order Butterworth low-pass filter (0.15 Hz cutoff) prior to connectivity computation (Allen et al., 2014; Allen et al., 2011; Damaraju et al., 2014; Lottman et al., 2017).

Functional connectivity was calculated as the pairwise Pearson correlation between RSN time courses, resulting in individual correlation matrices for each subject's fMRI data and MEG frequency bands. Time courses were extracted from the independent components obtained by ICA. Correlations were based on the average time course across voxels within each independent component. Individual correlation matrices were averaged across subjects in each group in order to generate group-level connectivity matrices. Within- and between-group differences in functional connectivity were evaluated on subject-level correlation matrices via paired and two-sample univariate t-tests, respectively.

In exploratory post hoc analyses, the relationship between functional connectivity and symptom scores including, BPRS, SAPS, and SANS was explored using Spearman correlations in order to examine the potential impact of symptoms on functional connectivity.

Results

Demographics

No significant differences in age, gender, parental SES, years of education, and smoking status (packs per day) were exhibited between patients and controls (Table 1). Of the patients enrolled in the study, most were treated with risperidone, one was unmedicated, two were treated with aripiprazole, one with clozapine, and one with a combination of ziprasidone and clozapine.

Table 1

Demographics and clinical assessments for MEG and fMRI connectivitya.

	HC (n = 24)	SZ (n = 21)b	t/χ2	p-value
Age (years)	24.13 ± 5.02	23.52 ± 4.64	0.415	0.680
Gender (male/female)	17/7	15/6	0.002	0.965
Parental SESc	3.46 ± 3.32	4.25 ± 4.38	6.272	0.617
Educationd	3.08 ± 0.58	2.85 ± 0.59	1.747	0.418
Smoking (packs per day)	0.01 ± 0.06	0.07 ± 0.15	−1.607	0.121
Duration of Illness (months)e	–	18.58 ± 26.68	–	–
BPRSf (n = 19)
Total score	–	32.26 ± 9.82
Positive symptom subscale	–	5.05 ± 2.95
Negative symptom subscale	–	5.74 ± 2.21
SANSg (n = 17)
Total composite scoreh	–	21.56 ± 17.81
Global summary score	–	6.47 ± 4.96
SAPSi (n = 17)
Total composite score	–	7.12 ± 11.77
Global summary score	–	2.47 ± 3.48
RBANSj
Total index	93.43 ± 8.66	74.35 ± 15.31	4.579	<0.001
Immediate memory	99.67 ± 12.11	79.65 ± 14.81	4.587	<0.001
Visuospatial	86.24 ± 14.98	79.94 ± 17.15	1.208	0.235
Language	103.71 ± 10.72	82.82 ± 11.17	5.862	<0.001
Attention	96.14 ± 13.54	76.47 ± 20.66	3.532	0.001
Delayed Memory	92.05 ± 8.27	79.88 ± 19.60	2.393	0.026

Abbreviations: HC, healthy control; SZ, schizophrenia; SES, socioeconomic status; Y, yes; N, no; BPRS, Brief Psychiatric Rating Scale; SANS, Scale for the Assessment of Negative Symptoms; SAPS, Scale for the Assessment of Positive Symptoms; RBANS, Repeated Battery for the Assessment of Neuropsychological Status.

Mean ± SD unless otherwise indicated.

SZ values for SES, education, and smoking (packs per day) calculated with n = 20 due to data missing for one patient.

SES ranks reported from Diagnostic Interview for Genetic Studies scale (1–18); high rank (lower numerical value) corresponds to high socioeconomic status.

Years of education reported from Diagnostic Interview for Genetic Studies scale.

Duration of illness (months) calculated with n = 19 due to missing data for 2 patients.

BPRS reported on 1–7 scale; positive (conceptual disorganization, hallucinatory behavior, and unusual thought content); negative (emotional withdrawal, motor retardation, and blunted affect).

SANS includes five subscales: affective flattening or blunting, alogia, avolition-apathy, anhedonia-asociality, and attention.

SANS average total composite score calculated with data from 16 SZ due to missing data from one SZ.

SAPS includes four subscales: hallucinations, delusions, bizarre behavior, and positive formal thought disorder.

RBANS data missing from 3 HC (n = 21) and 3 SZ (n = 17).

As illustrated in Table 1, patients demonstrated significantly lower immediate memory, language, attention, delayed memory, and total RBANS scores in comparison to controls. Patient symptoms scores can also be found in Table 1. No significant differences in mean framewise displacement (t = −0.78, p = .44) between controls and patients were found.

Independent components

fMRI

Of the 50 extracted ICs, 30 were identified as RSNs as opposed to artifacts based on the criteria described in the methods (Fig. 1). Identified RSNs were organized into six different networks including subcortical (1 RSN), auditory (1 RSN), sensorimotor (7 RSNs), visual (5 RSNs), task positive (9 RSNs), and default-mode (7 RSNs) networks (Fig. 1; Supplementary Table 1).

MEG

Of the 50 extracted MEG ICs, 33 were identified as RSNs (Fig. 1) and organized into auditory (3 RSNs), sensorimotor (10 RSNs), visual (6 RSNs), task positive (9 RSNs), and default-mode (5 RSNs) networks (Fig. 1; Supplementary Table 2).

fMRI functional connectivity

Mean functional connectivity for controls and patients are illustrated in the top row of Fig. 2. Patients demonstrated significantly higher connectivity (p < .001) between subcortical (IC 49) and auditory (IC 18) RSNs compared to controls (Fig. 2, Fig. 3). However, patients also exhibited significantly lower (p < .001) connectivity between two sensorimotor RSNs (i.e., IC1 and IC2) when compared to controls (Fig. 2, Fig. 3). Connectivity results are corrected for multiple comparisons using the false discovery rate (FDR) and reported at pFDR < 0.05.

Fig. 2

Resting-state fMRI and MEG functional connectivity. Functional connectivity for resting-state fMRI (top) and MEG averaged across frequency bands (bottom). Connectivity for controls (left column), first-episode patients (middle column), and the significant difference between groups (right column) are shown. FDR-corrected significant differences (p < 0.05) are indicated with • and uncorrected significant differences (p < 0.05) are indicated with *.

Fig. 3

fMRI connectogram and symptom correlations. (A) fMRI connectogram depicting significant (p < 0.05) connectivity differences between patients and controls. Warm colors represent decreased and cool colors represent increased connectivity in patients. (B) Correlation between significant subcortical-auditory connectivity values and patients BPRS positive scores.

Exploratory post hoc analyses revealed a significant relationship (rho = 0.59; p = .013, confidence interval: 0.1175 to 0.8533), using Spearman correlations over 1000 bootstrap, between patient BPRS positive symptom score and functional connectivity between subcortical and auditory RSNs (Fig. 3). This indicates that as BPRS positive symptoms increase, connectivity between subcortical and auditory networks in patients also increases. However, given the number of correlations conducted, including with the BPRS negative subscale, the SAPS and SANS total scores (all p > .05), this correlation would not survive multiple comparison correction.

MEG functional connectivity

Mean functional connectivity for the average of all frequency bands is shown for controls and patients in the bottom row of Fig. 2. In the delta frequency band, patients demonstrated significantly lower connectivity (p < .001) between sensorimotor (IC 5) and task positive (IC 34) RSNs (Fig. 4, Fig. 5). Patients also exhibited significantly lower connectivity (p < .001) between sensorimotor (IC 25) and task positive (IC 34) RSNs in the delta band (Fig. 4, Fig. 5). Connectivity results are corrected for multiple comparisons and reported at pFDR < 0.05. No other significant group differences corrected for multiple comparisons were exhibited in the remaining frequency bands. Frequency band-specific functional connectivity for each group and group differences (puncorrected < 0.05) are illustrated in Supplementary Fig. 2, Supplementary Fig. 3, Supplementary Fig. 4, Supplementary Fig. 5.

Fig. 4

Resting-state delta frequency band MEG functional connectivity. Functional connectivity for resting-state MEG delta band. Connectivity for controls (left column), first-episode patients (middle column), and the significant difference between groups (right column) are shown. FDR-corrected significant differences (p < 0.05) are indicated with • and uncorrected significant differences (p < 0.05) are indicated with *.

Fig. 5

MEG delta frequency band connectogram. MEG delta band connectogram depicting significant (p < 0.05) connectivity differences between patients and controls. Warm colors represent decreased and cool colors represent increased connectivity in patients.

Supplementary Fig. 2

Resting-state theta frequency band MEG functional connectivity. Functional connectivity for resting-state MEG theta band. Connectivity for controls (left column), first-episode patients (middle column), and the significant difference between groups (right column) are shown. FDR-corrected significant differences (p < 0.05) are indicated with • and uncorrected significant differences (p < 0.05) are indicated with *.

Supplementary Fig. 3

Resting-state alpha frequency band MEG functional connectivity. Functional connectivity for resting-state MEG alpha band. Connectivity for controls (left column), first-episode patients (middle column), and the significant difference between groups (right column) are shown. FDR-corrected significant differences (p < 0.05) are indicated with • and uncorrected significant differences (p < 0.05) are indicated with *.

Supplementary Fig. 4

Resting-state beta frequency band MEG functional connectivity. Functional connectivity for resting-state MEG beta band. Connectivity for controls (left column), first-episode patients (middle column), and the significant difference between groups (right column) are shown. FDR-corrected significant differences (p < 0.05) are indicated with • and uncorrected significant differences (p < 0.05) are indicated with *.

Supplementary Fig. 5

Resting-state gamma frequency band MEG functional connectivity. Functional connectivity for resting-state MEG gamma band. Connectivity for controls (left column), first-episode patients (middle column), and the significant difference between groups (right column) are shown. FDR-corrected significant differences (p < 0.05) are indicated with • and uncorrected significant differences (p < 0.05) are indicated with *.

Exploratory post hoc analyses did not indicate a significant relationship between patient symptom scores (e.g., BPRS, SANS, SAPS) and significant functional connectivity in the delta frequency band.

Spatial maps

Spatial maps for each modality were evaluated as a one-sample t-test of subject-specific back-reconstructed spatial maps across groups (Fig. 1). In accordance with previous studies (Houck et al., 2017), spatial correlation of the rsfMRI and MEG group-level (rather than subject back-reconstructed) spatial map components was evaluated. Spatial correlations between modality components and >0.50 were exhibited in sensorimotor, visual, and default-mode networks and illustrated in Supplementary Fig. 6, Supplementary Fig. 7, Supplementary Fig. 8.

Supplementary Fig. 6

Sensorimotor networks spatial correlations. Spatial correlation (r) between resting-state fMRI and MEG components in the sensorimotor network. Only correlations of r > 0.5 are illustrated.

Supplementary Fig. 7

Visual networks spatial correlations. Spatial correlation (r) between resting-state fMRI and MEG components in the sensorimotor network. Only correlations of r > 0.5 are illustrated.

Supplementary Fig. 8

Default-mode networks spatial correlations. Spatial correlation (r) between resting-state fMRI and MEG components in the default-mode network. Only correlations of r > 0.5 are illustrated.

Discussion

To our knowledge this is the first study examining resting-state functional connectivity in first-episode patients with schizophrenia using both 7 T fMRI and MEG. We describe fMRI connectivity alterations in subcortical, auditory, and sensorimotor networks and MEG connectivity alterations in sensorimotor and task positive networks in patients. As hypothesized, group differences in MEG functional connectivity were identified in a low frequency band. The identified connectivity alterations are all regions implicated in the disorder of schizophrenia and lend further support to the dysconnectivity hypothesis of schizophrenia (Friston and Frith, 1995; Pettersson-Yeo et al., 2011; Stephan et al., 2009).

Functional connectivity analysis of fMRI identified hyperconnectivity between the subcortical network (i.e., caudate) and auditory network (i.e., superior temporal gyrus) in first-episode patients (Fig. 2). Exploratory post hoc exploration of this subcortical-auditory relationship indicates a significantly positive correlation (r = 0.59; p = .01) between connectivity strength and BPRS positive scores in patients (Fig. 4B). These findings may be indicative of the link both regions have been shown to exhibit with positive symptoms (e.g., hallucinations) in schizophrenia. A relationship between hallucinations and structural and functional alterations in the superior temporal gyrus has been reported in a number of studies (Allen et al., 2008). In addition, these results are also indicative of the importance of the caudate in schizophrenia and its role in treatment response as it is highly innervated by dopaminergic neurons (Hutcheson et al., 2014). Future studies examining unmedicated first-episode patients would provide greater insight into the relationship between the caudate, temporal lobe, and positive symptoms.

Decreased fMRI functional connectivity between two sensorimotor network components was also demonstrated in patients (Fig. 2). Specifically, this hypoconnectivity demonstrated in patients is exhibited between the sensorimotor network components comprised of the right precentral gyrus (IC 1) and left postcentral gyrus (IC2). Further examination of these sensorimotor components reveals they share interhemispheric homotopic or bilaterally symmetric regions as both IC1 and IC2 comprise regions of the precentral and postcentral gyrus. This altered sensorimotor connectivity is consistent with previous studies demonstrating decreased interhemispheric homotopic connectivity in first-episode patients with schizophrenia (Guo et al., 2014b; Li et al., 2015). Decreased connectivity between homotopic motor regions are thought to be related to cognitive dysfunction associated with motor function in schizophrenia (Guo et al., 2014b). Moreover, decreased interhemispheric connectivity in the precentral gyrus has also been shown to be negatively related to positive (Guo et al., 2014b) and negative symptomatology (Li et al., 2015).

Group differences in the MEG delta band indicated hypoconnectivity between sensorimotor and task positive network connections among patients. Specifically, the precuneus (IC 5) and superior parietal lobule (IC 25) connections with the insula (IC34) are altered in the delta frequency band in patients. Alterations in low frequency neural oscillations at rest have been implicated in previous studies in schizophrenia (Moran and Hong, 2011; Uhlhaas and Singer, 2010). Supporting our results, previous studies have demonstrated altered delta band activity in frontal, temporal and parietal areas in schizophrenia (Alamian et al., 2017). In addition, our results are consistent with the literature indicating slow oscillations are more prone to travel longer distances and connect relatively distant regions of the brain (Moran and Hong, 2011). In fact, previous studies have suggested functional connectivity is affected by anatomical distance in schizophrenia (Guo et al., 2014a; Guo et al., 2017). Accordingly, our results indicate long range functional connections between the task positive and parietal sensorimotor networks are altered in schizophrenia. These results may be reflective of positive symptoms due to faulty sensory information processing characteristic of schizophrenia (Alamian et al., 2017; Wylie and Tregellas, 2010).

Visual inspection of the patterns of group differences (puncorrected < 0.05) indicate patients tend to exhibit instances of hypoconnectivity across all networks with the exception of hyperconnectivity across subcortical network connections. However, patients generally demonstrated instances of hypoconnectivity within and between sensorimotor network connections in MEG in the overall mean and individual frequency bands. Although future studies would be necessary to determine if these differences survive multiple comparison correction in a larger sample size, altered sensorimotor connectivity has been reported in schizophrenia (Kaufmann et al., 2015).

Our results using group independent component analysis as well as functional connectivity measures of both resting-state fMRI and MEG indicate the importance of using these complementary imaging modalities to capture a more comprehensive picture of disorder-related abnormalities. In particular, using both fMRI and MEG allows us to investigate both spatial and temporal characteristics of resting-state network functional connectivity that would have been missed by examining a single modality alone (Cetin et al., 2016; Houck et al., 2017). In fact, a previous study of chronic schizophrenia patients indicated that the combination of fMRI and MEG functional connectivity data more accurately classifies patients than either modality alone (Cetin et al., 2016). Furthermore, both modalities provide common and unique findings characterizing schizophrenia. As illustrated in Supplementary Fig. 6, Supplementary Fig. 7, Supplementary Fig. 8, similar and spatially overlapping RSNs were generated from separate group independent component analyses of fMRI and MEG further confirming the electrophysiological underpinnings of BOLD fMRI. While functional connectivity of fMRI and MEG both indicated connectivity alterations involving the sensorimotor networks (although the spatial overlap between fMRI and MEG sensorimotor networks alterations is not clear), unique findings of subcortical-auditory and delta band sensorimotor-task positive connectivity were exhibited in fMRI and MEG, respectively.

When interpreting our findings, several strengths and limitations should be considered. To minimize variance in the data, subjects were matched on factors including age, gender, parental SES, and smoking status. As it has been shown that the early stages of schizophrenia are critical in determining the overall course and outcome of the disorder (Larsen et al., 1996), a sample of medicated first-episode patients with schizophrenia were enrolled in the study. However, our sample size was limited and replication studies with larger sample size will be needed. Although patients were examined early in the course of the disorder and minimally treated, antipsychotic medication effects on brain function cannot be differentiated from intrinsic characteristics of the disorder. Patients in this study were receiving treatment with antipsychotic medication. Our previous fMRI studies have shown that these drugs can modify functional connectivity patterns (Hadley et al., 2014; Kraguljac et al., 2016). Further studies enrolling unmedicated or medication naïve patients are needed. While rigorous preprocessing steps were implemented on both the fMRI and MEG data, heart rate and breathing were not directly controlled for during data acquisition. Although our preprocessing steps and ICA analysis attempted to control for these potential sources of noise, our results may still be impacted by artifacts due to motion or physiological processes. Subjects self-reported that they did not fall asleep during the MEG scan, but it is still possible that some may have briefly fallen asleep and impacted our results. The number of significant functional connectivity group differences reported are much smaller than typically reported in the literature. We believe our findings may be impacted by the modest sample size. In addition, we examine functional connectivity for multiple regions independently comprising a network rather than using a single network template. Therefore, as multiple networks are examined in a functional connectivity analysis, our results may be impacted by the number of comparisons being made between functional networks (e.g., over 450 and 500 comparisons for fMRI and MEG, respectively). However, we believe that the compatibility of our results with the schizophrenia literature indicate the importance of our multiple comparison corrected and uncorrected functional connectivity findings. Finally, as previous chronic schizophrenia studies have explored (Cetin et al., 2016), future studies would benefit from utilizing classification analyses to quantitatively determine if examining functional connectivity using both fMRI and MEG provides better classification of first-episode patients with schizophrenia than using either modality alone.

Conclusion

In conclusion, our results support the dysconnectivity hypothesis of schizophrenia as patients demonstrated altered connectivity between fMRI subcortical and auditory networks, as well as fMRI and MEG altered intra- and inter- sensorimotor connectivity in both fMRI and MEG. These findings illustrate the strength in combining complementary modalities like fMRI and MEG to examine resting-state brain function in schizophrenia. More specifically, this study further highlights the importance of examining multimodal neuroimaging data when studying schizophrenia as critical information relating to different aspects of the brain (e.g., electrophysiological, hemodynamic, structural, etc.) may be missing in a single modality.

The following are the supplementary data related to this article.

Supplementary material 1

Supplementary material 2

Declaration of Competing Interest

Kristin Lottman, Timothy Gawne, Nina Kraguljac, Jeffrey Killen, Meredith Reid, and Adrienne Lahti declare no potential conflicts of interest.