Cross-frequency synchronization connects networks of fast and slow oscillations during visual working memory maintenance

Neuronal activity in sensory and fronto-parietal (FP) areas underlies the representation and attentional control, respectively, of sensory information maintained in visual working memory (VWM). Within these regions, beta/gamma phase-synchronization supports the integration of sensory functions, while synchronization in theta/alpha bands supports the regulation of attentional functions. A key challenge is to understand which mechanisms integrate neuronal processing across these distinct frequencies and thereby the sensory and attentional functions. We investigated whether such integration could be achieved by cross-frequency phase synchrony (CFS). Using concurrent magneto- and electroencephalography, we found that CFS was load-dependently enhanced between theta and alpha–gamma and between alpha and beta-gamma oscillations during VWM maintenance among visual, FP, and dorsal attention (DA) systems. CFS also connected the hubs of within-frequency-synchronized networks and its strength predicted individual VWM capacity. We propose that CFS integrates processing among synchronized neuronal networks from theta to gamma frequencies to link sensory and attentional functions. DOI: http://dx.doi.org/10.7554/eLife.13451.001


Introduction
Working memory (WM) has a limited capacity of 3À4 objects (Luck and Vogel, 1997) and is comprised of sensory storage and central executive control for manipulating the stored information and supporting the preparation of a contingent response (Baddeley, 1996;Miller and Cohen, 2001;Sreenivasan et al., 2014). During visual WM (VWM) maintenance, these functions are achieved in visual and fronto-parietal (FP) brain areas, which exhibit enhanced neuronal activity levels during memory maintenance in both monkey single-and multi-unit recordings Goldman-Rakic, 1995;Siegel et al., 2009) and functional MRI (fMRI) of humans (Munk et al., 2002;Pessoa et al., 2002;Todd and Marois, 2004). While the visual cortex is responsible for the processing of visual information and its maintenance in VWM (Riggall and Postle, 2012;Emrich et al., 2013;Kravitz et al., 2013), lateral prefrontal cortex (LPFC) is thought to regulate VWM maintenance, represent goals and rules, and govern response selection (Miller and Cohen, 2001;Rowe et al., 2000;Sreenivasan et al., 2014;Markowitz et al., 2015). Recent studies show that even in primate prefrontal cortex, these VWM functions are achieved in multiple dissociable networks (Markowitz et al., 2015;Lundqvist et al., 2016).
Phase synchronization (PS) is a mechanism for integrating anatomically distributed processing and regulating the communication between distant assemblies (Singer, 1999;Fries, 2005Fries, , 2015. Macaque cortical recordings have shown that in beta- (b,(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) and gamma-(g, 30À90 Hz) frequency bands, inter-areal PS is indeed enhanced by the maintenance of visual information in VWM and it has been suggested to support the maintenance of VWM content (Siegel et al., 2009;Salazar et al., 2012). Using MEG and EEG, source-localization, and network analyses (Bullmore and Sporns, 2009), we observed that VWM involves large-scale b-and g-band PS in the visual and parietal regions as in invasive LFP recordings from monkeys, but also high-alpha (a)-band (11-14 Hz) synchronization in fronto-parietal areas (Palva et al., 2010). This anatomical pattern suggests a putative functional division where b-and g-band synchronization underlies the maintenance and integration of visual information while the high-a-band synchronization supports the coordination of processing among brain regions serving attentional and central executive control. Likewise in MEG, also human intracranial EEG (iEEG) data have revealed long-range phase-synchronization in theta (q, 4-8 Hz), b, and g bands in hippocampal memory circuits (Fell et al., 2001;Axmacher et al., 2008;Burke et al., 2013).
Similar conclusions have also been made in studies addressing the relationship of VWM performance and oscillation amplitude or power modulations that reflect local cortical synchronization. The amplitudes of both b-and g-band oscillations are positively correlated with the task performance and individual capacity limits (Sauseng et al., 2009;Roux et al., 2012;Honkanen et al., 2015) and thought to support the neuronal representations of VWM contents (Tallon-Baudry et al., 1998;Honkanen et al., 2015). In contrast, the plateau of load-dependent strengthening in a-frequency band amplitudes during VWM is negatively correlated with the individual VWM capacity  and a oscillations are stronger for distractors than memorized items (Bonnefond and Jensen, 2012;Park et al., 2014). Together with a large body of other evidence (Palva and Palva, 2007;Buffalo et al., 2011;Saalmann et al., 2012;Haegens et al., 2015), these findings suggest that q and/or a oscillations underlie the attentional functions of VWM by facilitating the processing of relevant and/or by suppressing irrelevant information .
Hence, a crucial challenge now is to reveal what kinds of cross-frequency coupling (CFC) mechanisms integrate the neuronal processing that during VWM maintenance is distributed to multiple frequency bands in anatomically distinct networks. These mechanisms would conceivably also support the emergence of subjectively coherent VWM from sensory and central executive functions. Phase amplitude coupling (PAC) is one form of CFC that has been suggested to underlie cross-frequency integration (Palva et al., 2005;Palva and Palva, 2007;Sirota et al., 2008;Canolty and Knight, 2010;Fell and Axmacher, 2011;Giraud and Poeppel, 2012;Lisman and Jensen, 2013). Prior studies have shown that during VWM maintenance, g oscillation amplitudes are coupled to q phase via PAC in cortical (Sauseng et al., 2009;Siegel et al., 2009) and hippocampal recordings (Axmacher et al., 2010;Park et al., 2016). However, PAC is independent of the phase of the fast oscillations and hence cannot promote a consistent temporal and spike-time relationship between the slow and fast oscillations (Palva et al., 2005;Fell and Axmacher, 2011;Palva and Palva, 2012), which are essential in the regulation of neuronal communication (Singer, 1999;Fries, 2015).
In contrast to PAC, n:m-cross-frequency phase synchrony (CFS) (Tass et al., 1998;Palva et al., 2005;Palva and Palva, 2007) is a form of CFC that indicates direct phase coupling of the slow and fast oscillations. CFS could thus mechanistically underlie the cross-frequency integration of synchronized networks at the temporal resolution of the faster oscillation (Fell and Axmacher, 2011;Palva and Palva, 2012) through coordination of consistent spike timing relationships between the oscillations. In CFS, the ratio of high (f high ) and low (f low ) frequencies is defined by the integers n and m so that n. f high = m. f low . In humans, CFS has been observed in resting-state (Nikulin and Brismar, 2006) and VWM tasks among EEG sensors (Schack et al., 2005;Hamidi et al., 2009;Sauseng et al., 2009) and in a mental calculation WM task among MEG sensors (Palva et al., 2005). However, being limited to the sensor level, these studies leave the cortical sources of CFS as well as its functional significance unclear. A recent study using intracranial EEG (iEEG) also reported behaviorally significant q:b and q:g CFS during VWM in human hippocampus using a serial Stenberg task (Chaieb et al., 2015). Additionally, CFS has also been observed in the neocortical (Roopun et al., 2008) and hippocampal (Belluscio et al., 2012) microcircuits in rats.
We postulate that CFS supports the integration of neuronal processing between synchronized networks in distinct frequency bands during VWM maintenance and thereby underlies the integration of sensory and attentional functions of VWM ( Figure 1). We use here concurrent magneto-and electroencephalography (M/EEG) with individual-anatomy-based source reconstruction and datadriven analyses to investigate the functional significance of CFS in WVM. The key predictions of our postulate are that CFS should be observable and load-dependent in a parametric VWM task (Palva et al., 2010), that it should connect the task-relevant within-frequency-synchronized networks, and finally that it should be predictive of individual behavioral VWM performance. We found solid support for each of these predictions, which shows that CFS may indeed be a mechanism for coordinating the neuronal processing and regulating communication among oscillatory networks at different frequencies.

Inter-areal CFS is a robust phenomenon during VWM retention period
We assessed here the functional role of CFS in VWM by analyzing M/EEG data acquired during a delayed match-to-sample VWM task (Luck and Vogel, 1997;Palva et al., 2010, which had earlier revealed concurrent large-scale 1:1 synchronization in a-, b-, and g-frequency bands (Palva et al., 2010). In each trial of the task, the subjects were presented a to-be-memorized 'Sample' stimulus containing one-to six-colored squares and, after a memory retention period of 1 s, a 'Test' stimulus, wherein 50% of trials one square was of a different color than in the Sample (Figure 2a). The subjects gave a forced-choice response about whether the Test was same or different than the Sample. The hit rate (HR) decreased as a function of VWM load but remained above chance level for the highest load of six objects ( . VWM modulates inter-areal cross-frequency phase synchrony (CFS) in human cortex. Example of prominent single-trial CFS during the VWM retention period. (a) Illustration of the delayed-match-to-sample experimental paradigm where the Sample stimulus containing 1-6 colored squares is presented for 0.1 s and then after a 1 s retention period a Test stimulus appears and the subject responds whether the Test contains an object with a color different from than in Sample. The flattened cortical surface shows an example of 1:4 CFS (purple line) transiently phase coupling a (13 Hz) oscillations in the right medial frontal sulcus (MFS) with g (54 Hz) oscillations in left occipital pole (Occ.Pole) during VWM retention. These a and g oscillations were concurrently also 1:1 synchronized in networks of other cortical areas (blue and red lines, respectively). The cortical surface is colored by brain systems as indicated in Figure 4. ( The following source data and figure supplements are available for figure 2: We estimated the phase time series of single-trial ongoing neuronal activity for all cortical grey matter sources by filtering artefact-cleaned M/EEG recordings with Morlet wavelets ranging from 3 to 90 Hz in frequency and using cortically-constrained minimum-norm estimates (Lin et al., 2006) for the inverse transform. These source estimates were collapsed with individually optimized collapse operators (Korhonen et al., 2014) to neuroanatomically labeled cortical parcellations of 400 parcels (see Figure 2-figure supplement 2).
To first re-illustrate with real data the concept (see Figure 1) of CFS connecting cortical parcels engaged in 1:1 synchronization in a and g bands, we selected two cortical parcels in a single subject: middle frontal sulcus (mFS) in lateral prefrontal cortex (LPFC) for a oscillations and occipital pole (early visual cortex) for g oscillations (Figure 2b), and estimated their phase time series (Figure 2c). We then quantified n:m CFS (Palva et al., 2005;Tass et al., 1998) and 1:1 synchronization with the phase-locking value (PLV) in sliding time windows ( Figure 2d) and found in single trials that CFS between these regions as well as 1:1 a and g band synchronization in separate networks were strengthened during VWM retention.
To characterize this phenomenon at the group level (N = 12), we quantified event-related CFS obtained across trials for each pair of frequency bands and each pair of the 400 cortical parcels. Mapping the correlation networks at the whole-brain level is crucial for understanding neuronal communication in cognitive architectures (Petersen and Sporns, 2015). As the n:m ratios, we used n = 1 and m = 2À9 and assessed CFS in two pre-stimulus and two VWM-retention period time windows (0.4À0.7 and 0.7À1 s from Sample onset). The retention period windows were free of stimulusevoked activity , which is a major confounder in analyses of both inter-areal PS and all forms of CFC (Palva and Palva, 2012). For robust statistical mapping of CFS at the group level, we collapsed the parcel-parcel adjacency matrices of CFS strength (PLV) to a coarser parcellation, the Destrieux atlas (Destrieux et al., 2010), where the cortical surface is divided into 148 parcels (see We first addressed whether the inter-areal CFS was strengthened or suppressed during VWM retention compared to the pre-Sample baseline (two-sided t-test, a = 0.05, corrected). We initially averaged the CFS estimates across the six memory-load conditions ('Mean condition') and represented significant differences in CFS between baseline and retention as graphs where cortical areas comprise the vertices and significant connections are edges (Bullmore and Sporns, 2009). These task effects on the CFS were indexed at the graph level by connection density, K, that was the proportion of statistically significant inter-areal interactions of all possible pair-wise interactions among the 148 brain areas. Visualizing K as a function of f low and the frequency ratio 1:m revealed that inter-areal CFS of high-a (11-14 Hz) with b-and g-band oscillations as well as CFS of high-q-band (5À7 Hz) with a-, b-and g-band oscillations were enhanced above baseline levels for most or all frequency ratios (Figure 3a, see external source data 1 for statistical information, Siebenhühner et al., 2016). In contrast, CFS of low-q (3-4 Hz) and low-a (7À10 Hz) oscillations with their corresponding higher frequency bands was concurrently suppressed below baseline levels ( Figure 3b; external source data 1, Siebenhühner et al., 2016). These data hence show that CFS among cortical oscillations is both up-and down-regulated dynamically by VWM maintenance. The horizontal concentration of positive and negative findings across the frequency ratios to consistent lower frequencies (see Figure 3a and b) suggested that the high-a and q oscillations were synchronized in a harmonic structure with higher frequency oscillations in aÀg bands. We quantified this by comparing the harmonic consistency of CFS across frequency ratios against surrogate data (see Materials and methods) and found this consistency to be significant for the enhanced CFS in high-q and high-a bands and for the suppressed CFS in the low-a band (Figure 3a and b, p < 0.05, FDR corrected).

Cons.
Cons. . K+ values (color scale) are shown for all studied pairs of 1:m ratios from 1:2 to 1:9 (x-axis) and lower frequencies from 3.3 to 45 Hz (y-axis) and represent data from an adjacency matrix for each frequency-pair. The grey line marks the boundary set by the highest investigated frequency (90 Hz). CFS of high-q and high-a with their upper frequencies was increased for essentially all ratios. Right: The brown line indicates the harmonic consistency of CFS across low frequencies. The dashed line denotes the 95%-ile confidence limit. CFS of high-q and high-a oscillations with their harmonics at higher frequencies is significantly consistent across ratios. (b) Fractions of inter-areal CFS connections (K-) that were suppressed below baseline levels during the retention period. Harmonic consistency (blue line) was estimated as in and shows that low-a consistent CFS was suppressed at all ratios.(c) Fractions of CFS connections that were significantly positively (K+) or (d) and negatively (K-) correlated with VWM load (Spearman Rank correlation test of CFS across the six VWM memory load conditions, p < 0.05, corrected). DOI: 10.7554/eLife.13451.007 The following figure supplements are available for figure 3:  These data thus suggest that q and high-a oscillations are systematically phase coupled with b and g networks for which they thus provide temporal framing across a range of frequency ratios.
To then assess the relationship of CFS with progressive VWM task demands, we asked whether the strength of CFS was systematically correlated with the memory loads from 1 to 6 objects ('Load condition'). We found that CFS between the high-a-band consistent frequency pairs was positively correlated with memory load (Spearman's rank correlation test, a = 0.05, false discovery rate (FDR) corrected) while for q-consistent CFS, VWM load-dependent modulations were weak. Further, low-a consistent CFS was negatively correlated with VWM load ( To consolidate these observations and ensure that they were not driven by outliers or attributable to effects that would be unreliably sampled at the current cohort size, we performed two lines of control statistics for positive observations in the mean and load conditions. First, using a leave-oneout approach, we created 12 datasets of 11 subjects and performed the group statistics as earlier.
We then examined each significant parcel-parcel CFS observation of the prior group statistics (see Figure 3a-d) and considered significant only those connections that were significant in every leaveone-out cohort. Second, we estimated the effect size that could be estimated at statistical power 1Àb of 0.8 with the original cohort size and then from the results of the original group statistics accepted as significant only observations exceeding this limit. Neither the leave-one-out (

CFS is not attributable to concurrent event-related amplitude changes
The results so far showed that CFS is a robust, dynamic, and task-dependent phenomenon in human cortical activity and is, both theoretically and in the light of the present empirical data, a plausible mechanism for the integration of neuronal processing distributed across frequencies. However, event-related changes in CFS could be attributable to changes in the signal-to-noise ratio (SNR) caused by task-induced changes in oscillation amplitudes. To rule out this possibility, we obtained source-space parcel-specific estimates of apparent SNR (aSNR) (Figure 3-figure supplement 3) and used the event-related and memory-load dependent changes in aSNR to predict the changes in CFS estimates that one would observe in the absence to true changes in phase coupling (Palva et al., 2010) see Materials and methods). However, we did not find any systematic relationship between SNR and CFS, and observed that as in similar prior analyses for 1:1 PS (Palva et al., 2010), the observed CFS changes were much larger than those predicted by changes in SNR in both the Mean ( To further corroborate that the observed CFS effects were not attributable to SNR changes in a case where the SNR would be overestimated, we quantified the correlation between the Mean-and Load-condition effects in PLV and oscillation amplitudes for each significant edge, for each ratio, and for each frequency-pair (see Materials and methods). In line with the prior findings, we found that the task effects on oscillation amplitudes were either insignificantly correlated with the CFS effects or, in the few ratio-frequency pairs where a significant correlation was observed, explained   . q-and high-a-consistent CFS is observed in fronto-parietal and visual networks. (a) Low-frequency CFS hubs (LF, left) and high-frequency CFS hubs (HF, right) and their most central connections during the VWM retention period for q-consistent CFS averaged across memory loads (Mean condition) and shared across frequency ratios from 1:2À1:9 (see Materials and methods for details) displayed on a flattened cortical surface. These networks and the most central low-and high-frequency hubs are presented here separately to illustrate the distinct neuroanatomical sources of the most central cortical q oscillations (left) and the aÀg oscillations (right) connected by this CFS coupling. The complete q-consistent CFS network is thus the combination of these graphs. Circles represent network vertices (cortical parcels), their size is proportional to their degree centrality, and both vertex and surface colors represent functional brain systems (see colors below). The lighter sector of each vertex represents the fraction of lowfrequency connections and the darker section the fraction of high-frequency connections this vertex has in the graph. The width of the edges represents the coupling strength between the two vertices and their color is a mixture of the brain system assignments of the vertices. The dotted lines at the end of the connections point towards the high-frequency (from LF hubs, left panel) or the low-frequency (from HF hubs, right panel) vertices for LF and HF networks, respectively. (b) As in (a), but for the high-a-consistent CFS shared across ratios from 1:3À1:9.  To further test whether event-related changes in non-sinusoidal/spike-like waveforms (Palva et al., 2005) could have led to artificial observations of CFS, we estimated cross-frequency amplitude-amplitude correlations. For CFS attributable to non-sinudoidal waveforms, these amplitude correlations should exhibit a frequency-ratio pattern similar to that of CFS. However, the observed pattern of CF amplitude correlations was saliently distinct from CFS (Figure 3-figure supplement 7; external source data 2, Siebenhühner et al., 2016) and within each frequency-ratio pair, the vast majority of significant CFS connections did not exhibit significant amplitude-amplitude coupling (Figure 3-figure supplement 8). These analyses thus rule out the possibility that changes in SNR or in waveform non-sinusoidalities could have contributed significantly to the observations of CFS in this study.

CFS during VWM retention links visual and attentional systems
If CFS supports the integration of spectrally distributed processing of attentional and representational functions during VWM, it should be localized to the task-relevant visual and attentional brain systems, and the strength of CFS in these regions should be predictive of behavioral performance. To illustrate the most robust CFS networks in these data, we identified the connectivity backbones of high-a-and q-consistent CFS networks across the frequency ratios 1:2À1:9 by pooling the adjacency matrices of significant connections, rejecting connections that were significantly attributable to any single ratio, and visualizing the most central connections and network hubs (see Materials and methods).
As CFS networks are asymmetric, the visualization was performed separately for the low-and high-frequency parts of CFS networks. In the visualization, we also accounted for the M/EEG signal mixing effects with a novel clustering approach that bundles the connections more likely to reflect a single true underlying neuronal interaction (see Materials and methods). The brain regions putatively supporting attentional and other functions were identified by co-localizing our parcellation with the We first visualized the most central hubs and connections of q-harmonic networks for ratios from 1:2 to 1:9 in the Mean condition separately for the low-frequency (LF, here q) and high-frequency (HF, here aÀg) parts of the CFS interactions. LF network hubs, i.e., areas wherein the q oscillations were coupled with aÀg oscillations elsewhere, were found mainly in the fronto-parietal (FP), dorsal (DA), and ventral (VA) attention systems and included the right-hemispheric middle frontal gyrus and sulcus (mFG, mFS), intraparietal sulcus (intPS), and bilateral inferior parietal gyri (iPG) (see large circles in Figure 4a, left; external source data 3, Siebenhühner et al., 2016). These q hubs were CFS-coupled with HF oscillations in visual, parietal, and frontal cortices (see dashed ends of connections). HF hubs, i.e., loci where aÀg oscillations were predominantly coupled with q elsewhere, were observed in the DA system, namely in the left inferior precentral sulcus (iprCS), right superior parietal gyrus (sPG), and intPS as well as in postcentral sulcus and gyrus (poCS/G) of the sensorimotor system (SM) and superior frontal gyrus (sFG) of the default-mode (DM) system (Figure 4a, right).
We next visualized the most central hubs and connections of the high-a-CFS network with the same approach. For the CFS between high-a and g, in networks with ratios 1:3À1:9, we observed LF (high-a) hubs in the right-hemispheric iPG, mFS, and iFS in the FP system (see blue hubs in Figure 4b, left panel) and in bilateral sFG, right superior frontal sulcus (sFS), and inferior temporal gyrus (iTG) in the DM system. Among other targets, the FP a-hubs were bilaterally and strongly connected to g oscillations in the visual cortex (see dashed line ends of red FP-Vis connections in Figure 4b, left panel). The HF (here, g) hubs included mFG and mFS (FP), bilateral iprCS and left Figure 6. CFS connects task-relevant functional brain systems. (a) Connection density K for q consistent CFS (ratios 1:2-1:9) in the Mean condition within and between functional brain systems of the Yeo parcellation. The x-axis shows the subsystems for the lower frequency and the y-axis for the higher frequency oscillations. K is shown for those brain system pairs where K exceeded the 95%-ile of the corresponding K values obtained with surrogate data. (b,c) Same as (a), but for high-a consistent CFS (ratios 1:3-1:9) in the Mean and Load conditions, respectively. (d,e) Same as (a), but for a:b CFS (ratio 1:2) in the Mean and Load conditions, respectively. (f) K of low-a consistent CFS decreases (negative Mean condition) within and between systems. DOI: 10.7554/eLife.13451.021 The following source data is available for figure 6:  intPS in the DA system and, importantly, the right-hemispheric cuneus (CN) and superior occipital gyrus (SOG) in the visual system. These visual system g hubs were connected to high-a-oscillations in LPFC and with left-hemispheric posterior parietal cortex (PPC) and visual system (Figure 4b, right). CFS between high-a and g hence connected the FP and DA internally and these attentional systems with the visual system. Moreover, in both the q-and high-a-consistent CFS networks (see Figure 4a and b), the LF hubs were more salient in right-hemispheric FP areas while HF hubs were more prominent in bilateral DA and visual areas. This tentative dichotomy is well in line with the idea that q/a oscillations serve higher attentional/executive functions (FP) while g oscillations serve visual representational and lower-level attentional (DA) functions. Both LF and HF hubs were also observed in DM and most commonly in sFG and sFS therein.
In contrast with the 1:3À1:9 a:g CFS, 1:2 CFS connections between high-a and b oscillations were saliently localized to the SM, DA, and VA systems (Figure 4-figure supplement 1; external source data 5, Siebenhühner et al., 2016) and were thus visualized separately from a:g CFS. The differences between a:b and a:g CFS suggest that they may underlie functionally distinct CFC.
To localize the suppression of CFS and address whether it was a characteristic of the task-relevant or -irrelevant circuitry, we visualized the hubs and connections for the low-a consistent CFS networks that exhibited the greatest suppression during VWM maintenance. We found that both LF and HF hubs were predominantly found in the DM and SM systems that are putatively task-irrelevant in VWM (Figure 4-figure supplement 2; external source data 6, Siebenhühner et al., 2016). However, a LF low-a hub was also observed in intPS (DA) and HF b / g hubs in iprCS (DA) and mFG (FP), which were connected to the nodes in the DM system.
We then visualized the most central connections and nodes of the high-a consistent networks that were positively correlated with VWM load ('Load condition'). High-a (LF) hubs were mostly found in the left-hemispheric visual system, right-hemispheric FP, including the mFS, and iPG, and in SM bilaterally ( Figure 5; external source data 4, Siebenhühner et al., 2016). The HF g hubs were observed in the areas overlapping with those of high-a hubs in the left-hemispheric visual system and in right-hemispheric mFG and sPG. Interestingly, mFS and iFS were key frontal HF nodes and were connected to the visual system similarly to what was observed in the Mean condition a(LF) networks. This suggests that the reciprocality of a and g oscillations in these areas is enhanced with increasing VWM load. Further, similarly to the Mean condition, also the load-dependent 1:2 CFS connections between high-a and b oscillations were localized to the SM, DA and VA networks (  with Benjamini-Hochberg FDR method). (b) Simplified illustration of CFS connecting 1:1 networks. The graph shows 15 brain areas with the greatest degree centrality in significantly 1:1 phase synchronized (Mean-condition) qband (green) and a-band networks (pink). Blue borders indicate the brain areas that belong to the most central low-frequency CFS-hubs while the violet borders indicate the brain areas that are among the most central highfrequency CFS hubs. Significant CFS coupling between these networks is shown in a separated manner for clarity: left column shows CFS coupling of the low-frequency CFS hubs with all of their targets in the 1:1 higher-frequency network. Conversely, the right column shows the CFS coupling of the high-frequency CFS hubs with all of their targets in the 1:1 lower-frequency network. Individual brain areas may appear both as low-and high-frequency hubs within the same graph. (c) CFS coupling between 1:1 high-a-band network (green) and the low-g-band network (pink). All visualization details as in (b). DOI: 10.7554/eLife.13451.023 The following source data and figure supplements are available for figure 7:

Brain-system-level connectivity analysis of the CFS effects
To further elucidate CFS connectivity among brain systems, we complemented the detailed network visualization of the strongest CFS connections (see Figures 4 and 5) by assessing the densities of all significant CFS connections, K, within and between the seven functional brain systems. These K values were compared against non-parametric statistical thresholds obtained by shuffling the adjacency matrices. The system-system connection density matrices largely corroborated the main phenomena found in the network visualizations and showed that VWM retention was associated with significantly denser CFS coupling of q-oscillations in DA, FP, visual, SM, and DM systems with aÀg oscillations in SM and DA than expected by chance (Figure 6a, see Figure 6-source data 1). q-oscillations in DA were also synchronized with aÀg oscillations in visual, FP, and limbic systems. In line with the visualization of the most central connections (See Figure 4b), high-a oscillations were synchronized with g oscillations within DA, and between DA, visual, SM, and FP systems ( Figure 6b). As a notable feature in both of these analyses (see Figure 6a and b, cf. Figure 4a and b), was that HF oscillations in the visual system were connected with LF oscillations in DA. However, in the Load condition, high-a oscillations in SM and visual systems were synchronized with g oscillations in visual, SM, and DA systems, and conversely g oscillations in FP with high-a oscillations in DA, limbic, and DM systems (Figure 6c).
Further, in line with the visualization of strongest connections, 1:2 a:b CFS was strengthened in Mean and Load condition among SM and attentional systems (Figure 6d and e, cf. Figure 4-figure supplement 1) as well as also between DA and FP. Low-a oscillations in SM, DA, FP and DM showed suppressed synchronization with b and g oscillations in DM as well as in DA, SM and FP systems ( Figure 6f).
Thus, these brain-system-level connectivity analyses showed that CFS connected high-a oscillations in the attentional FP and DA systems with the g oscillations in the attentional as well visual and SM systems in the Mean condition. In the Load condition, CFS was additionally strengthened within the visual system. For the q-consistent CFS, the connectivity pattern was more widespread and q oscillations in many brain systems were connected to attentional networks.
CFS connects the hubs of within-frequency (1:1) synchronized networks If CFS mediates the integration of slow and fast 1:1 phase-synchronized networks, it should connect the most central hubs of these corresponding within-frequency (1:1) synchronized networks. To assess this hypothesis, we tested whether the most central areas of the low-(f low ) and high-frequency The following source data is available for figure 9: Source data 1. Statistical (f high ) CFS networks, i.e., the CFS network hubs, were co-localized with the hubs of the corresponding 1:1 PS networks at f low and f high. (See Materials and methods). We found that for both the high-q (5À7 Hz) and high-a (11À14 Hz) consistent CFS networks, the node degree centrality values were correlated with the corresponding centrality values of within-frequency synchronized networks, (p < 0.05, corrected, Figure 7a, see Figure 7-source data 1). Hence, the low-frequency CFS hubs were significantly co-localized with the hubs of he slow 1:1 networks. The corresponding correlations of high-frequency CFS hubs with those of the fast 1:1 networks were weaker but still significant for several frequency pairs (Figure 7a). These findings thus support the notion that CFS connects the 1:1 networks.
We then used simplified graphs to illustrate the strongest hubs of the fast and slow 1:1 networks, their within-frequency 1:1 connectivity, and how these two networks were mutually connected by CFS. We first visualized the main hubs and their connections of 1:1 q and a networks and the 1:2À1:3 q:a CFS linking these networks by illustrating separately how the 1:1 q network was connected with targets in the 1:1 a network (Figure 7b, left column) and how the 1:1 high-a network was connected to targets in 1:1 q network (Figure 7b, right column). This showed that for LF-CFS connections, 1:1 q and LF CFS hubs were co-localized in mFG and iFG of the LPFC and in CS (Figure 7b, left column) while iPG and mFG were hubs for both the 1:1 a network and HF-CFS ( Figure 7b, right column). For 1:4À1:5 q:b CFS and for 1:6À1:9 q:g CFS, shared hubs for 1:1 q network and LF-CFS network were observed in the poCS and LPFC while shared HF-CFS and fast 1:1 hubs were observed in posterior parietal cortex (PPC) and LPFC and for or q:g CFS also in the visual system (Figure 7-figure supplement 1).
For high-a consistent CFS, we first examined the coupling of the high-a band 1:1 network with the low-g band 1:1 network (Figure 7c), as well as with b-band and high-g band networks (Figure 7figure supplement 2). This revealed a very similar pattern of CFS connections than that presented in Figure 4. Shared key hubs in both 1:1 a and LF-CFS networks were mFS, sFG, both angular and supramarginal parts of iPG, and precentral gyrus (prCG) which were connected to their g network targets in the visual system through CFS (Figure 7c, left column). Hubs shared by low-g band and HF-CFS were SOG and inferior temporal gyrus (iTG), which were connected with their a-network targets is LPFC and PPC (Figure 7c, right column). Thus as proposed by our hypothesis, high-a band network included co-localized CFS and 1:1 hubs in the FP and DA attention networks while the gband network included key hubs in the visual system and posterior brain areas. A similar pattern of co-localized 1:1 and CFS hubs was observed also for 1:2 high-a:b (Figure 7a and Despite both low-and high-g networks being characterized by hubs in the visual system, they were not mutually significantly correlated in parcel centralities (Figure 7-figure supplement 3), which implies that these g oscillations were relatively narrow-band and can not particularly be explained by broad-band g activities. Overall, these data show that, as hypothesized, CFS connects 1:1 slow and fast networks through common hubs that are most prominent in PPC and LPFC as suggested also by data presented in Figure 4.
Phase-amplitude coupling (PAC) and CFS are distinct phenomena PAC has been suggested to underlie cross-frequency coupling of slow and fast oscillations and prior studies have found PAC, e.g., between the q phases and g amplitudes during VWM (Axmacher et al., 2010;Siegel et al., 2009). To address the possibility that CFS and PAC would underlie similar cross-frequency integration, we tested whether also PAC was correlated with the VWM demands in the present data. We computed PAC between all parcels and for the same frequency pairs used in the CFS analysis. Somewhat surprisingly, as the strongest effect, this analysis revealed that q-and low-a-band phases were coupled with a-and b-band amplitudes at a 1:2 ratio. We also found weak coupling of the g oscillation amplitudes with the phase of q oscillations at frequency ratios up to 1:9 (Figure 8a; external source data 7, Siebenhühner et al., 2016). Importantly, in the high-a band, where prominent CFS was observed, PAC was suppressed (Figure 8b). Loaddependent modulations of PAC showed a similar spectral profile (Figure 8c-d). These data show that PAC and CFS are clearly distinct phenomena.

The strength of CFS predicts VWM accuracy and individual capacity
Finally, if CFS was functionally significant in the integration of neuronal processing across q-to gband synchronized networks, the strength of CFS should be predictive of individual VWM task performance. We thus quantified the predictive correlation between inter-individual variability in the strength of VWM-retention period CFS and the inter-individual variability in VWM capacity across subjects. CFS network strength, S, was measured for each subject, VWM load, and low frequency of low (1:2-1:5) and high (1:6-1:9) CFS frequency ratios, where S was the sum of phase-locking values of all significant connections in the mean condition (positive). VWM capacity C was obtained separately for each subject and VWM load (See Materials and methods). We estimated the correlation (Pearson correlation test) between S and C and found that for low ratios the strength of high-a consistent CFS networks predicted (r = 0.343, p = 0.003) VWM capacity (Figure 9a, see Figure 9source data 1). For the high (1:6À1:9) ratios, the individual capacity was correlated (r = 0.376, p = 0.001) with high-q consistent CFS (see Figure 9a). These values remained significant when testing all low frequencies for such correlations and correcting for multiple comparisons (Benjamini-Hochberg method). We also tested whether the high-a low ratios observation was attributable only to 1:2 a:b coupling but found both the high-a:b (ratio 1:2, r = 0.247, p = 0.036) and high-a:low-g (ratios 1:3À1:5, r = 0.232, p = 0.0496) to be significant. We also estimated (see Figure 9b, Figure 9-source data 2) the CFS-capacity correlation of all individual ratios to illustrate their systematic contribution to the correlations of their mean (see Figure 9a) in the high-a and q frequency bands.
We then estimated correlations of PAC network strength with individual VWM capacity in the same manner (Figure 9c,d). However, no p values exceeded the corrected significance threshold. These data thus indicate that cortical CFS is functionally significant in the human VWM maintenance and also further reinforce the view that CFS and PAC are functionally distinct phenomena.

Discussion
WM is comprised of sensory storage and attentional control for manipulating the stored sensory information (Baddeley, 1996;Miller and Cohen, 2001;Sreenivasan et al., 2014). While the maintenance and integration of sensory information in WM is tentatively supported by b and g phase-synchronization, q and a phase-synchronization appears to play a role in the coordination of attentional functions. We hypothesized that concurrently synchronized assemblies in different neuronal circuits and frequency bands could be integrated into coherent VWM by cross-frequency phase synchrony (CFS).
We used M/EEG and a VWM task to investigate whether CFS would connect networks synchronized in q-, a-, b-, and g-frequency bands and integrate the associated neuronal processing in visual and attentional systems during VWM maintenance. We demonstrate here that human VWM maintenance is indeed characterized by large-scale networks of functionally significant CFS among cortical oscillations. We found that CFS was correlated with task demands, predicted individual VWM capacity and connected 1:1-within-frequency synchronization among task-relevant visual and fronto-parietal regions. Importantly, CFS could not be explained by changes in the signal-to-noise ratio (SNR), concurrent amplitude-amplitude correlations, or non-sinusoidal waveforms, which supports the interpretation that CFS indeed reflects dynamic cross-frequency coupling of narrow-band neuronal oscillations. The present findings thus show that CFS reflects true dynamic cross-frequency coupling of neuronal oscillations among task-relevant cortical circuits and may thus constitute a so far little acknowledged computational mechanism for integrating and coordinating neuronal processing across distinct frequency bands. Localization of the cortical sources and finding evidence for the functional significance of CFS in the present study thus extend the scarce prior observations of CFS in neuronal circuits, which have so far revealed CFS in sensor-level MEG or EEG analyses (Hamidi et al., 2009;Nikulin and Brismar, 2006;Palva et al., 2005;Tass et al., 1998), in human (Chaieb et al., 2015) and rat (Belluscio et al., 2012) hippocampal circuits, and in the neocortical (Roopun et al., 2008) microcircuits in rats.

CFS connects visual and attentional systems and within-frequency synchronized networks therein
We found that in narrow q-and high-a-frequency bands and their harmonics, CFS was positively correlated with VWM task demands. CFS hence connected the q oscillations with those in a, b, and lowg bands as well as the high-a with b and g oscillations. Several studies suggest that the contents of VWM are maintained in the visual cortex (Riggall and Postle, 2012;Emrich et al., 2013;Kravitz et al., 2013;Sreenivasan et al., 2014;Honkanen et al., 2015), while LPFC is thought to underlie the attentional and executive regulation of the VWM maintenance (Rowe et al., 2000;Miller and Cohen, 2001;Pessoa et al., 2002;Curtis and D'Esposito, 2003;Voytek and Knight, 2010;Sreenivasan et al., 2014;Markowitz et al., 2015) in multiple functionally dissociable networks (Markowitz et al., 2015;Lundqvist et al., 2016). A complementary view is obtained from large-scale connectivity analyses of fMRI (Corbetta et al., 2000;Corbetta and Shulman, 2002;Szczepanski et al., 2013), MEG (Spadone et al., 2015), and intracranial EEG (iEEG) (Daitch et al., 2013) recordings which show that in attention tasks, functional connectivity between LPFC and PPC, i.e., within the FP and DA systems, plays a fundamental functional role in attentional control and regulation. Indeed, studies both in human fMRI (Munk et al., 2002;Mohr et al., 2006) and monkey LFP (Salazar et al., 2012) recordings have observed strengthened activity in both frontal and parietal areas also during VWM tasks.
Hence, if CFS is a functionally significant mechanism for the integration of representational and central executive functions of VWM, it should specifically connect synchronized networks across visual and attentional systems. Moreover, as the integration in brain networks has been shown to be dependent on key brain areas that are the most central nodes, i.e., the hubs in these networks (de Pasquale et al., 2016), CFS should connect the key hubs of the within-frequency synchronized networks among task-relevant, here attentional and visual, systems.
We found two lines of evidence for CFS playing such a role. First, the key hubs of CFS networks were largely localized to these brain systems and CFS connectivity linked these systems extensively. Both in the q and high-a band consistent CFS networks, the low frequency oscillations had hubs predominantly in FP and to a lesser extent in DA, and significant connectivity of these hubs to faster oscillations in the visual system and DA. Conversely, the higher frequency oscillations had hubs predominantly in DA, and in the g band also in visual cortex. Second, we found that the most central low-and high-frequency hubs of the CFS networks during VWM maintenance were anatomically colocalized with the hubs of the slow and fast, respectively, 1:1 networks and that CFS was salient among the most central nodes of the within-frequency synchronized networks between visual system and attentional networks.
The present data thus show that CFS would be well positioned to underlie the integration and regulation spectrally distributed neuronal processing between the fronto-parietal attention networks and the visual system. Thereby CFS could underlie the integration and coordination of attentional and representational functions suggested to be carried out by a and g band synchronization, respectively (Tallon-Baudry et al., 1998;Sauseng et al., 2009;Palva et al., 2010Bonnefond and Jensen, 2012;Roux et al., 2012;Park et al., 2014;Honkanen et al., 2015).

The roles of CFS in sensorimotor and default mode systems
CFS among high-a and b oscillations was distinct in network topography from higher-ratio CFS of high-a and g oscillations. a:b CFS had large hubs in SM system, strong within-SM connectivity, while still having connections with both LPFC and visual system. fMRI studies have shown that during VWM, neuronal activity is also strengthened in the SM system wherein the dorsal premotor cortex (PMC) has been suggested to underlie the maintenance and manipulation of spatial features (Munk et al., 2002;Mohr et al., 2006). Thus 1:2 high-a and b CFS in SM could reflect these functions during VWM. However, a:b CFS could also be biased by the classical sensory-motor mÀrhythm (Hari and Salmelin, 1997) that is comprised of strongly 1:2 phase coupled a and b oscillations (Nikulin and Brismar, 2006). While a:g CFS was largely distinct from a:b, also a:g CFS involved SM hubs. Together with the prior fMRI data, this supports the notion that the SM system could be taskrelevant in VWM, but further studies are needed to disentangle their functional roles.
Interestingly, both high-a-and q-consistent CFS was found also in the DM system, including superior frontal areas of LPFC as well as inferior and superior temporal areas and the angular part of iPG.
Despite DM being often categorized as a task-negative or -irrelevant system, recent studies have shown that neuronal activity in the angular part of iPG is stronger for successfully memorized information (Hutchinson et al., 2009;Lee and Kuhl, 2016;Spaniol et al., 2009) while the superior frontal cortex is partly functionally similar to mFG (Fuster, 2008), and thus both these regions could be task-positive in the present VWM task despite belonging to DM during resting state. sTS, on the other hand, is related to multisensory processing (Beauchamp et al., 2004) and verbal WM functions (Tanabe et al., 2005) and its presence in the CSF networks might be attributable to verbal rehearsal during the memory maintenance even though the present cohort was instructed not to verbalize the task. Overall, given that CFS was predominantly observed among the visual and attentional systems, we suggest that these findings reflect task-positive roles of these brain regions rather than CFS being a phenomenon characteristic to task-negative processing during VWM.

Suppression of CFS in task-irrelevant networks
Intriguingly, the positive CFS effects were paralleled by a suppression of CFS between d and low-a bands with their harmonic frequency bands. The suppression involved network hubs predominantly in the SM and DM systems and at the systems-level, disengaged DA and SM from DM and FP. The DM system is widely thought to support non-task-oriented processing (Fox et al., 2009), which suggests that the suppression of CFS there could represent the disengagement of task-irrelevant brain regions and networks, even though some parts of DM could be task-positive in VWM (see above). Low-a consistent suppression of CFS with the concurrent increase in the high-a consistent CFS among visual and FP systems may also reflect the disengagement of neuronal processing in a taskirrelevant frequency band from that in a task-relevant band. Overall, the task-negative and -positive patterns of CFS hierarchies over the low-and high-a bands, respectively, suggest that a oscillations support both the coordination of behaviorally relevant (Palva and Palva, 2007;Saalmann et al., 2012) and the disengagement of task-irrelevant neuronal processes (Klimesch et al., 2007;Jensen and Mazaheri, 2010) as suggested previously. These results suggest that CFS can underlie both the integration as well as segregation of neuronal processing suggested to be implemented via hierarchical network interactions (Sporns, 2013;Deco et al., 2015;Spadone et al., 2015).

CFS but not PAC of and high-a consistent oscillations predicts VWM performance
We hypothesized that CFS would regulate neuronal communication between fast and slow oscillatory networks. To obtain evidence for this, we tested whether CFS is correlated with task demands and whether it is predictive of the individual the behavioral performance. We found that CFS of q and high-a oscillations with their harmonics in b and g bands was strengthened above baseline levels during VWM maintenance and was positively correlated with the VWM load. As key evidence for the functional significance of this CFS, we observed that inter-individual variability of CFS network strength predicted inter-individual variability in VWM capacity. Intriguingly, the strength of high-a consistent CFS predicted individual VWM capacity at small frequency ratios while the strength of high-q consistent CFS was correlated with individual capacity at large frequency ratios. In contrast, we found no correlations between PAC and individual capacity.
As, synchronization in q and high-a bands is thought to support the attentional and executive functions of VWM (Palva et al., 2010Zanto et al., 2011;Bonnefond and Jensen, 2012;Park et al., 2014) and synchronization in g band the contents of VWM (Tallon-Baudry et al., 1998;Siegel et al., 2009;Salazar et al., 2012;Honkanen et al., 2015), the correlation of CFS between these bands with individual capacity is well in line with the suggestion that VWM capacity is set by the cross-frequency interactions among sensory and attentional rhythms (Lisman and Idiart, 1995;Lisman and Jensen, 2013;Roux and Uhlhaas, 2014). However, in contrast with prior proposals based on PAC, we found only CFS to be predictive of VWM capacity.
PAC and CFS have not earlier been directly compared in cognitive tasks-the present data thus provide evidence for that these two forms of CFC are both phenomenologically and functionally distinct. The critical difference between CFS and PAC is that CFS couples the slow and fast oscillations at a temporal accuracy of the faster oscillation. CFS, but not PAC, hence can achieve consistent spike timing between these oscillations, which could reciprocally regulate neuronal communication between the frequency bands (Palva et al., 2005). The present data thus show that the classical 1:1 'within-frequency' phase synchronization, which has been suggested to provide a mechanistic basis for regulating neuronal communication (Singer, 1999;Fries, 2005;Bastos et al., 2015;Fries, 2015), is complemented by cross-frequency phase-phase synchronization integrating cortical assemblies into hierarchical multi-scale networks.

Conclusions
We found that CFS connects slow and fast cortical oscillations across the visual, FP, and DA systems during VWM retention. The strength of CFS was dynamically up-and down-regulated, correlated with VWM task demands, and predictive of individual VWM capacity, which show that CFS is a functionally significant cross-frequency interaction mechanism during VWM maintenance. CFS hierarchies can thus be dynamically recruited to meet cognitive processing demands. The present data are in line with the notions that CFS integrates neuronal processing and regulates neuronal communication among spectrally and anatomically distributed cortical networks and that CFS could thereby support the integration of attentional and representational functions of VWM.

Materials and methods
The workflow for all steps leading to the results in Figures 3-9 is also shown schematically in

Task and recordings
The VWM task and MEG-EEG recordings of the data analyzed here are described in detail by Palva et al. (Palva et al., 2010. In the task, 12 right-handed subjects (age 28 ± 3, mean ± SD, 4 females) memorized a Sample stimulus comprising one to six randomly located squares of different colors and after a 1 s retention period answered whether the test display was the same or different as the sample stimuli and whether they were sure about the responses. Subjects were instructed not to use verbal rehearsal during memorization. A total of~1800 trials (~300 per memory load of 1À6 objects) were collected from each subject. Cortical activity was measured with concurrent EEG (60channel) and MEG data (204 planar gradiometers and 102 magnetometers, Elekta Neuromag Ltd., Helsinki, Finland) recordings at 600 Hz sampling rate. Ocular signals were measured with electrooculogram (EOG) and the behavioral thumb-twitch responses with electromyography (EMG). Only artifact-free trials with valid behavioral responses were included in the analysis. One subject was excluded from the analysis due to too large differences in head position between recording sessions. T1-weighted anatomical MRI scans for cortical surface reconstruction models were obtained at a resolution of 1x1x1 mm with a 1.5 T MRI scanner (Siemens, Germany). This study was approved by the ethical committee of Helsinki University Central Hospital and was performed according to the Declaration of Helsinki. Written informed consent was obtained from each subject prior to the experiment.

Behavioral performance
Behavioral VWM performance was measured with Hit Rate (HR) that was given by the behavioral accuracy in each VWM load condition, i.e., by the proportion of trials with correct responses of all trials that were accepted into M/EEG analyses after blink rejection.

Data pre-processing
We used the Maxfilter software (Elekta Neuromag Ltd., Helsinki, Finland) to suppress extra-cranial noise and to co-localize the recordings in signal space (Figure 3-figure supplement 1a). Independent component analysis (ICA, Matlab toolbox Fieldtrip, http://fieldtrip.fcdonders.nl (Oostenveld et al., 2011) was used to extract signal components both for the EEG and MEG data and to exclude components that were correlated with ocular artefacts identified by EOG or with heartbeat for which the reference signal was estimated from MEG magnetometers.

Filtering
Preprocessed M/EEG data were filtered (Figure 3-figure supplement 1b) into 25 approximately logarithmically spaced frequencies from f min = 3 Hz to f max = 90 Hz using Morlet Wavelets so that the filtered time series X ðt;f Þ were obtained by convolution of the original time series x ðtÞ with Morlet wavelets wðt; f Þ : X t; f ð Þ ¼ x t ð Þ w t; f ð Þ for each wavelet frequency f min f f max , where: , i is the imaginary unit, and the Morlet parameter m was m ¼ f sf ¼ 5 ( Tallon -Baudry et al., 1996).

Preparation of forward and inverse operators
Source reconstruction was performed using FreeSurfer software (http://surfer.nmr.mgh.harvard.edu/ ) for the volumetric segmentation of the MRI data, surface reconstruction, flattening, cortical parcellation, and neuroanatomical labeling with the Destrieux atlas (Dale et al., 1999;Destrieux et al., 2010;Fischl et al., 2002, see Figure 2-source data 1 for a list of parcels). MNE software (http://www.nmr.mgh.harvard.edu/martinos/userInfo/data/sofMNE.php) was used to create three-layer boundary element head conductivity models and cortically constrained source models as well as for the M/EEG-MRI co-localization and for the preparation of the forward and inverse operators (Figure 3-figure supplement 1c, Hämäläinen and Sarvas, 1989;Hämäläinen and Ilmoniemi, 1994). MEG and EEG data (204 MEG planar gradiometers, 102 MEG magnetometers, and 60 EEG electrodes) were hence integrated at the inverse transform stage of data analysis (Hämäläinen and Sarvas, 1989;Hämäläinen and Ilmoniemi, 1994). The source models had dipole orientations fixed to the pial surface normals and a 7 mm inter-dipole separation throughout the cortex, which yielded models containing 7000-8000 source vertices.

Inverse transform and vertex collapsing
A single inverse operator was used for all wavelet frequencies and was obtained by using broadband-(1-40 Hz) filtered pre-S1 baseline periods of all trials for estimating the noise covariance matrix. To reconstruct ongoing cortical phase-time series, the filtered complex single-trial M/EEG time series were inverse transformed (Figure 3-figure supplement 1d) to source-vertex time series that were then collapsed (Figure 3-figure supplement 1e) into time series of cortical parcels in a 400-parcel collection (see Figure 2-figure supplement 2). This parcellation was obtained by iteratively splitting the largest parcels of the Destrieux atlas along their most elongated axis using the same parcel-wise splits for all subjects (Palva et al., 2010;Rouhinen et al., 2013;Honkanen et al., 2015). Using neuroanatomical labeling as the anatomical 'coordinate system' made inter-subject morphing in group-level analyses unnecessary, which would have compromised the individual anatomical accuracy. Collapsing the source time-series into time-series of cortical parcels was performed with optimized collapse operators where only the source vertices with greatest source reconstruction accuracy were used (Korhonen et al., 2014). For creating the collapse operators, we used forward and inverse modeling of simulated data to assess the source reconstruction accuracy for each source vertex and cortical parcel separately for each subject. In these simulations, Morlet-filtered white noise was created independently for each cortical parcel and these time series together were forward modeled to simulate the M/EEG recording of ongoing brain activity with known source dynamics. Inverse modeling and parcel time series collapsing, carried out identically to the inverse modeling of real data, was then applied to the simulated M/EEG data to reconstruct the spatial correlation patterns attributable to signal mixing in M/EEG data acquisition and source reconstruction. Comparison of the original parcel time series with the forward-inverse modeled source-vertex time series was then used to identify for each parcel the collection source vertices that most accurately reconstruct the original parcel time series (Korhonen et al., 2014).

Removal of low-fidelity vertices
In phase synchrony analyses of collapsed parcel time series data, one major confounding factor is spurious edges resulting from signal mixing between neighbouring brain regions in data acquisition and source reconstruction (Palva and Palva, 2012). We assessed the reliability of data based on phase correlations between real and simulated data. Parcel fidelity (f) was defined as the phase correlation between the original true parcel time series and the forward-inverse modeled parcel time series (Korhonen et al., 2014). Parcel cross-talk (c), conversely, was taken as the phase correlation between the forward-inverse modeled parcel time series with the original true time series of all other parcels. For each parcel, we calculated c as the average of f values with all other parcels. Fidelity hence measures the parcel reconstruction accuracy and parcel cross-talk summarizes the degree of mixing a given parcel has with other parcels. Parcel spread (s) for each parcel p was calculated as the mean f of all other parcels with the original time series of parcel p and hence reveals parcels generating signals that are spuriously picked up by many other parcels. To decrease the probability of spurious synchronization and anatomically misplaced connections, for both statistical analyses and graph visualization in this study, we removed parcels with f cÁs <3000 (a.u.), so that~10% of the cortical surface was excluded. These parcels were located mostly deep and/or inferior sources, which generate the least detectable signals in M/EEG and hence are most likely to incorrectly reflect signals generated elsewhere.

1:1 and CF synchrony analysis
We estimated 1:1 phase-synchrony within frequencies and m:n phase-synchrony between frequencies from source-reconstructed phase-time series of M/EEG data (Figure 3-figure supplement 1f). The filtered time series Xðt; f Þ are complex and can be expressed as amplitude and a phase time series A x ðt; f Þ and x ðt; f Þ, respectively, so that We estimated n : m cross-frequency phase synchrony (CFS) and 1:1 within-frequency phase synchrony for each pair of cortical parcels p and q and between pairs of low and high frequencies f high and f low . We used the phase-locking value (PLV) to compute synchronization so that PLV p;q;n:m;flow;fhigh ¼ where iis the imaginary unit, the integers nand mdefine the frequency ratio so that n Á f high ¼ m Á f low , where N r is the number of trials r and N t is the number of samples t within a time window (Palva et al., 2005). For CFS we used n = 1 and m e {2,3,4,5,6,7,8,9} while for 1:1 within-frequency phase synchrony m = 1. Phase-amplitude coupling (PAC) was also quantified for the same frequency pairs as CFS. We estimated PAC by computing the 1:1 PLV between the phase of the slow oscillation and the phase of the amplitude envelope of the fast oscillation (Vanhatalo et al., 2004). The phase of the amplitude envelope was estimated by filtering it at f low . PAC was thus defined as: where E ðt; f low ; f high Þ is the phase of the filtered amplitude envelope time series Eðt; f low ; f hight Þ that was obtained by filtering Aðt; f high Þ with the Morlet wavelet wðt; f low Þ: Cross-frequency amplitude correlations were estimated for the same frequency pairs as the CFS by using the Pearson correlation coefficient CC of the amplitude time series A(t) so that CC p;q;flow;fhigh ¼

Single-trial analysis
To assess 1:1 synchrony and 1:4 CFS at the single-trial level (see Figure 2), we computed the 1:1 and 1:4 PLV for a single trial between parcels observed to exhibit clear 1:4 CFS in the group analyses (see below and Figure 4) in sliding 300 ms wide time windows at 50 ms distance. While 1:1 PLV estimates may be biased by residual linear mixing after source reconstruction, neither the long-range couplings nor the strengthening of 1:1 PLVs during the retention period can be attributed to linear mixing. To assess the statistical significance, we computed the 99th% -ile of the PLV distribution of 1000 realizations of surrogate data. The surrogate PLVs were assessed identically to those of real data and for the same sliding time windows, but so that one of the time series was shifted randomly by 0À300 ms. This method of surrogate generation is specific for measuring the significance of interareal correlations because all local features that bias the PLV estimates, such as temporal auto-correlations, are preserved.

Group-level statistics
We then computed m:n CFS, PAC, and 1:1 phase-synchronization across trials separately for all subjects (N = 12) between all cortical parcels, combinations frequency pairs of f low and f high , and for four time-windows; two before stimulus onset and for early (400-700 ms from Sample stimuli) and late (700-1000 ms from Sample stimuli) retention periods and for all memory loads l (N = 6).
Since PLV is sensitive to the number of trials, we equalized the number of trials between conditions (memory loads) within each subject to the lowest number of trials among conditions. Prior to statistics, the CFS adjacency matrices were collapsed (Figure 3-figure supplement 1g) to a coarser parcellation, the Destrieux atlas (Destrieux et al., 2010) of 148 parcels (see Figure 2-figure supplement 2).
We then estimated (Figure 3-figure supplement 1h) the statistical significance across the entire subject population of each pairwise interaction for each of the 1:m ratios. We first identified interactions that were significantly stronger or weaker in the retention period time windows than in the first pre-stimulus baseline by using a two-sided t-test with a statistical significance level p<0.05, and corrected for multiple comparisons for data averaged across the six memory loads (Mean condition). The effect sizes (Mean difference PLV ret -PLV BL divided by the standard deviation of the difference) and p-values for all frequencies and ratios are provided in external source data 1 (Siebenhühner et al., 2016). To correct for multiple comparisons, we removed as many of the least significant edges as predicted by the p-value (False discovery rate correction). For the effect of memory load (Load condition), we identified interactions whose strength was positively or negatively correlated with the memory load with a Spearman rank correlation test (p<0.05; FDR corrected). In the Load condition, the effect size is given by the correlation coefficients that together with the p values are provided in external source data 1 (Siebenhühner et al., 2016).
We then used graph theoretical notation (Bullmore and Sporns, 2009) to characterize the statistical CFS adjacency matrices at the graph level. Here, brain areas are the nodes, and connections of synchrony, the edges. We used connection density, K, that was defined to indicate the proportion of statistically significant inter-areal interactions from all possible interactions to indicate the extent of statistically significant synchronization separately for each frequency ratios of CFS (Figure 3), CF amplitude-amplitude correlations (Figure 3-figure supplement 7) and PAC (Figure 8). The effect sizes and p values for CFS, amplitude-amplitude correlations and PAC are provided in external source data 1, 2 and 7, respectively (Siebenhühner et al., 2016).
To compare the overlap between CFS and CF amplitude-amplitude coupling, we then computed the fraction of edges that were found significant in both CFS and CF amplitude-amplitude correlations to those that were significant in CFS only (Figure 3-figure supplement 8).

Estimation of consistency across ratios
To quantify how consistently the connection density K of CFS was increased or decreased across the 1:m frequency ratios for each f low < 18Hz, we introduced (Figure 3-figure supplement 1i) a consistency measure C: C l;flow ¼ K l;flow Á exp À P N m K m;l;flow À K l;flow N Á K l;flow ! where l denotes the memory load or test condition (Mean or Load), K m;i;flow is the connection density for CFS low: high at load l and K is the mean value for all K, N is the number of ratios (here 8). The C values were considered significant if they exceeded 95% -ile confidence limits estimated from 2000 realizations of surrogate data that were constructed by shuffling the Kvalues over frequencies. The same procedures described in this paragraph were also used for estimation of K and C for amplitude-amplitude (CC) and phase-amplitude (PAC) couplings.

Construction of summary graphs across ratios
We then performed (Figure 3-figure supplement 1j) graph visualization for CFS graphs in two stages. In the first stage, we created graphs that summarize the connections shared among harmonically consistent networks across frequency ratios. These summary graphs are provided in external source data 3-6 (Siebenhühner et al., 2016). In the second stage, the graphs were visualized so that the most central edges were selected and hierarchically clustered to attenuate the effect of spurious edges.
In the first stage, in order to identify CFS edges that were shared across ratios, edges from all ratios of interest were clustered into hyperedges where each hyperedge could contain edges from the adjacency matrices of multiple ratios. This edge clustering was performed so that an edge-edge adjacency matrix was created and the adjacency was defined to be the product of corresponding maximum parcel cross-talk c (see previous paragraph). Thus each element in the edge-edge adjacency matrix indicates the proximity of any two observed edges in terms of linear mixing attributable to data acquisition and inverse modeling on the group level. Here high adjacency indicates a high probability that these edges actually represent the same underlying cortico-cortical interaction. Hierarchical clustering of the edge-edge adjacency matrix thus groups all observed edges into a set of hyperedges that are mutually dissociable in the MEG linear mixing space. To achieve a representation of connectivity shared by multiple frequency ratios, we used permutation statistics to test for each hyperedge whether it could be significantly (p<0.05, uncorrected) attributable to edges from any single ratio. Then those hyperedges for which this was the case were excluded from further visualization. The remaining hyperedges were then 'flattened' into a regular adjacency matrix by summing the shared edges across ratios.

Bundling of edges into hyperedges
In the second stage ( Figure 3-figure supplement 1k), in order to visualize these graphs for the Mean and Load conditions we picked the 500 and 250 most central edges, respectively, i.e., edges that were connected to the most central low-or high-frequency hubs. For f low , we ranked the edges by the low-frequency degree centrality of the vertices to which their low-frequency end was connected. Conversely for f high , we ranked edges by the high-frequency degree centrality of the vertices that their high-frequency end was connected to.
We then re-used the edge-adjacency based hierarchical clustering described above in order to bundle the edges into clusters that are likely to each represent a single underlying true interaction. A bundle per se indicates the residual uncertainty about the anatomical localization of the interaction. As large bundles are much more likely to reflect major true neuronal interactions than small bundles, we excluded bundles of less than 4 edges and discarded edges with low centrality (<50% ) within a bundle.
We created low-and high-frequency graphs, i.e., graphs that showed the most central LF or HF hubs and their connections to other vertices (parcels), respectively. Circles were used to denote the vertices, i.e., parcels, so that their size was proportional to their degree in the graph. The vertices were coloured by the functional brain system (Yeo et al., 2011) they belonged to. As many vertices would be both at the high and low-frequency 'end' of connections within a graph, we chose to indicate their fractions of LF and HF degree by the light and darker shaded sectors. For example, a vertex with 3 connections where it was at the LF 'end' and one where it was at the HF 'end' would be divided into 75% area of lighter and 25% area of darker shade. CFS connections were indicated by lines where dashed ends pointed towards the high-frequency vertices in LF graphs and towards lowfrequency vertices in HF graphs. The neuroanatomical spread of the ends of each bundle illustrates the level of spatial smearing in M/EEG source reconstruction used here. The graphs were overlaid on flattened maps of the complete cortical surfaces of left and right hemispheres. These maps indicate the cortical gyration so that sulci are darker than gyri. Color in flattened brains identifies the 7 brain systems of the Yeo parcellations (Yeo et al., 2011) and the thin white or grey lines the 148 parcels of the Destrieux parcellation (Destrieux et al., 2010). These graphs are shown in Figures 4, 5 and their supplements.

Assessing interactions within and between functional brain systems
To estimate the interactions between functional subsystems, we morphed ( Figure 3-figure supplement 1l) the ratio-collapsed graphs graphs from 148x148 parcel interaction matrices into 7x7 subsystem interaction matrices, based on the Yeo 2011 parcellation scheme (Yeo et al., 2011) and computed for each element the connection density K. We then computed 5000 randomized versions of the same 148x148 matrix, keeping the number of edges constant. K values of the true subsystem interaction matrix were reported as significant (see Figure 6, Figure 6-source data 1), if they exceeded the 95th percentile of K values in randomized graphs.
Correlation between networks of CFS and 1:1 phase synchrony To assess whether the CFS connected 1:1 synchronized networks, we obtained (Figure 3-figure supplement 1) Mean condition 1:1 networks for all analyzed frequencies with the statistical approach identical to what was used for CFS. With these networks, we tested whether the vertex centralities (represented by vertex degree) of the low-frequency (f low ) CFS network vertices were correlated with centrality values for the vertices of the corresponding slow 1:1 networks. Conversely, we tested whether the vertices of the high-frequency (f high ) part of the CFS networks were correlated with the vertices of fast 1:1 networks. Finding such correlations would indicate that CFS connects the major nodes of spectrally distinct within-frequency synchronized networks. We used the vertex degree to quantify the centrality and Pearson correlation to quantify the correlations. The Pearson correlation p values were corrected for multiple comparisons with the Benjamini-Hochberg FDR method (Figure 7a

Degree correlation between frequencies
To assess the parcel centrality similarity between narrow-band 1:1 networks, we computed the vertex degree for each parcel in the Mean condition (positive tail) networks of all 25 narrow-band frequencies. We then evaluated, for each pair of frequencies, the correlation of vertex degrees across all parcels with the Pearson's correlation coefficient. The results are presented in Figure 7-figure supplement 3.

Correlation of CFS and PAC with the behavioral performance
To assess whether CFS and PAC were predictive of the subjects' behavioral VWM performance, we estimated ( Figure 3-figure supplement 1n) the correlation between the network strength of CFS networks and individual capacity for each memory load. Individual CFS strength was obtained as the sum of individual connection strengths across connections significant in the positive tail of Mean condition group statistics. Individual VWM capacity was defined for each VWM load as C(l) = HR(l) . l, where HR(l) is the hit rate at load l.
In detail, the interaction matrices for the Mean condition were used to create binary masks M. If the interaction at the 1:m-ratio m, at low frequency f low , in time window t was found significant between parcels pand q in the group-level analysis, Mðp; q; m; t; f low Þ was set to 1, otherwise to 0. For each subject, we multiplied CFS and PAC adjacency matrices with these masks and then summed over all parcel pairs. Subjects' individual network strength S for each m, memory load l, and frequency f low was thus calculated as: S 0 l; s ð Þ ¼ 1 Sabs S l; s ð Þ À S L s ð Þ À S pop l ð Þ À Á .
Identical detrending and normalization was applied to the VWM capacity values, i.e., for obtaining C 0 ðs; lÞ from Cðs lÞ. We then averaged S 0 separately across low (1:2À1:5) and high (1:6À1:9) CFS ratios and then computed for each f low the correlation of S' and C' across subjects and VWM loads with the Pearson correlation test (see Figure 9a and c, Figure 9-source data 1). Correction for multiple comparisons across all f low ðn ¼ 12Þ was performed using the Benjamini-Hochberg method. For dissecting post hoc the role of 1:2 CFS (see Results), we further also evaluated these correlations for the ratios 1:2 and 1:3À1:5 for f low ¼ 13 Hz. For visualization without statistics (see Figure 9b and d, Figure 9-source data 2), we also computed correlation coefficients separately for individual ratios.

Leave-one-out and effect size analyses
We performed leave-one-out and effect-size analyses to confirm that the findings of increased CFS (Mean condition) and of CFS positively correlated with VWM Load (Load condition) were not affected by the small sample size of our cohort. In the leave-one-out analyses, we created 12 datasets of 11 subjects, where in each one a different subject had been excluded and repeated the statistical analyses underlying Figure 2a and c. We then retained those edges that were found significant in the original analyses and had a p-value<0.05 in each leave-one-out dataset. These results are presented in Figure 3-figure supplement 2a and d for Mean and Load conditions, respectively. We then estimated that with 12 subjects in our cohort, to get alpha level 0.05 with statistical power of 0.8, our data needed to have an effect size of minimum of 0.9 for the Mean condition and 0.35 for the Load condition. In the effect size analyses, we thus retained those edges which were found significant in the original analyses and which had an effect size > 0.9 (Mean condition) or a correlation coefficient > 0.35 (Load condition). (Figure 3-figure supplement 2b and e). Finally, we show the combined results from these masks in Figure 3-figure supplement 2c and f, where only the edges significant in both the leave-one-out and effect-size analyses are retained. The same analysis was also done for PAC, with the results shown is Comparison of observed changes in PLV with changes predicted from changes in SNR or oscillation amplitudes Changes in phase synchrony may be partially caused by changes in oscillation amplitude. In particular, an increase in signal-to-noise ratio (SNR) can improve phase estimation and therefore also phase-locking values (PLV). However, we showed in our earlier analysis of the same dataset that changes in SNR could only explain a small fraction of the observed changes in PLV (Palva et al., 2010). We here extend the approach to the case of cross-frequency phase synchrony.
We simulated four time series X, Y, N x and N y as white noise with 10 6 samples in the range of -1 to 1 that was filtered with a Morlet wavelet at f 1 for X and N x and at f 2 for Y and N y , where f 2 ¼ n Á f 1 , obtaining complex time series that can be written in the form X ¼ A x Á exp i Á ' x ð Þandanalogously. We then obtained coupled time series as follows: Y', N x and N y were then normalized as in X 0 n ¼ X 0 =jX 0 j and so on, and the noise added to the coupled time series series: X 00 ¼ X 0 n þ s x Á N x;n and Y 00 ¼ Y 0 n þ s y Á N y;n , where s ¼ 1=SNR. Then, we estimated CFS, using PLV, between X'' and Y' for a range of values each for s x and s y , which were spaced so that at SNR i1 ¼ 1:3 Á SNR i . We estimated the apparent SNR (aSNR) by comparing the mean amplitude levels for experimental data across object loads and empty-room recordings that were obtained before or after each recording session. Empty-room data were inverse modeled with the inverse operator derived from the experimental data onto 400 parcels, and we defined the aSNR as follows: aSNR ¼ where A exp; s; p denotes the mean amplitude for experimental data and A ER; s; p the mean amplitude for empty room recordingsfor subject s (here N s ¼ 6) and parcel p. The values for aSNR are shown in Figure 3-figure supplement 3. The (non-linear) function linking aSNR and SNR was estimated numerically. We created white noise time series, noise and signal plus noise time series as above and obtained aSNR ¼ ASN ÀAN AN , where A SþN is the mean amplitude of the signal plus noise time series and A N the mean amplitude of the noise time series. By interpolating this function, we estimated SNR values for experimental data from aSNR values.
We used these results to predict how large changes in PLV would be caused by changes in amplitude in the low and high frequency band for selected frequency ratios (those that showed the largest K values). For each of these, we identified the 200 parcel pairs with the largest observed change in PLV in order to focus the analysis on the strongest effects in this dataset. If PLV were significantly biased by SNR changes, this interaction should be most salient in the strongest PLV findings.
To obtain the changes of PLV in experimental data in the Mean condition, we computed for each parcel pair j, k, the difference of the mean PLV between initial and modulated condition. For the Mean effect, we computed the difference of the mean over subjects and object loads l in the retention period with that in the baseline window: DPLV exp;j;k;f1;f2 ¼ 1 N s Á N l X s X l PLV s;j;k;f1;f2;ret À PLV s;j;k;f1;f2;BL For the Load effect, we computed instead the difference of the 6 object and 1 object conditions in the retention period: DPLV exp;j;k;f1;f2 ¼ 1 N s X s PLV s;j;k;f1;f2;ret l ¼ 6 ð Þ À PLV s;j;k;f1;f2;ret l ¼ 1 ð Þ We further estimated the SNR for these conditions and parcel pairs from the aSNR and used the SNR levels and PLV for the initial condition to estimate the coupling factor c from the simulated time series, for which: PLV exp;init f 1 ; f 2 ð Þ ¼ PLV sim;f1;f2 c in ; SNR 1;init ; SNR 2;init À Á We then estimated the predicted PLV from the simulated time series with the SNR levels in the modulated condition: PLV pred f 1 ; f 2 ð Þ ¼ PLV sim;f1;f2 c in ; SNR 1;mod ; SNR 2;mod À Á and thus obtained the change in PLV that could be expected from amplitude changes alone: The results for the Mean and Load conditions are shown in Figure 3-figure supplements 4 and 5, respectively. To assess whether PLV and amplitude effects were correlated directly, we calculated the Meanand Load-condition effects identically to what was done in the main analyses, for each ratio, frequency and parcel pair. In the Mean condition, we computed the average PLV difference between retention period and baseline time windows over subjects and conditions. In the load condition, we correlated the PLV values in the retention period with the object load across subjects. We then morphed these data from the 400-parcel to the coarser 148-parcel Destrieux parcellation. In the same manner, we computed for each frequency and ratio the Mean and Load effects for the measured parcel amplitudes, averaged these amplitude effects for each parcel pair, and morphed to the Destrieux parcellation. For both Mean and Load condition, we then computed the correlation coefficient r between the PLV and amplitude effect for all parcel pairs that were found to have a significant PLV effect in the corresponding group analysis. We compared these r against 1000 surrogate datasets where the amplitude values were randomly shifted across edge pairs. We considered the PLVeffect vs. amplitude-effect correlation significant if the r was in the highest or lowest 2.5%-ile of the surrogate values. No multiple comparison corrections were applied across the studied frequencyratio pairs. In Figure 3-figure supplement 6, we illustrate the fraction of variance in PLV effects explained by the amplitude effects with r 2 values for those ratios and frequencies that had a significant correlation between the PLV and amplitude effects and initially had K > 0.1% in the group analysis. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Additional information
Author contributions FS, SHW, Analysis and interpretation of data, Drafting or revising the article; JMP, Conception and design, Analysis and interpretation of data, Drafting or revising the article; SP, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article Author ORCIDs Satu Palva, http://orcid.org/0000-0001-9496-7391

Ethics
Human subjects: This study was approved by the ethical committee of Helsinki University Central hospital and was performed according to the Declaration of Helsinki. Written informed consent was obtained from each subject prior to the experiment.

Major datasets
The following dataset was generated: