Fifty-three adults with a primary diagnosis of major depressive disorder (MD), as assessed by a psychiatrist with the Structured Clinical Interview for the Diagnostic and Statistical Manual of Mental Disorders-Fourth Edition (Text Revision) DSM-IV-TR [SCID-IV-TR] (29), participated in this study, as previously described (19). Briefly, patients were excluded if they had any other Axis I disorder (except for anxiety disorders), recent (< 6 months ago) problems with substance abuse/dependence, an unstable medical condition, significant suicide risk, seizure history, or if they had been previously treated for their current depressive episode with the study medications. Medicated patients underwent a supervised washout period prior to study commencement (>5 weeks for fluoxetine, 1 week for other medications, consistent with wash-out protocols). As part of a clinical trial conducted between August 2007 and March 2012 (28, 30), patients received either escitalopram (ESC) and placebo, bupropion (BUP) and placebo, or a combination of the two medications for 12 weeks. Assignment to a specific treatment regimen was randomized (double blind).

Depressive symptoms were assessed using the MADRS (31, 32), every week during the first 4 weeks, and biweekly for the remaining 8 weeks. Dosage was increased if tolerated and remission was not yet reached (average dose at 12 weeks for the current sample: dual treatment: ESC = 32 mg, BUP = 379 mg; monotherapy: ESC = 34 mg, BUP = 425 mg). All patients had a baseline MADRS score ≥ 22. Patients whose MADRS scores improved ≥50% from baseline to 12 weeks were considered responders (R), while those who improved <50% were considered non-responders (NR). Due to participant drop-out and issues with EEG data quality, 10 participants were excluded, leaving 43 participants for data analysis (i.e. complete datasets at baseline, week 1 and 12). Of these, 25 were responders and 18 were non-responders. Demographic and clinical characteristics can be found in Table 1 (statistically compared on pertinent variables). All participants provided written informed consent and were reimbursed $30 CAD/testing session. This study was approved by the Royal Ottawa Health Care Group and University of Ottawa Research Ethics Boards.

We conducted power analyses in G*Power (33) to assess the statistical power our sample afforded us. As no established way exists to conduct power analyses for our multivariate statistical approach (partial least squares - see below), we estimated power based on the univariate model that most closely matched our experimental design, namely a 2 (responders vs non-responders) x 3 (baseline, week 1, week 12) mixed ANOVA. We determined that our sample size is sufficient to detect small within-subject and moderate between-subject effects at >80% power, and these calculations are in line with sample size guidelines proposed by Cohen and colleagues (34).

Resting state EEG recordings were collected before the start of treatment (baseline), 1 week and 12 weeks after treatment initiation at the Royal Ottawa Mental Health Centre. Participants abstained from caffeine and nicotine >3 hours prior to testing, and did not take any psychotropic medications, other than the prescribed antidepressants. Two 3-minute resting-state EEG recordings were collected, one with eyes open (EO) and one with eyes closed (EC), while participants sat in a temperature- and light-controlled testing chamber. Ip and colleagues (35) showed that 3-minute EEG recordings are enough to extract reliable EEG characteristics in the theta, alpha and beta bands using a test-retest design. Gudmundsson and colleagues (36) furthermore, examined reliability in EEG spectral power, complexity and connectivity, and found decent reliability with 40 seconds of data included in analysis. The order of EO and EC testing was counterbalanced between participants and sessions. The rationale for these short eyes open/eyes closed recordings was to mitigate the possibility of sleepiness contaminating the data. Furthermore, most data collection was done during the day (i.e., not first thing in the morning or in the evening), though variability existed. EEG was recorded using 32 Ag/AgCl (silver chloride) electrodes embedded in a cap (EasyCap, Inning am Ammersee, Germany), with electrodes positioned according to a variant of the 10-20 system (37). AFz served as the ground, and the average of the two mastoid channels (Tp9/Tp10) was used as the reference. Four additional channels were placed outside the left and right eye canthi, and above and below one eye, to monitor electrooculographic (EOG) activity. Data was sampled at 500 Hz, and impedance was <5KΩ (BrainVision Recorder, Gilching, Germany).

