# Detecting event-related changes of multivariate phase coupling in dynamic brain networks

## Abstract

Oscillatory phase coupling within large-scale brain networks is a topic of increasing interest within systems, cognitive, and theoretical neuroscience. Evidence shows that brain rhythms play a role in controlling neuronal excitability and response modulation (Haider B, McCormick D. *Neuron* 62: 171–189, 2009) and regulate the efficacy of communication between cortical regions (Fries P. *Trends Cogn Sci* 9: 474–480, 2005) and distinct spatiotemporal scales (Canolty RT, Knight RT. *Trends Cogn Sci* 14: 506–515, 2010). In this view, anatomically connected brain areas form the scaffolding upon which neuronal oscillations rapidly create and dissolve transient functional networks (Lakatos P, Karmos G, Mehta A, Ulbert I, Schroeder C. *Science* 320: 110–113, 2008). Importantly, testing these hypotheses requires methods designed to accurately reflect dynamic changes in multivariate phase coupling within brain networks. Unfortunately, phase coupling between neurophysiological signals is commonly investigated using suboptimal techniques. Here we describe how a recently developed probabilistic model, phase coupling estimation (PCE; Cadieu C, Koepsell K *Neural Comput* 44: 3107–3126, 2010), can be used to investigate changes in multivariate phase coupling, and we detail the advantages of this model over the commonly employed phase-locking value (PLV; Lachaux JP, Rodriguez E, Martinerie J, Varela F. *Human Brain Map* 8: 194–208, 1999). We show that the N-dimensional PCE is a natural generalization of the inherently bivariate PLV. Using simulations, we show that PCE accurately captures both direct and indirect (network mediated) coupling between network elements in situations where PLV produces erroneous results. We present empirical results on recordings from humans and nonhuman primates and show that the PCE-estimated coupling values are different from those using the bivariate PLV. Critically on these empirical recordings, PCE output tends to be sparser than the PLVs, indicating fewer significant interactions and perhaps a more parsimonious description of the data. Finally, the physical interpretation of PCE parameters is straightforward: the PCE parameters correspond to interaction terms in a network of coupled oscillators. Forward modeling of a network of coupled oscillators with parameters estimated by PCE generates synthetic data with statistical characteristics identical to empirical signals. Given these advantages over the PLV, PCE is a useful tool for investigating multivariate phase coupling in distributed brain networks.

- neuronal oscillations
- phase-locking value
- functional networks
- multivariate analysis

the brain exhibits network structure across a variety of different spatial scales, from the local synaptic connectivity of nerve cells within a cortical column (White and Keller 1989) to the large-scale networks of interconnected brain areas responsible for language, attention, and memory (Sporns 2011). Furthermore, goal-directed complex behavior such as reading a sentence requires the dynamic coordination of multiple distinct subnetworks, with computation and communication occurring via the transient activation and deactivation of appropriate functional networks. However, structural or anatomical connectivity is relatively static at the time-scale of perception and action. What mechanism could provide the fast dynamic regulation of active brain networks required to survive in an unpredictable world?

One possibility is that a hierarchy of neuronal oscillations could provide the dynamic scaffolding required to modulate network activity at a variety of different scales (Lakatos et al. 2005). In particular, the communication through coherence hypothesis (Fries 2005) proposes a mechanistic solution for dynamically regulating the strength of connectivity between brain areas via oscillatory phase coupling. In this view, oscillations within a given brain area are rhythmic variations in neuronal excitability that influence the effective gain of spiking input to an area as well as the probability of spike output from an area; that is, neuronal oscillations define recurring windows of communication for a cortical region. Therefore, if two interconnected brain areas are oscillating at a given frequency, then the relative phase difference between them may influence their effective connectivity. For example, given short axonal conduction delays, if two areas have a relative phase lag of 0 radians, then their open windows of communication will occur simultaneously, increasing effective connectivity, whereas a relative phase lag of π radians misaligns the windows of communication in the two areas, such that spikes leaving one area will arrive at the least effective time to elicit spiking in the other area. Finally, dynamic regulation of network activity can be accomplished via fast changes of the relative phase offset between areas. Proper empirical evaluation of the communication through coherence hypothesis thus requires a technique for tracking transient phase coupling in brain activity.

### Existing Methods for Estimating Phase Coupling

#### Linear coherence and its limitations.

A number of techniques exist for estimating phase coupling from brain activity. Such activity is typically recorded from microelectrodes as local field potentials (LFPs), from subdural macroelectrodes as the electrocorticogram (ECoG), from macroelectrodes on the scalp as the electroencephalogram, or from measurement of magnetic activity outside the head as the magnetoencephalogram. For each of these types of signals, the activity in a given frequency band can be described in terms of the amplitude and phase of filtered signals. One popular technique for investigating frequency-specific dependence between two signals is coherence (Carter 1987; Nunez et al. 1997). Given two wide-sense stationary signals *x* and *y*, we can compute the cross power spectrum
*R*_{xy}(τ) = *E*[*x*(*t*)*y*(*t* − τ)] is the cross-correlation function and *E*[.] is the expectation operator. The autospectra for *x* and *y* can be computed similarly. From these, we can calculate the complex-valued coherence as
*x* and *y* before computing the MSC allows one to track changes in coherence over time. However, while the MSC has seen wide use in neuroscience, two aspects of the coherence method make it less than ideal for evaluating event-related phase coupling in neuronal signals.

