Functional Connectivity States of Alpha Rhythm Sources in the Human Cortex at Rest: Implications for Real-Time Brain State Dependent EEG-TMS

Alpha is the predominant rhythm of the human electroencephalogram, but its function, multiple generators and functional coupling patterns are still relatively unknown. In this regard, alpha connectivity patterns can change between different cortical generators depending on the status of the brain. Therefore, in the light of the communication through coherence framework, an alpha functional network depends on the functional coupling patterns in a determined state. This notion has a relevance for brain-state dependent EEG-TMS because, beyond the local state, a network connectivity overview at rest could provide further and more comprehensive information for the definition of ‘instantaneous state’ at the stimulation moment, rather than just the local state around the stimulation site. For this reason, we studied functional coupling at rest in 203 healthy subjects with MEG data. Sensor signals were source localized and connectivity was studied at the Individual Alpha Frequency (IAF) between three different cortical areas (occipital, parietal and prefrontal). Two different and complementary phase-coherence metrices were used. Our results show a consistent connectivity between parietal and prefrontal regions whereas occipito-prefrontal connectivity is less marked and occipito-parietal connectivity is extremely low, despite physical closeness. We consider our results a relevant add-on for informed, individualized real-time brain state dependent stimulation, with possible contributions to novel, personalized non-invasive therapeutic approaches.


