Quantification of Kuramoto Coupling Between Intrinsic Brain Networks Applied to fMRI Data in Major Depressive Disorder

Organized patterns of system-wide neural activity adapt fluently within the brain to adjust behavioral performance to environmental demands. In major depressive disorder (MD), markedly different co-activation patterns across the brain emerge from a rather similar structural substrate. Despite the application of advanced methods to describe the functional architecture, e.g., between intrinsic brain networks (IBNs), the underlying mechanisms mediating these differences remain elusive. Here we propose a novel complementary approach for quantifying the functional relations between IBNs based on the Kuramoto model. We directly estimate the Kuramoto coupling parameters (K) from IBN time courses derived from empirical fMRI data in 24 MD patients and 24 healthy controls. We find a large pattern with a significant number of Ks depending on the disease severity score Hamilton D, as assessed by permutation testing. We successfully reproduced the dependency in an independent test data set of 44 MD patients and 37 healthy controls. Comparing the results to functional connectivity from partial correlations (FC), to phase synchrony (PS) as well as to first order auto-regressive measures (AR) between the same IBNs did not show similar correlations. In subsequent validation experiments with artificial data we find that a ground truth of parametric dependencies on artificial regressors can be recovered. The results indicate that the calculation of Ks can be a useful addition to standard methods of quantifying the brain's functional architecture.

The relevance of the dependence of Ham-D on head movement is discussed in section 4.1 of the Supplement. Fig. 1 shows the 20 IBNs detected in the empirical data set. Fig. 2 shows the 28 IBNs detected in the test data set.

Metastability Analysis
In order to measure metastability, the averaged synchronization of extracted IBNs was measured using the Kuramoto order parameter R (Deco et al., 2017). All IBN time courses were band pass filtered to 6 frequency bins using a 7th order Butterworth filter: 0.01 Hz -0.025 Hz, 0.025 Hz -0.05 Hz, 0.05 Hz -0.075 Hz, 0.075 Hz -0.1 Hz, 0.1 Hz -0.125 Hz. R(t) across all filtered IBN time courses was calculated per band. R(t) measures global level of synchronization of n oscillating Figure 1: IBN layouts of the 20 independent components (ICs) identified in the exploratory data set, obtained from spatial group-ICA, plotted as z-scores, thresholded at z>1, and displayed at the three most informative slices (MNI space). Networks are categorized into groups according to anatomical and functional properties. IC numbers refer to the matching networks in Allen et al., 2011. Figure 2: IBN layouts of the 28 independent components (ICs) identified in the test data set, obtained from spatial group-ICA, plotted as z-scores, thresholded at z>1, and displayed at the three most informative slices (MNI space). Networks are categorized into groups according to anatomical and functional properties. IC numbers refer to the matching networks in Allen et al., 2011. Figure 3: Median metastability measures per participant group and frequency bin with 50% confidence intervals. Additionally, the overall means and standard deviations across all patients and all controls, respectively, are indicated as line plot. Relevant significant figures marked. Analysis on the exploratory data set revealed a significant metastability reduction at frequencies 0.05 Hz to 0.075 Hz, which was tested and confirmed in the test data set.
signals. When a lack of synchronization is present the n phases are randomly distributed equating to R valued at near 0. If all n oscillating signals are synchronized, this will result in R = 1 (full synchronization). The mathematical expression of the Kuramoto order parameter is given by: where φ k (t) is the instantaneous phase of each narrow-band BOLD signal at node k. The instantaneous phase φ k (t) of each narrowband signal was computed using the Hilbert transform (as implemented in matlab, www.mathworks.com). Values of R(t) were streched from the range of 0 to 1 to the range of -1 to 1 and Fisher Z transformed to the full range of real numbers in order the measure the standard deviation across time. In taking the standard deviation across time, the number of transitions from synchronized to desynchronized states can be quantified, representing metastability. Statistical analysis was first performed on the exploratory data set. Metastability values per bin were tested for normality by a Lilliefors test implemented in matlab (www.mathworks.com). Distributions proved not to be normal. Therefore, statistical group comparisons between patients and controls per each of the nine frequency bins were performed as Wilcoxon rank-sum tests implemented in matlab. Here, group differences were regarded as significant with at P c < 0.05, Bonferroni corrected for the 6 tests resulting from the splitting into 6 frequency bands. The detected significant group differences provided the a very specific a priori hypothesis for the independent test data set regarding a metastability reduction in frequencies of 0.05 Hz to 0.075 Hz. Based on this, the group difference in the test sample was tested for significance by a right-tailed Wilcoxon rank-sum testing without Bonferroni correction for number of frequency bins and proofed to be significant (see Fig. 3).
In consequence we performed our analysis of K in the frequency bin of 0.05 Hz to 0.075 Hz.
2 Calculation of Kuramoto coupling parameters K

