Heterogeneity in clinical presentation and mechanisms of major depressive disorder (MDD) likely contributes to limited responses to current treatments in many patients. Identifying biologically meaningful phenotypes would constitute a major step towards the development of personalized treatment approaches. Brain-activity-based phenotyping offers a promising route toward this goal. In particular, brain oscillations—rhythmic patterns of neural activity that support information processing—have been implicated in depression, but have not previously been used to define biological phenotypes of the disorder. Yet, no studies have used brain oscillations to identify biological depression phenotypes. Here we report data-driven identification of oscillation phenotypes for MDD. We conducted a cross-sectional study and collected resting-state magnetoencephalography (MEG), structural MRI, and clinical symptom data from 263 patients with MDD and 75 healthy controls. We assessed oscillation-based functional connectivity from source-reconstructed MEG data with two coupling-mode measures and computed their low-dimensional brain-symptom associations to obtain latent components. Using clustering methods on these components, we identified five depression phenotypes which were characterized by distinct spectral and spatial patterns and differentiated clinically unique symptom profiles. These findings suggest that MEG-based oscillatory connectivity captures clinically relevant heterogeneity in MDD and provides candidate mechanistic phenotypes for future validation and treatment-stratification studies.

Symptom measures

All participants were required to complete 10 self-report symptom scales to assess various aspects of mental health and functioning. These scales include Patient Health Questionnaire (PHQ-9) measuring depressive symptoms (range: 0–27), Quick Inventory of Depressive Symptomatology 16-item (QIDS) evaluating depressive symptom severity (range: 0–27), Generalized Anxiety Disorder 7-item (GAD) Scale assessing anxiety symptoms (range: 0–21), Ruminative Responses Scale Short-version (RRS) evaluating ruminative thought patterns (range: 8–32), Sheehan Disability Scale (SDS) evaluating functional impairment (range: 0–30), Brief Experiential Avoidance Questionnaire (BEAQ) measuring experiential avoidance (range: 15–90), PTSD Checklist for DSM-5 (PCL) assessing post-traumatic stress disorder symptoms (range: 0–80), Alcohol, Smoking, and Substance Involvement Screening Test Lite (ASSIST) assessing substance use involvement (range: 0–20), Positive Valence Systems Scale (PVSS) examining positive affect and reward processing (range: 21–189), and WHO-5 Well-Being Index (WHO) measuring subjective well-being (range: 0–25).

Neuroimaging data acquisition

Fifteen minutes eyes-open resting-state MEG data were recorded with a 306-channel MEG system (TRIUX or TRIUXneo, MEGIN Oy, Helsinki, Finland; 204 planar gradiometers and 102 magnetometers). Bipolar horizontal and vertical electrooculography (EOG) and electrocardiography (ECG) were recorded for detection of eye movements and cardiac artefacts. Participants were instructed to sit in a dimly lit room and to focus on a fixation cross. T1-weighted anatomical MRI scans were obtained with a 3-tesla whole-body MRI scanner (Magnetom Skyra, Siemens, Erlangen, Germany) at a resolution of 0.8×0.8×0.8 mm, repetition time of 2530 ms, and echo time of 3.42 ms.

MEG data preprocessing and source modeling

Temporal signal space separation (tSSS) in the Maxfilter software (Elekta Neuromag) was used to suppress extracranial noise from MEG sensors, interpolate bad channels, and compensate for head motions. 50-Hz line noise and its harmonics were removed with a Finite Impulse Response (FIR) notch filter. Independent components analysis was used to remove ocular, heartbeat, and muscle artifacts. We used the FreeSurfer software (https://surfer.nmr.mgh.harvard.edu/) for volumetric segmentation of MRI data, surface reconstruction, flattening, and cortical parcellation. Source reconstruction was performed with minimum norm estimation (MNE) with the MNE software (https://mne.tools/stable/index.html). A surface-based source space with 5-mm spacing and 1-layer (inner skull) symmetric boundary element method (BEM) was used in computing the forward operator. Noise covariance matrices (NCM) were obtained from preprocessed data filtered to 151–249 Hz. We then estimated vertex fidelity by applying forward and inverse operators to complex white-noise time series and computing the correlation between original and forward-inverse-modeled time series. The obtained fidelity-weighted inverse operators were used for collapsing vertex time series into the 200 parcels of the Schaefer atlas. Broadband parcel time series were then filtered into narrowband time series with 32 Morlet wavelets with center frequencies spanning from 2.1 to 59.3 Hz in log-linear space.

Phase synchrony and amplitude correlations

Phase synchrony (PS) and amplitude correlations (AC) are two intrinsic modes of oscillation-based functional connectivity. We applied the weighted phase lag index to assess PS and the orthogonalized correlation coefficient to measure AC, which are maximally insensitive to false-positive interactions due to source leakage.