EEG data were preprocessed using EEGLAB v13.4.4b (38) in MATLAB 2014 (The MathWorks, Inc., Natick, Massachusetts). Raw EEG data were bandpass filtered (0.5-55 Hz; slope: 12 dB/octave) using ERPlab’s IIR butterworth filter, notch filtered at 60 Hz (lower and upper edge: 55-65 Hz) using EEGlab’s basic FIR filter and segmented into 2s epochs. We used an independent component analysis (ICA) to identify and eliminate noise and ocular artifacts. Channels with excessive noise or drift were excluded from the ICA procedure, and subsequently interpolated using EEGlab’s spherical spline interpolation function. The average number of interpolated channels was 0.11, and no more than two channels were interpolated for each participant and session (34). Epochs were visually inspected following ICA, and those with remaining artifacts were manually rejected. An average of 85.7 (range: 63-113) artifact free epochs were obtained per participant, state (EC/EO) and session (baseline, week 1 & 12), including 28 electrodes (Fp1/2; F3/4; F7/8; FC1/2; FC5/6; C3/4; CP1/2; CP5/6; P3/4; P7/8; T7/8; O1/2; Fz/Cz/Pz/Oz). The two reference channels and two additional channels that provided low data quality for most participants (FT9/FT10) were excluded from analysis. There was no statistical difference between responders and non-responders in number of artifact-free epochs per session or state (p-values:.09-.9).

Functional connectivity, as quantified by the weighted phase lag index (WPLI), was calculated for each unique combination of the 28 channels (378 pairs) using the open source Fieldtrip toolbox (39) in MATLAB. WPLI is a modified version of the phase lag index (PLI), which was first described in 2007 by Stam and colleagues (40). It estimates connectivity by calculating the phase angle difference between EEG signals from two channels for each time point, and determining the consistency in these phase lags over time. As such, if the difference in phase between two channels is similar over time, the PLI will be high, indicating high connectivity between two channels. An advantage of the PLI compared to other EEG connectivity measures is that it is less sensitive to volume conduction, because it disregards any phase lags of 0 and π. The WPLI also takes into account that phase lags can easily turn into leads and vice versa (e.g. a slightly positive phase angle difference can turn into a slightly negative phase angle difference). While the PLI is sensitive to such small disturbances in phase lags, the WPLI resolves this issue by giving greater weight to angle differences around 0.5π and 1.5π (41). The result is a value between 0 and 1, with higher values indicating stronger connectivity. WPLI was calculated in two ways, first across epochs, as is commonly used and enables direct comparisions with previous findings, and then within epochs, which is less commonly used, but enables single-participant analyses.

Multiscale entropy (MSE) was used to quantify brain signal complexity. An advantage of MSE over other measures of complexity is that it incorporates multiple time scales. This feature is important because it differentiates between signals that are purely random (e.g., white noise) and those comprised of both random and deterministic components (e.g., 1/f or coloured noise). Signals that are purely random show a rapid decline in the MSE curve with increasing scale whereas those with temporal inter-dependencies will have a more gradual shift in the MSE curve (42, 43). A detailed description and theoretic background for MSE is outlined in Costa and colleagues (11). In short, MSE estimates the regularity of a signal by evaluating the ratio of similar patterns of different lengths repeating over several time scales. It is calculated in two steps. First, the raw signal is resampled several times to create data sequences that represent different temporal scales. Essentially, an increasing number of non-overlapping data points are averaged into one new data point. The first timescale is the (cleaned) raw time series. With a sample rate of 500Hz in the current study, time scale 1 had a temporal resolution of 2 milliseconds between data points. For time scale 2, two consecutive data points were averaged, yielding a temporal resolution of 4 milliseconds; for time scale 3, averaging occurs over three time points yielding a temporal resolution of 6 milliseconds, and so on. The coarsest scale used in this study was 20 (temporal resolution of 40 milliseconds), to ensure a sufficient number of data points (minimum 50) for the sample entropy calculation.