Kuramoto model
Synchronization denotes the phenomenon of several processes with individual rhythms coming into co-occurrence with a common rhythm after interacting dynamically. This effect has been observed and described in many different fields of research such as sociology, biology, technology and ecology (Arenas et al., 2008). The first attempt to describe synchronization in mathematical terms was done in 1967 by Winfree (Winfree, 1967) as he viewed the interacting units as phase oscillators. Since then, many synchronization-based models have been established of which the Kuramoto model (Kuramoto, 1975) is one of the most widely investigated (Acebrón et al., 2005;Bocharova, 2014). In the classical so called mean field coupling form it can be written aṡ where the interaction between oscillators ϕ i (t), i = 1, . . . , r is determined by the constant coupling strength C. We allow individual coupling by replacing C by coefficients K i,j . Furthermore, model (2) incorporates that each oscillator rotates at its own natural frequencies ω i which is taken from a distribution g(ω). While the coupling strength forces the frequencies to assimilate, the natural frequencies work as the opposite force. The larger the variance of g(ω), the less likely will the phases of the oscillators synchronize. The time courses in our experiments are filtered to a very narrow frequency band. Therefore, we chose to model the eigen-frequencies as the mean frequency of the respective frequency band. Another option would be to estimate each eigen-frequency from the largest peak in the frequency profile of each time course via fast Fourier transform (FFT). Our codes provide both options.

Euler's Method
Given the ODEφ Euler's method (Epperson, 2013) for a step-wise approximation of the solution reads as followŝ where h is a fixed step size. The smaller h, the more accurate is the approximationφ to the real solution ϕ. If we tried to solve our system of ODEs with this method, we would have to address stability questions, as for stiff ODEs, the Euler method is only stable for extremely small step sizes. In this work, however, we start with the exact values of the time courses and calculate the coefficients of the equation. Therefore, in the end, the step size h will only affect the calculations as a scaling factor for the coefficients K i,j , which will become apparent, when we rearrange the equations to the LES (see Sec. 2.3). Furthermore, we are only considering a very short time span, and only a relatively small number of steps, which alleviates concerns about stability. In future work investigating with other methods such as Runge-Kutta procedures would of course be possible as well.

Derivation of the LES
We derive the linear equation system (LES) matrix and inhomogeneity, starting from the following equation, which was derived by plugging in the phase values into Euler's method. As a first step, we bring ϕ i (s) and ω i to the left side of the equation and multiply the whole expression by 20. The factor 20 could be omitted, since it will only lead to an equal scaling of the coefficients. Nevertheless, we included it in our calculations to account for the 20 brain regions.
We then rewrite the sum on the right hand side as a vector-vector multiplication.

vector of unknown coefficients
For each s = 0, . . . , 298 we obtain one entry for b i and one corresponding row for S i . This leads to the following matrix-vector form: Note, that the matrix S i does not contain the entry for t = 299. This is due to the fact, that we only obtain 299 entries for b i , since it contains differences of consecutive ϕ i values. We have to omit the last entry in order to have matching dimensions of the system matrix and the inhomogeneity vector.
Here, we also see, that the matrix S i will contain a zero column. To be precise, the i-th column will only contain zero values, as the sine of zero is zero. This results in a singular matrix, for which we cannot calculate a solution vector. Therefore, we have to eliminate this column and the corresponding unknown variable K i,i , which we simply set to a fixed value of 1.