First, linear coherence is dependent on amplitudes as well as phases, that is, the coherence method is sensitive to correlated amplitude fluctuations and may exhibit significant event-related changes even if the phases of the two signals are statistically independent (Nunez and Srinivasan 2006). While the instantaneous value of amplitude or power in an area is also likely to influence effective connectivity between areas, it would be preferable to employ methods that can distinguish the influence of amplitude and phase separately, such as phase-only coherence (Nunez and Srinivasan 2006).

Second, event-related coherence methods tend to employ a sliding-window approach that estimates coherence over short time intervals. In this approach, the same window length (temporal duration) is employed for all frequencies from very low to very high frequencies. However, the duration of neuronal oscillations does not appear to be the same for all frequencies but is rather a function of the mean oscillatory frequency (Durka 2007). That is, the duration of an oscillatory “wavelet” or burst of rhythmic activity will be shorter for high frequencies than the duration of oscillatory activity occurring at low frequencies. Similar to techniques that employ the natural statistics of signals for compact, high signal-to-noise decompositions (Smith 2006), it would be preferable to investigate phase coupling using a method that employs optimal window sizes for the different frequencies of interest.

### Phase-Locking Value and Its Limitations

One popular technique that addresses both of these issues is the phase-locking value (PLV) method (Lachaux et al. 1999). First, this technique investigates one frequency at a time, permitting investigators to employ whatever filtering/windowing parameters are most appropriate for the question of interest. Second, the PLV method normalizes the filtered (complex valued) signals to have unit amplitude, so that only instantaneous phases are involved in the computations. Specifically, one formulation of the PLV method employs three steps: first, given two raw (real-valued) time series *x*_{1} and *x*_{2}, the instantaneous phases of both signals are extracted to generate two (complex-valued) time series *z*_{1} and *z*_{2}. These new signals *z*_{1} and *z*_{2} are normalized such that the modulus (amplitude) of each sample point *z*_{1}(*t*) = *z*_{2}(*t*) = 1 for all sample points *t*. This phase extraction step can be accomplished via convolution of the raw signal with an appropriate basis function [time-frequency signal atom; see Grochenig (2001) for more details]. Second, for event-related data, epochs or windows around the trial events of interest are extracted from each signal and the modulus (amplitude) of the trial-averaged phase differences is computed:
*z* = exp[iθ] ∈ C, |.| indicates complex modulus, |*z*| = 1, *i* = *k* indexes the trial number, and *m* and *n* index the first and second channels of interest, respectively. This trial-averaged modulus is the PLV and takes on values between 0 and 1. A PLV of 1 can only occur if the relative phase difference is identical from trial to trial. If the relative phase difference between signals is independent from one trial to the next and uniformly distributed over all possible phase values, then the PLV will be close to 0. Finally, the statistical significance of the PLV actually observed can be estimated through a randomized permutation test. That is, for some large number of surrogate trials (say, *N*_{SURROGATE} = 1000 surrogate runs), one computes a set of surrogate PLV {*PLV*_{SUR}}. For each surrogate run, the ordered list of trial epochs is preserved for signal 1, but is randomly permuted for signal 2. Rather than pairing trial *k* of *z*_{1} to trial *k* of *z*_{2}, now trial *k* of *z*_{1} is paired with trial *h* of *z*_{2}, where *h* is sampled without replacement from the set of integers {1, 2, …,*N*_{TRIALS}}. The trial-shuffling procedure preserves event-related phase changes caused by trial events alone and discounts any nontrial specific coupling between the two signals (which is the effect of interest). Statistical significance can be assessed by comparing the actual PLV to the set of surrogate PLV {*PLV*_{SUR}}; if an event-related increase in PLV is expected and only *s* surrogate PLV are larger than the actual PLV, this corresponds to *P* value of

The PLV method is an excellent technique for estimating the strength and significance of instantaneous phase coupling between two brain signals, and it has been used to investigate a multitude of different functional brain networks. Nonetheless, the PLV method suffers from two fundamental limitations that motivate the development of alternative methods. First, the output of the PLV is a single scalar number or index rather than a probability distribution, and second, the PLV is an inherently bivariate measure.

Because the PLV is an index and not a distribution, it cannot be directly used to generate simulated data with identical phase coupling characteristics. However, the set of single-trial relative phase differences used to calculate the PLV can also be used to fit a parametric statistical model. One of the simplest and most tractable statistical models over circular or phase variables is the von Mises distribution:
*I*_{0}(.) is the modified Bessel function of order 0. Once the von Mises parameters μ and κ have been fit for a set of relative phase differences (where ϕ = θ_{m} − θ_{n}) it is straightforward to generate synthetic data with the same phase coupling properties as the empirically recorded data; that is, synthetic data generated in this way will yield the same pair-wise PLVs as the empirical data. The generation of synthetic data is beyond the scope of this study [see Strogatz (1994) and Press et al. (2007) for additional details and standard numerical techniques]. In our experience working with the subdural electrocorticogram in humans and macaque LFPs, the von Mises distribution often provides an excellent model for the empirical distribution of phase differences.

The second limitation of the PLV is that it is an inherently bivariate measure. That is, if *N* different LFP channels are simultaneously recorded during an experiment, the PLV method requires that all possible pairs of channels be examined separately, which can produce misleading results. For example, suppose we have a three-node network of coupled oscillators including nodes *A*, *B*, and *C*. If nodes *A* and *B* are both phase coupled to node *C*, then phase value θ_{A} from node *A* and phase value θ_{B} from node *B* may nonetheless appear dependent even if there is no direct influence between nodes *A* and *B*. That is, it may be the case that
*p*(*x|y*) denotes the probability of observing an event *x*, given that *y* is known. That is, *Eq. 6* says that the phases of oscillators *A* and *B* are statistically dependent (if you discover the value of the phase of *A*, then you have also reduced the uncertainty of your knowledge of the phase of *B*), but at the same time *Eq. 7* says that the phases of *A* and *B* are conditionally independent, given the phase of oscillator *C* (that is, if you discover the phases of both *A* and *C*, you still have the same amount of uncertainty as before about the phase of *B*). Therefore, examining nodes *A* and *B* in isolation and ignoring the influence of node *C* may produce incorrect estimates of the network structure. The ideal case, given *N-*coupled oscillators, would be to fit a probability density function (pdf) describing the likelihood of seeing every possible *N*-dimensional vector of phases. In contrast to the PLV, such a multivariate pdf would permit one to draw samples from it to run numerical simulations as well as avoiding the confounding issues of conditional dependence inherent in applying bivariate measures to a highly-interconnected network of coupled oscillators.

## METHODS

Here we describe a method for estimating task- and event-related changes in phase coupling that improves upon these existing methods. In particular, this method termed phase coupling estimation (PCE) exhibits three key advantages over the widely used PLV method: *1*) it is a fully multivariate approach that can distinguish direct coupling from indirect coupling between nodes in a large network; *2*) it is a probabilistic method that permits resampling from the best-fit distributions, aiding in-depth simulation studies; and *3*) it is a dynamical systems approach where the model parameters have a clear physical interpretation as the coupling terms describing interactions within a network of coupled oscillators (Cadieu and Koepsell 2010).

### PCE

The probability model used in multivariate phase coupling estimation is the maximum entropy distribution for *N*-dimensional phase vectors given pair-wise statistics. A number of maximum entropy distributions are used throughout the science and engineering communities. In the real-valued case the multivariate Gaussian distribution, and in the binary case the Ising model, serve as widely used multivariate maximum entropy distributions consistent with second order statistics. For a one-dimensional phase variable, the maximum entropy distribution given the first circular moment is the von Mises distribution (Kotz et al. 1985). Through a change of variables
*m* and *n* using the von Mises pdf:
_{m}, θ_{n}, ϕ_{mn}, and μ_{mn} ∈ [−π, π) and κ_{mn} ∈ (0, +∞). The phase vector θ = (θ_{1}, θ_{2}, …, θ_{N}) from *N* simultaneously recorded channels is a point on the *N*-torus. We can generalize the von Mises pdf to *N* dimensions, thus capturing the full multivariate dependence structure. Given the full set of pair-wise phase statistics, there is a unique maximum entropy distribution that reproduces the first and second statistics of these measurements (Cadieu and Koepsell 2010). Here we refer to this multivariate distribution of phase vectors, and its efficient estimation via score-matching (Cadieu and Koepsell 2010), as PCE.

For multivariate phases, the first circular moment is a measurement between two phases, θ_{m} and θ_{n}, and is defined as the complex quantity: 〈exp[*i*(θ* _{m}* − θ

*)]〉, where 〈.〉 is the average over samples. The real and imaginary parts are given as:*