Next, sample entropy is calculated at each time scale. Sample entropy determines the natural logarithm of the ratio of patterns of length m over patterns of length m+1 repeated within one epoch. This gives a value between 0 and 1, with higher numbers indicating a less predictable/more variable signal (i.e. fewer patterns of length m+1 compared to the number of patterns of length m). In line with previous studies, (e.g. 19, 44) and guidelines outlined by Richman and Moorman (45), we set parameter m to 2, while the similarity criterion r, which determines which points in the time series are considered to be ‘the same’, was set to 0.5 (i.e. two data points were treated as indistinguishable if their amplitudes differed <50% of the standard deviation of the time series). MSE was calculated for each epoch and electrode at each time scale, using the algorithm available at www.physionet.org/physiotools/mse/. Single epoch MSE data were used for statistical analyses at the individual level. MSE values were also averaged over epochs to provide one MSE value for each electrode and time scale per participant, session (baseline, week 1 and 12) and state (EC and EO), which were used for group level analyses.

To control for differences in age between responder and non-responder groups (see Table 1), and because both brain signal complexity and connectivity have been observed to change with age (19, 46–49), age was regressed out of the data before the statistical group comparisons using an in-house MATLAB script (50, 51).

Consistent with MSE and connectivity differences between EO and EC states observed in previous studies (19, 52–55), we found strong EO/EC effects in our analyses that masked changes occurring over assessment sessions (see Supplementary Figures S1, S2). Therefore, we performed analyses on EO and EC data separately. To maximize the chance of replicating group findings at the individual level, we performed single-participant analyses on the data showing the strongest effects (EC for WPLI, EO for MSE), and present those group findings below (other group findings are presented in Supplementary Materials [Figures S3, S4]).

Partial least squares with singular value decomposition (PLS-SVD) is a multivariate statistical approach that can detect condition- and/or group-related differences in whole-brain variables (56, 57). Briefly, PLS-SVD calculates the normalized average across participants within each condition (in this case, responder status and assessment sessions) for each element in the brain characteristic matrices (in this case, frequencies and electrode pairs for WPLI, or time scales and electrodes for MSE). Then, this mean-centered matrix is decomposed using SVD into orthogonal latent variables (LVs) that account for most of the covariance between groups/conditions and brain characteristics, revealing the optimal associations between specific groups/conditions and spatiotemporal patterns in the brain. LVs contain several components. One is the singular value, which indicates the strength of the effect the LV represents. Another component holds the condition contrast, which reveals the linear combination of weights capturing differences between conditions (groups and assessment sessions in our data). The third component contains the element loadings (singular image), which represent the pattern of the specific data elements (in this study frequencies and electrode pairs for WPLI, and time scales and electrodes for MSE) that show the given contrast. These element loadings are used to compute brain scores: the dot product of the element loadings with each participants’ data for each assessment session. Brain scores represent the extent to which each participant expresses the given contrast in a single number per participant and condition (see Figure 1C).