Optimized Solution -Gaussian Normal Equations
In order to obtain the coefficients K i,j , i, j = 1, . . . , 20, i = j we have to solve an over-determined LES, which can be written in a general form as: where the system matrix S ∈ R n×m with n ≥ m, the solution vector x ∈ R m×1 and the inhomogeneity b ∈ R n×1 . This means we have more equations than unknowns. We proposed to accept an approximated solutionx for which we obtain the least squared error. Thus, if we define the error of x in the k-th row, or equation respectively, we want x to be chosen such that This is an equivalent formulation to findingx with since the sum in equation (4) is precisely the squared 2 norm in equation (5). For a minimum of E(x) the derivative of E(x) with respect to x has to vanish In matrix form E(x) can be written as The derivative with respect to x is then Plugging in this derivative into equation (6) and dividing by 2 yields Thus, by solving the so called Gaussian normal equations (7) we obtain an approximated solutionx for the original problem (3), which is optimal in the sense of least squared errors.

Kuramoto coupling coefficients
The calculated Kuramoto coupling coefficients (K res , see Sec. 2.2 in the main paper) of one subject can be nicely visualized with a heatmap, where we can directly compare estimated coupling strengths between different IBNs. See Fig. 4 as an example for one subject.
3 Artificial data simulation 3.1 Solving the simulation model to generate phase courses To obtain simulated phase courses, we have to solve the system of ODEs given in Appendix A, Equ.
(15) of the main paper. As the numerical solver of the ODE system, we choose not to work with the same as when calculating the coefficients (i.e., Euler's method). This prevents to simply get out what we put in. The solver used for the simulations is the classical Runge-Kutta algorithm (Schwarz and Köckler, 2011) also called RK-4 with step size 1. One iterations step of the described simulation is given by the following formula ϕ i (s + 1) = ϕ i (s) + 1 6 · (g 1 + 2 · g 2 + 2 · g 3 + g 4 ) + n · ε i (s).
We define a function g as with θ, f and w real numbers, a, l and ∈ R r . The terms g 1 to g 4 in formula (8) are then given by Providing initial values ϕ i (0), i = 1, . . . , r, formula (8) iteratively yields the phase courses. Note, that we switched from t to s as the function variables, since we are now considering iteration steps instead of time steps. k i and m i denote the i−th row of the matrices K, M, respectively and ϕ(s) is the vector of all r phase values at iteration step s.

Technical parameters
In this section, we provide more details about the chosen parameters of the artificial data simulation (see Appendix A of the main paper). The size of our data set is chosen as s × r × T = 24 × 20 × 300, which represents s = 24 instances (representing the subjects) with r = 20 phase courses (representing the IBN's activity phase courses) with T = 300 measure points, where each measure point corresponds to an activity value every two seconds. This is the same size as the exploratory data set and corresponds to a typical clinical study setting. We set the number of iteration steps in the RK-4 procedure to 3000 and afterwards cut off the first 300 values to eliminate possible transient phenomena. We then down sample by choosing every 9 th value to get a time series of length 300 as for the empirical data.The initial phases for each region of each subject are randomly chosen within [−π, π) and the eigenfrequencies are drawn from a Gaussian distribution N (0.0625, 0.0125). This corresponds to the frequency band with the statistically highly significant finding for the empirical data. The frequency band needs to be chosen rather narrow since this is a prerequisite of the Kuramoto model. For each of the 24 subjects the 20 × 20 coupling coefficients matrix K o is randomly initialized with values ranging from approximately −20 to 20 as this is the range for most of the calculated Kuramoto coefficients derived with our method for data in the same frequency band. Each subject is randomly assigned a score value (representing the Ham-D) between 0 and 40 resulting in the vector s.
To maintain a certain balance between synchronization and de-synchronization, the weight factor d in equation (8) needed to be set such that the oscillators deviate from their respective eigenfrequencies, but at the same time the influence needs to be damped to prevent highly noisy and strongly oscillating phase courses. Starting above a d value of 0.6, coupling leads to severe impact on the phase time courses, to an extent that the overall spectral properties of the time courses are distorted. As this does not reflect a realistic scenario we limited our simulations to d values within the range of [0.01, 0.6].
Gaussian noise was added with intensity levels within the range [0.5, 3] resulting in 5% to 30% noise, which falls into the regime of stochastic resonance in modelling of fMRI time courses (Deco et al., 2009). Fig. 5 shows a typical example of a synthetically generated IBN phase course within freq3. The according procedure is described in Appendix A of the main paper.

Rediscovery analysis
We analyzed in how far the exact couplings parametrically depending on the score s could be recovered by the analysis pipeline. We, therefore, separated couplings into three categories being the couplings actually parametrically manipulated within the closed cluster in the ground truth data (amounting to 156 couplings, upper left part of the matrix in Fig. 6, couplings 'within the patterns' N s and N as respectively, termed IC), those couplings between an IBN from within the closed correlations pattern and an 'outside' IBN (amounting to 182 couplings, bridging couplings termed BC), and those only connecting IBNs from outside the set of manipulated coefficients (amounting to 42 couplings, lower right part of the matrix in Fig. 6, reference couplings termed RC). In this last group of reference couplings the detection of a significant parametric dependence of the K res Figure 5: Typical simulated phase course. The figure shows an (already down sampled) phase course as typically produced by the simulation process described in Appendix A of the main paper. We used the convention to state phase values within the interval [−π, π). Generated values (circles) are linearly interpolated for better visualization. parameters on the score s was regarded as false positive. In the group of 'bridging couplings' such a dependence would highlight at least one region actually affected by true parametric dependence on s. We evaluated the rediscovery of significant clusters with two versions.
Metric 1. The numbers of significant coefficients in the respective regions IC, BC and RC are divided by the total number of significant coefficients, i.e. within the significant set.
Metric 2. The numbers of significant coefficients in the respective regions IC, BC and RC are divided by the total number of significant coefficients, i.e. within the significant set, and it is additionally scaled by the size of the respective region (i.e., 156 for IC, 182 for BC and 42 for RC).
Note that the total number of coefficients in the cluster can vary largely across runs and parameter settings. Metric 1 expresses IC, BC and RC shares of the number of significant coefficients in the cluster and, therefore only accounts for the size of the cluster, but not the sizes of the IC, BC and RC areas. Metric 2 accounts for both. Since metric 1 does not account for the sizes, it gives some insight irrespective of the choice of the set of coefficients with induced correlations (i.e., manipulated coef-ficients), which influences the rediscovery values for IC, BC and RC due to the scaling with their sizes.
We focused on scenarios for which the significance of retrieving a significant coupling set by the analysis pipeline was ≤ 0.05. We, therefore, evaluated the two metrics in two versions. 4 Additional empirical data results

Dependence of K on head movement and cardiac rate
In the exploratory data set, patients with greater disease severity exhibited more motion during the scan (P < 0.034, Spearman correlation). When testing for a dependence of K on head movement, we also found a significant relation ('all': P = 0.028, 'pos': P = 0.014, 'neg': P > 0.05). To analyse the contributions of movement and Ham-D on the sets of Ks we followed the Baron and Kenny's steps for mediation (1986). The direct effect (not mediated by Ham-D) of movement on K was calculated by regressing Ham-D out of the movement parameters and using the resulting regressor in the set-level statistical analysis and vice versa. Results of the mediations are shown in Fig. 7. While Ham-D almost completely mediates the dependence of motion on K, motion does not mediate the dependence of Ham-D on K.
In the test data set a bite bar had been used. There was neither a dependence of K on head movement nor of Ham-D on head movement in the test patient data set.
There was no dependence of K on the cardiac rate in the exploratory or test data set.

Additional mask matrices
In the main paper in Appendix A we presented cluster permutation test results for the two correlation patterns N s and N as and three mask matrices. We performed simulations with three additional mask matrices. All six masks are depicted in Fig. 9. Fig. 10 and 11 show the collective results for all masks for N s and N as , respectively.

Parameters determined by the analysis pipeline deviated from zero to a very limited extent
For all runs of all parameter settings we applied a Wilcoxon signed rank test for all coefficients K res i,j (Bonferroni corrected for multiple testing with 380), if coefficients were significantly different from zero. There was no substantial evidence for the detected coefficients to significantly deviate from zero (see Appendix A of the main paper). Table 2 shows the collective Wilcoxon signed rank test results for the coefficients K res for all correlation patterns and mask scenarios.

Results for random initial matrices
We further present in Fig. 12 the median P -value results for the random initial matrices K o (see Appendix A in the main paper).

Recovered sets of couplings are shifted to joining regions
For each mask and correlation pattern we calculated the two metrics (do not account for area size vs account for size) for IC, BC and RC, respectively and evaluated both versions A and B (consider all significant (n, d) settings across runs vs consider only settings which are significant in the median of runs). Table 3 summarizes the results. For metric 1 we can see for both versions, that for all scenarios BC > IC > RC with three exceptions, where IC > BC > RC. When accounting for the sizes of the areas, the ratio of false positives increases. However, IC and BC combined still contain the most significant coefficients. Furthermore, also the SD increases disproportionately for RC compared to IC and BC, which makes the rations for recovered false positives less stable.
We conclude, the reliability of significant coefficients can be mainly addressed on an IBN level, but not necessarily for the specific coefficient. In general, couplings of group BC had the highest likelihood of being included into a significant set.
In Table 3 all results for both metrics and both evaluation versions are presented for the all correlation pattern and mask settings.  . It only occurred twice that for a single (n, d) combination more than one coefficient resulted in a positive test. This was only the case in run 2 for the setting N s with M 0.5/0 and in run 5 for the setting N as with M 1/0 , for which in one (n, d) combination, 3, respectively 2, coefficients were seen as significantly different from 0 by the test. The cluster permutation test in these cases yielded a significant result, however, the coefficients with a positive signed rank test were all not part of this significant set.  IC 0.38 ± 0.16 0.34 ± 0.08 0.25 ± 0.11 0.22 ± 0.05 M 0.5/0 IC 0.27 ± 0.16 0.25 ± 0.08 0.17 ± 0.10 0.16 ± 0.05 BC 0.51 ± 0.15 0.54 ± 0.05 0.28 ± 0.08 0.30 ± 0.03 BC 0.60 ± 0.16 0.60 ± 0.08 0.33 ± 0.09 0.33 ± 0.05 RC 0.11 ± 0.06 0.12 ± 0.05 0.26 ± 0.14 0.29 ± 0.12 RC 0.14 ± 0.10 0.15 ± 0.07 0.32 ± 0.24 0.36 ± 0.17 M 1/0 IC 0.34 ± 0.19 0.23 ± 0.12 0.22 ± 0.12 0.15 ± 0.08 M 1/0 IC 0.27 ± 0.15 0.27 ± 0.08 0.17 ± 0.10 0.17 ± 0.05 BC 0.52 ± 0.17 0.58 ± 0.11 0.28 ± 0.09 0.32 ± 0.06 BC 0.59 ± 0.17 0.59 ± 0.10 0.32 ± 0.09 0.32 ± 0.05 RC 0.15 ± 0.12 0.19 ± 0.08 0.35 ± 0.28 0.46 ± 0.18 RC 0.14 ± 0.09 0.14 ± 0.05 0.33 ± 0.21 0.33 ± 0.13 M 1/0.5 IC 0.47 ± 0.14 0.47 ± 0.03 0.30 ± 0.09 0.31 ± 0.02 M 1/0.5 IC 0.41 ± 0.21 0.37 ± 0.00 0.26 ± 0.13 0.24 ± 0.00 BC 0.44 ± 0.13 0.44 ± 0.02 0.24 ± 0.07 0.24 ± 0.01 BC 0.49 ± 0.20 0.50 ± 0.00 0.27 ± 0.11 0.28 ± 0.00 RC 0.08 ± 0.05 0.08 ± 0.01 0.19 ± 0.11 0.19 ± 0.03 RC 0.10 ± 0.10 0.12 ± 0.00 0.  Table 3: Retrieved cluster coefficients. For each correlation pattern and mask matrix combination, the mean values and standard deviations for the two metrics (1 and 2) and two evaluation versions (A and B) are given. We can observe for metric 1 where we do not consider the size of the areas IC, BC and RC, that the BC area had, in general, the highest likelihood of containing significant cluster coefficients. For most cases, false positives (i.e., RC) had the lowest likelihood. In specific cases, e.g. N s and M 1/0.5 , the IC had the highest likelihood. * Values of metric 2 are naturally very small since they are not only divided by the total number of coefficients with significant correlations, but additionally by the size of the respective area. For better readability and comparability these values were multiplied by 100. Here, we can observe, that the proportion of false positives increases in comparison to the other areas, especially as significance (low P -values) in the median is more prominent (see Fig. 10 and 11). However, SD also increases disproportionately for RC. Further, IC and BC combined remain the largest portion of recovered coefficients. combinations, K c with the correlation pattern N s . Each subplot shows the results for one of the 6 mask matrices. Cells with significant P -values are colored in black. It shows that increasing d leads to a better capacity of the analysis pipeline to retrieve the significant dependence on the sore s from the artificial data, when the contrast in dependence of the K c coefficients outside the set of coefficients with strong correlations are set to 0 ((b) and (c)) or have very weak influence (f). When the contrast in dependence of the K c coefficients outside the set of coefficients with strong correlations to within this set is weak to moderate ((d) and (e)), or even non present (a) significance is only retrieved randomly. Added noise ε in the coupling process has a detrimental effect on recovering the dependence on s.  Figure 11: Results for N as . The above figures show the median P -values over the 6 runs for all (n, d) combinations, K c with the correlation pattern N as . Each subplot shows the results for one of the 6 mask matrices. Cells with significant P -values are colored in black. It shows that increasing d leads to a better capacity of the analysis pipeline to retrieve the significant dependence on the sore s from the artificial data, when the contrast in dependence of the K c outside the set of coefficients with induced strong correlations are set to 0 ((b) and (c)) or have very weak (f) to moderate (e) influence. When the contrast in dependence of the K c coefficients outside the set of coefficients with strong correlations to within this set is weak (d) or non present (a) significance is only retrieved randomly. Added noise ε in the coupling process has a detrimental effect on recovering the dependence on s. Significant effects can be retrieved more prominently for this correlation pattern than for the symmetric correlation pattern N s (see Fig. 10).
(a) Results for K o and M 1/1 . The median P -values show no dependence on (n, d). significant settings only appear randomly.
(b) Results for K o and M 1/0 . The median P -values show no dependence on (n, d) even if a large portion of coefficients is excluded in the simulation process. Figure 12: Results on random initial matrices. Cells with significant P -values are colored in black.
(a) If we simulate phase courses where the coefficients do not show a significant correlation with the the score s (i.e., K o and mask M 1/1 ), the significance in the statistical analysis shows a random behavior. Some parameter combinations indicate the finding of a network cluster associated with s, but this does not correspond to the intensity of the coefficients in the simulation or the noise level. In general, fewer and less significant results are found. (b) To determine, if the results are associated with the correlation structure rather than with the applied masks M 1/0 or M 0.5/0 , we repeat the experiment for the original random correlated matrices K o , but set the same coefficients to zero as in the cases of N s and N as . Again the results show a random behavior, which indicates a strong impact from the correlations and a weaker impact of the exclusion of certain coefficients by setting them to zero on the results.