_{n}Given these statistics, it can be shown (Cadieu and Koepsell 2010) that the corresponding maximum entropy distribution is given as:
_{1}, θ_{2}, …, θ_{N}) is the *N* dimensional set of phases and **K** specifies the parameters of the distribution. We have used trigonometric identities to combine the sine and the cosine of the differences of the phase pairs into one term for each pair of phases. The terms κ_{mn} and μ_{mn} are the coupling between phases θ_{m} and θ_{n} and the phase offset between phases θ_{m} and θ_{n}, respectively. The term *Z*(**K**) is the normalization constant and is dependent on the parameters of the distribution. Given phases from *N* different LFP channels, we can estimate the probability of observing a particular *N*-dimensional vector of phases using a multivariate phase distribution. An equivalent but more compact expression for the probability distribution given in *Eq. 12* is,
*N*-dimensional vector of phase variables *z* as a vector of unit length complex variables *z*_{m} = exp[*i*θ_{m}) and θ_{m} is an element of the real-valued interval [−π, π]. The *N* × *N* coupling matrix **K** is Hermitian and traceless. The elements of K encode the coupling parameters between channels; e.g., **K**_{mn} encodes the coupling between the *m*-th and *n*-th phase variables. Each element of K is a complex number **K**_{mn} = κ_{mn} exp[*i*μ_{mn}], where the modulus κ_{mn} encodes coupling strength and the angle μ_{mn} denotes the preferred phase offset between channels. The diagonal elements of **K** are zero (**K**_{ii} = 0), but nonuniform univariate phase distributions can be modeled by augmenting the observed matrix of phase variables with an additional variable of fixed phase, resulting in a (*N* + 1) × (*N* + 1) coupling matrix. The normalization constant *Z*(**K**) is a function of the coupling matrix and in general cannot be computed analytically. Note that *Eqs. 12* and *13* are equivalent, but *Eq. 13* uses complex notation. Given an observed set of phase variables, we then estimate the parameters of the distribution using an efficient technique derived in (Cadieu and Koepsell 2010). The lack of a closed form to the partition function *Z*(**K**) makes standard maximum-likelihood estimators computationally expensive and prone to convergence problems. The estimator derived in (Cadieu and Koepsell 2010) is a linear system of equations using the measurements of the phase variables. This linear system of equations can be solved using standard techniques; code to estimate the distribution using score-matching (Hyvarinen 2005) is available at https://github.com/koepsell/phase-coupling-estimation. This estimator has been shown to correspond to the maximum-likelihood estimate and performs well in high dimensions and with limited data (Cadieu and Koepsell 2010). As with other statistical measures, we will be less certain of PCE parameter values when we have few sample points per dimension for high-dimensional data sets, but as shown in Fig. 3 of Cadieu and Koepsell (2010), accurate parameter estimates can be obtained for as few as 100 sample points per dimension.

Interestingly, the estimated parameters of this maximum entropy multivariate phase distribution have a physical interpretation in a dynamic system of coupled oscillators. In fact, the parameters in the phase distribution are identical to the interactions between the oscillators. This implies that rather than starting with a set of pair-wise phase statistics and deriving the maximum entropy distribution consistent with these statistics, we can instead derive the multivariate phase distribution from a dynamical systems model of coupled oscillators. Given the dynamical system,
*Eqs. 12* and *13*, save for the introduction of a multiplicative factor β within the exponential to account for the variance of the noise terms *v*_{m}(*t*). Thus the parameters of the matrix **K** estimated from observed phase data may be interpreted as the interaction terms between a physical system of coupled oscillators. For simplicity, in the remainder of this work, we ignore this parameter β [see Cadieu and Koepsell (2010) for more details].

In summary, the multivariate phase distribution of *Eq. 13* provides the most parsimonious statistical model of the joint multivariate phase distribution given only pair-wise phase measurements. Maximum entropy solutions make the fewest assumptions required to satisfy a given set of constraints (here, matching the first and second moments of the data) and are preferred to more structured statistical models when the true joint distribution is unknown. PCE is a technique that can be used to estimate the parameters of the joint multivariate phase distribution from measurements. The coupling parameters in the estimated multivariate phase distribution can be interpreted as interactions in a physical system of coupled oscillators. For further details on this distribution, the PCE technique, and the interpretation of this distribution as a physical system of oscillators, we refer the reader to Cadieu and Koepsell (2010).

### Relationship Among the PLV, the Von Mises Distribution, and PCE

How is the bivariate PLV related to the multivariate PCE? While the PLV and von Mises distribution are closely related to each other, they are only indirectly related to the PCE coupling parameters. This is important because the PCE coupling parameters can be used to distinguish direct from indirect coupling between pairs of oscillators. That is, in the context of brain networks, two cortical areas *A* and *B* can influence each other directly, presumably through monosynaptic anatomical corticocortical projections, or they can influence each other indirectly via indirect polysynaptic chains that involve one or more intermediate cortical areas. In particular, the oscillatory signals recorded from areas *A* and *B* may exhibit statistical dependence even if there is no direct link between *A* and *B*. This is a severe limitation of the *PLV*, which cannot distinguish direct from indirect coupling. In this section, we show exactly how the *PLV* relates to PCE, and how PCE can distinguish direct from indirect coupling.

Recall that the PLV is the amplitude of the first circular moment of the measured phase difference between two phases:
*m* and *n*, we can see the relationship between the phase-locking value *PLV*_{mn} and the PCE coupling parameter **K**_{mn} by examining the marginal distribution of phase differences. The marginal distribution of phase differences between channels *m* and *n* indicates how strongly coupled the two channels are, including all possible network paths between *m* and *n* (that is, including both direct and indirect or intermediate connections). Strong coupling will show up in the marginal distribution as a narrow distribution over phases, while weak coupling will be indicated by the presence of a broad distribution over phases. In other words, given the phase of one channel (say, channel *m*), the marginal distribution of phase differences between *m* and *n* tells us how much uncertainty remains about the phase value of channel *n*. The marginal distribution is defined as:
_{k} with *k* ≠ *m*, *n*, which can be either the first or second variable in the cosine. Intuitively, the effect of the integration in *Eq. 16* is to condense the influence of all possible network paths between *m* and *n* into a single “virtual” link between channels *m* and *n*, now represented by the one-dimensional probability distribution. After applying the variable substitution θ_{k} = ∼θ_{k} + θ_{n}, all terms in *Eq. 16* either depend on the phase difference θ_{m} − θ_{n}, or are independent of θ_{m} and θ_{n}. The independent terms integrate to a constant and the remaining terms combine to a distribution in the pair-wise phase difference given by
_{mn} and concentration parameter γ_{mn}. *Equation 17* states that the (unnormalized) distribution on the left-hand-side is proportional to the expression on the right-hand-side, but one of the axioms of probability theory is that the integral of a continuous probability distribution must be equal to unity (Jaynes 2003). Solving for the normalization term that makes this true, we find that
_{mn} for a pair of phases the *phase concentration*. Recall that *Eq. 4* relates the PLV to the *first circular moment* such that 〈exp[*i*(θ* _{m}* − θ

*)]〉 =*

_{n}*PLV*exp[

_{mn}*i*Δ

*]. This indicates how we can compute the parameters of the distribution in*

_{mn}*Eq. 18*: the mean phase Δ

_{mn}is the angle of the first moment and the concentration parameter γ

_{mn}can be obtained by numerically solving the equation

*I*

_{0}(

*x*) and

*I*

_{1}(

*x*) denote the modified Bessel functions of zeroeth and first order, respectively (Evans et al. 2000).

Therefore, there is a nontrivial relationship between the PLV and the PCE coupling parameters. There is a simple one-to-one mapping between the phase locking value *PLV*_{mn} and the von Mises concentration parameter γ_{mn} through *Eq. 19*, but the value of the (bivariate) von Mises concentration parameter γ_{mn} is related to the (multivariate) coupling parameters **K** through the nonintuitive *Eq. 16* together with *Eqs. 17* and *18*. Thus the *PLV*_{mn} is related to the full set of PCE coupling parameters **K** through *Eqs. 16*, *17*, *18*, and *19* (see also Fig. 1 for a schematic overview of these different measures).

One intuitive way to make sense of these nontrivial relationships is to recall the dynamical system interpretation of PCE. Under this interpretation of the probability distribution, the interaction between two oscillators *m* and *n* is given by the complexvalued coupling parameter **K**_{mn} = κ_{mn} exp[*i*μ_{mn}], where κ_{mn} represents the phase concentration due to *direct* coupling alone, and μ_{mn} is the mean phase offset due to direct coupling alone (e.g., *m* ↔ *n*). However, in a complex network there will also exist many *indirect* paths linking *m* and *n* (e.g., *m* ↔ *k* ↔ *n*, *m* ↔ *k* ↔ *l* ↔ *n*, etc.), and these indirect paths will also influence the value of PLV and the marginal (von Mises) distributions. In general there is no simple relationship between the PCE coupling parameters and the measured PLV or phase concentration. However, unlike bivariate approaches that aggregate direct and indirect interactions into one lumped *total* interaction term, the PCE approach permits investigators to infer the direct interactions between the oscillators separately from the indirect interactions mediated through the rest of the network. This separation is essential for accurately capturing dynamic activity in complex networks. We make explicit these direct and indirect interactions in the next section.

### Empirical, Isolated, and Network Distributions

In terms of marginal phase distributions, we can separate the influence of direct and indirect modes of coupling. That is, given a set of phase measurements we can directly compute the marginal distribution of the phase difference between a specific pair of phases. We call the marginal distribution computed from the difference of phase measurements of θ_{m} and θ_{n} the *empirical* distribution *p*(θ_{m} − θ_{n}). In a network of many oscillators, the empirical distribution is determined by a direct interaction between nodes *m* and *n* and an indirect interaction through the rest of the network. However, this empirical distribution can be decomposed into an *isolated* distribution, which captures the direct interaction, and a *network* distribution, which captures the interaction through the network.

For a given set of oscillators and coupling parameters the empirical distribution is given as:
*Eq. 16* but with the product containing only one term for each pair of oscillators. The integration is over all phases θ_{k} with *k* ≠ *m, n*. Because the integration is over phases from all channels except *m* and *n*, we can factor out the terms containing the coupling parameters between nodes *m* and *n*:
_{k} = ∼θ_{k} + θ_{n}, and all terms in *Eq. 21* either depend on the phase difference θ_{m} − θ_{n} or are independent of θ_{m} and θ_{n}. The independent terms integrate to a constant and the remaining terms combine to a von Mises distribution in the pair-wise phase difference. Therefore, the empirical distribution *p*(θ_{m} − θ_{n} **K**), written as *p*_{emp}(θ_{m} − θ_{n} **K**) for clarity, is proportional to a product of two von Mises distributions:
**K̄**_{mn} is the set of parameters excluding the direct coupling parameters κ_{mn} and μ_{mn}.

The (indirect) concentration κ̄_{mn} and (indirect) phase offset μ̄_{mn} are determined through the integral in *Eq. 21* and depend on all of the terms in the coupling matrix **K** excluding the direct coupling term **K**_{mn} = κ_{mn} exp[*i*μ_{mn}]. We refer to the distribution that contains the direct coupling parameters the isolated distribution *p*_{iso}(θ_{m} − θ_{n} κ_{mn}, μ_{mn}) because it is the distribution that would be measured if these two oscillators were isolated from the rest of the network and their interaction would only be determined by the direct interaction. We refer to the distribution that contains the network effects on the empirical distribution the network distribution *p*_{net}(θ_{m} − θ_{n} **K̄**_{mn}) because it is the distribution that would be measured if there were no direct interaction between the nodes and only indirect interactions through the network were present.

