Abstract
The amplitude of high broadband activity in human cortical field potentials indicates local processing and has repeatedly been shown to reflect motor control in the primary motor cortex. In a group of male and female subjects affected by essential tremor and undergoing deep brain stimulation surgery, ventral intermediate nucleus low-frequency oscillations (<30 Hz) entrain the corticomotor high broadband activity (>40 Hz) during rest, relinquishing that role during movement execution. This finding suggests that there is significant cross-rhythm communication between thalamocortical regions, and motor behavior corresponds to changes in thalamocortical phase-amplitude coupling profiles. Herein, we demonstrate that thalamocortical coupling is a crucial mechanism for gating motor behavior.
SIGNIFICANCE STATEMENT We demonstrate, for the first time, how thalamocortical coupling is mediating movement execution in humans. We show how the low-frequency oscillation from the ventral intermediate nucleus, known as the motor nucleus of the thalamus, entrains the excitability of the primary motor cortex, as reflected by the phase-amplitude coupling between the two regions. We show that thalamocortical phase-amplitude coupling is a manifestation of a gating mechanism for movement execution mediated by the thalamus. These findings highlight the importance of incorporating cross-frequency relationship in models of motor behavior; and given the spatial specificity of this mechanism, this work could be used to improve functional targeting during surgical implantations in subcortical regions.
Introduction
Coupling of local field oscillations in distal regions of the brain may play an important role in the communication and control dynamics between brain structures (Fries, 2005; Canolty and Knight, 2010; Miller, 2010; Canolty et al., 2012; Voytek et al., 2013; Maris et al., 2016). Current evidence associates lower frequency oscillations (LFOs; <30 Hz) with subcortical regions projecting onto the cortex (Osipova et al., 2008; Roux et al., 2013; Malekmohammadi et al., 2015) and with widespread distributed cortical modulations (Osipova et al., 2008; Hindriks et al., 2014), suggesting that LFOs can be a network-scale driver for neural excitability (Canolty et al., 2012; Voytek et al., 2013; Maris et al., 2016). In contrast, high broadband activity (HBA; >40 Hz) is associated with local cortical processing (Miller et al., 2007, 2009, 2010b; Jiang et al., 2015) and has been shown to be correlated with local neural firing rates (Ray et al., 2008; Manning et al., 2009; Miller, 2010). Moreover, previous studies have highlighted the role of LFOs, in specific sub-bands, such as theta (Mormann et al., 2005; Canolty et al., 2006; Sirota et al., 2008; He et al., 2010; Voytek et al., 2010; Lega et al., 2016), alpha (and mu) (Osipova et al., 2008; Voytek et al., 2010; Yanagisawa et al., 2012), and beta (de Hemptinne et al., 2015; Malekmohammadi et al., 2015; Wang et al., 2016); and in cross-frequency coupling as a communication mechanism for cognitive processing during visual, motor, attention, memory, or sleep tasks (Mormann et al., 2005; Canolty et al., 2006; Osipova et al., 2008; Cohen et al., 2009; Canolty and Knight, 2010; Miller et al., 2010a; Yanagisawa et al., 2012; Roux et al., 2013; Lega et al., 2016), within the same local neuronal cluster (Canolty et al., 2006; Yanagisawa et al., 2012; Shreve et al., 2017) or across different regions (Roux et al., 2013; Shimamoto et al., 2013; Malekmohammadi et al., 2015, 2018). Herein, we focus on cross-frequency coupling between LFOs and HBA representing a mechanism by which slow oscillations influence the firing rate of distal neurons (Volgushev et al., 1998; Canolty et al., 2006; Osipova et al., 2008; Roux et al., 2013), through a gating mechanism (Miller et al., 2012; Yanagisawa et al., 2012; Pascucci et al., 2018) that mediates task performance. We are particularly interested in the phase-amplitude coupling (PAC) relationship between the LFO and HBA of two distal nodes in the thalamocortical network, respectively, the ventral intermediate nucleus (VIM) of the thalamus and the primary motor cortex (M1). While it has been previously shown that PAC is involved in motor execution within the somatosensory cortex (Miller et al., 2012; Yanagisawa et al., 2012) and between the GPi and M1 (AuYong et al., 2018), we show that this relationship further expands in the motor network, involving the VIM, referred to as the motor nucleus of the thalamus (Molnar et al., 2005).
We hypothesize that the VIM LFO acts as a gating mechanism for movement execution by selectively modulating the motor cortical HBA, which in turn facilitates the execution of motor commands (Yanagisawa et al., 2012). To test this, we recorded local field potentials (LFPs) from cortical (M1) and subcortical (VIM) regions in human subjects affected by essential tremor undergoing deep brain stimulation (DBS) surgery. The subjects enrolled in our cohort exhibited intraoperative suppression of tremor symptoms following DBS electrode insertion (Morishita et al., 2010). We found that VIM-M1 PAC was present at rest, while it was suppressed during movement execution.
Materials and Methods
Participants.
Eleven subjects (7 male and 4 female) were enrolled at the University of Florida Fixel Center for Neurological Disorders with medically refractory ET scheduled for DBS surgery. Informed consent was obtained at least 1 d before the surgery. Participation in the study entailed the temporary implantation of an additional intraoperative monitoring subdural electrode strip, used in epilepsy and tumor resection surgeries to perform functional mapping, which was placed above the hand motor cortex and removed before the end of the surgery. The study protocol was in agreement with the Declaration of Helsinki, approved by the University of Florida Institutional Review Board, and did not interfere with the outcome of the clinical procedure. The subjects were enrolled in the study only for the duration of the surgery.
Surgical procedure.
The electrode implants consisted of 2 leads: (1) Medtronic model 3387 4-contact DBS lead placed in the VIM contralateral to the limb with dominant tremor as part of standard clinical care; and (2) 4- or 6-contact platinum intraoperative monitoring subdural strip placed acutely over hand motor cortex (model 2110–06-077, PMT). To localize the VIM target for this study, T1 structural MRI was fused to a stereotactic head CT and coregistered to a University of Florida-developed 3D deformable brain atlas, which was adapted to match the subject brain to predict the location of the VIM region as accurately as possible (Sudhyadhom et al., 2012). After choosing a safe trajectory, a burr hole was drilled in the skull at the trajectory entry site. Anesthetic was delivered locally at the site of the implantation, allowing the subject to be awake during the surgery. A microelectrode was advanced through the brain along the planned trajectory using a micropositioner to allow for intermittent single-unit recording with an FHC 4000 LP+ system (FHC) and to precisely identify the VIM/nucleus ventralis caudalis boundary. Before the initiation of any functional mapping and recording, the subject sedatives were stopped. At each stop along the trajectory, a trained neurologist evaluated the neural firing in response to stimuli/commands (touch, volitional movements, light cues) to verify whether the brain region activity matched the expected firing patterns. When the target was localized, the microelectrode was substituted with the DBS lead. Targeting was verified through intraoperative macrostimulation for effectiveness of tremor suppression and side effects (e.g., paresthesia). The cortical strip was inserted from the same burr hole made for the VIM lead. The cortical strip was placed proximal to the hand knob (Yousry et al., 1997; Panov et al., 2017), over the hand motor cortices, based on MRI-CT fusion images and the execution of a modified real-time functional mapping protocol (Hill et al., 2012), which verified correct cortical positioning. The methodology checks for any difference in the HBA power compared with baseline (using the coefficient of determination, r2) (Schalk et al., 2004). The contacts that show the highest r2 during a specific task, compared with baseline, are classified to be related to the task itself, given that HBA provides accurate functional localization (Miller et al., 2007, 2009, 2010b; Jiang et al., 2015). At the end of the research recordings, the subject was sedated, the subdural strip was removed (http://www.clinicaltrials.gov/: NCT02712515), and standard surgical care resumed. Two of the subjects received chronic subdural strips connected to an FDA IDE-approved bidirectional neurostimulator (Activa PC+S, Medtronic; http://www.clinicaltrials.gov/: NCT02649166).
Experimental design.
The subjects were instructed on how to execute the tasks before the surgery. During surgery, the subjects were awake for clinical purposes and remained awake during the research protocol. The subjects were in a lying down position, from which they could look at a monitor that displayed the task instructions/cues. The subjects performed 3 tasks: (1) baseline, in which the subject was asked to remain still in a resting position (4 ± 1 min), without contracting muscles to maintain any posture; (2) cued-go hand opening and clasping, in which the subject was cued to their left or right hand in a random order. When cued, the subject had to prepare to move a specific hand (CUE, for 5 ± 1 s). Next, the subject was cued to perform the movement (GO CUE, for 6 ± 1 s). Each moving phase was stopped by a rest cue, where the subject returned to the resting position used during the baseline task for 7 ± 2.5 s (see Fig. 1B); (3) cued cup grasping, in which the subject was asked to reach a cup (CUE every 12 ± 4 s), bring it toward their mouth and bring it back to the hand of the operator/researcher, in the original position. Only one hand was used at a time. No DBS was delivered during the execution of these tasks (see Fig. 1A); a subset of subjects executed the hand task (2) with the same structure as (3), without the cue for which hand to use, only the movement CUE (see Fig. 1A).
The tasks were performed after the placement of both electrode leads. Both leads were connected to a g.HiAmp external bioamplifier (Guger Technologies). The recordings were stored in a monopolar configuration by using as reference and ground the tissue near the burr hole. EMG and inertia (accelerometer, gyroscope, and magnetometer) data were acquired, respectively, at 1111.11 and 148.25 Hz with the Trigno Wireless EMG System (Delsys). The EMG/inertia sensors were placed on the upper limbs of the subject (hand dorsum, extensor carpi ulnaris, flexor carpi radialis, bicep brachii on both arms) to detect movement onset. This dataset also verified the suppression of tremor in our subjects after DBS lead implantation. Video data were acquired at 24 Hz to capture behavioral activity. All the different sources were aligned at the beginning of the recording session with train of square pulses (“sync trigger”) delivered to a dedicated shared channel. Each acquisition source has one, while the camera has an LED in front of it, matching the “sync trigger” status. Signals from different systems were aligned based on these square pulses offline.
Analysis.
All analyses were executed within the MATLAB environment (The MathWorks). To speed up the processing, the algorithm was scaled on the UF Hipergator HPC (Research Computing, University of Florida, Gainesville, Florida) and on a GPU (Nvidia, CUDA; https://developer.nvidia.com/accelerated-computing-training).
Signal processing of neural signals.
As a first preprocessing step, the LFP channels (stored in a monopolar configuration) were inspected for artifacts and SNR through visual inspection of the raw time series, power spectral densities (PSDs) across tasks and spectrograms. The raw channels were then recombined for monopolar and bipolar analyses. For monopolar analysis, a common average reference filter was applied to each lead for denoising. For bipolar analyses, combinations with one neighborhood and two neighborhood Laplacian referencing were performed on each lead. Line interference at 60 Hz was rejected only for the monopolar analysis with a notch filter. An additional high-pass filter at 1 Hz was used to eliminate the time series drift and to suppress the pressure artifact present only in the thalamic LFP, which is caused by blood flow variations due to heartbeats. No channels were dropped from the subsequent analyses. PSDs and spectrograms were computed by using an autoregressive model (Schalk et al., 2004; Stoica and Moses, 2005), with a frequency resolution of 1 Hz.
Movement events.
Movement events within each block were marked via visual inspection of accelerometer EMG and video data, obtaining 6 total event types: baseline, baseline during task (intermovement intervals), cue for right hand movement, cue for left hand movement, right hand movement, and left hand movement. The same subdivision was used for both cued hand opening and clasping and cued cup grasping.
Imaging.
MRI data were collected 1 d before surgery with a 3T Siemens Magnetom Prisma scanner and a 3T Siemens Verio scanner. Cortical reconstruction and volumetric segmentation were performed with the Freesurfer (Fischl, 2012) image analysis suite after surgery, which is documented and freely available for download online (http://surfer.nmr.mgh.harvard.edu/). This allowed to obtain the anatomical positioning of the cortical strip contacts. For the subcortical structures (white matter and deep gray matter), the LeadDBS v2 (Horn et al., 2019) toolbox was used, which is documented and freely available for download online (http://www.lead-dbs.org/). This allowed to normalize and bring the lead contacts in the Human Motor Thalamus Atlas (Ilinsky et al., 2018). The normalization was executed with the ANTs package (Avants et al., 2008, 2011).
PAC.
Coupling within the same contact and across contact combinations (and therefore across structures) was investigated with PAC by using the Tort method, a Kullback-Lieber-based method (Tort et al., 2010). Phase and amplitude were extracted from one (within PAC) or two channels (across PAC), respectively, defined as phase and amplitude time series. The raw time series selected for PAC were band-passed within the ranges of frequencies specified for phase and amplitude (e.g., respectively, 4–70 Hz, with 2 Hz bin, and 4–300 Hz, with 2 Hz bin) with a linear FIR filter (Delorme and Makeig, 2004). This generated two sets of band-passed signals (one for phase and one for amplitude). Hilbert transform was applied to both sets. The phase for the phase set was obtained by calculating the phase angle from the Hilbert-transformed phase set. The amplitude for the amplitude set was obtained by calculating the absolute value from the Hilbert-transformed amplitude set. The modified Kullback-Lieber method was used to compute the modulation index (MI) for each frequency for phase and frequency for amplitude combination. For statistical analysis, a surrogate population of MIs was created by randomly time shifting the frequency for amplitude vector 400 times. The p value was calculated by considering the percentage of surrogate values above the “true” MI. A single thalamic–cortical combination that coupled the most (highest MI) was chosen for each subject. In addition, the frequency for phase that coupled the most with HBA was computed. The instantaneous PAC was computed by estimating the PAC within a time window, specifically between the previously mentioned frequency for phase and the frequencies for amplitude (e.g., 4–300 Hz). No differences were noted between PAC computed from monopolar and bipolar combinations (comparison not shown). Bipolar combinations were preferred and used throughout this work to avoid common reference contamination.
Statistical analysis.
The MATLAB statistical package was used (R2018a). Unless indicated otherwise, all correlations are Pearson and all p values are two-tailed. Wilcoxon signed-rank test was used in paired tests. Unless indicated otherwise, a p value <0.05 was considered statistically significant. For the PAC calculations, a Benjamini and Hochberg false discovery rate (FDR) procedure (Benjamini and Hochberg, 1995; Benjamini et al., 2001; Genovese et al., 2002) was used.
Results
Electrophysiology data were recorded in a total of 11 subjects affected by essential tremor intraoperatively during DBS surgery while performing tasks (Fig. 1A). Data were recorded during rest (Fig. 1B) and motor execution (Fig. 1C). Clinical characteristics of the subjects are provided in Table 1. Raw traces show the difference between the thalamic and cortical signal amplitudes, approximately by one order of magnitude. Spectral analyses were consistent with the human motor cortex literature (Freeman et al., 2000; Miller et al., 2007, 2010b) and showed an increase in HBA during upper limb movements contralateral to the site of implantation and a decrease in the mu and beta LFOs (∼8–30 Hz) range during ipsilateral or contralateral movements (Fig. 1C, left), with a more prominent decrease during contralateral movements. In the VIM recordings (Brücke et al., 2013), significant changes were observed only in the lower bands (4–30 Hz, LFO), with a less prominent power decrease during movements ipsilateral to the side of implantation compared with contralateral movements (Fig. 1C, right).
Task timeline, raw and spectral analyses. A, Task timeline of the recordings. The top timeline did not have a preparation cue. Cued cup reaching task was administered only in the top timeline format. Boxes with a dashed border represent the uncued actions the subject took during the task. B, C, E, Data collected in an ET DBS implantation surgery along with contralateral and ipsilateral hand accelerations. The subject was asked to clench his/her hand repetitively when CUED (yellow). The GO CUE event (green) represents when the subject was asked to move their hand, and the MOVING event (violet) represents when the subject moved their hand. B, Baseline/rest data were collected in a separate run. D, Normalized PSD for both contralateral and ipsilateral movement for LFO (8–30 Hz) and HBA (>40 Hz) in motor cortex (Wilcoxon signed-rank test, LFO p = 0.001 and HBA p = 0.003 for baseline vs contralateral movement, LFO p = 0.002 and HBA p = 0.413 for baseline vs ipsilateral) and VIM (Wilcoxon signed-rank test, LFO p = 0.001 and HBA p = 0.001 for baseline vs contralateral movement, LFO p = 0.002 and HBA p = 0.067 for baseline vs ipsilateral). E, Zoomed section of LFP from C during contralateral hand movement. F, Thalamocortical coherence in baseline/rest on the right, normalized coherence for both contralateral and ipsilateral movement on the left (Wilcoxon signed-rank test, p = 0.002 for baseline vs contralateral movement, p = 0.0137 for baseline vs ipsilateral). G, Thalamocortical PLV in baseline/rest on the right, normalized PLV for both contralateral and ipsilateral movement on the left (Wilcoxon signed-rank test, p = 0.001 for baseline vs contralateral movement, p = 0.003 for baseline vs ipsilateral). *p < 0.01.
Subject demographic in this study populationa
Thalamocortical coherence and phase synchrony changes during movement execution
A large decrease in thalamocortical coherence was observed during movement in the LFO range (in respect to baseline, Fig. 1F) (Wilcoxon signed-rank test, p = 0.002 for baseline vs contralateral movement, p = 0.0137 for baseline vs ipsilateral). Coherence was spectrally concentrated in the most desynchronized band in PSD analysis in both M1 and VIM. To verify whether the aforementioned coherence modulation could be explained by phase synchrony, we computed the phase locking value (PLV) (Lachaux et al., 1999; Mormann et al., 2000) between the VIM and M1 for these spectral bands, as PLV disregards amplitude information. A significant decrease in thalamocortical PLV was elicited by movement execution (Fig. 1G; Wilcoxon signed-rank test, p = 0.00098 for baseline vs contralateral movement, p = 0.0029 for baseline vs ipsilateral). The frequency range that yielded highest phase locking matched the coherence findings.
PAC coupling during rest and suppression during movement execution
The mechanism of interest in this work is the modulation between thalamocortical nodes, in the form of cross-frequency coupling. Specifically, we focused on PAC pattern changes during hand movements compared with rest. The relationship between thalamic (VIM) LFO and cortical (M1) HBA is highlighted by the instantaneous PAC (Fig. 2C,H), wherein we observed decreased coupling during ipsilateral and contralateral movements compared with baseline. This suggests that PAC is the manifestation of a gating mechanism for movement. Across subjects (n = 11), we observed consistent coupling between thalamic LFO phase (9–15 Hz) and cortical HBA amplitude (>40 Hz) during rest (Fig. 3A–D,I–L), with significant phase-amplitude combinations across at least the 70% of the subjects (p < 0.01). The PAC p values were adjusted with Benjamini and Hochberg FDR procedure (Benjamini and Hochberg, 1995; Benjamini et al., 2001; Genovese et al., 2002) to address the multiple-comparison problem. PAC within the motor cortex was also consistent with previous findings (Miller et al., 2012; Murta et al., 2017). Both corticocortical and thalamocortical PAC were absent during movement of either hand at a group level; however, in 3 of 11 subjects, ipsilateral coupling decreased without being completely suppressed (Fig. 3E,F,M,N). The absence of statistically significant cortical LFO to thalamic HBA coupling supports the unidirectionality of the thalamic (LFO) to cortical (HBA) mechanism. Hence, these findings support modulation of local motor processing through thalamocortical coupling, with the cortical motor areas kept in a “dynamically suppressed” state (Miller et al., 2010a, 2012; Bonnefond et al., 2017; Pascucci et al., 2018), which can quickly revert to a “processing state” to responsively execute any motor command. To support this hypothesis, we show that the HBA amplitude was the highest following the trough/peak of LFO bursts (Fig. 3G,H) (Miller et al., 2010a).
Thalamocortical spectral and PAC modulations during a cued hand movement task. Data collected in an ET DBS implantation surgery along with (A) contralateral and (E) ipsilateral hand accelerations. The subject was asked to clench his/her hand repetitively when CUED (yellow). The MOVING event (green) represents when the subject moved their hand. The normalized spectrograms of (B) hand motor cortex (MC) and (D) VIM show movement-related spectral modulations, with increased HBA in the MC during contralateral movements, and decreased LFO amplitudes in MC and VIM during both hand movements. C, Instantaneous PAC between the phase of a subset of VIM LFO (12–20 Hz) and the amplitude of a range of oscillations (4–200 Hz) of the MC. The subset of LFO for the phase was chosen based on the frequencies for phase that coupled the most during the baseline task. The normalized ERD-ERS spectrograms aligned to contralateral movement onset and end show the increased HBA in MC (F) and decreased LFO amplitudes in both MC (F) and VIM (H). G, Event-related PAC of the instantaneous PAC of C, aligned to contralateral movement onset and end. The thalamocortical PAC decreases during movement. PSD of ipsilateral MC (I) and VIM (J), showing clear desynchronization for both ipsilateral and contralateral movement (VIM at ∼8–50 Hz, MC at ∼8–20 Hz). There is less modulation for ipsilateral movement. K, Coherence between VIM and ipsilateral M1. There are clear coherence peaks during baseline (10, 20 Hz) and ipsilateral movement (20 Hz).
A–D, Columns represent PAC patterns between different brain regions: P, Phase; A, amplitude. PAC combinations that are statistically significant (p = 0.01, Benjamini and Hochberg FDR-adjusted) across at least 70% of the subjects are contoured in black. Different task conditions are presented in the rows. A, B, Columns represent corticocortical and thalamocortical PAC, respectively, during baseline/rest, which decreases during ipsilateral and contralateral movements. E, F, Columns represent the corticocortical and thalamocortical coupling, respectively, in a subset of patients (3 of 11) with PAC during ipsilateral movement. G, Single-subject phase coupling palette is constructed to show the coupling at different frequencies for intraregions and interregions (e.g., thalamocortical). The frequency for phase selected is the one that coupled the most within thalamocortical PAC profiles (15 Hz in this subject). The trough of the LFO is at 0 radians. A consistent spectral change (HBA > 40 Hz from cortical contacts) is phase-locked with the thalamus LFO, with the increase in HBA power after the trough of the LFO (at 0 phase). H, Polar representation of the PAC palette that shows the most coupling (>40 Hz). I–L, Columns represent PAC combinations that are statistically significant (p = 0.01, Benjamini and Hochberg FDR-adjusted) in the A–D columns. The combinations that are significant across at least 70% of the subjects are contoured in black. Different task conditions are presented in the rows. M and N show respectively the cortico-cortical and thalamocortical significant coupling combinations in the subset of patients from E–F columns.
Spectral and PAC relationship shows how the LFO modulator mediates motor execution
There was a clear match between the VIM LFO frequency with the largest decrease in power during contralateral movement execution (normalized to baseline) and the frequency for phase of the highest thalamocortical PAC (VIM LFO to M1 HBA) during baseline (method at Fig. 4, results at Fig. 5A, top, B). A similar relationship involving the coherence and PAC can also be seen (Fig. 5A, bottom, C). These findings reinforce the mechanistic relationship and demonstrate how the thalamic activity modulates the cortex through coupling at a specific phase frequency for each subject (<30 Hz). The PAC at these specific phase frequencies significantly decreased during contralateral and ipsilateral movement for each subject (Fig. 5D, Wilcoxon signed-rank test, p = 0.002 for baseline vs contralateral movement, p = 0.004 for baseline vs ipsilateral) and did not occur at any other point in the PAC comodulogram (Fig. 6A,B, green dashed line).
Analysis of link between spectral content and PAC. A, Raw signals from cortex and VIM are processed. B, PSD is computed for VIM. The first frequency of interest is the one with the highest decrease in power during contralateral movement (green arrow). C, PAC is computed between VIM LFO and cortical HBA during a baseline/rest task. D, The average with respect to the frequency for amplitude (range 40–250 Hz) is z-normalized and used to find the frequency for phase with the highest PAC (green arrow), which will be the second frequency of interest. E, This procedure is repeated for each subject to the obtain the scatterplot of PAC max peak frequency and PSD max power drop/decrease, and the spread of their difference.
Relationships between spectra and PAC. A, Top, Group analysis shows relationship between frequency with maximum PSD drop (LFO) in VIM during contralateral movement execution compared with baseline, and the phase frequency of the highest thalamocortical PAC (VIM LFO to cortical HBA) during baseline. Bottom, Group analysis shows relationship between frequency with maximum/peak thalamocortical coherence (LFO) and the phase frequency of the highest thalamocortical PAC (VIM LFO to cortical HBA), both during baseline. B, The highly limited spread of the difference shows how the frequency that decreases the most in the VIM during movement matches the frequency that couples the most with the broadband high-frequency cortical activity. C, The highly limited spread of the difference shows how the frequency that has the most thalamocortical coherence in the LFO range during baseline matches the frequency that couples the most with the broadband high-frequency cortical activity. D, Thalamocortical PAC during movement normalized to baseline presents a statistically significant decrease (Wilcoxon signed-rank test, p = 0.002 for baseline vs contralateral movement, p = 0.004 for baseline vs ipsilateral). The comparison is made at the phase frequency that couples the most during baseline/rest. *p < 0.01.
Spatial specificity of the thalamocortical PAC. A, Exemplary normalized PSD profiles for contralateral and ipsilateral movement tasks (B) with their respective PAC profiles: P, Phase; A, amplitude. PAC between the VIM and cortex is specific to the precentral gyrus contacts (B, third column, first row). Posterior cortical channels (channels 1–2, 2–3) do not exhibit local or distributed PAC. Green dashed line that overlays A, B indicates the relationship between the frequency of maximum PSD drop (in LFO) and the phase frequency of the PAC maximum peak. C, Image reconstruction showing the subdural cortical electrode positioning (4 contact configuration) and (D) DBS lead positioning (model 3387). Red represents the VIM electrodes used for the coupling in B. C, D, Based on presurgical MRI and postoperative CT. E, Group-level PAC patterns between postcentral gyrus (PostC) and VIM (combination specified with phase [P]and amplitude [A]). PAC combinations that are statistically significant (p = 0.01) are FDR (Benjamini and Hochberg) adjusted. Different task conditions are presented in the rows. F, PAC combinations statistically significant across subjects. VIM-PostC (somatosensory) combination did not show a consistent group-level significance. G, PAC modulation normalized to baseline did not show statistically significant difference across different task conditions (Wilcoxon signed-rank test, p = 0.770 for baseline vs contralateral movement, p = 0.375 for baseline vs ipsilateral).
Thalamocortical PAC is spatially selective for motor cortices
Moreover, thalamocortical motor modulation was a spatially specific mechanism, as it existed between VIM and motor areas (precentral gyrus), but not between VIM and the primary somatosensory (S1) cortex (postcentral gyrus) (Fig. 6B, single subject; Fig. 6E,F, group level), even though both cortical areas of interest (M1 and S1) revealed spectral modulation during movement execution (Fig. 6A). Group-level PAC modulation normalized to baseline between VIM and postcentral gyrus (S1) did not yield a statistically significant difference across different task conditions (Wilcoxon signed-rank test, p = 0.770 for baseline vs contralateral movement, p = 0.375 for baseline vs ipsilateral), contrary to what was shown for VIM and precentral gyrus (M1) (Figs. 3A,B, 5D). The corticocortical coupling observed in Figure 3A could be explained by mirroring of the thalamic LFO phase in M1 as a mechanism to suppress local activity, as there is high thalamocortical phase locking as shown by coherence and PLV, in the LFO range (Fig. 2K). The corticothalamic (Fig. 3C) and thalamothalamic (Fig. 3D) PAC did not reveal any group-level relationship.
PAC modulation shows no difference between motor behaviors
Group analysis of thalamocortical coupling (Fig. 7A,B) showed statistically significant modulation during upper limb movement both with cued-go hand opening and closing (Wilcoxon signed-rank test, p = 0.008 for baseline vs contralateral movement, p = 0.008 for baseline vs ipsilateral), and with cued cup reaching movement (Wilcoxon signed-rank test, p = 0.008 for baseline vs contralateral movement, p = 0.008 for baseline vs ipsilateral). The two movement-related PAC modulations (Fig. 7C) did not show any significant difference between the two motor tasks (hand movement and cup grasping) tested in this study (Spearman r = 0.810, p = 0.022).
PAC modulation shows no difference between motor behaviors. A, Group analysis of thalamocortical PAC during hand movement normalized to baseline presents a statistically significant decrease (Wilcoxon signed-rank test, p = 0.008 for baseline vs contralateral movement, p = 0.008 for baseline vs ipsilateral). B, Group analysis of thalamocortical PAC during cup reaching normalized to baseline presents a statistically significant decrease (Wilcoxon signed-rank test, p = 0.008 for baseline vs contralateral movement, p = 0.008 for baseline vs ipsilateral). A, B, The comparisons are made at the phase frequency that couples the most during baseline/rest. The group analysis is executed in the subject that executed both tasks. C, Spearman correlation shows similarity between the PAC profiles in two different motor behaviors (r = 0.80952, p = 0.022). *p < 0.01.
Waveshape analysis does not indicate spurious thalamocortical PAC
Previous work from other groups (Cole and Voytek, 2017; Cole et al., 2017) has underlined how the shape of the oscillation/waveform can influence the presence of PAC. An excessive inequality between trough and peak sharpness could lead to spurious PAC, especially when the phase and amplitude are computed from the same channel (e.g., corticocortical PAC in Fig. 3A). Our work focused on internodal PAC (thalamus and cortex), which is less affected by this issue (Cole and Voytek, 2017). The analysis of the wave shape showed that there was no statistical difference between trough and peak sharpness in either the M1 or VIM in any subject (Fig. 8 shows results for the same subject in Fig. 6).
Single-subject wave shape analysis of alpha-beta oscillations. A, Filtered signal (bandpass filtered 6–40 Hz) parsed to find peaks (red crosses), throughs (blue crosses), and cross points (blue and green squares). The highlighted blocks (blue, violet, red) represent the found bursts based on cycle thresholding (amplitude >0.1 of z-normalized values). Top, Burst/cycle analysis from M1. Bottom, From VIM. B, Sharpness distribution for trough and peak of alpha-beta oscillations. C, Distributions of rise-decay and peak-trough period symmetry of alpha-beta oscillations. Although the peak-trough symmetry median for M1 deviated from 0.5 (by 0.023242), this shift does not affect the presence of PAC, but only the phase at which PAC occurs. *p < 0.01.
Discussion
The mechanism by which PAC acts as a communication modality between distal nodes likely relies on LFOs, which have been shown to better propagate through tissue, to influence the membrane potential of the “receiver” neural cell and consequently to affect its excitability (Volgushev et al., 1998; Canolty et al., 2006; Osipova et al., 2008; Roux et al., 2013). We demonstrate in an ET population how this distal communication affects the network involving VIM and M1. The role of other regions, such as cerebellum, which directly projects onto the thalamus and is important for aspects of motor control, remains unexplored. Herein we did not explore PAC as a descriptor of the pathophysiology of ET. Instead, we focused on PAC as the expression of the underlying gating mechanism for movement execution.
Our findings are compatible with previously reported coupling between LFO and HBA within the same cortical node (Canolty et al., 2006; de Hemptinne et al., 2013, 2015). The thalamocortical coupling mechanism presented herein, together with the phase locking between the thalamus and cortex LFOs, explains this PAC between the receiving node LFO phase and its own HBA amplitude. This finding could also represent an additional responsive mechanism to suppress motor execution that is triggered by thalamic input (Bonnefond et al., 2017). All subjects had a complete absence of corticothalamic coupling (M1 phase-coupled with VIM amplitude), showing that the mechanism is not symmetric. Other studies have reported a similar mechanism between GPi-sensorimotor cortex (Malekmohammadi et al., 2018), and STN-sensorimotor cortex (Shimamoto et al., 2013).
Ipsilateral movement led to complete PAC suppression in 8 of 11 subjects, whereas the remainder had significant coupling (p = 0.01, Benjamini and Hochberg FDR-adjusted), resulting in a PAC profile similar to baseline PAC (Fig. 3E,F, second row). These findings support the theory that the release of the gating mechanism is a necessary process for cortical excitability, but it is not sufficient for movement execution. The lack of PAC between the thalamus and somatosensory cortex reinforces this mechanism involvement solely in the execution of the movement, disregarding sensory input. In addition, the lack of significant differences between the PAC decrease across different movement tasks suggests that this mechanism is unrelated to the type or complexity of the task at hand.
In addition, other ongoing studies in our group (Fig. 9) have shown the lack of phase-amplitude coupling in the centromedian thalamus, which is the DBS target used in subjects affected by Tourette syndrome, even though it is adjacent to the VIM. Hence, considering the spatial selectivity of this mechanism, specifically only between the primary motor cortex and the VIM (and not other thalamic structures), this work would contribute to the improvement of functional targeting during surgical implantation in subcortical regions.
PAC patterns at rest from a single subject affected by Tourette and implanted in M1 and centromedian (CM) thalamus. PAC combinations that are statistically significant (p = 0.025, Benjamini and Hochberg FDR-adjusted) are contoured in black. While there is consistent corticocortical PAC (A), the thalamocortical PAC (B) between CM-thalamus and M1 has a clear different pattern (no consistent PAC at any phase for frequency) compared with the VIM-M1 PAC.
Previous research on thalamocortical connectivity has shown low-frequency phase-to-high-frequency amplitude coupling during rest, such as theta phase-beta amplitude (Malekmohammadi et al., 2015) and alpha phase-broadband gamma amplitude (with magnetoencephalography) (Roux et al., 2013). However, wider frequency spectra and behavioral related changes have yet to be explored in detail. Our study demonstrates that (1) cortical HBA is modulated by a VIM LFO (but not the opposite); (2) VIM-cortex PAC is likely a gating mechanism for movement execution; (3) thalamocortical PAC is suppressed during movement; (4) PAC is spatially selective, even more than spectral power profiles suggest (somatosensory nodes with LFO desynchronization during movement do not necessarily elicit PAC); and (4) ipsilateral and rest PAC reflect motor laterality (in a subset of subjects).
Limitations
There were limitations in the extent of the possible recordings due to temporal and clinical constraints during DBS surgical implantation. The baseline recordings were 3–4 min, whereas movement recordings consisted only of 12–18 trials per task, lasting on average for 6 min. Data presented in this study were collected after insertion of DBS lead in the VIM, which is known to cause temporary symptom relief (Koop et al., 2006; Morishita et al., 2010) and was experienced by our subjects as well.
Footnotes
This work was supported by National Institutes of Health, National Institute of Neurological Disorders and Stroke Grant UH3NS09553 to K.D.F. and A.G., De Luca Foundation to S.C., and National Institutes of Health, National Institute of Neurological Disorders and Stroke Grant T32NS082128 to S.C. We thank Robert Eisinger, Jackson Cagle, Sarah Long, and Brandon Parks for helpful comments on the manuscript; and Julie Segura for administrative support.
The authors declare no competing financial interests.
- Correspondence should be addressed to Aysegul Gunduz at agunduz{at}ufl.edu
