## 2 Materials and Methods
### Reproducibility between Measurements: Characteristics and Statistical Tests
To check this, we computed and plotted matrices of Spearman correlations between comodulograms from the first half-session of every subject with the second half-session of every other subject. In the event of the existence of reproducible and unique patterns in the comodulograms of each subject, these matrices should have large diagonal values (within-subject correlations) and smaller off-diagonal values (between-subject correlations). The reproducibility between the half-sessions of the phase and amplitude frequencies of maximal TGC was characterized by two quantities: the mean absolute difference (change) and the intraclass correlation coefficients (ICC). The mean difference characterizes the absolute variability in Hz, while the ICC measures the reliability relative to the between-subject variance. We calculated ICC (C,1) in the terminology of [61] and tested the hypothesis \"ICC = 0\" using the MATLAB package \"Intraclass Correlation Coefficient (ICC)\" [62]. This test determines whether there is a significant subject effect on the frequencies (the \(p\)-values are actually the same as for the subject effect in a two-way ANOVA with factors: subject and half-session). Thus, this test shows if there are reliably inferred nonzero differences between mean frequencies for different subjects. This is important in the context of brain stimulation for the individualization procedures to be meaningful. To compare the methods in terms of absolute variability of peak frequencies, we performed [PERSON]'s tests for the effect of the algorithm on the absolute frequency changes, followed by post-hoc tests using the Tukey-Kramer method. To test for a possible effect of the half-session (due to fatigue, learning, or other causes), we performed [PERSON]'s signed-rank tests on the peak frequencies from the two half-sessions for each PAC calculation method. To illustrate the degree of reproducibility of the comodulograms, we included their plots and scatter plots between half-sessions for a single subject as an example. We selected this subject to have average peak variability. To this end, we calculated a composite measure of variability derived from both the phase and amplitude peak frequencies. This was performed by dividing the absolute frequency change by the length of the corresponding frequency range and averaging the resulting two values:  
\[\Delta f_{composite}=\frac{1}{2}\bigg{(}\frac{\Delta f_{p}}{r_{p}}+\frac{ \Delta f_{a}}{r_{a}}\bigg{)}\]  
where \(\Delta f_{p}\) is the absolute change between half-sessions of the phase frequency of maximum PAC, and similarly for \(\Delta f_{a}\); \(r_{p}\) and \(r_{a}\) are the lengths of the theta and gamma frequency ranges (\(8-4=4\) Hz and \(60-30=30\) Hz, respectively). Thus, this composite measure indicates by what fraction of the whole frequency interval the peak has shifted. The subject for illustrations was chosen to have the composite variability (averaged by all the algorithms) closest to the group median. We also compared the algorithms in terms of the central tendency of the frequencies of maximal TGC to see if any of the methods produced consistently lower or higher peak frequencies. For this, we performed [PERSON]'s tests separately for the phase and amplitude peak frequencies, followed by post-hoc tests using the Tukey-Kramer method. This analysis was not concerned with reproducibility, and hence we performed it for comodulograms computed from whole sessions rather than half-sessions.