## RESULTS

### PCE Validation Via Simulations of Coupled Oscillator Networks

The ability to separate direct from indirect interactions has important consequences for accurate estimation of network structure and explains the differences observed between PLV and the coupling parameters estimated from the PCE approach. Here we illustrate these differences using simulations of simple networks of coupled oscillators, where “ground truth” is known.

For the first simple network, we consider three nodes and are interested in determining the interaction between nodes *A* and *B* (see Fig. 2*A*). A third node (unlabeled) is coupled to both nodes *A* and *B* (κ = 1.1 and μ = 0 to node *A*, and κ = 0.9 and μ = 0 to node *B*), but nodes *A* and *B* have no direct coupling (κ = 0). We simulate this network using a discrete time simulation of equation *Eq. 14*, a sampling rate of 1,000 Hz, an oscillator frequency, ω, of 10 Hz, and sample *v* to be independent Gaussian noise with zero mean and a variance of 0.002. This simulation generates a phase time-series for each oscillator in the network. We then used the phase time series from oscillators *A* and *B* (ignoring the third unlabeled oscillator) to compute the bivariate measure *PLV*_{AB}. The phase concentration parameter γ_{AB}, describing the width of the distribution of A–B phase differences, can be found from *PLV*_{AB} via *Eq. 19* and is represented as the length of the black arrows in Fig. 2*B*. The von Mises distributions of *A-B* phase differences, generated from the phase concentration values via *Eq. 18*, are shown in black in Fig. 2*C*. Clearly, using the PLV estimation results in the estimation of a *spurious coupling* between nodes *A* and *B*.

We also estimated the PCE parameters and isolated distributions for this network (red in Fig. 2, *B* and *C*), the critical difference being that PCE uses the phase time series from all three oscillators in the network to estimate the coupling between A and B, rather examining the phases of *A* and *B* alone. Figure 2*B* in red shows the estimated phase coupling and Fig. 2*C* in red shows the estimated isolated distribution, corresponding to *Eq. 23*. For this case, the PCE parameters clearly show no interaction between oscillators *A* and *B*, which correctly matches the network simulation. Therefore, in this example, where PLV produces a spurious coupling, the PCE correctly estimates the lack of coupling between nodes *A* and *B*. PCE estimates a model in which the interactions between nodes *A* and *B* through the third node explain the correlations observed between *A* and *B*. No direct interaction between *A* and *B* is necessary to explain this observed correlation.

We illustrate two additional networks: *missing coupling* (Fig. 2, *D-F*) and *incorrect phase offset* (Fig. 2, *G–I*). In the missing coupling network with four nodes the coupling between nodes *A* and *B* is κ = 0.5 with μ = 3π/4 (135°), the coupling strength between the unlabeled nodes and *A* and *B* is κ = 0.75 with μ = 3π/4 (135°) and μ = −π/2 (270°) between one node and *A* and *B*, respectively, and with μ = −π/2 (270°) and μ = 3π/4 (135°) between the other unlabeled node and *A* and *B*, respectively. In the incorrect phase offset network with three nodes the coupling between nodes *A* and *B* is κ = 0.3 with μ = 3π/4 (135°), between node *A* and the unlabeled node κ = 0.9 with μ = −3π/4 (225°), and between node *B* and the unlabeled node the coupling is κ = 1.1 with μ = 3π/4 (135°). In each case, *PLV* and/or phase concentration does not reflect the true direct interaction between the indicated oscillators. In contrast, inferring the parameters of the full probabilistic distribution using PCE correctly recovers the true coupling between the indicated oscillators and all other pairs. That is, for these simple networks the estimated PCE parameters reflect the “ground-truth” connectivity of each network, while the bivariate measures do not.

### Multivariate Phase Coupling in Nonhuman Primate LFP Signals

Importantly, these examples of spurious coupling, missing coupling, and incorrect phase offset, due to the inability of bivariate methods to factor the empirical distribution into a product of isolated (direct) and network (indirect) distributions, are not limited to toy examples of simulated oscillator networks. Similar differences between the bivariate approach and the multivariate approach can be seen in real neurophysiology data. For example, Figs. 3 and 4 provide examples of coupling differences seen in LFP data recorded from the macaque motor system. These LFP signals were sampled at 1000 Hz from chronically implanted tungsten microwire arrays [see Ganguly and Carmena (2009) for more details of the task, the surgery, and recording methods used]. Conducted procedures were in compliance with the National Institutes of Health's *Guide for the Care and Use of Laboratory Animals* and approved by the University of California (Berkeley) Institutional Animal Care and Use Committee. To extract estimates of instantaneous phase, LFPs were filtered with a complex Gabor atom (Gaussian envelope) with a center frequency of 32 Hz and a frequency-domain SD of 4.4 Hz, corresponding to the observed motor beta-activity for this subject (Canolty et al. 2010). Figure 3 shows examples of coupling estimates for pairs of channels that appear to exhibit cases of spurious coupling (*A*, *D*, and *G*), missing or reduced coupling (*B*, *E*, and *H*), and incorrect phase offset (*C*, *F*, and *I*). While ground truth is not know with certainty in these cases, the similarity to the simulation cases suggests the need for caution when interpreting the results of bivariate methods. Importantly, the network structure inferred from the bivariate and multivariate methods is quite different.

