Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Conf Proc IEEE Eng Med Biol Soc. Author manuscript; available in PMC 2013 August 26.
Published in final edited form as:
PMCID: PMC3753031

A Beamforming Approach to Phase-Amplitude Modulation Analysis of Multi-Channel EEG

A.L. Sampson,1,* B. Babadi,1,2 M.J. Prerau,1,2 E.A. Mukamel,3 E.N. Brown, Fellow, IEEE,1,2,4 and P.L. Purdon, Member, IEEE1,2,*


Phase-amplitude modulation is a form of cross frequency coupling where the phase of one frequency influences the amplitude of another higher frequency. It has been observed in neurophysiological recordings during sensory, motor, and cognitive tasks, as well as during general anesthesia. In this paper, we describe a novel beamforming procedure to improve estimation of phase-amplitude modulation. We apply this method to 64-channel EEG data recorded during propofol general anesthesia. The method improves the sensitivity of phase-amplitude analyses, and can be applied to a variety of multi-channel neuroscience data where phase-amplitude modulation is present.


Oscillations are thought to underlie many aspects of brain function, but the mechanisms by which these oscillations organize neural activity across different temporal and spatial scales remains an area of active investigation. Recently, cross-frequency coupling has been observed where the phase of theta oscillations (4-8 Hz) modulates the amplitude of gamma oscillations (> 30 Hz) [1]. Similar phase amplitude relationships have been observed during different sensory, motor, and cognitive tasks [2]. These phase-amplitude relationships are usually estimated from single channel data, even when multi-channel data are acquired simultaneously, using non-parametric models of the relationship between phase and amplitude (e.g., constructing a histogram of amplitudes across discrete phase bins). The efficiency of these analyses could be improved substantially if data could be incorporated across multiple channels, and if parametric representations could be used to model the phase-amplitude relationship.

In recent work, we studied phase-amplitude modulation during general anesthesia. We examined how the phase of slow oscillations (< 1 Hz) influenced the amplitude of alpha-band oscillations (8-14 Hz), and found distinct patterns of modulation corresponding to different levels of consciousness under general anesthesia [3]. In this work, we describe a novel beamforming procedure to improve estimation of phase-amplitude modulation. In this method, we model the phase-amplitude relationship parametrically using a low-order Fourier series. We then identify the optimal linear combination of channels to minimize the quadratic error for this phase-amplitude modulation curve. We apply this method to 64-channel EEG data acquired during induction of general anesthesia with the drug propofol. The method improves the sensitivity of phase-amplitude analyses, and can be applied to a broad range of multi-channel neuroscience data, such as EEG or local field potentials (LFP), where phase-amplitude modulation is present.


A. EEG Recordings

We induced and maintained general anesthesia in healthy volunteers using the intravenous anesthetic propofol. The anesthetic induction was carried out by increasing the targeted effect-site propofol concentration to levels of 0, 1, 2, 3, 4, and 5 μg/ml every 14 minutes with a computer controlled infusion pump [4, 5]. We recorded 64-channel EEG continuously during this time (BrainAmp MRPlus, BrainProducts, GMBH). These studies were approved by the Massachusetts General Hospital Human Research Committee. In this paper, we analyze a subset of the data including n = 2 subjects to demonstrate our method for multi-channel estimation of phase-amplitude modulation between the slow (0.1-1 Hz) and alpha (8-14 Hz) bands.

B. Beamforming Formulation of Phase-Amplitude Coupling

Let x(t) := [x1(t), x2(t), ···, xN(t)]T denote the EEG time-series corresponding to EEG channels for time 0 ≤ tT. Let α(t) := [α1(t), α2(t), ···, αN(t)]T and s(t) := [s1(t), s2(t), ···, sN(t)]T denote the alpha rhythm and slow oscillation time-series which are obtained by band-pass filtering x(t) in the frequency bands of 8-14 Hz and 0-1Hz, respectively [3]. Let αH(t)α(t)+iH(α(t)) and sH(t)s(t)+iH(s(t)), where H() is the discrete-time Hilbert transform. In [3] we showed that the amplitude of the alpha rhythm is modulated by the phase of the slow oscillation during general anesthesia, based on an analysis of single-channel Laplacian-derived EEG. Assuming that the phase-amplitude modulation arises from a unified and possibly spatially localized mechanism in the brain, the problem reduces to reconstructing a single phase-amplitude modulation relationship based on the observation through the multi-channel array of EEG sensors.

