Can polygenic-informed EEG biomarkers predict differential antidepressant treatment response? An EEG stratication marker for rTMS and sertraline

The treatment of major depressive disorder (MDD) is hampered by low chances of treatment response in each treatment step, which is partly due to a lack of rmly established outcome-predictive biomarkers. Here, we hypothesize that polygenic-informed EEG biomarkers may help predict differential antidepressant treatment response. Using a polygenic-informed electroencephalography (EEG) data-driven, data-reduction approach, we identify a functional brain network that is sex-specically associated with polygenic risk for MDD in psychiatric patients (N=1,123). Subsequently, we demonstrate the utility of this network in predicting response to transcranial magnetic stimulation (TMS) and antidepressant medication in two independent datasets (N=196 and N=1,008). A simulation aimed at stratifying patients to TMS, sertraline or escitalopram/venlafaxine based on only this EEG component yields up to >30% improved remission rates. Overall, our ndings highlight the power and utility of a combined polygenic and neurophysiological approach in the search for clinically-relevant biomarkers in psychiatric disorders.


Introduction
Major depressive disorder (MDD) is a common psychiatric disorder with a complex etiology that is generally explained from a biopsychosocial model, in which multiple biological, psychological, and social factors are all considered important contributors. 1,2 Furthermore, genetic risk factors of MDD overlap with other psychiatric disorders. 3 It is assumed that this multifactorial model for MDD underlies its heterogeneous symptomatology and variable treatment e cacy, 4,5 as well as the higher prevalence and different presentation of MDD in women compared to men. 6,7 In line with the biological heterogeneity of MDD that in turn may be related to treatment outcome, pharmacogenomic studies have focused on genetic biomarkers of antidepressant treatment response in MDD. Genome wide association studies (GWASs) have identi ed several (common) genetic variants that are associated with antidepressant e cacy, but clinically-relevant and converging evidence has remained elusive. [8][9][10][11][12][13][14][15] For repetitive transcranial magnetic stimulation (rTMS) responsiveness, to our knowledge, only one GWAS at relatively limited power is available. 16 Thus, antidepressant treatment outcome is likely a complex trait and explained by several loci of small effect, 17 with recent evidence indeed suggesting that antidepressant response is polygenic. 18 Consequently, a polygenic instead of single gene or locus approach, by calculation of the individual's polygenic risk score (PRS) seems valuable to associate genetic risk with treatment (non)response. 19 However, despite recent evidence showing the power of PRS of MDD for a range of MDD-related phenotypes, 20 at present evidence for the out-of-sample value of polygenic risk approaches in the prediction of treatment outcome is limited. 18,[21][22][23][24] A proposed strategy to effectively predict therapeutic outcomes for clinically prognostic purposes, is to integrate PRS with other predictors, such as neuroimaging and clinical characteristics. 25 Electroencephalography (EEG) is a non-invasive neuroimaging technique that is used to quantitatively analyze oscillatory brain activity of neurons with high temporal resolution. 26 Several EEG patterns are heritable, in particular characteristics within the alpha frequency band and EEG power across the power spectrum. [27][28][29][30] Some studies have also demonstrated heritability for functional connectivity and ('small world') network organization parameters, which have been linked to pathological states of the brain. [31][32][33] EEG biomarker research for treatment prediction within MDD has shown that certain EEG patterns or abnormalities are differentially associated with drug-speci c or drug-class speci c antidepressant treatment effects [34][35][36] , as well as rTMS outcome. [37][38][39][40] Such studies have also demonstrated qualitative sex differences in topographic distribution of EEG activity and sex-speci c predictors of treatment response of alpha asymmetry, 35 EEG connectivity, 41 and event-related potentials. 42 Until recently, it was concluded that studies are insu ciently validated and replicated, and do not yet support the use of EEG for clinical decision making. 43 However, two recent studies using machine-learning algorithms applied to resting-state EEG features identi ed predictive signatures for sertraline, a selective serotonin-reuptake inhibitor, that related differentially to rTMS response. 44,45 This nding is of clinical relevance as it suggests that EEG signatures may be useful as a clinical tool to stratify patients to one of two evidencebased antidepressant treatments (rTMS vs. antidepressant medication), aiming to increase initial treatment response, without the requirement to consider off-label prescriptions or simply 'withhold' treatment due to a biomarker predicting low likelihood of response. 46 Here, our aim was to predict treatment outcome in MDD based on an EEG biomarker using a polygenicinformed EEG data-driven, data-reduction approach: selection of functional brain networks for subsequent response prediction was guided by PRS-MDD, thus combining genetics with neurophysiology approaches. To that end, we conducted a functional independent component analysis (fICA) using LORETA (Low Resolution Brain Electromagnetic Tomography), producing independent spectral-spatial components (i.e. functional brain networks), in a large dataset of severely ill psychiatric patients and healthy controls (dataset 1). In a prior study, this fICA method was tested and validated, 47,48 and demonstrated to reliably identify the default mode network (DMN) and task-positive network (TP) in a sample of 1,397 subjects, replicated in an independent ADHD sample. 47 We then found that one functional network was signi cantly associated with polygenic liability to MDD. Subsequently, validation analyses in two large independent datasets (datasets 2 and 3) demonstrated how this EEG biomarker is differentially associated with antidepressant treatment to rTMS and antidepressant medication. Finally, in simulation approaches we show that this biomarker increases prediction accuracy of antidepressant treatment response and remission, by strati cation to one of two treatments by the EEG biomarker.