Statistical testing occurs at two levels in PLS-SVD analyses. First, the overall significance of the LV is determined using permutation tests. In each permutation, the data are randomly shuffled between conditions (within participants) and between groups, and PLS-SVD analysis is performed on the shuffled data just as on the actual data. LVs are considered significant when their singular value is more extreme than 95% of the singular values calculated from the randomly shuffled data (corresponding to p <.05). In the current study, 500 permutations were performed for each analysis. Second, the stability of the identified pattern across participants is established through bootstrap resampling. In essence, the PLS-SVD analysis is repeated with different subsamples of participants, to see how consistently each electrode pair/electrode and frequency/time scale display the identified pattern of differences across the whole sample. This consistency is quantified as a bootstrap ratio (BSR), which is calculated by dividing the element loadings by the standard error of the created bootstrap distribution for each element. In addition to determining the stability of the pattern, bootstrap resampling also protects against the influence of outliers, as subsamples with and without the outlier would produce different outcomes, thereby decreasing the consistency of the findings (i.e. the bootstrap ratio). In practice, this means that effects that are driven largely by an outlier get attenuated. Bootstrap ratios are similar to z-scores, with absolute values ≥ 3.1 corresponding to ~99% confidence interval. In this study, bootstrap resampling was performed 200 times. As each statistical test is computed in one mathematical step, no correction for multiple comparisons is necessary (56). P-values indicating significance levels, and percentage of crossblock covariance explained (PCCE) are reported for each LV of interest. PLS-SVD analyses were applied both at a group and individual level.

Both groups (responders vs. non-responders) and all sessions (baseline, 1 & 12 weeks of treatment) were entered in four PLS-SVD analyses: two for connectivity (EO/EC states separately) and two for complexity (EO/EC). As all showed interaction effects between groups and assessment sessions, two additional analyses were run for each analysis for responders and non-responders separately, again including all sessions. The p-values of these follow-up analyses were corrected for multiple comparisons using the Bonferroni method.

The input data consisted of across-epoch WPLI/averaged MSE values, organized into 2D matrices with n * k rows, and m * t columns, with n being the number of participants (R: 25; NR:18) and k the number of conditions (assessment sessions: 3). M and t represent the spatiotemporal elements, namely the number of electrode pairs (378) and frequencies (99) for the WPLI analyses and the number of electrodes (28) and timescales (20) for the MSE analyses. The same procedure was followed for the averaged, single-epoch WPLI data (Methods – Connectivity analyses).

Non-rotated (hypothesis-driven) PLS analyses were performed for each individual, using single-epoch EC WPLI data, and single-epoch EO MSE data. Non-rotated PLS was chosen because it allows one to determine whether and to what extent a specific, predefined contrast is present in the data (57). Using the contrasts (patterns of differences in connectivity and complexity between groups and across sessions) found in the group analyses, non-rotated PLS was used to test whether each individual followed the pattern of change observed in responder and non-responder groups. Single epoch WPLI/MSE values were organized into similar 2D matrices as described for the group analyses, only now each participant had their own matrix, with the n dimension representing the number of epochs instead of participants (57). No correction for multiple comparisons was applied, as these analyses aimed to replicate group findings in separate datasets for each individual.

The similarity of the individual PLS outcomes to the group PLS-SVD outcomes was quantified in two independent ways to balance the advantages and suitability of quantitative and qualitative measures (see Figure 1D). First, the similarity was estimated quantitatively, by correlating the stable (|BSR| > 2, corresponding to ~95% confidence interval) element loadings (i.e. the spatiotemporal brain pattern) of the group results with the element loadings of each participants’ individual analysis in MATLAB. For connectivity, the element loadings from the group analyses on averaged single-epoch WPLI were used for this correlation procedure (Methods – Connectivity analyses). For participants whose PLS analysis was non-significant (p >.05) or did not match the predefined contrast (responder/non-responder), indicating a different timing and/or direction of the change highlighted by element loadings, the element loadings were not correlated with the group results and were included as ‘showing no correlation with the group pattern’ in the summaries. The percentage of participants showing moderate-strong correlations (r≥.4; 58) with their own group outcome was considered as a quantitative indicator of how well group results translated to individuals. Second, to balance arbitrary cut-offs, significant individual outcome patterns were visualized and classified by two independent raters, blind to response status, as being similar to either or both the responder or non-responder group patterns, or neither. Important responder and non-responder features were selected based on visual inspection of the most consistent changes across time (i.e. those with |BSR| > 3.1) in the group analyses, and used to construct the rating criteria that both independent raters used to determine similarity. The percentage of participants being classified as conforming to their own group pattern exclusively was determined as a qualitative indicator of the replicability of the group patterns at the individual level.