This problem is well-studied in array signal processing, and a viable solution is given by beamforming [6]. The idea of beamforming is to form a scalar signal based on the array observations in order to minimize an appropriate cost function representing the underlying system model. Let w := [w1, w2, ···, wN]T denote a weight vector (beamforming vector) and consider the corresponding projection of the alpha rhythm and slow oscillation time-series given by αw(t) := wT αH(t), and sw(t) := wT sH(t), respectively. The amplitude of αw(t) and the phase (argument) of sw(t)are given by




respectively. Suppose that for a given value of the phase of sw(t), denoted by θ, the amplitude of the alpha rhythm Aw(t) has a distribution given by the density pw(A; θ). Then, the phase-amplitude modulation relation is defined as


where the ensemble averaging Epw is with respect to the density pw(A; θ). The function Aw(θ; t) is clearly periodic with the full period defined as [–π, π]. We further assume that Aw(θ; t) is stationary during the observation period [0, T] and hence drop the dependence on t. Assuming that the function Aw(θ) has sufficient smoothness properties, it can be represented in the Fourier basis as follows:


where μ, ak, and bk denote the expansion coefficients. A suitable model for estimating Aw(θ) is given by its truncated Fourier expansion to the first L terms, with L ≤ 3. This reduced-order model enforces a smooth phase-amplitude modulation relation, which is consistent with empirical observations in [3]. A suitable cost function for estimating Aw(θ) is given by the following quadratic form:


Since the densities pw(A; θ) are unknown, it is not possible to compute Aw(θ) = Epw{Aw(t)|θ}. Hence, we resort to an empirical quadratic cost function for estimating Aw(θ) by substituting the ensemble averaging operator Epw by the corresponding temporal averaging as follows:


where tw(θ) denotes the inverse function of θw(t), p(θ) is the prior distribution of the slow oscillation phase and Et denotes temporal averaging. Note that we are implicitly assuming the ergodicity of the underlying processes during the observation period [0, T] and hence are replacing ensemble averaging by temporal averaging. Since the prior p(θ) is unknown, the cost function can be further approximated by substituting the ensemble averaging over θ by the corresponding temporal averaging as follows:


For a given beamformer w, it is possible to minimize the cost function over the parameters μ, ak, and bk. Then, we can choose the best such beamformer by minimizing the resulting cost function over w. This, in fact, corresponds to a cost minimization formulation for estimating the reduced-order phase-amplitude modulation relation that is most consistent with the data (in the sense of the above quadratic cost function). Assuming that the beamformer elements are bounded as wwkw, for some constants w and w, the overall optimization procedure can be expressed as:


The inner minimization can be easily carried out by linear regression and the resulting solution can be expressed explicitly in terms of Aw(t) and θw(t). The outer minimization can be performed using standard optimization routines. In particular, since the constraints on wk form a convex set, we have employed the interior point method for the outer minimization stage.

C. Data Analysis

For this analysis, the data from two subjects was used to compute optimal weight coefficients w. In our earlier studies, we showed that the phase-amplitude modulation of frontal EEG under GA [3] undergoes two different patterns of modulation, corresponding to depth of anesthesia. The first pattern, occurring before and after the point of loss of consciousness, consists of maximum alpha amplitude occurring at the trough (surface-negative) of the slow oscillation, which we refer to as the “trough-max” pattern. Under a deeper level of GA, the relationship reverses and maximum alpha amplitude occurs at the peak (surface-positive) of the slow oscillation, which we refer to as the “peak-max” pattern. In order to compute the electrode weights that would show both modes of the phase-amplitude modulation most clearly, equal-length segments of data from both modes were chosen and used to compute the optimal weights for each mode. These trough-max and peak-max data for the two subjects was were used to perform the averaging described in Equation (7). The data used in the optimization consisted of four-minute segments, chosen as periods during which the phase-amplitude modulation was relatively constant, based on phase-amplitude histograms computed using Laplacian-referenced data.