Results
An overview of the baseline demographic characteristics and response and remission rates per dataset after EEG preprocessing (online Methods) can be found in Table 1. The analysis procedure that was performed in this study is visualized in Figure 1.

ICA analysis identi es 58 components
Of the 1,195 participants enrolled in dataset 1, the nal sample for the LORETA-fICA analysis after quality control (online Methods) consisted of 1,061 hospital-admitted psychiatric patients and 62 healthy controls (N=1,123; dataset 1). The appropriate dimensionality of the data was established using sphericity test which indicated 58.2 dimensions; hence the LORETA-fICA analysis was constrained to 58 components that explained 96.0% of the total variance in EEG power (see Figure 1: discovery).
Relating components to polygenic risk for MDD Of the 1,123 participants, PRS-MDD association analysis was performed using the data of 722 participants remaining after EEG pre-processing and genetic quality control (Supplementary information). The individual loading score on 1 out of 58 independent components correlated negatively (r=-0.182, R 2 =3.3%) at p<0.001 with PRS-MDD at P T <5.0×10 -8 , after controlling for age and the rst ve genetic principal components (PCs), in women only (see Supplementary information for more details about the PRS model t). This was followed up by an exploratory SNP analysis separately for men and women, which revealed that -for ve variants -the component loading was signi cantly different between homozygotes or heterozygotes for the alternate allele compared to homozygotes for the reference allele (see Supplementary information for all results). Notably, the loading was signi cantly different in subjects with a homozygote variant on the SGIP1 (SH3 domain GRB2 like endophilin interacting protein 1) gene (rs6656912) compared to those homozygotes for the reference allele, which was more pronounced and in opposite direction in men (Cohen's d, d=-0.435; p=0.007) compared to women (d=0.310; p=0.041). Figure 2 shows the component that was associated with PRS-MDD in women, representing jointly deactivation and activation of neural activities coming from sets of regions that form functional spatialspectral networks (see Supplementary information for a non-scaled color map of the networks). Most prominent were frontal alpha power, mainly seen at the left dorsolateral prefrontal cortex (DLPFC), inversely associated with delta and theta power in the right anterior portion of the PFC and delta power seen at a region surrounding the left lateral sulcus, mainly including somatosensory-motor cortices. We will refer to this component as the prefrontal and sensorimotor (PF-SM) network.
The individual loadings on the PF-SM network as visualized in Figure 2, did not correlate with baseline characteristics and a sensitivity analysis revealed that the displayed results cannot be explained by frontal alpha asymmetry, as earlier reported by Arns et al. 35 on the same data (see Supplementary information for the correlation sensitivity analysis), thus ruling out that only the asymmetric alpha was predictive and therefore this represents a novel EEG biomarker.

Relating the PRS-informed EEG component to antidepressant treatment outcomes
The primary outcome for validation analysis (see Figure 1: validation) was dimensional improvement of depressive symptoms using linear regression, based on self-report questionnaire scores at baseline and after rTMS (dataset 2) or medication treatment (dataset 3). All data were normally distributed. The secondary analysis was focused on categorical improvement: response (de ned as ≥50% reduction of baseline severity score) and remission (de ned as a score of ≤12 on the Beck Depression Inventory II, BDI-II, or ≤5 on the Quick Inventory of Depressive Symptomatology, QIDS).

Relating the PRS-informed EEG component to rTMS outcome (dataset 2)
Of the 196 participants, data of 186 were usable (receiving rTMS 1 Hz right DLPFC or 10 Hz left DLPFC, clean EEG and all channels available).
Primary linear regression analysis of ∆BDI-II on individual EEG component loading with age as covariate yielded an R 2 of 7.9% (p=0.007) in women, and R 2 of 8.0% (p=0.005) in men, and R 2 of respectively 6.7% (p=0.009) and 6.4% (p=0.008) when baseline BDI-II score was also added as covariate.  Figure 3B).
A sensitivity analysis con rmed that EEG component loading alone signi cantly contributed to the rTMS response and remission prediction models as outlined above (see Supplementary information for discriminant sensitivity analysis).
To explore if ndings were driven by one of the two rTMS protocols (1 Hz R-DLPFC vs. 10 Hz L-DLPFC rTMS) a sensitivity analysis was performed. For response, running the ANCOVAs as above, adding rTMS protocol as xed factor, yielded a signi cant main effect for 10 Hz rTMS in men only (d=0.963; F(1,34)=9.752; p=0.004), but not for 1 Hz rTMS, albeit the effect was in the same direction (d=0.295). This indicates the effect on response in men was mostly attributable to 10 Hz rTMS. For remission, no signi cant interactions with rTMS protocol were found (see Supplementary information for detailed outline of the interactions and results).

Relating the PRS-informed EEG component to antidepressant medication outcome (dataset 3)
Of the 1,008 participants, data of 535 were usable (treated per protocol, clean EEG and all channels available).
Primary linear regression analysis of ∆QIDS on individual EEG component loading with age as covariate yielded an R 2 of 3.1% (p=0.015) in all subjects (men and women together) receiving sertraline treatment, and R 2 of 2.4% (p=0.019) when baseline QIDS score was also added as covariate. No signi cant associations were found within subjects receiving escitalopram or venlafaxine.
Second, ANCOVA with EEG component loading as dependent variable and response and sex as xed factors, and age as covariate yielded a signi cant response x treatment interaction (F(2,489)=3.278; p=0.039), but no interaction with sex. Repeating the analysis for all three antidepressants separately without treatment as xed factor resulted in a main effect of response for sertraline (d=-0.309; F(1,177)=4.316; p=0.039). A discriminant analysis resulted in a signi cant contribution of age to the prediction of sertraline response in women and men together (Wilk's Lambda=0.970; Chi-Square=5.431; p=0.020). Running the same analysis including the EEG component loading resulted in an improved model (Wilk's Lambda=0.948; Chi-Square=9.618; p=0.008), with a PPV and NPV of 0.60 and 0.55, respectively. The AUC of the ROC analysis was 0.634 (see Figure 3A).
ANCOVA with EEG component loading as dependent variable and remission and sex as xed factors, and age and baseline QIDS score as covariates yielded no signi cant interactions, suggesting the increased network activity is predictive for sertraline response, but not remission when controlled for baseline symptom severity.
A sensitivity analysis con rmed that EEG component loading alone signi cantly contributed to the sertraline prediction model for response as outlined above (see Supplementary information for discriminant sensitivity analysis).

Strati cation demonstrates likelihood of remission is increased when using the EEG component
Based on the merged (signi cant) changes of self-report scores against the PF-SM network loading ( Figure 4A), a simulation was conducted to assess the improvement in treatment outcome, using a cutoff of zero ( Figure 4B). To that end, we calculated what the effect would have been (see Figure 4C) if women with greater (>0) and men with lower (<0) loadings than zero had been treated with rTMS, and if men or women with greater loading (>0) than zero had been prescribed sertraline (see Supplementary information for details of this simulation). The remission rate to rTMS (10 Hz L-DLPFC or 1 Hz R-DLPFC) would have improved from 55% to 69% in women (+26%) and from 56% to 62% in men (+11%). The improvement was largest for 10 Hz rTMS: from 58% to 75% in women (+29%) and from 59% to 70% in men (+18%). The remission rate to sertraline would have improved from 35% to 46% in men and women (+34%). The rTMS response rate in men would have improved from 61% to 68% in men (+11%), and was largest for 10 Hz rTMS, in line with our rTMS protocol sensitivity analysis outlined above: from 68% to 80% (+18%). The response rate to sertraline would have improved from 53% to 69% (+30%). Women with a loading smaller than zero would respond less to both sertraline and rTMS and treatment with escitalopram or venlafaxine would result in slightly improved remission and response rates for them (+5% and +4%, respectively).

