## Abstract

Although the thalamus is believed to regulate and coordinate cortical activity both within and across functional regions, such as motor and visual cortices, direct evidence for such regulation and the mechanism of regulation remains poorly described. Using simultaneous invasive recordings of cortical and thalamic electrophysiological activity in 2 awake and spontaneously behaving human subjects, we provide direct evidence of thalamic regulation of cortical activity through a mechanism of phase-amplitude coupling (PAC), in which the phase of low frequency oscillations regulates the amplitude of higher frequency oscillations. Specifically, we show that cortical PAC between the theta phase and beta amplitude is spatially dependent on and time variant with the magnitude of thalamocortical theta coherence. Moreover, using causality analysis and MR diffusion tractography, we provide evidence that thalamic theta activity drives cortical theta oscillations and PAC across structures and that these thalamocortical relationships are structurally constrained by anatomic pathways. This relationship allows for new evidence of thalamocortical PAC. Given the diffuse connectivity of the thalamus with the cerebral cortex, thalamocortical PAC may play an important role in addressing the binding problem, including both integration and segregation of information within and across cortical areas.

## Introduction

Just like a conductor works reciprocally with an orchestra to receive, integrate, and coordinate components, the brain can be hypothesized to require a conductor to coordinate multiple inputs and outputs, cognitive processes, and attentional factors. The diffuse yet specific patterns of structural connectivity of the thalamus with the cerebral cortex suggests that thalamocortical connectivity could play an important role in regulating cortical activity. A regulatory role of the thalamus is strongly suggested by the coherence of cortical sleep spindles and alpha activity with thalamocortical interplay activity and intrathalamic feedback (McCormick and Huguenard 1992; Golomb et al. 1994; Contreras et al. 1996; Bazhenov et al. 1999; Goncalves et al. 2006; Liu et al. 2012). In particular, studies using EEGfMRI and fMRIfMRI resting state correlation analyses have implicated the pulvinar nucleus of the thalamus in generating and modulating alpha rhythms in the occipital lobe, between which there are extensive reciprocal thalamocortical connections (Goncalves et al. 2006; Liu et al. 2012). Other studies have utilized theoretical methods to indicate the presence of specific alpha rhythm generators at the thalamic level (McCormick and Huguenard 1992; Destexhe et al. 1993; Golomb et al. 1994; Contreras et al. 1996; Bazhenov et al. 1999). The role of the pulvinar in regulating cortical function was further elucidated in a nonhuman primate study which demonstrated that the pulvinar synchronizes activity between interconnected cortical areas according to attentional allocation, suggesting a critical role for the thalamus not only in attentional selection but more generally in regulating information transmission across visual cortices (Saalmann et al. 2012). Staudigl et al. (2012) investigated thalamocortical communication during human long-term episodic memory retrieval and observed the impact of such thalamocortical communication on local frontal networks is expressed via a modulation of gamma power by the thalamic beta (∼20–23 Hz) phase. And in a more recent study, FitzGerald et al. (2013) used data from epileptic patients undergoing thalamic deep brain stimulation and reported existence of PAC both within thalamus and prefrontal cortex (PFC) and between them. Thalamocortical regulation, however, is not limited to the pulvinar. For example, the ventrolateral and mediodorsal nuclei of the thalamus are believed to be important components of the frontal cortical-basal ganglia–thalamic circuits mediating motivation and emotional drive, planning, and the expression of goal-directed behaviors (Haber and Calzavara 2009).

A putative electrophysiological mechanism for thalamocortical regulation is cross-frequency coupling (CFC), and more specifically phase-amplitude coupling (PAC), in which the phase of a low frequency rhythm from one signal regulates the power of higher frequency activity (either from the same or another signal). PAC has been described extensively as an inherent property of cortical electrophysiology and is postulated to play a role in regulation of an array of neural networks, including memory and learning, attention, and in sensory and motor processing (Bhattacharya 2001; Fries 2005; Lakatos et al. 2005). While many works describe the existence of PAC within a cortical region, the literature has rarely described the spatial specificity of this phenomenon, the time-variant nature of these complex relationships, nor the factors that may regulate such variation (Lakatos et al. 2005; Henriksson et al. 2009; Tort et al. 2009; Axmacher et al. 2010; Canolty and Knight 2010; Lopez-Azcarate et al. 2010; Pienkowski and Eggermont 2010). Specifically, the origin of the phase-encoding frequency remains speculative. PAC provides a plausible physiological explanation for such subcortico-cortical mechanisms that require dynamic coordination of different frequencies with different spatial properties such as lower frequencies like theta or alpha rhythms (from subcortical region) and higher frequencies like beta or gamma rhythms (cortical sites) (Canolty and Knight 2010). Miller et al. (2010) provide the most comprehensive report to date of temporal modulation of cortical PAC that is behaviorally-dependent. The authors identify the thalamus as a putative synchronizing locus that could provide the phase-encoding low frequency rhythms that modulate cortical activity. Given differences in the spatio-temporal profiles of lower and higher frequency electrophysiological rhythms (e.g., it has been suggested that lower frequencies like theta modulate activity over large spatial regions and in long temporal windows while higher frequencies modulate activity over small spatial regions and short temporal windows activity) (Canolty and Knight 2010), the possibility of distinct sources for these signals is plausible and in fact probable. Changes in thalamic neuronal activity between tonic (Henning Proske et al. 2011) and bursting (Llinas and Jahnsen 1982; McCormick and Huguenard 1992) modes of activity are thought to create electrophysiologic oscillations that spread widely through both intrathalamic and thalamo-cortical connections, creating characteristic thalamocortical rhythms that could modulate cortical activity. Pathological thalamocortical rhythms, or thalamo-cortical dysrhythmia (TCD), have been hypothesized to be central to electrophysiological changes in some functional brain disorders (Llinas et al. 1999; Sarnthein et al. 2005, 2006; Walton and Llinas 2010). TCD posits that abnormal internally generated low frequency oscillations (mainly in theta range) in the thalamocortical network disrupts the normal state dependent flow between thalamus and cortex resulting in a broad range of functional disorders, including Parkinson's disease, depression, and chronic pain (depending on the thalamocortical network involved) (Llinas et al. 1999; Llinas and Steriade 2006; Kane et al. 2009; Jones 2010). Still, direct evidence for a subcortical source for low frequency rhythms is lacking to date and therefore remains speculatory.

Multimodality analyses including both structural and functional data within subjects provide unique opportunities to elucidate the nature of thalamocortical relationships. These analyses, however, hinge on the assumption that structural connectivity analyses, especially that of thalamocortical relationships, are functionally valid. In a seminal study using magnetic resonance diffusion tensor imaging (DTI) and probabilistic tractography, Johansen-Berg and colleagues reported anatomically specific and histologically concordant segmentation of the human thalamus based on each thalamic subregion (or nucleus) having a distinct pattern of cortical connectivity (Behrens et al. 2003). Subsequent studies functionally validated this approach, demonstrating that thalamic functional activations during a motor task (detected using functional magnetic resonance imaging) colocalize with the thalamic regions with the highest probability of connectivity with motor and prefrontal cortical areas, respectively (Behrens et al. 2007). Direct comparisons between structure and function within subjects however provide much stronger validation. This opportunity is uniquely afforded by neurosurgical procedures that provide access to invasive local field potential (LFP) recordings from multiple sites within the brain in patients who have also undergone DTI. Using this approach, we recently showed that probabilistic tractography can be used to segment the thalamus and identify specific thalamic subregions to be targeted for deep brain stimulation that will result in tremor suppression, suggesting that using this method to segment brain structures is in fact functionally meaningful and reliable across subjects (Pouratian et al. 2011). In a separate study of patients undergoing intracranial EEG monitoring with intrathalamic depth and subdural cortical strip electrodes, we validated thalamocortical tractography by tracking somatosensory evoked potentials (SSEP) through the thalamus and to the cortex in a manner that was concordant with predictions based on tractography (Elias et al. 2012). While these correlations are suggestive, the thalamic subregions identified by tractography do not necessarily correspond to thalamic nuclei or subnuclei.

While electrophysiological activity and functional relationships should theoretically correspond well with structural connectivity, the strength or robustness of an anatomical connection is not always reflective of its functional significance. This is particularly true when investigating the potential role of one region in modulating or regulating another. To better elucidate the putative role of PAC as a mechanism of thalamocortical regulation, we analyzed and compared simultaneous cortical and thalamic electrophysiological recordings in 2 awake and spontaneously behaving human subjects who were undergoing invasive neurophysiological monitoring for medically refractory epilepsy who had had preoperative DTI. We first extensively characterize cortical PAC to illustrate the time-variant nature and spatial specificity of these signals. The time and spatial variance of cortical PAC are compared with simultaneously recorded thalamic LFP power spectra to evaluate the relationship between PAC and thalamocortical coherence of the phase-encoding rhythm. We subsequently use causality analysis and probabilistic tractography analyses to assess the thalamocortical flow of phase-encoding rhythms and the structurally constraints on this flow of information.

## Materials and Methods

### Signals and Recording

We studied 2 subjects with intractable epilepsy who underwent invasive monitoring to identify epileptogenic foci. Subdural electrocorticographic strips, each containing four to eight 6 mm diameter contacts with 10 mm spacing (AD-Tech, Racine, WI), were implanted through standard frontal and parietal burr holes, per clinical protocol. After obtaining informed consent from the study participants, depth electrodes were stereotactically implanted subacutely in the thalamus based on a research protocol approved by the University of Virginia Institutional Review Board originally intended to understand the role of thalamus in seizure propagation (Bertram et al. 2001). Each depth electrode consisted of tubing (1.1-mm outer diameter) with 10 cylindrical platinum contacts of 2.3 mm length with an interelectrode distance of 5 mm (AD-Tech). All recordings were done with 200 Hz sampling rate with patients spontaneously behaving (based on simultaneous video recordings). One hour epochs of data were sampled from each subject. Data were extracted from seizure free periods (based on expert neurologist interpretation). Likewise, episodes with noise contamination were identified by visual inspection of raw traces demonstrating non-physiologic patterns and abnormally high amplitudes and spliced out of the data. To minimize effects of volume conduction, the data were used in a bipolar montage, with the time courses from adjacent electrodes subtracted from one another.

### Power Spectral Analysis

Before further processing 60 Hz line noise was removed from the data and all following steps were done using a moving window of 10 s length and 50% overlap. To estimate the power spectral density of both cortical and thalamic signals, we used the multitaper method implemented in Chronux (Bokil et al. 2010). To circumvent the limitations of conventional Fourier analysis in estimation of power spectra (which introduces undesirable bias in its estimate in the setting of noise), the multitaper method uses mutually orthogonal tapers (which are multiplied element-wise by signals) providing multiple independent estimates of spectra (called tapered spectra). The final spectrum is obtained by averaging these tapered spectra.

### Phase-Amplitude Coupling

PAC we estimated using Tort's relative entropy method (Tort et al. 2010), which we refer to as Modulation Index (MI). Given signals $Xraw,1(t)$ and $Xraw,2(t)$ containing phase and amplitude components of interest, respectively, we used wavelet transform (Morlet packet with width 7) to extract the instantaneous phase and amplitude signals, respectively, providing $φph(t)$ and $Aamp(t)$. $φph(t)$ phases are then binned and the mean of amplitude ($Aamp(t)$) over each bin is calculated and normalized and referred to as $P(j)$ (for $j=1,2,…,N$ where $N$ is number of bins). $P$ has the characteristics of a probability density function and can be referred to as an “amplitude distribution” when plotted as a function of the phase bins. When there is strong coupling/modulation between phase and amplitude signals this distribution deviates from the uniform distribution. Therefore we define MI to be the distance between uniform distribution and.$P$. Mathematical formulation for MI would be as follows:

$D(P,U)=∑j=1N⁡P(j)logP(j)U(j),$

where $U$is the uniform distribution defined over the same phase bins. $D(P,U)$ is the KL distance between $P$ and uniform distribution. Because

$D(P,U)=log⁡(N)+∑j=1N⁡P(j)log⁡(P(j)).$