By design, responders had lower MADRS scores at week 12, but not at baseline or week 1 (Table 1). Apart from responders being younger than non-responders, the two groups did not differ statistically in clinical and demographic characteristics (Table 1). We accounted for the age difference by regressing age effects out of our data before running the statistical tests at the group level.

The PLS-SVD analysis including both groups and all sessions identified one significant LV (p <.001, PCCE = 35.07%). As this LV presented an interaction effect between groups and sessions, two additional analyses for each group separately were run, with the statistical significance threshold corrected to α <.025. These analyses revealed a complex, opposite pattern of change from weeks 1 to 12 in responders (p=.024, PCCE=55.8%), but only approached significance in non-responders after Bonferroni correction (p = .032, PCCE=56.6%; Figure 2). The most prominent frequencies for each group are highlighted by red boxes in Figure 2: Non-responders showed a widespread increase in alpha connectivity from weeks 1 to 12 (10Hz), while responders exhibited an extensive increase in beta connectivity (22Hz). Considering the same frequencies in the opposite groups (e.g. alpha in responders; highlighted by blue boxes) revealed more spatially contained changes in the opposite direction: Responders showed a decrease in connectivity at 10Hz, while non-responders showed a decrease at 22Hz. In both groups, changes in alpha connectivity were most pronounced at interhemispheric frontal-to-occipito-parietal electrode pairs but involved.

additional electrode pairs in non-responders. The most consistent beta changes occurred in left intra-hemispheric connections in both groups, but also included right central and parietal channel pairs in responders (Figure 3).

The group PLS-SVD analyses performed on the averaged single-epoch WPLI data showed a similar pattern of change in connectivity in responder and non-responders (Supplementary Material and Figure S5). However, the results spread over multiple frequencies (e.g. from 8-14 Hz instead of dominantly at 10 Hz), which is unsurprising, considering the reduced spectral resolution associated with sliding window approaches (34).

The PLS-SVD analysis examining changes in MSE over time in responders and non-responders identified one significant LV (p <.001, PCCE = 82.71%), which revealed an interaction effect. The analyses exploring changes for each group separately each found one significant LV (responders: p = .006, PCCE = 93.85%, non-responders: p = .02, PCCE = 86.9%, significant at α <.025). Both groups showed a decrease in coarse scale complexity from baseline to 12 weeks, but the timing and extent of change differed. Responders showed an early (starting at week 1) and widespread decrease in coarse scale complexity, while non-responders showed a later (only present at week 12) decrease in coarse scale complexity in limited electrodes (Figures 4 & 5). Additionally, fine scale complexity increased only in non-responders.

The pattern of change across assessments identified by the group PLS-SVD examining connectivity was similar regardless of the approach used to calculate WPLI. Specifically, it consisted of a change in connectivity from week 1 to week 12 in both non-responders and responders (Figure 2; Supplementary Figure S5). Therefore, we applied this contrast in the non-rotated single-participant PLS-SVD analyses (Figure 1C). As the non-responder and responder patterns only differed in the direction of change (i.e., increase or decrease in WPLI), only one contrast was defined for each analysis (0 1 -1). This contrast examines changes in WPLI from week 1 to week 12 but leaves the direction of change and at which frequencies this occurs to be determined by the data. All non-responders and 22/25 responders exhibited the predefined pattern of change at an uncorrected significance level (all p <.05), of whom 33 (19R/14NR) survived Bonferroni correction (p <.001).

The similarity of the pattern of connectivity across channel pairs and frequencies between the individual and group-level results was examined in two ways (see section 2.10 and Figure 1D for details). First, each individual’s connectivity pattern was correlated with the connectivity pattern identified in the group analyses. This procedure showed that 60.5% of individual patients exhibited moderate-strong positive correlations (i.e., r ≥ .4; 58) between their individual and group outcomes. Another 9.3% showed weak positive correlations (i.e., .1 < r < .4), while 14.0% revealed negative correlations between individual and group PLS outcomes. The remaining 16.3% of patients’ individual analyses either correlated negligibly (-.1 < r < .1) or did not reach significance and were therefore not correlated (Figure 6).

Second, the main characteristics differentiating responders and non-responders were defined and their presence in individual results were visually inspected and rated by two independent, blind raters. Following the direction of change found at 10Hz (alpha) and 22Hz (beta) for the responder and non-responder groups, individuals showing a pattern of a meaningful decrease in alpha (~8-14Hz) and/or increase in beta (~18-30Hz) WPLI from week 1 to 12 were considered to fit the responder pattern, while patients showing a pattern of a meaningful increase in alpha and/or decrease in beta WPLI from week 1 to 12 were considered to fit the non-responder pattern. Two authors (K.C. and G.W.) examined each individual outcome pattern and independently decided whether they conformed to the outlined definitions. Inter-rater reliability was quantified using Cohen’s Kappa: κ=0.76 for rating whether patients fit the responder pattern, and κ=0.57 for rating whether patients fit the non-responder pattern. Raters discussed discrepancies until consensus, which indicated that 39.5% of patients fit the pattern of their own group exclusively, 34.9% fit both patterns, 14.0% only showed the pattern of the opposite groups, and the remaining 11.6% did not conform to either pattern (7% did not reach significance).

Table 2 shows the characteristics of patients who were rated as fitting their own WPLI pattern exclusively versus the other categories. Aside from individual responders who conformed to the responder group pattern showing higher correlations with their own WPLI group pattern compared to responders assigned to other categories (p = .005, did not survive Holm’s sequential Bonferroni test for multiple comparisons), there were no obvious clinical and demographic differences between individuals in the different responder and non-responder categories. We also illustrate the correlations between individual WPLI patterns and the WPLI pattern of participants’ own groups along the four categories in Table 2, colour coded according to treatment regimen and sex in Figure 7 (top row). No clear patterns of individual variation in WPLI related to sex or treatment regimen emerged.

Non-rotated PLS was performed for all participants with two predefined contrasts: a linear contrast (1 0 -1) as the responder pattern (change across the three time points, top of Figure 4) and 1 1 -2 as the non-responder pattern (no change from baseline to week 1, change from weeks 1 to 12, top of Figure 4). While all individual analyses revealed significant results for both predefined LVs (all p <.001), 16% of patients did not show the same pattern of change over assessment sessions as defined in their own group contrast.

The similarity of the pattern of complexity across channels and time scales between the individual and group-level results was examined in two ways (see section 2.10 and Figure 1D for details). First, each individual’s complexity pattern was correlated with the complexity patterns identified in the group analyses. This procedure showed that 53.5% of individual outcome patterns correlated positively with moderate-high strength (r ≥ .4; 58) to that of their groups. Another 9.3% showed a weak (.1 < r < .4) positive correlation between their individual outcome and that of their group. Of the remaining individuals’ analyses, 11.6% yielded negative correlations, while 25.6% either showed negligible correlations (-.1 < r < .1) or did not show the predefined contrast and were therefore not correlated with group patterns (Figure 6).

Second, the main characteristics differentiating responders and non-responders were defined and their presence in individual results were visually inspected and rated by two independent, blind raters. Participants were characterized as fitting the responder pattern if they showed a meaningful decrease in coarse scale MSE from baseline to week 1 in their first LV, and characterized as fitting the non-responder pattern if they demonstrated no meaningful change in coarse scale MSE from baseline to week 1, and any change in coarse scale MSE from week 1 to week 12 in their second LV. Again, two raters (K.C. and G.W.) examined the significant patterns found in each individual analysis and independently decided whether they conformed to these definitions. Inter-rater reliability yielded κ=0.91 for rating whether patients fit the responder pattern, and κ=0.96 for ratings on whether patients fit the non-responder pattern. Based on consensus, 46.5% of patients exclusively showed the pattern of their own group, 25.6% fit both patterns, 14.0% exclusively showed the opposite pattern, and 14.0% showed neither. Patients in this last group, including 4 responders and 2 non-responders, showed an early increase in coarse scale complexity instead.