Discussion
Here, we aimed to elucidate whether polygenic-informed EEG biomarkers may help predict differential antidepressant treatment response. Using a polygenic risk score-informed data-driven, data-reduction approach applied to resting-state EEG in a large set of hospital-admitted psychiatric patients and healthy controls (dataset 1), we were able to identify one spectral-spatial independent component ('functional network'). Given psychological measures mapping poorly on neurobiology and cognizant of the scarce diagnostic and prognostic biomarkers in MDD, [49][50][51] we have here taken a novel, genetics-informed approach. We thus uncovered a functional network biomarker that in turn was differentially associated with two evidence-based antidepressant treatments in independent datasets.
Visualizing our functional network (Figure 2), we found 1) prefrontal jointly left-sided alpha power (mainly DLPFC) that was inversely associated with right-sided slow delta and theta power (mainly in the anterior portion of the PFC); 2) slow delta power surrounding the left lateral sulcus, including the somatosensorymotor and auditory cortex; 3) asymmetrical alpha activity in the visual cortex. The individual strength of this PF-SM network was differentially associated with treatment outcomes to rTMS in a sex-speci c manner, and to sertraline (similar for men and women), but no such associations for escitalopram or venlafaxine were detected.
In women, but not in men, increased neural activity of this network was related to lower PRS-MDD. The predictive value of the network with regards to treatment outcome was validated in MDD patients receiving TMS treatment (dataset 2) and randomized antidepressant treatment (dataset 3). Primary analysis showed that increased network PF-SM strength was dimensionally associated with response to rTMS in women but to non-response in men, while it was non-sex-speci cally associated with sertraline response. Secondary analysis con rmed these results for remission from MDD. Subsequent discriminant analysis suggested that the PF-SM network loading improved the basic model including the clinical variables age and baseline severity symptom score for the prediction of remission and response. Lastly, based on the relative change on self-reporting questionnaires, a clinical cut-off individual PF-SM network loading of 0 was established. The results of the simulation indicated an improvement in the number of predicted rTMS remitters with 11%-26% and of sertraline remitters with 34%. Rest-EEG recordings and subsequent calculation of PF-SM network loading in treatment-naive MDD patients before treatment inception is likely relatively economical and non-invasive. Such a simulation may thus in future provide a useful construct for treatment strati cation, thereby enhancing chances of initial remission (and response), thus limiting the relative ine ciency of the current stepped-care, 'trial-and-error' approach. For example, if the loading value for a given patient is >0, a clinician may decide to prescribe sertraline and alternatively advise rTMS for female MDD patients. If the loading value is <0, they may choose rTMS for a male patient, but escitalopram or venlafaxine for a female patient. Given that e cacy of antidepressant treatment in the general MDD population is moderate, [52][53][54] and antidepressant discontinuation and switching rates are high, [55][56][57] only slightly increased response chances may reduce disease burden and duration.
Several hypotheses might explain the predictive value of the PF-SM network for antidepressant treatment outcomes in MDD. Abnormalities of the PFC as a network node are known to be implicated in the etiology of MDD and have previously been associated with treatment outcome. 58 TMS applied to the PFC, however, results in transsynaptic activation of deeper areas such as the sgACC, 59 and the frontal-vagal pathway. 60 It is plausible that, by modulating neural activity at the stimulation site, TMS synchronically activates remote cortical areas and thereby modulates dysfunctional functional connectivity between areas of the PF-SM network in a cross-frequency manner. An organized neural circuit between the PFC and motor and somatosensory cortices has been described. 61 Also, TMS induces anticorrelations between the DLPFC and medial prefrontal areas of the default mode network. 62 Notably, the effect we detected was mostly driven by high frequency (10 Hz) stimulation at the left DLPFC, which was previously associated with a network synchronizing effect in the alpha frequency band. 40 This protocolspeci c effect was most predominantly found in men. Furthermore, while increasing PF-SM network activity was signi cantly related to improvement of depressive symptoms after rTMS in women, the reverse effect was found in men. These ndings hint at different underlying mechanisms of action of TMS on neural activity in men relative to women, and are supported by previous studies reporting sexspeci c differences in TMS response. [63][64][65] On a SNP level, the PF-SM network loading was also different in subjects homozygous for the rs6656912 variant (SGIP1), compared to reference allele homozygotes, which was more pronounced and in opposite direction in men compared to women. SGIP1 is expressed predominantly in the brain and encodes a protein required for the neuronal regulation of energy homeostasis. 66 Recently, it was hypothesized that SGIP1 may be involved in processes regulated by Wnt signaling, together with other protein interactors of Wntless (Wl, essential protein for the secretion of multiple Wnt proteins 67 ) that are expressed in the brain, such as the dopamine transporter. 68 Differences in Wnt signaling may partly underlie the reversed treatment outcomes found men and women. Further research is warranted to investigate these sex-speci c mechanisms.
In both sexes, increased PF-SM network activity was antidepressant-speci cally related to sertraline response, but not to escitalopram or venlafaxine. Sertraline, in contrast to the other two antidepressants, is also a synaptic dopamine reuptake inhibitor that increases extracellular levels of dopamine. 69,70 Potential involvement of dopamine in the PF-SM network is supported by the nding that rTMS of the left PFC can induce dopamine release in the striatum. 71 Dopamine plays a prominent role in many functions that are impaired in MDD, such as execution of movement, executive cognitive functioning, and the ability to experience pleasure. 72 Downregulation of the dopamine system, leading to dysfunction of neural circuits such as PFC-amygdala functional connectivity, has been implicated in the pathophysiology of MDD. 73 Evidence also suggests involvement of dopaminergic systems in the modulation of sensorimotor gating, 74 and indicates that MDD is characterized by dysregulation of sensorimotor processes that modulate depressive symptoms via fronto-limbic circuits. 75,76 Thus, the PF-SM network we uncover here may partly re ect a disrupted dopamine system, and therefore allows sertraline outcome prediction.
Cross-validation using three large, independent datasets is an important strength of this study. In addition, the fICA-LORETA method is applicable to all EEGs independently of apparatus, electrode con guration or number of electrodes since it is derived from the voxel-level rather than the electrode level. Furthermore, to allow for future clinical translation of our ndings we have highlighted several clinically intuitive outcome measures that indicate clinical relevance of the EEG component we retrieve.
Nonetheless, limitations of our study include the lack of a placebo control arm, precluding analyses that parse placebo effects. On the other hand, the opposite effects for men and women for rTMS argue against notion of none-speci c effects. Since all patients received psychotherapy concurrent with rTMS, we could not rule out that the EEG component was predictive for rTMS only. Furthermore, for visualization of neural activity, the fICA-LORETA method calculates power on a categorical scale (i.e. frequency bands) instead of a continuous scale (i.e. power spectrum), thereby limiting the interpretation of the functional networks that are obtained by fICA. Finally, while for our strati cation model we relied on an EEG biomarker, future studies should aim to further optimize strati cation by also including other baseline variables, which are likely to further improve the clinical response.
In conclusion, we show for the rst time how a genetics-informed data-driven, data-reduction approach identi es an EEG functional brain network that increases response prediction to two treatments in MDD.
Our method highlights the clinical applicability of such an approach and sets the stage for future strati ed psychiatry research.