Introduction
The origin of alpha waves and the function they subserve constitute long-lasting scientific issues in neuroscience. Already by 1929, Berger had managed to isolate alpha waves by means of a pioneering EEG set-up using scalp electrodes and described this rhythm as the most prominent in the human electroencephalogram [1]. Recent advances due to brain source reconstruction and invasive electrophysiological recordings have provided evidence that alpha waves originate from several cortical and subcortical sites, with direct evidence suggesting both thalamus and diverse cortical areas as possible origins of such rhythm [2][3][4].
The functional role of alpha oscillations also is still far from being completely understood. The "pulsed inhibition hypothesis", for instance, proposes that alpha oscillations actively inhibit neuronal firing in a phasic manner, opening and closing interleaved periods of "high" and "low" excitability of the cortex by cyclically producing bouts of inhibition [5,6]. This hypothesis is in line with the idea behind brain-state dependent stimulation, that the outcome of an electro/magnetic perturbation depends on the instantaneous phase state of a specific brain rhythm in a given area. In fact, studies investigating alpha in the occipital cortex demonstrate that the alpha phase at the instant of a visual stimulus predicts high or low probability of detection [7].
However, it is still not clear which is the exact role of alpha oscillations in opening states of cortical excitability. For example, the aforementioned "pulsed inhibition hypothesis" does not seem to apply to findings of brain-state dependent stimulation in the parietal cortex, where the same directionality of inhibition vs. excitation as in the visual cortex does not emerge [8][9][10]. In this regard, studies investigating µ-alpha in the parietal cortex show that a principle of pulsed facilitation rather than inhibition would better explain the role of alpha in the hand-knob of the sensorimotor cortex. In fact, it has been shown by means of real-time EEG phase-dependent transcranial magnetic stimulation (TMS) that stimulation at µ-alpha troughs results in facilitated motor evoked responses (MEPs, [8]), while at peaks of the same oscillation, inhibition does not seem to occur in a significant proportion of instances [9].
Even if there is no consensus as to whether alpha shapes neuronal recruitment by determining windows of lower and higher excitability, or rather only by opening windows of greater excitability, alpha oscillations appear to have a role in shaping neural recruitment.
Most of the protocols conceived to modulate cortical excitability by means of TMS use a predefined stimulus sequence irrespective of the instantaneous brain-state, as opposed to a real-time brain-state dependent stimulation, which delivers the stimulus in determined phases of the alpha cycle [8]. If the phase of alpha reflects a phasic increase and decrease of cortical excitability (as reported in the occipital cortex by [11], we should find different outcomes depending on the local alpha phase at the instant of stimulation. First attempts in this direction have tried to post-hoc determine the phase of alpha waves at the moment of the stimulus. This has been mainly tested in the occipital cortex, where the conscious visual percept has been linked post-hoc to the alpha phase at the instant of stimulus presentation [7,12]). However, the effects were generally assessed by statistically estimating the probability that the stimulus could be delivered at a given phase. Differently, phasestate dependent stimulation leverages instantaneous EEG phase-states to trigger the TMS pulse exactly at the phase of interest, without post-hoc tracing of it back to the moment of analysis. Therefore, studies using this novel approach have led to more consistent and reproducible results. For example, EEG-TMS has been used to investigate motor excitability depending on mu-alpha in the parietal cortex. This line of research is relatively recent, but several consistent pieces of evidence show that participants with a detectable mu-alpha rhythm show larger motor evoked potentials (MEPs) when TMS stimulation at the motor spot is delivered at mu troughs or at the early rising phase, compared to the conditions where positive peaks or random phases are targeted [8,13,14]. These results are in line with the hypothesis of the alpha rhythm as a mechanism modulating cortical excitability. Attempts to link specific alpha phases to enhanced cortical excitability have also been made in studies targeting the dorsolateral frontal cortex (DLPFC), where alpha-synchronized rTMS at troughs appears to provide for a local alpha power decrease in patients with drug resistant depression disorder. The same result is not obtained by intermittent theta burst or random phase stimulation [15].
To sum up, the role of alpha phase has been investigated in at least three different cortical areas and a relationship between alpha phase and cortical excitability has been found in frontal, parietal and occipital regions [8,11,16]. However, there is still no clear and organic consensus on a generalized functional role of the alpha phase. For example, the alpha phase in the occipital cortex seems to elicit opposite effects with respect to the other cortical areas: trough targeting has been proposed to decrease the possibility of perceiving the stimulus [7]. Moreover, it has been shown that the conscious perception of a visual stimulus not only depends on alpha phase in the occipital cortex but also on that in fronto-central regions at the moment of stimulation [7,12] this would be in line with the idea of the frontal control network allowing access to a conscious perception. This evidence, together with the assumption that coherence between brain regions underlies integration of information [17], suggests that knowing the functional coupling between stimulated (or post-hoc investigated) regions and other areas potentially connected with the stimulation site is crucial in designing experiments aiming at a trial-based selective perturbation of a given "brain state". In this light, considering the functional coupling patterns of the stimulated site is important in explaining EEG-TMS results and is even more crucial when EEG-TMS protocols are intended, not only in terms of a single stimulation site, but also as a technique for pathway-specific modulation targeting multiple functional hubs: the next horizon of brain-state dependent stimulation. In this regard, this technique has been shown to have the potential to modulate specific pathways [18], for example, by coupling the stimulation to an activity state of that pathway (e.g., [13]). Therefore, reliable metrics are required that can extract extended network states through long-range connectivity. Such solid brain state landmarks would be certainly useful (before and after a neuro-modulatory intervention) to assess whether the brain-stimulation has exerted the desired network state modulation. Most importantly, however, they would be even more crucial when time-resolved connectivity state estimates, describing the connectivity pattern of the network, are used as real-time trigger condition. Nevertheless, gold-standard metrics for connectivity generally consist of pseudo-statistics regarding different phase consistency over several tens of trials and a consensus for a conceptual definition of singletrial instantaneous connectivity is still missing. For this reason, we here try to open the way to addressing the need for brain-state dependent protocols to take into account not only the local state which usually triggers the stimulus, but also a general sensitivity state of the system being modulated in its excitability.
Therefore, this study aims to identify suitable functional pathways between wellknown cortical alpha generators, whose connectivity state can be assessed with MEG/EEG at rest. Here, we present a pipeline that effectively determines potential connections from resting MEG/EEG data using Weighted Pairwise Phase Consistency (WPPC; [19]) and Weighted Phase Lag Index (WPLI; [20]). It is worth noticing that the choice of a connectivity metric has to take into account its advantages and disadvantages in the context of application [21,22]. Here we used WPPC and WPLI, exploiting their partial complementarity, for the data under investigation. Both metrics depend on the consistency between the phases of the signal of interest, and are not biased by sample size. WPPC is based on the pairwise phase difference, thus it is still affected by zero-lag correlations introduced by the spatial spread of the inverse solution. As a matter of fact, WPPC also provides positive and statistically significant values when the phase difference is exactly zero. This is a result that does not reflect real synchronization, because the finite (but not null) propagation time of the nervous signal on the pathway connecting the two areas is not taken into account. In contrast, WPLI does not have this problem because it is computed from the sign of the imaginary part of the coherence, which vanishes for zero-lag correlations. However, the WPLI signal to noise ratio is optimal for a phase delay of π/2, which corresponds at a typical frequency of 10 Hz to an absolute time delay of 25 ms. Since we are investigating the temporally resolved cortico-cortical connectivity, this is a large delay if compared to the typical brain conduction delays within the same hemisphere. Thus WPLI, despite being a less sensitive measure of phase consistency between the brain regions under investigation, is useful to confirm that the connectivity already detected by the more sensitive WPPC is real and not due to the spatial spread of the inverse solution. In this sense, the two metrices are partially complementary. With this methodological set-up, we investigated 203 healthy subjects with several minutes of resting MEG data and found alpha connectivity stronger for parietal-prefrontal areas rather than for occipito-parietal areas. The results of this study will be relevant to design real-time brain state-dependent EEG-TMS experiments: connectivity priors [23] will be highly relevant as off-line acquired a priori information for the development of real-time connectivity state estimation algorithms.