As with the PLV, PCE can be used to track transient, event-related changes in phase coupling during behavior. Figure 4 shows the chronic microelectrode arrays (Fig. 4*A*) implanted in multiple areas (Fig. 4*B*) of the macaque motor system while subjects engaged in 3,730 trials of a delayed center-out reach task over 19 days (Fig. 4*C*) [see Ganguly and Carmena (2009) for details]. From these electrodes, we selected 20 channels from three different brain areas (bilateral primary motor cortex and left dorsal premotor cortex) for further analysis here. A set of *N* channels has *N*(*N* − 1)/2 unique (unordered) channel pairs; thus, 20 channels yield 190 unique channel pairs. As shown earlier, there is a one-to-one correspondence between the PLV index and the von Mises concentration parameter; after computing the PLV for each channel pair as described above, we can convert this value into a corresponding von Mises parameter to compare it with the PCE concentration parameters.

Critically, the previous examples considered above do not explicitly model a change in PCE parameter values over the time course of an experimental trial. With enough sample points per dimension, however, we can simply perform an independent PCE fitting for an ordered set of time points within the trial (cross-trial fitting). That is, we can perform a PCE fitting at the time of the trial go-cue (0 ms), then a second, independent PCE fitting 10 ms after the go-cue, and a third fitting at 20 ms after the go-cue, and so on. In particular, when fitting the 0 ms PCE (for example), we will include one sample point from each trial for each of the 20 channels for all of the 3,730 trials; the input for the PCE function will be a 20 × 3730 (channels × trials) matrix of phases. As shown in Cadieu and Koepsell (2010), excellent parameter estimates can be obtained for as few as 100 sample points per dimension; here we have 3,730 trials and 20 channels for 186.5 sample points per dimension for each of these independent PCE fittings over the course of the experimental trial. For time-resolved estimates of PCE parameter values in studies with fewer trials, we will need to either include more sample points per trial or reduce the number of channels considered. Alternatively, we can shift the focus from event-related PCE changes, as in this example, to experimental-condition-related PCE changes, as shown for the human ECoG language study discussed in the next section.

This event-locked PCE approach results in several differences from PLV-based methods. Figure 4, *D* and *E*, shows the event-related changes in phase coupling associated with movement onset (0 ms) in the center-out reach task. First, note that the overall strength of coupling is lower for the multivariate PCE method than for the bivariate von Mises method. That is, while both methods explain the observed statistics of the data, the PCE parameters are sparser than the von Mises parameters. Second, both methods show a drop in phase coupling associated with movement, consistent with prior results showing a drop in beta-power and phase coupling during movement (Sporns 2011). That is, even given a relatively small number of trials per dimension (Cadieu and Koepsell 2010), PCE can track rapid changes in instantaneous phase coupling in dynamic functional networks. Finally, Fig. 4*F* shows that the bivariate and multivariate methods differ in their statistical specificity. With the use of the PLV resampling procedure described above, 98.4% (187/190) electrode pairs exhibited sustained phase coupling during reaching. In contrast, under PCE only 49.5% (94/190) electrode pairs exhibited statistically significant phase coupling. As in the simulation studies, it is likely that this sparse pattern of coupling is due to PCE explicitly accounting for indirect network effects, unlike the PLV method. Again, these striking differences in the significance and strength of phase coupling will have a strong effect on the interpretation of dynamic network activity.

### Multivariate Phase Coupling in Human ECoG Signals

While the example above focused on time-dependent event-related changes in coupling, it is possible to focus on task- or stimulus-related differences as well. For example, Fig. 5 shows phase coupling data from a human intracranial study on language comprehension and target detection. The study protocol, approved by the University of California, San Francisco, and University of California, Berkeley, Committees on Human Research, did not interfere with the ECoG recording made for clinical purposes and presented minimal risk to the participating subjects. Briefly, patients undergoing surgical treatment for epilepsy were implanted with a grid of subdural electrodes for clinical purposes. During their hospital stay, some patients volunteered to participate in a simple auditory target detection task. Subjects listened to a string of action verbs (e.g., throw, blow, etc.), acoustically matched but scrambled nonwords (unintelligible), or people's names (e.g., Bob, Mary, etc.). The task was to respond with a button press to the rarely presented target items (people's names) and ignore all other stimuli. This task elicits complex patterns of dynamic activity across space, time, and frequency. Here we focus on changes in theta (4–8 Hz)-phase coupling within the widespread network of electrodes of one subject that show strong target-selectivity (red electrodes in Fig. 5*B*). Target-selective electrodes were defined as those electrodes showing a significant increase in high gamma (80–150 Hz)-power following presentation of targets vs. verbs [for details, see Canolty et al. (2007)]. The mean strength of phase coupling within this network is dependent on the type of stimulus presented, with the strongest coupling occurring after target presentation and the weakest occurring after the presentation of scrambled nonwords (Fig. 5*A*; verbs > nonwords, *P* < 0.01; targets > verbs, *P* < 0.01; randomized permutation test).