Online Methods
Participants, dataset 1 The rst dataset was used for functional independent component analysis (fICA). EEG recordings of participants potentially eligible were collected from September 2013 until September 2018 at Ziekenhuis Netwerk Antwerpen (ZNA), a large community hospital in Antwerp, Belgium. The study was approved by the Institutional Review Board of ZNA. We abided by the principles of the Declaration of Helsinki. A total of 1,195 adult participants -1,132 psychiatric patients and 63 healthy controls -were included and provided written informed consent. Exclusion criteria for all participants were age < 18 years, inability to give informed consent for whatever reason, restlessness that could interfere with the EEG, and inability to sit still. Only patients who were hospital-admitted for stabilization and/or treatment of a psychiatric disease (no further selection was made) were included to allow collection of a representative and heterogenous cohort of (severely ill) psychiatric patients. Healthy controls were de ned as having no current psychiatric episode and never been treated by a mental health service. After preprocessing, the total sample for fICA consisted of 1,123 (1,061 patients and 62 healthy controls).

Participants of the rTMS study, dataset 2
The second dataset was used for validation purposes and the evaluation of treatment effects. It consisted of 196 patients, diagnosed with non-psychotic MDD or dysthymia and BDI-II ≥14 at baseline, receiving protocolized rTMS treatment concurrent with psychotherapy. All participants provided written informed consent. Only participants receiving high-frequency TMS (10 Hz left dorsolateral prefrontal cortex, DLPFC) or low-frequency TMS (1 Hz right DLPFC) were included; participants receiving both 1 Hz and 10 Hz sequentially were excluded since this would make interpretation of results di cult. All patients completed at least 10 sessions of treatment, and lled in the BDI-II at baseline and at the last session (on average session 21). Details about this sample are described elsewhere. 49,77 Participants of the medication study, dataset 3 The third dataset used for validation purposes and the evaluation of treatment effects was an international multi-center, randomized, prospective open-label trial (phase-IV clinical trial), or iSPOT-D sample (International Study to Predict Optimized Treatment in Depression). This study consisted of 1,008 patients diagnosed with non-psychotic MDD who were subsequently randomized to escitalopram, sertraline, or venlafaxine. All participants provided written informed consent and this study was approved by the institutional review boards at all of the participating sites and this trial was registered with ClinicalTrials.gov under id NCT00693849. At baseline and after 8 weeks of treatment patients lled in the Quick Inventory of Depressive Symptomatology (QIDS). Only data from participants who completed 8 weeks of randomized medication treatment ('per protocol' sample) were included. Details about this sample have been published elsewhere. 35,78 DNA isolation and genotyping (dataset 1) DNA was isolated and genotyped for 887 participants who provided informed consent for DNA extraction and analyses. One 10 ml ethylenediaminetetraacetic acid (EDTA) tube was lled during standard blood draws at the ward. DNA was extracted in the clinical molecular genetics laboratory of the University Medical Center Utrecht (UMCU). Samples were brought to a DNA concentration of 50 ng/μl with a total concentration of 200 ng DNA per participant. Subsequently, samples were sent in two batches to the Human Genotyping Facility of Erasmus Medical Center (Erasmus MC) Rotterdam for Global Screening Array v.1 (GSA) by Illumina, Santa Clara (California), USA, that has excellent validity and reliability. 79 Genetic quality control (dataset 1) Quality control (QC) and principal component analysis (PCA) were done with PLINK 1.9, 80 and performed on two batches separately (see Supplementary information). Pre-imputation involved the creation of a superset with the highest quality SNPs for subsequent sample QC. The superset of SNPs was created by excluding those with genotype call rates <0.01, minor allele frequencies (MAF) <0.1, Hardy-Weinberg equilibrium (HWE) <10 -4 , and linkage disequilibrium (LD) r 2 >0.2, with a window size of 50 and window shifting of a step size of 5. Using the superset, subjects were removed who: 1) had a mismatch in their sex between reported and genotyped; 2) were too extremely hetero-or homozygous (their F-values differed ≥3 SDs from the mean in the whole cohort); 3) were related (their pi-hat was above 0.1: one of each pair was randomly excluded); and 4) were cohort outliers (had values for the rst two ancestry principal components (PCs) that deviated ≥3 SDs from the mean of the whole cohort).
This was followed by a regular SNP-level QC for exclusion of ill performing SNPs: variants with genotyping rate <0.01; MAF < 0.01; HWE p-value <10 -5 were thus excluded. European ancestry was checked by comparing with the HapMap population: individuals were removed who deviated ≥3 SDs from the rst and second genetic ancestry PCs from the Northern and Western European (CEU) population.
Lastly, before imputation, genotypic data was compared with the Haplotype Reference Consortium panel. Strands, alleles, positions and frequency differences were checked. Chromosomes were pre-phased and imputed using the Michigan Imputation Server. 81 Post-imputation QC was performed to include reliable SNPs: variants that had a MAF >0.05 and LD r 2 ≥0.8 were included, resulting in 5,211,700 SNPs available to generate PRS.
Polygenic risk score calculation (dataset 1) The polygenic risk score (PRS) for MDD per individual was generated for those SNPs and individuals remaining after QC. The summary statistics of MDD (meta-analysis of PGC and UK Biobank 82 ) were used to generate PRSs. 83 That dataset has no sample overlap with the current dataset as Belgian individuals had not been included in the large MDD GWAS. If only odds ratios (ORs) were reported in the summary statistics, ORs were log-converted to beta values as effect sizes. To that end, the beta values, effective allele, and p-values were extracted from all summary statistics.
SNPs that overlapped between the summary statistics GWASs (training datasets), 1,000 genomes (reference), and our dataset (target) were extracted. Then, insertions or deletions, and ambiguous SNPs, were excluded. To account for complicated LD structure of SNPs in the genome, these SNPs were clumped in two rounds using PLINK 1.90b3z, 84 according to previously established methods; 85,86 round 1 with the default parameters (physical distance threshold 250 kb and LD threshold (r 2 ) 0.5); round 2 with a physical distance threshold of 5,000 kb and LD threshold (r 2 ) 0.2. Additionally, we excluded all SNPs in genomic regions with strong or complex LD structures. Sample overlap between training datasets with our target samples is unlikely since all samples belong to different cohorts and to our knowledge no Belgians had been included in the aforementioned GWASs.