Based on the definition of Shannon entropy we can define $H(P)$ to be:

$H(P)=−∑j=1N⁡P(j)log⁡(P(j)).$

Therefore $D(P,U)$can be rewritten as:

$D(P,U)=log(N)−H(P).$

And finally we put:

$MI=D(P,U)log⁡(N)=1−H(P)log⁡(N).$

Because $H(U)=log⁡(N)$, when there is no coupling and $P$ is in fact a uniform distribution, $H(P)$ would be equal to $log⁡(N)$ and $MI=0$.

Where there is stronger coupling $H(P)$ becomes smaller than $log⁡(N)$and MI increases. We calculated MI using frequencies ranging 1–70 Hz in 1-Hz steps for the phase-encoding signals and 2-Hz steps for amplitude (power) signal.

### Surrogate Data Analysis

A surrogate data analysis using a shuffling procedure was used to evaluate the significance of derived MI values. For each signal pair ($φph(t)$, $Aamp(t)$), we generated 100 temporally shuffled versions of amplitude signals and calculated MI values for each. We then compared the true MI relative to the 100 surrogate MI values and obtained a Z-score. Only MI values with corresponding Z-scores > 1.96 (corresponding to a P-value < 0.05) were maintained for subsequent consideration. MI values that were non-significant (i.e., Z-value < 1.96) were not included in MI maps.

### Cross Coherence Analysis

After obtaining power spectra ($Sx(f)$ and $Sy(f)$) and also cross spectra ($Sxy(f)$ which is a measure of joint power of 2 signals x and y per unit frequency) using the multitaper method, we estimated cross coherence $Cxy(F)$ as:

$Cxy(f)=|Sxy(f)|Sx(f)Sy(f).$

We calculated the coherence for pairs of signals in thalamocortical network and derived cross coherence as a measure of frequency, which describes degree of covariability between the 2 signals over different frequency ranges. To find the time–frequency representation of coherence for the same pair of signals we used a moving window approach to find the coherence at each frequency over time. This form of time–frequency representation allows for observing pattern of change in coherence over time.

### Granger Causality and Directed Transfer Function Analysis

Directional transfer function (DTF) is a derivation of Granger causality (GC), which is a data-driven approach to assess the causal relationship (directional causal interaction) between 2 time series (Kaminski and Blinowska 1991). This method has been widely used in the analysis of both LFP and intracranial recordings (Brovelli et al. 2004; Chen et al. 2006; Bressler et al. 2007; Wang et al. 2007, 2008; Bollimunta et al. 2008; 2008; Gow et al. 2009; Tass et al. 2010; Ding et al. 2011; Zhang et al. 2012). However, GC limits analysis to 2 time series. We therefore used the adapted multivariate version of it, the directed transfer function to further characterize thalamocortical functional relationship. Causality was evaluated using eConnectome (Electrophysiological Connectome) (He et al. 2011) toolbox, an open source Matlab software package that has been specially developed to investigate directional interactions between multiple electrophysiological signals. DTF takes into account all signals simultaneously and makes possible estimation of activity flow in a given direction as a function of frequency. DTF is robust in respect to noise and constant phase disturbances; in particular it discriminates against volume conduction, which propagates with the zero phase (Kaminski et al. 2001). Normalized version of DTF (which has been used here) is calculated as previously reported (Kaminski et al. 2001):

$DTFa→b2(f)=|Hba(f)|2∑i=1ch⁡|Hbi(f)|2.$

In which $Hba(f)$ is the $(a,b)$element in the MultiVariate AutoRegressive (MVAR) model Matrix for the system of signals and $DTFa→b2(f)$ is a number between 0 and 1 producing the ratio between inflow from channel $a$ to channel $b$ to all the inflows to channel $b$.

### Diffusion Tractography Probabilistic Segmentation