Dataset and Acquisition
We analyzed magnetoencephalographic (MEG) resting state recordings of 203 healthy participants (age range 18-57 years; see Figure 1a for age distribution) from the Cambridge Centre for Aging and Neuroscience (CamCAN) public dataset. Prior to inclusion in the dataset for neuroimaging measurements, all participants were tested for cognitive decline (MMSE > 24; [24]), for matching vision, hearing and English language inclusion criteria and for absence of serious neurological and psychiatric conditions (see [25] for details about the dataset and acquisition parameters). For each participant, resting state activity (eyes closed; 8 min and 40 s) and empty-room noise background (3 min) were recorded using a 306-channel VectorView MEG system (102 magnetometers; 204 first order planar gradiometers; sampling rate = 1000 Hz; high pass filter 0.03 Hz; low pass filter 330 Hz). Anatomical landmarks (nasion, left and right pre-auricular points) were registered, as well as at least 75 additional (isotrak) points, to model the head surface. For human resting state data, continuous monitoring of the head position (cHPI), electro-ocular (EOHG, EOVG) and electro-cardiac (ECG) recordings were available. Additionally, T1 weighted anatomical data were collected using a 3T Siemens TIM Trio scanner (MPRAGE, TR = 2250 ms, TE = 2.99 s; FOV = 256 × 256 × 192; voxel size 1 × 1 × 1 mm). All the subsequent analysis was performed only on planar gradiometers using Fieldtrip ( [26]), SPM ( [27]), CAT12 (http: //www.neuro.uni-jena.de/cat/ (accessed on 1 October 2021)) and custom MATLAB code. (c) Distribution of all condition numbers from SVD, aiming to define optimal dipole direction for spectral data.

MEG Data Preprocessing
Collected MEG datasets (resting state and empty-room) were inspected and bad channels showing high noise levels and/or SQUID jumps were marked and removed from the data. External magnetic source nuisance was removed by means of temporal Signal Space Separation (tSSS, [28,29]) with a correlation threshold of 0.98 and a sliding time window of 10 s. Line noise contribution was removed by applying a 5th-order Butterworth two-pass filter centered at the line frequency and its first 3 harmonics, and spanning an interval of 2 Hz. Head position of the resting state data was corrected every 200 ms using cHPI data and rereferenced to a common head position. Human MEG data were cleaned from physiological artifacts using an automated procedure. Muscular activity was detected by filtering data between 100 and 140 Hz, transforming each channel data into a z-score relative to the whole channel time series and rejecting segments where the z-score averaged across channels was greater than 5 for at least 200 ms. Then we ran an extended Infomax Independent Component Analysis (ICA; [30,31]) on data filtered between 0.5 and 125 Hz (Butterworth 4th order, two-pass) and resampled to 250 Hz. Data resampling and digital filtering commute with the ICA unmixing matrix estimation, provided the artifactual activity of interest lies in the spared frequency band [31]. For this reason, and for computational efficiency, we computed ICA on resampled data and applied unmixing matrices to unfiltered data at the original sampling rate of 1000 Hz. Correlation coefficients of each extracted component with electrophysiological channels (EHOG, EVOG and ECG) were computed and rejected using a recursive z-score based procedure that robustly discarded all components with correlation more than 2 standard deviations from the coefficient distribution mean of all components. Data cleaned from artifacts were then filtered and segmented. Empty-room data, later used for noise covariance matrix estimation, were filtered broadly between 0.5 and 140 Hz (4th-order Butterworth two-pass filter) and segmented in 1 s epochs. Human resting state recordings were instead filtered around the alpha frequency band (4th-order Butterworth two-pass filter between 5 and 16 Hz) and split into 2 s segments. Finally, all empty-room and resting state segments exceeding a threshold of 10 standard deviations with respect to the channel average, were further discarded in order to deal with potential residual SQUID jumps and/or filter border effects.

Anatomical Data Processing
Structural T1w MRI scans were processed with the aim of extracting cortical surfaces for forward and inverse model calculation and for common space mapping. We processed T1w images using the CAT 12 toolbox pipeline (http://www.neuro.uni-jena.de/cat/ (accessed on 1 October 2021)). After bias normalization, denoising and skull stripping, volumetric data were automatically co-registered to a template space (IXII 555 MNI space; www.brain-development.org, (accessed on 1 October 2021)) Different tissue was then segmented in native space in order to extract surface models of the head, the brain enclosing surface and the cortical mantle. We modeled the cortical mantle with a tessellation (20,484 vertices) of the mid thickness surface, i.e., the surface in between the pial and the white/gray matter interface. The surface enclosing the brain and a surface model of the head were also modeled as 20,000 vertices meshes. In addition to the extracted surfaces in native subject space, co-registered spheres for projection on the FreeSurfer Average (5th order icosahedron "fsaverage5"; [32]) superficial template space were computed, as well as the corresponding interpolation matrices between template and native coordinate spaces. Interpolation coefficients for each point were defined as inversely weighted average of first neighborhood. The mapping of the native space cortical surface onto the FreeSurfer average allowed us to identify, in template space, vertices belonging to 360 regions of interests of the 'state of the art' multimodal parcellation from the Human Connectome Project [33]. This atlas was generated by combining structural, diffusion and resting state fMRI data from 210 healthy young adults. Being interested in phase consistency between frontal, parietal and occipital areas, we defined three corresponding brain sectors, for each hemisphere, by pooling together correspondent ROIs from the atlas. We selected the ROIs in order to achieve a sufficiently large coverage of the three sectors of interest, while avoiding excessive proximity that might lead to spurious connectivity due to potential spread of the inverse solution. Finally, this led to 16, 7 and 18 regions of interest for the occipital, parietal and frontal sectors, respectively. For the sake of computational efficiency, we refined the sectors by trimming the borders in order to have the same number of dipoles per sector (D = 930). A list of corresponding regions of interest, with MNI coordinates, can be found in Table 1. All the subsequent analyses were performed only on these ROIs.

MEG Source Reconstruction Based on Individual Anatomies
We used Minimum Norm Estimation [34] to solve the inverse problem and compute the projection matrix from sensors to source space. Native space surfaces were co-registered to MEG coordinates in a first step, using correspondence between anatomical landmarks as recorded during the MEG session and in the MRI anatomical image. Second, co-registration was refined by aligning additional head surface points registered during MEG acquisition to the tessellation representing the head surface. In this procedure, MEG isotrak points anterior to the nasion were discarded since the anatomical MRI images were defaced prior to being publicly available. We then computed a forward model solution for planar gradiometers using the singleshell method and the brain enclosing surface as computed before, and depth normalizing lead fields with a factor of 0.5. The source model for MNE was defined as a free orientation set of dipoles uniformly distributed on the cortical surface (mean spacing 3.1 mm) and the inverse solution was computed, with a 1% noise covariance regularization. This procedure leads to a set of three time series of estimated cortical activations for each resting state epoch and for each vertex of the source model mesh, in the three cartesian coordinate system. MEG is blind to the radial component of magnetic field generated by a current dipole in the source model [34,35]. Thus, even for realistic forward solutions, the estimated current in the most radial direction, with respect to the head surface, is almost zero. Therefore, we performed a Singular Value Decomposition (SVD) of the three source time series at each vertex and for each epoch, retaining only the components associated with the first two singular values, while the third was always almost zero. In this way we reduced the current estimate to its projection on the plane orthogonal to the radial direction, resulting in two source activity time series for each vertex and for each dipole represented hereafter by the two-dimensional vectors x(t).

Spectral Analysis
For each epoch and vertex of the cortical mesh, we computed Fourier coefficients X(f) from the time dependent activity vectors x(t) in a frequency interval f = [8:14] Hz using a Hanning window taper. Having thus two Fourier coefficients at each frequency, dipole and epoch, we decided to reduce spectral data defining an optimal dipole orientation as follows: we computed the cross spectral density matrix between the two Fourier components, performed an SVD and keeping only the direction relative to the first singular value. This resulted in an optimal dipole orientationû that depends on the dipole position, the epoch and the frequency of interest, thus optimizing the detection of the brain signal at the frequency of interest. We reduced accordingly the Fourier coefficient vector to a scalar one defined as X( f ) = X( f ) ·û. As a measure of the quality of the procedure, we pooled together all condition numbers resulting from SVDs on all subjects: the average condition number was 1.5 × 10 3 . This means that, on average, the dipole orientations were, in time, almost fixed (see also [36]) and thus our optimization procedure captures most of the spectral content of the data. The distribution of the logarithm of all condition numbers can be inspected in Figure 1c, showing the reliability of the assumption. Fourier coefficients in the optimized direction X(f) were then used to compute power spectral density at each vertex and for each epoch. Pooling together all spectral densities and detecting the peak in the frequency band of interest, we defined, for each subject, an Individual Alpha Frequency (IAF); a distribution of all the 203 IAFs is reported in Figure 1b. All the subsequent connectivity analysis was then performed at the individual alpha frequency. For this reason, the frequency dependence of Fourier coefficients X(f) will be dropped hereafter in the notation.

Connectivity Analysis and Group Statistical Validation
We are interested in connectivity and phase relationships between region of interest belonging to the frontal, parietal and occipital areas within each hemisphere. For this reason, we computed two different spectral based connectivity metrics between all the combinations of areas belonging to the three different sectors. In particular, given the symmetricity of the connectivity metrics we use, we crossed occipital areas with parietal and frontal, and parietal region of interests with frontal ones. For each resulting combination, we computed the Weighted Pairwise Phase Consistency (WPPC; [20]) and the Weighted Phase Lag Index (WPLI; [19]). The WPPC is based on the distribution of the phase differences between all the pairs of observations (resting state 2 s epochs in our case). WPPC is a robust non-biased measure of phase consistency of brain signals, but it is still affected by zero-lag artificial correlations induced by the imperfection of the inverse model solution [21,22]. Complementary to WPPC, we estimated WPLI as a connectivity measure not affected by the artificial zero-lag correlations. Being based on the imaginary part of the coherence, WPLI is immune to zero-lag connectivity but it has another disadvantage: it achieves maximal Signal to Noise Ratio (SNR) when the two brain signals are in a π/2 relationship. For the frequency band of interest (around 10 Hz) this means a time delay of~25 ms, a long time when compared to typical brain signal propagation time. However, both metrics strongly depend on a consistent phase relationship between signal of interest, and then they can complementarily provide information about the phase consistency of alpha rhythmicity between the sectors we are investigating.
Given the combination of two brain sectors (A and B) from the ones defined above (O = occipital; P = parietal; F = frontal) we computed whole sector connectivity matrices, using both WPLI and WPPC, as follows. Named for brevity as X = X A (d a ) and Y = Y B (d B ), the Fourier coefficients at each vertex d a and d b of the sectors A and B, respectively, we computed a regional dipole-wise connectivity matrix C M AB (d a , d b ) ∈ M(D × D), where M represents the metrics (M = WPPC or WPLI). Henceforth, for the sake of clarity, subscripts and superscripts will be omitted when not necessary. In addition, we computed, for each metrics, correspondent null connectivity matrices C A (d a , d r ) and C B (d b , d r ) from sector of interest to a set of D dipoles {d r } randomly chosen within each subject space among the set of dipoles not belonging to any sector of interest. These null connectivity matrices were then used for group statistical validation, under the null hypothesis that connectivity to random dipoles, differently chosen for each subject, will be distributed according only to the metrics' bias and sensitivity. To this aim we performed a bootstrap procedure at the group level by comparing actual and random connectivity matrices: for each combination (A, B) and for each subject, the actual connectivity matrix C AB and the random ones C A and C B were permuted 10,000 times in order to estimate the empirical distribution of the Fisher regularized difference of connectivity dC ≡ tan −1 (C AB ) − tan −1 C A/B . The comparison between the real dC with its null empirical distribution provided, for each dipole, a bootstrap p-value. (only positive tailed comparisons were considered). Resulting D 2 p-values were then FDR corrected (q = 0.05) and only significant elements of C_AB were considered in the subsequent analysis. Furthermore, we reduced the information in each connectivity matrix by summarizing connectivity between the two sets of ROIs {r a } and {r b } from the atlas [33] belonging to sectors A and B, respectively. To this aim we defined, for all connectivity metrics, the following two quantities: where # [·] and fˆsg [·] represents the count and the distribution of significant connectivity values, respectively, while P_95 [·] represents the 95% percentile. The quantities Γ and ∆ were computed between regions for each combination {r_a,r_b} and from a single region of interest to all the target sector {r_a,B}. Finally, bi-hemispheric results were collapsed by averaging the contribution of both hemispheres. We defined the quantities ∆ and Γ to summarize the connectivity between ROIs and/or sectors, given that inspecting he whole dipole by dipole connectivity matrices would have been confusing and not clear to the reader. The two quantities are conceptually derived from standard graph theory analysis. The connectivity degree ∆(r_a,r_b) simply counts, from the connectivity matrix, the number of statistically significant connections between dipoles of the two ROIs, disregarding the strength of the connectivity and normalizing the count to the total number of dipoles in the sector. This is conceptually analogous to the node degree in the context of network analysis. The mean significant connectivity Γ(r_a,r_b) provides a refinement of the connectivity degree, including the information about the connectivity (or phase consistency) strength. This is achieved by selecting from the connectivity matrix only the statistically significant connections between the dipoles of the two ROIs, extracting the 95% percentile of the resulting distribution instead of just counting them. It is worth noticing that choosing the 95% percentile is a conservative choice, given that usually connectivity values bounded between 0 and 1, even when Fisher regularized, can give rise to non-symmetrical distributions.

Results
Degree and mean connectivity between each ROI combination (∆(r a , r b ) and Γ(r a , r b )) and from each ROI to the other sectors as a whole (∆(r a , B) and Γ(r a , B)) for WPPC and WPLI are shown in Figures 2 and 3 as color coded source maps and connecto-grams. All extracted values are listed in Tables 2 and 3. In general, we can notice that the connectivity estimated from both metrics, each one differently depending on a consistent phase relationship, is predominant in the parietal-frontal connection. While the connection between occipital and frontal sectors appears to be still relevant, the weakest phase consistency has been found between the occipital and the parietal areas. The predominance of phase consistency in parieto-frontal connections is also confirmed by the value of Γ(P, F) = 0.06, to be compared with Γ(O, P) = 0.03 and Γ(O, F) = 0.03 obtained by WPPC considering the whole area as a single ROI in the analysis. When considering the same values but extracted from WPLI, the pattern is less evident (Γ(O, P) = 0.012; Γ(O, F) = 0.015 and Γ(P, F) = 0.013), since the mean connectivity between sectors as detected by the imaginary part of the coherence is almost the same. However, this has to be interpreted in the light of the different sensitivity of the two metrics with respect to the absolute value of phase shift at IAF between the source activity in the two areas, which is ultimately related to the propagation time of the nervous signal on the pathway connecting the ROIs under consideration. As a matter of fact, comparing values of connectivity parameters as extracted from WPPC and WPLI, we can see that the former are in general larger than the first ones. The pattern of a predominant fronto-parietal connectivity with a less consistent pareto-occipital connectivity is also less marked in the WPLI case. We interpret this, methodologically, by considering the different sensitivity of WPPC and WPLI with respect to phase delay absolute values. Being based on the imaginary part of the coherence, WPLI is more sensitive to phase consistent signals whose time delay corresponds to values close to φ = π/2, which translate in the alpha band into a propagation time delay of signals of~25 ms. This time is relatively large if compared to typical brain conduction delays, therefore we expect the phase delay value giving rise to significant connectivity measures to be more toward a phase difference of φ = 0, thus expecting WPPC to be more sensitive. However, because of the finite conduction delays in the white matter, and in general in the brain, we expect more distant areas having larger propagation time. Thus, we expect phase consistency between distant areas to be more evident in the WPLI.  Γ(r a , B). In the inset the color code is explained: assigning red/green/blue colors to occipital/parietal/frontal sectors, respectively, the intensity of each area represents the overall magnitude of the parameter of interest, while the color encodes the proportion of the magnitude due to connectivity to the other sectors. So, for example, a parietal area more towards blue has more consistent phase relationship to the frontal sectors, while, if more towards red, the most relevant connection is to the occipital sector.  Γ(r a , B). In the inset the color code is explained: assigning red/green/blue colors to occipital/parietal/frontal sectors, respectively, the intensity of each area represents the overall magnitude of the parameter of interest, while the color encodes the proportion of the magnitude due to connectivity to the other sectors. So, for example, a parietal area more towards blue has more consistent phase relationship to the frontal sectors, while, if more towards red, the most relevant connection is to the occipital sector.  Inspecting more in detail ∆(r a , B) and Γ(r a , B)), as we can see from Figure 2a,b, the overall higher values for both the quantities extracted from WPPC correspond to the paretofrontal connections, giving frontal area 8C the overall highest values (∆(8C, P) = 0.1 and Γ(8C, P) = 0.06)). This predominance can be also appreciated in the connecto-grams (Figure 2c), where connections between frontal areas 8C and 8av appear to be the highest in this sector's combination and with respect to the other sector combinations. This result is further confirmed by values obtained from WPLI, as shown in Figure 3. In this case also, referring to the parietal-frontal connection, area 8C is the most connected, with ∆(8C, P) = 0.0009 and Γ(8C, P) = 0.013. We can thus conclude that this pattern is not due to the proximity of area 8C with respect to the frontal sector, WPLI being insensitive to zero-lag connections by design. However, in this case, the degree and the mean connectivity of areas 8c and 8av are not the highest with respect to other sector combinations, as can be appreciated in the connecto-grams of Figure 3c and in the color-coded source maps of mean connectivity in Figure 3b.
As regards the occipito-frontal connections, the highest WPPC connectivity and degree values has been found from frontal area 9a to the occipital sector, being ∆(9a, O) = 0.  Figure 3c for WPLI, where, in the latter case, connectivity from frontal areas to occipital is more evident than in the WPPC case. As we pointed out before, we interpret this as a byproduct of the different sensitivity of WPPC and WPLI to coherent signals whose phase delay is more towards φ = 0 or φ = π/2. Finally, also in the single ROI connectivity values inspection, the occipito-parietal connectivity seems to be less relevant. This can be appreciated in the connecto-grams for WPPC in Figure 2c, where the highest values, in this case, can be found in frontal area 2 (Γ(2, O) = 0.034 and in occipital area V6 (Γ(V6, F) = 0.033) for the WPPC. This is confirmed also from WPLI results (Γ(2, O) = 0.013 and (Γ(V6, F) = 0.012). However, in this case the predominance of these two areas in the occipito-parietal connectivity is less marked as can be noticed from the connecto-grams in Figure 3c.

Discussion
Our aim was to provide the brain-state dependent EEG-TMS community with a work describing functional phase-dependent relationships between frontal, parietal and occipital areas. In these regions, alpha phase-dependent cortical excitability has been studied both in healthy controls and patients via brain-state dependent stimulation [8,[13][14][15]. (Since the attempts at manipulating the effects of stimulation triggered by the phase detected on a region different from the target have no clear effectiveness [8], we believe that taking into account off-line phase-dependent connectivity patterns between the areas of interest would be of substantial help in predicting results from these studies. Furthermore, consideration of the general phase-dependent connectivity state of the target region, (either trigger-region or not), could be essential for explaining changes in connectivity patterns for the time after application of plasticity protocols up to the moment where the behavioral measure has returned to baseline. Considering not only the TMS trigger-state (phase of the frequency of interest) but also the general connectivity state of the brain between the regions of interest would in fact be useful both for predicting plasticity protocols' efficacy and for interpreting results. For these reasons, we investigated spectral-resolved connectivity metrics that depend on a consistent phase relationship at IAF between frontal, parietal and occipital regions at rest. Results showed greater functional coupling between frontal and parietal areas, compared to very low functional coupling values between frontal and occipital and parietal and occipital pairs ( Figure 2). These results are consistent with previous evidence of a fronto-parietal network at rest (e.g., [37]). We used MEG data and three different phase-coherence measures to show that a clear coupling emerges at IAF between frontal and parietal parcels at rest.
The ability to assess the strength of frontoparietal coupling from resting state data before and after an intervention makes this a suitable pathway for the investigation of pathway-specific plasticity. The alpha frequency band is relevant because this is the frequency band where phase-specific modulation of cortical excitability was previously demonstrated both in the motor system [8,13,14] and in the prefrontal cortex [15]. Therefore, in frontal and sensorimotor regions, state-dependent stimulation based on the phase of alpha frequency appears effective in modulating cortical excitability at rest. However, when TMS is applied to the motor cortex in response to the phase of alpha recorded from the occipital cortex at rest, no modulation of corticospinal excitability was found [8].
For future development of real-time estimators of instantaneous connectivity states, we propose that connectivity metrics which depend on phase relationship between areas, and therefore describe whether two regions share trial-averaged functional coupling, can provide off-line connectivity priors in order to enable more accurate time-resolved connectivity estimates. This broad state of functional coupling will, in fact, constitute a novel level of a priori information to assess whether real-time EEG-TMS paradigms will be effective in modulating cortical excitability in one area when the trigger is recorded from another area. In this light, for example, our finding of the limited functional coupling at rest between occipital and parietal areas can help elaborate on the reason why no effect was found when the hand knob of the motor cortex was stimulated on the basis of the occipital alpha phase [8].
Therefore, we propose that, in order to improve the efficacy of brain-state dependent stimulation protocols, the definition of "brain state" should not only refer to the instantaneous phase of the considered frequency which is triggering the stimulus, but also to the predisposition of the network to be globally perturbed by a stimulus triggered by a given phase of a determined frequency. In this regard, measures of functional coupling of regions at rest or during tasks that involve different functional coupling patterns would allow the obtaining of "off-line" connectivity priors.
In essence, beside the targeting of a local oscillation phase, brain-state dependent stimulation might leverage the off-line acquired, a priori knowledge regarding functional patterns in order to causally test the functional coupling between the brain regions involved. As an example of potential relevance of this approach, it has been shown that a network of prefrontal and occipito-parietal areas is involved in visual target detection and that alpha phase in this network is essential for a correct visuospatial processing [38,39]. Standard brain-state dependent stimulation could be used in order to further test the relationship between these cortical areas. In fact, brain-state dependent stimulation could potentially modulate cortical excitability of one area given the phase of alpha from another region. However, according to our approach, this brain-state dependent stimulation could lead to more effective results only when the two regions are embedded in a functional pattern that describes a state of sensitivity to be modulated [13,21]. A relevant part of cortical excitability could derive therefore from functional coupling of the target areas at rest (or during a task) at a given frequency.
This new perspective would also be useful for developing new and more effective individualized approaches for rehabilitation via brain-state dependent stimulation. In this regard, individualized stimulation rehabilitative protocols have been proposed in stroke patients for whom specific adaptive connectivity patterns between specific cortical regions should be reinforced in order to achieve recovery [40]. In this light, our approach combining brain-state dependent stimulation with phase-dependent connectivity measures would nicely suit this purpose.
Finally, it must be noticed that our work is the first attempt at measuring coherence at IAF during rest between these three regions of the cerebral cortex. Most of the resting state literature consists of fMRI studies, which lack the temporal resolution that MEG measures provide [41][42][43]. It is also worth noticing that, depending on the connectivity metric employed, we obtain slightly different results in the coupling between frontal and occipital areas, with the WPLI showing some functional coupling between the two regions compared to other measures. The WPLI is a measure based on the imaginary part of coherence. This metric is not sensitive to zero-lag spurious correlations. However, its maximum SNR spuriously coincides with π/4. Since we are investigating IAF which, on average, refers to around 10 Hz, the maximum SNR of WPLI measures has a phase-lag of about 25 ms. For these reasons, the functional relationship between distant areas like the frontal and occipital cortices might be captured thanks to this longer phase-lag that allows pick-up of the coupling between signals with longer propagation time.
To summarize, our findings suggest that brain-state dependent stimulation could benefit from taking into account a broader concept of "state" of the system. Depending on whether brain-state dependent stimulation is practiced at rest or during a task, one should consider the functional coupling patterns between the areas of interest in the frequency band whose phase is used as a trigger. Furthermore, we suggest that the physical distance between the regions could also be taken into account when choosing the coherence metric to assess functional coupling based on a consistent phase relationship in the frequency band of interest.

Conclusions
Brain-state dependent brain-stimulation has the potential to reliably modulate specific neural pathways using a "Closed-Loop" approach [44]. How to achieve a time-resolved realtime estimation of EEG-derived brain connectivity states remains a key challenge. In this study, we addressed the off-line estimation of phase-based connectivity patterns between different regions of interest, presenting a pipeline to determine candidate pathways for real-time paradigms, for instance where the trigger-signal is extracted from one cortical site, and the stimulation targets a different site, or considering the target region and downstream effects to distal cortical regions.
In summary, offline connectivity measures have a twofold purpose: first, the estimate of connectivity patterns before and after can represent a biomarker of cortico-cortical changes after single and repetitive state-dependent TMS. In fact, for mu-alpha, stimulation at troughs was found to enhance Transcranial evoked potentials (TEPs) in both ipsi and contralateral hemisphere at 100 ms after stimulus even when the stimulus was provided at 90% of motor threshold [10]. Moreover, the same kind of stimulation has been found to be responsible for an enhancement of connectivity between the two sensorimotor cortices [45]. Connectivity joined with TEP analysis could be even more relevant when analyzing nonmotor areas such as the dorsolateral prefrontal cortex, where one cannot rely on as reliable and consistent an outcome as that of the MEPs.
Second, an off-line connectivity analysis can be used as prior information for future online estimates, to determine a correlation between areas oscillating at the same frequency or in an anticorrelation, a phase state which was previously suggested as a possible marker of functional segregation [46], or no correlation at all. The offline measures can serve as a benchmark for determining the optimal trade-off between the temporal resolution (windowlength) and the estimator error variance for a given experimental condition, noting that the window-length must be short enough to capture physiological transitions [47]. Institutional Review Board Statement: Data analyzed in this work were obtained from the CamCAN repository. Data were collected in compliance with the Helsinki Declaration (Cambridgeshire 2 Research Ethics Committee; reference: 10/H0308/50). See [24] for details.
Informed Consent Statement: Data analyzed in this manuscript are publicly available from of the Cambridge Centre for Aging and Neuroscience (CamCAN) public dataset [25]. No data were collected by the authors for this study. For information about original informed consent see the CamCAN dataset documentation [25].