EEG recordings and preprocessing
During EEG recordings, subjects were seated in a sound and light attenuated room that was controlled at an ambient temperature of around 22°C. The participants were instructed to sit still for the duration of the recording without thought instructions. The operator did not intervene when drowsiness patterns were observed in the EEG.
Subsequently, the following steps were taken in the EEG preprocessing and artefact rejection procedure using Brain Vision Analyzer 2.0 (Brain Products, Germany): 1) data ltering: 0.5-90 Hz (dataset 1) or 0.3-100 Hz (dataset 2 and 3), and notch lter; 2) removal and spherical spline interpolation of noisy signals or at lines; 3) electro-oculography (EOG) correction, using a regression-based technique 88 ; 4) segmentation in 4-second epochs; and 4) artefact-rejection using an automatic procedure (criteria: maximal allowed difference of 150 µV peak-to-peak). This resulted in a minimum of one-minute data per subject.

LORETA-fICA model
The EEG was used for estimating the cortical source distribution of electric neuronal activity by means of LORETA (low-resolution electromagnetic tomography; free academic software available at https://www.uzh.ch/keyinst/loreta). This method weights minimum norm inverse solution, and localization inference is based on the standardized estimates of the current density. 89 The following analysis steps were performed using the collection of 4-second artefact-free epochs obtained from dataset 1. In the rst step, each EEG recording was transformed to the frequency domain, using the discrete Fourier transform. The cross-spectral matrices were obtained for six frequency bands, de ned as: delta (1.5-3.5 Hz), theta (4-7.5 Hz), alpha (8-13 Hz), beta (14.5-30 Hz), low-gamma (31-47 Hz), and high-gamma (>70 Hz). Aiming to eliminate the notch bands used at different sites in the EU and US, the 48-69 Hz range was excluded. In the second step, from data of each cross-spectrum matrix, the spectral density was computed for each cortical voxel, sampled at 5 mm resolution in a realistic head model, using the MNI152 template. 48 In the third step, the spectral-spatial data of all subjects was concatenated, and ICA was performed on this data, aiming to identifying independent spectral-spatial components (i.e. functional networks). This method was recently validated in Aoki et al. and Gerrits et al. and reliable identi ed DMN (default mode network) and TP (task-positive) networks. 47,48 The typical ICA model assumes that the source signals are not observable, statistically independent and non-Gaussian, with an unknown, but linear, mixing process, 90 and is described by the following formula: x=As where x, A and s represent matrices. In our case, these three matrices consisted of the following data: 1. Matrix x with 1,123 rows corresponding to all subjects of dataset 1, and the data per subject consists of 37,434 (6,239x6) columns corresponding to the spectral power at 6,239 cortical voxels for the six frequency bands. This approach, using a priori determined frequency bands, is a unique feature of the method used. 91 2. Matrix s with 58 rows corresponding to the number of statistically independent components (i.e. functional networks), and 37,434 columns. In this way, each functional network contains 6 spatial images corresponding to the neural activity of each frequency band (i.e. in a cross-frequency manner).
3. Matrix A with 1,123 rows and 58 columns. Thus, what remains of this data reduction for every subject is the amount of each component that was used for that subject. This amount is expressed as a loading (i.e. signed weight or score) per functional network for each subject.

Independent components
Each independent cross-frequency spectral-spatial functional network (fICA network or EEG component) represents sets of brain regions that are consistently activated or deactivated together within and across a given frequency band. The number of EEG components here was estimated from a measure related to Wackermann's Omega Complexity, 92 indicating 58.2 dimensions, hence the LORETA-fICA analysis was constrained to 58 components explaining 96.0% of total EEG power variance, which is a great reduction of data.
To visualize the functional networks (i.e. correlation of brain regions that are consistently activated or deactivated), a threshold at 3 z-values was set. Individual loadings per fICA network were obtained for each subject, corresponding to the strength of that network for a given individual subject.
The functional networks that were established based on the rst dataset, were prospectively applied to dataset 2 and 3. Likewise, for each subject in each dataset, EEG component loadings were obtained per network. These loadings were used in the statistical analysis.

Outcome measures
For biomarker identi cation (discovery, Figure 1 main text) PRS-MDD was examined for association with independent EEG components (dataset 1). The primary outcome for the prediction analysis (validation, Figure 1 main text) was dimensional improvement of depressive symptoms, based on the BDI-II for rTMS sample (dataset 2) and the QIDS for iSPOT-D sample (dataset 3), which are both self-report questionnaires. In a secondary analysis we focused on categorical improvement de ned as response or remission. Self-reports were taken at intake and after treatment completion (on average at session 21 for rTMS and at week 8 for antidepressants). Response was de ned as ≥50% reduction of baseline score, while remission was de ned as a score of ≤12 on the BDI-II 93  A priori, analyses were strati ed by sex, since previous iSPOT-D studies reported sex-speci c predictors of treatment outcome, 35,41,42,78 and this would enable us to identify strati cation biomarkers. Hence, sex was included as main factor or women and men were analyzed separately if the analysis could not accommodate sex as main factor, rather than handled as covariate since covariation can only resolve quantitative -not qualitative -sex differences. If no sex interaction was found or the effect for both sexes was the same direction, men and women were analyzed together.
First, a discovery analysis was performed to examine if there was an association between one or more fICA components and PRS-MDD (dataset 1). To that end, a partial correlation analysis -controlling for age and the rst ve genetic ancestry principal components (PCs) on men and women separately -was run between all EEG fICA components loadings and the most stringent PRS-MDD p-value thresholds (P T =5.0×10 -8 to P T =0.05) in order to choose the optimal P T , which is unknown a priori. 83 The signi cance level was set to α<0.001, to obtain a highly signi cant result and reduce the likelihood of false positives and non-reproducibility of ndings. 95 Second, validation analyses were performed (dataset 2 and 3) to examine if the EEG component that was found to be signi cantly associated with the PRS-MDD in the discovery analysis, was predictive of treatment outcome (dataset 2 and 3). The signi cance cut-off for the analyses following up was set at conventional α=0.05 as these analyses were intended for validation of the ndings in the discovery analysis. We investigated possible associations between individual EEG component loading and absolute changes in BDI-II score (dataset 2) and QIDS score (dataset 3). The absolute change (∆) was de ned as the symptom severity score difference between baseline and at treatment completion. Therefore, for men and women separately, ∆BDI-II and ∆QIDS were regressed on EEG component loading, adding age and baseline self-report score as covariates. If signi cant associations were found, factorial ANCOVAs were run to investigate if the individual loadings on the EEG component were signi cantly different in remitting and responsive patients relative to non-remitters and non-responders, respectively. In dataset 2, treatment response/remission and sex were added as xed factors; in dataset 3, treatment response/remission, sex and antidepressant subtype were added as xed factors. Age was added as covariate in all models. For remission, baseline BDI-II in dataset 2 and QIDS in dataset 3 (i.e. clinical score) were added as additional covariate, since baseline severity was related to treatment remission, while response was not. 49,78 Subsequently, to assess the predictive value of the EEG component, a discriminant analysis on treatment outcome was performed for datasets 2 and 3. Prior studies had already tested several psychological (personality, anxiety etc.), demographic and behavioral measures and their ability to predict response or remission in these samples, and failed to nd robust and clinically relevant predictors. 49,51 The baseline model consisted of 1) age predicting response and 2) age and baseline severity predicting remission.
Then we tested whether the model performance improved when EEG component loading was added as a predictor. Also, a receiver operating curve (ROC) of both models was regressed on response and remission to establish the predictive value of the discriminant analysis. The positive predictive value (PPV) and negative predictive value (NPV) for predicting response and remission were assessed.
Finally, we attempted to make a construct, based on the data of both dataset 2 and dataset 3. In order to merge both samples, Z-scores of the ∆BDI-II and ∆QIDS were calculated. A clinical cut-off loading was visually determined based on a graph with dataset 2 and 3 merged, representing standardized ∆BDI-II and ∆QIDS against the EEG component loading. Based on this cut-off, a simulation was built in order to evaluate the clinical usefulness, by calculating the percentage response/remission improvement if only subjects above or below this cut-off were assigned to rTMS or sertraline treatment.
Declarations the rst draft. All authors were involved in the writing and revision of further drafts, have approved the submitted version and agreed to be personally accountable for their own contributions.
Competing Interests statement MA is unpaid chairman of the Brainclinics Foundation, a minority shareholder in neuroCare Group (Munich, Germany), and a co-inventor on 4 patent applications related to EEG, neuromodulation and psychophysiology, but receives no royalties related to these patents; Research Institute Brainclinics received research funding from Brain Resource (Sydney, Australia), Urgotech (France) and neuroCare Group (Munich, Germany), and equipment support from Deymed, neuroConn and Magventure. EG is founder and receives income as CEO and chairman for Brain Resource Ltd. and he has stock options in Brain Resource Ltd. Additional funding for this project was obtained through a personal UMC Utrecht Brain Center Rudolf Magnus Young Talent Fellowship (H150) to JL. The remaining authors, HM, BL, GvW, DD, BdW, JvH, PN and, KvE, declare no competing nancial or non-nancial interest.

Data Availability statement
The data that support the ndings of this study are available from the corresponding author, MA, upon reasonable request.