Interestingly, while the macaque motor results discussed above show a drop in both beta-band power and beta-phase coupling, here theta-power and phase coupling disassociate. Specifically, in this task theta-power remains steady for nonwords, increases for verb distracters, and decreases for target presentation (Canolty et al. 2007). However, while short-range (<2 cm) theta-phase coupling drops for target presentation, many long-range electrode pairs exhibit an increase in theta-phase coupling (e.g., Fig. 5*C*). This shows that power and phase coupling can change independently, at least over large spatial distances. Critically, these stimulus-specific phase coupling results are only apparent using PCE; the PLV and von Mises approaches do not reveal these changes (Fig. 5*C*). Similarly, Fig. 5*D* shows a case where PCE finds a stimulus-dependent phase shift for targets but not distracters (as suggested by the communication through coherence hypothesis), while the bivariate methods do not. As before, this is likely due to the inability of the PLV to separate direct and indirect influences on phase coupling, an interpretation bolstered by ground-truth simulation studies that reach similar findings (Cadieu and Koepsell 2010).

## DISCUSSION

### Summary of Advantages Over Existing Methods

The PCE method, employing a multivariate probabilistic approach, exhibits several advantages over the commonly used PLV method. Importantly, we show that PCE accurately captures both direct and indirect (network mediated) coupling between network elements, and we show that the direct coupling terms have a straightforward physical interpretation, representing interaction terms in a dynamical system of coupled oscillators. It is straightforward to model this dynamical system using standard techniques and thus generate synthetic data with second-order statistical properties identical to the recorded data. The PCE model parameters also describe a probability distribution that can be employed for resampling or simulations. We have found empirically that PCE output is sparser than bivariate approaches and generates models with fewer interactions and thus a more parsimonious interpretation of the data. In addition, using the PCE coupling parameters, forward modeling via simulations generates synthetic data with statistical characteristics identical to empirical signals.

A primary purpose of this study is to show how the multivariate phase distribution derived in Cadieu and Koepsell (2010) can be used to estimate task- and event-related changes in phase coupling in neurophysiological signals such as intracortical LFPs, subdural ECoG, and scalp-recorded signals such as the electroencephalogram or magnetoencephalogram and to highlight the advantages of this approach over bivariate methods. However, PCE has proven useful in addressing a variety of other questions. For example, PCE has been employed to characterize the dependence of neuronal spiking upon large-scale patterns of LFP phase coupling (Canolty et al. 2010) and has also been used to estimate multivariate phase-amplitude cross-frequency coupling (Canolty et al. 2012). Any question that requires the examination of phase coupling in multivariate signals will likely benefit from employing PCE, for the reasons detailed above.

### Limitations of the PCE Method and Future Directions

Despite the advantages over the PLV approach described earlier, the PCE method does not fully describe coupling within brain networks. In particular, three limitations of the PCE method described here point toward directions for future development. First, the PCE approach cannot model phase-phase coupling between multiple frequency bands. The approach described here targets single-frequency phase coupling between *N* simultaneously recorded signals. More complex networks involving coupling between *M* distinct frequencies active within *N* channels will require a probabilistic model of different form.

Second, the PCE method does not incorporate information about the power of ongoing oscillations. This is an important limitation in that transient changes in oscillatory power may act to change the coupling strength between nodes. Intuitively, brain areas that exhibit higher oscillatory power may be relatively “louder” than areas with lower oscillatory power, exerting more influence on downstream nodes. Clearly, it is important to keep the effects of phase and power separate, in fact, this was our primary criticism of the linear coherence approach to estimating phase coupling. However, a modeling approach that incorporates both power and phase as separate and distinct elements will likely produce a better fit to observed data than will a phase-only approach. Additionally, when band-limited power is high, we can have more confidence in the stability of phase estimates compared with times of low power. Unfortunately, it is not immediately obvious how to incorporate amplitude or power into a multivariate model similar to PCE; this is an issue that will require future research.

Third, if some oscillators or signals are not included in the estimation procedure, then PCE may not reveal the “ground-truth” coupling between all oscillators. However, PCE will still generate the most parsimonious model of the available data that matches second-order statistics. That is, the rationale for employing maximum-entropy (minimum assumption) models will still hold in this case, and to our knowledge no other method currently available can recover or estimate the coupling to unobserved (hidden) variables.

### Conclusion

Using ground-truth simulations and empirical data sets, we showed how PCE can be used to investigate changes in multivariate phase coupling and detailed the advantages of this multivariate model over bivariate approaches. Given these strong advantages, we contend that PCE is a useful tool for investigating multivariate phase coupling in distributed functional brain networks.

## GRANTS

This work was supported by DARPA Contract N66001-10-C-2008 (to J. M. Carmena), National Science Foundation Grant CBET-0954243 (to J. M. Carmena), National Science Foundation Grant IIS-0917342 (to C. F. Cadieu and K. Koepsell), and National Institute of Neurological Disorders and Stroke Grant NIH-R01-NS21135 (to R. T. Knight).

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).