Abstract
Spatially distributed phaseamplitude coupling (PAC) is a possible mechanism for selectively routing information through neuronal networks. If so, two key properties determine its selectivity and flexibility, phase diversity over space, and frequency diversity. To investigate these issues, we analyzed 42 human electrocorticographic recordings from 27 patients performing a working memory task. We demonstrate that (1) spatially distributed PAC occurred at distances >10 cm, (2) involved diverse preferred coupling phases, and (3) involved diverse frequencies. Using a novel technique [Nway decomposition based on the PARAFAC (for Parallel Factor analysis) model], we demonstrate that (4) these diverse phases originated mainly from the phaseproviding oscillations. With these properties, PAC can be the backbone of a mechanism that is able to separate spatially distributed networks operating in parallel.
Introduction
In a fastchanging complex environment, it is essential that the brain can selectively route information through multiple networks operating in parallel. Oscillatory coupling provides the temporal and spatial dynamics necessary to implement this. Oscillations reflect rhythmic modulations of neuronal excitability, affecting the efficacy of incoming EPSPs and the probability of spike output. It has been proposed that coherently oscillating networks create selective windows of communication between neuronal groups by synchronizing their periods of maximum excitability (Fries, 2005). We investigated spatially distributed phaseamplitude coupling (PAC), a phenomenon that may emerge in spatially distributed oscillating networks. We observed two key properties, phase diversity and frequency diversity, that allow spatially distributed PAC to flexibly and selectively route information through distributed neuronal networks.
PAC describes the coupling between the phase of a slow oscillation and the amplitude of a fast oscillation, with the highest amplitude occurring at the socalled preferred coupling phase (Jensen and Colgin, 2007; Lakatos et al., 2008; Canolty and Knight, 2010). PAC has been observed in multiple species, including rats (Chrobak and Buzsáki, 1998; Sirota et al., 2008; Tort et al., 2008), monkeys (Lakatos et al., 2005, 2008), and humans (Schack et al., 2002; Bruns and Eckhorn, 2004; Mormann et al., 2005; Canolty et al., 2006; Cohen, 2008; Osipova et al., 2008; Miller et al., 2010; Voytek et al., 2010; Maris et al., 2011). A recent study has shown that PAC in human electrocorticography (ECoG) is widely spatially distributed (Maris et al., 2011), which is a key requirement for routing information through distributed networks. Using a novel decomposition technique, this study showed that amplitude and phaseproviding oscillations occurred at broadly distributed sites. However, the decomposition only allowed for diversity over space in the preferred coupling phases of the amplitudeproviding oscillations and not of the phaseproviding oscillations. Selective routing of information could greatly benefit from phase diversity in the phaseproviding oscillation, because its phase could be used to select neuronal populations for interactions. Such phase diversity has not been shown so far. Frequency diversity is another key property determining the flexibility of PAC in selective routing of information. Although theta–gamma PAC dominates the literature, recent reports (He et al., 2010; Miller et al., 2010; Maris et al., 2011) have shown that PAC occurs at many different frequencies.
We analyzed ECoG recordings from 27 patients to investigate the phase diversity and frequency diversity in spatially distributed PAC. We show that PAC occurred over distances that exceed 10 cm, that there was strong phase diversity, and that it involved diverse frequencies. Using a modified version of Nway decomposition based on the PARAFAC (for Parallel Factor analysis) model (Maris et al., 2011), we were able to show that spatially distributed phaseproviding oscillations were the main source of phase diversity. These oscillations showed large and consistent phase diversity over space. In contrast, the amplitudeproviding oscillations showed bursts that were much more synchronized. This phase and frequency diversity are two important attributes that determine how flexibly and selectively spatially distributed PAC can route information through distributed networks operating in parallel.
Materials and Methods
Subjects.
Twentyseven patients (12 female, 15 male) with pharmacoresistant epilepsy were implanted with subdural grid, strip, and depth electrodes before resective surgery. Patients were selected from a large pool of datasets if they had >15 electrodes and >70 trials per recording session after artifact rejection. Informed consent was obtained from the patients or their guardians if they were underage. The research protocol was approved by the appropriate institutional review boards at the Hospital at the University of Pennsylvania (Philadelphia, PA), Children's Hospital (Philadelphia, PA), University Clinic (Freiburg, Germany), Children's Hospital (Boston, MA), and Brigham and Women's Hospital (Boston, MA). Some of the datasets have been analyzed previously, but the analyses presented here are novel (Rizzuto et al., 2003; Raghavachari et al., 2006; Jacobs and Kahana, 2009; van Vugt et al., 2010) or complementary (Maris et al., 2011).
Experimental paradigm.
Recordings were obtained from patients performing a Sternberg working memory task (Sternberg, 1966). Patients were presented with a series of letters (from one to six) on a computer screen that they had to remember. At the beginning of each trial, a fixation cross was presented, followed by 700 ms of letter presentation, and then by 275–350 ms (uniformly distributed) of blank screen. The last letter was followed by a retention interval of 425–575 ms (uniformly distributed), after which a probe letter was presented. Patients had to indicate by key press whether the probe letter was part of the previous letter series. After the key press, visual feedback was given and the patient could initiate the next trial by another key press. The main purpose of our study was to characterize fundamental properties of PAC (spatial distribution, phase diversity, and frequency diversity), and therefore we did not investigate any behavioral contrasts (e.g., correct vs incorrect) or stimulus type contrasts (e.g., number of letters). We only analyzed the period between the fixation cross and the onset of the probe letter, during which patients were actively engaged.
Recordings and preprocessing.
Electrophysiological recordings were obtained from subdural grid, strip, and depth electrodes. Recordings were sampled at 256–1024 Hz, depending on the hospital, and were referenced to a common average reference. Note that using a nearestneighbor bipolar referencing scheme did not substantially change the electrode pairs that showed strong PAC or the level of diversity in the preferred coupling phases. Only recordings from grid and strip electrodes were analyzed. Artifact rejection was performed by visual inspection. All trials and/or electrodes contaminated by epileptiform activity were removed. To remove power line noise, we bandstop filtered the data with 1 Hz windows at 50 and 60 Hz (depending on continent) and at other frequencies containing line noise. All recordings were bandpass filtered between 0.01 and 100 Hz. All filters were fourthorder Butterworth.
Electrode locations were determined by first coregistering a postoperative computed tomography scan with a higherresolution preoperative magnetic resonance image. All patients' brains were normalized to Talairach space (Talairach and Tornoux, 1988), and coordinates were subsequently computed. All preprocessing and the first step of our spectral analyses were performed using the FieldTrip opensource MATLAB toolbox (Oostenveld et al., 2011) developed at the Donders Institute for Brain, Cognition, and Behavior (http://www.ru.nl/neuroimaging/fieldtrip).
Amplitudeweighted phaselocking factor.
To quantify PAC, we calculated amplitudeweighted phaselocking factors (wPLFs). These coefficients were calculated using the output of a timeresolved spectral analysis. This spectral analysis was performed by convolving the data with complexvalued wavelets, one for every frequency of interest. All wavelets were obtained from an elementwise multiplication of a threecycle complex exponential and a Hanning taper of equal length. For every given sampling rate, we only used frequencies for which the corresponding wavelet has an integer number of samples per cycle. Under this constraint, we selected frequencies between 2 and 67 Hz in steps of ∼1 Hz. This resulted in one complexvalued time series per trial per frequency bin, called the wavelet transform, describing the timevarying amplitudes and phases.
A wPLF is a complexvalued number representing the relation averaged over time between the phase of one oscillation (obtained from electrode k at frequency m) and the amplitude of another (obtained from electrode j at frequency l). Thus, a wPLF is indexed by an electrode pair (indices k and j) and a frequency pair (indices m and l). A wPLF is normalized, with magnitude ranging from 0 to 1. The magnitude of a wPLF measures the consistency, over trials, of the phase of the phaseproviding oscillation at which amplitude increases of the amplitudeproviding oscillations occur. The angle of a wPLF indicates this phase, called the preferred coupling phase. A wPLF is amplitude weighted because trials with high amplitudes in one oscillation have a bigger influence than trials with low amplitudes.
The calculation of the wPLFs can be expressed as follows: in which A(x_{j}, f_{l}) denotes the meancentered absolutevalue of the wavelet transform at frequency l of the raw signal x of electrode j, W(x_{k}, f_{m}) denotes the wavelet transform at frequency m of the raw signal x of electrode k, 〈 , 〉 denotes the inner product, and ‖ ‖ denotes the norm. wPLFs were computed for all possible electrode and frequency pairs. Thus, each electrode in the dataset provides amplitude and phase information for PAC with all other electrodes and for all estimated frequencies. This results in one 4way array of wPLFs for each dataset. The dimensions correspond to (1) amplitudeproviding electrodes (of size J), (2) phaseproviding electrodes (of size K), (3) amplitudeproviding frequencies (of size L), and (4) phaseproviding frequencies (of size M). We show a schematic of the construction of this 4way array in Figure 1, A and B.
Selecting significant wPLFs.
As a part of the analysis of the 4way arrays of wPLFs, statistically significant wPLFs were selected. Statistical significance was assessed by comparing every wPLF to a datasetspecific reference distribution obtained under the null hypothesis that the timevarying amplitudes and phases are uncorrelated. Reference distributions were created by randomly pairing the amplitudes of one trial with the phases of another trial. This was repeated 50 times for each dataset, providing 50 random wPLFs for each electrode pair and frequency pair. A normal probability density function was then estimated for every wPLF, using the mean and SD of the magnitude of these 50 random wPLFs. wPLFs were selected if their magnitude surpassed the 99th percentile of this estimated probability density function. Apart from this selection based on a statistical threshold, we also removed all wPLFs in which the phaseproviding frequency is higher than or equal to the amplitudeproviding frequency.
Evaluating the reliability of preferred coupling phases.
To evaluate the reliability of the preferred coupling phase of our significant wPLFs, we used a splithalf procedure. This involved a random split of the trials of each dataset in two partitions, followed by constructing a 4way array of wPLFs for both partitions. In this way, we obtained two independent estimates of every wPLF. The more the preferred coupling phase is influenced by random noise, the larger the phase difference will be between the two estimates. Based on these splithalf wPLFs, we calculated a splithalf reliability coefficient: In this formula, we first take the difference between the preferred coupling phase of the wPLFs of the first partition θ_{1j} and of the second partition θ_{2j}. We do this for all J significant wPLFs, as determined above (using the 4way array of wPLFs based on all trials). These phase differences are then expressed as unitmagnitude complex numbers and averaged over all significant wPLFs (indexed by j). The splithalf reliability coefficient is then attained by taking the absolute value of the resulting complex number, also known as the mean resultant vector of phase differences. Such a coefficient was calculated for each dataset.
Nway decomposition based on the PARAFAC model.
We use Nway decomposition to refer to the decomposition of an Nway array with more than two dimensions. Unlike 2way decompositions, such as principal and independent component analysis, Nway decomposition has not been used extensively in neuroscience (for exceptions, see Beckmann and Smith, 2005; Mørup et al., 2006). Nway decomposition decomposes an Nway array into several components, each consisting of N loading vectors, one corresponding to each dimension of the original array. Every component describes one aspect of the array, and the original array can be reconstructed from its components. There are multiple models for decomposing an Nway array, but here we will only describe and use the most parsimonious decomposition, which is based on the PARAFAC model (Harshman, 1970), also known as CANDECOMP (for Canonical Decomposition) (Carrol and Chang, 1970). This decomposition can be derived from a few plausible assumptions about the spatiospectral characteristics of the sources that are involved in PAC (Maris et al., 2011). Crucially, for Nway arrays with more than two dimensions (N > 2), Nway decomposition based on the PARAFAC model is unique up to scaling and permutation, which are two transformations that do not affect the interpretation of the components. The Nway PARAFAC algorithm we used (see below) will be implemented in the FieldTrip opensource MATLAB toolbox (Oostenveld et al., 2011) developed at the Donders Institute for Brain, Cognition and Behavior (http://www.ru.nl/neuroimaging/fieldtrip).
Nway decomposition of 4way arrays of wPLFs into two complexvalued spatial maps and two frequency profiles.
We used Nway decomposition based on the PARAFAC model to decompose each 4way array of wPLFs into one or more components. We show this schematically in Figure 1C. Each component consists of four loading vectors, one for each dimension. Because the first two dimensions of a 4way array of wPLFs correspond to amplitude and phaseproviding electrodes, the corresponding loading vectors in a component describe spatial locations. We denote these two loading vectors as amplitudeproviding and phaseproviding spatial maps. Furthermore, because the last two dimensions of a 4way array of wPLFs correspond to amplitude and phaseproviding frequencies, we denote these as amplitude and phaseproviding frequency profiles. Each component thus describes a PAC pattern that is characterized by an amplitude and phaseproviding spatial map and an amplitude and phaseproviding frequency profile.
The decomposition of a single wPLF can be expressed in a formula involving elementwise multiplication: A wPLF_{jklm} is described as the sum, over components f, of the product of the loadings a_{jf}, b_{kf}, c_{lf}, and d_{mf}. These loadings are organized in the loading matrices A, B, C, and D, respectively. Matrices A and B contain as columns the amplitude and phaseproviding spatial maps, and matrices C and D contain as columns the amplitudeand phaseproviding frequency profiles. The spatial maps A and B are complex valued, whereas frequency profiles C and D are real valued. This differs from the previous approach (Maris et al., 2011), in which only the amplitudeproviding spatial map (A) was complex valued, reflecting the assumption that there are no betweenelectrode phase differences in the phaseproviding oscillation over electrodes (except for phase differences of exactly ±π, which are translated into loadings that have different signs). For our current approach, investigating phase diversity in PAC, it is essential that the phaseproviding spatial map B is complex valued as well.
Decomposition of preferred coupling phases in PAC into relative phases in two spatial maps.
Our Nway decomposition decomposes all preferred coupling phases in spatially distributed PAC into two complexvalued spatial maps. As such, phase diversity in PAC is fully explained by the phase relations within the two spatial maps. Phase diversity in the phaseproviding spatial map reflects consistent phase differences of the phaseproviding oscillation over electrodes. Phase diversity in the amplitudeproviding spatial map reflects time delays between amplitude increases of the amplitudeproviding oscillation. The exact time delay depends on the cycle length of the phaseproviding oscillation. In Figure 2, we show a schematic of this decomposition. This schematic shows PAC between five electrodes (Fig. 2A) and their decomposition into an amplitudeproviding (Fig. 2B) and a phaseproviding (Fig. 2C) spatial map. The phase relations within the amplitudeproviding spatial map reflect betweenelectrode time delays between bursts of the amplitudeproviding oscillations (Fig. 2B). The phase relations within the phaseproviding spatial map reflect the betweenelectrode phase differences in the phaseproviding oscillations (Fig. 2C). Note, we cannot distinguish between (1) the case in which every cycle of the phaseproviding oscillation shows a burst of the amplitudeproviding oscillation (Fig. 2A, first row of table) and (2) the case in which only some cycles show such a burst (Fig. 2A, second row of table).
Indeterminacies of the PARAFAC model.
The spatial maps and frequency profiles can only be determined up to scaling and permutation. Because of permutation indeterminacy, the order of components is irrelevant, and because of scaling indeterminacy, any loading vector of the same component can be multiplied with any number, as long as another loading vector of the same component is multiplied with the inverse of this number. Moreover, because the two spatial maps are complex valued, there is phase indeterminacy. If one spatial map is multiplied with e^{−iθ} and the other is multiplied with e^{i}^{θ}, then all phases are shifted by θ, yet the decomposition remains exactly the same. Note that this does not affect the phase differences within a component. Because of the above, the components in our decomposition are sorted by explained variance, all loading vectors are normalized to have a norm of 1, and all spatial maps have an average magnitudeweighted phase of 0. This means that absolute phases inside a spatial map cannot be interpreted. To stress that, within a component, only betweenelectrode phase differences can be interpreted, the phase of the spatial maps will be denoted as relative phases. Analogously, the magnitude of each electrode in a spatial map and the value of each frequency in the frequency profiles can only be interpreted relative to the other electrodes in the map and the other frequencies in the profile, respectively.
Reconstructing 4way arrays of wPLFs and evaluating their accuracy.
Using the extracted components, we can reconstruct each 4way array of wPLFs. We computed reconstructions based on all components to evaluate the accuracy of the Nway decomposition for every dataset. We also computed reconstructions based on a single component to select significant electrodes in a spatial map. We show both reconstructions schematically in Figure 1D. The reconstruction of wPLF x_{jklm} is denoted by x̂_{jklm} and, when based on all components, it is calculated as follows: Thus, the wPLF at amplitude and phaseproviding electrodes j and k and at amplitude and phaseproviding frequencies l and m can be reconstructed by taking the product a_{jf}b_{kf}c_{lf}d_{mf} and summing over the components. When based on a single component f, x̂_{jklm} is equal to the product a_{jf}b_{kf}c_{lf}d_{mf}.
We evaluated the accuracy of reconstructed wPLFs (based on all components) by a coefficient comparing the reconstructed wPLFs to the observed wPLFs. This coefficient was calculated as follows: We took the inner product 〈 , 〉 of the complex conjugate of the vectorized 4way array of observed wPLFs v̅e̅c̅(̅X̅)̅ and the vectorized 4way array of reconstructed wPLFs vec(X̂). This inner product is then normalized by the product of the vector norms, and the absolute value is taken. This number ranges from 0 to 1. It can be interpreted as a generalized cosine of the angle between two complexvalued vectors.
Reporting on the results of the Nway decomposition.
To report the decomposition results, we selected electrodes on the basis of componentspecific reconstructed wPLFs. This selection is necessary because we are interested in the phases of the two types of spatial maps, and these phases can only be reliably estimated for electrodes that are involved in the component. We performed our electrode selection by first reconstructing a 4way array of wPLFs on the basis of a single component (see Materials and Methods). Next, this array of reconstructed wPLFs was compared with the reference distribution from the same dataset (see above, Selecting significant wPLFs). When more than one component was extracted from an array of wPLFs, the same reference distribution was used multiple times. An electrode in the amplitudeproviding spatial map was selected if any of the reconstructed wPLFs that have this electrode as the amplitudeproviding electrode exceeded the 99th percentile of the reference distribution. The same criterion was applied to the phaseproviding spatial maps.
We report on phase differences within each spatial map by computing the difference between pairs of selected electrodes. Among others, we show phase differences occurring between 0 and ±π. To demonstrate that these phase differences are reliable (i.e., reflecting true phase differences between 0 and ±π), we calculated the splithalf reliability of our decomposition results. This involved randomly splitting the trials of each dataset and calculating two 4way arrays of wPLFs, one for each of the two sets of trials. We then decomposed both arrays into the same number of components and computed the phase differences between pairs of selected electrodes (see above). In addition to this selection based on statistical significance, we selected electrode pairs with phase differences in the intervals from −2π/3 to −π/3 and from π/3 to 2π/3. For these electrode pairs, we calculated the following splithalf reliability coefficient: In this formula, we first take the difference between the two independent estimates of the phase difference for electrode pair j, one obtained from the first (θ_{1j}) and the other from the second (θ_{2j}) partition. We do this for all J selected electrode pairs. All phase differences are then expressed as unitmagnitude complex numbers and averaged, producing a mean resultant vector. The splithalf reliability coefficient is then obtained by taking the magnitude of this mean resultant vector.
An alternating leastsquares algorithm for Nway decomposition of a 4way array of wPLFs.
Nway decomposition according to the PARAFAC model can be performed using an alternating leastsquares (ALS) algorithm that has been implemented for realvalued arrays (Bro, 1998) and for complexvalued arrays (Sidiropoulos et al., 2000). The algorithm for complexvalued arrays produces only complexvalued components. In contrast, in our application, we decompose a complexvalued array into components that consist of two complexvalued spatial maps and two realvalued frequency profiles. We now describe the algorithm and how we adapted it for complexvalued arrays for our application.
The ALS algorithm is an iterative algorithm with, per iteration, as many steps as the number of different loading matrices. In our case, in every stage of an iteration, a loading matrix leastsquares estimate is calculated while keeping the other loading matrices constant. The algorithm continues until an iteration does not provide an increase in fit over and above the previous iteration. All loading matrices are initialized by random starting values, which are orthogonal over components. The algorithm can converge to a suboptimal solution, which is a local minimum of the leastsquares objective function that we want to minimize. This is undesirable but can be controlled for by running the algorithm many times with random starting values. If the algorithm converges multiple times to the same solution using different random starting values and this solution also achieves the smallest objective function value, then it is assumed to have converged to the global minimum. It is crucial to detect and discard degenerate models that occur when component pairs are nearly identical but negatively correlated (Bro, 1998). To perform an Nway decomposition based on the PARAFAC model, it is also necessary to estimate the number of components, or the socalled rank of the array. We determined this rank using a splithalf procedure, identical to the procedure by Maris et al. (2011).
To describe our ALS algorithm, it is convenient to make use of the Khatri–Rao product ⊗, which is defined as follows: This applies to any matrix A and B with an equal number of columns F. The Khatri–Rao product is defined as the concatenation of the Kronecker tensor products ⊗ of column 1 to F of A and B. Using the Khatri–Rao product, we can express the 4way PARAFAC model as follows: In this formula, X^{J·KLM} denotes a 2way array that is obtained by unfolding the 4way array X along its last three dimensions. The formula expresses that X^{J·KLM} is the sum of a model term A(D ⊗ C ⊗ B)^{T} and an error term E^{J·KLM}. The model term is a function of the loading matrices A, B, C, and D with dimensions J × F, K × F, L × F, and M × F, respectively. F, the number of columns in each loading matrix, denotes the number of components being extracted. In our application, J = K (the number of electrodes), and L = M (the number of frequencies). The error term E^{J·KLM} is necessary to express the fact that the observed wPLFs may differ from the model wPLFs (which are determined by the loadings) as a result of sampling error.
Each leastsquares estimate is calculated using the following equation. Because of the symmetry between the four loading matrices, we only present the estimation equations for one loading matrix, which we denote by A. A single iteration of the algorithm estimates all loading matrices once and then determines the fit. By keeping the loading matrices B, C, and D fixed, the leastsquares estimation for loading matrix A is the following: where Z = (D⊗C⊗B), Z̄ denotes the complex conjugate of Z, Z* denotes the complex conjugate transpose of Z, and + denotes the Moore–Penrose pseudoinverse. The leastsquares estimate of a realvalued loading matrix (matrices C and D) is obtained by replacing Z and X by the realvalued matrices Z = [Re(Z)Im(Z)], which is the rowwise concatenation of Re(Z) and Im(Z) (the real and the imaginary parts of Z) and X = [Re(X) Im(X)], which is the columnwise concatenation of Re(X) and Im(X).
Results
We analyzed PAC in human ECoG recordings from 42 datasets obtained from 27 patients performing a working memory task (see Materials and Methods). We investigated the phase diversity of spatially distributed PAC by means of wPLFs. These wPLFs are complexvalued association measures, quantifying coupling between the phase of one oscillation and the amplitude of another, averaged over time (see Materials and Methods). Our wPLF is a correlational and not a causal measure. To reflect this, we use the causally neutral terms “phaseproviding” and “amplitudeproviding” to denote the first and the second oscillation, respectively. wPLFs were computed for all amplitude and phaseproviding electrodes and frequencies. This results in a 4way array of wPLFs. As an example, we show one slice of such a 4way array (Fig. 3A), containing the wPLFs for all frequency pairs and a single electrode pair. For this electrode pair, the strongest coupling is between the phase of a theta oscillation (center frequency, 6 Hz) and the amplitude of a beta/lowgamma oscillation (center frequency, 23 Hz). The preferred coupling phase is π, which corresponds to the trough of the theta oscillation.
We used two approaches to investigate phase diversity in spatially distributed PAC, one based on selecting significant wPLFs from the 4way arrays (one array for every dataset) and the other based on a decomposition of each of these arrays. We first report on the results obtained by selecting significant wPLFs, showing that PAC occurred over long distances with substantial phase diversity. Next, we report on the decomposition results, showing that this phase diversity originated mainly from the spatially distributed phaseproviding oscillations.
PAC occurred over long distances, involved diverse preferred coupling phases, and involved many frequencies
We selected statistically significant wPLFs from each of the 42 datasets. Significance was assessed by comparing every wPLF with a reference distribution obtained under the null hypothesis of independence of phases and amplitudes (see Materials and Methods). wPLFs were selected if their magnitude exceeded the 99th percentile of this distribution. On average, 15.7 ± 9.0% (SD) of the wPLFs were selected from each of the 42 datasets. These were combined into one large data array used for all analyses on significant wPLFs. The contribution of each of the 42 datasets to this data array was on average 2.4 ± 2.3% (SD), and the contribution of each of the 27 patients was on average 3.7 ± 3.5% (SD).
To investigate the spatial extent of the observed PAC, we computed the Euclidian distance (using Talairach coordinates) between all electrode pairs involved in the significant wPLFs. We constructed the density of wPLFs as a function of their strength (horizontal axis) and the distance between electrodes within a pair (vertical axis) (Fig. 3B). We observed that (1) PAC occurred predominantly at distances ∼6 cm, (2) PAC occurred at distances as large as 14 cm, and (3) PAC strength decreased with distance.
Next, we investigated the diversity in preferred coupling phases. We constructed the density of significant wPLFs as a function of their phase and their strength (Fig. 3C). We make the following three observations: (1) PAC occurred with diverse preferred coupling phases, (2) for weak coupling, phases were clustered around the peak of the phaseproviding oscillation (phase = 0), and (3) for strong coupling, phases were clustered both around the peak and the trough (phase = ±π). The observed diversity in preferred coupling phase is not produced by sampling error resulting from unreliable phase estimates (Fig. 4).
We obtained peaks of frequency profiles of the amplitude and phaseproviding frequencies of the significant wPLFs for each dataset by counting the significant wPLFs in the three other dimensions. We constructed a scatter plot (Fig. 3D) and observed that PAC involved many frequencies. The peak phaseproviding frequencies showed a substantial spread, ranging from delta (2 Hz or lower) to alpha (12 Hz), and so do the amplitudeproviding frequencies, ranging from alpha (15 Hz) to gamma (67 Hz or higher).
From our analysis of significant wPLFs, we conclude that PAC (1) occurred over long distances, (2) showed substantial diversity in preferred coupling phases, and (3) involved many frequencies. However, analyzing significant wPLFs does not inform us about the spatial distribution of PAC or the origin of the phase diversity. More specifically, we do not know (1) whether the observed longdistance PAC is generated by multiple small spatially separated sources or by one very large source, (2) whether the phase diversity is attributable to phase differences within a source or between sources, or (3) how the spatial distribution of the phaseproviding electrodes is related to that of the amplitudeproviding electrodes. To investigate these issues, we decomposed each 4way array of wPLFs into sets of two spatial maps and two frequency profiles. Importantly, these spatial maps provide information about the spatial distribution of PAC and the origin of the phase diversity.
Nway decomposition reveals the spatial distribution of PAC in sets of two spatial maps and two frequency profiles
To analyze the spatial distribution and phase diversity of PAC, we used Nway decomposition based on the PARAFAC model (see Materials and Methods). This method has been used previously (Maris et al., 2011) but in a version that was unable to reveal the phase diversity that we identified (see Materials and Methods). Each 4way array of wPLFs was decomposed into one or more components. Every component characterizes one PAC pattern and consists of an amplitudeproviding and a phaseproviding spatial map and an amplitudeproviding and a phaseproviding frequency profile. Because a dataset may involve multiple PAC patterns, the decomposition can extract multiple components. To illustrate the decomposition, we show an example component of a representative subject (Fig. 5). The spatial maps are shown as grids on a template brain (Fig. 5A,B; not all ECoG grids are shown), and the frequency content is shown in the frequency profiles (Fig. 5C). With respect to the example component, we observed that (1) both the amplitudeproviding (Fig. 5A) and the phaseproviding (Fig. 5B) spatial map had a wide spatial distribution over cortex, (2) the phaseproviding spatial map had a wider spatial distribution than the amplitudeproviding spatial map, and (3) there was more phase diversity in the phaseproviding than in the amplitudeproviding map (Fig. 5A). These observations are representative for all 42 datasets. Note that the phaseproviding spatial map in the example component shows a spatial structure similar to that of traveling waves. This was the case for 42 of 85 components.
The example component (Fig. 5A–C) reflects the main pattern in the original 4way array of wPLFs. To show this, we selected two electrode pairs that share the same amplitudeproviding electrode (electrode 43) but have different phaseproviding electrodes (electrodes 29 and 57). We show the wPLFs for both pairs, for all frequency pairs (Fig. 5D). The frequencies that exhibit strong coupling closely match the frequency profiles from the decomposition (Fig. 5C). Moreover, the phase difference between the two electrode pairs closely match the phase differences in the phaseproviding spatial map (Fig. 5B). [Note that, in contrast to betweenelectrode phase differences, absolute phases cannot be interpreted (see Materials and Methods).] It is important to note that relative phases in amplitude and phaseproviding spatial maps reflect different properties of spatially distributed PAC (Fig. 2; see Materials and Methods). From the 42 datasets, we extracted 85 components, explaining on average 50.7 ± 20.1% (SD) of the variance of significant wPLFs. We evaluated the accuracy of the reconstruction of all preferred coupling phases. For the phases of all significant wPLFs, we compared their reconstructed and their corresponding observed values. All magnitudes were set to 1. For every dataset, we calculated a coefficient that quantified reconstruction accuracy (see Materials and Methods), which ranges from 0 to 1. This coefficient of reconstruction was on average 0.70 ± 0.15 (SD), indicating that the observed wPLFs can be accurately reconstructed from the decompositions.
We now show two sets of aggregated results obtained from all components. First, we investigated the spatial extent of the spatial maps. We constructed a scatter plot of the mean betweenelectrode Euclidian distance of each spatial map per component (Fig. 6A). We selected electrodes in each spatial map by comparing the componentspecific reconstructed wPLFs to the 99th percentile of the reference distribution used to select significant wPLFs (see Materials and Methods). The mean distance was on average higher for the phaseproviding than for the amplitudeproviding spatial maps (paired samples t test; t_{(84)} = −5.56, p < 1e6). This shows that our PAC was not generated by distributed sharpedged waveforms because both maps would be equally large. This issue has been discussed previously (Maris et al., 2011).
Second, we investigated the frequency profiles of the amplitude and phaseproviding oscillations of all components. We constructed a scatter plot of the peaks of the amplitude and phaseproviding frequency profiles (Fig. 6B). Peaks of the phaseproviding frequency profiles were spread out, ranging from delta (∼2 Hz) to alpha (∼16 Hz). Peaks of amplitudeproviding frequency profiles were spread out even more, ranging from theta (∼5 Hz) to gamma (∼67 Hz or higher). Note the difference with Figure 3D. The peak frequencies determined from the decomposition (Fig. 6B) showed a much larger spread than to those from the 4way array of wPLFs (Fig. 3D). This most likely results from the fact that each 4way array of wPLFs carries several PAC patterns.
Phaseproviding and amplitudeproviding spatial maps show different phase configurations
To investigate the phase configuration within both types of spatial maps, we calculated the phase differences for all possible pairs of electrodes selected from a spatial map (Fig. 7). As before, we selected electrodes based on the componentspecific reconstructed wPLFs (see Materials and Methods). To show the aggregate phase diversity in both spatial maps, we show the density of electrode pairs as a function of their strength and their phase difference (Fig. 7A,B). We observe (1) that both amplitude and phaseproviding spatial maps showed phase differences clustered around 0, (2) that phaseproviding spatial maps also showed phase differences clustered around ±π, and (3) that both maps showed phase differences between 0 and ±π. This indicates that the amplitudeproviding spatial maps mainly showed synchronous phase configurations. Conversely, the phaseproviding spatial maps also showed socalled antiphasic configurations, with two groups of electrodes having small withingroup but large betweengroup phase differences. Importantly, both maps showed phase differences between 0 and ±π. Modulation of phase relations across the complete circle (not limited to synchrony/antiphase) may have important consequences for the computational mechanism that PAC reflects. More specifically, if PAC reflects selective routing of information by modulating excitability of neuronal groups, then phase determines its flexibility.
From this perspective, it is crucial to establish that the observed phase differences between 0 and ±π are not estimation errors but true phase differences. To demonstrate this, we evaluated the splithalf reliability of these phase differences. For each dataset, we randomly partitioned trials in two sets and calculated 4way arrays of wPLFs for both. We then decomposed each array into the same number of components and computed phase differences as above. We show the splithalf reliability as the density of betweenelectrode phase differences from the first set of trials plotted against those of the second set (Fig. 7C,D). We observe that all phase differences were highly similar between the two sets of trials, and this holds for both spatial maps. Thus, phase differences in the decomposition reflect true phase differences. As an additional quantification, we calculated differences between the two sets of phase differences (Fig. 7E,F). We selected phase differences between −2π/3 to −π/3 and π/3 to 2π/3, as estimated using the first set of trials. For both types of spatial maps, betweenset phase differences cluster around 0 (Fig. 7E,F). This indicates that, even for this narrow range, the phase differences reflect true phase differences. To support this quantitatively, we computed a reliability coefficient, ranging from 0 to 1 (see Materials and Methods). For both spatial maps, the coefficient r was high for both intervals (Fig. 7E,F).
Different connectivity structure for components with different phaseproviding frequencies
To reveal the connectivity structure of PAC in our 85 extracted components, we aggregated over all pairs of selected electrodes (Fig. 8). As before, we selected electrodes based on the componentspecific reconstructed wPLFs (see Materials and Methods). To visualize the connectivity structure, we downsampled the anatomical locations of electrodes (using Talairach coordinates) to 23 locations on the left and 23 locations on the right hemisphere. We did not observe all possible connections based on these downsampled locations: of the 2116 possible connections, we observed 1698 (80.3%) in our data. However, as will become clear, our quantifications and subsequent comparisons were not biased by this incompleteness. For every observed connection, we calculated the proportion selected electrode pairs within a downsampled connection. This proportion estimates the probability that PAC is found between the downsampled locations. The number of electrode pairs over which this proportion was calculated differed greatly across location pairs. This was especially the case for the small number of contralateral versus the large number of ipsilateral electrode pairs (Fig. 8C).
We use connectograms to show the connectivity structure of PAC (Fig. 8A). Whereas no clear structure is revealed by the connectogram for all components together (Fig. 8A, top connectogram), much structure is revealed when separate connectograms are made for components with different phaseproviding frequencies (Fig. 8A, three bottom connectograms). In these connectograms, line color indicates the amplitudeproviding frequency of each connection; pie charts show the percentage of nonzero contralateral and ipsilateral connections. Because of spatial downsampling, every connection consists of many electrode pairs from multiple components. The frequencies used for the connectograms were obtained by taking the mode over the peaks of all their frequency profiles. Components that differ in their phaseproviding frequency differ greatly in their connectivity pattern, including their degree of crosshemispheric lateral connections. We observe that (1) delta oscillations show mostly ipsilateral connections, (2) theta oscillations show both ipsilateral and contralateral connections, and (3) alpha oscillations show predominantly contralateral connections. The imbalance of contralateral versus ipsilateral connections in the top connectogram is partly the result of the relatively small number of patients with bilateral recordings sites (11 of 27). This source of the imbalance also affects the bottom three connectograms, whose combination forms the top connectogram. Importantly, because this source of imbalance affects all bottom three connectograms to the same extent, their number of contralateral connections can be safely compared with each other.
Our findings were obtained by investigating the number of connections that exist between downsampled locations. We also investigated the coupling strength between these locations, calculated as the median connection strength over all electrode pairs belonging to that connection. For each electrode pair, connection strength was calculated by selecting their componentspecific reconstructed wPLFs and taking a weighted average over the two frequency dimensions. We now report on this, with separate averaging over (1) contralateral and ipsilateral and (2) within and betweenlobe connections (Fig. 8B). We observe that (1) there is only a small difference in connection strength between lobes compared with within lobes, and (2) the contralateral and ipsilateral connection strengths did not differ greatly: the average strength of the contralateral connections was 69.7% of the average strength of the ipsilateral connections (Fig. 8B). We also investigated the above patterns separately for each of the three phaseproviding frequencies (delta, theta, and alpha), but we found no substantial differences (figure not shown).
Discussion
We provided evidence for three key properties that could allow PAC to flexibly and selectively route information through distributed neuronal networks. (1) We showed that PAC was widely spatially distributed. From our analyses of significant wPLFs, we found that PAC occurred at distances >10 cm. Using a decomposition based on the PARAFAC model, we showed that this PAC was generated by spatially distributed phase and amplitudeproviding oscillations, of which the phaseproviding oscillations were more spread out. The spatial distribution of PAC is required to be able to route information through distributed networks. (2) We showed that, over these spatially distributed networks, there was great phase diversity. The phase diversity we observed was mainly explained by phase diversity in the phaseproviding oscillation, showing phase differences over space, across the whole circle. This phase diversity can determine the flexibility of PAC in selecting neuronal populations for interaction. (3) We showed that PAC occurred between oscillations of many different frequencies (He et al., 2010; Miller et al., 2010; Maris et al., 2011). Amplitudeproviding frequencies ranged from theta to gamma, and phaseproviding frequencies ranged from delta to alpha. This frequency diversity can determine the flexibility of PAC in separating neuronal networks operating in parallel.
Besides providing evidence for these three key properties, we also made a first step toward identifying different roles for phaseproviding oscillations at different frequencies. We showed that delta oscillations establish mostly ipsilateral connections, theta oscillations establish both ipsilateral and contralateral connections, and alpha oscillations establish predominantly contralateral connections.
Although we find PAC that is widely spatially distributed, most reports so far have shown local PAC (Chrobak and Buzsáki, 1998; Lakatos et al., 2005; Mormann et al., 2005; Canolty et al., 2006; Cohen, 2008; Lakatos et al., 2008; Osipova et al., 2008; Penny et al., 2008; Axmacher et al., 2010; He et al., 2010; Miller et al., 2010; Voytek et al., 2010). The evidence for crossarea PAC is much less abundant (Sirota et al., 2008; Tort et al., 2008; Maris et al., 2011). This is surprising given that, at least for theta, oscillations are measurable over broad regions (von Stein and Sarnthein, 2000; He et al., 2008).
PAC is a form of oscillatory synchronization that involves two different frequencies. A more familiar form of oscillatory synchronization involves only a single frequency: phase consistency between oscillations. This type of oscillatory synchronization was originally put forward as a mechanism to bind different features of an object, encoded in different neuronal populations, and was confined to gamma (Singer and Gray, 1995; Singer, 1999). More recently, it has been proposed as a mechanism of communication between neuronal groups (Fries, 2005). Central to this proposal is that neurons preferably fire in a specific phase of the gamma cycle, implying a temporal relation between spikes of coherently oscillating neurons. This allows neurons to synchronize their periods of maximum excitability and communicate effectively. An expanded mechanism, involving PAC, would be that the phaseproviding oscillation modulates when neuronal populations engage in oscillatory phase synchronization. Such a mechanism would require phase consistency during highfrequency bursts. Although we did not investigate this, longrange phase consistency at high frequencies has been reported (Gregoriou et al., 2009). If it is indeed the case that phaseproviding oscillations select neuronal populations for communication, then such crossfrequency interaction could provide dynamic gating of information (Vogels and Abbott, 2009; Akam and Kullmann, 2010). The phaseproviding oscillation as a selector would also greatly benefit from a process that creates substantial phase diversity across sites, which we have observed, in which sitespecific phases could function as a selection variable. Besides the separation of neuronal activity via phase diversity (i.e., “phase multiplexing”), the separation of multiple networks operating in parallel could also be supported by frequency diversity (i.e., “frequency multiplexing”).
Although we provide evidence for three key properties that makes PAC a likely candidate for routing for information, we do not have a mechanistic neurophysiological model to explain our observations, and providing such an explanation is a crucial challenge for future research. Obviously, neuronal spiking is the signal for targeted communication between neurons. Therefore, we face the challenge to link our findings using ECoG to neuronal signals with a very different spatial specificity. Inevitably, any attempt to provide such a link will involve some speculation. However, it is justified to the extent that it can be related to relevant findings in the literature. We hypothesize that PAC may reflect an interaction between slow and fast rhythmic synaptic input streams and that the efficacy of this interaction can be modulated by adjusting the timing between these rhythms. Furthermore, we hypothesize that the slow rhythm, the phaseproviding oscillation, can be used to segregate information streams by variations in its phase over space. These hypotheses have links with the substantial literature on up and down states (Steriade et al., 1993; Destexhe et al., 2003; Cash et al., 2009; Haider and McCormick, 2009).
Up and down states are brain states that can be characterized at many levels, ranging from intracellular recordings to macroscopic electroencephalography (Steriade et al., 2001; Destexhe et al., 2003; Volgushev et al., 2006). Up and down states are most easily identified during slowwave sleep and anesthesia, when they alternate rhythmically (Hô and Destexhe, 2000; SanchezVives and McCormick, 2000; Shu et al., 2003; Hasenstaub et al., 2005; Haider et al., 2006). This typically occurs at a frequency <1 Hz, producing the socalled cortical slow oscillation. During wakefulness and shallow sleep, isolated down states called Kcomplexes occur in a nonrhythmic manner, sometimes preceded by an up state (Amzica and Steriade, 1997; Cash et al., 2009). In both rhythmic and nonrhythmic cases, up and down states are spatially distributed (Volgushev et al., 2006). Furthermore, recordings in animals (SanchezVives and McCormick, 2000; Volgushev et al., 2006) and humans (Massimini et al., 2004) have shown that the cortical slow oscillation is phase diverse over space. In ECoG recordings, up and down states can be identified as positive and negative deflections (Cash et al., 2009). These originate from a sourcesink configuration, with a source in layer II/III and a sink in layer I (Cash et al., 2009). The neocortical up state is a network phenomenon, characterized by a balance in excitatory and inhibitory input (Haider et al., 2006). During an up state, neurons have an increased membrane potential, bringing them closer to their firing threshold (Haider et al., 2007). Animal studies have shown that, during an up state, there is an increase in spiking activity and highfrequency (30–100 Hz) local field potential (LFP) fluctuations (Nowak et al., 1997; Haider and McCormick, 2009). A down state provides the opposite: the membrane potential is hyperpolarized, and there is a strong decrease in spiking and highfrequency LFP fluctuations.
The rhythmic alternation of up and down states has similarities to our PAC results: (1) it involves a slow oscillation with highfrequency oscillations occurring at specific phases, (2) the slow oscillation is spatially distributed, and (3) the slow oscillation is phase diverse over space. Based on these similarities, we propose that phaseproviding oscillations of PAC could affect neuronal populations in a similar way. The phaseproviding oscillation could, like the cortical slow oscillation, modulate the basal membrane potentials of neuronal groups, which could provide joint windows of communication during certain phases. This could allow neuronal groups to exchange information during bursts of highfrequency LFP fluctuations. This information exchange could involve coherent oscillations. Importantly, there are a number of differences between up and down states, and our PAC results complicate this comparison. First, slow cortical oscillations are typically <1 Hz, whereas we identified phaseproviding oscillations between 2 and 16 Hz (lower boundary is restricted by epoch length). Second, the phase diversity we observed across space is much larger than the phase diversity in the slow oscillation. Last, highfrequency LFP fluctuations during up states are locked to the peak of the slow oscillation, whereas we report strong diversity in preferred coupling phases of PAC. The reported phase of the slow oscillation, however, depends on where in the neuropil the signal is recorded: because the up and down states are characterized by a sourcesink configuration between layer I and layer II/III (Cash et al., 2009), the polarity of the slow oscillation would reverse if one would record from layer II/III. Assuming that the neurophysiological makeup of the brain allows for different sourcesink configurations across layers, it should be possible to generate PAC with diverse phases as measured on the surface of the brain.
More evidence is needed to show that spatially distributed PAC in ECoG signals reflect a neurophysiological mechanism that also modulates spiking activity. At least part of the evidence must come from in vivo experiments involving simultaneous recordings from multiple sites in the neuropil. The present study can guide the analysis of these recordings. We have demonstrated that PAC is a pervasive phenomenon that has a wide spatial distribution, a strong diversity in preferred coupling phases, and involves oscillations at many frequencies. With these properties, PAC is a plausible candidate for supporting selective neuronal communication.
Footnotes

We gratefully acknowledge the support of the BrainGain Smart Mix Programme of the Dutch Ministry of Economic Affairs and the Dutch Ministry of Education, Culture, and Science.

The authors declare no competing financial interests.
 Correspondence should be addressed to Eric Maris, Radboud University Nijmegen, Donders Institute for Brain Cognition and Behaviour, 6525 HR Nijmegen, The Netherlands. e.maris{at}donders.ru.nl