Probabilistic diffusion tractography was used to define patterns of structural connectivity between regions of interest, using methods previously described in detail using FSL tools (FMRIB's Diffusion toolbox (FDT); http://www.fmrib.ox.ac.uk/fsl) (Johansen-Berg et al. 2005). FDT (BEDPOSTX) uses Bayesian techniques to estimate a probability distribution function (PDF) on the principal fiber direction at each voxel, accounting for the possibility of crossing fibers within each voxel. Two fibers modeled per voxel, a multiplicative factor (i.e., weight) of 1 for the prior on the additional modeled fibers, and 1000 iterations before sampling (Behrens et al. 2007). Eddy current correction was used to apply affine registrations to each volume in the diffusion dataset to register it with the initial reference B0 volume prior to performing tractography. Skull stripping was performed using the brain extraction tool (BET). Using these PDFs and PROBTRACKX, we could then determine the probability of connection between seed voxels (in the desired cortical strip) and the predefined thalamic targets chosen to cover the inter-space between adjacent thalamic contacts (using 5000 samples, a 0.2 curvature threshold, and loopcheck termination). The cortical seed (corresponding to the location of ECoG recordings) was then segmented into distinct regions on a voxel-by-voxel basis based on the thalamic target with which each cortical voxel was most dominantly connected.

## Results

Simultaneously acquired ECoG and thalamic depth recordings of LFP in 2 awake and spontaneously behaving patients undergoing invasive electrophysiological monitoring for medication-refractory epilepsy were evaluated. Multiple thalamocortical pairs were evaluated in each subject with similar patterns observed for each thalamocortical pair studied. For descriptive and analytic purposes, detailed results are provided for a single thalamocortical pair in one subject. Results of the other thalamocortical pairs and the other subject are provided in Supplementary Figures 1 and 2, respectively.

### Temporal Dynamics and Spatial Variability of Cortical Power Spectra and PAC

Cortical ECoG signals demonstrated time-variant spectral power (Fig. 1A) with pronounced variability occurring in the theta band (19.92 ± 3.54 dB, Fig. 1B). In some contacts (see Discussion of spatial variability below), significant PAC between the phase of theta frequencies and the amplitude of beta rhythms was noted (P < 0.05, Fig. 1C,D). The phase-encoding frequency corresponded to the band containing the most time-dependent variability in power and the peak frequency within this band. At contacts with significant theta-beta PAC, the magnitude of the PAC was temporally dynamic (3.54 × 10−4 ± 3.62 × 10−4, Fig. 1E). While PAC was noted in several cortical contacts, it was not uniformly observed across all contacts (Fig. 2). In fact, the presence of PAC was related to the presence of a theta peak in the PSD of the ECoG signals at each contact (Fig. 2A–C, contacts 3–5).

Figure 1.

(A) Time–frequency representation of power spectral density for the cortical contact 4 from Figure 2A. Color bar represents power in dB. (B) Time-dependent variability of distinct frequency bands, demonstrating variability in theta power. (C,D) Assessment of PAC within cortical contacts consistently revealed significant PAC between the phase of theta frequencies and the amplitude of beta rhythms (color bar represents MI (unitless)), with peak beta amplitudes occurring contemporaneously with theta troughs (E). Theta-beta PAC, like theta power, was temporally dynamic, with as much as 3.54 × 10−4 ± 3.62 × 10−4 variability over time.

Figure 1.

(A) Time–frequency representation of power spectral density for the cortical contact 4 from Figure 2A. Color bar represents power in dB. (B) Time-dependent variability of distinct frequency bands, demonstrating variability in theta power. (C,D) Assessment of PAC within cortical contacts consistently revealed significant PAC between the phase of theta frequencies and the amplitude of beta rhythms (color bar represents MI (unitless)), with peak beta amplitudes occurring contemporaneously with theta troughs (E). Theta-beta PAC, like theta power, was temporally dynamic, with as much as 3.54 × 10−4 ± 3.62 × 10−4 variability over time.

Figure 2.

(A) From left to right: electrode placement for thalamic and cortical contacts (B) Power Spectral density for cortical contacts placed according to panel A. (the presence of theta peak is indicated by the downward arrow) (C) Intracortical PAC within same cortical contacts as in panel B in which theta-beta PAC is observed in cortical signals with theta peak in their power spectral density (as in panel B). (D) coherence between thalamic (columns) and cortical (rows) pairs. (E) PAC between phase of thalamic signals and amplitude of cortical signals. Note the peak of coherence in theta frequency range between thalamic signals and cortical signals with significant theta-beta PAC (red box) (using the same scale as panel (C) for PAC).

Figure 2.

(A) From left to right: electrode placement for thalamic and cortical contacts (B) Power Spectral density for cortical contacts placed according to panel A. (the presence of theta peak is indicated by the downward arrow) (C) Intracortical PAC within same cortical contacts as in panel B in which theta-beta PAC is observed in cortical signals with theta peak in their power spectral density (as in panel B). (D) coherence between thalamic (columns) and cortical (rows) pairs. (E) PAC between phase of thalamic signals and amplitude of cortical signals. Note the peak of coherence in theta frequency range between thalamic signals and cortical signals with significant theta-beta PAC (red box) (using the same scale as panel (C) for PAC).

### Temporally Dynamic and Spatially Specific Thalamocortical PAC

To elucidate relationships between thalamocortical activity and cortical electrophysiological dynamics, we initially evaluated thalamocortical coherence as a measure of thalamocortical functional connectivity between the cortical contacts showing or being adjacent to those with strong PAC and simultaneously recorded LFP from throughout the thalamus. Indeed, inspection of the data suggested cortical contacts demonstrating theta peaks (Fig. 2B) and intracortical PAC (Fig. 2C, contacts 3–5) seemed to demonstrate increased thalamocortical theta coherence with certain thalamic contacts, suggesting a link between thalamic and cortical theta activity. The relationship between thalamocortical theta coherence and cortical theta power and intracortical PAC is analyzed more quantitatively in subsequent analyses (see next paragraph and Fig. 3). Furthermore, thalamocortical pairs demonstrating significant thalamocortical PAC (Fig. 2E, red box), in which the phase of the theta frequency band from thalamic signals modulates the amplitude of the beta band activity recorded from cortical contacts, were also found to have strong thalamocortical theta coherence (corresponding to the same contacts in which we observed intracortical PAC previously described in Fig. 2C). The degree of thalamocortical coherence was in fact significantly correlated with the strength of thalamocortical PAC across thalamocortical pairs (R2 = 0.4544 with P = 2.2 × 10−4).

Figure 3.

(A) Average theta band thalamo-cortical coherence over time. (B) Scatter plot between average theta coherence and average cortical theta power (R2 = 0.3379, P = 4.0898 × 10−15). (C) Scatter plot between thalamo-cortical theta coherence and thalamo-cortical PAC (R2 = 0.4138, P = 3.9989 × 10−19). (D) Average cortical theta power over time. (E) Average Thalamo-cortical PAC over time. (F) Scatter plot between average cortical theta power and cortico-cortical PAC (R2 = 0.4988, P = 2.8901 × 10−24). (G) Scatter plot between thalamo-cortical and cortico-cortical PAC (R2 = 0.4138, P = 8.4986 × 10−14). (H) Average cortico-cortical PAC over time.

Figure 3.

(A) Average theta band thalamo-cortical coherence over time. (B) Scatter plot between average theta coherence and average cortical theta power (R2 = 0.3379, P = 4.0898 × 10−15). (C) Scatter plot between thalamo-cortical theta coherence and thalamo-cortical PAC (R2 = 0.4138, P = 3.9989 × 10−19). (D) Average cortical theta power over time. (E) Average Thalamo-cortical PAC over time. (F) Scatter plot between average cortical theta power and cortico-cortical PAC (R2 = 0.4988, P = 2.8901 × 10−24). (G) Scatter plot between thalamo-cortical and cortico-cortical PAC (R2 = 0.4138, P = 8.4986 × 10−14). (H) Average cortico-cortical PAC over time.

Like cortical theta power and PAC (Fig. 1B and E), the spatially specific and concordant thalamocortical theta coherence and theta-beta PAC were temporally dynamic (Fig. 3A,E) and significantly correlated in time with one another (Fig. 3C, with R2 = 0.4138 and P = 3.9989 × 10−19). Thalamocortical theta coherence was also highly correlated with cortical theta power (Fig. 3A,B,D, R2 = 0.3379, P = 4.0898 × 10−15), which itself was highly correlated with intracortical theta-beta PAC (Fig. 3D,H, with R2 = 0.4988 and P = 2.8901 × 10−24). Finally, thalamocortical and intracortical theta-beta PAC were likewise highly correlated (Fig. 3G, with R2 = 0.3209 and P = 8.4968 × 10−14). The time-variant relationships of thalamic theta power was also compared with cortical theta power (R2 = 0.0128 with P = 0.0157), thalamocortical theta coherence (R2 = 0.018 with P = 0.02), and thalamocortical theta-beta PAC (R2 = 0.0261 with P = 5.53 × 10−4). While significant in each case, R-values were much lower than other relationships described above and illustrated in Figure 3. Moreover, the relationship between thalamic theta power and cortical theta power (R2 = 0.0147 with P = 0.0743) as well as between thalamic theta power and thalamocortical theta coherence (R2 = 0.024 with P = 0.34) were not found to be significant in the second subject.

### Theta Rhythms Flow from Thalamus to Cortex

We used causality to evaluate which area is driving (causing) activity in the other one, with a particular focus on the phase-encoding theta band in order to understand if in fact thalamic signals are driving cortical PAC patterns. Causality analysis between thalamic and cortical contacts found to have high theta flow during a 1-min period of high PAC demonstrated temporally dynamic patterns as illustrated in Figure 4A-F with cortico-cortical, thalamo-thalamic, and thalamo-cortical influences but little to no cortico-thalamic flow of signals in the theta band. On average (Fig. 4G), while there is flow of theta signals within thalamus and within cortex, the direction of flow of theta signals between thalamic and cortical signals is from the thalamus to the cortex.

Figure 4.

(A–F) Average theta DTF for thalamocortical grid every 10 s during a 1-min period of high PAC (C indicates cortical and T indicates thalamic). (G) Average theta DTF of panels A–F, demonstrating flow of theta signals within thalamic and within cortical contacts, but only from thalamus to cortex.

Figure 4.

(A–F) Average theta DTF for thalamocortical grid every 10 s during a 1-min period of high PAC (C indicates cortical and T indicates thalamic). (G) Average theta DTF of panels A–F, demonstrating flow of theta signals within thalamic and within cortical contacts, but only from thalamus to cortex.

### Thalamocortical Dynamics are Constrained by Structural Connectivity

The variability of theta power across cortical contacts (Fig. 2A) as well as the variability of coherence and PAC across thalamocortical pairs (Fig. 2D,E) suggest a spatial specificity to the observed variance. Further analyses were therefore done to elucidate the potential contributions of structural connectivity to the observed thalamocortical functional phenomena.

In order to confirm the tight relationship between functional and structural thalamocortical connectivity, we performed probabilistic tractography to determine the relative strength of connectivity between cortical and thalamic recording sites (Fig. 5). The strongest probabilistic structural connectivity was identified between the thalamic and ECoG contacts displaying thalamocortical coherence and PAC (Fig. 4). Specifically, using the entire ECoG strip as a seed mask, we found that the part of ECoG strip that is most strongly structurally connected to the thalamic contacts of interest (defined as the thalamic contacts that demonstrated thalamocortical coherence and PAC [see Fig. 2], as indicated by the target mask as in panels A and B,) corresponds to the 3 cortical contacts demonstrating the time-variant thalamocortical theta coherence and thalamocortical theta-beta PAC detailed in Figures 2 and 3. We calculated the average number of probabilistic tractography “hits” within spherical masks around cortical LFP contacts seeded from thalamic contacts indicated in Figure 5A and used this as a measure of strength of connectivity to those thalamic contacts. Fig 5F shows that the strength of thalamocortical connectivity between the thalamic contacts of interest and each cortical contact strongly correlated with average thalamocortical theta-beta PAC (R2 = 0.9166, P = 0.01).

Figure 5.

(A) Sagittal view of high-resolution structural T1 image illustrating the thalamic target contacts used for probabilistic tractography, including the first 3 thalamic contacts (red ellipse). (B) View of cortical EcoG strip indicated by yellow ellipse along with 3 lines showing where the transverse views of probabilistic tractography maps in panels (C–E) are located. (C–E) Three transverse views indicated by lines 1–3 in panel (B) showing the cortical area in the cortical strip most connected to the thalamic target mask (thalamic contacts 1–3), corresponding to the same thalamocortical pairs found to have high thalamocortical coherence and PAC described in Figures 2 and 4. (F) Scatter plot showing the correlation between the strength of connectivity between thalamic contacts and each cortical contact (as numbered in panel (B)) and the average thalamocortical PAC (R2 = 0.9166 with P = 0.01).

Figure 5.

(A) Sagittal view of high-resolution structural T1 image illustrating the thalamic target contacts used for probabilistic tractography, including the first 3 thalamic contacts (red ellipse). (B) View of cortical EcoG strip indicated by yellow ellipse along with 3 lines showing where the transverse views of probabilistic tractography maps in panels (C–E) are located. (C–E) Three transverse views indicated by lines 1–3 in panel (B) showing the cortical area in the cortical strip most connected to the thalamic target mask (thalamic contacts 1–3), corresponding to the same thalamocortical pairs found to have high thalamocortical coherence and PAC described in Figures 2 and 4. (F) Scatter plot showing the correlation between the strength of connectivity between thalamic contacts and each cortical contact (as numbered in panel (B)) and the average thalamocortical PAC (R2 = 0.9166 with P = 0.01).

### Generalizability to Other Thalamocortical Pairs in Same and Other Subjects

Identical thalamocortical functional and structural connectivity patterns were observed for other pairs of thalamic and cortical electrodes in this same subject (Supplementary Fig. 1). Likewise, consistent patterns of time variant, spatially specific, and structurally constrained thalamocortical functional coupling and modulation were found in the second subject (Supplementary Fig. 2).

## Discussion

We sought to use this unique and difficult to acquire dataset of simultaneously recorded ECoG and thalamic LFP recordings in humans to provide new evidence of thalamic modulation of cortical electrophysiological activity via PAC. In order to do this, we first highlight the rarely discussed point that cortical PAC is in fact not a universal phenomenon, but a spatially specific and temporally dynamic one. The variable nature of PAC across time and across cortical recording sites and thalamocortical pairs is an essential prerequisite to positing that PAC serves as a dynamic regulatory mechanism of cortical activity. This is consistent with the hypothesized function of PAC in regulation of cortical activity, as it pertains to memory, attention, sensory, and motor programming (Bhattacharya 2001; Fries 2005; Lakatos et al. 2005). While PAC has been recognized as a putative regulation mechanism by which subcortical structures may regulate cortical activity (Tort et al. 2009; Lopez-Azcarate et al. 2010; Pienkowski and Eggermont 2010), the etiology and regulation of the phase-encoding frequency has remained theoretical. Our analyses provide the first direct electrophysiological evidence of thalamocortical coupling with causality analysis that implicates thalamic contributions to the modulation of spontaneous cortical activity. This study supports the central role of the thalamus, at least in part, in regulating cortical activity and affirms conclusions of other studies based on indirect evidence that thalamus is in fact modulating cortical activity. While the results of the current address the exact purpose or function of PAC, one can hypothesize that PAC could mediate the binding problem, by coordinating activity in distinct cortical areas in order to both integrate information across cortical areas (addressing the combination problem) as well as to dynamically segregate distinct information processing streams. Studies in nonhuman primates in fact support the potential role of the thalamus in regulating activity across distinct visual cortices, showing that the pulvinar coordinates activity between interconnected cortical areas in an attention-dependent manner (Saalmann et al. 2012). Further task-specific studies in humans are necessary to further characterize the function of PAC.

As would be predicted, thalamic modulation of ongoing spontaneous cortical activity is regulated by structural constraints imposed by direct anatomic connectivity, as defined by MR diffusion connectivity analyses. We did not observe any cases of modulation by second order connections implying the regulatory role of the thalamus is limited to first order connections. This is in fact similar to what Saalmann et al. (2012) found in nonhuman primates indicating pulvinar synchronizing attentional activity between interconnected cortical regions that are also directly connected to the pulvinar thalamus. The thalamus and its subregions are therefore most likely only regulating and coordinating between cortical regions with which there is direct anatomic connectivity and not exerting a regulatory role via distant synapses. Coordination of distant functional cortices may in fact be mediated by intrathalamic connectivity and interplay.

The coordinated dynamic nature of thalamic and cortical signals and the causality analysis presented are strongly suggestive of a regulatory role of thalamocortical PAC. In these spontaneously behaving subjects, significant time-dependent variability is seen in theta power, beta power, thalamocortical theta coherence, intracortical PAC, and thalamocortical PAC. Despite the significant time-dependent variability, these measures demonstrate remarkable covariance and a consistent relationship, suggesting a uniform process of regulation and modulation. Failure of such coupling may account for pathophysiology of various diseases, as has been suggested by the concept of thalamocortical dysrhythmia in which disruption of the normal flow between thalamus and cortex has been contributed to different functional disorders (Henning Proske et al. 2011; Llinas et al. 1999; Llinas and Steriade 2006; Kane et al. 2009; Jones 2010).

While these studies have provided significant insight, there are limitations to the current dataset and analyses. As data were primarily acquired for clinical purposes, sampling frequency was limited to 200 Hz which precludes assessment of coupling of high gamma band frequencies, which are an area of increasing interest. With the current insight, future studies should further investigate cortical–subcortical coupling phenomenon using very high sampling rates to better assess high gamma (>70 Hz) and very high gamma (>200 Hz) activity. Moreover, as a consequence of using data primarily acquired for clinical purposes, current analyses are done in awake and spontaneous behaving patients with epilepsy, rather than with specific tasks. While this limits the extensibility of the conclusions, it does provide insight into normal brain function in an unconstrained system. Task-based analyses will enable investigators to tease out the components of coupling and better understand causality of thalamocortical relationships in a more controlled setting. Because of the need for invasive recordings, there is no foreseeable way to circumvent recordings in diseases patients; these limitations must always be considered in interpreting the generalizability of the current results.

Nevertheless, having established and provided direct evidence for the first time of thalamocortical regulation in spontaneously behaving humans, future studies must further elucidate the regulatory role and dynamic nature of thalamocortical PAC. Questions remain as to the precise role of this regulation. Drawing from studies of PAC in Parkinson's disease (PD), it seems that PAC and therefore thalamocortical PAC and modulation is likely inhibitory in nature, with studies suggesting excess PAC in the motor cortex of patients with PD that improves with therapeutic intervention (de Hemptinne et al. 2013). Moreover, the role of intrathalamic connectivity in coordinating regulation across remote brain regions remains to be better elucidated.

## Conclusion

In accordance with the diffuse connectivity of the thalamus with the cerebral cortex, the current work provides direct causal evidence that the thalamus exerts regulatory control over ongoing cortical activity, at least in part, through a mechanism of PAC, in which low frequency thalamic rhythms, such as theta, directly regulate higher frequency cortical activity in the beta range and possibly beyond, although that remains speculative based on the current data and analysis. Such modulation is spatially specific and structurally constrained, as evidenced by within subject diffusion tractography analysis.

## Supplementary Material

Supplementary Material can be found at http://www.cercor.oxfordjournals.org/online.

## Funding

This study was supported in part by the UCLA Neurosurgery Visionary Ball Fund, the UCLA Scholars in Translational Medicine Program (NP); NIBIB K23EB014326 (NP).

## Notes

Conflict of Interest: None declared.

## References

Axmacher
N
Henseler
MM
Jensen
O
Weinreich
I
Elger
CE
Fell
J
.
2010
.
Cross-frequency coupling supports multi-item working memory in the human hippocampus
.
.
107
:
3228
3233
.
Bazhenov
M
Timofeev
I
M
Sejnowski
TJ
.
1999
.
Self-sustained rhythmic activity in the thalamic reticular nucleus mediated by depolarizing GABAA receptor potentials
.
Nat Neurosci
.
2
:
168
174
.
Behrens
TE
Berg
HJ
Jbabdi
S
Rushworth
MF
Woolrich
MW
.
2007
.
Probabilistic diffusion tractography with multiple fibre orientations: what can we gain?
Neuroimage
.
34
:
144
155
.
Behrens
TE
Johansen-Berg
H
Woolrich
MW
Smith
SM
Wheeler-Kingshott
CA
Boulby
PA
Barker
GJ
Sillery
EL
Sheehan
K
Ciccarelli
O
et al
2003
.
Non-invasive mapping of connections between human thalamus and cortex using diffusion imaging
.
Nat Neurosci
.
6
:
750
757
.
Bertram
EH
Mangan
PS
Zhang
D
Scott
CA
Williamson
JM
.
2001
.
The midline thalamus: alterations and a potential role in limbic epilepsy
.
Epilepsia
.
42
:
967
978
.
Bhattacharya
J
.
2001
.
Reduced degree of long-range phase synchrony in pathological human brain
.
Acta Neurobiol Exp
.
61
:
309
318
.
Bokil
H
Andrews
P
Kulkarni
JE
Mehta
S
Mitra
PP
.
2010
.
Chronux: a platform for analyzing neural signals
.
J Neurosci Methods
.
192
:
146
151
.
Bollimunta
A
Chen
Y
Schroeder
CE
Ding
M
.
2008
.
Neuronal mechanisms of cortical alpha oscillations in awake-behaving macaques
.
J Neurosci
.
28
:
9976
9988
.
Bressler
SL
Richter
CG
Chen
Y
Ding
M
.
2007
.
Cortical functional network organization from autoregressive modeling of local field potential oscillations
.
Stat Med
.
26
:
3875
3885
.
Brovelli
A
Ding
M
Ledberg
A
Chen
Y
Nakamura
R
Bressler
SL
.
2004
.
Beta oscillations in a large-scale sensorimotor cortical network: directional influences revealed by Granger causality
.
.
101
:
9849
9854
.
Canolty
RT
Knight
RT
.
2010
.
The functional role of cross-frequency coupling
.
Trends Cogn Sci
.
14
:
506
515
.
Chen
Y
Bressler
SL
Knuth
KH
Truccolo
WA
Ding
M
.
2006
.
Stochastic modeling of neurobiological time series: power, coherence, Granger causality, and separation of evoked responses from ongoing activity
.
Chaos
.
16
:
026113
.
Contreras
D
Destexhe
A
Sejnowski
TJ
M
.
1996
.
Control of spatiotemporal coherence of a thalamic oscillation by corticothalamic feedback
.
Science
.
274
:
771
774
.
De Hemptinne
C
Ryapolova-Webb
ES
Air
EL
Garcia
PA
Miller
KJ
Ojemann
JG
Ostrem
JL
Galifianakis
NB
Starr
PA
.
2013
.
Exaggerated phase-amplitude coupling in the primary motor cortex in Parkinson disease
.
.
110
:
4780
4785
.
Destexhe
A
McCormick
DA
Sejnowski
TJ
.
1993
.
A model for 8–10 Hz spindling in interconnected thalamic relay and reticularis neurons
.
Biophys J
.
65
:
2473
2477
.
Ding
M
Mo
J
Schroeder
CE
Wen
X
.
2011
.
Analyzing coherent brain networks with Granger causality
.
Conf Proc IEEE Eng Med Biol
.
2011
:
5916
5918
.
Elias
WJ
Zheng
ZA
Domer
P
Quigg
M
Pouratian
N
.
2012
.
Validation of connectivity-based thalamic segmentation with direct electrophysiologic recordings from human sensory thalamus
.
Neuroimage
.
59
:
2025
2034
.
FitzGerald
THB
Valentin
A
Selway
R
Richardson
MP
.
2013
.
Cross-frequency coupling within and between the human thalamus and neocortex
.
Front Hum Neurosci
.
7
:
84
.
Fries
P
.
2005
.
A mechanism for cognitive dynamics: neuronal communication through neuronal coherence
.
Trends Cogn Sci
.
9
:
474
480
.
Golomb
D
Wang
XJ
Rinzel
J
.
1994
.
Synchronization properties of spindle oscillations in a thalamic reticular nucleus model
.
J Neurophys
.
72
:
1109
1126
.
Goncalves
SI
de Munck
JC
Pouwels
PJW
Schoonhoven
R
Kuijer
JPA
Maurits
NM
Hoogduin
JM
Van Someren
EJW
Heethaar
RM
da Silva
FHL
.
2006
.
Correlating the alpha rhythm to BOLD using simultaneous EEG/fMRI: inter-subject variability
.
Neuroimage
.
30
:
203
213
.
Gow
DW
Jr
Keller
CJ
Eskandar
E
Meng
N
Cash
SS
.
2009
.
Parallel versus serial processing dependencies in the perisylvian speech network: a Granger analysis of intracranial EEG data
.
Brain Lang
.
110
:
43
48
.
Haber
SN
Calzavara
R
.
2009
.
The cortico-basal ganglia integrative network: the role of the thalamus
.
Brain Res Bull
.
78
:
69
74
.
He
B
Dai
Y
Astolfi
L
Babiloni
F
Yuan
H
Yang
L
.
2011
.
eConnectome: a MATLAB toolbox for mapping and imaging of brain functional connectivity
.
J Neurosci Methods
.
195
:
261
269
.
Henning Proske
J
Jeanmonod
D
Verschure
PF
.
2011
.
A computational model of thalamocortical dysrhythmia
.
Eur J Neurosci
.
33
:
1281
1290
.
Henriksson
L
Hyvarinen
A
Vanni
S
.
2009
.
Representation of cross-frequency spatial phase relationships in human visual cortex
.
J Neurosci
.
29
:
14342
14351
.
Johansen-Berg
H
Behrens
TE
Sillery
E
Ciccarelli
O
Thompson
AJ
Smith
SM
Matthews
PM
.
2005
.
Functional-anatomical validation and individual variation of diffusion tractography-based segmentation of the human thalamus
.
Cereb Cortex
.
15
:
31
39
.
Jones
EG
.
2010
.
Thalamocortical dysrhythmia and chronic pain
.
Pain
.
150
:
4
5
.
Kaminski
M
Ding
MZ
Truccolo
WA
Bressler
SL
.
2001
.
Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance
.
Biol Cybern
.
85
:
145
157
.
Kaminski
MJ
Blinowska
KJ
.
1991
.
A new method of the description of the information-flow in the brain structures
.
Biol Cybern
.
65
:
203
210
.
Kane
A
Hutchison
WD
Hodaie
M
Lozano
AM
Dostrovsky
JO
.
2009
.
Enhanced synchronization of thalamic theta band local field potentials in patients with essential tremor
.
Exp Neurol
.
217
:
171
176
.
Lakatos
P
Shah
AS
Knuth
KH
Ulbert
I
Karmos
G
Schroeder
CE
.
2005
.
An oscillatory hierarchy controlling neuronal excitability and stimulus processing in the auditory cortex
.
J Neurophys
.
94
:
1904
1911
.
Liu
ZM
de Zwart
JA
Yao
B
van Gelderen
P
Kuo
LW
Duyn
JH
.
2012
.
Finding thalamic BOLD correlates to posterior alpha EEG
.
Neuroimage
.
63
:
1060
1069
.
Llinas
R
Jahnsen
H
.
1982
.
Electrophysiology of mammalian thalamic neurones in vitro
.
Nature
.
297
:
406
408
.
Llinas
RR
Ribary
U
Jeanmonod
D
Kronberg
E
Mitra
PP
.
1999
.
Thalamocortical dysrhythmia: a neurological and neuropsychiatric syndrome characterized by magnetoencephalography
.
.
96
:
15222
15227
.
Llinas
RR
M
.
2006
.
Bursting of thalamic neurons and states of vigilance
.
J Neurophys
.
95
:
3297
3308
.
Lopez-Azcarate
J
Tainta
M
Rodriguez-Oroz
MC
Valencia
M
Gonzalez
R
Guridi
J
Iriarte
J
Obeso
JA
Artieda
J
Alegre
M
.
2010
.
Coupling between beta and high-frequency activity in the human subthalamic nucleus may be a pathophysiological mechanism in Parkinson's disease
.
J Neurosci
.
30
:
6667
6677
.
McCormick
DA
Huguenard
JR
.
1992
.
A model of the electrophysiological properties of thalamocortical relay neurons
.
J Neurophys
.
68
:
1384
1400
.
Miller
KJ
Hermes
D
Honey
CJ
Sharma
M
Rao
RPN
den Nijs
M
Fetz
EE
Sejnowski
TJ
Hebb
AO
Ojemann
JG
et al
2010
.
Dynamic modulation of local population activity by rhythm phase in human occipital cortex during a visual search task
.
Front Hum Neurosci
.
4
:
197
.
Pienkowski
M
Eggermont
JJ
.
2010
.
Nonlinear cross-frequency interactions in primary auditory cortex spectrotemporal receptive fields: a Wiener–Volterra analysis
.
J Comput Neurosci
.
28
:
285
303
.
Pouratian
N
Zheng
Z
Bari
AA
Behnke
E
Elias
WJ
Desalles
AA
.
2011
.
Multi-institutional evaluation of deep brain stimulation targeting using probabilistic connectivity-based thalamic segmentation
.
J Neurosurg
.
115
:
995
1004
.
Saalmann
YB
Pinsk
MA
Wang
L
Li
X
Kastner
S
.
2012
.
The pulvinar regulates information transmission between cortical areas based on attention demands
.
Science
.
337
:
753
756
.
Sarnthein
J
Morel
A
von Stein
A
Jeanmonod
D
.
2005
.
Thalamocortical theta coherence in neurological patients at rest and during a working memory task
.
Int J Psychophysiol
.
57
:
87
96
.
Sarnthein
J
Stern
J
Aufenberg
C
Rousson
V
Jeanmonod
D
.
2006
.
Increased EEG power and slowed dominant frequency in patients with neurogenic pain
.
Brain
.
129
:
55
64
.
Staudigl
T
Zaehle
T
Voges
J
Hanslmayr
S
Esslinger
C
Hinrichs
H
Schmitt
FC
Heinze
HJ
Richardson-Klavehn
A
.
2012
.
Memory signals from the thalamus: early thalamocortical phase synchronization entrains gamma oscillations during long-term memory retrieval
.
Neuropsychologia
.
50
:
3519
3527
.
Tass
P
Smirnov
D
Karavaev
A
Barnikol
U
Barnikol
T
I
Hauptmann
C
Pawelcyzk
N
Maarouf
M
Sturm
V
et al
2010
.
The causal relationship between subcortical local field potential oscillations and Parkinsonian resting tremor
.
J Neural Eng
.
7
:
16009
.
Tort
AB
Komorowski
R
Eichenbaum
H
Kopell
N
.
2010
.
Measuring phase-amplitude coupling between neuronal oscillations of different frequencies
.
J Neurophys
.
104
:
1195
1210
.
Tort
AB
Komorowski
RW
Manns
JR
Kopell
NJ
Eichenbaum
H
.
2009
.
Theta-gamma coupling increases during the learning of item-context associations
.
.
106
:
20942
20947
.
Walton
KD
Llinas
RR
.
2010
.
Central pain as a thalamocortical dysrhythmia: a thalamic efference disconnection?
In:
Kruger
L
Light
AR
, editors.
Translational pain research: from mouse to man
.
Boca Raton, FL
.
Wang
X
Chen
Y
Bressler
SL
Ding
M
.
2007
.
Granger causality between multiple interdependent neurobiological time series: blockwise versus pairwise methods
.
Int J Neural Syst
.
17
:
71
78
.
Wang
X
Chen
Y
Ding
M
.
2008
.
Estimating Granger causality after stimulus onset: a cautionary note
.
Neuroimage
.
41
:
767
776
.
Zhang
L
Chen
G
Niu
R
Wei
W
Ma
X
Xu
J
Wang
J
Wang
Z
Lin
L
.
2012
.
Hippocampal theta-driving cells revealed by Granger causality
.
Hippocampus
.
22
:
1781
1793
.