The characteristics of patients rated as conforming to their own group pattern exclusively versus those who did not are presented in Table 3. Again, the most notable difference was that responders exclusively showing the responder pattern exhibited higher correlations between their individual and group MSE patterns compared to responders assigned to the other categories; this comparison did not survive Holm’s sequential Bonferroni test to correct for multiple comparisons (p = .006). The correlations between individual MSE patterns and the MSE pattern of participants’ own groups along the four categories in Table 3, colour coded according to treatment regimen and sex are presented in Figure 7 (bottom row). We observed no noticeable patterns of individual variation in MSE related to treatment regimen or sex.

Despite the novelty of the presented work, certain limitations exist. First, we grouped patients receiving different treatment regimens (i.e. escitalopram, bupropion or both) to ensure sufficient statistical power for our main analyses. Similarly, as the clinical trial was designed to compare different treatment arms, we were not able to differentiate between treatment and placebo effects, a noteworthy future direction. Our sample was also largely White, and from WEIRD (Western, Educated, Industrialized, Rich and Democratic; 82) societies, and thus findings cannot be generalized outside of such populations. While our sample was fairly balanced in terms of sex, we were also unable to test for sex differences due to power issues. In addition, our sample included patients with different MD symptom subtypes (e.g. melancholic, atypical) and several patients had comorbid anxiety disorders (given the high co-occurrence of anxiety in depressed individuals these individuals were not excluded). Although this heterogeneity in the sample complicates statistical analyses, it does give an accurate representation of the population seeking treatment for MD. At the same time, the treatment provided was not representative of what people with comorbid conditions may receive in clinical practice, as participants were not allowed to take psychotropic medications outside of those prescribed for the study. While heterogeneity may also have contributed to the individual variation we found in the individual analyses, we did not find any indication that individual variation in EEG connectivity and complexity patterns could simply be explained by differences in treatment regimen, sex, clinical profiles or other clinical and demographic characteristics (see Tables 2, 3; Figure 7). That being said, future studies with larger samples should be conducted to explore the influence of heterogeneity of treatment and patient characteristics in more detail. They might also collect data during more narrow time windows during the day, or control for sleepiness/drowsiness in other ways (e.g., as a covariate in the analyses), which we were unable to do here. Similarly, it could be relevant to include more follow-up measurements across longer time periods, other types of individual information (e.g., neuroimaging) and outcome metrics (e.g., functional measures) in such studies, as well as explore alternative analysis methods (e.g., classification methods using Bayesian or machine learning approaches).

In addition, there are numerous ways to quantify similarity between individual and group patterns. The correlation procedure was objective, but also led to arbitrary limits for categorizing who did and did not match the group patterns (i.e. r ≥.4). The independent ratings avoided setting such arbitrary limits. While inter-rater reliability was high for ratings of MSE (κ = 0.91-0.96), there was less agreement for the ratings of WPLI (κ = 0.57-0.76), highlighting the subjectivity of this method. The lower inter-rater reliability for WPLI was likely due to the larger number of elements included in this analysis (378*99 compared to 28*20 for MSE), or could reflect more variety in WPLI response patterns compared to MSE patterns. Finally, we applied no correction for multiple comparisons in the individual analyses, as this would have made it more difficult to find the group patterns in individuals. However, lowering the significance threshold to α = .001 did not alter our MSE findings, and only changed 16% of the individual WPLI analyses (7 patients) from being significant to being insignificant. Overall, both measures of similarity showed that there was substantial individual variation in the connectivity and complexity group patterns.