The alpha amplitude as a function of slow wave phase was fit for each of these data sets using the Fourier model in Equation (4), as well as a non-parametric model using 100 phase bins (Figure 1A). The Fourier model order was chosen using an F-test for inclusion of successive terms. Based on this test, a first order model was used everywhere except in the peak-max data segment from subject 2, where a second order model provided a better fit to the data.

Fig. 1
For trough-max modulation (A), and peak-max modulation (B), estimates of the empirical phase/amplitude relationship (blue) are shown along with first and second order model fits. The optimal weights for each electrode are shown to the right as estimated ...

The results of the optimization in all cases tended to give the greatest weight to one or two electrodes, and included smaller contributions from the other channels. The values of the weights for each electrode for subject 2 are shown in Figure 1, using the first Fourier harmonic for the trough-max period, and the first and second Fourier harmonics for the peak-max period.

To assess the efficacy of the beamformer weighting, we computed time-varying non-parametric phase-amplitude histograms using: 1) a single frontal channel with bipolar reference, 2) a single frontal channel with Laplacian reference (i.e., using the average of neighboring electrodes as the reference, as in [3]), and 3) the proposed method with optimal weights. These phase-amplitude histograms were computed using 100 phase bins, averaging over 1-minute windows with 30 seconds overlap, normalized by average alpha amplitude within each window. Figure 2 shows these time-varying phase-amplitude histograms for both subjects during periods of both trough-max and peak-max modulation. The modulation relationship appears clearest when the optimal weights w are used. To characterize this quantitatively, we computed the modulation depth for each of the methods by taking the average phase-amplitude histogram in each period (trough-max and peak-max), and taking the difference between the maximum and minimum points of the histogram. The beamforming method produced the largest modulation depth, followed by the Laplacian method, with bipolar referencing showing the lowest modulation depth in both regimes.

Fig. 2
(A) The targeted effect-site propofol concentration is shown, along with the time-varying phase-amplitude histogram computed based on Laplacian-referenced data for a period spanning the trough-max and peak-max regimes of modulation. (B) Phase-amplitude ...


Phase-amplitude modulation has been observed in a number of oscillatory neural systems sensory, motor, or cognitive tasks [2]. The modulation of alpha wave amplitude by slow wave phase is present in the EEG during general anesthesia, and could be an important variable for assessing brain activity and level of consciousness during general anesthesia. In EEG data, the ability to detect cross-frequency coupling is strongly influenced by the electrode reference scheme, and is often difficult to detect with a standard bipolar reference. The beamforming method presented here provides a means to obtain electrode weights that minimize the least-squares error in a parametric sinusoidal model of the phase-amplitude relationship. This optimal weighting of EEG electrodes allows for improved detection of phase-amplitude modulation across time and subjects. This method could be useful in studies of phase-amplitude modulation in the EEG under general anesthesia, as well as other conditions where this phenomenon might arise.

Modulation Depth For Different Methods


Research supported by NIH grants DP2-OD006454 (Purdon), K25-NS057580 (Purdon), DP1-OD003646 (Brown), and R01-EB006385 (Brown).


1. Canolty RT, Edwards E, Dalal SS, Soltani M, Nagarajan SS, Kirsch HE, Berger MS, Barbaro NM, Knight RT. High gamma power is phase-locked to theta oscillations in human neocortex. Science. 2006 Sep 15;313:1626–8. [PMC free article] [PubMed]
2. Canolty RT, Knight RT. The functional role of cross-frequency coupling. Trends in cognitive sciences. 2010 Nov;14:506–15. [PMC free article] [PubMed]
3. Mukamel EA, Wong KF, Prerau MJ, Brown EN, Purdon PL. Phase-based measures of cross-frequency coupling in brain electrical dynamics under general anesthesia. IEEE EMBS. 2011 Aug;:1981–4. 2011. [PMC free article] [PubMed]
4. Shafer SL, Gregg KM. Algorithms to rapidly achieve and maintain stable drug concentrations at the site of drug effect with a computer-controlled infusion pump. J Pharmacokinet.Biopharm. 1992;20:147–169. 04. [PubMed]
5. Schnider TW, Minto CF, Shafer SL, Gambus PL, Andresen C, Goodale DB, Youngs EJ. The influence of age on propofol pharmacodynamics. Anesthesiology. 1999;90:1502–1516. 06. [PubMed]
6. Van Trees HL. Optimum Array Processing. Wiley; New York: 2001.